Nanoconfined electrolyte

Molecular dynamics simulation of an aqueous NaCl solution confined between two solid walls


Responsive image


The objective of this tutorial is to use molecular dynamics and simulate an electrolyte confined between two rigid carbon walls. There are several main parts to this tutorial:

  • System generation - First, the electrolyte and the two walls are generated using LAMMPS internal commands.
  • Minimisation - Seconds, an energy minimisation is performed to bring the system in a more favorable configuration.
  • System equilibration - Third, the system is equilibrated at a pressure of 1 atm and temperature of 300 K.
  • Going further with exercice - Finally, some exercices are proposed (just send me an email if you want the solutions, or register here).

If you are new to LAMMPS, I recommend you to follow this simpler tutorial first.

Click here if you are looking for help with your project, if you want to support me (for free or not), or if you have any suggestions for these tutorials.



System generation

Create a new folder that will contain all the files. Open a blank page using a text editor, and save it in the new folder using the name 'input.01.lammps'. Note that the '.lammps' extension has no importance. Copy the following lines into input.01.lammps:

# Initialisation
units		real
atom_style	full
bond_style	harmonic
angle_style	harmonic
pair_style	lj/cut/tip4p/long 1 2 1 1 0.1546 12.0
kspace_style	pppm/tip4p 1.0e-4

Explanations: There are many differences with respect to tutorial 01. With the unit style 'real', masses are in grams per mole, distances in Ångstroms, time in femtoseconds, energies in Kcal/mole. With the atom style 'full', each atom is a dot with a mass and a charge. In addition, each atom can be linked by bonds, angles, dihedrals and impropers potentials (for example to form molecules). The 'bond_style' and 'angle_style' commands define the styles of bond and angle used in the simulation, respectively, and the 'harmonic' keyword imposes the potential to use.

Note: with the rigid water model that is used here, the bond and angle styles are not relevant, but LAMMPS requires them to be specified.

With the 'pair_style' named 'lj/cut/tip4p/long', atoms interact through both a Lennard-Jones (LJ) potential and through Coulombic interactions. This style is specific to four points water model, and automatically accounts for the additional massless site. The six numbers are, respectively,

  • 1 - the atom type for the oxygen O of the tip4p water,
  • 2 - the atom type for the hydrogen H of the tip4p water,
  • 3 - the OH bond type,
  • 4 - the HOH angle type,
  • 5 - the distance from O atom to massless charge,
  • 6 - the cutoff in Ångstroms.

The cutoff applies to both LJ and Coulombic interactions, but in a different way. For LJ 'cut' interactions, atoms interact with each other only if they are separated by a distance smaller than the cutoff. For Coulombic 'long', interaction between atoms closer than the cutoff are computed directly, and interaction between atoms outside that cutoff are computed in reciprocal space. Finally the kspace defines the long-range solver for Coulombic interactions. The pppm style refers to particle-particle particle-mesh. From Luty and van Gunsteren' paper:

The PPPM method is based on separating the total interaction between particles into the sum of short-range interactions, which are computed by direct particle−particle summation, and long-range interactions, which are calculated by solving Poisson's equation using periodic boundary conditions (PBCs).

Let us create the box: copy the following lines into input.01.lammps, after the 'kspace_style' command:

# ------------- System definition
lattice 	fcc 4.04
region		box block -5 5 -5 5 -13 13
create_box	6 box &
		bond/types 1 &
		angle/types 1 &
		extra/bond/per/atom 2 &
		extra/angle/per/atom 1 &
		extra/special/per/atom 2

Explanations: The 'lattice' command defines the unit cell. Here 'fcc' with a scale factor of 4.04 has been chosen for the future positioning of the carbon atoms in a cubic crystal lattice. The 'region' command defines a geometric region of space, and by choosing 'xlo=-5' and 'xhi=5', and because we have previously chosen a lattice with scale factor of 4.04, the region 'box' extends from -20.2 to 20.2 Ångström (units of lengths are in Ångstrom because we have chosen the unit style 'real'). Finally, the 'create_box' command creates a simulation box with 6 types of atoms in the simulation: oxygen and hydrogen atoms of the water, Na, Cl, and the carbon atoms of the top and bottom walls, respectively. Note that the carbon atoms of the top and bottom walls will be identical (same mass, same pairwise interaction, etc.), and the assignment of two different numbers is for practical purpose.

