PME with CustomNonbondedForce and 14scale fcator

[CHANG YUN SON]

Hi,

Now I'm trying to implement a force field developed by our group into openmm.
This requires quite many different functions.
1. Drude polarization
2. Custom nonbonded interaction with 1-4 interaction scaling factor(0.5)
3. PME for electrostatics and drude

I figured out how to use the drude plugin in openmm (thanks to Peter).
But so far, I could not find out simple and efficient way to implement those task # 2 and 3.

One naive idea I can think of is a redundant definition of nonbonded interaction.
1. define regular NonbondedForce with PME and LJ, then set all sigma and epsilon to be 0.
2. define a Custom Nonbonded Interaction which is scaled by .5 with bondCutOff set to be 2
3. define a Custom Nonbonded Interaction which is scaled by .5 with bondCutOff set to be 3

In this way I think openmm will be able to calculate the nonbonded interactions as I want to be.
Though I feel this is too redundant and inefficient, and wasting lots of computation time to do either trivial calculation of 0 LJ or doing exact same force calculation twice.

Would there be better way to implement this in openmm?
Should I go into the source codes and develop my own plugin for this?

ChangYun

[Peter Eastman]

Hi ChangYun,

1-4 interactions are basically a bonded force. It's a fixed list of atom pairs chosen based on the bonded topology. So you can use a single CustomBondForce to implement them with whatever function you need.

1. define regular NonbondedForce with PME and LJ, then set all sigma and epsilon to be 0.

That's actually a special case that it checks for. If it sees that epsilon is 0 for every particle, it #ifdefs out the LJ calculation from the nonbonded kernel so you don't pay any cost for calculating an interaction you aren't using.

Peter

[CHANG YUN SON]

Thank you Peter !

That's actually a special case that it checks for. If it sees that epsilon is 0 for every particle, it #ifdefs out the LJ calculation from the nonbonded kernel so you don't pay any cost for calculating an interaction you aren't using.

That's great !
I actually hoped it might be the case

Thank you for your advice for the 14 interaction.
I was obsessed on this CustomNonbondedForce term so didn't think about defining it in bonded term.
Will see how I do with this implementation and will update if I make it to work or need additional help

Best.
ChangYun

[CHANG YUN SON]

Hi Peter,

I still have problem with this 14 interaction.
If I use the CustomBondForce to calculate this 14 interaction,
I think I need to define those 14 pairs as bond in residue topology.
If that's the case, that will mess the conectivity of the residue,
and so it will increase the risk of adding unintended angle forces or torsional forces along with removing incorrect atom pairs from nonbonded interactions due to the exclusion of nonbonded force.

Do you have any different class to define the 14 type of pairs?
Or am I misunderstanding something here?

I think one solution for this could be adding 14 pairs manually within simulation python script, but not sure how easy it is to do.
Do you have any another idea or suggestions?

ChangYun

[Jason Swails]

You are mistaken here. Every potential energy term must be explicitly added. So adding a CustomBondForce between two particles will not introduce additional Angle or Torsion forces between atoms that will form an Angle when accounting for the bond you just added.

The only thing you may want to be careful about is if you use the option to generate exclusions and exceptions in the nonbonded list from bonds and angles (although that should still work, since you are using the CustomBondForce to implement a nonbonded exception).

[CHANG YUN SON]

Thank you for the reply.

I understand that I can explicitly add every angle and torsional forms based on atomtype rather than class. Though I sill think there might be problem with the exclusions in the nonbonded list.

I want to use CustomNonbondedForce and exclude bonded pairs and 13 pairs, and sclae the nonbonded interaction to be half for 14 pairs. I assumed if I use the exclusion option in CustomNonbondedForce will try to find 12 or 13 pairs based on the bond definition of residue. Though since I need to define the 14 pairs as a bond in residue, topology will treat 15 pairs as 13 pairs and so on. Is this not correct?

If that's the case, the only option I can think of is to explicitly define CustomBondForce for all nonbonded interaction pairs in a residue. Though it seems to hectic if the molecule gets larger.
Is there any other way?

------
I think there could be a less hectic option. To define bond for every 13 pairs as well, but not define any CustomBondForce for those pairs, and use the exclusion=1 option for CustomNonbondedForce. I'll try to use this option and see how it goes.

