Project

General

Profile

Bug #49

Wrong image convention for evaluating dummy atoms.

Added by Ramon Garcia almost 14 years ago. Updated over 13 years ago.

Status:
Closed
Priority:
High
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

After running mdrun_d with a water system with van der Eerden model, mdrun_d
crashed. (Use attached files to reproduce the problem). After debugging, I
discovered that the reasons was that the coordinates of a dummy atom were wrong.
When they were evaluated, the positions of the atoms have been shifted twice
from the original ones. As the coordinates of dummy atoms are a linear
combination of other atoms, if the image conventions of them are inconsistent,
the dummy can have very wrong values. This caused a crash in the calculation of
forces.

Looking at the source code, it looks like the logic of shifting before virtual
site evaluation is wrong

if (bVsites) {
if (graph) {
/* Following is necessary because the graph may get out of sync * with the coordinates if we only have every N'th coordinate set
*/
if (bRerunMD || bExchanged)
mk_mshift(log,graph,state->box,state->x);
shift_self(graph,state->box,state->x);
}
construct_vsites(log,state->x,&mynrnb,inputrec->delta_t,state->v,
&top->idef,graph,cr,fr->ePBC,state->box,vsitecomm);
if (graph)
unshift_self(graph,state->box,state->x);
}

This code is written as if the particles were not shifted before, and so they
are shifted and then unshifted to restore the previous state. But previously
there has been a call to do_pbc_first and therefore the positions have been
already shifted. I changed the code above to:

if (bVsites) {
if (graph) {
/* Following is necessary because the graph may get out of sync * with the coordinates if we only have every N'th coordinate set
*/
if (bRerunMD || bExchanged) {
unshift_self(graph, state->box, state->x);
mk_mshift(log,graph,state->box,state->x);
shift_self(graph,state->box,state->x);
}
}
construct_vsites(log,state->x,&mynrnb,inputrec->delta_t,state->v,
&top->idef,graph,cr,fr->ePBC,state->box,vsitecomm);
}

which seems to be correct for me. Please review the proposed change. I have only
tested with the system I am working with.

sim.tgz (83.9 KB) sim.tgz Files used for reproducing the problem. Ramon Garcia, 02/18/2006 05:52 PM
md.c.diff (937 Bytes) md.c.diff Proposed fix Ramon Garcia, 02/18/2006 05:54 PM
after_em.pdb (141 KB) after_em.pdb Corrected coordinate file David van der Spoel, 02/24/2006 10:36 PM
md.log (13.7 KB) md.log Log file generated by a simulation with LINCS, otherwise with the original files posted. Ramon Garcia, 02/25/2006 05:58 PM

History

#1 Updated by Ramon Garcia almost 14 years ago

Created an attachment (id=15)
Files used for reproducing the problem.

#2 Updated by Ramon Garcia almost 14 years ago

Created an attachment (id=16)
Proposed fix

Go to the directory where md.c is, and do patch < md.c.diff

#3 Updated by David van der Spoel over 13 years ago

Hi Ramon,

tahnk you for a detailed bug report with proposed solution. I have reproduced
the problem, but I wonder whether it is triggered by broken molecules. i.e. the
starting molecules not being continuous. Gromacs can not handle this... This is
a bug in itself. If you have a coordinate file that is correct it might not
trigger this bug. I will investigate it further.

#4 Updated by David van der Spoel over 13 years ago

I have implemented your fix but it still crashes with a LINCS error.

#5 Updated by David van der Spoel over 13 years ago

Created an attachment (id=18)
Corrected coordinate file

Try running with this input file. It works fine for me. It was created by doing
an energy minimization of the system using flexible bonds.

#6 Updated by Ramon Garcia over 13 years ago

Sure the problem happens with "broken molecules", that is, different atoms of
the same molecule having different image conventions. That happens because the
coordinates come from DL_POLY.

I am not sure why Gromacs does not handle broken molecules. Gromacs has already
all the logic to do it correctly in mshift.c (I had to go through it while
debugging this problem).

For me Gromacs ran fine after the fix (I am using SHAKE rather than LINCS). And
I am sure that broken molecules where corrected by Gromacs (except that before
the fix, the shift was repeated and so the molecules where broken again).

I will see why it breaks with LINCS.

The configuration is already relaxed. One should not need any kind of minimization.

#7 Updated by Ramon Garcia over 13 years ago

I can't reproduce your LINCS warnings.

I am running exactly the same files as I uploaded, but with the constraint
algorith set to LINCS, and have no trouble. I killed Gromacs so that the output
was flushed. Constraint deviations are not particulary large.

I am running gromacs from the CVS tag release-3-3-patches, and with the
additional change described in this bug report. It is compiled with double
precision.

#8 Updated by Ramon Garcia over 13 years ago

Created an attachment (id=19)
Log file generated by a simulation with LINCS, otherwise with the original
files posted.

#9 Updated by Ramon Garcia over 13 years ago

Ping.

van der Spoel? are you there?

#10 Updated by Ramon Garcia over 13 years ago

ping again...

Please can you answer?

#11 Updated by David van der Spoel over 13 years ago

I have discussed this with Berk Hess and the concensus is that the shifting code
is fine in 3.3. The problem is that the broken molecules are not handled
correctly. Once this is fixed the 3.3 code should work without your patch.

#12 Updated by Ramon Garcia over 13 years ago

I am not sure if I understand. So you are saying that it was a chance that my
system worked after my patch and that there are other reasons why broken
molecules are not well supported? I don't have full knowledge of Gromacs source
code. All I can say is that my system didn't work before and after my patch it
worked fine.

If there are other deep reasons why broken molecules are not supported, I
suggest to test it (perhaps in grompp) and to warn the user. If the user knows
the problem he can workaround it, but with no more information than a segfault
he will be totally lost.

#13 Updated by David van der Spoel over 13 years ago

What needs to be done is to fix broken molecules once at the start of mdrun and
then use the code as is in 3.3.

#14 Updated by Erik Lindahl over 13 years ago

Fixed. We now call the routine rm_pbc() to make all molecules in mdrunner() before actually calling any of
the main routines.

Also available in: Atom PDF