[gmx-users] Problem in Free Energy Calculations
shagun krishna
krishna.shagun123 at gmail.com
Tue Jan 5 09:14:53 CET 2016
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.
; 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.
>>
>
>
More information about the gromacs.org_gmx-users
mailing list