[Peter Eastman]

Oh, right. If you're doing this from a force field XML file, you need to cheat a little bit. The XML file needs to contain a <script> tag. That contains a Python script that gets executed after the System has been built, and can do any final processing on it to implement things that don't fit in the standard logic.

Take a look at charmm_polar_2013.xml. It uses a <script> tag for exactly this sort of thing. It actually needs a pretty complicated script with a bunch of hardcoded lookup tables, because that's a really complicated force field with lots of unusual features. Yours can be a lot simpler. But like you, it uses NonbondedForce for Coulomb interactions and CustomNonbondedForce for vdW interactions.

Here's how it does this. In the <NonbondedForce> tag, it specifies sigma=1, epsilon=1 for every atom type. That isn't because we actually want those values. One of the things the script does is to set epsilon to 0 for every atom:

Code: Select all

for i in range(nonbondedForce.getNumParticles()):
    (q, sig, eps) = nonbondedForce.getParticleParameters(i)
    nonbondedForce.setParticleParameters(i, q, 1, 0)

The reason for this is so we can distinguish which exceptions are 1-4 interactions (since they'll have epsilon=1), and which are 1-2 or 1-3 interactions (since they'll have epsilon=0). It then loops over all the exceptions, picks out the 1-4 interactions, and modifies their parameters based on values in a set of lookup tables:

Code: Select all

for i in range(nonbondedForce.getNumExceptions()):
    (p1, p2, q, sig, eps) = nonbondedForce.getExceptionParameters(i)
    class1 = typeToClass[data.atomType[data.atoms[p1]]]
    class2 = typeToClass[data.atomType[data.atoms[p2]]]
    index = class1*numAtomClasses+class2
    if eps._value != 0.0:
        eps = epsilon14[index]
    nonbondedForce.setExceptionParameters(i, p1, p2, q, sigma14[index], eps)
    customNonbondedForce.addExclusion(p1, p2)

You can do something similar. Loop over them and set epsilon to 0 for all the exceptions (since you don't want NonbondedForce calculating them). Call addExclusion() on the CustomNonbondedForce so it won't calculate them either. And then, if it's a 1-4 interaction, call addBond() on your CustomBondForce.

Peter

[CHANG YUN SON]

Hi Peter,

I'm still confused by this script.

How does this script preserve the information regarding 1-4 untouched?
After setting all the eps values to be 0 for all atoms, why only 1-4 interactions are having non-zero eps values?

I also posted another question regarding the implementation of my FF,
and I think if there is simple example simulation script of this charmm-polarizable FF, it would be really beneficial for me.

Could you send or show me any examples using this FF?

Thank you for your help !
Best,
ChangYun

[Peter Eastman]

An exception is an interaction that is calculated based on different parameters from the ones associated with the particles. So setting epsilon to 0 for all particles will remove all interactions that are not exceptions, but it has no effect on exceptions. Those are calculated based on the parameter values specified in the exception itself, not based on parameters for the particles.

Exceptions are automatically added for 1-2, 1-3, and 1-4 pairs. The 1-2 and 1-3 pairs have epsilon set to 0, but the 1-4 pairs do not. That's how it tells them apart.

Peter

[Joe Napoli]

I am trying to use CustomNonbondedForce to implement a nonbonded interaction whose parameters depend on the classes of the two atoms in the interaction. My understanding is that I need to define interaction groups. This is a water model, so the groups would be H-O, O-O, and H-H, and I'd like the parameters epsilon, gamma, and sigma below to take on specific values according to which interaction pair is being evaluated. To be clear, I'd like to be able to define something like the following in my XML:

Code: Select all

<CustomNonbondedForce energy="2*epsilon/(1-3/(gamma+3))*(sigma^6/(sigma^6+r^6))*(3/(gamma+3)*exp(gamma*(1-r/sigma))-1)" bondCutoff="0">
  <InteractionGroup class1="O" class2="O" sigma="0.36174" epsilon="1.0334" gamma="13.256"/>
</CustomNonbondedForce>

What is the most straightforward way to achieve this? From the API docs I can see how one does it in a script, but in the XML case it is less clear to me whether a <Script> tag is necessary.

Thank you very much!

-Joe