[gmx-users] Question about system equilibration

[Justin Lemkul]

On 10/27/16 1:43 PM, José Manuel Suárez Sierra wrote:
> Hi!
>
> I have been trying to equillibrate my system using first this NVT script:
>
> title		= OPLS NVT equilibration
> define		= -DPOSRES	; position restrain the protein
> ; Run parameters
> integrator	= md		; leap-frog integrator
> nsteps		= 50000		; 2 * 50000 = 100 ps
> dt		    = 0.002		; 2 fs
> ; Output control
> nstxout		= 500		; save coordinates every 1.0 ps
> nstvout		= 500		; save velocities every 1.0 ps
> nstenergy	= 500		; save energies every 1.0 ps
> nstlog		= 500		; update log file every 1.0 ps
> ; Bond parameters
> continuation	        = no		; first dynamics run
> constraint_algorithm    = lincs	    ; holonomic constraints
> constraints	            = all-bonds	; all bonds (even heavy atom-H bonds) constrained
> lincs_iter	            = 1		    ; accuracy of LINCS
> lincs_order	            = 4		    ; also related to accuracy
> ; Neighborsearching
> cutoff-scheme   = Verlet
> ns_type		    = grid		; search neighboring grid cells
> nstlist		    = 10		; 20 fs, largely irrelevant with Verlet
> rcoulomb	    = 1.0		; short-range electrostatic cutoff (in nm)
> rvdw		    = 1.0		; short-range van der Waals cutoff (in nm)
> ; Electrostatics
> coulombtype	    = PME	; Particle Mesh Ewald for long-range electrostatics
> pme_order	    = 4		; cubic interpolation
> fourierspacing	= 0.16	; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl		= V-rescale	            ; modified Berendsen thermostat
> tc-grps		= Protein Non-Protein	; two coupling groups - more accurate
> tau_t		= 0.1	  0.1           ; time constant, in ps
> ref_t		= 300 	  300           ; reference temperature, one for each group, in K
> ; Pressure coupling is off
> pcoupl		= no 		; no pressure coupling in NVT
> ; Periodic boundary conditions
> pbc		= xyz		    ; 3-D PBC
> ; Dispersion correction
> DispCorr	= EnerPres	; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel		= yes		; assign velocities from Maxwell distribution
> gen_temp	= 300		; temperature for Maxwell distribution
> gen_seed	= -1		; generate a random seed
>
>
> After this I tryid to equillibrate the system using this NPT script:
>
> title		= OPLS NPT equilibration
> define		= -DPOSRES	; position restrain the protein
> ; Run parameters
> integrator	= md		; leap-frog integrator
> nsteps		= 250000		; 2 * 50000 = 100 ps
> dt		    = 0.002		; 2 fs
> ; Output control
> nstxout		= 500		; save coordinates every 1.0 ps
> nstvout		= 500		; save velocities every 1.0 ps
> nstenergy	= 500		; save energies every 1.0 ps
> nstlog		= 500		; update log file every 1.0 ps
> ; Bond parameters
> continuation	        = yes		; Restarting after NVT
> constraint_algorithm    = lincs	    ; holonomic constraints
> constraints	            = all-bonds	; all bonds (even heavy atom-H bonds) constrained
> lincs_iter	            = 1		    ; accuracy of LINCS
> lincs_order	            = 4		    ; also related to accuracy
> ; Neighborsearching
> cutoff-scheme   = Verlet
> ns_type		    = grid		; search neighboring grid cells
> nstlist		    = 10	    ; 20 fs, largely irrelevant with Verlet scheme
> rcoulomb	    = 1.0		; short-range electrostatic cutoff (in nm)
> rvdw		    = 1.0		; short-range van der Waals cutoff (in nm)
> ; Electrostatics
> coulombtype	    = PME		; Particle Mesh Ewald for long-range electrostatics
> pme_order	    = 4		    ; cubic interpolation
> fourierspacing	= 0.16		; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl		= V-rescale	            ; modified Berendsen thermostat
> tc-grps		= Protein Non-Protein	; two coupling groups - more accurate
> tau_t		= 0.1	  0.1	        ; time constant, in ps
> ref_t		= 300 	  300	        ; reference temperature, one for each group, in K
> ; Pressure coupling is on
> pcoupl		        = berendsen	    ; Pressure coupling on in NPT
> pcoupltype	        = isotropic	            ; uniform scaling of box vectors
> tau_p		        = 0.5		            ; time constant, in ps
> ref_p		        = 1.0		            ; reference pressure, in bar
> compressibility     = 4.5e-5	            ; isothermal compressibility of water, bar^-1
> refcoord_scaling    = com
> ; Periodic boundary conditions
> pbc		= xyz		; 3-D PBC
> ; Dispersion correction
> DispCorr	= EnerPres	; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel		= no		; Velocity generation is off
>
>
> First I used Parrinello-Rahman, but system does not converges to 1 atm, I tried it in several scales of time until 500 ps, but system does not reach equilibrium, so I tried Berendesen, but same happened, it does not converges to 1 atm, I obtanied -2 atm in the best case. Temperature, volume and density converged easily, but pressure does not. I need help with this. ¿Can you help me?
>

http://www.gromacs.org/Documentation/Terminology/Pressure

See also the same exact question posted and answered within the last few days.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================