I am simulating a system of electrolytes where charge regulation occurs so the total particle number is changing during the simulation. However, in the production run process of my script, an error of "Particle node not found" is observed. I confirm if I remove the lines below, the script works well:
So I am confused about why these lines lead to an error. What is the correct way to output the configuration so that one can analyze the configuration? My script is attached below, thank you!
from __future__ import print_function
from __future__ import division
import espressomd
from espressomd import code_info
from espressomd import analyze
from espressomd import integrate
from espressomd.interactions import *
from espressomd import reaction_ensemble
from espressomd import interactions
from espressomd import electrostatics
from espressomd import electrostatic_extensions
from espressomd.shapes import Wall
import sys
import numpy as np
import time
box_lxy = 20.0
box_lz= 10.0 #distance between planar interfaces
fraction_empty = 1.5
system = espressomd.System(box_l=[box_lxy, box_lxy, (1+fraction_empty)*box_lz],periodicity = [1,1,1],time_step=0.01)
l_bjerrum = 3.0
temperature = 1.2
system.thermostat.set_langevin(kT=temperature, gamma=100.0)
system.set_random_state_PRNG()
np.random.seed(seed=system.seed)
system.cell_system.skin = 2
system.cell_system.max_num_cells = 2744
# Particle setup
N0 = 10 # total number of monomers=N0*N0*2
nNaOH = 0 # number of initial Na+OH- (additional OH-'s)
nHCl = 0 # number of initial H+Cl- (additional H+'s)
type_HA = 0 # type 0 = HA
type_A = 1 # type 1 = A-
type_H = 2 # type 2 = H+
type_OH = 3 # type 3 = OH-
type_Na = 4 # type 4 = Na+
type_Cl = 5 # type 5 = Cl-
type_sub = 6
charges = {}
charges[type_HA] = 0
charges[type_A] = -1
charges[type_H] = 1
charges[type_OH] = -1
charges[type_Na] = 1
charges[type_Cl] = -1
charges[type_sub] = 0
# Walls
system.constraints.add(shape=Wall(dist=0, normal=[0, 0, 1]), penetrable=False, particle_type=type_sub)
system.constraints.add(shape=Wall(dist=-box_lz, normal=[0, 0, -1]), penetrable=False, particle_type=type_sub)
c = 0
for i in range(N0):
for j in range(N0):
system.part.add(pos=[i*box_lxy/N0, j*box_lxy/N0, 0.5], id=c,
type=type_A, q=charges[type_A], mass=1, fix=[1,1,1])
c = c + 1
for i in range(N0):
for j in range(N0):
system.part.add(pos=[i*box_lxy/N0, j*box_lxy/N0, box_lz-0.5], id=c,
type=type_A, q=charges[type_A], mass=1, fix=[1,1,1])
c = c + 1
for i in range(N0*N0*2):
system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy, 1+np.random.random()*(box_lz-2)],
id=c, type=type_H, q=charges[type_H], mass=1)
c = c + 1
# setting up other ions
for i in range(nNaOH):
system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy,
np.random.random()*box_lz],
id=c, type=type_OH, q=charges[type_OH], mass=1)
c = c + 1
for i in range(nNaOH):
system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy,
np.random.random()*box_lz],
id=c, type=type_Na, q=charges[type_Na], mass=1)
c = c + 1
for i in range(nHCl):
system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy,
np.random.random()*box_lz],
id=c, type=type_H, q=charges[type_H], mass=1)
c = c + 1
for i in range(nHCl):
system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy,
np.random.random()*box_lz],
id=c, type=type_Cl, q=charges[type_Cl], mass=1)
c = c + 1
#h5 = h5md.H5md(filename="trajectory.h5", write_pos=True, write_vel=True)
# setting up LJ-interactions
lj_eps = 1.0
lj_sig = 1.0
lj_cut = 1.122462048
lj_shift = 0
for i in range(6):
for j in range(6):
system.non_bonded_inter[i, j].lennard_jones.set_params(epsilon=lj_eps,
sigma=lj_sig, cutoff=lj_cut, shift=lj_shift)
for i in range(6):
system.non_bonded_inter[i, type_sub].lennard_jones.set_params(epsilon=lj_eps,
sigma=lj_sig*0.5, cutoff=lj_cut*0.5, shift=lj_shift)
# Setting up reaction
mode = "constant_pH_ensemble" #"reaction_ensemble"
K_diss = 10.0**(-7)
RE = None
if(mode == "reaction_ensemble"):
RE = reaction_ensemble.ReactionEnsemble(temperature=temperature, exclusion_radius=1)
elif(mode == "constant_pH_ensemble"):
RE = reaction_ensemble.ConstantpHEnsemble(temperature=temperature, exclusion_radius=1)
RE.constant_pH = 7
# HA <--> A- + H+
RE.add_reaction(gamma=K_diss, reactant_types=[type_HA], reactant_coefficients=[1], product_types=[
type_A, type_H], product_coefficients=[1, 1],
default_charges={type_HA: charges[type_HA], type_A: charges[type_A], type_H: charges[type_H]})
print(RE.get_status())
RE.set_volume (box_lxy*box_lxy*(box_lz-1))
RE.set_wall_constraints_in_z_direction (0.5,box_lz-0.5)
system.setup_type_map([type_HA, type_A, type_H, type_OH, type_Na, type_Cl])
system.force_cap = 500
for i in range(50):
system.integrator.run(200)
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0)*(len(system.part.select(type=2))))
print(i, system.number_of_particles(type=type_HA),
system.number_of_particles(type=type_A),
system.number_of_particles(type=type_H),
system.number_of_particles(type=type_OH),
system.number_of_particles(type=type_Cl),
system.number_of_particles(type=type_Na),
energy['total'], energy['coulomb'], energy["kinetic"], temp_measured,system.analysis.min_dist())
system.force_cap = 0
p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=1e-5)
system.actors.add(p3m)
elc = electrostatic_extensions.ELC(gap_size=box_lz*fraction_empty, maxPWerror=1e-5, delta_mid_top=0.9, delta_mid_bot=0.9)
system.actors.add(elc)
for i in range(5):
RE.reaction()
system.integrator.run(200)
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0)*(len(system.part.select(type=2))))
print(i, system.number_of_particles(type=type_HA),
system.number_of_particles(type=type_A),
system.number_of_particles(type=type_H),
system.number_of_particles(type=type_OH),
system.number_of_particles(type=type_Cl),
system.number_of_particles(type=type_Na),
energy['total'], energy['coulomb'], energy["kinetic"], temp_measured)
print("\n<production run starts>\n")
f_config = open("config.txt", 'w+')
for i in range(10):
RE.reaction()
system.integrator.run(200)
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0)*(len(system.part.select(type=2))))
number_of_particles = system.number_of_particles(type=type_HA) + system.number_of_particles(type=type_A) + system.number_of_particles(type=type_H) + system.number_of_particles(type=type_OH) + system.number_of_particles(type=type_Cl) + system.number_of_particles(type=type_Na)
print(i, number_of_particles, system.number_of_particles(type=type_HA),
system.number_of_particles(type=type_A),
system.number_of_particles(type=type_H),
system.number_of_particles(type=type_OH),
system.number_of_particles(type=type_Cl),
system.number_of_particles(type=type_Na),
energy['total'], energy['coulomb'], energy["kinetic"], temp_measured)
print("frame ID", file=f_config)
print(i, file=f_config)
print("number of particles", file=f_config)
print(number_of_particles, file=f_config)
print("type", file=f_config)
for k in range(number_of_particles):
print(system.part[k].type, file=f_config)
print("position", file=f_config)
for k in range(number_of_particles):
print(system.part[k].pos, file=f_config)