<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tutorial 2
Step1:  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
Step2:  We save the force field parameters in python dictionaries, now with parameters for the walls
Step3:  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.
Step4:  In the next snippet, we add the walls to our system. Our constraint takes two arguments
Step5:  Now we place the particles at random position without overlap with the walls
Step6:  The scheme to set up the Lennard-Jones interaction is the same as before, 
extended by the Electrode-Ion interactions
Step7:  Next is the Lennard-Jones Equilibration
Step8:  As described, we use P$^3$M in combination with ELC to account for the 2D-periodicity. 
ELC is also added to the <tt>actors</tt> of the system and takes gap size and maximum
pairwise errors as arguments.
Step9:  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
Step10:  This is followed by our standard temperature equilibration
Step11:  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.
Step12:  Finally, we calculate the average, normalize the data with the bin volume and save it to
a file using NumPy's <tt>savetxt</tt> command.
Step13:  Finally we can plot the density of the ions.
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function
from espressomd import System, assert_features, electrostatics, electrostatic_extensions
from espressomd.shapes import Wall
import espressomd
import numpy
Explanation: 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:
End of explanation
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
Explanation: 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:
End of explanation
# 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}
Explanation: We save the force field parameters in python dictionaries, now with parameters for the walls:
End of explanation
# 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)
Explanation: 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.
End of explanation
# 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"])
Explanation: In the next snippet, we add the walls to our system. Our constraint takes two arguments: 
First the <tt>shape</tt>, in our case a simple plane defined by its normal vector and the distance from the origin, 
second the <tt>particle_type</tt>, which is used to set up the interaction between particles and constraints.
End of explanation
# 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"])
Explanation: Now we place the particles at random position without overlap with the walls:
End of explanation
# 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")
Explanation: The scheme to set up the Lennard-Jones interaction is the same as before, 
extended by the Electrode-Ion interactions:
End of explanation
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'])
Explanation: 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.
End of explanation
# 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)
Explanation: As described, we use P$^3$M in combination with ELC to account for the 2D-periodicity. 
ELC is also added to the <tt>actors</tt> of the system and takes gap size and maximum
pairwise errors as arguments.
End of explanation
for p in system.part:
    p.ext_force = [0,0,Ez * p.q]
Explanation: 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:
End of explanation
# 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)
Explanation: This is followed by our standard temperature equilibration:
End of explanation
# 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
Explanation: 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.
End of explanation
# 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)")
Explanation: Finally, we calculate the average, normalize the data with the bin volume and save it to
a file using NumPy's <tt>savetxt</tt> command.
End of explanation
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()
Explanation: Finally we can plot the density of the ions.
End of explanation
<END_TASK>