Temperature and pressure values are extremely high during equilibration

User discussions
1

GROMACS version: 2023.5 (single precision)
GROMACS modification: Yes/No

The output pressure.xvg file from the following commands appears to have very high values ranging from 4,000-6,000 bars.

gmx grompp -f amber/npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -maxwarn 1000
gmx mdrun -ntmpi 1 -v -deffnm npt -nb gpu -gpu_id 0
gmx energy -f npt.edr -o pressure.xvg -xvg none

npt.mdp

title                   = Amber NPT equilibration  
define                  = -DPOSRES  ; position restrain the protein

; Parameters describing what to do, when to stop and what to save
integrator              = md        ; leap-frog integrator
dt                      = 0.001     ; 0.001 = 1 fs
nsteps                  = 100000    ; nsteps = 100 ps / dt (ps) ; aim for 100 ps​
nstenergy               = 500       ; save energy and temperature every 1.0 ps

continuation            = yes      

; Pressure coupling is on
pcoupl                  = C-rescale             ; Pressure coupling on in NPT
pcoupltype              = isotropic             ; uniform scaling of box vectors
tau_p                   = 5.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

; Keep system temperature fluctuating physically correct
tcoupl                  = V-rescale           ; modified Berendsen thermostat
tc-grps                 = system  ; coupling groups 
tau_t                   = 1.0     ; time constant, in ps
ref_t                   = 300     ; reference temperature, one for each group, in K

; Settings that make sure we run with parameters in harmony with the selected force-field
constraints             = h-bonds   ; bonds involving H are constrained
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
vdw-modifier            = Potential-shift-Verlet ; Amber specific
DispCorr                = EnerPres  ; account for cut-off vdW scheme
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
fourierspacing          = 0.125     ; grid spacing for FFT

; Settings for LINCS
lincs-order             = 8 ; default = 4
lincs-iter              = 4 ; default = 1

A similar situation with temperature with values hovering at around 17,000k

nvt.mdp

title                   = Amber NVT equilibration 
define                  = -DPOSRES  ; position restrain the protein

; Parameters describing what to do, when to stop and what to save
integrator              = md        ; leap-frog integrator
dt                      = 0.001     ; 0.001 = 1 fs
nsteps                  = 100000    ; nsteps = 100 ps / dt (ps) ; aim for 100 ps​
nstenergy               = 500       ; save energy and temperature every 1.0 ps

; Keep system temperature fluctuating physically correct
tcoupl                  = V-rescale           ; modified Berendsen thermostat
tc-grps                 = system  ; coupling groups 
tau_t                   = 1.0     ; time constant, in ps
ref_t                   = 300     ; reference temperature, one for each group, in K

; Pressure coupling is off
pcoupl                  = no

; Velocity generation
gen_vel                 = yes                 ; assign velocities from Maxwell distribution
gen_temp                = 300                 ; temperature for Maxwell distribution

; Settings that make sure we run with parameters in harmony with the selected force-field
constraints             = all-bonds ; h-bonds = bonds involving H are constrained, all-bonds = all bonds are constrained
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
vdw-modifier            = Potential-shift-Verlet ; Amber specific
DispCorr                = EnerPres  ; account for cut-off vdW scheme
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
fourierspacing          = 0.125     ; grid spacing for FFT

; Settings for LINCS
lincs-order             = 8 ; default = 4
lincs-iter              = 4 ; default = 1

The simulation itself can be completed without issue, there are no observable stereochemical issues throughout the process and I am able to obtain seemingly normal results. It is possible that I may be interpreting something incorrectly as I am a bit of a gromacs newbie so any guidance would be greatly appreciated.

2

A high pressure when you run with fixed volume is not surprising. Running NpT will give you the pressure you ask for.

But a temperature of 17000 K when the thermostat is set to 300 K sounds very strange. I would think this is impossible. If something is wrong, I would expect the simulation to crash, not to have an extremely high temperature deviation. Is the temperature fluctuating around 17000 K or shooting up rapidly?

3

Thank you for your response. In that case it is only the temperature that is incorrect. In all the simulations I have done with this configuration the temperature has always consistently fluctuated around some arbitrarily high temperature. There are no major fluctuations or rapid increases / decreases of temperature at the start or end.

These are the commands I use to obtain temperature:

gmx grompp -f amber/nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -maxwarn 1000
gmx mdrun -ntmpi 1 -debug -v -deffnm nvt -nb gpu -gpu_id 0
gmx energy -f nvt.edr -o temperature.xvg -xvg none -b 20

Am I perhaps misinterpreting the second column as the temperature in kelvin?

4

Could you post the list of output data you get from the gmx energy and what you chose to export? It might help if you present the output you get from the command as well.

You can also find the average temperature of your simulation near the end of the md.log file. If that is different from what you get using gmx energy, you are probably outputting the wrong energy data.

5

I just figured out the issue, essentially the wrong option was being selected after running the command. The output is now normal.

Thanks for the assistance both of you and my apologies for the issue stemming from an oversight on my end.