normal free energy perturbation mdp

  • atsmonraz
  • atsmonraz's Avatar Topic Author
  • Visitor
6 years 2 months ago #7484 by atsmonraz
free energy perturbation mdp was created by atsmonraz
Hi everyone,

I’ve been trying to get 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, Yoav

Please Log in or Create an account to join the conversation.

More
6 years 2 months ago #7486 by riccardo
Replied by riccardo on topic free energy perturbation mdp
Just to make sure I understand what you are trying to do. You have a dipeptide and you want to perturb one of its residues, LEU, into ALA - and get the free energy difference. Is that correct? Is the dipeptide then just solvated in water?

Please Log in or Create an account to join the conversation.

  • atsmonraz
  • atsmonraz's Avatar Topic Author
  • Visitor
6 years 2 months ago #7487 by atsmonraz
Replied by atsmonraz on topic free energy perturbation mdp
i actually have a pentapeptide and i'm trying to find the free energy difference in bulk water for perturbing a single residue. In the current system, the peptide is indeed solvated just in water.

Please Log in or Create an account to join the conversation.

More
6 years 2 months ago #7493 by bart
Replied by bart on topic free energy perturbation mdp
I know it is a nasty question, but are you sure your simulations somehow reached convergence? I couldn't tell you without testing some things, but I could imagine a pentapeptide does not sample all its configurations in water within 1 mus. The value indeed seems rather high in my perception, however I have never performed such an alchemical transformation.

Cheers,

Bart

Please Log in or Create an account to join the conversation.

More
6 years 2 months ago #7494 by riccardo
Replied by riccardo on topic free energy perturbation mdp
I don't see any trivial mistake in the mdp.
Besides what Bart suggested, you could also try to reproduce some published free energy difference for a similar perturbation. I can't think of any paper with a similar system right now but I'm sure there must be quite a few in the literature. Like that you'd make sure your procedure is correct.

Please Log in or Create an account to join the conversation.

  • atsmonraz
  • atsmonraz's Avatar Topic Author
  • Visitor
6 years 2 months ago - 6 years 2 months ago #7506 by atsmonraz
Replied by atsmonraz on topic free energy perturbation mdp
Thanks for your responses bart and riccardo , I'm actually following up on the protocol that was described by Singh and Tieleman in JCTC 2011, 7, 2316-2324
In that work, they actually used a 0.5 us/lambda while i'm using 1us/lambda with a preliminary 100ns NPT and 100ns NVT equilibration steps. However, I do indeed have some concerns regarding the convergence, although that is not the main issue for me right now. As I've mentioned in my post, my biggest concern for the time being is if the mdp is actually doing what its intended for. From what I understand with couple-lamdba0=vdw and couple-lambda1=none i'm actually phasing out all vdw interactions of state 0, so by state 1, their are no vdw interactions for Protein_A so I'm now trying to run the same script only with couple-lambda0=vdw-q and couple-lambda1=vdw-q so it will transition from state 0 to 1 while accounting for the interactions of the two states while exchanging between them.

Thanks in advance
Last edit: 6 years 2 months ago by atsmonraz.

Please Log in or Create an account to join the conversation.

Time to create page: 0.104 seconds