Re: Python API to evaluate ewald Coulomb force
Posted: Fri Jun 28, 2024 3:44 am
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
A short text to describe your forum
https://simtk.org/plugins/phpBB/
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)