Project

General

Profile

Bug #2942

Pressure not converging with semiisotropic scaling and wall restraints

Added by Marvin Bernhardt 14 days ago. Updated 8 days ago.

Status:
Fix uploaded
Priority:
Normal
Assignee:
Category:
documentation
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I am trying to keep a mixture of methanol and water in a slab, that is surrounded by water. The methanol has a virtual site at its COM, which is restrained with a flat-bottomed harmonic potential in z-direction.

I tried different barostats with isotropic and semiisotropic scaling. Strangely, the pressure is not correct, when I use semiisotropic box scaling (xy or z). This is consistent with different barostats (tested with Berendsen, Parrinello-Rahman).

I am not completely sure this is a bug, since the setup is a bit exotic. But It is very weird to me, that with semiisotropic scaling the pressure is constantly off and the box not scaled further, while it works with isotropic scaling.

Tested on
24 core, AMD Opteron(tm) Processor 6174
4 core, Intel(R) Core(TM) i5-6300U

grompp-iso.mdp (584 Bytes) grompp-iso.mdp Marvin Bernhardt, 05/07/2019 05:13 PM
grompp-xy.mdp (596 Bytes) grompp-xy.mdp Marvin Bernhardt, 05/07/2019 05:13 PM
conf.gro (5.88 MB) conf.gro Marvin Bernhardt, 05/07/2019 05:13 PM
ffbonded.itp (205 KB) ffbonded.itp Marvin Bernhardt, 05/07/2019 05:13 PM
ffnonbonded.itp (61.5 KB) ffnonbonded.itp Marvin Bernhardt, 05/07/2019 05:13 PM
forcefield.itp (862 Bytes) forcefield.itp Marvin Bernhardt, 05/07/2019 05:13 PM
MET.itp (962 Bytes) MET.itp Marvin Bernhardt, 05/07/2019 05:13 PM
SOL.itp (1.09 KB) SOL.itp Marvin Bernhardt, 05/07/2019 05:13 PM
topol.top (157 Bytes) topol.top Marvin Bernhardt, 05/07/2019 05:13 PM
restraint.gro (3.84 MB) restraint.gro Marvin Bernhardt, 05/07/2019 05:13 PM
conf.gro (34.5 KB) conf.gro Marvin Bernhardt, 05/08/2019 10:41 AM
grompp-iso.mdp (585 Bytes) grompp-iso.mdp Marvin Bernhardt, 05/08/2019 10:41 AM
grompp-xy.mdp (597 Bytes) grompp-xy.mdp Marvin Bernhardt, 05/08/2019 10:41 AM
grompp-z.mpd (595 Bytes) grompp-z.mpd Marvin Bernhardt, 05/08/2019 10:41 AM
restraint.gro (22.5 KB) restraint.gro Marvin Bernhardt, 05/08/2019 10:41 AM
topol.top (697 Bytes) topol.top Marvin Bernhardt, 05/08/2019 10:41 AM

Associated revisions

Revision f63083b6 (diff)
Added by Berk Hess 7 days ago

Correct user guide postion restraint virial

Corrected comment in the user guide note on virial and pressure
with position restraints and refcoord-scaling none or com.

Refs #2942

Change-Id: I3cdb5a01640c388f1c4222b1a43b342ac55fd882

History

#1 Updated by Marvin Bernhardt 13 days ago

I made a much simpler example with two types of LJ-particles, where the first is restrained into a slab.

With isotropic scaling pressure converges to 1000 bar.
With xy scaling the pressure converges to ~850 bar.
With z scaling the pressure converges to ~1200 bar.

This is not be expected behavior, is it?

#2 Updated by Berk Hess 13 days ago

Difficult to say. For your xy case, the pressure in x and y is 1000, as you requested. In z it is 500, whereas you requested 0, but I think it might be impossible to achieve 1000 in x/y and 0 in z in a LJ system.

#3 Updated by Marvin Bernhardt 13 days ago

Thanks for the answer. I did not have the idea to look at pressure components yet. Neither to set the target pressure in z direction to 1000 bar.

However, now that i tried it, things become a bit more clear.

