|
From: | Konrad Breitsprecher |
Subject: | Re: [ESPResSo-users] electrostatic layer correction in espressomd |
Date: | Thu, 4 Apr 2019 15:40:19 +0200 |
Dear Konrad,Hello! Thank you so much for the suggestion. I also anticipate this is due to some particles move out of the ELC region. Now I use the below script (I try a smaller system for a quick checking) where I first equilibrate the system using LJ potential and after that turn on electrostatics. The simulation is running normally when I set both dielectric mismatches to be 0, but If I change one of them to be even 0.01 (so 0.001 and 0), then again the error of the non-neutral system. I do not think a dielectric mismatch 0.001 will induce too large image force? I am confused and I am not really sure if the below script is modified according to your ideas.Best,Jiaxing==================================from __future__ import print_functionfrom __future__ import divisionimport espressomdfrom espressomd import code_infofrom espressomd import analyzefrom espressomd import integratefrom espressomd.interactions import *from espressomd import reaction_ensemblefrom espressomd import interactionsfrom espressomd import electrostaticsfrom espressomd import electrostatic_extensionsfrom espressomd.shapes import Wallimport sysimport numpy as npimport timebox_lxy = 10box_lz= 5 #distance between planar interfacesfraction_empty = 1.0Acc = 1e-5system = espressomd.System(box_l=[box_lxy, box_lxy,(1+fraction_empty)*box_lz],periodicity = [1,1,1],time_step=0.005)#system.cell_system.set_domain_decomposition(use_verlet_lists=False)#system.cell_system.set_layered(n_layers=1)l_bjerrum = 3.5temperature = 1.0system.thermostat.set_langevin(kT=temperature, gamma=100.0)system.set_random_state_PRNG()np.random.seed(seed=system.seed)system.cell_system.skin = 0.5system.cell_system.max_num_cells = 274400# Particle setupN = 10 # number of salt moleculetype_cat = 0type_ani = 1type_sub = 2charges = {}charges[type_cat] = 1charges[type_ani] = -1charges[type_sub] = 0# Wallssystem.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 ionsc = 0for 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)],id=c, type=type_cat, q=charges[type_cat], mass=1)c = c + 1for 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)],id=c, type=type_ani, q=charges[type_ani], mass=1)c = c + 1# setting up LJ-interactionslj_eps = 1.0lj_sig = 1.0lj_cut = 1.122462048lj_shift = 0.0for 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])system.force_cap = 500for i in range(1500):system.integrator.run(200)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)p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=Acc)system.actors.add(p3m)elc = electrostatic_extensions.ELC(gap_size=box_lz*fraction_empty, maxPWerror=Acc, delta_mid_top=0, delta_mid_bot=0)system.actors.add(elc)system.force_cap = 0for i in range(1000):system.integrator.run(200)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)print("\n<production run begins>\n")f_config = open("config7.txt", 'w+')for i in range(100000):system.integrator.run(200)energy = system.analysis.energy()temp_measured = energy['kinetic'] / ((3.0 / 2.0)* (len(system.part.select(type=0))+len(system.part.select(type=1))))number_of_particles = system.number_of_particles(type=0) + system.number_of_particles(type=1)print(i, 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)On Wed, Apr 3, 2019 at 7:48 PM Konrad Breitsprecher <address@hidden> wrote:Hello Jiaxing,in your script you do the warmup (with force capping) after setting up the electrostatics. That is generally a bad strategy, you get a more stable behaviour by equilibrating your LJ interaction first, add the electrostatics and equilibrate again. I could cure your script (at least no immediate error as before) by simply running the warmup with the force capping before adding p3m.The error was a bit misleading indeed and is a result of the order of the sanity checks in espresso. Most likely, a particle slipped through your constraint into the ELC-Gap region and the actual (valid) system volume became non-neutral (which is not allowed with dielectric contrast != +/- 1). The ELC sanity check fired before the actual error (particle X violated constraint) could be recognised.Hope this helps,KonradAm Mi., 3. Apr. 2019 um 08:28 Uhr schrieb Jiaxing Yuan <address@hidden>:Dear all,
I am simulating a system of electrolytes (a set of cations and anions) between two planar dielectric interfaces. I want to use the electrostatic layer correction to deal with the surface polarization. The below is my script for this system, the very strange thing is if I set the dielectric mismatch to be the same, say 0.939, everything is working fine. But if I change one of the dielectric mismatches to be slightly different, say 0.93 and 0.939, then an error comes to me: ELC does not work for non-neutral systems and non-metallic dielectric contrast. Why does this happen? Thank you for any suggestions.
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 = 100
box_lz= 50 #distance between planar interfaces
fraction_empty = 3.0
Acc = 1e-5
system = espressomd.System(box_l=[box_lxy, box_lxy,(1+fraction_empty)*box_lz],periodicity = [1,1,1],time_step=0.01)
#system.cell_system.set_domain_decomposition(use_verlet_lists=False)
#system.cell_system.set_layered(n_layers=1)
l_bjerrum = 3.5
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.3
system.cell_system.max_num_cells = 274400
# Particle setup
N = 218 # 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)
# setting up ions
c = 0
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)],
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)],
id=c, type=type_ani, q=charges[type_ani], mass=1)
c = c + 1
# setting up LJ-interactions
lj_eps = 1.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])
p3m = electrostatics.P3M(prefactor=l_bjerrum * temperature, accuracy=Acc)
system.actors.add(p3m)
elc = electrostatic_extensions.ELC(gap_size=box_lz*fraction_empty, maxPWerror=Acc, delta_mid_top=0.93, delta_mid_bot=0.939)
system.actors.add(elc)
system.force_cap = 100
for i in range(1000):
system.integrator.run(200)
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)
system.force_cap = 0
for i in range(1000):
system.integrator.run(200)
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)
print("\n<production run begins>\n")
f_config = open("config7.txt", 'w+')
for i in range(100000):
system.integrator.run(200)
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0)* (len(system.part.select(type=0))+
len(system.part.select(type=1))))
number_of_particles = system.number_of_particles(type=0) + system.number_of_particles(type=1)
print(i, 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)
[Prev in Thread] | Current Thread | [Next in Thread] |