Viscosity of an ionic liquid with APPLE&P¶

Creating the system¶

Let’s start by creating the atomic structure with the Builder.

PYR14 = from_smiles("CCCC[N+]1(CCCC1)C")
TFSI = from_smiles("C(F)(F)(F)S(=O)(=O)[N-]S(=O)(=O)C(F)(F)F")

box = appleandp_packmol(
    molecules=[PYR14, TFSI],
    n_molecules=[30, 30],
    density=0.5,
    region_names=["cation", "anion"],
    forcefield_file=forcefield_file
)

We intentionally chose not to fill the box to the experimental density of 1367 kg/m³, we will adjust the density later.

Setting up the calculation¶

Next we need to set up the APPLE&P engine.

GUIPython (PLAMS)

1. Select the ForceField engine by switching to the Force Field panel: →

2. Set Periodicity → Bulk.

3. Set Type → APPLE&P.

4. Click on to go to the Model → APPLE&P panel.

5. Set formal charges for each molecule found in the structure as -1 for C2F6NO4S2 and +1 for C9H20N.

6. Click the Generate button which will run a tool to generate a system-specific set of APPLE&P parameters.

Pre-optimization¶

It’s always a good idea to perform a quick optimization of the system.

GUIPython (PLAMS)

From the Main tab.

1. Select Task → Geometry Optimization.

2. Click on to go to the Details → Geometry Optimization panel.

3. Enable the Optimize lattice check-box.

4. Set the Convergence to Basic.

You can now save the input as appViscosity_preopt (for example), and run the pre-optimization.

Adjusting the simulation box¶

We would like now to adjust the simulation box dimensions in order to obtain a density close to the experimental density of 1367 kg/m³. Click View → Properties to display various information about the system. The current density is approximately 491 kg/m³. In order to reach a reasonable density, we need to shrink the simulation box. We will therefore perform molecular dynamics under deformation to reduce the simulation box dimensions to the target size. Let’s first define the MD settings:

GUIPython (PLAMS)

1. From the Main tab, select Task → Molecular Dynamics.

2. Click on to go to the Model → Molecular Dyanmics panel.

3. Set the total Number of steps to 10000 and the Time step to 0.5 fs.

4. We also initialize the velocities to the Initial temperature of 333 K and set the Checkpoint frequency to 10000.

Note

With APPLE&P it is recommended to use a timestep smaller than 1 fs.

We now define the deformations. The deformation procedure will linearly adjust the simulation box to the target dimension during the molecular dynamics simulation. A bit of algebra gives us the equation for the final cell dimension \(L_{f} = \sqrt[3]{\frac{\rho_i}{\rho_f}}L_i \approx 25\) Å.

Note

We do not need to be very precise in the final volume as we will later perform an isobaric simulation and let the density of the system equilibrate.

1. Select Model → MD… → Deformations.

2. Add a deformation by clicking the button

3. Set the Target length to 25 25 25 Å.

Save your calculation as appViscosity_scandensity, and run it. It should take about 10 minutes on a modern laptop. Click yes when asked to update the coordinates. You can appreciate the volume change from SCM → Movie. You can check that the density of the system is now close to the experimental value from View → Properties.

Thermalization¶

We will now perform a short molecular dynamics simulation in the canonical ensemble (NVT) at 333 K to thermalize the system.

GUIPython (PLAMS)

1. From the Main panel, select Task → Molecular Dynamics.

2. Click on to go to the Model → MD panel.

3. We will keep the MD settings from the previous deformation run.

4. Click the to go to the thermostat settings.

5. Click the button to add a thermostat.

6. Select the NHC Thermostat, with the Temperature of 333 K, and a Damping constant of 100 fs.

7. Before running the thermalization we need to remove the deformation: go to Model → MD… → Deformations and click the button left to Deformation: #1.

You can now save your input as appViscosity_NVT and run the calculation. Check the box Use MD velocities from result from the pop-up Results found, and click yes to update the coordinates and bonds.

We can now equilibrate the simulation box under isothermal-isobaric conditions (NPT ensemble).

