espressomd-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [ESPResSo-users] electrostatic layer correction in espressomd


From: Jiaxing Yuan
Subject: Re: [ESPResSo-users] electrostatic layer correction in espressomd
Date: Sat, 6 Apr 2019 16:19:43 +0800

Dear Konrad.

Yes, it works now! Thanks a lot!

Best,
Jiaxing 

On Thu, Apr 4, 2019 at 9:40 PM Konrad Breitsprecher <address@hidden> wrote:
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 Breitsprecher

Am 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_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 = 10
box_lz= 5 #distance between planar interfaces
fraction_empty = 1.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.005)
#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.5
system.cell_system.max_num_cells = 274400
# Particle setup
N = 10 # 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])
system.force_cap = 500
for 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 = 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)

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,
Konrad

Am 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)

reply via email to

[Prev in Thread] Current Thread [Next in Thread]