Observations:
1. xy-scaling, no restraint, p_tgt_z = 0 -> Pres-ZZ = 1000
2. xy-scaling, no restraint, p_tgt_z = 1000 -> Pres-ZZ = 1000
3. xy-scaling, slab restraint, p_tgt_z = 0 -> Pres-ZZ = 500
4. xy-scaling, slab restraint, p_tgt_z = 1000 -> Pres-ZZ = 500

Conclusions:
1. p_tgt_z does not play a role, which is expected, since compressibility_z = 0.
2. The xy-semiisotropic barostat only matches the x and y pressure component. This makes sense and should, in the retrospect, have been clear to me, since the z pressure target can be independently set.

For my slab system with xy-scaling this means, that the overall x- and y-pressure is correct. Because of the walls the pressure in the slab will be higher than outside (osmotic pressure), e.g. 1500 and 500 bar. That still means, that also in z direction, if the virial would ignore the restraints, the pressure should be 1000 bar, right? Appereantly it isn't and i think the virial includes the forces from the restraints. That is why the total pressure reported does not match 1000 bar, even though i would argue it is 1000 bar.

For the slab system with z-scaling this means, that the overall x- and y-pressure are off, since the virial takes into account the restraint forces and thereby the pressure in the system is to high.
That is also what i observe:
Pres-XX 1502.12
Pres-YY 1505.27
Pres-ZZ 993.386

So i still think this is a bug, even though I (think I) understand it better now and can work around it.

#4 Updated by Berk Hess 13 days ago

This is a complex issue. I think it depends on what you define should happens to the restraints when you scale the volume. If the position restraint reference coordinates scale along, which would be the case with the option refcoord-scaling=all, I think the restraints should contribute to the virial. If you do not scale the position restraint coordinates, the virial and pressure might be ill defined.

#5 Updated by Marvin Bernhardt 12 days ago

You are right, this is a complex issue and in the end is a question of definition.

The reason, why I was puzzled in the beginning was that i thought of the virial pressure as a sum over pair interactions only. Therefor i thought the wall restraints had no direct effect on the virial, since it is not really a pair interaction, but more a spring to a fixed position.

I think it would be more consistent if forces exerted by fixed positions (which have no counter force) would not contribute to the virial. But that's my view and it might break something somewhere else. If you decide to not change this, i think there should be a (preprocessor) warning that the virial pressure is not representing the "real" pressure, if there are restraints to fixed positions. I guess there is only an effect if there are multiple particles restrained in a way that they constantly are pushed together or pulled apart.

A related issue was found by a colleague of mine. With a periodic polymer in z-direction and a xy-scaling barostat, the total pressure is always below the expected result, and I guess it is because of the polymer pulling together in z-direction, thereby affecting the virial.

#6 Updated by Marvin Bernhardt 12 days ago

Another system one might consider is a cuboid of fluid surrounded by vacuum as explained in the restraints help1.

Without taking the wall forces into the virial and allowing for refcoord-scaling, this would allow the user to control the pressure of the fluid within the cuboid, right?

[1] http://manual.gromacs.org/current/reference-manual/functions/restraints.html

#7 Updated by Berk Hess 12 days ago

The periodic polymer case is correct. If the polymer can not deform sufficiently, you simply can not get the requested pressure.

If am not sure if setting the virial constributions of restraints with refcoord-scaling=none to zero gives consistent results. If it does, we should simply remove the contribution. You can achieve this with the following changes to the code:

index 0fa2a2fed..76530666c 100644
--- a/src/gromacs/listed-forces/position-restraints.cpp
+++ b/src/gromacs/listed-forces/position-restraints.cpp
@ -304,7 +304,10 @ real fbposres(int nbonds,
}
}

- forceWithVirial->addVirialContribution(virial);
+ if (refcoord_scaling erscALL)
+ {
+ forceWithVirial->addVirialContribution(virial);
+ }

return vtot;
}
@ -386,7 +389,7 @ real posres(int nbonds,
}
}

- if (computeForce)
+ if (computeForce && refcoord_scaling erscALL) {
forceWithVirial->addVirialContribution(virial);
}

#8 Updated by Marvin Bernhardt 11 days ago

Thanks for pointing out, where in the code the virial for the restraints are computed.

