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:
ESPResSo features a number of basic shapes like cylinders, walls or spheres to simulate confined systems. Here, we use two walls at $z = 0$ and $z = L_z$ for the parallel plate setup ($L_z$: box length in $z$-direction)
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.
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.
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 [1].
Let's walk through the python script. We need additional imports for the wall shapes and the ELC algorithm:
from espressomd import System, 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:
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
k_B = 1.380649e-23 # units of [J/K]
q_e = 1.602176634e-19 # units of [C]
epsilon_0 = 8.8541878128e-12 # units of [C^2/J/m]
coulomb_prefactor = q_e**2 / (4 * numpy.pi * epsilon_0) * 1e10
l_bjerrum = 0.885**2 * coulomb_prefactor / (k_B * temp)
wall_margin = 0.5
Ez = 0
num_steps_equilibration = 200
num_configs = 50
integ_steps_per_config = 100
['ADDITIONAL_CHECKS', 'AFFINITY', 'BMHTF_NACL', 'BOND_CONSTRAINT', 'BUCKINGHAM', 'COLLISION_DETECTION', 'CUDA', 'DIPOLAR_BARNES_HUT', 'DIPOLAR_DIRECT_SUM', 'DIPOLES', 'DP3M', 'DPD', 'EK_BOUNDARIES', 'ELECTROKINETICS', 'ELECTROSTATICS', 'ENGINE', 'EXCLUSIONS', 'EXPERIMENTAL_FEATURES', 'EXTERNAL_FORCES', 'FFTW', 'GAUSSIAN', 'GAY_BERNE', 'GSL', 'H5MD', 'HAT', 'HERTZIAN', 'LANGEVIN_PER_PARTICLE', 'LB_BOUNDARIES', 'LB_BOUNDARIES_GPU', 'LB_ELECTROHYDRODYNAMICS', 'LENNARD_JONES', 'LENNARD_JONES_GENERIC', 'LJCOS', 'LJCOS2', 'LJGEN_SOFTCORE', 'MASS', 'MEMBRANE_COLLISION', 'METADYNAMICS', 'MMM1D_GPU', 'MORSE', 'NPT', 'P3M', 'PARTICLE_ANISOTROPY', 'ROTATION', 'ROTATIONAL_INERTIA', 'SMOOTH_STEP', 'SOFT_SPHERE', 'TABULATED', 'THOLE', 'VIRTUAL_SITES', 'VIRTUAL_SITES_INERTIALESS_TRACERS', 'VIRTUAL_SITES_RELATIVE', 'WCA']
We save the force field parameters in python dictionaries, now with parameters for the walls:
# 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.
# 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 = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3
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.
# 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"])
<espressomd.constraints.ShapeBasedConstraint at 0x1540740b71d0>
Now we place the particles at random position without overlap with the walls:
# 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:
# 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, followed by the thermostat:
energy = system.analysis.energy()
print("Before Minimization: E_total = {:.3e}".format(energy['total']))
system.minimize_energy.init(f_max=10, gamma=10, max_steps=2000,
max_displacement=0.01)
system.minimize_energy.minimize()
energy = system.analysis.energy()
print("After Minimization: E_total = {:.3e}".format(energy['total']))
# Set thermostat
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)
Before Minimization: E_total = 3.272e+15 After Minimization: E_total = -7.134e+05
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.
# 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.66872e+00 7.78356e+01 7.071e-03 7.784e+01 accuracy not achieved
12 7 4.90514e-01 7.66872e+00 1.70554e+01 7.071e-03 1.706e+01 accuracy not achieved
16 7 4.90514e-01 7.66872e+00 9.28496e-01 7.071e-03 9.285e-01 accuracy not achieved
18 7 4.90514e-01 7.66872e+00 2.79314e-01 7.071e-03 2.792e-01 accuracy not achieved
22 7 4.90514e-01 7.66872e+00 4.24312e-02 7.071e-03 4.184e-02 accuracy not achieved
26 7 4.90514e-01 7.66872e+00 1.21072e-02 7.071e-03 9.828e-03 accuracy not achieved
28 7 4.76143e-01 7.90432e+00 9.97014e-03 7.071e-03 7.029e-03 16.95
28 6 4.90514e-01 7.66872e+00 2.74754e-02 7.071e-03 2.655e-02 accuracy not achieved
32 7 4.19415e-01 8.99349e+00 9.99017e-03 7.071e-03 7.057e-03 14.69
32 6 4.76143e-01 7.90432e+00 1.42755e-02 7.071e-03 1.240e-02 accuracy not achieved
34 7 3.97298e-01 9.50319e+00 9.91693e-03 7.071e-03 6.953e-03 18.69
34 6 4.19415e-01 8.99349e+00 2.31351e-02 7.071e-03 2.203e-02 accuracy not achieved
resulting parameters: mesh: (32 32 44), cao: 7, r_cut_iL: 4.1942e-01,
alpha_L: 8.9935e+00, accuracy: 9.9902e-03, time: 14.69
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:
for p in system.part:
p.ext_force = [0, 0, Ez * p.q]
This is followed by our standard temperature equilibration:
# 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("progress={:.0f}%, t={:.1f}, E_total={:.2f}, E_coulomb={:.2f}, T={:.4f}"
.format(i * 100. / int(num_steps_equilibration / 100 - 1), system.time,
energy['total'], energy['coulomb'], temp_measured), end='\r')
system.integrator.run(100)
print()
progress=100%, t=0.2, E_total=-18474879.29, E_coulomb=-32691922.01, T=7146.8495
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.
# 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):
print('progress: {:>3.0f}%'.format(i * 100. / num_configs), end='\r')
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
print('progress: 100%')
progress: 100%
Finally, we calculate the average, normalize the data with the bin volume and save it to a file using NumPy's savetxt command.
# 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.
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$-axis $(\\mathrm{\\AA})$', fontsize=20)
plt.ylabel('Density $(\\mathrm{u\\AA}^{-3})$', 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).
[1] 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