|
| From: | Konrad Breitsprecher |
| Subject: | Re: [ESPResSo-users] ICC* |
| Date: | Fri, 12 Jul 2019 12:55:19 +0200 |
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
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,Jiaxingfrom __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):
#system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy, 1+np.random.random()*(box_lz-2)],
system.part.add(pos=[0.2692+3,-1.9375+3,2.4000+12],
id=c, type=type_cat, q=charges[type_cat], mass=1)
c = c + 1
for j in range(N):
#system.part.add(pos=[np.random.random()*box_lxy, np.random.random()*box_lxy, 1+np.random.random()*(box_lz-2)],
system.part.add(pos=[0.0586+3,-0.7526+3,-2.4000+12],
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.actors.add(p3m)
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.actors.add(icc)
# 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)
| [Prev in Thread] | Current Thread | [Next in Thread] |