[gmx-users] Simulated annealing in implicit solvent

[Ольга Кононова]

Hello,
I got a problem, trying to make preparation (in particular, heating from 0 to 300 K) of my system in implicit solvent (GBSA OBC model). For gradual heating I use simulated annealing (see mdp below) and sd integrator, which allows me do not use thermostat. I also set up gen_vel = yes and gen_temp = 0 (I am not sure, is it correct way or not). The problem appears: on the 0th step in output line temperature is ~305 K and it drops to ~8 K during 0.7 ps of simulation after that it starts to fluctuate around constant value (~ 10 K). At the beginning in output file there is line 'Initial temperature: 0 K', but started mdrun shows temperature 305 K. At the same time 'Current ref_t for group System' changes correctly step by step.
I did try different parameters (turn on thermostat, turn off annealing, set gen_vel 0 or 300 K) - the same story. My guess is the problem with minimization. Algorithm "steepest descent" had convergent for 300 steps (in spite of very low emtol which I set up) and total potential energy dropped to constant minimum, I recompiled gromacs with double precision but minimization result was the same. Looks like my system (it's ~400000 atoms) fell in some local minimum which gives me high initial temperature.
What is wrong? Or it's proper behavior? Does it matter that I use cut-offs?
Here is mdp file:
define = -DFLEXIBLE
constraints = none
integrator = sd
dt = 0.001; 1fs
nsteps = 300000 ; 300 ps
nstcomm = 1
nstxout = 1000 ; frequency to write coordinates to output trajectory
nstvout = 0 ; frequency to write velocities to output trajectory; the last velocities are always written
nstfout = 0 ; frequency to write forces to output trajectory
nstlog = 10 ; frequency to write energies to log file
nstenergy = 100 ; frequency to write energies to edr file
nstcalcenergy = 100
vdwtype = cut-off
coulombtype = cut-off
pbc = no
table-extension = 20.0
nstlist = 100
ns_type = grid
rlist = 1.0
rcoulomb = 1.1
rvdw = 1.1
comm-mode = angular
comm-grps = system
optimize_fft = yes
;heating
annealing = single
annealing_npoints = 2
annealing_time = 0 300
annealing_temp = 0 300
 
ld_seed = 8072012
 
; V-rescale temperature coupling is on
Tcoupl = no
tau_t = 0.02
tc_grps = system
ref_t = 0
; Pressure coupling is off
Pcoupl = no
; Generate velocites is on
gen_vel = yes
gen_temp = 0
gen_seed = 8042012
; Implicit solvent
implicit_solvent = GBSA
gb_algorithm = OBC; Still ; HCT ; OBC
nstgbradii = 1
rgbradii = 1.0 ; [nm] Cut-off for the calculation of the Born radii. Currently must be equal to rlist
gb_epsilon_solvent = 80 ; Dielectric constant for the implicit solvent
sa_algorithm = Ace-approximation
sa_surface_tension = -1
 
I will be very appreciated for any help.
 
- Olga Kononova