CustomHbondForce in XML File?

[Mary Van Vleet]

I'm looking to use CustomHbondForce in my force fields. For consistency, I would prefer to add a CustomHbondForce tag to a pre-existing xml file that describes the other forces present in my force field.

Looking at the documentation,

http://docs.openmm.org/7.0.0/userguide/ ... e-xml-file

however, there is no example listed for how to use CustomHbondForce in the xml file. Is it possible to include this force in the xml, or are we required to use the python layer to define the functions and parameters for CustomHbondForce?

If using CustomHbondForce in the xml file is possible, an example (either on the forum or in the OpenMM documentation itself) would be very helpful.

Thanks in advance,
Mary

[Peter Eastman]

There's no XML tag for CustomHbondForce. Implementing that in a general way would be a bit more complicated than some of the other forces due to the difficulty of specifying donors and acceptors. With, say, CustomBondForce it's easy. You just specify the types of two atoms, and it adds a term for every pair of bonded atoms that have those types. But with CustomHbondForce, a donor or acceptor can involve up to three atoms that don't need to be bonded in any particular way.

That said, though, it doesn't sound like a very hard problem. For example, we could just require that atom1 be bonded to atom2, and atom2 to atom3. I can't imagine a case where the atoms wouldn't all be bonded to each other, so really all that does is restrict how you order the atoms when defining the force. So we really should add that tag.

Peter

[Mary Van Vleet]

I can't think of an exception either to the strategy you mentioned, and I imagine that requiring bonded acceptor and donor groups would allow you to use a nearly-identical set-up as CustomAngleForce (except with Donor and Acceptor tags rather than Angle tags).

Let me know if you plan to add this in the near future; otherwise, assuming only forcefield.py needs to be modified, I can probably write the CustomHbondforce xml parser myself. Either way, thanks for the help and the quick response!

[Peter Eastman]

I'm working on it now! Hopefully I'll have something to show in the next day or two.

Peter

[Peter Eastman]

I sent a pull request with the changes: https://github.com/pandegroup/openmm/pull/1593. Try it out and see how it works for you.

Peter

[Mary Van Vleet]

Peter,

I didn't see your post until just now, so sorry for the delay.

I've downloaded the changes and ran the example you listed in the documentation; that particular example works fine for me. However, I notice a couple of bugs (either that or I'm doing something wrong) when using different combinations of parameters for particlesPerDonor and particlesPerAcceptor. I've attached my tests using a benzene dimer as a trial system; a summary of the tests results is as follows:

particlesPerDonor="3" particlesPerAcceptor="3"
Works as expected.

particlesPerDonor="3" particlesPerAcceptor="2"
Works as expected.

particlesPerDonor="3" particlesPerAcceptor="1"
Program doesn't crash, but the number of acceptors is listed as zero (I anticipated 12 acceptors), and thus the energy from CustomHbondForce is zero.

particlesPerDonor="2" particlesPerAcceptor="3"
Program crashes with the following error message:

Code: Select all

  File "run_benzene.py", line 27, in <module>
    ewaldErrorTolerance=0.0005)
  File "/srv/home/mvanvleet/anaconda2/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 1220, in createSystem
    force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
  File "/srv/home/mvanvleet/anaconda2/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 2817, in createForce
    force.addDonor(bond[0], bond[1], -1, self.donorParamValues[i])
TypeError: '_BondData' object does not support indexing

particlesPerDonor="2" particlesPerAcceptor="2"
Program crashes with the same error message as for the particlesPerDonor="2" particlesPerAcceptor="3" case.

particlesPerDonor="2" particlesPerAcceptor="1"
Program crashes with the same error message as for the particlesPerDonor="2" particlesPerAcceptor="3" case.

particlesPerDonor="1" particlesPerAcceptor="3"
Program doesn't crash, but the number of donors is listed as zero (I anticipated 12 donors), and thus the energy from CustomHbondForce is zero.

particlesPerDonor="1" particlesPerAcceptor="2"
Program doesn't crash, but the number of donors is listed as zero (I anticipated 12 donors), and thus the energy from CustomHbondForce is zero.

particlesPerDonor="1" particlesPerAcceptor="1"
Program doesn't crash, but both the number of acceptors and donors are listed as zero (I anticipated 12 for each), and thus the energy from CustomHbondForce is zero.


Hope that helps. Let me know if you need more information from my end or if I can be helpful in the debugging process.

Cheers,
Mary

[Peter Eastman]

Thanks for testing it! I just sent an update with some fixes: https://github.com/pandegroup/openmm/pull/1605. See if it works now.

Peter

[Mary Van Vleet]

Peter,

All of my example test cases now run without crashing, and the number of acceptor/donor particles is correct.

I still think there's a bug in the code whereby the acceptor/donor labels don't match between the energy expression (a1, d1, etc.) and the Donor/Acceptor tag (class1, class2, etc.):

To illustrate this, I've recopied the 3 Donor/2 Acceptor and 3 Donor/3 Acceptor test cases below. Because the energy expression doesn't depend on a3/Acceptor class3,

