Issues with xml forcefield file for AmoebaMultipoleForce

[Mary Van Vleet]

As detailed below, I have three issues (two coding-related, one science-related) regarding the implementation of the AmoebaMultipoleForce and its xml force field interface. For context, I'm trying to develop force field that contains both induced dipoles and permanent multipoles, and am currently in the process of testing the electrostatics of the model on a series of water trimers. I included some of my input files (which can be run by executing "bash run_all.sh") in case it's helpful.

#1: Possible bug in AmoebaMultipoleForce xml file


I'd like to test the effect of 1-4 scaling factors on my force field, as the defaults from AMOEBA may not be the best accuracy for my model. As far as I can tell, the permanent multipole and induced dipole scaling factors should be in the header of the xml file:

Code: Select all

<AmoebaMultipoleForce  direct11Scale="0.0"  direct12Scale="1.0"  direct13Scale="1.0"  direct14Scale="1.0"  mpole12Scale="0.0"  mpole13Scale="0.0"  mpole14Scale="0.4"  mpole15Scale="0.8"  mutual11Scale="1.0"  mutual12Scale="1.0"  mutual13Scale="1.0"  mutual14Scale="1.0"  polar12Scale="0.0"  polar13Scale="0.0"  polar14Intra="0.5"  polar14Scale="1.0"  polar15Scale="1.0"  >

However, changing these scale factors in the xml file seems to have no effect on the resulting trimer energies. My C++ skills and knowledge of the OpenMM code base is rather poor, but based on this file it looks like the scale factors are hard-coded in. Assuming this is true, is it possible to fix the AmoebaMultipoleForce xml interface such that changing the xml tags changes the scale factors used in computing energies and forces?

#2: Clarification question on polarization groups

I parameterized the induced dipoles to already account for intramolecular polarization, and want to make sure that I choose the correct scale factors and polarization groups such that, for each water molecule, the permanent multipole moments on the molecule *do not* polarize induced dipoles on the same water molecule, but *do* polarize induced dipoles on other water molecules. If I'm understanding the language of the AMOEBA2013 model correctly, (doi:10.1021/ct4003702) this means that both the hydrogen and oxygen atomtypes should be in the same polarization group, and I wanted to check (since there's no documentation on this in the OpenMM manual) that I achieve this by setting up my polarization groups like in the amoeba2013.xml file (type 247 is oxygen on water, type 248 is hydrogen on water):

Code: Select all

  <Polarize type="247" polarizability="0.000837" thole="0.3900" pgrp1="248" />
  <Polarize type="248" polarizability="0.000496" thole="0.3900" pgrp1="247" />
  
   <Polarize type="85" polarizability="0.00175" thole="0.3900" pgrp1="83" pgrp2="86" pgrp3="87" />
  <Polarize type="86" polarizability="0.000696" thole="0.3900" pgrp1="85" />
  <Polarize type="87" polarizability="0.00175" thole="0.3900" pgrp1="85" pgrp2="88" />

#3: Issue with the Thole damping factors


This one's part coding, part science question, and I'll try and ask the AMOEBA folks directly about any parts where the OpenMM team can't help. In line with the findings of the recent AMOEBA papers, the 2-body water energies are supposed to be relatively insensitive to the choice of Thole damping factors, whereas the choice of damping factor should matter for the many-body (3-body and up) energies. In manually testing a variety of Thole damping factors (from 0.01 to 1.0), however, I see that the 3-body energies can vary by over 100 kJ/mol, which seems highly unusual to me. This leads me to the following questions:

1. Is this magnitude of variance in the 3-body energies (as a function of Thole screening factor) normal??
2. Is there a way to test (preferably with minimal changes to the OpenMM code) whether or not a polarization catastrophe is occurring for a particular choice of damping factor?
3. Is there a known upper limit to the damping factor that will ensure (for normal molecular systems) no polarization catastrophe occurs?
4. (I'll probably have to ask the AMOEBA people this) Assuming that the many-body energies are this sensitive to the Thole damping factors, and other than comparing against benchmark 3-body energies, is there a way to determine optimal Thole damping factors for a particular system?


As always, many thanks to the OpenMM team for your help on these issues.

Cheers,
Mary

[Peter Eastman]

based on this file it looks like the scale factors are hard-coded in.

Correct, the scale factors are hard coded. This isn't a bug; it's a feature that doesn't exist. I just checked the documentation and didn't see anything suggesting they could be modified?

<Polarize type="247" polarizability="0.000837" thole="0.3900" pgrp1="248" />
<Polarize type="248" polarizability="0.000496" thole="0.3900" pgrp1="247" />

That's correct. In fact, those are precisely the lines that appear in amoeba2013.xml.

#3: Issue with the Thole damping factors

I'm afraid these questions are beyond my knowledge. You can call getInducedDipoles() on the AmoebaMultipoleForce to make sure it's producing reasonable values, but I'm not sure what the expected effect is from changing the dampling factor.

Peter