Python API to evaluate ewald Coulomb force
Re: Python API to evaluate ewald Coulomb force
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
- Posts: 2551
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Python API to evaluate ewald Coulomb force
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.
Re: Python API to evaluate ewald Coulomb force
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
- Posts: 2551
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Python API to evaluate ewald Coulomb force
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!
Re: Python API to evaluate ewald Coulomb force
Sorry about being confusing, please see below,
Initialize the positions
PME coulombic forces and LJ forces using Openmm,
From this code, I can get a forces array seem meaningful (my main purpose in this setting is to get force value)
However, when I changed the integrator, (I want to validate with the results by langevin)
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?
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)
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)
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))
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)
- Peter Eastman
- Posts: 2551
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Python API to evaluate ewald Coulomb force
Code: Select all
simulation.step(num_steps)
state.getPositions(asNumpy=True)