Python API to evaluate ewald Coulomb force

[JZ Most]

Thanks for replying. The values of positions are from 0 to 1e-8m, I transfered them to the unit of angstrom. Thats why I multiple them with 1e+10

[Peter Eastman]

Can you provide a complete, working example that demonstrates the problem? The code you posted before doesn't correspond to what you described since it doesn't use a real integrator, doesn't run any dynamics, and never retrieves the positions.

[JZ Most]

Thanks for replying. The physical system I tried to simulate is 'a cubic system (L = 1e-8m), with 2 layer of fixed charges at x = 0 and x = L/2, and certain number of free ions '. There are two approaches I implemented, one is to use openmm to calculate the PME coulombic and LJ forces and use python random number to simulate the random diffusion, so I can update the positions of free ions every time step. Another is to just use the Langevin integrator to run the simulation to a certain timepoint. The reason for why I did not use a real integrator is that I only want to know the force value. However, no matter which integrator I used (include real integrator), I still have my position array became 0 after 1 time step. So I am a bit confused.

[Peter Eastman]

Please provide a complete example I can use to reproduce the problem. Otherwise, I don't really know what you're doing or why it isn't working!

[JZ Most]

Sorry about being confusing, please see below,
Initialize the positions

Code: Select all

import numpy as np
from simtk import openmm as mm
from simtk.openmm import app
from simtk.unit import nanometer, picosecond, angstrom, kilojoule_per_mole, kelvin, kilocalorie_per_mole
import openmm.unit as unit
import time
from matplotlib import pyplot as plt

length = 1e-8
num_particles = 200
rho_perm_pos = np.random.rand(36, 3) * length ##fixed charge at x = 0
rho_perm_pos[:,0] = 0

rho_perm_neg = np.random.rand(36, 3) * length ##fixed charge at x = L/2
rho_perm_neg[:,0] = 0.5*length
np.random.seed(42)
positions_pos = np.random.rand(num_particles, 3) * length
positions_neg = np.random.rand(num_particles, 3) * length
positions = np.random.rand(num_particles*2, 3) * length

positions = np.concatenate((positions_neg, rho_perm_neg, positions_pos, rho_perm_pos), axis=0)

PME coulombic forces and LJ forces using Openmm,

Code: Select all

box_size = 100*angstrom

# define a new system
system = mm.System()
system.setDefaultPeriodicBoxVectors((box_size, 0, 0), (0, box_size, 0), (0, 0, box_size))

# define the ion number and charge
num_Na = num_particles+36 ## free ions and 36 fixed
num_Cl = num_particles+36
charge_Na = +1.0 # +
charge_Cl = -1.0 # -
mass_Na = 22.98976
mass_Cl = 35.453

# charge array
charges = np.concatenate([np.full(num_Cl, charge_Cl), np.full(num_Na, charge_Na)])

# add ions to the system
for _ in range(num_Cl):
    system.addParticle(mass_Cl) 
for _ in range(num_Na):
    system.addParticle(mass_Na)
    
# NonbondedForce, Particle Mesh Ewald
nonbonded = mm.NonbondedForce()
nonbonded.setNonbondedMethod(mm.NonbondedForce.PME)
nonbonded.setEwaldErrorTolerance(1.0e-4)
nonbonded.setCutoffDistance(50*angstrom)

alpha = 0.052/nanometer
# Set PME parameters (alpha, reciprocal space cutoff)

grid_size_x = 17  # Example FFT grid size for x dimension
grid_size_y = 17  # Example FFT grid size for y dimension
grid_size_z = 17  # Example FFT grid size for z dimension

# Set PME parameters (alpha, reciprocal space cutoff, FFT grid sizes)
nonbonded.setPMEParameters(alpha, grid_size_x, grid_size_y, grid_size_z)

# add the charge to the system
for _ in range(num_Cl):
    nonbonded.addParticle(-1, 0.227 * nanometer, 0.15*4.184 * kilojoule_per_mole)
for _ in range(num_Na):
    nonbonded.addParticle(+1, 0.141 * nanometer, 0.0469*4.184 * kilojoule_per_mole)

# add the force to the system
system.addForce(nonbonded)

# create simulation
dummy_integrator = mm.CustomIntegrator(0)
simulation = app.Simulation(topology=None, system=system, integrator=dummy_integrator, platform=platform)
simulation.context.setPositions(positions*1e+10*angstrom)
state = simulation.context.getState(getForces=True)
forces = state.getForces(asNumpy=True)

From this code, I can get a forces array seem meaningful (my main purpose in this setting is to get force value)

Code: Select all

Quantity(value=array([[ -619.46679688,   468.09777832,   -35.26379395],
       [ -345.91641235,   -78.44809723,    90.83336639],
       [-1286.58483887, -1409.69042969,  3650.68798828],
       ...,
       [ -240.84684753,    78.62600708,   530.60791016],
       [  148.15328979,   503.44345093,   212.35394287],
       [ -301.76553345,  3715.84716797,  -208.17672729]]), unit=kilojoule/(nanometer*mole))

However, when I changed the integrator, (I want to validate with the results by langevin)

Code: Select all

dummy_integrator = mm.LangevinIntegrator(300 * unit.kelvin, 0.1 / unit.picosecond, 0.01 * unit.picoseconds)
simulation = app.Simulation(topology=None, system=system, integrator=dummy_integrator, platform=platform)
simulation.context.setPositions(positions * 1e+10 * angstrom)
num_steps = 1000  # Number of steps to run
simulation.step(num_steps)
state.getPositions(asNumpy=True)

The positions will not be change (but the forces value is the same to Customintegrator(0)) until the 'Particle coordinate is NaN.' Does that mean I did not get the correct value?

[Peter Eastman]

Code: Select all

