[gmx-users] Simulated Annealing

[Martina Bertsch, PhD]

Greetings.

I would like to perform the following simulated annealing protocol on a 
GPCR model in vacuum, while keeping the C-alpha atoms of the 
transmembrane (TM) helices restrained (fixed?):

1) Quick heating from 300-600 K in 1 ps with a 1 fs step.
2) Equilibration at 600 K for 100 ps.
3) Slow cooling to 500 K.
4) Equilibration at 500 K for 200 ps.
5) Slow cooling to 400 K.
6) Equilibration at 400 K for 200 ps.
7) Slow cooling to 300 K.
8) Equilibration at 300 K for 200 ps with gradual release of the 
restrained (fixed?) TM C-alpha atoms.
9) Minimization of the resulting structure.

My questions are:

a) Is there a way to program the whole protocol in a single mdp? Which 
mdp options do I use?
In the past, I used a simple SA cooling from 600 to 0 K in 1 ns (see 
enclosed mdp), but the papers suggest that a more gradual cooling is 
preferred.

2) During this process, I would like to keep just the transmembrane 
C-alpha atoms fixed, e.g., C-alpha in residues 6-30, 43-62, etc. Should I:
make_ndx -f GPCR_postsd.gro -o index.ndx
 > r 6-30 & a CA
 > r 43-62 & a CA
...
 >q
I seem to recall reading that an atom can only belong to one group. Does 
that mean that I have to split the original [ C-alpha ] into the [ 
TM_C-alpha ]
and the [nonTM_C-alpha] and delete the original [ C-alpha ] group? 
Wouldn't that also imply that I should get rid of the [ backbone ] group 
-- that sounds a bit suspicious?

3) Once I obtain the [ TM_C-alpha ] group, how do I fix it in the mdp? 
Should I just define:

freeze_groups [TM_C-alpha]
freeze_dim y y y

Or, should I prepare a restraint file TM_posres.itp and define 
-DPOSRES_TM in the mdp?

4) How would I encode the gradual release of the [TM_C-alpha] group at 
the last stage of the protocol?

5) If I did not want to fix any parts of the structure, how would I 
obtain the energy and rmsd for parts of the structure (with respect to 
the beginning structure), e.g., for TM helix 5, TM helix 7, etc?
I know that I could define groups TM1-TM7 in the index.ndx, and then for 
the rmsd calculation:
g_rms -s gpcr_SA.tpr -f gpcr_SA.trr -o TM7.xvg
and select the desired group at the prompt.
I am not sure how to obtain energy profiles for certain groups.

The "simple" SA protocol:
title = MT1 Simulated Annealing
cpp = /lib/cpp
constraints = all-bonds
integrator = md
dt = 0.001 ; ps !
nsteps = 1000000 ; total 1 ns.
;nstcomm = 1
nstxout = 2000
nstvout = 2000
nstfout = 2000
nstlog = 2000
nstenergy = 2000
nstxtcout = 2000
nstlist = 10
ns_type = grid
rlist = 1.0
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
pbc = xyz
comm_mode = linear
nstcomm = 1
fourierspacing = 0.15
pme_order = 4
optimize_fft = yes
comm_grps = protein
; Berendsen temperature coupling is on in one group
Tcoupl = berendsen
tau_t = 0.1
tc-grps = protein
ref_t = 600
; Energy monitoring
energygrps = protein
; Pressure coupling is not on
Pcoupl = no
pcoupltype = isotropic
tau_p = 5
compressibility = 4.5e-5
ref_p = 1.0
; Simulated annealing
annealing = yes
zero_temp_time = 1000000 ; ps
; Generate velocites is on at 600 K.
gen_vel = yes
gen_temp = 600.0
gen_seed = 280572

Regards,

Martina Bertsch