The create box command extends over 6 lines thanks to the '&' character. The second and third lines are used to specify that the simulation with contain 1 type of bond and 1 type of angle (for the water molecule). The parameters of these bond and angle constraints will be given later. The three last lines are for memory allocation.

We can now add atoms to the system. First, we create two sub-regions corresponding respectively to the two solid walls. Then we create the atoms of type 5 and 6 within the two regions, respectively.

# ------------- create the walls
region		rbotwall block -5 5 -5 5 -12 -10
region		rtopwall block -5 5 -5 5 10 12

create_atoms 	5 region rtopwall
create_atoms 	6 region rbotwall

In order to add the water molecules, we first need to download the file TIP4P2005.txt and place it in the same folder. It contains all the necessary information about the water molecule, such as positions, bonds, and angle. Then, add the following lines to input.01.lammps:

# ------------- create the fluid
region		rliquid block -5 5 -5 5 -9 9
lattice		sc 4.04
molecule 	h2omol TIP4P2005.txt
create_atoms 	0 region rliquid mol h2omol 482793

Explanations: With the last four lines, a region used to deposit the water is created on the last defined lattice, which is 'fcc 4.04'. Then, on the next line, we define a new simple cubic lattice in order to position the water molecules on it, with a distance of 4.04 Ångstroms between each water molecule (note: the new lattice replaces the previous one, LAMMPS reads a script from top to bottom). Note that 4.04 Ångstroms is larger than the typical equilibrium distance between water molecules in a liquid, but this will allow us to insert ions more safely. Then, the 'molecule' command opens the 'TIP4P2005.txt' file, and name the associated molecule 'h2omol'. Finally, molecules are created on the sc lattice by the 'create_atoms' command. The first parameter is '0' because we use the atom id from the 'TIP4P2005.txt' file. The number '482793' is a seed that is required by LAMMPS, it can be any positive integer. Finally, let us deposit 20 ions (10 Na, 10 Cl) in between the water molecules by adding these two lines to input.01.lammps:

fix		mydep1 all deposit 10 3 1 56513 region rliquid near 0.3
fix		mydep2 all deposit 10 4 1 58613 region rliquid near 0.3

Each 'fix deposit' will add one ion at a random position within the 'rliquid' region every timestep. So we can just make a very short simulation of 10 timesteps and export the generated configuration. Feel free to increase the salt concentration by increasing both number of steps and number of desired ions in the 'deposit' fixes. Note that 'mydep1' and 'mydep2' are the name I choose to give to the fix.

We need to define the parameters of the simulation: the mass of the 6 atoms (O, H, Na, Cl, Al, Al), the pairwise interaction parameters (here the parameters for the Lennard-Jones potential), and the bond and angle parameters. Copy the following line into input.01.lammps:

# ------------- settings
include		PARM.lammps

Create a new text file, call it 'PARM.lammps', and copy it in the same folder where 'input.01.lammps' is. Copy the following lines into PARM.lammps:

mass 		1 15.9994
mass 		2 1.008
mass		3 28.990
mass		4 35.453
mass		5 26.9815
mass		6 26.9815

pair_coeff 	1 1 0.185199 3.1589
pair_coeff 	2 2 0.0 0.0
pair_coeff	3 3 0.04690 2.4299
pair_coeff	4 4 0.1500 4.04470
pair_coeff	5 5 11.697 2.574
pair_coeff	6 6 11.697 2.574

bond_coeff 	1 0 0.9572

angle_coeff 	1 0 104.52

Explanations: The parameters for water (mass 1, mass 2, pair_coeff 1 1, pair coeff 2 2, bond_coeff 1 and angle_coeff 1) are given by the TIP4P/2005 force field, the parameters for Na\(^+\) and Cl\(^-\) (mass 3, mass 4, pair_coeff 3 3, pair_coeff 4 4) are given by the CHARMM-27 force field, and the parameters for the wall (mass 5, mass 6, pair_coeff 5 5, pair coeff 6 6) are parameters for Al atoms. Each 'mass' command assigns a mass in grams/mole to an atom type. Each 'pair_coeff' assigns respectively the depth of the potential in Kcal/mole, and the distance at which the particle-particle potential energy in Ångstrom. We have only assigned pairwise interaction between atoms of identical type. By default, LAMMPS calculates the pair coefficients for the interactions between atoms of type i and j by using the geometrical rule. Other rules can be set with the 'pair_modify' command, but for the sake of simplicity, we are going to keep the default option here.

