AmoebaReferenceMultipoleForce Scaling Parameters

[MSchauperl]

I am a new Postdoc in Mike Gilson’s lab and, as @jchodera mentioned in an earlier post #2033 , I will work on the development of polarizable force fields. We would like to work with OpenMM and have already had an intensive discussion, whether or not OpenMM supports everything we need. It seems that the current version does not have all the features (on the API level) we would like to use.

Especially, we would be interested in using 1-2, 1-3, 1-4 and 1-5 universal exclusion and scaling factors for the permanent electrostatics and the induced dipoles for our calculations, instead of polarization groups. It seems that these factors (mScale, uScale, pScale, dScale) are already included in AmoebaReferenceMultipoleForce.cpp with a fixed value, but there is no way to change the factors via the API. Is it possible to make these factors changeable through the C++/Python API?

@jchodera and @davidlmobley suggested that the best way to proceed would be to post a request here and get @peastman 's opinion on this issue.

[peastman]

I'm marking this as a feature request. Currently those factors are computed automatically. Allowing them to be set manually would require changing the API, figuring out how to store and access them efficiently, and of course changing all kernels that use them.

[davidlmobley]

To be clear, @peastman , a big group of us would very much like this as it is necessary for a good chunk of science we'd like to explore within the Open Force Field initiative. What can we do to help make it easy or, if not easy, at least something that can be done in a reasonable time without too much pain and suffering?

[peastman]

Pretty much any change to the AMOEBA electrostatics code involves pain and suffering. That's its nature, I'm afraid! :/

[jchodera]

Especially, we would be interested in using 1-2, 1-3, 1-4 and 1-5 universal exclusion and scaling factors for the permanent electrostatics and the induced dipoles for our calculations, instead of polarization groups. It seems that these factors (mScale, uScale, pScale, dScale) are already included in AmoebaReferenceMultipoleForce.cpp with a fixed value, but there is no way to change the factors via the API. Is it possible to make these factors changeable through the C++/Python API?

It looks like these are indeed hard-coded constants in the Reference platform AmoebaReferenceMultipoleForce.cpp, and would just require exposing these through the API and passing them along to the force kernels.

For the CUDA kernels, these functions return the appropriate scaling factors. It shouldn't be too painful to have the kernels take these constants during their initialization and use them in these functions.

[jchodera]

But yes, in general, I can see how it's going to be painful to eventually overhaul all that code...

[peastman]

It shouldn't be too painful to have the kernels take these constants during their initialization and use them in these functions.

If I were to implement it that way, it would significantly degrade the performance. It would mean transferring more data to the kernels and storing more data in registers and/or shared memory. Implementing it in a way that doesn't hurt performance will be a big challenge.

[jchodera]

Ah! Sorry---hadn't realized that was the case.

We don't expect any of these scaling factors to be changeable after Context creation. For many of the other forces, you create the kernel code on the fly and JIT compile it when the Context is created with the constants embedded in the kernel source. Could that work here? Or is the problem that Mark didn't write the AMOEBA code with the infrastructure in place to dynamically generate kernel sources?

[peastman]

The problem is that this would be a fundamentally more complicated calculation involving more data. Right now, for any pair of atoms we just need two bits: are they in the same covalent group, and are they in the same polarization group? From those two bits, we can calculate all the different scale factors. With this change we would need four floating point numbers (128 bits) for each atom pair instead of just two bits.

[jchodera]

What about at least exposing the exiting hard-coded constants through the API? There should be no performance hit for that if the kernels are dynamically generated.

Also, I wonder if we might be able to draw @andysim into helping rework the AMOEBA kernels for CUDA and OpenCL in the longer term!

[peastman]

What about at least exposing the exiting hard-coded constants through the API?

@MSchauperl would that meet your needs? For reference, here is the code used to compute the scale factors. covalent and polarizationGroup are bit flags indicating whether the two atoms are in the same covalent group and the same polarization group. Could you achieve what you want just by changing the constants in those functions?

__device__ real computeDScaleFactor(unsigned int polarizationGroup, int index) {
    return (polarizationGroup & 1<<index ? 0 : 1);
}

__device__ float computeMScaleFactor(uint2 covalent, int index) {
    int mask = 1<<index;
    bool x = (covalent.x & mask);
    bool y = (covalent.y & mask);
    return (x ? (y ? 0.0f : 0.4f) : (y ? 0.8f : 1.0f));
}

