[gmx-users] Calculating Binding Affinity between Protein and Ligand using FEP

[Justin A. Lemkul]

sunita gupta wrote:
> Hello all
> 
> Earlier also I posted a query regarding FEP for protein ligand complex, 
> but I dint not help much.
> Again I would like to share detailed information regarding the protocol 
> I am following for calculation the binding free energy between protein 
> and ligand using FEP(free energy perturbation) method. Please correct me 
> if I am wrong anywhere....as the values of *dVpot/dlambda  
> dEkin/dlambda  dG/dl constr *are continuously coming zero(0).
> 
> I followed the gromacs tutorial of protien-ligand complex 
> (http://cinjweb.umdnj.edu/~kerrigje/pdf_files/trp_drug_tutor.pdf) for 
> preparing the coordinate and topology file for the whole system.
> For Free energy Calculation I followed the tutorial 
> (http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial) 
> for lambda value ranging from zero(0) to 1 and setup 11 independent job 
> for each lambda value for 5 ns.
> But the* dVpot/dlambda  dEkin/dlambda  dG/dl constr* values in the *.log 
> are continuously coming zero.
> Any help will be highly appreciated.
> 

What version of Gromacs are you using?  If it's 3.3.x, you need to define a 
B-state for the molecule that you're trying to decouple.  Furthermore, you 
should be decoupling charges and vdW parameters in separate steps; I believe 
that is mentioned in the tutorial material you quote above, and if not, it has 
surely been posted to this list several times.

If you're using version 4.0.x, note that the free energy options have changed 
somewhat.  Changes to the topology are not necessary, but different options must 
be set in the .mdp file.  See the manual.

-Justin

> For convenience I am also pasting my pro_constV.mdp (please let me know 
> if any parameter is wrong of missing...which is leading to such problem)
> 
> ; RUN CONTROL PARAMETERS =
> integrator               = md
> ; start time and timestep in ps =
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 5000000
> ; number of steps for center of mass motion removal =
> nstcomm                  = 100
> ; OUTPUT CONTROL OPTIONS =
> ; Output frequency for coords (x), velocities (v) and forces (f) =
> nstxout                  = 50000
> nstvout                  = 50000
> nstfout                  = 0
> ; Output frequency for energies to log file and energy file =
> nstlog                   = 500
> nstenergy                = 500
> energygrps               = protein non-protein
> ; Output frequency and precision for xtc file =
> nstxtcout                = 5000
> xtc-precision            = 1000
> ; This selects the subset of atoms for the xtc file. You can =
> ; select multiple groups. By default all atoms will be written. =
> ;xtc_grps                 =
> ; NEIGHBORSEARCHING PARAMETERS =
> ; nblist update frequency =
> nstlist                  = 10
> ; ns algorithm (simple or grid) =
> ns_type                  = grid
> ; Periodic boundary conditions: xyz or none =
> ;pbc                      = xyz
> ; nblist cut-off         =
> rlist                    = 1.0
> ;domain-decomposition     = no
> ; OPTIONS FOR ELECTROSTATICS AND VDW =
> ; Method for doing electrostatics =
> coulombtype              = pme
> ;rcoulomb-switch          = 0
> rcoulomb                 = 1.0
> ; Dielectric constant (DC) for cut-off or DC of reaction field =
> epsilon-r                = 1
> ; Method for doing Van der Waals =
> vdw-type                 = cut-off
> ; cut-off lengths        =
> ;rvdw-switch              = 0.8
> rvdw                     = 1.4
> ; Apply long range dispersion corrections for Energy and Pressure =
> DispCorr                  = EnerPres
> ; Spacing for the PME/PPPM FFT grid =
> fourierspacing           = 0.1
> ; FFT grid size, when a value is 0 fourierspacing will be used =
> fourier_nx               = 0
> fourier_ny               = 0
> fourier_nz               = 0
> ; EWALD/PME/PPPM parameters =
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
> optimize_fft             = yes
> ;restraints
> ;dihre=yes
> ;dihre-fc=1
> ;nstdihreout=1000
> ;disre=simple
> ;disre_fc=1
> ; Berendsen temperature coupling is on
> Tcoupl = berendsen
> tau_t = 0.1 0.1
> tc_grps = protein non-protein
> ref_t = 300 300
> ;OPTIONS FOR PRESSURE COUPLING
> Pcoupl                   = berendsen
> tau_p                    = 0.5
> compressibility          = 4.5e-05
> ref_p                    = 1.0
> ; Free energy control stuff
> free_energy              = yes
> init_lambda              = 0.0
> delta_lambda             = 0
> sc_alpha                 = 0.5
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> ; GENERATE VELOCITIES FOR STARTUP RUN =
> gen_vel                  = yes
> gen_temp                 = 300
> gen_seed                 = 173529
> ; OPTIONS FOR BONDS     =
> constraints              = hbonds
> ; Type of constraint algorithm =
> constraint-algorithm     = Lincs
> ; Do not constrain the start configuration =
> unconstrained-start      = no
> ; Relative tolerance of shake =
> shake-tol                = 0.0001
> ; Highest order in the expansion of the constraint coupling matrix =
> lincs-order              = 12
> ; Lincs will write a warning to the stderr if in one step a bond =
> ; rotates over more degrees than =
> lincs-warnangle          = 30
> 
> 
> 
> -- 
> Best Regards
> SUNITA GUPTA
> Member Research Team
> LeadInvent Technology
> TBIU, IIT Delhi, India
> Email- sunita at leadinvent.com <mailto:sunita at leadinvent.com>
> Ph- +9111 26581524 (Ex-6)
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================