[gmx-users] Free energy perturbation parameters sc-alpha sc-power sc-sigma
Francesca Vitalini
francesca.vitalini11 at gmail.com
Wed Jan 22 17:10:09 CET 2014
Dear Gromacs users,
I am trying to perform a series of simulations where I gradually transform
a valine capped amino acid (Ac-VAL-NHMe) into an alanine capped amino acid
(Ac-ALA-NHMe), ie I have to turn on the HG and transform CG into HB. I have
defined 6 lambda values at which I would like to equilibrate my system
(lambda=0,0.2,0.40.6,0.8,1)
I am using GROMACS.4.5.5. and OPLS-AA force field.
I have been following the tutorial by Justin Lemkul
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html
I have problems with the nvt equilibration with lambda=0.8 and lambda=1 as
my systems systematically crashes. I believe that the reason is in the
tuning of the parameters sc-alpha sc-power sc-sigma in the mdp file.
I understand from the manual that sc_sigma should not be the default option
as I have 6 hydrogens disappearing, and consequently the sc-alpha sc-power
parameters will be effected as the interactions of the hydrogens in the
soft core state will get stronger.
However I have no experience in how to select these parameters and would
like to ask, if anyone has more experience, which could be plausible values.
Thank you very much,
Francesca
For completeness I attach here a copy of my topology and mdp files.
TOPOLOGY
; Include forcefield parameters
#include "oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Protein_chain_A 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB
; residue 1 ACE rtp ACE q 0.0
1 opls_135 1 ACE CH3 1 -0.18 12.011
opls_135 -0.18 12.011
2 opls_140 1 ACE HH31 1 0.06 1.008
opls_140 0.06 1.008
3 opls_140 1 ACE HH32 1 0.06 1.008
opls_140 0.06 1.008
4 opls_140 1 ACE HH33 1 0.06 1.008
opls_140 0.06 1.008
5 opls_235 1 ACE C 2 0.5 12.011
opls_235 0.5 12.011
6 opls_236 1 ACE O 2 -0.5 15.9994
opls_236 -0.5 15.9994
; residue 2 VAL rtp VAL q 0.0
7 opls_238 2 VAL N 3 -0.5 14.0067
opls_238 -0.5 14.0067
8 opls_241 2 VAL H 3 0.3 1.008
opls_241 0.3 1.008
9 opls_224B 2 VAL CA 3 0.14 12.011
opls_224B 0.14 12.011
10 opls_140 2 VAL HA 3 0.06 1.008
opls_140 0.06 1.008
11 opls_137 2 VAL CB 4 -0.06 12.011
opls_137 -0.18 12.011
12 opls_140 2 VAL HB 4 0.06 1.008
opls_140 0.06 1.008
13 opls_135 2 VAL CG1 5 -0.18 12.011
opls_140 0.06 1.008
14 opls_140 2 VAL HG11 5 0.06 1.008
opls_140 0.00 0.000
15 opls_140 2 VAL HG12 5 0.06 1.008
opls_140 0.00 0.000
16 opls_140 2 VAL HG13 5 0.06 1.008
opls_140 0.00 0.000
17 opls_135 2 VAL CG2 6 -0.18 12.011
opls_140 0.06 1.008
18 opls_140 2 VAL HG21 6 0.06 1.008
opls_140 0.00 0.000
19 opls_140 2 VAL HG22 6 0.06 1.008
opls_140 0.00 0.000
20 opls_140 2 VAL HG23 6 0.06 1.008
opls_140 0.00 0.000
21 opls_235 2 VAL C 7 0.5 12.011
opls_235 0.5 12.011
22 opls_236 2 VAL O 7 -0.5 15.9994
opls_236 -0.5 15.9994
; residue 3 NAC rtp NAC q 0.0
23 opls_238 3 NAC N 8 -0.5 14.0067
opls_238 -0.5 14.0067
24 opls_241 3 NAC H 8 0.3 1.008
opls_241 0.3 1.008
25 opls_242 3 NAC CH3 9 0.02 12.011
opls_242 0.02 12.011
26 opls_140 3 NAC HH31 9 0.06 1.008
opls_140 0.06 1.008
27 opls_140 3 NAC HH32 9 0.06 1.008
opls_140 0.06 1.008
28 opls_140 3 NAC HH33 9 0.06 1.008
opls_140 0.06 1.008
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
5 6 1
5 7 1
7 8 1
7 9 1
9 10 1
9 11 1
9 21 1
11 12 1
11 13 1
11 17 1
13 14 1
13 15 1
13 16 1
17 18 1
17 19 1
17 20 1
21 22 1
21 23 1
23 24 1
23 25 1
25 26 1
25 27 1
25 28 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
1 8 1
1 9 1
2 6 1
2 7 1
3 6 1
3 7 1
4 6 1
4 7 1
5 10 1
5 11 1
5 21 1
6 8 1
6 9 1
7 12 1
7 13 1
7 17 1
7 22 1
7 23 1
8 10 1
8 11 1
8 21 1
9 14 1
9 15 1
9 16 1
9 18 1
9 19 1
9 20 1
9 24 1
9 25 1
10 12 1
10 13 1
10 17 1
10 22 1
10 23 1
11 22 1
11 23 1
12 14 1
12 15 1
12 16 1
12 18 1
12 19 1
12 20 1
12 21 1
13 18 1
13 19 1
13 20 1
13 21 1
14 17 1
15 17 1
16 17 1
17 21 1
21 26 1
21 27 1
21 28 1
22 24 1
22 25 1
24 26 1
24 27 1
24 28 1
[ angles ]
; ai aj ak funct c0 c1
c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
1 5 6 1
1 5 7 1
6 5 7 1
5 7 8 1
5 7 9 1
8 7 9 1
7 9 10 1
7 9 11 1
7 9 21 1
10 9 11 1
10 9 21 1
11 9 21 1
9 11 12 1
9 11 13 1
9 11 17 1
12 11 13 1
12 11 17 1
13 11 17 1
11 13 14 1
11 13 15 1
11 13 16 1
14 13 15 1
14 13 16 1
15 13 16 1
11 17 18 1
11 17 19 1
11 17 20 1
18 17 19 1
18 17 20 1
19 17 20 1
9 21 22 1
9 21 23 1
22 21 23 1
21 23 24 1
21 23 25 1
24 23 25 1
23 25 26 1
23 25 27 1
23 25 28 1
26 25 27 1
26 25 28 1
27 25 28 1
[ dihedrals ]
; ai aj ak al funct c0 c1
c2 c3 c4 c5
2 1 5 6 3
2 1 5 7 3
3 1 5 6 3
3 1 5 7 3
4 1 5 6 3
4 1 5 7 3
1 5 7 8 3
1 5 7 9 3
6 5 7 8 3
6 5 7 9 3
5 7 9 10 3
5 7 9 11 3
5 7 9 21 3
8 7 9 10 3
8 7 9 11 3
8 7 9 21 3
7 9 11 13 3 dih_VAL_chi1_N_C_C_C
7 9 11 17 3 dih_VAL_chi1_N_C_C_C
21 9 11 13 3 dih_VAL_chi1_C_C_C_CO
21 9 11 17 3 dih_VAL_chi1_C_C_C_CO
7 9 11 12 3
10 9 11 12 3
10 9 11 13 3
10 9 11 17 3
21 9 11 12 3
7 9 21 22 3
7 9 21 23 3
10 9 21 22 3
10 9 21 23 3
11 9 21 22 3
11 9 21 23 3
9 11 13 14 3
9 11 13 15 3
9 11 13 16 3
12 11 13 14 3
12 11 13 15 3
12 11 13 16 3
17 11 13 14 3
17 11 13 15 3
17 11 13 16 3
9 11 17 18 3
9 11 17 19 3
9 11 17 20 3
12 11 17 18 3
12 11 17 19 3
12 11 17 20 3
13 11 17 18 3
13 11 17 19 3
13 11 17 20 3
9 21 23 24 3
9 21 23 25 3
22 21 23 24 3
22 21 23 25 3
21 23 25 26 3
21 23 25 27 3
21 23 25 28 3
24 23 25 26 3
24 23 25 27 3
24 23 25 28 3
[ dihedrals ]
; ai aj ak al funct c0 c1
c2 c3
1 7 5 6 1 improper_O_C_X_Y
5 9 7 8 1 improper_Z_N_X_Y
9 23 21 22 1 improper_O_C_X_Y
21 25 23 24 1 improper_Z_N_X_Y
; Include Position restraint file
#ifdef POSRES
#include "Ac_V_NHMe/Ac_V_NHMe.itp"
#endif
; Include water topology
#include "oplsaa.ff/tip3p.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "oplsaa.ff/ions.itp"
[ system ]
; Name
ACETYL-VALINE-METHYLAMIDE in water
[ molecules ]
; Compound #mols
Protein_chain_A 1
SOL 682
NVT_VDW-Q.mdp
here I intend to switch off the coulomb interactions.
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 10000 ; 20 ps
nstcomm = 100
; Output control
nstxout = 10000
nstvout = 10000
nstfout = 0
nstlog = 1000
nstenergy = 500
nstxtcout = 500
xtc-precision = 500
xtc_grps = Protein
; Neighborsearching and short-range nonbonded interactions
nstlist = 5
ns_type = grid
pbc = xyz
rlist = 1.0
; Electrostatics
coulombtype = PME
rcoulomb = 1.0
; van der Waals
vdw-type = switch
rvdw-switch = 0.8
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
optimize_fft = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 300
; Pressure coupling is off for NVT
Pcoupl = No
tau_p = 0.5
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda = 0.8
delta_lambda = 0
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.005
couple-moltype = Protein_chain_A ; name of moleculetype to
decouple
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = vdw ; turn off everything, in this case only
vdW
couple-intramol = no
; Generate velocities to start
gen_vel = yes
gen_temp = 300
gen_seed = -1
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
NVT_VDW.mdp
here I intend to switch off the Van der Waals interactions. The only
difference with the previous file is in the Free energy parameters, so I
will report only those.
; Free energy control stuff
free_energy = yes
init_lambda = 0.8
delta_lambda = 0
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
couple-moltype = Protein_chain_A ; name of moleculetype to
decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case
only vdW
More information about the gromacs.org_gmx-users
mailing list