CHARMM Development Project Forums User Discussion & Questions Other User Discussion and Questions Different MSD from COOR ANALYSIS and CORREL

Hello, I have a question.

How can I match the time step (at x-axis of graph) ?
I did the MSD calculation but I could not match it.

My system information is:
time step: 2 fs
number of steps: 500,000
skip: 1000
Hence, 500 files are exist at a 1 dcd.
I wrote the files using a below script.

write r@n dumb time unit 10@n
* dummy
*

After I did calculations, I got below data:
If I did the MSD calculation with 1 ns, number of data are total 128 and
it starts from 0.000000 to 254.000000 (time step).

If I did the MSD calc. with 10 ns, number of data are total 2048 and
it starts from 0.000000 to 4.094000 (time step).

I read correl doc. and searched but I could not find why those values like that...
I want to know 2 things:

1. How those data show like that
2. How can I plot a graph related to time vs. msd (x vs. y plot).

Thanks a lot.

Your question about matching the timestep is unclear.
You also have not shown us what you actually did; this is required, since your description of the input file is less accurate than the input file itself.

This should have been a new post; now it is.

Due to lack of detail posted, the questions are difficult to impossible to answer. Please review the READ BEFORE POSTING topic for guidelines on how to post a question in these forums.

Dear Lennart and Rick,
Thanks for your answers.
I added details.
This script is for the lateral diffusion calc.

----------------------------------------------------------------------
! Define correlation function
correl maxt 40000 maxs 1024 maxa 40000 noupdate

! Define timeseries
set n 1
set nresid 128
label enterloop
enter X@n atom x sele segid MEMB .and. resid @n end
enter Y@n atom y sele segid MEMB .and. resid @n end
incr n by 1
if @n .le. @nresid goto enterloop

! Read trajectories
open unit 10 file name unfold_1-10.dcd
traj first 10 nunit 1 nocheck

! Write timeseries
set n 1
set nseries 128

label writeloop
corfun X@n X@n diff
enter s@n dupl corr
corfun Y@n Y@n diff
enter t@n dupl corr

enter r@n dupl s@n
mantime r@n add t@n

open unit 10@n write card name traj_10/msd_@n.dat

write r@n dumb time unit 10@n
* dummy
*

incr n by 1
if @n .le. @nseries goto writeloop

stop
--------------------------------------------------------------------

My approach to this would to obtain the center of mass time series for each lipid, set the z component to zero, and then use "CORFUN name name DIFF" on each of the modified time series.

I'm guessing from your first post that there may some issues with the time base of the DUMB TIME files you are writing. This can be dealt with by using the EDIT subcommand, via the DELTA, SKIP, and OFFSET keywords.

Dear Rick,

Thanks a lot for you kind answers.
As you mentioned, I tried all of them.
Here is my edited script (I modified a trajectory by setting the z values to 0):

--------------------------------------------------------------------------------------
! Define correlation function
correl maxt 40000 maxs 1024 maxa 40000 noupdate

! Define timeseries
set n 1
set nresid 128
label enterloop
enter M@n atom xyz sele mass segid MEMB .and. resid @n end
incr n by 1
if @n .le. @nresid goto enterloop

! Read trajectories
open unit 10 file name 1ns_2d.dcd
traj first 10 nunit 1 nocheck

! Write timeseries
set n 1
set nseries 128

label writeloop
corfun M@n M@n diff
enter r@n dupl corr

edit r@n veccod 1 delta 0.002 skip 1 offset 1

open unit 10@n write card name msd_@n.dat
write r@n unit 10@n dumb time
* dummy
*

increase n by 1
if @n .le. @nseries goto writeloop

stop
--------------------------------------------------------------------------------------

However, unfortunateIy, I could not see the exact accordance of time step.
As before, the data start from 1.000000 and the data end to 1.2540000.
Total number of data are 128 (from the output data) and there are 500 trajectories in a 1 dcd file (my trajectory file).

