Tutorial 2: A Simple Charged System, Part 2

7 2D Electrostatics and Constraints

In this section, we use the parametrized NaCl system from the last task to simulate a molten salt in a parallel plate capacitor with and without applied electric field. We have to extend our simulation by several aspects:

Confinement

ESPResSo features a number of basic shapes like cylinders, walls or spheres to simulate confined systems. Here, we use to walls at $z = 0$ and $z = L_z$ for the parallel plate setup ($L_z$: Box length in z-direction)

2D-Electrostatics

ESPResSo also has a number of ways to account for the unwanted electrostatic interaction in the now non-periodic z-dimension. We use the 3D-periodic P$^3$M algorithm in combination with the Electrostatic Layer Correction (ELC). ELC subtracts the forces caused by the periodic images in the z-dimension. Another way would be to use the explicit 2D-electrostatics algorithm MMM2D, also available in ESPResSo.

Electric Field

The simple geometry of the system allows us to treat an electric field in z-direction as a homogeneous force. Note that we use inert walls here and don't take into account the dielectric contrast caused by metal electrodes.

Parameters

For our molten NaCl, we use a temperature $100 \ \mathrm{K}$ above the melting point ($1198.3 \ \mathrm{K}$) and an approximated density of $\rho = 1.1138 \ \mathrm{u \mathring{A}}$$^{-3}$ found in Janz, G. J., Thermodynamic and Transport Properties of Molten Salts: Correlation Equations for Critically Evaluated Density, Surface Tension, Electrical Conductance, and Viscosity Data, J. Phys. Chem. Ref. Data, 17, Suppl. 2, 1988.

Let's walk through the python script. We need additional imports for the wall shapes and the ELC algorithm:

In [1]:
from __future__ import print_function
from espressomd import System, assert_features, electrostatics, electrostatic_extensions
from espressomd.shapes import Wall
import espressomd
import numpy

If we target a liquid system, we should not set up the particles in a lattice, as this introduces unwanted structure in the starting configuration. We define our system size by the number of particles and the density. The system parameters lead to the following values:

In [2]:
required_features = ["EXTERNAL_FORCES", "MASS", "ELECTROSTATICS", "LENNARD_JONES"]
espressomd.assert_features(required_features)
print(espressomd.features())

# System parameters
n_part = 1000
n_ionpairs = n_part/2
density = 1.1138
time_step = 0.001823
temp = 1198.3
gamma = 50
#l_bjerrum = 0.885^2 * e^2/(4*pi*epsilon_0*k_B*T)
l_bjerrum = 130878.0 / temp
wall_margin = 0.5
Ez = 0

num_steps_equilibration = 3000
num_configs = 200
integ_steps_per_config = 100
['ELECTROSTATICS', 'EXTERNAL_FORCES', 'FFTW', 'LATTICE', 'LB', 'LB_BOUNDARIES', 'LENNARD_JONES', 'MASS', 'P3M', 'ROTATION', 'ROTATIONAL_INERTIA', 'SWIMMER_REACTIONS', 'VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']

We save the force field parameters in python dictionaries, now with parameters for the walls:

In [3]:
# Particle parameters
types       = {"Cl":          0, "Na": 1, "Electrode": 2}
numbers     = {"Cl": n_ionpairs, "Na": n_ionpairs}
charges     = {"Cl":       -1.0, "Na": 1.0}
lj_sigmas   = {"Cl":       3.85, "Na": 2.52,  "Electrode": 3.37}
lj_epsilons = {"Cl":     192.45, "Na": 17.44, "Electrode": 24.72}

lj_cuts     = {"Cl":        3.0 * lj_sigmas["Cl"], 
               "Na":        3.0 * lj_sigmas["Na"],
               "Electrode": 3.0 * lj_sigmas["Electrode"]}

masses      = {"Cl":  35.453, "Na": 22.99, "Electrode": 12.01}

To finally calculate the box size, we take into account the diameter of the electrode interaction. Additionally, ELC needs a particle-free gap in the z-direction behind the wall.

In [4]:
# Setup System
box_l = (n_ionpairs * sum(masses.values()) / density)**(1. / 3.)
box_z = box_l + 2.0 * (lj_sigmas["Electrode"]+wall_margin)
elc_gap = box_z * 0.15
system = System(box_l=[box_l,box_l, box_z+elc_gap])
system.seed=42
box_volume = numpy.prod([box_l,box_l, box_z])

