[gmx-users] rare huge forces in constraint simulation with direction_periodic

[Jochen Hub]

Hi,

I am running a constraint simulation, with 12 molecules restrained with 
respect to the center of mass of a lipid membrane. During 5ns of 
simulation, I get a few frames (e.g. 10 or so) with huge forces on all 
molecules, whereas the forces in all other frames are fine. When looking 
at the simulation, nothing looks spurious. There are no LINCS warnings 
or so. Here is an example. Note the large forces at t=385.32. In 
pullx.xvg, the molecules are not set to the correct constrained position 
in the problematic frame:

pullf (12 pull groups, 4 pull shown here):
    385.12   6.65942   325.524  -122.863  -156.645 ...
    385.16   312.699   498.147  -282.365   66.9014 ...
     385.2   271.748   581.749   42.3289  -48.3744 ...
    385.24   256.169   481.929  -147.006  -163.566 ...
    385.28   360.289   435.808  -128.435   126.654 ...
    385.32   27214.7   27461.9   26817.9   27025.9 ...
    385.36   140.386   129.482   359.494  -5.16442 ...
     385.4 -0.0698604  -457.767    -68.48    173.24 ...
    385.44  -323.406  -161.396  -187.093   396.581 ...
    385.48  -427.383   19.5112  -144.655  -95.9228 ...
    385.52  -102.406   61.0266  -758.817   35.7622 ...

pullx (12 pull groups, reference group + 4 pull groups shown here). The 
four molecules shown here are restrained at z=-1.831 with respect to the 
lipid membrane:
    385.12    3.5217    -1.831    -1.831    -1.831    -1.831 ...
    385.16    3.5217    -1.831    -1.831    -1.831    -1.831 ...
     385.2    3.5216    -1.831    -1.831    -1.831    -1.831 ...
    385.24   3.52142    -1.831    -1.831    -1.831    -1.831 ...
    385.28   3.52231    -1.831    -1.831    -1.831    -1.831 ...
    385.32    3.5222  -1.83208  -1.83208  -1.83208  -1.83208 ...
    385.36   3.52209    -1.831    -1.831    -1.831    -1.831 ...
     385.4   3.52192    -1.831    -1.831    -1.831    -1.831 ...
    385.44   3.52175    -1.831    -1.831    -1.831    -1.831 ...
    385.48   3.52165    -1.831    -1.831    -1.831    -1.831 ...
    385.52   3.52156    -1.831    -1.831    -1.831    -1.831 ...


I used gmx 4.5.5-dev, the relevant mdp options are (pressure coupling in 
z is switched off):
pull                     = constraint
pull_geometry            = direction_periodic
pull_dim                 = N N Y
pull_constr_tol          = 1e-06
pull_ngroups             = 12
pull_group0              = refgroup
pull_weights0            =
pull_pbcatom0            = 54
pull_group1              = Solute_Layer0_Channel0
pull_weights1            =
pull_pbcatom1            = 0
pull_vec1                = 0.0 0.0 1.0
pull_init1               = -1.8310
pull_rate1               = 0
...

Apart from this, I use the md integrator, LINCS for all bonds, V-rescale 
T-coupling (tau=5ps), pressure coupling in the xy plane only 
(compressibility = 4.5e-5 0 with Pcoupltype = semiisotropic), LJ cutoff 
at 1nm, PME, so normal settings.

Did anyone see a similar behavior? Or any hints what one could try? Bug 
or feature of the constraint algorithm?

Many thanks in advance,
Jochen



-- 
---------------------------------------------------
Dr. Jochen Hub
Computational Molecular Biophysics Group
Institute for Microbiology and Genetics
Georg-August-University of Göttingen
Justus-von-Liebig-Weg 11, 37077 Göttingen, Germany.
Phone: +49-551-39-14189
http://cmb.bio.uni-goettingen.de/
---------------------------------------------------