CustomNonbondedForce and exceptions/exclusions

[Michal Krompiec]

Hi,
I'm developing a force field which has "standard" Coulomb forces (like in OPLS) and a custom dispersion/repulsion term. I'd like the Coulombic forces and CustomNonbondedForce to be scaled by 0.5 for 1-4 interactions.
Initially, I tried to use <NonbondedForce> with epsilon=sigma=0 and the charges from my model, <CustomNonbondedForce bondCutoff=4> and <CustomBondForce> for all 1-3 interactions, but I got this error Exception: All Forces must have identical exceptions.
To get rid of this exception, I set bondCutoff for the custom NB to 3, and subtracted half of the non-bonded force using CustomBondForce, but this is a dirty workaround which I really don't like.
I would be very grateful if someone could suggest a better solution. For example, is it possible to change the bondCutoff for Coulomb to 4, and then add the missing 1-3 Coulomb via CustomBondForce?
Thanks,
Michal

[Peter Eastman]

Code: Select all

<CustomNonbondedForce bondCutoff=4>

Are you sure that's really what you wanted? bondCutoff=4 means to exclude all pairs separated by up to 4 bonds, which is to say, 1-5 interactions. If you only want to exclude 1-4 interactions, you specify bondCutoff=3.

[Michal Krompiec]

Thank you very much for your reply. Oh yes, you are right. Indeed, bondCutoff should be set to 3, and I will add scaled 1-4s via CustomBondForce. Do I understand correctly that NonbondedForce with zeroed epsilons includes scaled 1-4 Coulomb? (assuming the scaling parameter is set correctly)
Thanks,
Michal

[Peter Eastman]

Correct. Also, if you set all the epsilons to 0 then it ifdefs out the parts of the kernels that compute LJ interactions, so there's no performance cost from them.

[Michal Krompiec]

I still can't get this right. At the moment, I have a CustomNonbondedForce with bondCutoff=3 and a NonbondedForce which handles Coulomb interactions (1-4 interactions are scaled by 0.5). Now, I would like to add 1-4 interactions to the CustomNonbondedForce (scaled by 0.5, or unscaled). Is it possible to do it in the xml (without modifying anything in python)? Such a solution would be preferred because I'd like to keep my FF compatible with ForceBalance.
If I do it by adding 1-4 interactions as bonds (in <Residue>) with a CustomBondForce, these interactions will become "bonded" in the topology and will disappear from NonbondedForce, and any neighboring atoms will become 1-3 interactions (instead of 1-5). So using this https://github.com/jmcdaniel43/Thermody ... ionpol.xml as an example won't help.

[Peter Eastman]

You actually can put Python code right in the XML file. It goes in a <Script> tag. That gives you a chance to post-process the System in whatever way you need, such as by adding a CustomBondForce. For an example of a force field that uses this mechanism, take a look at https://raw.githubusercontent.com/openm ... r_2013.xml. Within the script, the variable "sys" refers to the System that was just created and "data" is a ForceField._SystemData object containing lots of useful information about bonds, angles, dihedrals, atom types, etc.

[Michal Krompiec]

Thank you, this is great! When exactly is this code executed?

[Peter Eastman]

It's the very last step, after all standard processing has been done and just before the System gets returned:

https://github.com/openmm/openmm/blob/m ... 1331-L1334