Ionic Conductivity

In this tutorial we will compute the ionic conductivity of sodium chloride ions in an aqueous solution by means of molecular dynamics simulations with the ForceField engine.

The system and workflow presented here is very similar to the one described in the publication Estimates of Electrical Conductivity from Molecular Dynamics Simulations: How to Invest the Computational Effort J. Phys. Chem. B 124, 9680−9689 (2020).

Ionic conductivity (S/m) of 1.21 mol/L NaCl in water at 300K

AMS

7.3

Ref. (MD)

5.9 \(\pm\) 2.9

The tutorial consists of the following steps:

  • Building a sodium chloride solution with AMBER95 atom types.

  • Equilibrating the system at 300 K

  • Running a simulation in the NVT ensemble

  • Obtaining a value for the ionic conductivity

  • Comparing that value to literature

If you wish, you can download the already-relaxed sodium chloride solution NaCl_wat.in and jump directly to the diffusion coefficient calculation

/scm-uploads/doc.trunk/Tutorials/_images/IonicConductivityNaCl.png

Building the NaCl solution

1. Open a new AMSinput window
2. Switch to ForceField: ADFPanel → ForceFieldPanel
3. Select Type: → Amber95

First we create the building blocks for our system from scratch.

3. Draw a single Na atom
4. In the main panel, Click on MoreBtn next to Atom types
5. In the column FF Type fill in the AMBER atom type IP
6. In the column FF Charge fill in 1.0
7. File → Export coordinates → .in and give it the name “Na.in”
/scm-uploads/doc.trunk/Tutorials/_images/Na.png

  • Repeat for a Cl atom (FF Type IM, FF Charge -1.0), and save as “Cl.in”.

  • Repeat for a single water molecule (FF Types O: OW, H: HW, FF Charges O: -0.8340, H: 0.4170), and save as “h2o.in”.

Now we are ready to build the full system.

1. Remove all atoms from the AMSinput window
2. Click on Edit → Builder
3. Make a box of 19x19x19 Angstrom by typing 19 into the diagonal elements of the Lattice vectors
4. Click on the folder icon next to copies of:, and import the previously saved file “Na.in”
5. Change the 100 copies into 5 copies
6. Click on AddButton Molecules: to add a new molecule
7. Import the previously saved file “Cl.in”, and change to 5 copies
8. Add a third molecule, import the previously saved file “h2o.in”, and change to 220 copies
9. Click the Generate Molecules button
10. Close the Builder tool by clicking the Close button

This should result in an aqueous NaCl solution with a density of 1.029 g/mL.

/scm-uploads/doc.trunk/Tutorials/_images/builder.png

Important

While building the system, make sure that the Amber95 force field is still selected. Changing the ForceField Type after building will result in all atom typing information being erased.

Equilibrating the system

First perform a geometry optimization.

1. In the main panel, select Task → Geometry Optimization
2. Periodicity → Bulk
3. File → Save As… and give it an appropriate name (e.g. “NaCl_optimization”)
4. Run the calculation
5. Click yes when asked if the structure in AMSinput should be updated

Equilibrate the resulting system at the desired temperature

1. In the main panel, select Task → Molecular Dynamics
2. Click on MoreBtn next to Task: Molecular Dynamics to go to the MD details
3. Set the number of steps to 20000
4. Set the Time step to 0.5 fs
5. Set the Initial temperature to 300 K
/scm-uploads/doc.trunk/Tutorials/_images/equilib12.png
6. Click on MoreBtn next to Thermostat to go to the thermostat details
7. Set the Type to Berendsen
8. Set the Temperature to 300 K
9. Set the damping constant to 100 fs
/scm-uploads/doc.trunk/Tutorials/_images/equilib3.png

You are now ready to run the simulation.

1. File → Save As… and give it an appropriate name (e.g. “NaCl_equilibration”)
2. File → Run
3. You can follow the progress of the calculation in AMSMovie and AMStail (SCM → Logfile)

AMSMovie can plot the potential energy of the system over the course of the simulation.

1. In AMSMovie select Properties → Potential Energy
2. In AMSMovie select Properties → Temperature
/scm-uploads/doc.trunk/Tutorials/_images/equilibrated.png

The potential energy of the system has converged to a stable value over the course of the simulation, and the temperature oscillates around 300K.

Running the Production Simulation

