[gmx-users] dihedral restrains fail in free energy calculation

[Justin Lemkul]

On 3/28/17 6:07 AM, Ahmet Yıldırım wrote:
> Dear users,
>
> I was trying to do the free energy calculation on GROMACS-2016.1 but I
> faced a problem. To apply intermolecular restraints I add the [
> intermolecular_interactions ] part at the end of complex topology file.
> First I turn off intermolecular interactions and then turn them on one by
> one. Intermolecular bond and angle restraints work together, but dihedral
> restraints fail when decoupling the ligand. That means every three
> restraints (bond, angle and dihedral) don't work together when ligand is
> decoupled from protein. The system gives lincs warning error in the NVT,
> NPT or md step for different ligands. The same systems work properly for
> the ligand bound state with the same intermolecular interactions. Could
> someone give me some pointers about
> why the dihedral restraints fail when decoupling the ligand?
>
> Some more details on simulation: I use sd integrator, AMBER99-ILDN
> force-field, tip3p water model, constraint on Hydrogen atoms, position
> restraint at NVT and NPT.
>
> The ID of atoms used for intermolecular_interactions are as follows. The
> backbone atoms of aminoacid that is the closest one to the ligand:
> 1793: C of the carboxyl group of amino acid
> 1778: C_alpha atom of amino acid
> 1776: N of amino acid
>
> The ID of heavy atoms of ligand:
> 3: O2
> 2: O1
> 1: S1
>
> Free energy control stuff used in the mdp files of the decoupled state is
> here:
> free-energy              = yes
> init-lambda              = 1
> delta-lambda             = 0
> sc-alpha                 = 0.3
> sc-power                 = 1
> sc-sigma                 = 0.25
> sc-coul                  = yes
> couple-moltype           = ligand
> couple-intramol          = no
> couple-lambda0           = vdw-q
> couple-lambda1           = none
> nstdhdl                  = 100
>
> [ intermolecular_interactions ] part in the complex topology identifies a
> suitable set of Boresch restraints. The last part of the complex.top is as
> follows:
>
> [ molecules ]
> ; Compound        #mols
> ligand               1
> protein             1
> SOL             4812
> NA                 30
> CL                 15
>
> [ intermolecular_interactions ]
> [ bonds ]
> ;    i     j  type     r0A     r1A     r2A    fcA    r0B     r1B     r2B
> fcB
>  1776     2 10     0.3198     0.3198    10.0000      0.000     0.3198
> 0.3198    10.0000   4184.000
>
> [ angle_restraints ]
> ;   ai    aj    ak    al  type    thA      fcA    multA  thB      fcB
> multB
>  1778  1776     2  1776          1   132.7383      0.000          1
> 132.7383     41.840          1
>  1776     2     1     2          1   157.5476      0.000          1
> 157.5476     41.840          1
>
> [ dihedral_restraints ]
> ;   ai    aj    ak    al  type    phiA     dphiA  fcA    phiB      dphiB
> fcB
>  1793  1778  1776     2          1    65.9452     0.0000      0.000
> 65.9452     0.0000     41.840
>  1778  1776     2     1          1    98.2017     0.0000      0.000
> 98.2017     0.0000     41.840
>  1776     2     1     3          1   179.1898     0.0000      0.000
> 179.1898     0.0000     41.840
>

It's really easy to choose "bad" atoms for restraining.  Consider the points in 
dx.doi.org/10.1021/ci300505n, and cited work by Boresch, Roux, etc. in 
determining the proper restraints.  I have found that virtual sites constructed 
in the middle of rings are very effective in accomplishing orientational 
restraints, even in very complex systems.  They run without a hitch if chosen 
judiciously.

-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

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