The bond coefficient (here for the O-H bond of the water molecule) sets respectively the energy of the harmonic potential and the equilibrium distance in Ångstrom. The value is '0' for the energy, because we are going to use a rigid model for the water molecule. The shape of the molecule will be preserved by the shake algorithm (see later). Similarly, the angle coefficient (here for the H-O-H angle of the water molecule) sets the energy of the harmonic potential (also 0) and the equilibrium angle is in degree.

Finally, add the following line to input.01.lammps:


dump 		mydmp all image 10 dump.*.jpg type type

# ------------- run
run		10

set 		type 3 charge 1.0
set 		type 4 charge -1.0

write_data	data.01.lammps

Explanations: The 'run 10' requires the simulation to run for 10 timesteps. The value of the timestep (1 fs by default) does not matter yet because the atoms are not moving. We also need to specify the charge of the newly added ions, which is done using the 'set' commands. The write 'data_file' finally creates a file named 'data.01.lammps' containing all the information required to restart the simulation from the final configuration generated by input.01.lammps. The 'dump' command with the 'image' option will save 'jpg' images of the system.

The first input script is ready to be executed with LAMMPS. In the terminal windows, you should see something like that:

LAMMPS (20 Nov 2019)
Lattice spacing in x,y,z = 4.04 4.04 4.04
Created orthogonal box = (-20.2 -20.2 -52.52) to (20.2 20.2 52.52)
  1 by 1 by 1 MPI processor grid
Created 1000 atoms
  create_atoms CPU = 0.00466238 secs
Created 1000 atoms
  create_atoms CPU = 0.000538191 secs
Lattice spacing in x,y,z = 4 4 4

(...)

Per MPI rank memory allocation (min/avg/max) = 35.92 | 35.92 | 35.92 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0            0 8.5279913e+20            0 8.5279913e+20 1.3643111e+21 
      10            0 8.5279913e+20            0 8.5279913e+20 1.3643111e+21 
Loop time of 1.65591 on 1 procs for 10 steps with 8917 atoms

Performance: 0.522 ns/day, 45.998 hours/ns, 6.039 timesteps/s
99.2% CPU use with 1 MPI tasks x no OpenMP threads

(...)

Neighbor list builds = 10
Dangerous builds = 0
Setting atom values ...
  10 settings made for charge
Setting atom values ...
  10 settings made for charge
System init for write_data ...

(...)

Total wall time: 0:00:01

If the simulation is successfully, a file named 'data01.lammps' should be in the folder, as well as the following snapshot of the system:



Responsive image





Energy minimisation

Why is this step necessary?

It is clear from the way the system has been created that the atoms are not at equilibrium distances from each others. Indeed, some of the ions added using the 'fix deposit' commands are too close to the water molecules. If we were to start a 'normal' molecular dynamics simulation now (i.e. solve the equations of motion) with a 'normal' timestep (1 or 2 femto-seconds), the atoms would exert huge forces on each others, as a consequence they would accelerate brutally, and the simulation would fail.

MD simulations failing due to overlapping atoms is extremely common. If it occurs, you can either (1) delete the overlapping atoms or (2) move the atoms to more reasonable distances (here we do the second one).

Here we need to find a way to move the atoms and place them in a more favorable position before starting the simulation. This step is called 'energy minimization', and is often necessary.

Implementation

To perform the energy minimization with our system, let us open a new blank sheet in the same folder, and rename it input.02.lammps. The first lines will be very similar to the previous input file:

# ------------- Initialisation
boundary	p p p
units		real
atom_style	full
bond_style	harmonic
angle_style	harmonic
pair_style	lj/cut/tip4p/long 1 2 1 1 0.1546 12.0
kspace_style	pppm/tip4p 1.0e-4

# ------------- System definition
read_data	data.01.lammps

# ------------- settings
include		PARM.lammps

The only difference with input.01.lammps is that, instead of creating a box and atoms, we open and read the previously created file 'data.01.lammps' which contains the definition of the simulation box and the positions of the atoms. Next, we are going to create some groups with the different atom types:

