Simple Binary Material

[Gideon Simpson]

I am interested in using OpenMM to simulate a simple binary material where I have three types of pair potentials for A-A, B-B, and A-B bonds. I've been looking at some of the documentation for this. Is the right approach here that I build an XML force file for this and a separate PDB file for the atomic positions? Or is there a way to do all/part of this "online" within the python script?

[Peter Eastman]

It's possible to do everything within the Python script, if that's what you prefer. If you describe a bit more about your system, I can suggest ways of doing it. What's the form of the interaction? Is it a bonded (only specific pairs of particles interact) or nonbonded (everything interacts with everything) interaction? How many particles do you have?

[Gideon Simpson]

In principle, this is a non-bonded system, with different pair potentials depending on specifies. The simplest example i would like to do is a dimer in an LJ fluid. Here, Atoms 1 and 2 interact with each other via a bistable potential, and they interact with atoms 2,3,...,N via an LJ pair potential. Atoms 2,3,...,N then interact with one another via an LJ pair potential (with possibly different parameters).

The more complicated problem would be to have N atoms, with approximately half species A and half species B. A atoms interact with one another through one type of pair potential, say LJ, B atoms interact with one another through a different pair potential, say Morse, and A-B interactions are through some third, similar, pair potential.

My main interest in being able to set this up in python is that I would like to be able to rapidly change system size (number of atoms) and their initial configurations.

[Peter Eastman]

What are the actual functional forms you want to use for the interactions? I don't just mean "for example", I mean what are the specific ones you want to use? Likewise, what is the specific range of system sizes (numbers of particles) you want to simulate? These details matter! I'm going to suggest different approaches depending on the answers.

[Gideon Simpson]

For a first problem, I want to study the solvated dimer problem from Example 1.3.2.4 of Lelievre, Rousset, and Stolz (if you have that book handy). In it, particles 1 and 2 interact through the pair potential

V(r) = h * (1 - (r-r0-w)^2/w^2)^2, r = distance between particles 1 and 2, and h, r0, and w are parameters

Particles 3-N interact with one another via LJ pair potentials with finite range. Particles 3-N also interact with particles 1 and 2 with the same LJ pair potentials. In the text I am referring to, LJ potential is only C0, but not C1, but this is not essential to me. To replicate the text examples, they use a total of N = 100 particles in 2D with periodic BCs. I would likely extend this, eventually, to 3D with N = 10^3 or 10^4 atoms.

For a second problem, I would look at a binary system of particles in a 2D/3D periodic domain with at least N = 10^2 and likely upwards of 10^4 particles. In this setting, the three types of pair interactions would be governed by standard LJ pair potentials with finite range interactions, but each with its own set of parameters (i.e., sigma and epsilon values for A-A, B-B, and A-B pairs). Roughly half would be A particles and half would be B particles.

I think maybe the point is that the sorts of potentials I am currently imagining are all posed in terms of operations found on a scientific calculator (polynomials, exp, log, trig functions), with a modest number of parameters (no more than 4).

[Peter Eastman]

In the first case, nearly all atoms interact with the same potential and same parameters, right? The only exception is that particles 1 and 2 use a different form? For that case, you can just go ahead and use a standard NonbondedForce. Then call addException() on it to tell it to skip that one particular pair. To implement the interaction between those two particles, use a CustomBondForce and give it the particular functional form you want.

For the second case, everything is LJ and the only difference is in the parameters, right? Are the A-B parameters derived from the A-A and B-B ones by a standard combining rule? If so, you can use NonbondedForce for that too. If not, use a CustomNonbondedForce with a lookup table for the parameters. I can point you to example code when you get to that point.

[Gideon Simpson]

I'm having some difficulty with the dimer model problem. I attempted to use the commands:

Code: Select all

# remove NonbondedForce interaction between 0 and 1
force.addException(0, 1, 0.0*charge*charge, sigma, 0.0*epsilon)
system.addForce(force)
# add custom force between 0 and 1 (simple version, for now)
dimer_force = openmm.CustomBondForce("0.5 * k *  (r - r0)**2")
dimer_force.addPerBondParameter("k")
dimer_force.addPerBondParameter("r0")
r0 = 1.5 * unit.nanometers
k = 1000 * unit.kilojoules_per_mole / unit.nanometers ** 2

dimer_force.addBond(0,1, [k, r0])
system.addForce(dimer_force)

This executes without issue, except when I get to:

Code: Select all

simulation = app.Simulation(topology, system, integrator)

I get the error:

Code: Select all

libc++abi.dylib: terminating with uncaught exception of type Lepton::Exception: Parse error in expression "0.5 * k *  (r - r0)**2": unexpected token: *
Abort trap: 6

[Peter Eastman]

Exponentiation is indicated with ^ not **. See http://docs.openmm.org/latest/userguide ... xpressions.

[Gideon Simpson]

Yup, that did it.

[Gideon Simpson]

I'm seeing some slightly unexpected behavior when using a CustomBondForce to represent the interaction of the two atoms making up the dimer. I'm trying to run on a periodic domain, but I see the two dimer atoms drifting out of the periodic box. I've tried using setUsesPeriodicBoundaryConditions(True)