- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Daniel Roe <daniel.r.roe.gmail.com>

Date: Thu, 5 Jan 2017 08:52:17 -0500

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
*

Date: Thu, 5 Jan 2017 08:52:17 -0500

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:

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.

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.

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

-- ------------------------- 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/amberReceived on Thu Jan 05 2017 - 06:00:03 PST

Custom Search