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? 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? 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). 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!! 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? 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
-- 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