Code: Select all

<CustomHbondForce particlesPerDonor="3" particlesPerAcceptor="3" bondCutoff="2"
           energy="scale*k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2">
        <GlobalParameter name="scale" defaultValue="1"/>
        <PerDonorParameter name="theta0"/>
        <PerAcceptorParameter name="k"/>
        <PerAcceptorParameter name="r0"/>
        <Donor class1="Hsp2" class2="Csp2" class3="Csp2" theta0="2.1"/>
        <Acceptor class1="Csp2" class2="Csp2" class3="Csp2" k="115.0" r0="0.2"/>
       </CustomHbondForce>

these two cases should give identical energies. However, I get the following output:

Code: Select all

3_donor_2_acceptor
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7fd11a02a930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7fd11a02a960> > 0.0 kJ/mol
----
3_donor_3_acceptor
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7f9998c51930> > 12200.575102 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f9998c51960> > 0.0 kJ/mol
----

which shows different energies for the two systems. If I replace all the a1 labels with a2 in the 3 Donor/3 Acceptor case, I then get consistent energies:

Code: Select all

3_donor_1_acceptor
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7f1e7ca6c930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f1e7ca6c960> > 0.0 kJ/mol
----
3_donor_2_acceptor
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7f6cd0385930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f6cd0385960> > 0.0 kJ/mol
----
3_donor_3_acceptor
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7f51cea79930> > 12200.575102 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f51cea79960> > 0.0 kJ/mol
----
3_donor_3_acceptor_switch_labels
num acceptors: 12
num donors: 24
num exclusions: 144
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7feb3676a930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7feb3676a960> > 0.0 kJ/mol
----

This behavior seems incorrect to me, and I'll look more into it tonight to see if I can find out more specifically why this inconsistency is happening (in particular, if there's a bug in the code or if I'm just running into some unexpected behavior). In the meantime, let me know what you think.

[Mary Van Vleet]

Follow-up: If a print out the parameters for each of the acceptor and donor pairs, I get the following:

Code: Select all

3_donor_2_acceptor
num acceptors: 12
Acceptor parameters (a1,a2,a3,params):
[0, 1, -1, (115.0, 0.2)]
[1, 2, -1, (115.0, 0.2)]
[2, 3, -1, (115.0, 0.2)]
[3, 4, -1, (115.0, 0.2)]
[4, 5, -1, (115.0, 0.2)]
[5, 0, -1, (115.0, 0.2)]
[12, 13, -1, (115.0, 0.2)]
[13, 14, -1, (115.0, 0.2)]
[14, 15, -1, (115.0, 0.2)]
[15, 16, -1, (115.0, 0.2)]
[16, 17, -1, (115.0, 0.2)]
[17, 12, -1, (115.0, 0.2)]
num donors: 24
Acceptor parameters (d1,d2,d3,params):
[0, 1, 7, (2.1,)]
[0, 5, 11, (2.1,)]
[1, 0, 6, (2.1,)]
[1, 2, 8, (2.1,)]
[2, 1, 7, (2.1,)]
[2, 3, 9, (2.1,)]
[3, 2, 8, (2.1,)]
[3, 4, 10, (2.1,)]
[4, 3, 9, (2.1,)]
[4, 5, 11, (2.1,)]
[5, 0, 6, (2.1,)]
[5, 4, 10, (2.1,)]
[12, 13, 19, (2.1,)]
[12, 17, 23, (2.1,)]
[13, 12, 18, (2.1,)]
[13, 14, 20, (2.1,)]
[14, 13, 19, (2.1,)]
[14, 15, 21, (2.1,)]
[15, 14, 20, (2.1,)]
[15, 16, 22, (2.1,)]
[16, 15, 21, (2.1,)]
[16, 17, 23, (2.1,)]
[17, 12, 18, (2.1,)]
[17, 16, 22, (2.1,)]
num exclusions: 144
energy expression: scale*k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7ffc35ee5930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7ffc35ee5960> > 0.0 kJ/mol
----
3_donor_3_acceptor
num acceptors: 12
Acceptor parameters (a1,a2,a3,params):
[0, 1, 2, (115.0, 0.2)]
[0, 5, 4, (115.0, 0.2)]
[1, 0, 5, (115.0, 0.2)]
[1, 2, 3, (115.0, 0.2)]
[2, 3, 4, (115.0, 0.2)]
[3, 4, 5, (115.0, 0.2)]
[12, 13, 14, (115.0, 0.2)]
[12, 17, 16, (115.0, 0.2)]
[13, 12, 17, (115.0, 0.2)]
[13, 14, 15, (115.0, 0.2)]
[14, 15, 16, (115.0, 0.2)]
[15, 16, 17, (115.0, 0.2)]
num donors: 24
Acceptor parameters (d1,d2,d3,params):
[0, 1, 7, (2.1,)]
[0, 5, 11, (2.1,)]
[1, 0, 6, (2.1,)]
[1, 2, 8, (2.1,)]
[2, 1, 7, (2.1,)]
[2, 3, 9, (2.1,)]
[3, 2, 8, (2.1,)]
[3, 4, 10, (2.1,)]
[4, 3, 9, (2.1,)]
[4, 5, 11, (2.1,)]
[5, 0, 6, (2.1,)]
[5, 4, 10, (2.1,)]
[12, 13, 19, (2.1,)]
[12, 17, 23, (2.1,)]
[13, 12, 18, (2.1,)]
[13, 14, 20, (2.1,)]
[14, 13, 19, (2.1,)]
[14, 15, 21, (2.1,)]
[15, 14, 20, (2.1,)]
[15, 16, 22, (2.1,)]
[16, 15, 21, (2.1,)]
[16, 17, 23, (2.1,)]
[17, 12, 18, (2.1,)]
[17, 16, 22, (2.1,)]
num exclusions: 144
energy expression: scale*k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7f6aed512930> > 12200.575102 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f6aed512960> > 0.0 kJ/mol
----
3_donor_3_acceptor_switch_labels
num acceptors: 12
Acceptor parameters (a1,a2,a3,params):
[0, 1, 2, (115.0, 0.2)]
[0, 5, 4, (115.0, 0.2)]
[1, 0, 5, (115.0, 0.2)]
[1, 2, 3, (115.0, 0.2)]
[2, 3, 4, (115.0, 0.2)]
[3, 4, 5, (115.0, 0.2)]
[12, 13, 14, (115.0, 0.2)]
[12, 17, 16, (115.0, 0.2)]
[13, 12, 17, (115.0, 0.2)]
[13, 14, 15, (115.0, 0.2)]
[14, 15, 16, (115.0, 0.2)]
[15, 16, 17, (115.0, 0.2)]
num donors: 24
Acceptor parameters (d1,d2,d3,params):
[0, 1, 7, (2.1,)]
[0, 5, 11, (2.1,)]
[1, 0, 6, (2.1,)]
[1, 2, 8, (2.1,)]
[2, 1, 7, (2.1,)]
[2, 3, 9, (2.1,)]
[3, 2, 8, (2.1,)]
[3, 4, 10, (2.1,)]
[4, 3, 9, (2.1,)]
[4, 5, 11, (2.1,)]
[5, 0, 6, (2.1,)]
[5, 4, 10, (2.1,)]
[12, 13, 19, (2.1,)]
[12, 17, 23, (2.1,)]
[13, 12, 18, (2.1,)]
[13, 14, 20, (2.1,)]
[14, 13, 19, (2.1,)]
[14, 15, 21, (2.1,)]
[15, 14, 20, (2.1,)]
[15, 16, 22, (2.1,)]
[16, 15, 21, (2.1,)]
[16, 17, 23, (2.1,)]
[17, 12, 18, (2.1,)]
[17, 16, 22, (2.1,)]
num exclusions: 144
energy expression: scale*k*(distance(a2,d1)-r0)^2*(angle(a2,d1,d2)-theta0)^2
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7faec108b930> > 12479.5535467 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7faec108b960> > 0.0 kJ/mol
----

