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

[Ahmet Yıldırım]

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

Cheers,

-- 
Ahmet Yıldırım