Core-Shell potential...

[Kasra Momeni]

Hello,

I wondered to know if I can implement core-shell inter-atomic potential using OpenMM?

regards,
Kasra.

[Peter Eastman]

Can you provide a link to the specific potential you're interested in using?

Peter

[Kasra Momeni]

Here is the title for the paper:

Piezoelectric constants for ZnO calculated using classical polarizable core–shell potentials

doi: 10.1088/0957-4484/21/44/445707

regards,
Kasra.

[John Chodera]

Hi Kasra,

I've just skimmed the reference, but it states:

"The core–shell potentials for ZnO that we consider in the present work take a Buckingham-type form"

and lists a Coulomb plus exponential-six Buckingham potential as Eq. 1. These terms should be easily implementable in OpenMM using the CustomNonbondedForce class.

Coulomb potentials (q_i q_j / r_{ij}) can already be handled by the NonbondedForce class, but you can implement just the Buckingham exponential-six part, or even the whole of Eq. 1, with a CustomNonbondedForce using code like the following:

// Implement just the Buckingham exponential-six part.
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("A*exp(-r/rho) - C/r^6;");
// more here...
system.addForce(customNonbonded)

Now, the difficulty is that we need to also specify the parameters A, rho, and C. If there were simple mixing rules, like A = 0.5*(A1+A2), where interacting atoms 1 and 2 had A parameter of A1 and A2, you could simply include those in the energy expression fed to CustomNonbondedForce. However, it looks like Table 1 lists the O-O, O-Zn, and Zn-Zn sets of parameters as all being different---is that right?

I don't think OpenMM's CustomNonbondedForce supports a lookup table for these parameters. In the past, I've had to get around this by defining per-particle parameters like "isO" which will be 1 if the particle is an oxygen and 0 if it is zinc; and "isZn" which will be 1 if the particle is a zinc and 0 if it is an oxygen. Then, the energy expression becomes much more complex:


// Implement just the Buckingham exponential-six part for Binks parameters.
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("A*exp(-r/rho) - C/r^6; A = 9547.96*isO1*isO2 + 529.70*(isO1*isZn2+isO2*isZn1); rho = 0.21916*isO1*isO2 + 0.3581*(isO1*isZn2+isO2*isZn1); C = 32.0*isO1*isO2;");

All of that new junk just selects the appropriate parameters from Table 1 depending on the identity of the atoms in the interacting pair. It's a bit unpleasant, but it at least works. You'll also have to change the numbers that appear from eV and A to the natural OpenMM units of kJ/mol and nm.

Peter: Any further suggestions?

- John

[Kasra Momeni]

Thanks John. I have already implemented that Buckingham code using custom force and coulomb force.
My main question is on implementation of Core-Shell model. In this case I have to model the nucleus as a massive particle and electron electron as a massless shell which interacts with the core. This interaction can be modeled with a harmonic potential.
I do not know how I should model this massive core- massless shell interactions using OpenMM.

I am still looking forward for your help.

With best regards,
Kasra.

[John Chodera]

Hi Kasra,

Ah! Apologies for not understanding your question earlier.

It seems there might be (at least) two ways to implement the mobile shell within OpenMM:

#1. If you want to add the shell center of mass as another pseudoparticle, you can add a harmonic potential that restrains the shell particle to the core particle with a spring constant K (the one appearing in Eq. 2). You'd have to then create a CustomNonbondedForce for describing the interactions of the shell with other shells and with cores, where you would have to work out the shell-shell and shell-point charge Coulomb energies. There would also have to be some cleverness in exclusions such that the shell does not interact with its own core electrostatically. Additionally, you couldn't make the shell massless in OpenMM---you'd have to use a small mass that would have to be tuned to the core mass and timestep, in a manner similar to Carr-Parinello dynamics. There might be all sorts of problems to watch out for here.

#2: The much simpler way is to exploit the isomorphism with having a single polarizable dipole of polarizability \alpha discussed in Section 2. OpenMM is now supporting the AMOEBA polarizable forcefield, which introduces polarizable point dipoles (as well as fixed multipoles). You could certainly make use of these polarizable point dipoles to implement the polarizable shells. These point dipoles *are* solved self-consistently each timestep, so will behave as if you had your massless shells.

Peter Eastman can advise about the status of AMOEBA. I think the capability you need is probably already complete in the development version, and will almost certainly be in the next release.

[Peter Eastman]

AMOEBA will be in the next release, which hopefully will be available quite soon (within a few weeks).

Peter

[Spela Ivekovic]

Hi,

I have a follow-up question on this.

I'd like to use this approach for a Lennard-Jones potential with custom sigma and epsilon values for each pair of species but I'll ask the question in terms of the potential discussed in this thread.

In order to identify the type of the particle at the time of initialisation, we do something like this:

force->addPerParticleParameter("isO");
force->addPerParticleParameter("isZn");
vector<int> params(2);
// Oxygen
params[0] = 1;
params[1] = 0;
force->addParticle(params);
// Zinc
params[0] = 0;
params[1] = 1;
force->addParticle(params);
...

However, at the time of force calculation, the following equations are used:

A*exp(-r/rho) - C/r^6;
A = 9547.96*isO1*isO2 + 529.70*(isO1*isZn2+isO2*isZn1);
rho = 0.21916*isO1*isO2 + 0.3581*(isO1*isZn2+isO2*isZn1);
C = 32.0*isO1*isO2;

Does the parser know that any parameter name with a number after it refers to one of the particles in the pair at the time of interaction or is this specified explicitly by the user somewhere? I.e., if both particles interacting are Oxygens, their parameters are isO=1 and isO=1 - do they become isO1 and isO2 when involved in a pairwise force calculation?

Thank you!

Spela

[Peter Eastman]

if both particles interacting are Oxygens, their parameters are isO=1 and isO=1 - do they become isO1 and isO2 when involved in a pairwise force calculation?

Exactly. For every per-particle parameter you define, CustomNonbondedForce creates two variables whose names are the parameter name followed by "1" and "2", containing the values of the parameter for the two interacting particles.

Peter

[Spela Ivekovic]

Hi Peter,

we've implemented a custom LJ potential through the procedure described above and we have found that a simple waterbox test case is around 13-times slower than when the default nonbonded force calculation is used. Is that to be expected?

Thanks,
Spela