GROMACS version: 2021.4
GROMACS modification: Yes/No
Here post your question : I want to simulate my system at different temperature, say 290K. So I have to change the ref_t from 300K to 290K in nvt.mdp, npt.mdp and md.mdp files. In the nvt.mdp, do I need to change the gen_t (Maxwell distribution) to 290K too? The following is the nvt.mdp file-
title = charmm36 CYHE NVT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
I’m pretty sure you should change it as well, gen_t decides what velocities to assign to each atom so that it matches the maxwell distribution of the temperature you assign and ref_t is the temperature that the algorithm will try to move the system towards during equilibration, keeping the temperature consistent should provide the best results and avoids unnecessary fluctuations
To add to Karis’s answer, when you are changing the temperature relatively little, like in this case, it might be better not to regenerate all velocities. You would certainly reach the target temperature quicker if you regenerate them, but the equilibration might take longer.
So, either should work fine. If it has already taken a long time to equilibrate the system I would probably avoid regenerating the velocities.
I changed both the temperatures (ref_t and gen_temp) to 290K for protein-water system. After NVT equilibration (100 ps) I got the following plot for temperature-
It’s difficult for us to say. But in general, if you are not sure if your system is equilibrated it might be worth equilibrating a bit longer. Are you equilibrating a protein-water system for 100 ps NVT and 100 ps NPT? That sounds short to me.
One more thing, why are you using such a small tau-t? I wonder if that’s the cause of your large temperature fluctuations. I would recommend 1ps or even 2ps.
Yes, I am equilibrating a protein-water system (protein is ~8.2 kDa, cubic box). I equilibrated for 100 ps for both NVT and NPT. Then can you please suggest for how much time step do I need to equilibrate my system at temperature 290K and 310K? Is tau_t, the only change should I make. I am attaching my .mdp files for your reference.
NVT.mdp-
title = charmm36 Protein-water system NVT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
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 = 290 290 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 290 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
NPT.mdp-
title = charmm36 Protein-water system NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = yes ; continuing from NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
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 = 290 290 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = C-rescale ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; 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 is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; velocity generation off after NVT
I can’t say for how long you should equilibrate your system.
As far as I can see the settings in general look OK (no guarantees - you must take the responsibility for your own settings). The only things I would change are tau-t (to 1 or 2) and tau-p (to 5 or 10). I’m curious to learn where these low coupling constants for v-rescale and c-rescale are recommended.
As I am using the charmm36 forcefield, so I used most of the input file from the protein-ligand MD tutorial (V-rescale constant was recommended there, and pressure constant was actually recommended for Parrinello-Rahman barostat, but as you suggested me to change it to C-rescale for cyclohexane box of biphasic system, I used the same barostat here as I have to compare the two systems).