system.periodicity = [1, 1, 1]
system.time_step = time_step
system.cell_system.skin = 0.3
system.thermostat.set_langevin(kT=temp, gamma=gamma)

In the next snippet, we add the walls to our system. Our constraint takes two arguments: First the shape, in our case a simple plane defined by its normal vector and the distance from the origin, second the particle_type, which is used to set up the interaction between particles and constraints.

In [5]:
# Walls   
system.constraints.add(shape=Wall(dist=wall_margin,
            normal=[0,0,1]),particle_type=types["Electrode"])
system.constraints.add(shape=Wall(dist=-(box_z-wall_margin),
            normal=[0,0,-1]),particle_type=types["Electrode"])
Out[5]:
<espressomd.constraints.ScriptInterfaceHelper at 0x110cbd910>

Now we place the particles at random position without overlap with the walls:

In [6]:
# Place particles
for i in range(int(n_ionpairs)):
    p = numpy.random.random(3)*box_l
    p[2] += lj_sigmas["Electrode"]
    system.part.add(id=len(system.part), 
                    type=types["Cl"],  pos=p, q=charges["Cl"], mass=masses["Cl"])
for i in range(int(n_ionpairs)):
    p = numpy.random.random(3)*box_l
    p[2] += lj_sigmas["Electrode"]
    system.part.add(id=len(system.part), 
                    type=types["Na"],  pos=p, q=charges["Na"], mass=masses["Na"])

The scheme to set up the Lennard-Jones interaction is the same as before, extended by the Electrode-Ion interactions:

In [7]:
# Lennard-Jones interactions parameters 

def combination_rule_epsilon(rule, eps1, eps2):
    if rule=="Lorentz":
        return (eps1*eps2)**0.5
    else:
        return ValueError("No combination rule defined")

def combination_rule_sigma(rule, sig1, sig2):
    if rule=="Berthelot":
        return (sig1+sig2)*0.5
    else:
        return ValueError("No combination rule defined")

for s in [["Cl", "Na"], ["Cl", "Cl"], ["Na", "Na"], ["Na", "Electrode"], ["Cl", "Electrode"]]:
        lj_sig = combination_rule_sigma("Berthelot", 
                lj_sigmas[s[0]], lj_sigmas[s[1]])
        lj_cut = combination_rule_sigma("Berthelot", 
                lj_cuts[s[0]], lj_cuts[s[1]])
        lj_eps = combination_rule_epsilon("Lorentz", 
                lj_epsilons[s[0]],lj_epsilons[s[1]])

        system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")

Next is the Lennard-Jones Equilibration: Here we use an alternative way to get rid of the overlap: ESPResSo features a routine for energy minimization with similar features as in the manual implementation used before. Basically it caps the forces and limits the displacement during integration.

In [8]:
energy = system.analysis.energy()
print("Before Minimization: E_total=", energy['total'])
system.minimize_energy.init(f_max = 10, gamma = 10, max_steps = 2000,
                            max_displacement= 0.1)
system.minimize_energy.minimize()
energy = system.analysis.energy()
print("After Minimization: E_total=", energy['total'])
Before Minimization: E_total= 1.2455422308e+14
After Minimization: E_total= 44312.6671337

As described, we use P$^3$M in combination with ELC to account for the 2D-periodicity. ELC is also added to the actors of the system and takes gap size and maximum pairwise errors as arguments.

In [9]:
# Tuning Electrostatics
p3m = electrostatics.P3M(prefactor=l_bjerrum*temp, 
        accuracy=1e-2)
system.actors.add(p3m)
elc = electrostatic_extensions.ELC(gap_size = elc_gap, 
        maxPWerror = 1e-3)