Thus the difference between the various test cases lay in the way the acceptor/donor groups were defined. In practical application, and in contrast to the current test system, the energy expression should be uniquely defined with respect to the different choices that OpenMM could make when assigning atom indices to a1,a2, and a3. In my case, having different particlesPerAcceptor led to forcefield.py picking different (but valid) indices for a1/a2/a3, which led to different energies.

To summarize, I don't think there's a bug that needs fixing, and you can probably merge the changes to CustomHbondForce, but the 'non-unique index' issue is certainly something users should be careful of.

[Mary Van Vleet]

Peter,

One last question on this topic (hopefully). I've started testing use of the CustomHbondForce, and, as a sanity check, wanted to make sure that CustomNonbondedForce and CustomHbondForce were giving me the same energies for simple energy expressions. I've copied my xml files and run scripts below, but the output for a linear function (energy = A*r) is below:

Code: Select all

CustomHbondForce
num acceptors: 1
num donors: 1
num exclusions: 0
energy expression: scale*A*r;     A=Aex;     Aex=(Aexch1*Aexch2);      r = distance(a1,d1)

CustomNonbondedForce
num particles: 2
num exclusions: 0
energy expression: A*r;      A=Aex;     Aex=(Aexch1*Aexch2)

<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7fe794a92a80> > 0.0 kJ/mol
<simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7fe794dbc8a0> > 98.64615798 kJ/mol
<simtk.openmm.openmm.CustomHbondForce; proxy of <Swig Object of type 'OpenMM::CustomHbondForce *' at 0x7fe794dbc7b0> > 98.64615806 kJ/mol

Though the differences are small, it appears that r (from CustomNonbondedForce) and distance(a1,d1) (from CustomHbondForce) return slightly different values. I'm guessing this is a difference in precision levels, but wanted to check:

1. Are the energy differences in fact due to different precision levels?
2. If so, is there a way to force both CustomForce routines to use the same precision levels throughout?

This will help me with debugging as I implement a few custom force fields.

As always, thanks for your time.
Mary