Hi,
For some reason your message came in all in a single line (i.e. there
were no newlines) which made reading your code a little difficult, but
I've spotted some problems:
On Thu, Jan 5, 2017 at 7:39 AM, wei <lwstudy.sina.com> wrote:
> Dear all,
> I want to verify that the translational diffusion constant I calculated for a sugar molecule is correct. I found that the magnitudes of mean square displacement (MSD) for one atom I computed using Ptraj/Cpptraj is fairly different from what I obtained from a small code I wrote (I just read in coordinates for all atoms directly then calculate MSDs as follows:
> open(11,file='image.crd')
> do
> read(11,*)
> read(11,'(10f8.3)') (x(i),y(i),z(i),i=1,natom)
> nframe=nframe+1
> cx(nframe)=x(1)
> cy(nframe)=y(1)
> cz(nframe)=z(1)
> enddo
> close(11)
This code is incorrect for reading Amber ASCII trajectories. The title
line only appears once (see
http://ambermd.org/formats.html#trajectory). You need to move the
'read(11,*)' statement to outside the do loop. Also, if your
trajectory has any box coordinates they need to be read after the
coordinates. A good check to see if your code is doing the right thing
is to try and write the data back out and see if you can visualize it.
> do t=0,2000
> do i=1,nframe-t
> dvx(t+1)=dvx(t+1)+(cx(i+t)-cx(i))**2
> dvy(t+1)=dvy(t+1)+(cy(i+t)-cy(i))**2
> dvz(t+1)=dvz(t+1)+(cz(i+t)-cz(i))**2
> enddo
> dvx(t+1)=dvx(t+1)/real(nframe-t)
> dvy(t+1)=dvy(t+1)/real(nframe-t)
> dvz(t+1)=dvz(t+1)/real(nframe-t)
>enddo).
This isn't the right calculation for mean-squared displacement. You
want the average of the delta of each from from the *initial* frame
because you want to know how far a particle moved (i.e. was displaced)
from it's initial position, not how far it is moving each frame on
average (see also
https://en.wikipedia.org/wiki/Mean_squared_displacement). It looks
like you're doing some kind of correlation calculation above.
> Ptraj input was wrote as:
> diffusion :1.C1 2 average output
> Then I got (MSD for a specific t):Mine code (239.17)Ptraj (936.91)
> Thus I am reading the code related to 'diffusion' command in Ptraj (please find it in AMBER14, 'actions.c'). I feel confused while seeing these lines of comments:
> /* * if the particle moved more than half the box, assume * it was imaged and adjust the distance of the total * movement with respect to the original frame... */ if ( state->box[0] > 0.0 ) { if ( delx > state->box[0]/2.0 ) diffusionInfo->deltax[currentAtom] -= state->box[0]; else if ( delx < -state->box[0]/2.0 ) diffusionInfo->deltax[currentAtom] += state->box[0]; if ( dely > state->box[1]/2.0 ) diffusionInfo->deltay[currentAtom] -= state->box[1]; else if ( dely < -state->box[1]/2.0 ) diffusionInfo->deltay[currentAtom] += state->box[1]; if ( delz > state->box[2]/2.0 ) diffusionInfo->deltaz[currentAtom] -= state->box[2]; else if ( delz < -state->box[2]/2.0 ) diffusionInfo->deltaz[currentAtom] += state->box[2]; }
This is the code in ptraj that attempts to account for imaging in
orthorhombic cells. As I wrote previously, if you don't 'unwrap' an
imaged trajectory prior to the diffusion calculation, you have to
correct for imaging "jumps", otherwise you get the wrong value for
diffusion. Cpptraj does this for both orthorhombic and
non-orthorhombic cells.
Hope this helps,
-Dan
> I am just trying to find the reason why I did not compute the same MSDs as Ptraj. I have two simple question:1. How to understand this description--"adjust the distance of the total movement with respect to the original frame"? Does 'total movement' refer to external degrees of freedom of the simulated box?2. Ptraj images atoms for all frames or just for those frames whose particle moved more than half the box?
> Any help would be highly appreciated and thanks in advance!
> Regards,Liu Wei--------------------------------
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
--
-------------------------
Daniel R. Roe
Laboratory of Computational Biology
National Institutes of Health, NHLBI
5635 Fishers Ln, Rm T900
Rockville MD, 20852
https://www.lobos.nih.gov/lcb
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jan 05 2017 - 06:00:03 PST
