[gmx-users] Simulated annealing in implicit solvent

[Justin Lemkul]

On 8/7/12 8:10 PM, Ольга Кононова wrote:
>
> 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

Sounds like a known issue:

http://redmine.gromacs.org/issues/757

> 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?

Your results will be very inaccurate and your simulations unstable.  I have 
found the only way to obtain a reliable trajectory using implicit solvent is to 
use the optimized all-vs-all kernels.  Set:

pbc = no
rlist = 0
nstlist = 0
rvdw = 0
rcoulomb = 0
rgbradii = 0

Pertaining to the issue you described, are you running on GPU or CPU?  If it's 
GPU, I find the same behavior - regardless of the target temperature, mdrun 
always reports a temperature 10 K higher than the target.  I cannot determine if 
it is an output bug or a calculation bug.  The behavior you describe is 
consistent with simulated annealing not working with implicit solvent (I don't 
know whether it should or not) and thus getting stuck at 0 K and thus reporting 
10 K throughout the simulation.

If you're running on CPU, ignore the previous paragraph ;)

-Justin

> 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
>

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================