[gmx-users] Question about drift in positions of frozen atoms

[Andrew DeYoung]

Hi, 

I am using freeze groups to hold the x, y, and z positions of specified
atoms fixed (to roughly simulate a fixed surface).  I would like to keep the
positions of atoms of residuetypes "AT1" and "AT2" fixed.  So, in my .mdp
file, I have the following keywords pertaining to "freeze":

freezegrps = AT1 AT2
freezedim = Y Y Y Y Y Y ; freeze x, y & z directions

When I run the simulation, I observe that the positions of the "frozen"
atoms change significantly.  For example, suppose I consider the x-position
of atom 2 (which is of residuetype "AT1" and is thus frozen).  In my input
file, conf.g96, I have the following entry for atom 2:

    1 AT1   CP         2    0.353750000    0.061271297   11.700000000

So, the x-position of atom 2 is 0.353750000 nm.

Then I run my simulation (in double precision) using grompp_d followed by
mdrun_d.  I save configurations every timestep (my timestep is 1 fs), using
nstxout = 1, and I use the NVT ensemble (using tcoupl = nose-hoover, tc_grps
= system, ref_t = 350, and tau_t = 1.0).  

Then, I use "g_traj_d -f my.trr -s my.tpr -n my.tpr -ox my.xvg -fp -noy
-noz" to extract only the x-position of atom2 and print these values to an
.xvg file.  I obtain the following:

 0       3.56560087860643e-01
 0.001   3.57711670093127e-01
 0.002   3.58265131434284e-01
...
 0.011   3.60305324395584e-01
 0.012   3.60525009903859e-01
...
 0.021   3.63643322576657e-01
 0.022   3.64169229293184e-01
...
 0.041   3.76755705048180e-01
 0.042   3.77036248430013e-01
...
 0.05    3.78067513842026e-01
 0.051   3.78112854943937e-01
...
 0.101   3.78302201368430e-01
 0.102   3.78302202249019e-01
...
 0.501   3.78302205990565e-01
 0.502   3.78302205990565e-01
...
 0.71    3.78302205990565e-01
 0.711   3.78302205990565e-01
...
 0.911   3.78302205990565e-01
 0.912   3.78302205990565e-01
...

and so on over time (the value converges to 3.78302205990565e-01 nm).  I
have two questions, if you have time:

(i) Why is the x-position of atom2 at time t = 0 ps "3.56560087860643e-01
nm" rather than something slightly closer to the value given in conf.g96
("0.353750000 nm")?  However, this discrepancy at time t = 0 ps is not all
that bad; the discrepancy is about 0.0028 nm.

(ii) Why does the x-position of atom2 drift drastically over time, until it
appears to converge to "3.78302205990565e-01 nm"?  This discrepancy seems
very significant (a drift of about 0.0245 nm from the initial x-position).
A drift of ~0.02 nm for my "frozen atoms," since the bond length between
adjacent frozen atoms is 0.1415 nm.  Why does this drift happen, and why
does it converge to a single (incorrect) value, even when I run several
nanoseconds?

Do you have any thoughts about why my atoms are not remaining frozen?  Is
there anything else I can try?  I checked the y- and z-components other
atoms, and their drifts are of the same order of magnitude.

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University