espressomd-users
[Top][All Lists]

## Re: [ESPResSo-users] ICC*

 From: Jiaxing Yuan Subject: Re: [ESPResSo-users] ICC* Date: Fri, 12 Jul 2019 23:28:50 +0800

Thank you for your help. Using the same script as yours, I get a completely different answer for the energy when I turn on ICC*. What version of EspressoMD are you using? If there any special points to note when compiling the ICC* into EspressoMD? Thank you very much in advance.

Coulomb Energy Ion Pair:  0.017193882552404505

Coulomb Energy Ion Pair + ICC:  -0.44437159139925037

Difference:  0.46156547395165487

Best,
Jiaxing

On Fri, Jul 12, 2019 at 6:55 PM Konrad Breitsprecher <address@hidden> wrote:
Hello Jiaxing,

how did you compare 'with and without' ICC? If I use the script below, I get the right value (0.01729) and no significant difference in energy for ion pair vs. ion pair + ICC. In the script, I use the following protocol:
setup ICC particles with zero charge -> setup ion pair -> activate p3m -> measure energy ion pair -> set icc start charges -> activate icc -> measure energy ion pair + ICC.

Cheers,

The test script:

import espressomd
from espressomd import electrostatics
from espressomd import electrostatic_extensions
import numpy as np

box_lxy = 6
box_lz= 6 #distance between planar interfaces
fraction_empty = 3.0
Acc = 1e-5
Acc_ICC = 1e-12
system = espressomd.System(box_l=[box_lxy, box_lxy,(1+fraction_empty)*box_lz],periodicity = [1,1,1],time_step=0.0000000000000000001)
# Set the ICC line density and calculate the number of
# ICC particles according to the box size
l = 0.25
nicc = int(box_lxy / l)
nicc_per_interface = nicc * nicc
nicc_tot = 2 * nicc_per_interface
iccArea = box_lxy * box_lxy / nicc_per_interface
l = box_lxy / nicc
l_bjerrum = 1.0
temperature = 1.0
#system.thermostat.set_langevin(kT=temperature, gamma=100.0)
system.set_random_state_PRNG()
np.random.seed(seed=system.seed)
system.cell_system.skin = 0.5
system.cell_system.max_num_cells = 274400
# Lists to collect required parameters
iccNormals = []
iccAreas = []
iccSigmas = []
iccEpsilons = []
# Particle setup
N = 1 # number of salt molecule
type_cat = 0
type_ani = 1
type_sub = 2
charges = {}
charges[type_cat] = 1
charges[type_ani] = -1
charges[type_sub] = 0

qs_icc = 0.0

c = 0
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, 9], q=qs_icc, fix=[1, 1, 1], type=type_sub)
c = c + 1
iccNormals.append([0,0,1])
# Right interface (normal [0,0,-1])
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, 15], q=-qs_icc, fix=[1, 1, 1], type=type_sub)
c = c + 1
iccNormals.append([0,0,-1])
# setting up ions
for i in range(N):
id=c, type=type_cat, q=charges[type_cat], mass=1)
c = c + 1
for j in range(N):
id=c, type=type_ani, q=charges[type_ani], mass=1)
c = c + 1

iccAreas.extend([iccArea] * nicc_tot)
iccSigmas.extend([0] * nicc_tot)
iccEpsilons.extend([1.00000100000] * nicc_tot)

p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=Acc)

system.integrator.run(0)
ECoulombIonPair = system.analysis.energy()['coulomb']

# Set icc start charges
c = 0
for xi in range(nicc):
for yi in range(nicc):
system.part[c].q=0.001
c = c + 1
for xi in range(nicc):
for yi in range(nicc):
system.part[c].q=-0.001
c = c + 1

icc = electrostatic_extensions.ICC(first_id=0, n_icc=nicc_tot, convergence=Acc_ICC,
relaxation=1.5, ext_field=[0, 0, 0], max_iterations=20000, eps_out=1.0, normals=iccNormals,
areas=iccAreas, sigmas=iccSigmas, epsilons=iccEpsilons)

system.integrator.run(0)
ECoulombIonPairICC = system.analysis.energy()['coulomb']

print('Coulomb Energy Ion Pair: ', ECoulombIonPair)
print('Coulomb Energy Ion Pair + ICC: ', ECoulombIonPairICC)
print('Difference: ', ECoulombIonPair - ECoulombIonPairICC)

