[gmx-users] Pressure very high with chloroform and gromacs-4.5.5

[Mark Abraham]

On 12/10/2011 8:04 PM, Nuno Azoia wrote:
> Thank you Mark for your answer!
>
> I agree with you when you say that the high pressure is not itself a
> problem. To avoid that problem I change my system from 1000 molecules
> to 27000 molecules, and the pressure change from several thousands to
> several hundreds bar.
> What I found very strange was the increasing volume (and pressure). My
> system is very unstable, I know, because the density is very low, and
> not the opposite, so I was expecting to see the box getting smaller.
> Looking to the trajectory I can see a box almost empty (with empty I
> mean empty space, with chloroform molecules aggregating in small
> droplets), and that's why I found the system behavior very strange.

Sounds like your initial density might be far too low. grompp reports 
it, so do check.

Mark
>
> I will follow your suggestion concerning nsttcouple and nstpcouple and
> the other initial conditions.
>
> Nuno Azoia
>
>
> On Wed, Oct 12, 2011 at 6:24 AM, Mark Abraham<Mark.Abraham at anu.edu.au>  wrote:
>> On 12/10/2011 2:22 AM, Nuno Azoia wrote:
>>> Hello!
>>>
>>> I found something very strange while making a CHCl3 box using
>>> gromacs-4.5.5.
>>> A look the mailing list, the manual and some release notes for
>>> gromacs-4.5 and I couldn't found the answer for my problem. It's
>>> possible that I'm doing something wrong, but I can not find what, so
>>> I'm describe my problem.
>>>
>>> I start a chloroform box from scratch, using genconf, and I get o
>>> chloroform box with 1000 molecules. I get energy minimization without
>>> problems.
>>> Then I've run some equilibration steps in a NVT ensemble and in the
>>> end I get pressure in the order of hundreds of bar.
>> The high pressure is not itself a problem - small numbers of molecules and
>> short simulation times lead to doubtful statistics for pressure.
>>
>>>   Then I change to
>>> NPT conditions and both the pressure and the volume keep increasing
>>> with time.
>> Be sure to visualize your trajectory to confirm its behaviour matches what
>> you expect from the trends in P and V.
>>
>>> To discard the possibility of a size problem, I repeat everything with
>>> a box of 27000 molecules, with a volume of ~13800 nm^3. The problem
>>> was the same. Very high pressures (150-200 bar) and very low densities
>>> (<    200 g/L) after 750 ps simulation time. And both volume and pressure
>>> increasing with time.
>>>
>>>
>>> I'm doing now the same procedure, but using gromacs-4.0.7 and I'm
>>> getting very different (and better) results. After energy minimization
>>> I run 50000 steps in a nvt ensemble and I got pressure around -30 bar
>>> (Ok for me). After that I start to run the simulations in npt ensemble
>>> and the pressure start to increase slowly, with negative values
>>> because the system have very low densities (~400 g/L), and the volume
>>> is decreasing. So I'm getting the normal reaction from the system.
>>>
>>>
>>>
>>> Where is the problem? There are some different parameters to set in
>>> the mdp file and I didn't realize that, or is this a problem in
>>> gromacs-4.5?
>> It seems you are generating some numerical instability with your choice of
>> initial and simulation conditions, and that gromacs-4.0 was luckier in
>> avoiding the problem. If this is a parallel simulation, then the way
>> communication is managed in order to organize T and P coupling changed going
>> up to 4.5. For equilibration, I would encourage setting nsttcouple and
>> nstpcouple to 1. Using a shorter time step can help with numerical
>> instability. Berendsen T coupling or a larger coupling time constant might
>> also help. After gentler NVT, then switch to NPT, let the box size
>> stabilize, and finally move to larger coupling periods, time steps and
>> v-rescale T coupling.
>>
>> Mark
>>
>>>
>>> In both cases I used this parameters:
>>> -------------------------------------------------------------------------
>>> integrator          =  md
>>> dt                    =  0.002
>>> nsteps              =  50000
>>> nstcomm          =  1
>>> nstxtcout           =  500
>>> xtc-precision      =  1000
>>> nstxout             =  0
>>> nstvout             =  0
>>> nstfout             =  0
>>> nstlog              =  500
>>> nstenergy           =  500
>>>
>>> nstlist               =  5
>>> ns_type             =  grid
>>>
>>> ;Reaction field
>>> rlist               = 0.8
>>> coulombtype         = Reaction-field
>>> rcoulomb            = 1.4
>>> epsilon_r           = 1.0
>>> epsilon_rf          = 4.8
>>> vdwtype             = cut-off
>>> rvdw                = 1.4
>>>
>>> ; temperature coupling
>>> Tcoupl              =  v-rescale
>>> tc-grps             =  CHCL3
>>> tau_t               =  0.05
>>> ld_seed             =  -1
>>> ref_t               =  300
>>>
>>> ; Energy monitoring
>>> energygrps          =  CHCl3
>>>
>>> ; Isotropic pressure
>>> Pcoupl              =  no
>>> Pcoupltype          =  isotropic
>>> tau_p               =  0.5
>>> compressibility     =  4.5e-5
>>> ref_p               =  1.0
>>>
>>> ; Generate velocites is on at 300 K.
>>> gen_vel             =  yes
>>> gen_temp            =  300
>>> gen_seed            =  -1
>>>
>>> constraint_algorithm=lincs
>>> lincs_order         =  4
>>> lincs-warnangle     = 90
>>> constraints         =  all-bonds
>>> ------------------------------------------------------------------------
>>>
>>> and of course, for the npt ensemble I just change
>>> _________________________________________
>>> Pcouple     = berendsen
>>>
>>>
>>> gen_vel      = no
>>> ------------------------------------------------------------------------
>>>
>>> Thank you in advance
>>>
>> --
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> Please don't post (un)subscribe requests to the list. Use the www interface
>> or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
>