Help with simulated annealing of a small peptide

User discussions
, ,
1

GROMACS version: 2024
GROMACS modification: No

Hi! I am trying to perform a simulation where I look into the stability of a small peptide in water with increasing temperature (i.e.: melting/denaturation curve).

For this simulation, I thought it would be useful to use the simulated annealing function. After looking a bit around, I have seen some suggestions to “turn on” simulated annealing during the equilibration phase by introducing the following lines into my npt.mdp file:

; Simulated annealing
annealing       = single single
annealing_npoints       = 6 6
annealing_time  = 0 500000 600000 1100000 1200000 1700000 0 500000 600000 1100000 1200000 1700000
annealing_temp  = 293 293 313 313 333 333 293 293 313 313 333 333

However, when I look into the temperature curve generated after NPT equilibration, the temperature does not change, and it is kept constant at 293 K. The number of steps for equilibration is 1700000, so the problem is not that not enough time was given for the temperature to change. Here is the rest of the npt.mdp file for context:

title                   = OPLS Simulated annealing - NPT equilibration
define                  = -DPOSRES  ; position restrain the protein
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 1700000     ; 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             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings
cutoff-scheme           = Verlet    ; Buffered neighbor searching
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)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; 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                   = 293     293           ; 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                   = 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
; Velocity generation
gen_vel                 = no        ; Velocity generation is off

Thus, I was wondering if someone knew or could provide any guidance as to what I am doing wrong? Why is the temperature not changing?

Another question I have is, when it comes to the time to do the simulation, should the md.mdp file also have simulated annealing turned on? Should it be the exact same lines as the npt.mdp file?

Is the any other (better) way to perform the simulation I had initially intended for?

Thanks for all the help!

2

Hi,
There is an error in your annealing input. The time units are in ps, not in steps. With your current setup you state that the first and second temperature coupling groups should be at 293K from 0 to 500ns (your simulation is 3.4ns).

Try:

; Simulated annealing
annealing       = single single
annealing_npoints       = 6 6
annealing_time  = 0 1000 1200 2200 2400 3400 0 1000 1200 2200 2400 3400
annealing_temp  = 293 293 313 313 333 333 293 293 313 313 333 333

I think that will fix your problems.

Regarding the production vs equilibration stages, it all depends on what you want to study. If you want to change the temperature during your production simulation you should include the annealing settings. Otherwise, you could run your production simulation starting from the output of the equilibration, but then you should almost certainly start at the final temperature in your annealing simulations (333K).

What is important is that you ask yourself what you are studying and why. From your description, I think it sounds like you should equilibrate (run the npt simulation) at a fixed temperature (293K?) to ensure a stable and equilibration starting configuration, and then use the annealing only in the production stage. Perhaps you want to run two production simulations, one with and one without annealing?