Since I have refcoord_scaling on, your proposed change does not change anything for my systems.

For the reasons I mentioned before I think restraints to fixed positions should not contribute to the virial at all. So I would boldly, for the sake of discussion, put forward the following changes, even though, I have not thought through, if there are other systems, where this leads to inconsistencies. With these changes i get the expected behavior (similar pressure in x, y, z, tot) for all three scaling options (iso, xy, z).

index 460b055d9..cbc31371b 100644
--- a/src/gromacs/listed_forces/position_restraints.cpp
+++ b/src/gromacs/listed_forces/position_restraints.cpp
@@ -215,7 +215,6 @@ real fbposres(int nbonds,

     rvec *f      = as_rvec_array(forceWithVirial->force_.data());
     real  vtot   = 0.0;
-    rvec  virial = { 0 };
     for (i = 0; (i < nbonds); )
     {
         type = forceatoms[i++];
@@ -299,13 +298,9 @@ real fbposres(int nbonds,
         for (m = 0; (m < DIM); m++)
         {
             f[ai][m]  += fm[m];
-            /* Here we correct for the pbc_dx which included rdist */
-            virial[m] -= 0.5*(dx[m] + rdist[m])*fm[m];
         }
     }

-    forceWithVirial->addVirialContribution(virial);
-
     return vtot;
 }

@@ -354,8 +349,6 @@ real posres(int nbonds,
         f              = as_rvec_array(forceWithVirial->force_.data());
     }
     real        vtot   = 0.0;
-    /* Use intermediate virial buffer to reduce reduction rounding errors */
-    rvec        virial = { 0 };
     for (i = 0; (i < nbonds); )
     {
         type = forceatoms[i++];
@@ -381,16 +374,10 @@ real posres(int nbonds,
             if (computeForce)
             {
                 f[ai][m]  += fm;
-                virial[m] -= 0.5*(dx[m] + rdist[m])*fm;
             }
         }
     }

-    if (computeForce)
-    {
-        forceWithVirial->addVirialContribution(virial);
-    }
-
     return vtot;
 }

#9 Updated by Berk Hess 11 days ago

Removing the virial contribution from scaling restraints is wrong. The pressure is the derivative of the free-energy with respect to the volume. So the restraint potential should contribute to the virial. You might not want that, but then you don't want the "real" virial and pressure, but a pressure like quantity that excludes this contribution.

#10 Updated by Marvin Bernhardt 11 days ago

That is a good argument. I took a second look at the manual and the code and I think my intuition was wrong. Thanks for clearing it up for me.

For weather your proposed change should be implemented:
I think it would be a counter-intuitive outcome for my slab system and xy-scaling. The pressure would then depend on the option refcoord-scaling, while the forces (the hamiltonian) would not be affected, because there is no scaling in z-direction.

So from my perspective it would be ok, to leave it as it is (and close this issue). Maybe you want to add a grompp note or something.
Thank you again for looking into this for me.

#11 Updated by Berk Hess 11 days ago

Yes, but that is how it is. If you would restrain two layers of atoms to form walls to put all other atoms in between and do not scale the position restraint reference coordinates with the volume scaling, the potential will not change. I am not sure that in general we should not add the virial contribution when we do not scale the reference coordinates.

At least the comment in the manual: "The reference coordinates for position restraints are not modified. Note that with this option the virial and pressure will depend on the absolute positions of the reference coordinates." is incorrect. The virial does not depend on absolute positions. I think this used to be the case several years ago.

#12 Updated by Berk Hess 8 days ago

  • Category set to documentation
  • Status changed from New to Fix uploaded
  • Assignee set to Berk Hess
  • Target version set to 2019.3

I concluded that one can generally not define a virial and pressure that is consistent with the free-energy with position restraints and refcoord-scaling is none or com. I changed the mdp user guide to say the pressure might be ill-defined in such cases. I left in the virial contribution from the position restraints to the virial in these cases, as this lowers the pressure and the pressure is even lower without restraints, so the results of an NPT equilibration run is likely closer to what the user would want. This also keeps the virial independent from the refcoord-scaling option.

#13 Updated by Marvin Bernhardt 8 days ago

Thanks!

Also available in: Atom PDF