order parameters

[Vasiliy Triandafilidi]

Dear MDanalysis community.

Are there are any written programms of the orientational order for polymers?
I am new to MDanalysis and I wrote my code but it looks ugly and I was wondering if anyone can share their experience.

Thank you,
Vasiliy.

[Tyler Reddy]

Are there are any written programms of the orientational order for polymers?

Is this substantially different from the S^2 order parameter (for bond vectors) frequently calculated using NMR spectroscopy? I know plenty of people have written code to calculate that, although I don't think there is anything built into MDAnalysis that does it at the moment.

I am new to MDanalysis and I wrote my code but it looks ugly and I was wondering if anyone can share their experience.

It might be useful to share the pertinent snippets of your code and provide a reference for the definition of the orientational order parameter (if it is different from the NMR version). If you'd like it to be built in to MDAnalysis you might want to file an Issue that describes these things, and includes an example file / data of some sort and the value of the order parameter (so it can be developed against a unit test, etc.). 

 

Thank you,
Vasiliy.

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at http://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

[Oliver Beckstein]

Am Aug 7, 2014 um 22:09 schrieb Tyler Reddy <tre...@dal.ca>:

> Is this substantially different from the S^2 order parameter (for bond vectors) frequently calculated using NMR spectroscopy? I know plenty of people have written code to calculate that, although I don't think there is anything built into MDAnalysis that does it at the moment.

The calculation of lipid order parameters was described in the MDAnalysis J Comp Chem paper and there is also an example script in the distribution.

Oliver

[Vasiliy Triandafilidi]

Thank you for your quick response,

Here is my first programm that gets the z-axis. The main part that bothers me in my code is how to get the bond information effectively ( i.e how to get the bond vectors, because Mdanalysis easily gives me the bond length but is reluctant to give my the bond vector)
Even for a system of 500 atoms my programm doesn't work that fast. I plan to work with 100 k atoms, and its seems like i will be stack

#PROGRAM 1 
# get the z-axis
#get the z-axis vector (along which everything is aligned)
#this is being done by tensor Qab = 1/2 < 3(r_i*r_i) - 1> where r_i*r_i is the diadic product ( np.outer )  ,averaged by all atoms and some time steps
#then I get the eigen vector corresponding to the largest eigenvalue - this will be the z axis

def qab(coor):
  bb_res = coor
  nq = bb_res.shape[0]
  out = np.zeros((3,3))
  dij = np.eye(3)

  for i in xrange(0,nq):
    out += np.outer(bb_res[i],bb_res[i])
  return (1.5*out ) / float(nq) - 0.5*dij



u = Universe("polymer.psf", "trajectory_newparam.xtc")
bb = u.selectAtoms("all")

# get the bond list 

count = 0.
Qab = np.zeros((3,3))
for ts in u.trajectory[200:-1]:
	print "the frame is = ",ts.frame
	bonds_list = []
	for j in xrange(0,501):
		b = u.bonds[j]
		a1 = b.atom1.pos  		# <- this looks ugly 
		a2 = b.atom2.pos
		b1 = a2-a1
		bonds_list.append(b1/np.linalg.norm(b1))
	bonds_np = np.array(bonds_list)
	Qab += qab(bonds_np)
	count += 1.
print Qab
print np.trace(Qab)
vals, vecs =  np.linalg.eig(Qab)
print vals[0], vecs[:,0]
Qab = Qab/count

# the eigen vector is
# # [-0.62762023  0.4279941  -0.6503183 ]

Here is my second programm which basically analyses how aligned bonds are with respect to the z_axis

from MDAnalysis import *
import numpy as np
from numpy import linalg as LA
from pylab import *

#PROGRAM 2
# this program uses the z_axis direction
# and then calculates the P2 parameter = 3<cos*theta - 1> averaged over all atoms
# cos theta = (r_i r_z)/ ( |r_i| |r_z| ) 
# the angle between the bond and the z_axis


def P2(coor, principal_vector):
	nq = coor.shape[0]
	p2 = 0.
	pv2 = np.linalg.norm(principal_vector)
	for i in xrange(0,nq):
		v1 = np.linalg.norm(coor[i])
		some = np.dot(coor[i], principal_vector)/(v1*pv2)
		p2 += (3.*some*some - 1.)/2.
	p2 = p2 / float(nq-1)
	return p2

u = Universe("polymer_initial.psf", "trajectory_relaxation.dcd")
bb = u.selectAtoms("all")


#here is the order vector obtained from the previous programm
vec = np.array([ 0.13781049,  0.93179896,  0.33579603])
f = open('workfile', 'w')


Order_array = []