job = AMSNPTJob.restart_from(
    job,
    name="NPT_eq",
    nsteps=10000,
    timestep=0.5,
    samplingfreq=100,
    thermostat="NHC",
    temperature=333,
    tau=100,
    barostat="MTK",
    barostat_tau=1000,
    equal='XYZ',
    writebonds=True,
    writevelocities=False,
    writemolecules=False,
)

Save the input as appViscosity_NPT and run the calculation.

Click SCM → Movie to analyze the last NPT run. Let’s visualize the density with Graph → Add Graph, followed by MD Properties → Density. Ideally one would want to pursue the NPT equilibration to confirm that the density has converged. We found densities between 1336 and 1346 kg/m³ approximately 2% smaller than the experimental density of 1367 kg/m³.

To reduce the stress on the simulation box we will now select a frame with a small pressure. We first add the pressure to the graph with MD Properties → Pressure and select a frame where the pressure is close to zero. Here you also want to select a late frame in order to ensure that the system is equilibrated. You can use the horizontal slide bar to adjust the frame selection. Once the cursor is on a frame with small pressure, click File → Update Geometry In Input. We are now ready to perform the production run.

Production run¶

For the production run we will perform an NVT simulation as follow:

GUIPython (PLAMS)

1. Set the Model → MD… → Barostat to None.

2. From the MD tab, increase the Number of steps to 4000000 (4 million steps = 2 ns)

3. Set the Sample Frequency to 1000 and Checkpoint Frequency to 400000.

4. From Initial velocities select Random and set the Initial temperature to 333 K.

5. Click on to go to the Trajectory Options and enable the Calculate pressure option.

6. From the Details → Expert AMS scroll down to the MolecularDynamics Binlog section and enable Time and Pressure tensor.

Note

The last binlog option will write the pressure tensor at every step of the simulation. It is very important that you enable this or you will not be able to extract the viscosity!

Save the input as appViscosity_prod and run the calculation. You can follow the simulation with SCM → Movie and verify that the initial pressure was close to zero. It can be useful to perform a moving average over the pressure plot to better appreciate the averaged value. To achieve this click Graph → Analysis…, select the Moving Average tab and click the + button below. Here you can modify the window size, for example, 20 steps, and click the OK button.

Computing viscosity¶

The viscosity of the electrolyte can be calculated using the Green-Kubo relation:

where \(\eta\) is the viscosity, \(V\) is the cell volume, \(k_{\rm{B}}\) is the Boltzmann constant, \(T\) is the temperature, \(P_{\alpha\beta}\) is an off-diagonal component of the pressure tensor, and \(t\) is time. The function \(\langle P_{\alpha\beta}(0)P_{\alpha\beta}(t) \rangle\) is an autocorrelation function with the average taken over all time origins and off-diagonal components (αβ is yz, xz, or xy).

We will use a Python script will read the pressure tensor from binlog and compute the viscosity from the integral of the off-diagonal pressure tensor autocorrelation function. The resulting viscosity value is interpreted as the limiting value A of a double exponential fit defined as:

Practically, to extract the viscosity, download the viscosity extractor script, modify the path to your ams.rkf file (corresponding to your production run). You can then run the script with the following command (read more about Python scripting):

$AMSBIN/amspython viscosity_extractor.py /path/to/ams.rkf

Based on the 2 ns production simulation we obtain a limiting value of \(\eta\) = 24 mPa s.

We obtain a good agreement with the experimental viscosity of the mixture equal to 21.0 mPa s (at 333 K, see Borodin 2009 and Tokuda et al.).

Note that the viscosity may take a long time to converge. You should always ensure that the value is converged!

#!/usr/bin/env amspython
"""
Usage: 
$AMSBIN/amspython viscosity_extractor.py /path/to/ams.rkf

Feel free to change the max_dt_fs to change the maximum correlation time!

The script will produce a file green_kubo_viscosity.txt and viscosity.pdf in the same directory as ams.rkf.
"""

from scm.plams import *
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import os
import sys


