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

[Konrad Breitsprecher]

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)