|
From: | Jiaxing Yuan |
Subject: | Re: [ESPResSo-users] electrostatic layer correction in espressomd |
Date: | Sat, 6 Apr 2019 16:19:43 +0800 |
Ok actually I was wrong. Due to a bug, the neutrality is checked for the sum of particle and induced charges.I send you instructions to quickly fix it yourself in your espresso code or wait for the updated development version.Regards,Konrad BreitsprecherAm Do., 4. Apr. 2019 um 14:01 Uhr schrieb Jiaxing Yuan <address@hidden>: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] |