[gmx-users] Dihedral parameters from QM and MM relaxed surface scans

Mon Apr 30 03:09:26 CEST 2018

On 4/29/18 7:16 PM, Eric Smoll wrote:
> Justin,
>
> I see. Unless there are other problems impacting my observations (e.g, a
> need for double-precision GROMACS for energy minimization, some sort of
> QM-software dependent implementation of the QM method, etc.), the QM
> dihedral PES does not match the unrestrained MM PES from the literature.

In that case, I'd suggest contacting the corresponding author and asking 
for inputs. Some programs, like CHARMM, are very easy to use in this 
context. Others, like GROMACS, require the user to jump through a few 
additional hoops to get to the same place, and it's easy to become 
inconsistent in the approach.

-Justin

> Thank you for the guidance.
>
> Best,
> Eric
>
> On Sun, Apr 29, 2018 at 4:05 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>> On 4/29/18 5:48 PM, Eric Smoll wrote:
>>
>>> Mark,
>>>
>>> *Thank you for the reply! *
>>>
>>> I understand your conclusion but I am trying to rule out workflow errors
>>> on
>>> my part before I make decisions based on the suitability of this
>>> forcefield. Let's say I am testing a QM-geometry optimized molecule with
>>> two identical (same parameters) and coupled dihedrals: phiA and phiB.
>>>
>>> Across a QM dihedral relaxed surface scan (RSS), phiA is held at the
>>> initial QM geometry optimized value and phiB is incremented across the 360
>>> degrees of its rotation.
>>> Thus, each step of the QM RSS has two dihedral restraints.
>>>
>>> Then, each step of the associated MM RSS should use the same two dihedral
>>> restraints relaxing all other DOF's. The results can be used to calculate
>>> a
>>> "restrained" MM PES for phiB.
>>>
>>> In this example, how should the "unrestrained" PES be computed with "mdrun
>>> -rerun"? Should both dihedral restraints be removed across each step of
>>> the
>>> MM PES scan? Or should the dihedral restrain on phiA stay?
>>>
>> What you need to demonstrate is that the potential energy of the
>> configuration using the MM energy is equivalent to that of the QM
>> calculation. That requires the removal of all restraints/biases when
>> evaluating the energy.
>>
>> -Justin
>>
>>
>> "If the normal FF was parameterized for the species and configurations of
>>> interest, then you'd observe that the unrestrained PES matched the
>>> restrained PES, ie that the contributions from the restraint were
>>> negligible."
>>>
>>> I see. I did not make that connection. Thank you.
>>>
>>> Best,
>>> Eric
>>>
>>> On Sun, Apr 29, 2018 at 3:13 PM, Mark Abraham <mark.j.abraham at gmail.com>
>>> wrote:
>>>
>>> Hi,
>>>> Your two PES scans are over the same coordinates, and by construction you
>>>> are balancing the normal FF contributions against an intended-as-strong
>>>> restraint. The only thing that is different between your two sets of PE
>>>> calculation is the contribution of the restraint. If the normal FF was
>>>> parameterized for the species and configurations of interest, then you'd
>>>> observe that the unrestrained PES matched the restrained PES, ie that the
>>>> contributions from the restraint were negligible. Clearly that isn't the
>>>> case, so the only available conclusion is that the normal FF isn't
>>>> suitable
>>>> for modelling this PES. Whether that is a useful conclusion about the
>>>> chances for quality observations of quantities that depend on the FES is
>>>> less clear.
>>>>
>>>> Mark
>>>>
>>>> On Sat, Apr 28, 2018 at 9:30 AM Eric Smoll <ericsmoll at gmail.com> wrote:
>>>>
>>>> Hello Gromacs Users,
>>>>> I am trying to reproduce some published forcefield development work. In
>>>>> this paper, dihedral parameters are derived from QM and MM relaxed
>>>>>
>>>> surface
>>>>
>>>>> scan data collected on individual small molecules in vacuum.
>>>>>
>>>>> My replication workflow is as follows. The MM relaxed surface scan
>>>>> trajectory was prepared by concatenating the final frame of a set of
>>>>>
>>>> energy
>>>>
>>>>> minimization calculations. Each energy minimization held the dihedral of
>>>>> interest at a desired angle with a 1000 kJ/mol dihedral restraint. Other
>>>>> DOF's were allowed to relax with the "steep" algorithm in
>>>>>
>>>> single-precision
>>>>
>>>>> GROMACS 2018.1 to machine precision. A "restrained" dihedral PES was
>>>>> calculated by collecting the molecule's potential as a function of
>>>>>
>>>> dihedral
>>>>
>>>>> angle *with the restraint on the dihedral of interest*.
>>>>>
>>>>> An "unrestrained" dihedral PES was also calculated by removing the
>>>>>
>>>> dihedral
>>>>
>>>>> restraint on the dihedral of interest and calculating the energy with
>>>>> the
>>>>> rerun feature of mdrun.
>>>>>
>>>>> I find that the "restrained" MM dihedral PES is in excellent
>>>>> quantitative
>>>>> agreement with QM relaxed surface scan data. The "unrestrained" MM
>>>>>
>>>> dihedral
>>>>
>>>>> PES is only in qualitative agreement with the QM results. Peaks and dips
>>>>>
>>>> in
>>>>
>>>>> the "unrestrained" PES are exaggerated when compared to the "restrained"
>>>>> PES.
>>>>>
>>>>> As far as I understand, the "unrestrained" PES should be more
>>>>> appropriate/accurate than the "restrained" PES. Am I correct? Is there a
>>>>> flaw in my workflow?
>>>>>
>>>>> Best,
>>>>> Eric
>>>>> --
>>>>> 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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Assistant Professor
>> Virginia Tech Department of Biochemistry
>>
>> 303 Engel Hall
>> 340 West Campus Dr.
>> Blacksburg, VA 24061
>>
>> jalemkul at vt.edu | (540) 231-3129
>> http://www.thelemkullab.com
>>
>> ==================================================
>>
>>
>> --
>> 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.
>>
>
> On Sun, Apr 29, 2018 at 4:05 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>> On 4/29/18 5:48 PM, Eric Smoll wrote:
>>
>>> Mark,
>>>
>>> *Thank you for the reply! *
>>>
>>> I understand your conclusion but I am trying to rule out workflow errors
>>> on
>>> my part before I make decisions based on the suitability of this
>>> forcefield. Let's say I am testing a QM-geometry optimized molecule with
>>> two identical (same parameters) and coupled dihedrals: phiA and phiB.
>>>
>>> Across a QM dihedral relaxed surface scan (RSS), phiA is held at the
>>> initial QM geometry optimized value and phiB is incremented across the 360
>>> degrees of its rotation.
>>> Thus, each step of the QM RSS has two dihedral restraints.
>>>
>>> Then, each step of the associated MM RSS should use the same two dihedral
>>> restraints relaxing all other DOF's. The results can be used to calculate
>>> a
>>> "restrained" MM PES for phiB.
>>>
>>> In this example, how should the "unrestrained" PES be computed with "mdrun
>>> -rerun"? Should both dihedral restraints be removed across each step of
>>> the
>>> MM PES scan? Or should the dihedral restrain on phiA stay?
>>>
>> What you need to demonstrate is that the potential energy of the
>> configuration using the MM energy is equivalent to that of the QM
>> calculation. That requires the removal of all restraints/biases when
>> evaluating the energy.
>>
>> -Justin
>>
>>
>> "If the normal FF was parameterized for the species and configurations of
>>> interest, then you'd observe that the unrestrained PES matched the
>>> restrained PES, ie that the contributions from the restraint were
>>> negligible."
>>>
>>> I see. I did not make that connection. Thank you.
>>>
>>> Best,
>>> Eric
>>>
>>> On Sun, Apr 29, 2018 at 3:13 PM, Mark Abraham <mark.j.abraham at gmail.com>
>>> wrote:
>>>
>>> Hi,
>>>> Your two PES scans are over the same coordinates, and by construction you
>>>> are balancing the normal FF contributions against an intended-as-strong
>>>> restraint. The only thing that is different between your two sets of PE
>>>> calculation is the contribution of the restraint. If the normal FF was
>>>> parameterized for the species and configurations of interest, then you'd
>>>> observe that the unrestrained PES matched the restrained PES, ie that the
>>>> contributions from the restraint were negligible. Clearly that isn't the
>>>> case, so the only available conclusion is that the normal FF isn't
>>>> suitable
>>>> for modelling this PES. Whether that is a useful conclusion about the
>>>> chances for quality observations of quantities that depend on the FES is
>>>> less clear.
>>>>
>>>> Mark
>>>>
>>>> On Sat, Apr 28, 2018 at 9:30 AM Eric Smoll <ericsmoll at gmail.com> wrote:
>>>>
>>>> Hello Gromacs Users,
>>>>> I am trying to reproduce some published forcefield development work. In
>>>>> this paper, dihedral parameters are derived from QM and MM relaxed
>>>>>
>>>> surface
>>>>
>>>>> scan data collected on individual small molecules in vacuum.
>>>>>
>>>>> My replication workflow is as follows. The MM relaxed surface scan
>>>>> trajectory was prepared by concatenating the final frame of a set of
>>>>>
>>>> energy
>>>>
>>>>> minimization calculations. Each energy minimization held the dihedral of
>>>>> interest at a desired angle with a 1000 kJ/mol dihedral restraint. Other
>>>>> DOF's were allowed to relax with the "steep" algorithm in
>>>>>
>>>> single-precision
>>>>
>>>>> GROMACS 2018.1 to machine precision. A "restrained" dihedral PES was
>>>>> calculated by collecting the molecule's potential as a function of
>>>>>
>>>> dihedral
>>>>
>>>>> angle *with the restraint on the dihedral of interest*.
>>>>>
>>>>> An "unrestrained" dihedral PES was also calculated by removing the
>>>>>
>>>> dihedral
>>>>
>>>>> restraint on the dihedral of interest and calculating the energy with
>>>>> the
>>>>> rerun feature of mdrun.
>>>>>
>>>>> I find that the "restrained" MM dihedral PES is in excellent
>>>>> quantitative
>>>>> agreement with QM relaxed surface scan data. The "unrestrained" MM
>>>>>
>>>> dihedral
>>>>
>>>>> PES is only in qualitative agreement with the QM results. Peaks and dips
>>>>>
>>>> in
>>>>
>>>>> the "unrestrained" PES are exaggerated when compared to the "restrained"
>>>>> PES.
>>>>>
>>>>> As far as I understand, the "unrestrained" PES should be more
>>>>> appropriate/accurate than the "restrained" PES. Am I correct? Is there a
>>>>> flaw in my workflow?
>>>>>
>>>>> Best,
>>>>> Eric
>>>>> --
>>>>> 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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Assistant Professor
>> Virginia Tech Department of Biochemistry
>>
>> 303 Engel Hall
>> 340 West Campus Dr.
>> Blacksburg, VA 24061
>>
>> jalemkul at vt.edu | (540) 231-3129
>> http://www.thelemkullab.com
>>
>> ==================================================
>>
>>
>> --
>> 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.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

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