def main():
    if len(sys.argv) > 1:
        path_to_rkf = sys.argv[1]
    else:
        path_to_rkf = "ams.rkf"

    assert os.path.exists(path_to_rkf), f"{path_to_rkf} does not seem to exist. It should be an ams.rkf file."
    max_dt_fs = 200000
    outfile = os.path.join(os.path.dirname(os.path.abspath(path_to_rkf)), "green_kubo_viscosity.txt")
    compute_pressure_tensor_acf(path_to_rkf, outfile=outfile, max_dt_fs=max_dt_fs)

    assert os.path.exists(outfile), f"{outfile} does not seem to exist!"
    imagefile = os.path.join(os.path.dirname(outfile), "viscosity.pdf")
    viscosity = do_fit(outfile, imagefile=imagefile)

    print(f"Final viscosity: {viscosity:.3f} mPa s")


# The double exponential function
def viscosity_double_exponential(x, A, lam, tau1, tau2):
    return A * (lam * (1 - np.exp(-x / tau1)) + (1 - lam) * (1 - np.exp(-x / tau2)))


# The main function
def compute_pressure_tensor_acf(path_to_rkf, outfile, max_dt_fs=100000, discard_first_fs=0):

    # Change this to the path of your ams.rkf file
    job = AMSJob.load_external(path_to_rkf)

    print(f"Computing viscosity for {path_to_rkf}")
    print("Note: this may take several hours if the trajectory is long!")

    # Here we get the off-diagonal pressure tensor from binlog
    xy = job.results.get_history_property("PressureTensor_xy", history_section="BinLog")
    xz = job.results.get_history_property("PressureTensor_xz", history_section="BinLog")
    yz = job.results.get_history_property("PressureTensor_yz", history_section="BinLog")

    minN = min(len(xy), len(xz), len(yz))
    pressuretensor = np.stack((np.zeros(minN), np.zeros(minN), np.zeros(minN), yz[:minN], xz[:minN], xy[:minN]), axis=1)

    # Volume, temperature, time
    volume = job.results.get_main_molecule().unit_cell_volume()
    temperature = np.mean(job.results.get_history_property("Temperature", history_section="MDHistory"))
    times = job.results.get_history_property("Time", history_section="BinLog")
    timestep = times[1] - times[0]

    # We skip first 10000 frames
    discard_first_frames = int(discard_first_fs / timestep)
    pressuretensor = pressuretensor[discard_first_frames:, :]

    # Time interval to compute the ACF
    # We found that max_dt_fs = 100000 fs is a good value
    max_dt = int(max_dt_fs / timestep)

    # Some info to print
    print(f"Averaged temperature: {temperature}")
    print(f"Timestep: {timestep}")
    print(f"Total simulation time: {times[-1]} fs")
    print(f"Total #frames: {minN}")
    print(f"Discarding first {discard_first_fs:.2f} fs as equilibration period")
    print(f"Discarding first {discard_first_frames} frames as equilibration period")
    print(f"Maximum correlation time: {max_dt_fs:.2f} fs")
    print(f"Maximum correlation time: {max_dt} frames")
    print(f"Unit cell volume: {volume:.3f} ang^3")

    print(f"Calculating the pressure tensor autocorrelation function (this may take a while....)")

    # Compute the viscosity based on GK relation
    x, y = AMSResults._get_green_kubo_viscosity(
        pressuretensor=pressuretensor, max_dt=max_dt, time_step=timestep, volume=volume, temperature=temperature
    )

    # We save the data
    print(f"Writing cumulative integral of pressure tensor autocorrelation function to {outfile}")
    A = np.stack((x, y), axis=1)
    np.savetxt(outfile, A, header="Time(fs) Viscosity(mPa*s)")


def do_fit(fname, imagefile):
    A = np.loadtxt(fname)
    x = A[:, 0]
    y = A[:, 1]

    # The fit
    f, p0 = viscosity_double_exponential, (20.0, 0.5, 1000, 10000)
    popt, _ = curve_fit(f, x, y, p0=p0)
    prediction = f(x, *popt)

    # We plot the results
    plt.title(f"Limiting value: {popt[0]:.2f} mPa s")
    plt.plot(x, y, color="tab:blue")
    plt.plot(x, prediction, color="tab:red")
    plt.xlabel("Correlation Time (fs)")
    plt.ylabel("η (mPa s)")
    plt.savefig(imagefile)
    plt.show()

    return popt[0]


if __name__ == "__main__":
    main()