system.actors.add(elc)
P3M tune parameters: Accuracy goal = 1.00000e-02 prefactor = 1.30878e+05 
System: box_l = 3.16243e+01 # charged part = 1000 Sum[q_i^2] = 1.00000e+03
mesh cao r_cut_iL     alpha_L      err          rs_err     ks_err     time [ms]
10   7   4.90514e-01 7.66873e+00 7.78359e+01 7.071e-03 7.784e+01 accuracy not achieved
12   7   4.90514e-01 7.66873e+00 1.70555e+01 7.071e-03 1.706e+01 accuracy not achieved
16   7   4.90514e-01 7.66873e+00 9.28500e-01 7.071e-03 9.285e-01 accuracy not achieved
18   7   4.90514e-01 7.66873e+00 2.79315e-01 7.071e-03 2.792e-01 accuracy not achieved
22   7   4.90514e-01 7.66873e+00 4.24314e-02 7.071e-03 4.184e-02 accuracy not achieved
26   7   4.90514e-01 7.66873e+00 1.21077e-02 7.071e-03 9.828e-03 accuracy not achieved
28   7   4.76143e-01 7.90433e+00 9.97031e-03 7.071e-03 7.029e-03 17.85   
28   6   4.90514e-01 7.66873e+00 2.74752e-02 7.071e-03 2.655e-02 accuracy not achieved
32   7   4.19415e-01 8.99349e+00 9.98926e-03 7.071e-03 7.056e-03 18.77   
32   6   4.76143e-01 7.90433e+00 1.42753e-02 7.071e-03 1.240e-02 accuracy not achieved
34   7   3.97298e-01 9.50320e+00 9.87072e-03 7.071e-03 6.887e-03 21.98   
34   6   4.19415e-01 8.99349e+00 2.31349e-02 7.071e-03 2.203e-02 accuracy not achieved

resulting parameters:
28   28   40   7   4.76143e-01 7.90433e+00 9.97031e-03 17.85   

For now, our electric field is zero, but we want to switch it on later. Here we run over all particles and set an external force on the charges caused by the field:

In [10]:
for p in system.part:
    p.ext_force = [0,0,Ez * p.q]

This is followed by our standard temperature equilibration:

In [11]:
# Temperature Equilibration
system.time = 0.0
for i in range(int(num_steps_equilibration/100)):
    energy = system.analysis.energy()
    temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
    print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f}, T={3:.4f}".format(system.time,
                                                                            energy['total'], energy['coulomb'], temp_measured), end='\r')
    system.integrator.run(100)
t=5.3, E_total=-33559397.48, E_coulomb=-37024055.11, T=1185.7435

In the integration loop, we like to measure the density profile for both ion species along the z-direction. We use a simple histogram analysis to accumulate the density data. Integration takes a while.

In [12]:
# Integration
bins=100
z_dens_na = numpy.zeros(bins)
z_dens_cl = numpy.zeros(bins)
system.time = 0.0
cnt = 0

for i in range(num_configs):
    energy = system.analysis.energy()
    temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
    system.integrator.run(integ_steps_per_config)

    for p in system.part:
        bz = int(p.pos[2]/box_z*bins)
        if p.type == types["Na"]:
            z_dens_na[bz] += 1.0
        elif p.type == types["Cl"]:
            z_dens_cl[bz] += 1.0
    cnt += 1

Finally, we calculate the average, normalize the data with the bin volume and save it to a file using NumPy's savetxt command.

In [13]:
# Analysis

# Average / Normalize with Volume
z_dens_na /= (cnt * box_volume/bins)
z_dens_cl /= (cnt * box_volume/bins)
z_values = numpy.linspace(0,box_l,num=bins)
res = numpy.column_stack((z_values,z_dens_na,z_dens_cl))
numpy.savetxt("z_density.data",res,
        header="#z rho_na(z) rho_cl(z)")

Finally we can plot the density of the ions.

In [15]:
import matplotlib.pyplot as plt
plt.ion()

plt.figure(figsize=(10,6), dpi=80)
plt.plot(z_values,z_dens_na, label='Na')
plt.plot(z_values,z_dens_cl, label='Cl')
plt.xlabel('$z$', fontsize=20)
plt.ylabel('Density', fontsize=20)
plt.legend(fontsize=16)
plt.show()

The resulting density plot is very noisy due to insufficient sampling, but should show a slight depletion of the smaller Na atoms at the walls. Now try to put in an electric field that represents an applied voltage of $15 \ \mathrm{V}$ between the walls and compare the results. The density data should show strong layering at the walls, decaying towards the system center. The complete script is at /doc/tutorials/02-charged_system/scripts/nacl_units_confined.py. In the interactive script nacl_units_confined_vis.py, you can increase/decrease the electric field with the keys u/j (at your own risk).

missing
Figure 4: Snapshot and densities along the z-Axis with applied electric field for the ion species.