From the previous issue (http://www.charmm.org/ubbthreads/ubbthre...=true#Post34584), I calculated lateral MSD using the coor analy command. It's input is:

---------------------------------------------------------------------------------------
open read unit 101 file name 1ns_2d.dcd
traj query unit 101
traj iread 101 nread 1 skip 1000
calc nfile = ?NFILE
traj read

open unit 90 write form name 1ns_2d.msd
coor anal sele segid MEMB .and. resn DOPC end -
firstu 101 nunit 1 skip 1 -
imsd 90 ncors @nfile -
rspin 0.0 rspout 999.9 -
xbox @FFTX ybox @FFTY zbox @FFTZ

stop
-------------------------------------------------------------------------------------

From the coor anal, I got the exact data that show exact time step matching.
However, I am not sure that the MSD and D values are exact things which are showed by a output file (I doubt whether diffusion coefficient is right or not) . So I want to check it by using the corral command.

I am looking forward your answers.

Thanks a lot.

The time series starts at 1.0 because you used OFFSET 1; for an MSD time series, it should equal to the time spacing between frames.

The product of SKIP*DELTA should also be the time interval between frames; if that is 1.0 ps, then SKIP should be 500

I used a combination of EDIT and MANTIME commands in CORREL to set the z component of the center of mass time series.

I've said it before in these forums, and I'll say it again; COOR ANAL is only suitable for small rigid molecules such as solvents, as there is no facility for using the center of mass for the selected atoms. I tried it for carbohydrates, and it gave very different results from using the more rigorous approach via CORREL. It's even worse for lipids, which are known to diffuse on two different time scales; COOR ANAL uses a linear fit to the MSD data, which cannot be correct. Caveat emptor.

Yes, It is start from 1.00000 because of OFFSET 1.
In my case, the time step is 2.0 fs, so I set the skip to 1.
Alhough I changed the skip, the number of data is not changed.
It just shows only change of the position of dot (e.g., 1.126000 -> 126.000).

If you have another sight about this, please let me know.

In addition, In order to check the time step matching, I did short simulation and compared it which is in charmm test (c25test/cortst.inp).

It's input is:
---------------------------------------------------------------------------------------
dynamics VERLET strt timestep 0.002 nstep 100 -
iprfrq 100 ihtfrq 0 ntrfrq 0 ieqfrq 100 -
iunrea -1 iunwri -1 iuncrd 10 kunit -1 iunvel -1 -
nprint 50 nsavv -1 nsavc 1 inbfrq 20 ihbfrq 0 -
firstt 285.0 finalt 285.0 teminc 0.0 isvfrq 0 -
iasors 0 iasvel 1 iscvel 0 ichecw 0 twindh 5.0 twindl -5.0

if ?nocorrel .eq. 1 then stop

open unit 10 read file name @9cortst.trj

CORREL MAXT 5000 MAXS 10

ENTER DD atom x mass sele segid MAIN .and. resid 7 .and. type N end

TRAJ FIRSTU 10 NUNIT 1 BEGIN 1 STOP 100 SKIP 1

EDIT ALL DELTA 0.001

SHOW ALL

!mantime DD daverage
!write DD unit 100 dumb time

corfun DD DD diff

open write unit 100 card name msd.dat
write corr dumb time unit 100

END
-----------------------------------------------------------------------------------

--------------------------- It's output: -----------------------------------------

CORREL> SHOW ALL
NSER: 1
NAMES: DD
TOTALS: 100
AVERAGE: -0.974176
FLUCT: 0.165124
VECCOD: 1
CLASS: ATOM
VELCOD: 0
SKIP: 1
DELTA: 0.001000
OFFST: 0.002000
GECOD: 1
VALUE: 0.000000
Atom pointers for all time series:
SERPT: 1
SERNQ: 1
QAT : 29

CORREL>

CORREL> !mantime DD daverage
CORREL> !write DD unit 100 dumb time
CORREL>

CORREL> corfun DD DD diff
NSER: 1
NAMES: CORR
TOTALS: 32
AVERAGE: 0.040162
FLUCT: 0.031170
VECCOD: 1
CLASS: CORR
VELCOD: 0
SKIP: 1
DELTA: 0.001000
OFFST: 0.000000
GECOD: 1
VALUE: 0.000000
-----------------------------------------------------------------------------------

I checked that the TOTALS is changed from 100 to 32 after the CORFU DIFF command.
In the it's data file, there are total 32 data and in my opinion, the number of data is related to those TOTALS.
However, I can not understand why the value goes down from 100 to 32.
I also calculated this with turned off and truned on the EDIT command but there are no alteration at time step matching what I want to know.

Thank you again.

From the CORFUN section in correl.doc:

TOTAL integer

The TOTAL value determines the number of points to keep in
the correlation function. The number of points may not be
grater than the number of points in the time series. A reasonable
value is about 1/4 to 1/3 the length of the time series.
Correlation function values near the end have little weight.
The default value is the nearest power of two less than half of
the time series length.

OMG...
I'm sorry that I missed in searching the TOTALS at a correl.doc
As your mention, I calculated MSD using the TOALS option, as a result, I can get the exact result what I want to get.

Thanks a lot for your comment, and time.
Have a nice day.

Matthew Kim