Hassan Shallal wrote:
Dear Gromacs users,
I have some concerns about the both the pressure fluctuations and
averages I obtained during the equilibration phase. I have already read
through several similar posts as well as the following link
http://www.gromacs.org/Documentation/Terminology/Pressure. I understand
the pressure is a macroscopic rather than instantaneous property and the
average is what really matters. I also found out through similar posts
that negative average pressure indicates the system tendency to contract.
In the above link, it mentioned that pressure fluctuations should
decrease significantly with increasing the system's size. In my cases, I
have a fairly big systems (case_1 with *17393* water molecules
and case_2 with *11946 *water molecules). However, the pressure still
has huge fluctuations (around 500 bars) from the reference value (1
bar). Here are the average pressure and density values resulting from
the equilibration phases of two cases, please notice the negative
average pressure values in both cases...
Case_1_pressure:
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Pressure *-2.48342* 0.92 369.709 -4.89668
(bar)
Case_1_density:
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Density 1022.89 0.38 3.8253 2.36724
(kg/m^3)
Case_2_pressure:
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Pressure *-8.25259* 2.6 423.681 -12.1722
(bar)
Case_2_density:
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Density 1034.11 0.37 2.49964 1.35551
(kg/m^3)
So I have some questions to address my concerns:
1- each of the above systems has a protein molecule, NaCl to give 0.15 M
system and solvent (water) molecules... Could that tendency to contract
be an artifact of buffering the system with sodium and chloride ions?
I suppose anything is possible, but given that these are fairly standard
conditions for most simulations, I tend to doubt it. My own (similar) systems
do not show this problem.
2- how to deal with the tendency of my system to contract? Should
I change the number of water molecules in the system?
or
Is it possible to improve the average pressure of the above systems by
increasing the time of equilibration from 100 ps to may be 500 ps or
even 1 ns?
3- Is there a widely used range of average pressure (for ref_p = 1 bar)
that indicates acceptable equilibration of the system prior to the
production?
To answer #2 and #3 simultaneously - equilibration is considered "finished" when
your system stabilizes at the appropriate conditions (usually temperature and
pressure). Your results indicate that your equilibrium is insufficient.
4- I can't understand how the system has a tendency to contract whereas
the average density of the solvent is already slightly higher than it
should be (1000 kg/m^3).
The contraction causes the density to rise. Pressure and density are not
independent; density is a result of pressure.
I would like to ignore the pressure based judgement of the above
equilibration given that the average density values are very close to
the natural value (1000 kg/m^3) (by the way I am using tip3p water model
with CHARMM27 ff) Any comment!!
It is not guaranteed that simulation water models will reproduce the real
(experimental) density of water. If memory serves, the expected density of
TIP3P should be ~0.98 g mL^{-1}, but I could be wrong.
5- Is the huge fluctuation of the pressure values of the above system
despite thier large sizes still acceptable? or large fluctuation is only
acceptable for small size systems and is unacceptable for large size
systems?
If it is unacceptable, any idea of how could it be alleviated or minimized?
The fluctuations seem reasonable. You might try equilibrating with the
Berendsen barostat first, then switching to Parrinello-Rahman. The P-R barostat
allows for wider fluctuations, so if the system is poorly equilibrated, your
system will be slower to converge.
-Justin
I am including the .mdp used in the above equilibration in case it is
needed.
Any feedback or response to the above questions is so much appreciated..
Great regards
Hassan
.mdp used for the above equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 25000 ; 4 * 25000 = 100 ps
dt = 0.004 ; 4 fs, virtual sites along with heavy hydrogens are used
; Output control
nstxout = 100 ; save coordinates every 0.2 ps
nstvout = 100 ; save velocities every 0.2 ps
nstenergy = 100 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 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 = 6 ; also related to accuracy, changed to 6 because of
using virtual sites along with a larger time step
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 20 fs
rlist = 1.2 ; short-range neighborlist cutoff, equal to rcoulomb to
allow for PME electrostatics (in nm)
; Lennard-Jones
vdwtype = switch ; VDW interactions are switched of
between 1 and 1.2
rvdw_switch = 1 ;
rvdw = 1.2 ; short-range vdw cutoff, optimal for
CHARMM27 ff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
rcoulomb = 1.2 ; short-range electrostatic cutoff, optimal for CHARMM27
ff (in nm)
; 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 = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 1 ; in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for switch vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list [email protected]
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 [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
