Ionic Conductivity
******************

ForceField Batteries Amber95 Ionic conductivity Molecular Dynamics
Aqueous solution

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

   [image]


Building the NaCl solution
==========================

   **1.** Open a new **AMSinput** window
   **2.** Switch to ForceField: [image: ADFPanel][image] **→** [image: ForceFieldPanel][image]
   **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 [image: MoreBtn][image] 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 (System) → .in** and give it the name "Na.in"

   [image]

* 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*, both H
  atoms: *HW*, **FF Charges** O: -0.8340, H: 0.4170), and save as
  "H2O.in".

-[ Na.in ]-

   System
     Atoms
       Na 0.0 0.0 0.0 ForceField.Type=IP ForceField.Charge=1.0
     End
   End

-[ Cl.in ]-

   System
     Atoms
       Cl 0.0 0.0 0.0 ForceField.Type=IM ForceField.Charge=-1.0
     End
   End

-[ H2O.in ]-

   System
     Atoms
       O 0.0 0.0 0.40468149 ForceField.Type=OW ForceField.Charge=-0.834
       H 0.78258644 0.0 -0.20234075 ForceField.Type=HW ForceField.Charge=0.417
       H -0.78258644 0.0 -0.20234075 ForceField.Type=HW ForceField.Charge=0.417
     End
     BondOrders
       1 2 1.0
       1 3 1.0
     End
   End

Now we are ready to build the full system.

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.

   **1.** Remove all atoms from the **AMSinput** window
   **2.** Double-check that the **Type** on the force field panel is set to **Amber95**
   **3.** Click on **Builders → Packed Molecules**
   **4.** Set **Density** to 1.029 g/cm³.
   **5.** Set **Box bounds** to "19" "19" "19" Å.
   **6.** Tick **Use number of molecules**
   **7.** Click on the folder icon next to **SMILES or file**, and import the previously saved file "Na.in"
   **8.** Set **N mols** to "5"
   **9.** Click on [image: AddButton][image] to add a new molecule
   **10.** Import the previously saved file "Cl.in", and set **N mols** to "5"
   **11.** Click on [image: AddButton][image] to add a new molecule
   **12.** Import the previously saved file "H2O.in", and leave **N mols** empty. The appropriate number of molecules will be automatically determined.
   **13.** Click the **Generate Molecules** button
   **14.** 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, containing 5 Na⁺, 5 Cl⁻, and 220 water molecules.

Note:

  *New in AMS2026*: The possibility to automatically determine the
  number of water molecules was introduced in AMS2026. For previous
  AMS versions, set the number of water molecules to 220 and do not
  set the density.

   [image]


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.** When the calculation is finished, siwtch to AMSinput and click on [image: UpdateMolecule][image]

Equilibrate the resulting system at the desired temperature

   **1.** In the main panel, select **Task → Molecular Dynamics**
   **2.** Click on [image: MoreBtn][image] 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

   [image]

   **6.** Click on [image: MoreBtn][image] 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

   [image]

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.

   In **AMSMovie** select **MD Properties → Temperature**

   [image]

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**:

   **1.** Switch back to **AMSinput**
   **2.** In the main panel, click on [image: UpdateMolecule][image]

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

   **1.** In the main panel, select **Task → Molecular Dynamics**
   **2.** Click on [image: MoreBtn][image] 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 [image: MoreBtn][image] 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**)

Wait for the calculation to finish (this may take around 45 minutes,
depending on your hardware).


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**

   [image]

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. Your results may differ by
  a factor of 3 to the results presented here.

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.

   [image]

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.** Tick **Use → All values**
   **4.** Set **Start time slope** to **10000.0**
   **5.** Click on **Replace MSD**

   [image]

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

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.
