Re: [AMBER-Developers] Origin of dihedral term

From: Ross Walker <ross.rosswalker.co.uk>
Date: Tue, 20 Oct 2015 07:26:52 -0700

> Hi Carlos,

What do you fit to then? The forces? Or the second derivatives?

The issue is that when you expand out the 1+cos(n.theta - phi) terms you get

ndih n[1] n[1]
sum ( sum Vn + sum Vn.cos(n.theta - phi) )
i=1 j=1 j=1

The issue here is that if your fitting algorithm varies Vn in order to fit the energy to a QM reference - even if you have an offset so you are are only fitting relative energies - so your mean energy of all your MM terms (essentially the offset from the QM energies) also changes.

I think Jason is right that one should also fit the n=0 term but I 'think', need to look more closely, that this still has the same problem. For now I am just discarding the 1+ term and fitting just to Vn.cos(n.theta - phi) under the premise that this fitted against relative energy differences yields Vn terms that one can just use in the standard force field equation given that the 1+ falls out in the derivative.

I am interested to hear what others do though.

All the best
Ross

> On Oct 20, 2015, at 03:18, Carlos Simmerling <carlos.simmerling.gmail.com> wrote:
>
> For the protein force fields we don't fit the absolute energies, so the
> offset doesn't matter.
> Carlos
> On Oct 19, 2015 7:30 PM, "Jason Swails" <jason.swails.gmail.com> wrote:
>
>> On Mon, Oct 19, 2015 at 6:25 PM, Ross Walker <rosscwalker.gmail.com>
>> wrote:
>>
>>> Hi All,
>>>
>>> I am wondering if anyone knows the origin or why the dihedral term we use
>>> in the AMBER force field has the 1+ term in it.
>>>
>>> I.e. we have:
>>>
>>> ndih n[i]
>>> sum sum( 1/2 x Vn(1+cos(n.theta - phi)) )
>>> i=1 j=1
>>>
>>> Specifically we have the 1+cos term in there which is I guess to make the
>>> cos term oscillate between 0 and 2 rather than -1 and +1 and then we have
>>> the 1/2 in front of the Vn to get rid of the 2 making it between 0 and
>> Vn.
>>> However, as far as I can tell this is purely cosmetic. Is that correct?
>>>
>>
>> ​I always thought of it more as making it *look* a little more like the
>> other bonded terms (which are harmonic, and have the 1/2 term in front of
>> the force constant, and as a result have a minimum energy at 0). But
>> obviously the 1+ has no effect on forces, and the 1/2 can be pulled into
>> the Vn term (which, in Amber, it already is).
>>
>> As in I could ditch the 1/2, and ditch the 1+ and just have
>>> V[i,j].cos(n.theta-phi). The question is if that is true why don't we do
>>> this - does anyone know?
>>>
>>
>> ​At this stage, it would be historical.​
>>
>>
>> The issue arrises not in MD but when we try and refit the torsion terms. If
>>> we try to fit energies against quantum energies we always have an offset
>> in
>>> the mean due to the origins not matching - that doesn't matter since it
>>> would be constant during an MD run. However, if we are fitting Vn terms
>> the
>>> 1+cos term here causes our mean to drift as we adjust Vn. This is a pain
>> in
>>> the butt when it comes to getting a good fit. Thus I propose to just fit:
>>> Vn.cos(n.theta-phi) which, I believe would give perfectly transferable
>>> parameters to the 1+cos equation.
>>>
>>
>> ​Lachele actually addressed this precise problem in her presentation at the
>> Amber developer's meeting last year. The solution is simple: fit the
>> zero-periodicity term. That *gives* you an arbitrary constant to improve
>> your fit -- it doesn't contribute to an overfitting problem and solves the
>> issue you're describing. Then you're free to simply throw that term away
>> when making the frcmod file since it has no effect on forces. It seems to
>> me you are rediscovering her problems :). I'm sure she'd be happy to share
>> her slides if she can find them and you wanted them. </throwing Lachele
>> under the bus>
>>
>> Am I missing anything here? - Is there any reason why we couldn't go the
>>> whole step and just ditch the 1/2 and 1+ terms from the entire force
>> field?
>>> - Seems to me all that would happen is for a given set of parameters (say
>>> FF14SB) we just shift the energy origin but everything else would still
>>> work.
>>>
>>
>> ​A couple things. If you do a QM scan, you're certainly not going to get
>> an energy of zero at the midpoint -- the zero point energy will be
>> arbitrary. The easiest thing to do with your QM scan is simply scale the
>> whole potential by the minimum energy value, which will give you the same
>> minimum as the 1+cos() series at 0 (maybe that's why it's done?).
>>
>> Y
>> ​ou can fit however you'd like, but changing how Amber computes this
>> internally would break from how every other package calculates proper
>> torsions. It would make cross-program conversion validation much harder
>> than it currently is. This is a big deal IMO.
>>
>> All the best,
>> Jason
>>
>> --
>> Jason M. Swails
>> BioMaPS,
>> Rutgers University
>> Postdoctoral Researcher
>> _______________________________________________
>> AMBER-Developers mailing list
>> AMBER-Developers.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber-developers
>>
> _______________________________________________
> AMBER-Developers mailing list
> AMBER-Developers.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber-developers


_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Tue Oct 20 2015 - 07:30:06 PDT
Custom Search