[ESPResSo-users] use espressomd to model chemical reaction

[Jiaxing Yuan]

Dear all,

I am simulating a confined electrolyte between two planar dielectric impenetrable interfaces. Specifically, the model contains two planar interfaces located at z = 0 and z = Lz and mobile ions confined between. Each interface is functionalized by uniformly distributed acidic groups (which are fixed during the simulation). The ionization reaction of an acid monomer HA is given by HA ⇐⇒ H ^+ + A^− , where HA and A − are the protonated (non-ionized) and deprotonated (ionized) forms of the weak acid. The monomers at the interface may either have a negative charge −e (A − ) or may be neutral (HA), and the number of protons H + in the solution is adjusted to retain neutrality. The combination of reaction ensemble method and image charge-ELC, which are both implemented in espressomd, would allow me to study dielectric effects on this chemical process. 

I have been checking that the system I create is indeed initially charged neutral, whereas the espressomd told me: ERROR: ELC does not work for non-neutral systems and non-metallic dielectric contrast. I am a new user of espressomd and have been trapped here---I would appreciate anyone who helps!

With kind regards,

Jiaxing 

Below is my script:

from __future__ import print_function

from __future__ import division

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 espressomd

f_config = open("config.txt", 'w+') 

box_lxy = 20.0

box_lz= 20.0 #distance between planar interfaces

system = espressomd.System(box_l=[box_lxy, box_lxy, 1.5*box_lz])

l_bjerrum = 2.0

temperature = 1.0

system.set_random_state_PRNG()

np.random.seed(seed=system.seed)

system.cell_system.skin = 2.5

system.cell_system.max_num_cells = 274400

system.thermostat.set_langevin(kT=temperature, gamma=100.0)

system.time_step = 0.001

# Particle setup

N0 = 2 # 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)

# setting up the A-

c = 0

for i in range(N0):

    for j in range(N0):

        system.part.add(pos=[i*box_lxy/N0, j*box_lxy/N0, 1.0], 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-1.0], id=c,

        type=type_A, q=charges[type_A], mass=1, fix=[1,1,1])

        c = c + 1

# setting up H+

for i in range(N0*N0*2):

    system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy, 0.5+np.random.random()*(box_lz-1)],

    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

# setting up LJ-interactions

lj_eps = 1.0

lj_sig = 1.0

lj_cut = 1.122462048

lj_shift = 1.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 electrostatics

p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=1e-4)

system.actors.add(p3m)

elc = electrostatic_extensions.ELC(gap_size=box_lz*0.5, maxPWerror=1e-4,

    delta_mid_top=0.8, delta_mid_bot=0.8)

system.actors.add(elc)

# Setting up reaction

mode = "reaction_ensemble" #"constant_pH_ensemble"

K_diss = 0.000000001 #0.002694 smaller K_diss, more HA

K_w = 10.**(-14)*0.02694**2

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 = 6

# 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]})

# water autoprotolysis

RE.add_reaction(gamma=(1.0/K_w), reactant_types=[type_H, type_OH], reactant_coefficients=[

                1, 1], product_types=[], product_coefficients=[], 

                default_charges={type_H: charges[type_H], type_OH: charges[type_OH]})

print(RE.get_status())

RE.set_volume (box_lxy*box_lxy*box_lz)

RE.set_wall_constraints_in_z_direction (0,box_lz)

system.setup_type_map([type_HA, type_A, type_H, type_OH, type_Na, type_Cl])

alpha = []

nHA = []

nA = []

nH = []

nOH = []

# thermalization

print("position\n", system.part[:].pos)

print ("charge\n", system.part[:].q)

print("ID\n", system.part[:].id)

print("type\n", system.part[:].type)

print("fix\n", system.part[:].fix)

for i in range(100):

    system.time_step = 0.001

    system.force_cap = 0.001

    RE.reaction()

    system.integrator.run(200)

    print(i, " HA", system.number_of_particles(type=type_HA), "A-",

        system.number_of_particles(type=type_A), "H+", 

        system.number_of_particles(type=type_H),'OH-', 

        system.number_of_particles(type=type_OH), 'Cl-', 

        system.number_of_particles(type=type_Cl), 'NA+', 

        system.number_of_particles(type=type_Na))

    print("position\n", system.part[:].pos)

    print ("charge\n", system.part[:].q)

    print("ID\n", system.part[:].id)

    print("type\n", system.part[:].type)

    print("fix\n", system.part[:].fix)

# production run begins here

system.time_step = 0.01

system.force_cap = 0.0

print("<production run begins>\n")

for i in range(100000):

    RE.reaction()

    system.integrator.run(200) 

    print(i, " HA", system.number_of_particles(type=type_HA), "A-",

        system.number_of_particles(type=type_A), "H+", 

        system.number_of_particles(type=type_H),'OH-', 

        system.number_of_particles(type=type_OH), 'Cl-', 

        system.number_of_particles(type=type_Cl), 'NA+', 

        system.number_of_particles(type=type_Na))

    print(i, "\n", file=f_config)

    print("type\n", system.part[:].type, file=f_config)

    print("position\n", system.part[:].pos, file=f_config)

    alpha.append(system.number_of_particles(type=type_A)/N0)

    nHA.append(system.number_of_particles(type=type_HA))

    nA.append(system.number_of_particles(type=type_A))

    nH.append(system.number_of_particles(type=type_H))

    nOH.append(system.number_of_particles(type=type_OH))

# data analysis

alpha_av = np.mean(alpha)

alpha_err = np.std(alpha) / np.sqrt(len(alpha))

print("\n<alpha> = {} (err = {})\n".format(alpha_av, alpha_err))

nHA_av = np.mean(nHA)

nA_av = np.mean(nA)

nH_av = np.mean(nH)

nOH_av = np.mean(nOH)

print("\n<nHA> = {} ".format(nHA_av))

print("\n<nA>  = {} ".format(nA_av))

print("\n<nH>  = {} ".format(nH_av))

print("\n<nOH> = {} ".format(nOH_av))