We are trying to reliably calculate a self-diffusion constant using the MSD and are having trouble getting consistent results. In our case we have a somewhat large periodic box of a bcc solid below its melting point with one introduced vacancy. We now have good energy conservation in a NVE ensemble, but even when we do long runs ~ 1 ps, the results for the slope of MSD versus time vary significantly with the placement of the vacancy (middle or corner) and system size when we use the compute msd program. (Our tests suggest that there is not a problem with the periodic boundary conditions.) A textbook by Mark Tuckerman suggests that one can calculate better MSD averages by using using the individual time steps to also average over different initial times. I guess it would be possible to do this with an appropriate dump statement followed by post processing. Would this be what you would recommend for reliably calculating the diffusion constant or perhaps there is an easier method? Thanks in advance for your advice.
Sincerely, Natalie Holzwarth
N. A. W. Holzwarth email: natalie@…3310…
Department of Physics web: http://www.wfu.edu/~natalie
Wake Forest University phone: 1-336-758-5510
Winston-Salem, NC 27109 USA office: Rm. 300 Olin Physical Lab
Seems to me that you might get “relatively” large fluctuations from the slope of the MSD for a solid if that’s what you are calculating, since the diffusion coefficient is nearly zero.
Comparing results for the VACF might give more insight for the short time dynamics.
A textbook by Mark Tuckerman suggests that one can calculate better MSD averages by using using the individual >time steps to also average over different initial times. I guess it would be possible to do this with an appropriate >dump statement followed by post processing. Would this be what you would recommend for reliably calculating the >diffusion constant or perhaps there is an easier method? Thanks in advance for your advice.
yes, that is a good way to get better statistics from a single long run, i.e. use multiple time windows.
We don’t try to do that within LAMMPS, b/c it would require saving (and communicating) many “time-0”
snapshots as the simulation is running. So post-processing is the way to go.
many, many moons ago, i did some systematic study on convergence and
finite size effects of properties like g(r), MSD and so on. there
isn't much new in it that would warrant publishing it as a paper, but
over the years, it has been useful to show people what kind of errors
they are dealing with. please check out: http://klein-group.icms.temple.edu/akohlmey/files/talk-trieste2004-water.pdf
this wasn't done with LAMMPS and used a special purpose
post-processing code that was tuned for accuracy and minimizing
statistical errors.
Thanks very much for your very helpful advice. On a related topic, we tried also to look into the velocity correlation method using the advice given in your fix vector method (note: the link to fix/vector is incorrectly set; the correct page appears to be http://lammps.sandia.gov/doc/fix_vector.html). We attempted to use the following commands:
compute 2 all vacf
fix 5 all vector 1 c_2[4]
variable diff equal dt*trap(f_5)
thermo_style custom step v_diff
However we received the error message in our log file:
compute 4 all vacf
fix 5 all vector 1 c_4[4]
ERROR: Invalid fix style (../modify.cpp:714)
We are using lammps-1Feb14. Perhaps the vector command is in a special package that we did not load?? This may not be relevant to our case, since the velocity correlation functions seem to be extremely noisy, but it would be good to understand how to compute the velocity correlation for future reference. Thanks again. Natalie
Thanks very much for your very helpful advice. On a related topic, we
tried also to look into the velocity correlation method using the advice
given in your fix vector method (note: the link to fix/vector is incorrectly
set; the correct page appears to be http://lammps.sandia.gov/doc/fix_vector.html). We attempted to use the
following commands:
compute 2 all vacf
fix 5 all vector 1 c_2[4]
variable diff equal dt*trap(f_5)
thermo_style custom step v_diff
However we received the error message in our log file:
compute 4 all vacf
fix 5 all vector 1 c_4[4]
ERROR: Invalid fix style (../modify.cpp:714)
We are using lammps-1Feb14. Perhaps the vector command is in a special
package that we did not load??
No. you have too old a version of LAMMPS. fix vector was added later.
This may not be relevant to our case, since
the velocity correlation functions seem to be extremely noisy, but it would
you mention that you do "long" runs of 1ps. that is not long at all.