Possible distance restraint bug in GMX 4.0.5
Created an attachment (id=394)
Files required to reproduce the possible bug
The application of distance restraint between water oxygen atom and a dummy site fixed at the center of a sphere work in GMX 3.3.3 version but not in 4.0.5.
Follow as attachment the .mdp, .gro, .top, .ndx and .itp files, the .tpr's generated with these files with 3.3.3 and 4.0.5 versions and the gmxdump output from these two .tpr's.
#1 Updated by David van der Spoel over 9 years ago
If I run your system esfera3.tpr in gromacs 3.3.3 it seems to work fine. If I run it in the latest git it also runs fine, but produces energies (Coulomb, Lennard Jones) that are not a number (nan). The dynamics does not seem to be affected though, so maybe this is a quirk with the new all-to-all routines?
If I run it in 4.0.5 everything seems fine?
Now if I run mdrun (4.0.5) with the esfera4.tpr it crashes in the distance restraint code, which it doesn't do in git head. It seems that your original problem has been resolved, and it is replaced by a new one.
Erik, any insights in all-to-all problems?
#2 Updated by Erik Lindahl about 9 years ago
I can confirm that this does indeed (1) seem to be caused by all-vs-all kernels, and (2) my fix last week doesn't solve this one.
First, you can easily disable the all-vs-all kernels by setting the environment variable GMX_NO_ALLVSALL (to any value, e.g. 1), but this should hopefully be fixed later today...
#3 Updated by Erik Lindahl about 9 years ago
This bug is caused by the usage of an extra 'atom' in all water molecules with zero charge and LJ parameters, but the atom has exactly the same coordinates for all water molecules.
The only reason this works with the neighborlist code is that we optimize the neighborlists to remove any interaction with zero parameters, which isn't possible in the all-vs-all kernels.
However, fundamentally it is an error to require evaluation of nonbonded interactions between pairs of atoms with distance=0. The inverse square root of that distance is NaN (not-a-number), and even if all parameters are zero, Nan*0 is an undefined number (i.e., still NaN).
Thus, this topology is incorrect, and it is only a coincidence it works in some cases. As a side effect, the present charge groups will also have huge radii. While that doesn't matter when the cutoffs are infinite, it's bad practice in general and will definitely fail with cutoffs.
The proper way to restrain water molecules to an arbitrary position in space is to add a position restraint, e.g. to the oxygen, and then provide a position restraint reference file where the oxygen atom has the reference coordinate to which you want to restrain it.