# Thermal conductivity from NEMD¶

In this tutorial the temperature-dependent thermal conductivity is calculated from a non-equilibrium molecular dynamics (NEMD ) trajectory. The AMS engine used is the UFF force field .

A heat flow is realized by running two local thermostats, heat source and heat sink, inside the same molecular structure. The heat source is set to be 10K above the target temperature, while the heat sink is kept 10K below the target temperature. From the energy collected at the heat sink over time and the cross-sectional area, the heat flux and thermal conductivity can be calculated. Both the molecular structure and the workflow were taken from the combined experimental and theoretical study E.M. Moscarello, B.L. Wooten, H. Sajid, L.D. Tichenor, J.P. Heremans, M.A. Addicoat, P.L. McGrier, ACS Appl. Nano Mater. 2022, 5, 10, 13787–13793.

The heat flux Q can be calculated as

where *S* is the cross-sectional area. The heat flux is related to the thermal conductivity via Fourier’s law

where *L* is the length of the conduction zone and \(\Delta T = T_{\textrm{source}} - T_{\textrm{sink}}\).

Combining these two equations gives the thermal conductivity *k* as

## Setup¶

Begin by downloading the input structure `BBO_COF1.xyz`

. The file contains the cartesian
coordinates of a benzobisoxazole (BBO)-linked COF for which the thermal conductivity at 80K will be calculated.

Note

The periodic structure has been relaxed with GFN1-xTB (including the lattice), as described in the paper .

**1.**Download the coordinates

`BBO_COF1.xyz`

and import them into AMSinput**3.**Choose

**Task → Molecular Dynamics**

Define a heat source and a heat sink by using regions .

**1.**Change to the panel

**Model → Regions**

**3.**Select the atoms of the source region (hold down SHIFT key for multiple selections)

**5.**Enter name |

`source`

for the region**6.**Repeat steps

**3 → 5**with atoms of the

`sink`

regionThe source and sink regions can now be connected to two independent thermostats

**1.**Change to the panel

**Model → MD… → Thermostat**

**3.**Choose

**Thermostat → NHC**

**4.**Set

**Termperature**to

`90`

**5.**Set

**Damping constant**to

`10`

**6.**Choose

**Atoms in region → source**

**7.**Repeat steps

**2 → 6**to generate an

**NHC**thermostat at

**70K**for

**Atoms in region → sink**

As the final step, define the settings for the molecular dynamics calculation

**1.**Change to the panel

**Model → MD**

**2.**Set

**Number of steps**to

`200000`

**3.**Set

**Time step**to

`1`

fs**4.**Set

**Sample frequency**to

`1000`

**5.**Set

**Initial temperature**to

`80`

K**6.**Set

**Checkpoint frequency**to

`200000`

. This prevents AMS from writing checkpoint files that needlessly take up disk space.To start the calculation, do

**File → Save**and

**File → Run**

On a recently modern 12-core desktop computer the calculation will take around 20 minutes.

## Results¶

Once the calculation has finished we can view the results in AMSmovie.

**2.**Make sure that the x-axis on the right hand side shows

**MD-Time (fs)**. If it instead shows

**Frame**, toggle

**Graph → Try to Use Time On X-Axes**and restart AMSmovie.

The relevant bits for the NEMD calculation are the energies of the two thermostats as well the time. To visualize these values:

**1.**Remove the curve of the total energy

**Graph → Delete Curve**

**2.**Add the energy of the second thermostat

**MDProperties → NHCTstat2Energy**

**3.**If the Energy is not shown in in eV, double-click on the

**Energy (eV)**y-axis and set the unit to eV.

After the simulation has equilibrated, this curve should increase linearly with time. In this tutorial, we ran a relatively short simulation of 200 ps but ideally it would be run for much longer.

Starting with AMS2023, the **slope of the curve** can easily be calculated
inside AMSmovie.

**1.**

**Graph → Analysis**

**2.**Switch to the

**Linear Regression**panel

**4.**In the

**Start**field, enter

`100000`

**5.**Click

**OK**

This will perform linear fit starting from 100,000 fs. Here, we consider the first 100,000 fs to be part of the equilibration period.

From the above figure, you see that the **slope = dE/dt = 4.41 × 10⁻⁵ eV/fs = 7.07 × 10⁻⁹ J/s**.

**S** is the cross-sectional area, i.e. the area from which heat enters the system, that can be calculated from the lattice vectors **a** and **c** of the system:

**1.**In the

**AMSinput**window, switch to the

**Model → Lattice**panel

**2.**Read off the area S = 281.97 Å².

**3.**The conduction length L is half the y-component of the

**b**lattice vector: L = 72.305 /2 = 36.15 Å.

Recall the thermal conductivity `k`

:

The factor 2 stems from the periodic boundary conditions because the heat can flow in two directions.

Inserting the above numbers (with ΔT = 90 - 70 = 20 K) yields **k = 0.23 W m⁻¹ K⁻¹** (experimental value: 0.195 W m⁻¹ K⁻¹)

Note

Due to the non-deterministic nature of the MD thermostats and the fitting procedure, the result of this tutorial can vary slightly. Still, the thermal conductivity should be reasonably close to the experimental value.

## Easier calculation with Python¶

You can also use for example the below Python script to automate the procedure of performing the linear fit and the unit conversions:

```
#!/usr/bin/env amspython
from scm.plams import *
import numpy as np
def main():
# Unit conversion factors
eV_to_J = Units.convert(1.0, "eV", "J")
fs_to_s = Units.convert(1.0, "fs", "s")
ang_to_m = Units.convert(1.0, "angstrom", "m")
# length L and cross-sectional area S
L = 36.15 * ang_to_m # m
S = 281.9 * ang_to_m**2 # m^2
# manually calculated slope and delta_T
dEdt = 4.41e-5 * eV_to_J / fs_to_s # J s^-1 = W
delta_T = 20 # K
# automatically calculated slope
# job = AMSJob.load_external("jobname.results")
# dEdt = get_slope(job, min_time_fs=100000)
# delta_T = get_delta_T(job)
# calculate thermal conductivity k
k = (dEdt * L) / (2 * S * delta_T)
print(f"Thermal conductivity: {k:.3f} W m^-1 K^-1")
def get_slope(job, min_time_fs=0):
"""
Calculates the (absolute value of the) slope of NHCTstat2Energy.
Discards the first min_time_fs from the trajectory.
Returns: float. The slope in W ( = J/s)
"""
from scipy.stats import linregress
time = job.results.get_history_property("Time", history_section="MDHistory")
time = np.array(time) # time in fs
nhctstatenergy = job.results.get_history_property("NHCTstat2Energy", history_section="MDHistory")
nhctstatenergy = np.array(nhctstatenergy) # energy in hartree
ind = time >= min_time_fs
time = time[ind]
nhctstatenergy = nhctstatenergy[ind]
result = linregress(time, nhctstatenergy)
# take absolute value, convert to watt (J/s)
slope = np.abs(result.slope) * Units.convert(1.0, "hartree", "J") / Units.convert(1.0, "fs", "s")
return slope
def get_delta_T(job):
"""
Assume there are two Thermostats and take the difference of their temperatures
Returns: float. The difference in temperature in Kelvin
"""
t = job.settings.input.ams.MolecularDynamics.Thermostat
assert isinstance(t, list) and len(t) == 2, f"get_delta_T can only be called if there are exactly two thermostats"
delta_T = abs(float(t[0].Temperature) - float(t[1].Temperature))
return delta_T
if __name__ == "__main__":
main()
```