simulation.step(num_steps)
state.getPositions(asNumpy=True)

You don't retrieve a new State object. You call getPositions() on some state you created earlier, though I don't know when or how you created it. That isn't shown in your code. If you want to get the current positions, you need to call getState(getPositions=True).

[JZ Most]

Hi Peter,

Thanks dor reply, I am sure I created a state before as

Code: Select all

state = simulation.context.getState(getForces=True)

Now I change the sample code according to your mention.

Code: Select all

length = 1e-8
num_particles = 50
rho_perm_pos = np.random.rand(36, 3) * length # 36 fixed ions at x = 0
rho_perm_pos[:,0] = 0

rho_perm_neg = np.random.rand(36, 3) * length # 36 fixed ions at x = L/2
rho_perm_neg[:,0] = 0.5*length
np.random.seed(42)
positions_pos = np.random.rand(num_particles, 3) * length # free +ve ions
positions_neg = np.random.rand(num_particles, 3) * length # free -ve ions

positions = np.concatenate((positions_neg, rho_perm_neg, positions_pos, rho_perm_pos), axis=0)

box_size = 100*angstrom

# define a new system
system = mm.System()

system.setDefaultPeriodicBoxVectors((box_size, 0, 0), (0, box_size, 0), (0, 0, box_size))

# define the ion number and charge
num_Na = num_particles+36
num_Cl = num_particles+36
charge_Na = +1.0 # +
charge_Cl = -1.0 # -
mass_Na = 22.98976
mass_Cl = 35.453

# charge array
charges = np.concatenate([np.full(num_Cl, charge_Cl), np.full(num_Na, charge_Na)])

# add ions to the system
for _ in range(num_particles):
    system.addParticle(mass_Cl) # assume the mass is 1 for simplicity
for _ in range(36):
    system.addParticle(0) # assume the mass is 1 for simplicity
for _ in range(num_particles):
    system.addParticle(mass_Cl) # assume the mass is 1 for simplicity
for _ in range(36):
    system.addParticle(0) # assume the mass is 1 for simplicity
    


# NonbondedForce, Particle Mesh Ewald
nonbonded = mm.NonbondedForce()
nonbonded.setNonbondedMethod(mm.NonbondedForce.PME)
nonbonded.setEwaldErrorTolerance(1.0e-4)
nonbonded.setCutoffDistance(50*angstrom)

epsilon = 1.0e-4  # Ewald error tolerance

alpha = 0.052/nanometer
# Set PME parameters (alpha, reciprocal space cutoff)

grid_size_x = 17  # Example FFT grid size for x dimension
grid_size_y = 17  # Example FFT grid size for y dimension
grid_size_z = 17  # Example FFT grid size for z dimension

# Set PME parameters (alpha, reciprocal space cutoff, FFT grid sizes)
nonbonded.setPMEParameters(alpha, grid_size_x, grid_size_y, grid_size_z)

# add the charge to the system and ignore VDW
for _ in range(num_Cl):
    nonbonded.addParticle(-1, 0.440 * nanometer, 0.4928 * kilojoule_per_mole)
for _ in range(num_Na):
    nonbonded.addParticle(+1, 0.325 * nanometer, 0.544 * kilojoule_per_mole)

# add the force to the system
system.addForce(nonbonded)

dummy_integrator = mm.LangevinIntegrator(300 * unit.kelvin, 1 / unit.picosecond, 0.01 * unit.picoseconds)
simulation = app.Simulation(topology=None, system=system, integrator=dummy_integrator, platform=platform)
simulation.context.setPositions(positions * 1e+10 * angstrom)
num_steps = 10000 # Number of steps to run
simulation.step(num_steps)

state = simulation.context.getState(getForces=True, getPositions=True)
ps = state.getPositions(asNumpy=True)

I set the timestep as 1e-2 picoseconds, but the simulation will crush due to some of the position become Nah. If I further decrease the time step to 1e-4 picoseconds, the simulation will not return a error, but mosst of the positions of free ions eventually accumulate at the fixed charge. I think I define something wrong and seems like the LJ forces did not apply.

[Peter Eastman]

Some of your settings are very usual. For example, you set the cutoff distance to 50 A, which is huge. Values around 10 A are more typical. A timestep of 0.01 ps is larger than people typically use, but 1e-4 ps is much smaller than normally used. Setting all the particles to random initial positions is also dangerous, since there could easily be two particles right on top of each other. That's probably the reason you need such a small step size to make it stable. It's better to spread them out more evenly. At the very least, call simulation.minimizeEnergy() to resolve the largest forces before you try to simulate.

There are two methods of setting Ewald parameters. One is to let it automatically pick them based on the cutoff distance and error tolerance. The other is to set them manually. You're trying to do both at once. You first call setEwaldErrorTolerance(), but then setPMEParameters() which overrides it and means the error tolerance you set is ignored.

How close to each other do the charges end up?

[JZ Most]

Thanks for replying. I just ended my travel. Your suggestions are very helpful, I changed the cutoff distance, and used latin hypercube for sampling the initial positions, before simulation I also minimize the energy. These changes allow the system to be simulated in 0.01 ps timestep. However, the results did not change, the +ve and -ve ions all accumulated in the fixed charged layer, the space between fixed charges are empty, seems like LJ forces did not play its role. The Charges ended up with a minimum Euclidean distance of ~0.2-0.3 nm, do you think this is too closed?

[Peter Eastman]

If the LJ force weren't being applied, the distance between opposite charges would shrink without limit until the massive forces caused the simulation to explode. You set sigma to 0.227 nm for Cl and 0.141 nm for Na, so distances in the range of 0.2 to 0.3 nm are about what you expect.

What were you expecting to happen? It doesn't seem surprising that the mobile charges are attracted to the fixed ones and end up close to them.