__device__ float computePScaleFactor(uint2 covalent, unsigned int polarizationGroup, int index) {
    int mask = 1<<index;
    bool x = (covalent.x & mask);
    bool y = (covalent.y & mask);
    bool p = (polarizationGroup & mask);
    return (x && y ? 0.0f : (x && p ? 0.5f : 1.0f));
}

[andysim]

Thanks for the ping, @jchodera; this is something I'm very much interested in for sure. I think that setting the constants themselves, as @peastman suggested above, should be easy to do and would give a lot of flexibility. I'm actually just finalizing an OpenMM plugin that's tangentially related. We've been working on a version of the CHARMM Drude FF that uses multipoles (MP) instead of lone pairs, and induced dipoles (ID) instead of Drude particles. The MPID force field is therefore a fully analytic (due to the extrapolated algorithm being used for IDs) analog of the Drude FF that can be formulaically mapped from Drude with no reparameterization.

MPID starts from the AMOEBA plugin but has a few key differences. First, there's only 1 set of induced dipoles (not the p and d variants found in AMOEBA) that obey a simple topological 1-2 and 1-3 exclusion rule in their response to the permanent field. Second, only heavy atoms are polarizable, but the polarizabilities are anisotropic to compensate for the missing degrees of freedom. Third, only lone-pair sites have multipoles, but they go up to octopoles.

The code will be permissively licensed and open to contributions - let me know if you think it’d be useful for you. I'm also very much interested in helping out with the AMOEBA code.

[MSchauperl]

I would like to exclude 1-2 and 1-3 direct polarization and scale 1-4 polarization by a factor instead of using the polarization groups.

It seems to me that only exposing these factors would not allow me to do that, if I understand the code correctly. I think that the computeDScaleFactor (I assume that is the factor for direct polarization) function is "only" checking if the atoms are in the same polarization group or not and setting a bit to 0 or 1 which I guess sets if the charge on one of the atoms is inducing a dipole on the other atom.
So perfect for me would be if it is possible to include the bonded information in the computeDScaleFactor function? So more or less doing the same as in computeMScaleFactor for DScale?

[andysim]

@MSchauperl , there are essentially two sets of induced dipoles in AMOEBA: the p and d subscripted quantities. The former obey topological exclusions, as you want, to decouple the polarization from the usual bonded terms. The d dipoles obey group exclusions, to increase transferability among different conformations and chemical compounds. The polarization energy can be equivalently computed as E_p dotted with mu_d or as E_d dotted with mu_p. In my pre-caffeinated state I think that removing the polarization groups, by simply dumping everything in the same group, would yield a zero energy because E_d (and consequently mu_d) would be zero, but I could be wrong about that.

The problem you have is exactly what I faced for the MPID code I mentioned above. There are two possible solutions: 1) replace each instance of dscale in the code with pscale instead, forcing both sets of dipoles to be equivalent to mu_p 2) go through and rip out all of the mu_d terms in the energy and force routines. I chose option 2), which is equal parts laborious and cathartic. The force expression resulting from the dual dipoles is very messy so the simplified code should be a decent amount faster when it's properly optimized (still a work in progress). The plugin I'm writing may be a bit too aggressively simplified for your purposes; it assumes 1-2 and 1-3 are absent from permanent and induced terms, while everything else is present without a scale factor; induced-induced terms are not excluded at all due to Thole damping.

[peastman]

I chose option 2), which is equal parts laborious and cathartic.

I can imagine! :) This should definitely be faster. If I recall correctly, each iteration of the mutual induced dipoles involves repeating the same calculation twice, once for each set of dipoles. So removing one set should produce close to a 2x speedup.

only heavy atoms are polarizable

What's the reason for that? It seems like it would complicate the code with very little performance benefit.

Third, only lone-pair sites have multipoles, but they go up to octopoles.

Ooo, that will be fun to implement!

I'm also very much interested in helping out with the AMOEBA code.

We would love your help on anything you're inclined to do!

[davidlmobley]

@andysim - I'll shoot you an e-mail soon to see about having a call, but after talking with Frank Pickard recently about what we're trying to do in the OpenFF initiative/with respect to OpenMM here, he thought you'd be very interested in this. The larger goal is to (a) get a lot of what we need to try a VARIETY of polarizable force field options into OpenMM, and (b) include infrastructure to make it really easy to use these options AND parameterize force fields that employ them into our OpenFF initiative.

[MSchauperl]

