[gmx-users] fep protocol for gromacs 2016
asaffarhi at post.tau.ac.il
asaffarhi at post.tau.ac.il
Tue Jan 9 16:23:28 CET 2018
Dear Yoav,
Hi, did you use a soft-core technique?
In which lambda range do you have large free energy differences?
Best,
Asaf
Quoting Yoav Atsmon-Raz <yoav.atsmonraz at ucalgary.ca>:
> 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
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
> or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list