It seems that cpptraj cannot do the following first and second time
autocorrelation functions:
<cos(phi(t)-phi(t+t0))>_t0
<(3*cos^2(phi(t)-phi(t+t0)) - 1)/2>_t0
using phi(t) from the output of dihedral command.
Thus I wrote a python script to do that, reading as input the output file
from dihedral command. Maybe someone else can profit from this.
#!/usr/bin/env python
import numpy
if __name__=='__main__':
import argparse
parser = argparse.ArgumentParser(description='calculate first and second
order correlations')
parser.add_argument('--infile',help='input file with dihedral angles')
parser.add_argument('--outfile',help='output file')
args = parser.parse_args()
phis=(numpy.pi/180.)*numpy.loadtxt(args.infile, comments='#',
usecols=(1,))
# First order Legendre Polynomial autocorrelation
<cos(phi(t)-phi(t+t0)>_t0
# We use cos(a-b)=cos(a)cos(b)+sin(a)sin(b)
cosphis=numpy.cos(phis)
sinphis=numpy.sin(phis)
c1=numpy.correlate(cosphis,cosphis,mode='same')+numpy.correlate(sinphis,sinphis,mode='same')
c1/=max(c1)
# Second order Legendre Polynomial autocorrelation
<0.5*(3*cos^2(phi(t)-phi(t+t0))-1)>_t0
# We use the relation cos^2(x)=(1+cos(2x))/2 and
cos(a-b)=cos(a)cos(b)+sin(a)sin(b)
cos_2phis=numpy.cos(2*phis) # cos(2*phi)
sin_2phis=numpy.cos(2*phis) # sin(2*phi)
c2=numpy.correlate(cos_2phis,cos_2phis,mode='same')+numpy.correlate(sin_2phis,sin_2phis,mode='same')
#<cos(2phi(t)-2phi(t+t0)>_t0
c2=(1+c2)/2 # <cos^2(phi(t)-phi(t+t0)>_t0
c2=(3*c2-1)/2
c2/=max(c2)
#ouput
n=len(phis)
nh=n/2
dt=1.0 #picoseconds
buf='# time c1 c2\n'
for i in range(nh,n):
buf+='{0:4.0f} {1:7.4f} {2:7.4f}\n'.format( dt*(i-nh), c1[i], c2[i] )
open(args.outfile,'w').write(buf)
On Wed, Jul 16, 2014 at 7:40 PM, Jose Borreguero <borreguero.gmail.com>
wrote:
> Dear Amber users,
>
> I want to calculate the time autocorrelation function of a dihedral angle.
> I've used dihedral command from cpptraj to get the angle. For instance:
>
> dihedral theta .Si6 .C18 .C20 .C24 out theta.dat
>
> Command timecorr requires a vector as input, so I guess my question is
> how do I translate dataset theta into vector set (cos(theta), sin(theta),
> 0) so that I can feed it to timecorr.
>
> Related question: with output file theta.dat I can create file vector.dat
> with columns cos(theta), sin(theta), and 0. Then is it possible to read
> file vector.dat into a vector set that I can pass to timecorr?
>
> Best regards,
> -Jose Borreguero
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jul 17 2014 - 12:00:02 PDT
