Thermal conductivity via non-equilibrium molecular dynamics

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 .

/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-thumbnail.png

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.

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
2. Change to ForceFieldPanel
3. Choose Task → Molecular Dynamics
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-1.png

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

1. Change to the panel Model → Regions
2. Switch on the periodic view, PeriodicViewTool
3. Select the molecule of the source region (hold down SHIFT key for multiple selections)
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-2A.png
4. Click on AddButton to generate a region
5. Enter name | source for the region
6. Repeat steps 3 → 5 with atoms of the sink region
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-2.png

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

1. Change to the panel Model → Thermostat
2. Click on AddButton to generate a 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
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-3.png

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 1000000
3. Set Time step to 1 fs
4. Set Initial temperature to 80 K
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-4.png

To start the calculation, do

File → Save and File → Run

On a recently modern 12-core desktop computer the calculation will take around 1.5 hours.

Results

Note

The calculation of the thermal conductivity requires the use of spreadsheet software such as Open Office Calc , or similar.

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

1. From SCMMenu select 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 potential energy Graph → Delete Curve
2. Add the energy of the first thermostat MDProperties → NHCTstat1Energy1
3. Add the energy of the second thermostat MDProperties → NHCTstat1Energy2
4. Add the time MDPropertiers → Time
5. Place the time on the X axes Graph → Curve on X Axes
/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-5.png

The energy of the sink thermostat is the one ascending over time, while the source is descending. To export these results into ASCII format from the GUI:

1. Export the graph via: Graph → Save as XY

Plot the sink energy over the simulation time, as opposed to the step(!) and calculate the heat flux (Q) via

\[Q = \frac{\left( \frac{dE}{dt} \right) /2 } {S} = 1174390731.54 Watt\]

where dE/dt is the slope of the graph sink energy [Hartree, a.u.] vs. time [fs]

/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-6.png

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

/scm-uploads/doc.2022/Tutorials/_images/T-NEMD-7.png

Tip

The lattice vectors can be found under Model → Lattice

The factor 2 stems from the periodic boundary conditions because the heat can flow in two directions. With the heat flux Q and the length of the conduction zone L, the thermal conductivity (k) can be calculated from Fourier’s law:

\[k_{\textrm{80K}} = \frac{Q*0.00436}{S} \frac{L}{\Delta T} = \textrm{ 0.176 W/mK (Exp.: 0.195 W/mK})\]

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.