for ts in u.trajectory[1:-1]:
	bonds_list = []
	for j in xrange(0,501):
		b = u.bonds[j]
		a1 = b.atom1.pos
		a2 = b.atom2.pos
		b1 = a2-a1
		bonds_list.append(b1/np.linalg.norm(b1))
	# print b1
	bonds_np = np.array(bonds_list)
	order = P2(bonds_np,vec)
	print "frame = ", ts.frame, "order = ", order
	f.write(str(ts.frame) + " " + str(order) + "\n")
	Order_array.append((ts.frame, order))
f.close()

Order_array = np.array(Order_array)


plot(Order_array[:,0], Order_array[:,1], 'r--', lw=2, label=r"$P_2 = \dfrac{3 cos^2 \theta - 1}{2}$")
xlabel("time (ps)")
ylabel(r"order parameter $P_2$ ($\AA$)")

savefig('order_param.pdf')

On Thursday, August 7, 2014 10:10:36 PM UTC-7, Tyler Reddy wrote:

Are there are any written programms of the orientational order for polymers?

Is this substantially different from the S^2 order parameter (for bond vectors) frequently calculated using NMR spectroscopy? I know plenty of people have written code to calculate that, although I don't think there is anything built into MDAnalysis that does it at the moment.

I am new to MDanalysis and I wrote my code but it looks ugly and I was wondering if anyone can share their experience.

It might be useful to share the pertinent snippets of your code and provide a reference for the definition of the orientational order parameter (if it is different from the NMR version). If you'd like it to be built in to MDAnalysis you might want to file an Issue that describes these things, and includes an example file / data of some sort and the value of the order parameter (so it can be developed against a unit test, etc.). 

 

Thank you,
Vasiliy.

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discussion+unsub...@googlegroups.com.

[Richard Gowers]

I think the way I'd get the bond vectors would be....

from MDAnalysis.core.AtomGroup import AtomGroup # import the AtomGroup class so we can make some

bondlist = u.bonds[:501] # Make a list of all the bonds you want to work on

# Make two atomgroups, each containing one atom from each bond
ag1 = AtomGroup([b.atom1 for b in bondlist])
ag2 = AtomGroup([b.atom2 for b in bondlist])

# ag1.positions returns an array of all the coordinates
vecs = ag1.positions - ag2.positions # so this gives the vectors between bonds as np array
vecs /= np.sqrt((vecs*vecs).sum(axis=1)) # normalise to unit vectors

It should be a bit easier to get though, maybe I could add a .vectors() method to a TopologyGroup of bonds.

If you're working on a whole trajectory, you could allocate the results array before entering the for loop, so something like

Order_array = np.zeros(len(u.trajectory))

for i, ts in u.trajectory:
    Order_array[i] = something

[Vasiliy Triandafilidi]

I tried to implement what you have written but it doesn't seem to work.

ag1 = AtomGroup([b.atom1 for b in bondlist])

this line doesn't work as written
So I try to localize the mistake
bondlist = u.bonds
bondlist
>[(0, 1),
......
 (499, 500),
 (500, 501)]

then I try to do
b = bondlist[1]
print b.atom1  

b.atom1
Traceback (most recent call last):

>  File "<ipython-input-91-6626cfa2f882>", line 1, in <module>
    b.atom1

AttributeError: 'tuple' object has no attribute 'atom1'

so I changed it to

print b[0] -  this line works

But when I write

 ag1 = AtomGroup([b[0] for b in bondlist])
Traceback (most recent call last):

  File "<ipython-input-92-2868256205e7>", line 1, in <module>
    ag1 = AtomGroup([b[0] for b in bondlist])

NameError: name 'AtomGroup' is not defined

How can one define an Atomgroup using the the selections by atom numbers
(something like u.selectAtoms("bynum 1:10") where instead of 1:10 will be a custom array)

Thank you in advance,
Vasiliy.

[Richard Gowers]

I'm a bit confused, I thought u.bonds was always a list of Bond instances, looks like you had that before with your previous code.

If you're working just with the atom numbers, then you can just slice u.atoms ie

u.atoms[[0, 1, 3]]

returns an AtomGroup with atoms 0, 1 & 3

So for what you're doing..

ag1 = u.atoms[[b[0] for b in bondlist]]
ag2 = u.atoms[[b[1] for b in bondlist]]

should do the trick.

[Vasiliy Triandafilidi]

Thanks to everyone who helped, especially to Richard Gowers.
Here is what I managed to get. The code is working and I benchmarked the order parameter code  with the gyration radius evolution and they matched.
The only thing one needs to change is the Universe - just plug in your files and you are good to go. Please see the attached files for the output.

# Program: bond_order.py
# Purpose: get the P2(cos\theta) parameter using Python
# Author:  Triandafilidi Vasiliy , MSc student at CHBE UBC, Vancouver
# Syntax:  python bond_order.py
# Requires: polymer.psf, trajectory.xtc 

