# Re: [AMBER-Developers] Origin of dihedral term

Date: Tue, 20 Oct 2015 13:27:32 -0400

Hi,

The extra sum_Vn is accounting for the fact that the ABSOLUTE value of
the QM energy is orders of magnitude larger than the MM energy, simply
duew to where our zeros of energy are (nuclei + electron separated for
QM, vs strain for MM).

We do not make a big deal when we note that our TOTAL MM energy is not
comparable to the TOTAL QM energy, so why are we making it a big deal in
dihedrals? We should not...

Now, one could make three possible choices for that shift (sum V_n)

1. Fit it as an extra parameter. This means fit the individual V_i's and
also a V_0, straight out.
2. Fit ONLY the V_is and DEFINE V_0 as sum (V_i), but do not fit it
separately.
3. Pre-fit V_0 as the shift between the average QM and the average MM
data (other choices also work) and then go to points 1 or 2 with smaller

For a complete sample of QM points, these are exactly equivalent. For an
incomplete sample they are not the same, but again, who cares ?

The forces of course are identical, and all delta Es also are.

On 10/20/15 1:06 PM, Carlos Simmerling wrote:
> Hi Ross,
> unless I'm misunderstanding, I'm not sure how this would relate to fitting
> relative energies. whatever the offset between QM and MM surfaces, it
> cancels when using the delta delta E. So we don't fit the n=0 term, but we
> also don't include it in the data, so it doesn't seem to matter.
> carlos
>
> On Tue, Oct 20, 2015 at 10:26 AM, Ross Walker <ross.rosswalker.co.uk> wrote:
>
>>> 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
>>
> _______________________________________________
> AMBER-Developers mailing list
> AMBER-Developers.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber-developers

```--