[gmx-users] dihedral MM scan

[Justin Lemkul]

On 1/18/17 9:51 PM, Mohsen Ramezanpour wrote:
> Dear Gromacs users,
>
> I want to do a MM scan on a dihedral and compare with the available QM
> profile using Gromacs 5.1.3
>
> I used the following (updated compared to my last email on parameters)
> em.mdp and md.mdp files for relaxation and zero-step md, respectively.
>
> *em.mdp:*
> define                     = -DPOSRESDIH
> integrator                 = steep
> emtol                       = 100.0
> nsteps                     = 50000
> vdwtype                   = Cut-off
> vdw-modifier            = none
> coulombtype             = Cut-off
> constraints                = h-bonds
> constraint_algorithm    = LINCS
> nstlist                        = 0
> ns-type                      = simple
> pbc                           = no
> cutoff-scheme           = group
> rcoulomb                   = 0
> rvdw                         = 0
> rvdw_switch             = 0
> rlist                            = 0
>
> Which, by -DPOSRES I am applying position restraint on specific dihedral of
> interest.
> The following lines were included at the end of molecule.itp
>
> #ifdef POSRES
>
> [ dihedral_restraints ]
>
>  1    18  8   9     1   120   0   10000
>
> #endif
>
>
>
> *as for md.mdp:*
> define                  =
> integrator              = md
> dt                        = 0.001
> nsteps                  = 0
> nstlog                  = 0
> nstxout                 = 0
> nstvout                 = 0
> nstfout                  = 0
> nstcalcenergy           = 0
> nstenergy                 = 0
> nstxout-compressed      = 100000
> compressed-x-precision  = 1000
> coulombtype                  = Cut-off
> vdwtype                         = Cut-off
> vdw-modifier                   = none
> constraints                     = h-bonds
> constraint_algorithm        = LINCS
> nstcomm                         = 100
> comm_mode               = linear
> comm_grps                 = system
> refcoord_scaling           = com
> nstlist                          = 0
> ns-type                        = simple
> pbc                            = no
> cutoff-scheme              = group
> rcoulomb                = 0
> rvdw                      = 0
> rvdw_switch            = 0
> rlist                         = 0
> continuation            = yes
>
> It seems to work fine, at least because grompp do not complain and it gives
> me the desired files.
>
> There is one problem:
> Looking at the Potential Energy which zero-step md for angle = 120 gives me
> a value of 384.682 KJ/mol (ca. 91 Kcal/mol) which is almost 20 times the
> value from QM profile (4.5 Kcal/mol) by GAAMP server (as I understood, the
> units for dihedral profiles by GAAMP are in Kcal/mol).
>
> I used the parameters GAAMP gave me (partial charges and dihedral
> parameters).
>
> Briefly, the same molecule which I uploaded to GAAMP was used for MM scan
> in vaccume. No box and infinite cutoff was used. I did an EM step by above
> em.mdp and let the system to relax. Then, I removed the restraints and did
> the zero-step md part and got this energy.
>
> It is worth mentioning that the molecule is positively charges. BUT, as I
> understood, this is what GAAMP also uses for the QM scan.
>
> What could be the problem?
>
> I highly appreciate your comments as expert in this.
>

All potential energy scans need to be offset to the global minimum of the scan. 
You can't take the energy of one configuration and compare its value to the QM. 
You have to do the whole scan and offset it.  The relative energy is what matters.

-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

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