[gmx-users] fep protocol for gromacs 2016
Mark Abraham
mark.j.abraham at gmail.com
Tue Jan 9 16:35:32 CET 2018
Hi,
You need to have VDW on any atom that has a charge, or opposite charges are
at risk of collison. I think your lambda schedule is the wrong way around.
Mark
On Tue, Jan 9, 2018 at 4:23 PM <asaffarhi at post.tau.ac.il> wrote:
> 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.
>
>
>
> --
> 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