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

[Mark Abraham]

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
>