[gmx-users] Re: Problem with pressure coupling

[David van der Spoel]

On 2010-09-29 21.15, Ondrej Marsalek wrote:
> Hi again,
>
> I have made available a testcase that reproduces the problem here:
>
> http://marge.uochb.cas.cz/~marsalek/tmp/test-volume-drift.tar.bz2
>
> To see the problem, use grompp
> from GMX 4.5.1 and run the resulting tpr. Then extract the volume:
>
> echo Volume | g_energy -f ener.edr -o volume
>
> and observe the drift. I would be very grateful if someone could help
> me confirm the behavior.
>
> I have found one more thing indicative of this being a bug in 4.5.1.
> If I use grompp from 4.0.7 and then run the tpr using 4.5.1 it runs
> fine, just as it runs in 4.0.7. The tpr produced by grompp from 4.5.1
> has the drift in box volume when run in 4.5.1. Of course, it will not
> run in 4.0.7, so there is no way of testing that.
>
> Thanks for any help on this,
> Ondrej

I have confirmed this as a bug.
In 4.0.5 (which I had handy) the density is constant
Density (SI)   999.554    3.82377    3.80651  0.0628345    1.25694

In 4.5.1 with 4.0.5 tpr file it is as you say:
Density        995.277        1.3    4.32625   -1.24517  (kg/m^3)

In 4.5.1 with 4.5.1 tpr file I get
Density        709.919         66    137.824   -467.285  (kg/m^3)


Now the differences in the tpr files are revealing:
[anfinsen:test-volume-drift] % gmxcheck -s1 topol405.tpr -s2 topol.tpr
<snip>
Note: tpx file_version 58, software version 73
Reading file topol.tpr, VERSION 4.5.1-%ci-21e8a32 (single precision)
comparing inputrec
inputrec->nstcalcenergy (1 - 10)
inputrec->rgbradii (2.000000e+00 - 1.000000e+00)
inputrec->sa_surface_tension (2.092000e+00 - 2.050160e+00)
inputrec->sc_sigma_min (0.000000e+00 - 3.000000e-01)
inputrec->nstdhdl (1 - 10)
inputrec->separate_dhdl_file (0 - 10)
<snip>

After setting nstcalcenergy = 1 I get
Density        996.894        1.2    3.85083    1.86718  (kg/m^3)

In other words the bug can be circumvented by setting this variable to 
1. I will file a bugzilla.

>
>
> On Thu, Sep 23, 2010 at 15:16, Ondrej Marsalek
> <ondrej.marsalek at gmail.com>  wrote:
>> Dear all,
>>
>> running a box of neat water with the below mdp file in GROMACS 4.5.1
>> results in a fast linear increase of the volume of the box, the water
>> getting torn into clusters. The same simulation in 4.0.7 (with
>> nstcalcenergy removed, of course) runs just fine, with the box volume
>> stable and fluctuating. In 4.5.1 the problem can be fixed by setting
>> nstlist = 1. I assume this would also be the case for other positive
>> numbers. The same behavior can be seen for the Parrinello-Rahman
>> barostat. I can of course provide full testcases ready to run, if
>> needed.
>>
>> This looks like a bug to me. Can anyone help me resolve this?
>>
>> Thanks,
>> Ondrej
>>
>>
>> ; RUN CONTROL PARAMETERS
>> integrator               = md
>> tinit                    = 0
>> dt                       = 0.002
>> nsteps                   = 5000
>> comm-mode                = Linear
>> nstcomm                  = 10
>>
>> ; OUTPUT CONTROL OPTIONS
>> nstxout                  = 0
>> nstvout                  = 0
>> nstfout                  = 0
>> nstlog                   = 100
>> nstcalcenergy            = 10
>> nstenergy                = 100
>> nstxtcout                = 100
>> xtc-precision            = 1000
>>
>> ; NEIGHBORSEARCHING PARAMETERS
>> nstlist                  = -1
>> ns_type                  = grid
>> pbc                      = xyz
>> periodic_molecules       = no
>> rlist                    = 1.3
>>
>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>> coulombtype              = PME-Switch
>> rcoulomb-switch          = 0.6
>> rcoulomb                 = 0.8
>> vdw-type                 = Shift
>> rvdw-switch              = 0.6
>> rvdw                     = 0.8
>> DispCorr                 = EnerPres
>> fourierspacing           = 0.12
>> pme_order                = 4
>> ewald_rtol               = 1e-05
>> ewald_geometry           = 3d
>> optimize_fft             = yes
>>
>> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
>> Tcoupl                   = v-rescale
>> tc-grps                  = system
>> tau_t                    = 0.5
>> ref_t                    = 300
>> Pcoupl                   = berendsen
>> Pcoupltype               = isotropic
>> tau_p                    = 1
>> compressibility          = 4.5e-5
>> ref_p                    = 1.0
>>
>> ; GENERATE VELOCITIES FOR STARTUP RUN
>> gen_vel                  = no
>>
>> ; OPTIONS FOR BONDS
>> constraints              = none
>>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se