group		gH2O type 1 2
group		gNa type 3
group		gCl type 4
group		gliquid type 1 2 3 4
group		gtopwall type 5
group		gbotwall type 6
group		gwall type 5 6
neigh_modify 	exclude group gwall gwall
dump 		mydmp all atom 1000 dump.02.lammpstrj

Creating groups allows us to apply different dynamics to the liquid and to the walls. The 'neigh_modify' command asks to LAMMPS not to evaluate interactions between the atoms of the walls. This command is used here because internal deformations of the walls will not be permitted (rigid walls). The 'neigh_modify' command is not necessary, but save computation time. The last line asks LAMMPS to generate a dump file with the atoms positions.

Let us force the carbon walls to remain rigid during the simulation:

# ------------- run
# ------------- walls
fix		mysf1 gtopwall setforce 0 0 NULL
fix		mysf2 gbotwall setforce 0 0 NULL
fix		myaf1 gtopwall aveforce NULL NULL 0
fix		myaf2 gbotwall aveforce NULL NULL 0

The first fix 'setforce' applies to the group 'gtopwall', which contains all of the atoms of type 5. This fix cancels the \(x\) and \(y\) components of the forces applied on the atoms of the group at each timestep. Therefore, if these atoms have no initial velocity (which is the case here), they wont move along \(x\) anf \(y\). These fixes do nothing to the \(z\) components thanks to the 'NULL' keyword. The first fix 'aveforce' applies to the group 'gtopwall', and averages all the force exerted on the atoms of the group over \(z\). As a consequence, the atoms of the group 'gtopwall' move as a block along \(z\) (and the same for the bottom wall).

Now, we can include the most important commands for the minimization:

fix		mynve all nve
compute		tliq gliquid temp
fix		myber gliquid temp/berendsen 1 1 1
fix_modify	myber temp tliq

The fix 'nve', which we apply to all atoms, performs constant NVE integration to update positions and velocities of the atoms at each timestep. The 'temp/berendsen' fix rescales the velocities of the atoms of the group liquid (ions + water) every timestep in order to reset the temperature. Since we want to perform a minimization step, both initial and final temperatures have been chosen equal to 1K. The third parameter is the damping factor, in time units, which determines how rapidly the temperature is relaxed. A damping factor of 1 fs would be too small for a molecular dynamics simulation, but is acceptable for a minimization step during which we just want the atoms to move slightly from their initial positions. The 'fix_modify' is used to assign the temperature of the group 'gliquid' as calculated by the compute 'tliq' to the Berendsen thermostat.

Without this 'fix_modify', LAMMPS would use the measured temperature of group 'all' as the reference temperature when rescaling the velocities of the atoms of group 'gliquid', which is inconsistent. The simulation would still work without error message, but the temperature of the atoms would be different from the desired value.

If we were to run the simulation as it is, it would fail because nothing maintains the shape of the water molecules (the bond and angle energies are equal to 0). Let us use the shake algorithm in order to apply bond and angle constraints to the water molecules. In addition, let us add a fix 'recenter' in order to maintain the system centered in the middle of the box in the \(z\) direction. The fix recenter has no influence on the dynamics.

fix		myshk gH2O shake 1.0e-4 200 0 b 1 a 1
fix		myrct all recenter NULL NULL INIT

Finally, let us choose a very small timestep (because we anticipate that the atoms are initially too close to each others) and run for 1000 timesteps (with the command thermo 1000, thermodynamic info are printed in the terminal every 1000 timesteps).

timestep	0.1
thermo		500
run		5000
write_data	data.02.lammps

Result

When running the input.02.lammps file, you should see that the total energy of the system (fifth column) decreases toward negative values:

Step Temp E_pair E_mol TotEng Press 
       0            0    105943.88            0    105943.88    178061.22 
     500     3.161202   -7033.8737            0   -6979.0416   -3364.4997 
    1000    2.2057861   -8866.4103            0   -8828.1502   -4212.2734 
    1500    1.7782626   -10038.926            0   -10008.081   -4763.0287 
    2000    1.6006587   -10847.999            0   -10820.235   -5152.8961 
    2500    1.5852741   -11446.485            0   -11418.988     -5447.37 
    3000    1.7011608   -11917.953            0   -11888.445   -5681.5579 
    3500    1.9695876   -12329.433            0    -12295.27   -5891.3157 
    4000     2.445561   -12722.314            0   -12679.894   -6098.0218 
    4500    3.2010334   -13037.534            0   -12982.011   -6250.5859 
    5000    4.3296591   -13339.374            0   -13264.275   -6393.1841 

