[gmx-users] ACETONITRILE with CHARM ff

[Sonia Milena Aguilera Segura]

> 
> Message: 3
> Date: Tue, 11 Jul 2017 12:57:05 -0400
> From: Justin Lemkul <jalemkul at vt.edu>
> > Moreover, I observed that this time I got lower values for P during the NPT
> > equilibration, but still is too far from 1 bar.  I really don't understand
> > why for the NVT simulation I get a T around 298, but as soon as I turn on
> > the pcoupl, the T rises to 300-301 K and the P gets average values of 7
> > and 4 bar (vs 8 and 14 for the previous simulations). Then at the end of
> > the 1-ns MD the temperature remains around 301 and the P is -1 and 2.7
> > bar. Considering the parameters I am using, is there anything I can change
> > to make the P coupling better? I am running a 3 nm box with 308 molecules.
> > This is the full mdp file:
> > 
> 
> http://www.gromacs.org/Documentation/Terminology/Pressure
> 
> Your box is very small and will be subject to large fluctuations.

Dear Justin, 

I increased the size of the simulation box to 4 nm. Indeed, the values of pressure improved (averages around -1 bar or 1.5, 2 bar, or so). However, the T keeps being overestimated (302 K). During NVT I got the right value, so I decided to run the 200 ns NPT equilibration with Berendsen barostat instead of P-R. I got the right T value, around 298.3 and a pressure of 1.5 more or less. Then I launched a continuation test with 200 ps MD with P-R (I couldn't use the -e option because it gives me the error Could not find energy term named 'Box-Vel-XX', moreover, I didn't use P-R before so I guess I shouldn't expect stored values for it?. But I am using -t ). Despite I got low pressure averages around 0.5 bar, the temperature raised again to 301-302. This happens very early in the simulation, which seems to indicate that for sd integrator/thermostat and P-R barostat there is something going on. I am getting practically the same density in between simulation, so I guess I can say that in term of P, the system has been equilibrated. However, what can I do to get the right T for this system? If I was able to get the right T with Berendsen barostat, I don't understand what's wrong when I change to P-R. 

Thank you in advance, 

Sonia Aguilera 
PhD student
ENSCM

> 
> -Justin
> 
> > ; Run control
> > integrator               = sd       ; Langevin dynamics
> > tinit                    = 0
> > dt                       = 0.0005
> > nsteps                   = 2000000   ; 1 ns
> > nstcomm                  = 100
> > ;energygrps		 = non-Water
> > ; Neighborsearching and short-range nonbonded interactions
> > cutoff-scheme            = verlet
> > nstlist                  = 20
> > ns_type                  = grid
> > pbc                      = xyz
> > rlist                    = 1.2
> > ; Electrostatics
> > coulombtype              = PME
> > rcoulomb                 = 1.2
> > ; van der Waals
> > vdwtype                  = cutoff
> > vdw-modifier             = force-switch
> > rvdw-switch              = 1.0
> > rvdw                     = 1.2
> > ; Apply long range dispersion corrections for Energy and Pressure
> > DispCorr                  = no
> > ; Spacing for the PME/PPPM FFT grid
> > fourierspacing           = 0.12
> > ; EWALD/PME/PPPM parameters
> > pme_order                = 6
> > ewald_rtol               = 1e-06
> > epsilon_surface          = 0
> > ; Temperature coupling
> > ; tcoupl is implicitly handled by the sd integrator
> > tc_grps                  = system
> > tau_t                    = 1.0
> > ref_t                    = 298.15
> > ; Pressure coupling is on for NPT
> > Pcoupl                   = Parrinello-Rahman
> > tau_p                    = 1.0
> > compressibility          = 4.5e-05
> > ref_p                    = 1.0
> > ; Do not generate velocities
> > gen_vel                  = no
> > ; options for bonds
> > constraints              = none  ; we only have C-H bonds here
> > ; Type of constraint algorithm
> > constraint-algorithm     = lincs
> > ; Constrain the starting configuration
> > ; since we are continuing from NPT
> > continuation             = yes
> > ; Highest order in the expansion of the constraint coupling matrix
> > lincs-order              = 12
> > 
> > 
> > Thank you very much,
> > 
> > Sonia Aguilera
> > PhD student
> > ENSCM
> >> ; Run control
> >> integrator               = sd       ; Langevin dynamics
> >> tinit                    = 0
> >> dt                       = 0.0005
> >> nsteps                   = 40000000   ; 20 ns
> >> nstcomm                  = 100
> >> ; Neighborsearching and short-range nonbonded interactions
> >> cutoff-scheme            = verlet
> >> nstlist                  = 20
> >> ns_type                  = grid
> >> pbc                      = xyz
> >> rlist                    = 1.2
> >> ; Electrostatics
> >> coulombtype              = PME
> >> rcoulomb                 = 1.2
> >> ; van der Waals
> >> vdwtype                  = cutoff
> >> vdw-modifier             = potential-switch
> >> rvdw-switch              = 1.0
> >> rvdw                     = 1.2
> >> ; Apply long range dispersion corrections for Energy and Pressure
> >> DispCorr                  = EnerPres
> > 
> > CHARMM uses a force switch, and only apply dispersion correction in the
> > case of
> > lipid monolayers.
> > 
> > http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM
> > 
> > -Justin
> > 
> > ==================================================
> > 
> 
> --
> ==================================================