Thanks Jason. I have used the quasi-harmonic facility with ptraj but I
wanted to compare that to what I would get from the Schlitter formula.
Lepka
On Wed, Jul 14, 2010 at 1:30 PM, Jason Swails <jason.swails.gmail.com>wrote:
> Hello,
>
> I believe this corresponds to the quasi-harmonic approximation to the
> entropy, right? I think if you do it correctly, ptraj will calculate the
> entropy for you. This is actually done in MMPBSA.py by setting entropy=1 in
> the input file. You don't have to use that script, but the code will show
> you how to do it.
>
> Good luck!
> Jason
>
> --
> Jason Swails
> Quantum Theory Project,
> University of Florida
> Ph.D. Graduate Student
> 352-392-4032
>
> On Jul 14, 2010, at 12:50 PM, Lekpa Duukori <duukori.gmail.com> wrote:
>
> > Hello,
> >
> > I am having some trouble with some analysis on mass weighted covariance
> > matrix from ptraj
> >
> > I am trying to get entropy estimates using the Schlitter formula (Chem
> Phys
> > Lett, 1993,215,617). I first do an rms fit to the first frame of my data
> > using the following ptraj commands ( I am using AmberTools 1.2)
> >
> > trajin model.mdcrd.nowaters
> > trajout model.mdcrd.fit
> > rms first mass out rms.out.fit
> >
> > Then do
> >
> > trajin model.mdcrd.fit
> > matrix mwcovar out mwcovar.out
> >
> > I then use the following python script to calculate the entropy
> >
> > #!/usr/bin/python
> > import numpy
> > sigma = numpy.loadtxt('mwcovar.out')
> > hbar = 1.054571628e-34 #Joules seconds
> > T = 300.0 # Kelvin
> > Kb = 1.3806504e-23 #Joules/Kelvin
> > Euler = numpy.exp(1)
> > constant = 10e-23 #Angstroms squared multiplied by g --> Kg (10^-3)
> >
> > unit_matrix = numpy.eye(len(sigma[0]),len(sigma[0]),dtype=float)
> >
> > matrix = unit_matrix + ((Kb*T*Euler*Euler)/(hbar*hbar))*sigma*constant
> > Entropy = 0.5*Kb*numpy.log(numpy.linalg.det(matrix))
> >
> > But I get an infinite determinant and a NaN entropy, the elements matrix
> are
> > very large (10^26 or so, see below)
> >
> > I have done some dimensional analysis and think my units should be ok,
> but
> > the values in the matrix indicates otherwise so I don't know exactly what
> is
> > going wrong here. Any comments welcome.
> >
> > Thanks.
> >
> > Lekpa
> >
> > mwcovar matrix
> > [[ 0.75 -0.541 0.269 ..., 0.232 -0.571 0.065]
> > [-0.541 0.894 -0.4 ..., -0.597 0.793 -0.113]
> > [ 0.269 -0.4 1.099 ..., 0.354 -0.537 0.167]
> > ...,
> > [ 0.232 -0.597 0.354 ..., 1.102 -0.755 -0.031]
> > [-0.571 0.793 -0.537 ..., -0.755 1.762 -0.103]
> > [ 0.065 -0.113 0.167 ..., -0.031 -0.103 0.491]]
> >
> >
> > Output of python script
> >
> > [[ 2.06396838e+26 -1.48880919e+26 7.40276658e+25 ...,
> 6.38454218e+25
> > -1.57136792e+26 1.78877259e+25]
> > [ -1.48880919e+26 2.46025031e+26 -1.10078313e+26 ..., -1.64291883e+26
> > 2.18230256e+26 -3.10971236e+25]
> > [ 7.40276658e+25 -1.10078313e+26 3.02440166e+26 ..., 9.74193074e+25
> > -1.47780136e+26 4.59576959e+25]
> > ...,
> > [ 6.38454218e+25 -1.64291883e+26 9.74193074e+25 ..., 3.03265754e+26
> > -2.07772817e+26 -8.53106929e+24]
> > [ -1.57136792e+26 2.18230256e+26 -1.47780136e+26 ..., -2.07772817e+26
> > 4.84894971e+26 -2.83451657e+25]
> > [ 1.78877259e+25 -3.10971236e+25 4.59576959e+25 ..., -8.53106929e+24
> > -2.83451657e+25 1.35121130e+26]]
> > -inf
> > Entropy = nan
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jul 14 2010 - 14:30:03 PDT
