http://redmine.gromacs.org/http://redmine.gromacs.org/favicon.ico?14024030812014-11-22T13:41:58ZGROMACS developmentGROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=100982014-11-22T13:41:58ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul><li><strong>File</strong> <a href="/attachments/1326/dimer.tgz">dimer.tgz</a> <a class="icon-only icon-download" title="Download" href="/attachments/download/1326/dimer.tgz">dimer.tgz</a> added</li></ul><p>Attached example shows that the Ecoul recip is the culprit. A water dimer is simulated for 0 steps. If the dimer is located inside the box group and verlet yield the same result, if one of the molecules is on the edge of the box the Ecoul recip gets a strange very high value after two iterations (note that the patch above takes care of differences due to the first do_force call).</p>
<p>Group:<br />LJ (SR) 0 -- 0 0 (kJ/mol)<br />Coulomb (SR) 0 -- 0 0 (kJ/mol)<br />Coul. recip. 0.0419108 -- 0 0 (kJ/mol)<br />Polarization 5.97909e-05 -- 0 0 (kJ/mol)<br />Potential 0.0419706 -- 0 0 (kJ/mol)<br />Kinetic En. 0.00202031 -- 0 0 (kJ/mol)</p>
<p>Verlet:<br />LJ (SR) 0 -- 0 0 (kJ/mol)<br />Coulomb (SR) <del>4.11256 -</del> 0 0 (kJ/mol)<br />Coul. recip. 110.869 -- 0 0 (kJ/mol)<br />Polarization 5.97914e-05 -- 0 0 (kJ/mol)<br />Potential 106.756 -- 0 0 (kJ/mol)<br />Kinetic En. 0.00202491 -- 0 0 (kJ/mol)</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=100992014-11-22T15:11:53ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul><li><strong>File</strong> <a href="/attachments/1327/dimer-nopol.tgz">dimer-nopol.tgz</a> <a class="icon-only icon-download" title="Download" href="/attachments/download/1327/dimer-nopol.tgz">dimer-nopol.tgz</a> added</li><li><strong>Subject</strong> changed from <i>Difference in energy for shell systems with Verlet scheme</i> to <i>Difference in energy systems with Verlet scheme</i></li><li><strong>Priority</strong> changed from <i>Normal</i> to <i>High</i></li></ul><p>The issue is not reserved to shell systems as we thought at first. Attached is an example without shells where the reciprocal energy is incorrect as well.</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101002014-11-22T15:19:09ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul><li><strong>Subject</strong> changed from <i>Difference in energy systems with Verlet scheme</i> to <i>Difference in energy for shell systems with Verlet scheme</i></li><li><strong>Priority</strong> changed from <i>High</i> to <i>Normal</i></li></ul><p>OK, sorry it happens only for shells anyway.</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101012014-11-22T17:06:13ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><p>It is somewhat more complicated. Due to some historic reason the mdp file has a surface correction for PME. This triggers an issue in the Verlet code, namely that the total system dipole needs to be computed. This only works correctly when molecules are intact. Calling put_atoms_in_box_omp breaks this. Removing the call to this function and having a good start structure fixes this for the first round.</p>
<p>In subsequent calls to do_force I can prevent that mu is computed again by changing the force_flags like:</p>
<pre><code>/* Try the new positions */<br /> do_force(fplog, cr, inputrec, 1, nrnb, wcycle,<br /> top, mtop, groups, state-&gt;box, pos[Try], &#38;state-&gt;hist,<br /> force[Try], force_vir,<br /> md, enerd, fcd, state-&gt;lambda, graph,<br /> fr, vsite, mu_tot, t, fp_field, NULL, bBornRadii,<br /> force_flags &#38; ~GMX_FORCE_NS &#38; ~GMX_FORCE_STATECHANGED);</code></pre>
<p>However, this means that the total system dipole is not updated during the minimization and that the total energy is not correct.</p>
<p>So how does one in the Verlet framework make molecules whole again?<br />If that is not possible the dipole will always be incorrect, right?</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101022014-11-23T08:35:22ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul><li><strong>File</strong> <a href="/attachments/1328/spc-surf.tgz">spc-surf.tgz</a> <a class="icon-only icon-download" title="Download" href="/attachments/download/1328/spc-surf.tgz">spc-surf.tgz</a> added</li><li><strong>Subject</strong> changed from <i>Difference in energy for shell systems with Verlet scheme</i> to <i>Difference in energy with Verlet scheme due to ewald_geometry = 3dc</i></li></ul><p>It seems there is a real issue with Verlet and PME. Since the molecules are broken the total system dipole is incorrect and hence the PME dipole correction can not be computed. The attached test system, spc-surf shows that for a system with ewald_geometry = 3dc<br />the dipole and energy of normal SPC water diverge quickly (within 20 steps in double precision). The dipole correction to the energy is way too high, e.g.<br />group/mdrun0.debug:dipcorr = 0.000 0.000 -22.673<br />verlet/mdrun0.debug:dipcorr = 0.000 0.000 -193.160</p>
<p>It should be noted that this is a pathetic case in some sense since there is no vacuum in the box and hence molecules are broken in the normal direction too which they would not normally be.<br />The dipoles in the X and Y direction as expected also different:<br />group/mdrun0.debug:mutot = -0.021 -0.495 -0.084<br />verlet/mdrun0.debug:mutot = -2.301 -1.260 -0.714</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101032014-11-23T08:46:38ZGerrit Code Review Botadmin@gromacs.org
<ul></ul><p>Gerrit received a related patchset '1' for Issue <a class="issue tracker-1 status-5 priority-4 priority-default closed" title="Bug: Difference in energy with Verlet scheme due to PME dipole correction (Closed)" href="http://redmine.gromacs.org/issues/1645">#1645</a>.<br />Uploader: David van der Spoel (<a class="email" href="mailto:davidvanderspoel@gmail.com">davidvanderspoel@gmail.com</a>)<br />Change-Id: I4c756e835f35834d474c6372ee942b2037c34e28<br />Gerrit URL: <a class="external" href="https://gerrit.gromacs.org/4232">https://gerrit.gromacs.org/4232</a></p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101042014-11-23T09:27:04ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul><li><strong>Subject</strong> changed from <i>Difference in energy with Verlet scheme due to ewald_geometry = 3dc</i> to <i>Difference in energy with Verlet scheme due to PME dipole correction</i></li></ul><p>Digging somewhat more I now found issues also with the default ewald_geometry = 3d and Verlet when the epsilon_surface is not 0 (infinity). E.g. when using 65 for SPC (which happens to be the magic SPC dielectric constant) one gets after some steps:<br />group/mdrun0.debug:Total dipole correction: Vdipole=0.294524<br />verlet/mdrun0.debug:Total dipole correction: Vdipole=11.3957</p>
<p>Since the dipole correction also influences the force this is a real problem, even though many would not use it. I guess we just have to turn off the dipole correction completely with Verlet.<br />We have:<br />group/mdrun0.debug:dipcorr = -0.061 -1.016 -0.423<br />verlet/mdrun0.debug:dipcorr = -6.345 -2.571 -0.424</p>
<p>and the code does the following:<br /> if (dipole_coeff != 0)
{<br /> for (j = 0; (j < DIM); j++)
{<br /> f[i][j] -= dipcorrA[j]*chargeA[i];<br /> }<br /> }<br />In other word large spurious forces will be added to the atoms.</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101052014-11-23T10:05:01ZBerk Hessgmx3@hotmail.com
<ul></ul><p>These are not Verlet scheme issues, but occur with non-neutral charge groups.</p>
<p>The 3dc geometry should only be used without pbc in z, we should at least add a warning for this. I also used it for water slab in the middle of the box with full pbc. Then you only have minor issues when a single water passes the border.</p>
<p>A non-zero surface epsilon could work with neutral molecules, but without charge groups we will need to correct for molecules split over pbc.</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101062014-11-23T11:02:09ZGerrit Code Review Botadmin@gromacs.org
<ul></ul><p>Gerrit received a related patchset '1' for Issue <a class="issue tracker-1 status-5 priority-4 priority-default closed" title="Bug: Difference in energy with Verlet scheme due to PME dipole correction (Closed)" href="http://redmine.gromacs.org/issues/1645">#1645</a>.<br />Uploader: David van der Spoel (<a class="email" href="mailto:davidvanderspoel@gmail.com">davidvanderspoel@gmail.com</a>)<br />Change-Id: I02e317cbddb47256f942312ec53c5bab2b13be2a<br />Gerrit URL: <a class="external" href="https://gerrit.gromacs.org/4233">https://gerrit.gromacs.org/4233</a></p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101082014-11-23T17:06:52ZBerk Hessgmx3@hotmail.com
<ul></ul><p>I would be nice to fix the actual issue. But to do this for the general case would require graph like code, which we already have, combined with domain decomposition. This is difficult and might be expensive. Doing it for molecules up to the size of the cut-off radius would not be too hard. We could add an "interaction"-type dipole to the DD reverse top which then automatically gets assigned to the domain where all coordinates of the molecule are available. The question is how often the surface_epsilon option is used (correctly).</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101092014-11-23T20:48:47ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><p>The bulk dipole correction drops off quickly with box size, therefore it can probably be neglected in most cases. The surface correction is more important but not a problem in practice either if the liquid is in the center of the box. So is it really worth the effort?</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=101232014-11-26T10:41:46ZBerk Hessgmx3@hotmail.com
<ul></ul><p>The dipole and surface epsilon correction are identical, except that surface epsilon works on all 3 dimensions.<br />The dipole correction is not really problematic, since molecules should not cross pbc in z.<br />The surface epsilon only makes sense when only molecules with no net charge are mobile and any charged molecules do not cross pbc. If there is a need for such cases, we could implement a correct treatment for the Verlet scheme.</p> GROMACS - Bug #1645: Difference in energy with Verlet scheme due to PME dipole correctionhttp://redmine.gromacs.org/issues/1645?journal_id=109942015-06-22T16:03:42ZErik Lindahlerik@kth.se
<ul><li><strong>Status</strong> changed from <i>New</i> to <i>Closed</i></li></ul><p>4233 has been merged and solved the main problem, and the remaining point appears to mostly be a theoretical issue. Closing for now, but if somebody working on shells wants to spend time on improving it we can reopen.</p>