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

[Nuno Azoia]

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.

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
>



-- 
Nuno Gonçalo Azoia Lopes

Rua Nascente, nº 25
4705-473 Escudeiros

Tel: 253 683 198
Móvel: 965 382 487