If you look at the trajectory using VMD, you will see some of the atoms (the one that where initially in a problematic position) slightly move from each others, as seen in this video.



System equilibration

Starting from the final configuration of the minimization step, we can now run the molecular dynamics simulation. Here we are going to let the system evolve freely until it reaches an equilibrium. Let us create a new input file and call it input.03.lammps. Add the following lines:

# ------------- Initialisation
boundary	p p p
units		real
atom_style	full
bond_style	harmonic
angle_style	harmonic
pair_style	lj/cut/tip4p/long 1 2 1 1 0.1546 12.0
kspace_style	pppm/tip4p 1.0e-4

# ------------- System definition
read_data	data.02.lammps

# ------------- Simulation settings
include		PARM.lammps
group		gH2O type 1 2
group		gNa type 3
group		gCl type 4
group		gliquid type 1 2 3 4
group		gtopwall type 5
group		gbotwall type 6
group		gwall type 5 6
dump 		mydmp all atom 1000 dump.03.lammpstrj

# ------------- Run

# ------------- To maintain the wall rigid
fix		mysf1 gtopwall setforce 0 0 NULL
fix		mysf2 gbotwall setforce 0 0 NULL
fix		myaf1 gtopwall aveforce NULL NULL 0
fix		myaf2 gbotwall aveforce NULL NULL 0
neigh_modify 	exclude group gwall gwall

# ------------- Dynamics
fix		mynve all nve
compute		tliq gliquid temp
fix		myber gliquid temp/berendsen 300 300 100
fix_modify	myber temp tliq
fix		myshk gH2O shake 1.0e-4 200 0 b 1 a 1
fix		myrct all recenter NULL NULL INIT

timestep	1.0
thermo		5000
thermo_modify 	temp tliq
run		50000
write_data	data.03.lammps

The two main differences with the previous step are the timestep (1 fs instead of 0.1 fs) and the thermostating. With an imposed temperature of 300 K, the fluid is expected to behave as a liquid. A timestep of 1 fs is classic for these force fields. Never use a larger timestep except if you have a good reason to do so. What you observe should resembles this video.

The 'thermo_modify' ensures that the printed temperature is the temperature of the liquid, and not the temperature of the whole system which is lower (because the walls are frozen and have an effective temperature of 0 K)

Running the last input script is likely to take time (more than 30 minutes). You can speed up the calculation by running LAMMPS on multiple cores using mpirun.

All the scripts and data files generated during this tutorial are available here.



Going further with exercises


Request the solutions by email, or register here and access all the solutions + additional LAMMPS content.


Exercise 1 : Extract data from the simulation

Starting from the equilibrated configuration (data.03.lammps), extract the density profiles (for Na\(^+\), Cl\(^-\), and water) along the axis normal to the solid walls.

Responsive image

Hint: data extraction can in general be done (1) using the internal LAMMPS command, and (2) using a post-processing analysis tool (e.g. Python).

Exercise 2 : Use fully flexible walls

Starting from the equilibrated configuration (data.03.lammps), or starting from scratch, do the same using fully flexible wall instead of frozen walls. Flexible walls can be thermalised, and are therefore more realistic.

What you observe should resembles this video.

Exercise 3 : Simulate the transport properties

Starting from the equilibrated configuration (data.03.lammps), induce a flow of the liquid in the direction tangential to wall, and measure the velocity profile.

Advice 1: when inducing a flow in a direction \(i\), one needs to exclude the degrees of freedom in the direction \(i\) from the thermalisation, otherwise the resulting temperature is likely to be wrong.
Advice 2: the forcing must be chosen with care. If too large, the system will be strongly out-of-equilibrium. If too small, no net velocity will be measured because of the noise.

Is the velocity profile you obtain consistent with Poiseuille's predictions?



Click here if you are looking for help with your project, if you want to support me (for free or not), or if you have any suggestions for these tutorials.