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?