from MDAnalysis import *
import numpy as np
from numpy import linalg as LA
from pylab import *

def get_qab(coor):
    #gets the Qab for all atoms and averages it
    nq = coor.shape[0]
    out = np.zeros((3,3))
    # print "nq = ", nq

    out = np.zeros((3,3))
    dij = np.eye(3)

    for i in xrange(0,nq):

        out += np.outer(coor[i],coor[i])

    return (1.5*out ) / float(nq) - 0.5*dij

def get_bond_list(bondlist):
    #gets the universe creates a normalized bonds list

    # Make two atomgroups, each containing one atom from each bond

    ag1 = u.atoms[[b[0] for b in bondlist]]
    ag2 = u.atoms[[b[1] for b in bondlist]]

    vecs = ag1.positions - ag2.positions # so this gives the vectors between bonds as np array

    norm = np.sqrt((vecs*vecs).sum(axis=1)) # normalise to unit vectors
    vecs /= norm[:,None] #the norm vector is a (nx1) and we have to create dummy directions -> (n,3)
    return vecs


def get_eigen_vector(u):
    #gets the universe and loops over the trajector
    print "getting the eigen_values and eigen_vectors"
    Qab = np.zeros((3,3))
    count =  0.
    for ts in u.trajectory[300:-1]:
        # print "frame is = ", ts.frame
        bondlist = u.bonds
        Qab += get_qab(get_bond_list(bondlist))
        count += 1.
    Qab = Qab/count
    

print Qab
    print np.trace(Qab)
    vals, vecs =  np.linalg.eig(Qab)
    print vals[0], vecs[:,0]

    return vecs[:,0]

def get_p2(coor, principal_vector):

    nq = coor.shape[0]
    p2 = 0.
    pv2 = np.linalg.norm(principal_vector)
    for i in xrange(0,nq):
        v1 = np.linalg.norm(coor[i])
        some = np.dot(coor[i], principal_vector)/(v1*pv2)
        p2 += (3.*some*some - 1.)/2.

    p2 /= float(nq-1)
    return p2


def get_order_param(u):
    #this function gets the universe
    #uses the get_eigen_vector(universe) to get the eigen_vector
    #then uses the get_bond_list(bondlist) to get the vectorized coor 
    #then uses the get_p2(coor, principal_vector) to get the p2
    print "running the get_order_param"
    principal_vector = get_eigen_vector(u)
    print "the principal vector is = ", principal_vector
    Order_array = np.zeros((len(u.trajectory),2))
    # Order_array = []

    for ts in u.trajectory[1:-1]:

        bondlist = u.bonds
        vecs = get_bond_list(bondlist)
        order = get_p2(vecs,principal_vector)
        # Order_array.append((ts.frame, order))
        Order_array[ts.frame,0] = ts.frame
        Order_array[ts.frame,1] = order
        print "frame is = ", ts.frame, "order = " , order 
    return Order_array


u = Universe("polymer.psf", "trajectory.xtc")
Order_array = get_order_param(u)

plot(Order_array[:,0], Order_array[:,1], 'r--', lw=2, label=r"$P_2 = \dfrac{3 cos^2 \theta - 1}{2}$")
xlabel("time (ps)")
ylabel(r"order parameter $P_2$ ($\AA$)")

savefig('order_param.pdf')

[Oliver Beckstein]

Hi Vasiliy,

Thank you for posting your solution.

We could add it as an example script in the MDAnalysis distribution if you like to contribute it (under examples – see https://code.google.com/p/mdanalysis/source/browse/package#package%2Fexamples ). All you'd need to do is to

- add a line to the comments licensing it under the GPL v2, e.g.

# Copyright (c) 2014 Vasiliy Triandafilidi
# Released under the GNU Public Licence, v2 or any higher version

- add a reference to the calculation of the order parameters and/or a short description including the actual equations used (you can use LaTeX typesetting if you're familiar with it)

We have also been thinking about a specific MDAnalysis.analysis.polymers module and if someone gets round to creating such a module then your code could also be incorporated there.

Cheers,
Oliver

On 13 Aug, 2014, at 12:16, Vasiliy Triandafilidi wrote:

> Thanks to everyone who helped, especially to Richard Gowers.
> Here is what I managed to get. The code is working and I benchmarked the order parameter code with the gyration radius evolution and they matched.
> The only thing one needs to change is the Universe - just plug in your files and you are good to go. Please see the attached files for the output.
>
> # Program: bond_order.py
> # Purpose: get the P2(cos\theta) parameter using Python
> # Author: Triandafilidi Vasiliy , MSc student at CHBE UBC, Vancouver
> # Syntax: python bond_order.py
> # Requires: polymer.psf, trajectory.xtc

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst * orbe...@gmail.com