CustomNonbondedForce with pair dependent constants

[Dennis Della Corte]

Hello,

I am looking for an elegant way to encode non-bonded interactions that are dependent on the atom-pair for an ions only simulation.

Example:
Atomtype A
Atomtype B
Atomtype C

Functional form:
V(r) = interaction_AB / r_AB

I would need to parametrize a forcefield to contain all the pairwise terms:

interaction_AB = 1
interaction_AC = 2
interaction_BC = 3

My application will have many atom types, and a more complicated functional form with ~8 pairwise parameters.
There is no simple averaging of per-atom parameters possible, as implemented for LJ interactions.

Any advice is greatly appreciated.

Best wishes,
Dennis

[Peter Eastman]

The best way of doing this is to assign a per-particle parameter called "type" which is just an integer 0, 1, 2, ... with the index of the atom type for each atom. You can then use those to look up the value in a Discrete2DFunction. This is how we handle Lennard-Jones interactions for CHARMM force fields that involve NBFIX terms. You can see the code at https://github.com/openmm/openmm/blob/d ... 2524-L2527:

Code: Select all

self.force = mm.CustomNonbondedForce('acoef(type1, type2)/r^12 - bcoef(type1, type2)/r^6')
self.force.addTabulatedFunction('acoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, acoef))
self.force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, bcoef))
self.force.addPerParticleParameter('type')

[Dennis Della Corte]

Thank you, Peter. I am struggling with addTabulatedFunction inside of the forcefield .html:

<CustomNonbondedForce energy="(Bij(lookup1,lookup2)/r)^2" bondCutoff="3">

<PerParticleParameter name="lookup"/>
<Atom type="0" lookup="0"/>
<Atom type="1" lookup="1"/>

<Function name="Bij" type="Discrete2DFunction" xsize="1" ysize="1">

</Function>

</CustomNonbondedForce>

Results in the error:
Parse error in expression "(Bij(lookup1,lookup2)/r)^12": unknown function: Bij

Also: system.getForces()[0].getNumTabulatedFunctions()
results in [0]

How do I correctly add the tabulated function inside of html definition?

[Peter Eastman]

See http://docs.openmm.org/latest/userguide ... -functions. It describes how to add a tabulated function in a force field.

[Dennis Della Corte]

I have been digging through the documentation, but maybe html syntax is the issue here. When it says "To define a function, include a <Function> tag inside the <CustomNonbondedForce> or <CustomGBForce> tag", I thought it means:

<CustomNonbondedForce energy="tabulatedfct(lookup1,lookup2)" bondCutoff="3">

...

<Function name="tabulatedfct" type="Discrete2DFunction" xsize="1" ysize="1">
10
</Function>

</CustomNonbondedForce>

Sorry for the noob question, but where should I place the <Function> tag, so that the definition is available in the energy namespace?

[Peter Eastman]

Yes, that looks correct. By the way, it should be an XML file, not an HTML one.

[Dennis Della Corte]

Ok, in that case it might be a bug? I tried it out on two different machines, both times resulting in the same error msg:

Exception: Parse error in expression "(tabulatedFct(lookup1,lookup2)/r)^2": unknown function: tabulatedFct

here the files to reproduce the error: https://byu.box.com/s/7vnrlpd86dww79my5r5ezrve195azpoz

[Peter Eastman]

You have the wrong value for the "type" attribute. It should be

Code: Select all

type="Discrete2D"

instead of

Code: Select all

type="Discrete2DFunction"

[Dennis Della Corte]

issue solved... thank you.

[Wei Lu]

Hello Peter,
I have a similar problem. but a little more complicate. The potential I want to add is not a LJ potential.(it's defined by tabulated continuous function) so I would hope to have a mixed Discrete-Continuous tabulated function. To be more specific, the input array will be 3D. The first 2 dimension is treated as discrete(indicate the type of atom), the last dimension should be treated as data for spline fit.
So I can write something like:
mm.CustomNonbondedForce("mixedPotential(type1, type2, r)")
self.force.addTabulatedFunction('mixedPotential', mm.DiscreteMixedContinuous3DFunction(numLjTypes, numLjTypes, n_points, mixedPotential))
self.force.addPerParticleParameter('type')

This is possible? or you have a better way to achieve this? Thanks.