@andysim Thanks for the clarification.
I think I still have a lack of understanding what pScale and dScale, respectively E_d and E_p are. I thought it divides direct and mutual polarization but that does not seem to be the case. Do I have two sets of induced dipoles on every atom? Is there any kind of documentation how these terms affect the energy function?
What would happen if every atom has its own polarization group?

[andysim]

@peastman you're right about the doubling of the work in the iterations, so I am expecting to see a speedup when I finally do benchmarking. I only polarize heavy atoms to be consistent with the drude FF. I think they do it more for physics (H atoms not very polarizable and fear of overfitting) reasons than for efficiency concerns. Right now the code just computes everything as if it were polarizable, with zeros in place for the H atoms; not such a terrible approach given that there's not much of a saving in the reciprocal space part anyway. The octopole code wasn't too bad to code up in quasi internal coordinates, thankfully.

@davidlmobley that sounds great! Coincidentally, I was just chatting with Frank this morning about your (and Peter's) direct chemical perception paper on bioRxiv. I stumbled across it this morning and think it's a really great idea. The OpenFF project will be incredibly helpful - perhaps my plugin will help to fill in some of the gaps for polarizable FFs; it could be useful for testing different flavors of electrostatics out because it supports both isotropic and anisotropic polarization, as well as multipoles up to octopoles. With a little bit of generalizing I can add back in the different permanent/induced exclusion rules if needed.

@MSchauperl sorry for the shorthand - E_p and E_d are the fields generated by topological- and group-based exclusion rules, respectively. Perhaps the best source of this stuff is the AMOEBA organic molecule paper. Essentially for each pair we generate the field and scale it by two different factors, resulting in two distinct fields; one scaled by dscale the other by pscale. For a 1-3 pair pscale would be zero (topologically excluded) and dscale would be 1 if they're in different groups or 0 if they're in the same group. For a hypothetical 1-6 pair in the same group, pscale would be 1 (too far apart to be excluded topologically) but dscale would be 0 (atoms in the same group don't interact). Therefore putting everything in its own group would result in the E_d (and mu_d, which is just E_d multiplied by the polarizability) being large because each atom would experience the field from all others. To my knowledge, there's no way to coerce this elaborate dual dipole scheme into a simple topological single dipole scheme like what you want. That's why I wrote a separate plugin (based on the existing AMOEBA code) rather than trying to extend the existing code. Feel free to ping me if this isn't clear - I've spent a lot of time staring at this stuff over the last few years but still can't claim to be an expert.

[davidlmobley]

Ha, @andysim - that's perfect timing, as when I get a moment to e-mail I was going to send you the preprint. @MSchauperl 's work with Mike Gilson is basically interfacing with the initiative -- trying to get a framework in place so ultimately we can try a variety of different approaches to polarizable FFs within OpenFF, while his work with Gilson will focus on one specific approach. I've suggested you and he touch bases directly as it sounds like you have the expertise we need to help us sort out what we need (and perhaps your plugin will actually provide the framework for us to try a lot of it).

Do you have any interest in coming on to our Slack where we're discussing the framework for this? Basically we're trying to sort out two things:

1. What different approaches people may reasonably want to try for polarizable FFs in the near-term, and
2. How to accommodate those within extensions of our SMIRNOFF format

Ultimately all of this should make it onto GitHub somewhere, but Slack can be a bit more convenient for batting around half-baked ideas/asking stupid questions, which has been a lot of what we're doing so far (well, speaking for myself on the `stupid questions' aspect).

[peastman]

The AMOEBA plugin code is really complicated. Partly that's just because of how it was written (not the way I would have done it), but partly it's because AMOEBA is complicated. The dual dipoles make a lot of things more convoluted than they would be otherwise. If we think most other multipole based force fields aren't going to use that, there will be a lot of value in having a separate plugin that avoids that complication but is otherwise more flexible.

Regarding the original request, possibly the thing to do is restructure the multipole calculation to work more like NonbondedForce does. That is, don't try to compute all the different scaled interactions in a single loop. Just have a single bit flag to mark any interaction that should be excluded, then compute the exceptions in a separate kernel. This would hopefully not affect performance much one way or the other. You need an extra kernel to compute the scaled interactions, but usually the number of them will be quite small (especially if your system is mostly water). And the main kernel for nonbonded interactions will hopefully get slightly faster: storing fewer bit flags will save a few registers, plus the occasional flop here and there can be eliminated. This would make it more reasonable to provide arbitrary scale factors for every interaction independently, just like NonbondedForce.addException().

[chemphys]

@MSchauperl , we also want to do something like this (see #2078 ), but is more complicated. In our case, we polarize all atoms (included 1-2, 1-3, 1-4) with different scaling factors. We currently have a plugin for OpenMM, but its eficiency is bad. It doesn't scale very well in CPU, and is relatively slower than our C++ code. I am in Francesco Paesani's lab, at UCSD. Sincerely, Peter is right. The AMOEBA plugin is complicated as hell, and minor changes to do what we need will be translated in changing the whole kernels and APIs. That's why we need some help from somebody in OpenMM if we want to do it, because I cannot definitely do it alone.

[davidlmobley]

@chemphys - sounds like you should be talking to @andysim . @mjschnie also has considerable interests in this area so hopefully the group can work together that will make something efficient that covers a lot of our major use cases.

[peastman]

Based on the description @chemphys posted in #2078, it sounds like there's an important additional requirement: he needs to be able to change the thole parameter used for particular pairs of atoms. So that also would be a configurable parameter for each exception.

Then there's the position dependent charges. That could be quite complicated depending on how it's defined.

[andysim]

@davidlmobley yes - Slack's a great medium for discussing these kinds of problems, so I'd be happy to join the conversation there. My handle's the same as my GitHub handle.

The Drude FF uses a different density for the Thole terms to that in AMOEBA. It sounds like the MB-Pol version is different still, so I'll see what I can do to add flexibility along those lines into the code I'm working on now and it's something that we should keep in mind for any general polarizable FF framework.

[chemphys]

@andysim , I am not very familiar with slack. Do you have a workspace there, or something like that that I can join? Or how can I contact you through there?

[davidlmobley]

@chemphys - we were discussing putting Andy onto the Slack connected with our open force field initiative; would be happy to put you there too if you can get me an e-mail address to send you login info.

[chemphys]

Please, @davidlmobley , go ahead. You can use my gmail email: marc.riera.1989@gmail.com
Thanks!

[davidlmobley]

Invites sent.

[chemphys]

Got it, thanks!

[mvanvleet]

A bit late to the discussion, but I wanted to mention/reiterate some of the possible polarization enhancements that would be highly useful to myself and others I know:

1. Flexible screening functions. As mentioned by @andysim in this thread and as discussed in detail by @chemphys in MB-pol Electrostatics in OpenMM #2078, polarization energies can depend significantly on one's choice of screening function. Currently, only the Amoeba-style Thole screening is implemented in OpenMM; however, other viable screening functions exist, and I can think of several groups (myself included) who are interested in studying how new and different screening functions will impact the resulting polarization energy. Thole's 1981 paper lists several possible screening functions, however this is not an exhaustive list.

2. Atom-specific Thole parameters with appropriate combination rules. For a given choice of screening function, each polarizable site requires a Thole screening parameter. There is some evidence to suggest that different atomtypes should use different Thole parameters. AmoebaMultipoleForce already allows for atomtype-specific Thole parameters, however the current "combination rule" for mixed interactions is simply to use the smaller of the two Thole parameters, which to my mind is somewhat over-simplified and arbitrary. I'm not aware of a general literature consensus as to an ideal combination rule for Thole parameters; nevertheless, based on the recommendations from this article and Eq. 18/19 in particular, one might surmise that a geometric combination rule would be most appropriate. Regardless, flexibility regarding combination rules for Thole parameters (and/or a better default combination rule) would be much appreciated.

3. Anisotropic polarization. As @andysim suggested, and particularly for systems with lone pairs and/or where H-atom polarizability is ignored, allowing for anisotropic induced dipole parameters would be highly useful, especially if this does not greatly increase the computational expense relative to isotropic parameters.

4. Higher-order induced multipoles. This has already been mentioned by several others, but I'll note that select force fields for water and other small molecules already take advantage of quadrupole polarizabilities. Including this functionality in OpenMM would allow these models and others to be simulated and tested.

Finally, thanks to everyone that has already contributed to this initiative. I'm very excited to see where this goes, and happy to help out however would be useful.

[aehogan]

@jchodera I think a flexible induced dipole polarizability force would be a great benefit to the community. Usage of AmoebaMultipoleForce for all induced dipole models leads to headaches because of how different the code is to the rest of OpenMM's forces. Many of the additions are already detailed above but I will point out that ripping out the dual dipoles can lead to performance / maintainability improvements. More importantly these options would allow for experimentation with polarizable forcefields options and models like MB-pol, POLIR, or our models to be simulated with OpenMM easily.