[gmx-users] Umbrella Sampling Pull code Problem

[Steinbrecher, Thomas (IPC)]

Dear Gromacs users,

I have encountered a strange problem when trying to set up umbrella sampling simulations using gromacs 4.5.3.

My system contains two groups with a COM distance of 1.95 nm (all distances measured by g_dist). 

Trying to set up the first US window, I use the following pull code, together with typical mdp parameters from the Gromacs US tutorial:

; pull options
pull            = umbrella
pull_geometry   = distance
pull_dim        = Y Y Y
pull_start      = yes
pull_init1      = 0.00
pull_nstxout    = 1
pull_nstfout    = 0
pull_group0     = Protein
pull_group1     = OHM
pull_rate1      = 0.0
pull_k1         = 10000


this results in a simulation in which the group distance remains close to 1.95nm, as expected.

However, when I make the following two changes in my input file:

pull_start      = no
pull_init1      = 1.95

which should (?) amount to an equivalent setup, a very different trajectory results in which the COM distance quickly increases to 2.7 nm and then appears to be restrained there. (Visualization confirms, in the first case, the groups remain in their starting conformation, in the second one, they are pushed appart)

Interestingly, the grompp output contains the following lines:

Pull group  natoms  pbc atom  distance at start     reference at t=0
       0       994       497 
       1         2     46347   1.224                 1.950              


in this case. Apparently, grompp (and mdrun thereafter) calculates the group COM distance differently from g_dist! I think this is not a PBC issue, every atomic distance within both groups is smaller than half the box size and the pbcatoms are close together. However, when I set pbcatom0 to various atom numbers, different 'distance at start' values are obtained, but never the correct COM distance. The two groups do not have overlapping atoms. I am sure I used the same group indexes for pulling and distance measurements.

This behaviour is so visibly wrong that I cannot believe this is a bug and rather think I am doing something incorrect. A search of the list revealed a somewhat similar report by Gavin Melaugh in 2011 which did not resolve the issue.

Any ideas on what might be the problem here?

I am willing to send around files to reproduce the problem of course.

Thomas

Dr. Thomas Steinbrecher
Institut für Physikalische Chemie, KIT
Kaiserstr. 12, 76131 Karlsruhe