Relevant output:

resulting parameters: mesh: (14 14 54), cao: 7, r_cut_iL: 4.1504e-01,

alpha_L: 7.2532e+00, accuracy: 9.9176e-06, time: 1.82

Coulomb Energy Ion Pair:  0.01729305574554396

Coulomb Energy Ion Pair + ICC:  0.017292442816883573

Difference:  6.129286603887008e-07

Am Fr., 12. Juli 2019 um 07:50 Uhr schrieb Jiaxing Yuan <address@hidden>:
Dear all,

I am using ICC* in EspressoMD for dielectric interfaces. However, I find a strange thing in the code. Using the attached script below (2 ions between two planar dielectric interfaces), if I just turn off the calculation due to dielectric mismatch by excluding the line "system.actors.add(icc)" and set those surface patches to have zero charge, EspressoMD gives a coulomb energy 0.0171938825524045. This coincides with other code. But when I turn on the calculation of ICC*, even if I set the dielectric mismatch to be extremely tiny (for instance iccEpsilons.extend([1.00000100000] * nicc_tot)), then the code gives Coulomb energy 0.013045563701536.
I don't think this is reasonable because physically the two energies should be as close as possible, given that I have set such a tiny dielectric mismatch and thus the surface polarization charge should be zero. Thank you so much for the clarification!

Best,
Jiaxing

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 = 6
box_lz= 6 #distance between planar interfaces
fraction_empty = 3.0
Acc = 1e-5
Acc_ICC = 1e-8
system = espressomd.System(box_l=[box_lxy, box_lxy,(1+fraction_empty)*box_lz],periodicity = [1,1,1],time_step=0.0000000000000000001)
# Set the ICC line density and calculate the number of
# ICC particles according to the box size
l = 0.25
nicc = int(box_lxy / l)
nicc_per_interface = nicc * nicc
nicc_tot = 2 * nicc_per_interface
iccArea = box_lxy * box_lxy / nicc_per_interface
l = box_lxy / nicc
l_bjerrum = 1.0
temperature = 1.0
#system.thermostat.set_langevin(kT=temperature, gamma=100.0)
system.set_random_state_PRNG()
np.random.seed(seed=system.seed)
system.cell_system.skin = 0.5
system.cell_system.max_num_cells = 274400
# Lists to collect required parameters
iccNormals = []
iccAreas = []
iccSigmas = []
iccEpsilons = []
# Particle setup
N = 1 # number of salt molecule
type_cat = 0
type_ani = 1
type_sub = 2
charges = {}
charges[type_cat] = 1
charges[type_ani] = -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)
# Add the fixed ICC particles:
# Left interface (normal [0,0,1])
c = 0
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, 9], q=0.001, fix=[1, 1, 1], type=type_sub)
c = c + 1
iccNormals.append([0,0,1])
# Right interface (normal [0,0,-1])
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, 15], q=-0.001, fix=[1, 1, 1], type=type_sub)
c = c + 1
iccNormals.append([0,0,-1])
# setting up ions
for i in range(N):
id=c, type=type_cat, q=charges[type_cat], mass=1)
c = c + 1
for j in range(N):
id=c, type=type_ani, q=charges[type_ani], mass=1)
c = c + 1
iccAreas.extend([iccArea] * nicc_tot)
iccSigmas.extend([0] * nicc_tot)
iccEpsilons.extend([1.00000100000] * nicc_tot)
p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=Acc)
icc = electrostatic_extensions.ICC(first_id=0, n_icc=nicc_tot, convergence=Acc_ICC,
relaxation=1.5, ext_field=[0, 0, 0], max_iterations=20000, eps_out=1.0, normals=iccNormals,
areas=iccAreas, sigmas=iccSigmas, epsilons=iccEpsilons)
# setting up LJ-interactions
lj_eps = 0.0
lj_sig = 1.0
lj_cut = 1.122462048
lj_shift = 0.0
for i in range(2):
for j in range(2):
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(2):
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)
system.setup_type_map([type_cat, type_ani, type_sub])
for i in range(1):
system.integrator.run(1)
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0)* (len(system.part.select(type=0))+
len(system.part.select(type=1))))
print(i, energy['total'], energy['coulomb'], energy["kinetic"], temp_measured)
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)