[gmx-users] fep protocol for gromacs 2016
Yoav Atsmon-Raz
yoav.atsmonraz at ucalgary.ca
Tue Jan 9 06:19:40 CET 2018
Hi everyone,
I’ve been trying to get a the free energy difference between a leucine residue vs an alanine residue by doing a set of FEP simulations in which I decouple the VDW and electrostatic interactions in the martini CG force field. In the particular example I’m using here, I’m aiming to decouple the Leucine sidechain and the topology I’m using is as follows:
--------------------------------------
[ moleculetype ]
; Name Exclusions
Protein_A 1
[ atoms ]
1 P5 1 TRP BB 1 0.0000 P5 0 72
2 SC4 1 TRP SC1 2 0.0000 SC4 0 72
3 SNd 1 TRP SC2 3 0.0000 SNd 0 72
4 SC5 1 TRP SC3 4 0.0000 SC5 0 72
5 SC5 1 TRP SC4 5 0.0000 SC5 0 72
6 Nda 2 LEU BB 6 0.0000 Nda 0 72
7 AC1 2 LEU SC1 7 0.0000 AC1 0 72
8 Nda 3 LEU BB 8 0.0000 Nda 0 72
9 AC1 3 LEU SC1 9 0.0000 DUM 0 72
10 Nda 4 LEU BB 10 0.0000 Nda 0 72
11 AC1 4 LEU SC1 11 0.0000 AC1 0 72
12 Qa 5 LEU BB 12 -1.0000 Qa -1 72
13 AC1 5 LEU SC1 13 0.0000 AC1 0 72
[ bonds ]
; Backbone bonds
1 6 1 0.35000 1250 ; TRP(E)-LEU(E)
6 8 1 0.35000 1250 ; LEU(E)-LEU(E)
8 10 1 0.35000 1250 ; LEU(E)-LEU(E)
10 12 1 0.35000 1250 ; LEU(E)-LEU(E)
; Sidechain bonds
1 2 1 0.30000 5000 ; TRP
6 7 1 0.33000 7500 ; LEU
8 9 1 0.33000 7500 ; LEU
10 11 1 0.33000 7500 ; LEU
12 13 1 0.33000 7500 ; LEU
; Short elastic bonds for extended regions
1 8 1 0.64000 2500 ; TRP(1)-LEU(8) 1-3
6 10 1 0.64000 2500 ; LEU(6)-LEU(10) 2-4
8 12 1 0.64000 2500 ; LEU(8)-LEU(12) 2-4
; Long elastic bonds for extended regions
1 10 1 0.97000 2500 ; TRP(1)-LEU(10) 1-4
6 12 1 0.97000 2500 ; LEU(6)-LEU(12) 1-4
[ constraints ]
2 3 1 0.27000 ; TRP
2 4 1 0.27000 ; TRP
3 4 1 0.27000 ; TRP
3 5 1 0.27000 ; TRP
4 5 1 0.27000 ; TRP
[ angles ]
; Backbone angles
1 6 8 2 134 25 ; TRP(E)-LEU(E)-LEU(E)
6 8 10 2 134 25 ; LEU(E)-LEU(E)-LEU(E)
8 10 12 2 134 25 ; LEU(E)-LEU(E)-LEU(E)
; Backbone-sidechain angles
2 1 6 2 100 25 ; TRP(E)-LEU(E) SBB
1 6 7 2 100 25 ; TRP(E)-LEU(E) SBB
6 8 9 2 100 25 ; LEU(E)-LEU(E) SBB
8 10 11 2 100 25 ; LEU(E)-LEU(E) SBB
10 12 13 2 100 25 ; LEU(E)-LEU(E) SBB
; Sidechain angles
1 2 3 2 210 50 ; TRP
1 2 4 2 90 50 ; TRP
[ dihedrals ]
; Backbone dihedrals
; Sidechain improper dihedrals
1 3 4 2 2 0 50 ; TRP
2 3 5 4 2 0 200 ; TRP
-------------------------------------------
The mdp I’m using goes like this:
-----------------------------------------
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.02
nsteps = 50000000 ; 1 us
nstcomm = 100
comm-grps = Protein_A Water
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxout-compressed = 1000
compressed-x-precision = 10
compressed-x-grps = Protein_A Water
energygrps = Protein_A Water
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
verlet-buffer-tolerance = 0.005
nstlist = 20
ns_type = grid
pbc = xyz
; Electrostatics
coulombtype = reaction-field
rcoulomb = 1.1
epsilon_r = 15
epsilon_rf = 0
; van der Waals
vdwtype = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = Protein_A Water
tau_t = 1.0 1.0
ref_t = 300 300
; Pressure coupling is on for NPT
Pcoupl = Berendsen
Pcoupltype = isotropic
tau_p = 5.0
compressibility = 3e-04 3e-04
ref_p = 1.0 1.0
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = -1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.20 0.40 0.60 0.80 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = Protein_A ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; van der Waals interactions and coul
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = yes
nstdhdl = 100
; Do not generate velocities
gen_vel = no
gen_temp = 300
gen_seed = 473529
; options for bonds
constraints = none ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
Can someone just tell me if I’m doing this correctly? The values I’m getting via MBAR in alchemical analysis seem to high for me (~170 kj/mol) and I want to be sure that its currently indeed addressing the switch from topology A to the mutant topology B. Also I’m using GMX 2016.3
Thanks in advance for everyone, Yoav
More information about the gromacs.org_gmx-users
mailing list