[gmx-users] Problem in Free Energy Calculations
Mark Abraham
mark.j.abraham at gmail.com
Tue Jan 5 10:20:45 CET 2016
Hi,
On Tue, Jan 5, 2016 at 9:15 AM shagun krishna <krishna.shagun123 at gmail.com>
wrote:
> Hiiii Justin,
>
> I have done the simulation of ligand in water. How ever while analyzing I
> found that the ligand comes out of the solvent box at the end of the
> simulation and during simulation too, it vibrates a lot. Please find my
> final md.mdp file and suggest me where am I doing wrong.
>
This is normal. See
http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions
PS Your PME and LINCS settings are very likely lots more expensive than
they need to be. No force field I've heard of was parameterised like that...
Mark
>
> ; Run control
> integrator = sd ; Langevin dynamics
> tinit = 0
> dt = 0.002
> nsteps = 5000000 ; 10 ns
> nstcomm = 100
> ; Output control
> nstxout = 500
> nstvout = 500
> nstfout = 0
> nstlog = 500
> nstenergy = 500
> nstxtcout = 0
> xtc-precision = 1000
> ; Neighborsearching and short-range nonbonded interactions
> nstlist = 10
> ns_type = grid
> pbc = xyz
> rlist = 1.7
> ; Electrostatics
> coulombtype = PME
> rcoulomb = 1.7
> ; van der Waals
> vdw-type = switch
> rvdw-switch = 0.8
> rvdw = 0.9
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr = EnerPres
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing = 0.12
> ; EWALD/PME/PPPM parameters
> pme_order = 6
> ewald_rtol = 1e-06
> epsilon_surface = 0
> optimize_fft = no
> ; Temperature coupling
> ; tcoupl is implicitly handled by the sd integrator
> tc_grps = system
> tau_t = 1.0
> ref_t = 300
> ; Pressure coupling is on for NPT
> Pcoupl = Parrinello-Rahman
> tau_p = 0.5
> compressibility = 4.5e-05
> ref_p = 1.0
> ; Do not generate velocities
> gen_vel = no
> ; options for bonds
> constraints = h-bonds ; 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
>
> Best regards,
>
> Shagun
>
> On Tue, Jan 5, 2016 at 12:04 PM, shagun krishna <
> krishna.shagun123 at gmail.com
> > wrote:
>
> > Hi Justin..
> > Thanks a ton to you. Finally I have completed the simulation of ligand in
> > solvent. In your earlier mail you have written that to calculate the
> > binding energy of protein-ligand I will need the nonbonded interaction
> > energy between protein-ligand coming from protein-ligand simulation and
> the
> > ligand-water interaction energy coming from the second simulation. For
> > final ΔG binding I should substract these two energies, if I am not
> wrong.
> > PLs suggest me.
> >
> > Thanx in advance..
> >
> > Best regards,
> >
> > SHAGUN
> >
> > On Mon, Jan 4, 2016 at 4:58 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
> >
> >>
> >>
> >> On 1/4/16 6:25 AM, shagun krishna wrote:
> >>
> >>> Hii Justin,
> >>>
> >>> Thanks for the suggestion. When I am performing the ligand simulation
> in
> >>> water I am getting the following error in nvt equilibration:
> >>>
> >>> Program mdrun, VERSION 4.6.5
> >>> Source code file: /build/buildd/gromacs-4.6.5/src/mdlib/domdec_top.c,
> >>> line:
> >>> 393
> >>>
> >>> Fatal error:
> >>> 56 of the 13118 bonded interactions could not be calculated because
> some
> >>> atoms involved moved further apart than the multi-body cut-off distance
> >>> (0.841804 nm) or the two-body cut-off distance (1.7 nm), see option
> -rdd,
> >>> for pairs and tabulated bonds also see option -ddcheck
> >>> For more information and tips for troubleshooting, please check the
> >>> GROMACS
> >>> website at http://www.gromacs.org/Documentation/Errors
> >>>
> >>>
> >>
> >>
> http://www.gromacs.org/Documentation/Terminology/Blowing_Up#Diagnosing_an_Unstable_System
> >>
> >> Your ligand topology is probably incorrect in some way. Turn off the
> >> free energy code, use sensible cutoffs as required by the parent force
> >> field, and follow the above link.
> >>
> >> -Justin
> >>
> >>
> >> My nvt.mdp file is as follows:
> >>>
> >>> ; Run control
> >>> integrator = sd ; Langevin dynamics
> >>> tinit = 0
> >>> dt = 0.002
> >>> nsteps = 50000 ; 100 ps
> >>> nstcomm = 100
> >>> ; Output control
> >>> nstxout = 500
> >>> nstvout = 500
> >>> nstfout = 0
> >>> nstlog = 500
> >>> nstenergy = 500
> >>> nstxtcout = 0
> >>> xtc-precision = 1000
> >>> ; Neighborsearching and short-range nonbonded interactions
> >>> nstlist = 10
> >>> ns_type = grid
> >>> pbc = xyz
> >>> rlist = 1.7
> >>> ; Electrostatics
> >>> coulombtype = PME
> >>> rcoulomb = 1.7
> >>> ; van der Waals
> >>> vdw-type = switch
> >>> rvdw-switch = 0.8
> >>> rvdw = 0.9
> >>> ; Apply long range dispersion corrections for Energy and Pressure
> >>> DispCorr = EnerPres
> >>> ; Spacing for the PME/PPPM FFT grid
> >>> fourierspacing = 0.12
> >>> ; EWALD/PME/PPPM parameters
> >>> pme_order = 6
> >>> ewald_rtol = 1e-06
> >>> epsilon_surface = 0
> >>> optimize_fft = no
> >>> ; Temperature coupling
> >>> ; tcoupl is implicitly handled by the sd integrator
> >>> tc_grps = system
> >>> tau_t = 1.0
> >>> ref_t = 300
> >>> ; Pressure coupling is off for NVT
> >>> Pcoupl = No
> >>> 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
> >>> foreign_lambda = 0.05
> >>> sc-alpha = 0.5
> >>> sc-power = 1.0
> >>> sc-sigma = 0.3
> >>> couple-moltype = 1M_R ; name of moleculetype to decouple
> >>> couple-lambda0 = vdw ; only van der Waals interactions
> >>> couple-lambda1 = none ; turn off everything, in this case
> >>> only vdW
> >>> couple-intramol = no
> >>> nstdhdl = 10
> >>> ; Generate velocities to start
> >>> gen_vel = yes
> >>> gen_temp = 300
> >>> gen_seed = -1
> >>> ; options for bonds
> >>> constraints = h-bonds ; we only have C-H bonds here
> >>> ; Type of constraint algorithm
> >>> constraint-algorithm = lincs
> >>> ; Do not constrain the starting configuration
> >>> continuation = no
> >>> ; Highest order in the expansion of the constraint coupling matrix
> >>> lincs-order = 12
> >>>
> >>>
> >>> On Mon, Jan 4, 2016 at 4:20 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
> >>>
> >>>
> >>>>
> >>>> On 1/4/16 1:09 AM, shagun krishna wrote:
> >>>>
> >>>> Hi Justin,
> >>>>> You mean to say that I don't need to follow your Free Energy
> >>>>> Calculations
> >>>>> :
> >>>>> Methane in Water.
> >>>>> Actually I want to calculate the binding energy of my ligands with my
> >>>>> protein. All I should do is to run the one simulation of my ligands
> >>>>> simply
> >>>>> in water and another is the protein-ligand complex simulation, both
> for
> >>>>> 10ns. I am really very confused because I have gone through
> literature
> >>>>> and
> >>>>> all of them have calculated the value of α, β and γ, which i guess
> >>>>> can be
> >>>>> calculated by using free energy codes.
> >>>>>
> >>>>>
> >>>>> I don't use the LIE algorithm, so I don't know all the details.
> AFAIK
> >>>> these are empirical parameters to weight various contributions to the
> >>>> interaction energy. LIE is ultimately an estimate of the free energy,
> >>>> and
> >>>> its accuracy depends on the quality of the fitted parameters.
> >>>>
> >>>> A rigorous free energy calculation is the best method for calculating
> >>>> free
> >>>> energy of binding. For instance, LIE (again, AFAIK) relies solely on
> >>>> enthalpic terms, so it does not account for entropy unless you do
> >>>> additional work to figure this out separately.
> >>>>
> >>>>
> >>>> -Justin
> >>>>
> >>>> --
> >>>> ==================================================
> >>>>
> >>>> Justin A. Lemkul, Ph.D.
> >>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
> >>>>
> >>>> Department of Pharmaceutical Sciences
> >>>> School of Pharmacy
> >>>> Health Sciences Facility II, Room 629
> >>>> University of Maryland, Baltimore
> >>>> 20 Penn St.
> >>>> Baltimore, MD 21201
> >>>>
> >>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> >>>> http://mackerell.umaryland.edu/~jalemkul
> >>>>
> >>>> ==================================================
> >>>> --
> >>>> 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.
> >>>>
> >>>>
> >> --
> >> ==================================================
> >>
> >> Justin A. Lemkul, Ph.D.
> >> Ruth L. Kirschstein NRSA Postdoctoral Fellow
> >>
> >> Department of Pharmaceutical Sciences
> >> School of Pharmacy
> >> Health Sciences Facility II, Room 629
> >> University of Maryland, Baltimore
> >> 20 Penn St.
> >> Baltimore, MD 21201
> >>
> >> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> >> http://mackerell.umaryland.edu/~jalemkul
> >>
> >> ==================================================
> >> --
> >> 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