When the equilibration calculation has finished, you may be asked whether to update coordinates and bonds. If not, close the NaCl_equilibration AMSinput window, select the job in AMSjobs, and re-open it in AMSinput:

Tick the use MD velocities from results for AMS MD restart box
Select Yes, new job

We will run the production simulation using a Nose-Hoover chain thermostat.

1. In the main panel, select Task → Molecular Dynamics
2. Click on MoreBtn next to Task: Molecular Dynamics to go to the MD details
3. Set the number of steps to 100000
4. Set the Sample frequency to 20

If you ran the equilibration yourself:

5. Set initial velocities to From File and set File to the ams.rkf file from the NaCl_equilibration job (if it isn’t set already)

Otherwise, if you simply imported the provided equilibrated structure:

5. Set initial velocities to Random and set Temperature to 300

Then continue with the thermostat settings:

6. Click on MoreBtn next to Thermostat to go to the thermostat details
7. Set the Type to NHC
8. Set the Temperature to 300 K
9. Set the damping constant to 100 fs

You are now ready to run the simulation.

1. File → Save As… and give it an appropriate name (e.g. “NaCl_production”)
2. File → Run
3. You can follow the progress of the calculation in AMSMovie and AMStail (SCM → Logfile)

Computing the Ionic Conductivity

The ionic conductivity is proportional to the slope of the charge weighted mean square displacement of the atoms.

\(\sigma = \lim_{t\to\infty} \frac{1}{2tdVk_BT} \sum_i^N q_i^2 \langle ({\bf{r_i}}(0) - {\bf{r_i}}(t))^2 \rangle\)

Here \(d\) is the dimensionality of the systems (usually \(d=3\)), and \(N\) is the number of ions.

To compute it, open the completed simulation in AMSMovie.

1. In AMSMovie select MD Properties → MSD
2. Select Property → Conductivity
3. Set Atoms to Na Cl
4. Click on Generate MSD
/scm-uploads/doc.trunk/Tutorials/_images/msd.png

Two new curves have appeared in the AMSMovie window. The bottom curve depicts the convergence of the ionic conductivity with the correlation time. The final value of the ionic conductivity is 7.34 S/m (or C²/Jms). The corresponding literature states that the value should fall between 1.2 and 13.5 S/m, with an average over many simulations of 5.9 S/m.

/scm-uploads/doc.trunk/Tutorials/_images/msd_graph.png

The top curve depicts the charge weighted mean square displacement curve, of which the ionic conductivity is the slope. The graph shows that at short correlation times the slope of the curve is large, and then gradually becomes constant. The point where the curve becomes linear roughly corresponds to the time it takes for the motions of the ions to decorrelate. AMSMovie attempts to estimate this point from the shape of the mean square displacement curve. Indeed, the bottom graph starts at around 7ps, which is AMSMovies guess. However, a brief look at the top curve suggests that a better guess for the decorrelation time would be around 10 ps. We can recompute the ionic conductivity with this new decorrelation time in mind.

1. In AMSMovie select MD Properties → MSD
2. Keep all settings the same as before
3. Set Start time slope to 10000.0
4. Click on Replace MSD
/scm-uploads/doc.trunk/Tutorials/_images/msd10000.png

The new value for the ionic conductivity is 7.29 S/m.

Important

The example simulation presented here is too short to obtain really reliable results. In practice it is advised to perform multiple simulations of several nanoseconds each.

Important

The publication J. Phys. Chem. B 124, 9680−9689 (2020) discusses how the exact, correlated, ionic conductivity would include off-diagonal contributions to the mean square displacement.

\(\sigma_{corr} = \lim_{t\to\infty} \frac{1}{2tdVk_BT} \left( \sum_i^N q_i^2 \langle ({\bf{r_i}}(0) - {\bf{r_i}}(t))^2 \rangle + \sum_{i \neq j}^N q_i q_j \langle ({\bf{r_i}}(0) - {\bf{r_i}}(t)) ({\bf{r_j}}(0) - {\bf{r_j}}(t)) \rangle \right)\)

In this example we used high NaCl concentrations to be able to compare with literature values. Indeed, at these concentrations, the off-diagonal contributions can change the value of \(\sigma\) by approximately a factor two. At lower ion concentrations the effect of the off-diagonal contributions will become negligible.

Important

AMSMovie uses the atomic charges as they are stored on the ams.rkf file (in this case the forcefield charges). In this example, the forcefield charges match the formal charges, but for larger ions, this may not be the case, and this will affect the results.