Thermal Expansion

Calculate the thermal expansion for a Morse Pair potential using the LAMMPS molecular dynamics simulation code. In the following three methods to calculate the thermal expansion are introduced and compared for a Morse Pair Potential for Aluminium.

As a first step the potential is defined for the LAMMPS molecular dynamics simulation code by specifying the pair_style and pair_coeff commands for the Morse Pair Potential as well as the Aluminium bulk structure:

from ase.build import bulk
import pandas

potential_dataframe = pandas.DataFrame(
    {
        "Config": [
            ["pair_style morse/smooth/linear 9.0", "pair_coeff * * 0.5 1.8 2.95"]
        ],
        "Filename": [[]],
        "Model": ["Morse"],
        "Name": ["Morse"],
        "Species": [["Al"]],
    }
)

structure = bulk("Al", cubic=True)

The pandas.DataFrame based format to specify interatomic potentials is the same pylammpsmpi uses to interface with the NIST database for interatomic potentials. In comparison to just providing the pair_style and pair_coeff commands, this extended format enables referencing specific files for the interatomic potentials "Filename": [[]], as well as the atomic species "Species": [["Al"]], to enable consistency checks if the interatomic potential implements all the interactions to simulate a given atomic structure.

Finally, the last step of the preparation before starting the actual calculation is optimizing the interatomic structure. While for the Morse potential used in this example this is not necessary, it is essential for extending this example to other interactomic potentials. For the structure optimization the optimize_positions_and_volume() function is imported and applied on the ase.atoms.Atoms bulk structure for Aluminium:

from atomistics.workflows import optimize_positions_and_volume

task_dict = optimize_positions_and_volume(structure=structure)
task_dict

{'optimize_positions_and_volume': Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])}

It returns a task_dict with a single task, the optimization of the positions and the volume of the Aluminium structure. This task is executed with the LAMMPS molecular dynamics simulation code using the evaluate_with_lammpslib() function:

from atomistics.calculators import evaluate_with_lammpslib

result_dict = evaluate_with_lammpslib(
    task_dict=task_dict,
    potential_dataframe=potential_dataframe,
)
structure_opt = result_dict["structure_with_optimized_positions_and_volume"]
structure_opt

/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/calculators/__init__.py:63: UserWarning: calc_static_with_qe(), evaluate_with_qe() and optimize_positions_and_volume_with_qe() are not available as the import of the module named 'pwtools' failed.
  raise_warning(module_list=quantum_espresso_function, import_error=e)

Atoms(symbols='Al4', pbc=True, cell=[4.047310585424964, 4.047310585424964, 4.047310585424964])

The result_dict just contains a single element, the ase.atoms.Atoms structure object with optimized positions and volume. After this step the preparation is completed and the three different approximations can be compared in the following.

Equation of State

The first approximation to calculate the thermal expansion is based on the Equation of State derived by Moruzzi, V. L. et al.. So in analogy to the previous example of calculating the elastic properties from the Equation of State, the EnergyVolumeCurveWorkflow is initialized with the default parameters:

from atomistics.workflows import (
    analyse_results_for_energy_volume_curve, 
    get_tasks_for_energy_volume_curve,
    get_thermal_properties_for_energy_volume_curve,
)

task_dict = get_tasks_for_energy_volume_curve(
    structure=structure_opt.copy(),
    num_points=11,
    vol_range=0.05,
    axes=["x", "y", "z"],
)
task_dict

{'calc_energy': {0.95: Atoms(symbols='Al4', pbc=True, cell=[3.9786988461213992, 3.9786988461213992, 3.9786988461213992]),
  0.96: Atoms(symbols='Al4', pbc=True, cell=[3.992610493736228, 3.992610493736228, 3.992610493736228]),
  0.97: Atoms(symbols='Al4', pbc=True, cell=[4.00642586504517, 4.00642586504517, 4.00642586504517]),
  0.98: Atoms(symbols='Al4', pbc=True, cell=[4.020146608667117, 4.020146608667117, 4.020146608667117]),
  0.99: Atoms(symbols='Al4', pbc=True, cell=[4.033774328510742, 4.033774328510742, 4.033774328510742]),
  1.0: Atoms(symbols='Al4', pbc=True, cell=[4.047310585424964, 4.047310585424964, 4.047310585424964]),
  1.01: Atoms(symbols='Al4', pbc=True, cell=[4.060756898772644, 4.060756898772644, 4.060756898772644]),
  1.02: Atoms(symbols='Al4', pbc=True, cell=[4.074114747931804, 4.074114747931804, 4.074114747931804]),
  1.03: Atoms(symbols='Al4', pbc=True, cell=[4.087385573728375, 4.087385573728375, 4.087385573728375]),
  1.04: Atoms(symbols='Al4', pbc=True, cell=[4.100570779804249, 4.100570779804249, 4.100570779804249]),
  1.05: Atoms(symbols='Al4', pbc=True, cell=[4.113671733924125, 4.113671733924125, 4.113671733924125])}}

After the initialization the generate_structures() function is called to generate the atomistic structures which are then in the second step evaluated with the LAMMPS molecular dynamics simulation code to derive the equilibrium properties:

result_dict = evaluate_with_lammpslib(
    task_dict=task_dict, potential_dataframe=potential_dataframe
)
result_dict

{'energy': {0.95: -14.609207927145922,
  0.96: -14.656740101454448,
  0.97: -14.692359030099395,
  0.98: -14.71688372487553,
  0.99: -14.731079276327012,
  1.0: -14.735659820057943,
  1.01: -14.731295089579728,
  1.02: -14.718611862249293,
  1.03: -14.69819671584233,
  1.04: -14.670598736769113,
  1.05: -14.636332030744796}}

While in the previous example the fit of the energy volume curve was used directly, here the output of the fit, in particular the derived equilibrium properties are the input for the Debye model as defined by Moruzzi, V. L. et al.:

structure_opt.get_volume()

66.29787349319821

fit_dict = analyse_results_for_energy_volume_curve(
    output_dict=result_dict,
    task_dict=task_dict,
)
fit_dict

{'b_prime_eq': 6.218602473906521,
 'bulkmodul_eq': 218.25686441091622,
 'volume_eq': 66.29753112311258,
 'energy_eq': -14.735788882589588,
 'fit_dict': {'fit_type': 'polynomial',
  'least_square_error': 1.2103551168057661e-08,
  'poly_fit': array([-3.72876215e-04,  8.44360952e-02, -6.27903075e+00,  1.39077950e+02]),
  'fit_order': 3},
 'energy': [-14.609207927145922,
  -14.656740101454448,
  -14.692359030099395,
  -14.71688372487553,
  -14.731079276327012,
  -14.735659820057943,
  -14.731295089579728,
  -14.718611862249293,
  -14.69819671584233,
  -14.670598736769113,
  -14.636332030744796],
 'volume': [62.98297981853827,
  63.645958553470244,
  64.30893728840229,
  64.97191602333424,
  65.63489475826624,
  66.29787349319821,
  66.96085222813018,
  67.62383096306218,
  68.28680969799419,
  68.94978843292616,
  69.61276716785807]}

import numpy as np

thermal_properties_dict = get_thermal_properties_for_energy_volume_curve(
    fit_dict=fit_dict,
    masses=structure.get_masses(),
    temperatures=np.arange(1, 1500, 50),
    output_keys=["temperatures", "volumes"],
)
temperatures_ev, volume_ev = (
    thermal_properties_dict["temperatures"],
    thermal_properties_dict["volumes"],
)

The output of the Debye model provides the change of the temperature specific optimal volume volume_ev which can be plotted over the temperature temperatures_ev to determine the thermal expansion.

Quasi-Harmonic Approximation

While the Moruzzi, V. L. et al. approach based on the Einstein crystal is limited to a single frequency, the quasi-harmonic model includes the volume dependent free energy. Inside the atomistics package the harmonic and quasi-harmonic model are implemented based on an interface to the Phonopy framework. Still the user interface is still structured in the same three steps of (1) generating structures, (2) evaluating these structures and (3) fitting the corresponding model. Starting with the initialization of the QuasiHarmonicWorkflow which combines the PhonopyWorkflow with the EnergyVolumeCurveWorkflow:

from atomistics.workflows import (
    get_tasks_for_quasi_harmonic_approximation,
    analyse_results_for_quasi_harmonic_approximation,
    get_thermal_properties_for_quasi_harmonic_approximation,
)

task_dict, qh_dict = get_tasks_for_quasi_harmonic_approximation(
    structure=structure_opt.copy(),
    num_points=11,
    vol_range=0.10,
    # fit_type='birchmurnaghan',
    interaction_range=10,
    displacement=0.01,
)
task_dict

{'calc_energy': {0.9: Atoms(symbols='Al108', pbc=True, cell=[11.7229062192894, 11.7229062192894, 11.7229062192894]),
  0.92: Atoms(symbols='Al108', pbc=True, cell=[11.80910715486485, 11.80910715486485, 11.80910715486485]),
  0.94: Atoms(symbols='Al108', pbc=True, cell=[11.894067681419225, 11.894067681419225, 11.894067681419225]),
  0.96: Atoms(symbols='Al108', pbc=True, cell=[11.977831481208684, 11.977831481208684, 11.977831481208684]),
  0.98: Atoms(symbols='Al108', pbc=True, cell=[12.060439826001351, 12.060439826001351, 12.060439826001351]),
  1.0: Atoms(symbols='Al108', pbc=True, cell=[12.141931756274893, 12.141931756274893, 12.141931756274893]),
  1.02: Atoms(symbols='Al108', pbc=True, cell=[12.222344243795412, 12.222344243795412, 12.222344243795412]),
  1.04: Atoms(symbols='Al108', pbc=True, cell=[12.301712339412747, 12.301712339412747, 12.301712339412747]),
  1.06: Atoms(symbols='Al108', pbc=True, cell=[12.38006930767338, 12.38006930767338, 12.38006930767338]),
  1.08: Atoms(symbols='Al108', pbc=True, cell=[12.457446749652004, 12.457446749652004, 12.457446749652004]),
  1.1: Atoms(symbols='Al108', pbc=True, cell=[12.533874715230777, 12.533874715230777, 12.533874715230777])},
 'calc_forces': {(0.9,
   0): Atoms(symbols='Al108', pbc=True, cell=[11.7229062192894, 11.7229062192894, 11.7229062192894]),
  (0.92,
   0): Atoms(symbols='Al108', pbc=True, cell=[11.80910715486485, 11.80910715486485, 11.80910715486485]),
  (0.94,
   0): Atoms(symbols='Al108', pbc=True, cell=[11.894067681419225, 11.894067681419225, 11.894067681419225]),
  (0.96,
   0): Atoms(symbols='Al108', pbc=True, cell=[11.977831481208684, 11.977831481208684, 11.977831481208684]),
  (0.98,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.060439826001351, 12.060439826001351, 12.060439826001351]),
  (1.0,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.141931756274893, 12.141931756274893, 12.141931756274893]),
  (1.02,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.222344243795412, 12.222344243795412, 12.222344243795412]),
  (1.04,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.301712339412747, 12.301712339412747, 12.301712339412747]),
  (1.06,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.38006930767338, 12.38006930767338, 12.38006930767338]),
  (1.08,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.457446749652004, 12.457446749652004, 12.457446749652004]),
  (1.1,
   0): Atoms(symbols='Al108', pbc=True, cell=[12.533874715230777, 12.533874715230777, 12.533874715230777])}}

In contrast to the previous workflows which only used the calc_energy function of the simulation codes the PhonopyWorkflow and correspondingly also the QuasiHarmonicWorkflow require the calculation of the forces calc_forces in addition to the calculation of the energy. Still the general steps of the workflow remain the same:

result_dict = evaluate_with_lammpslib(
    task_dict=task_dict,
    potential_dataframe=potential_dataframe,
)

The structure_dict is evaluated with the LAMMPS molecular dynamics simulation code to calculate the corresponding energies and forces. The output is not plotted here as the forces for the 108 atom cells result in 3x108 outputs per cell. Still the structure of the result_dict again follows the labels of the structure_dict as explained before. Finally, in the third step the individual free energy curves at the different temperatures are fitted to determine the equilibrium volume at the given temperature using the analyse_structures() and get_thermal_properties() functions:

eng_internal_dict, phonopy_collect_dict = analyse_results_for_quasi_harmonic_approximation(
    qh_dict=qh_dict,
    output_dict=result_dict,
    dos_mesh=20,
    number_of_snapshots=None,
)
thermal_properties_dict_qm = get_thermal_properties_for_quasi_harmonic_approximation(
    eng_internal_dict=eng_internal_dict,
    task_dict=task_dict,
    qh_dict=qh_dict,
    fit_type="polynomial",
    fit_order=3,
    temperatures=np.arange(1, 1500, 50),
    output_keys=["free_energy", "temperatures", "volumes"],
    quantum_mechanical=True,
)
temperatures_qh_qm, volume_qh_qm = (
    thermal_properties_dict_qm["temperatures"],
    thermal_properties_dict_qm["volumes"],
)

The optimal volume at the different temperatures is stored in the volume_qh_qm in analogy to the previous section. Here the extension _qm indicates that the quantum-mechanical harmonic oszillator is used.

thermal_properties_dict_cl = get_thermal_properties_for_quasi_harmonic_approximation(
    eng_internal_dict=eng_internal_dict,
    task_dict=task_dict,
    qh_dict=qh_dict,
    fit_type="polynomial",
    fit_order=3,
    temperatures=np.arange(1, 1500, 50),
    output_keys=["free_energy", "temperatures", "volumes"],
    quantum_mechanical=False,
)
temperatures_qh_cl, volume_qh_cl = (
    thermal_properties_dict_cl["temperatures"],
    thermal_properties_dict_cl["volumes"],
)

For the classical harmonic oszillator the resulting volumes are stored as volume_qh_cl.

Molecular Dynamics

Finally, the third and most commonly used method to determine the volume expansion is using a molecular dynamics calculation. While the atomistics package already includes a LangevinWorkflow at this point we use the Nose-Hoover thermostat implemented in LAMMPS directly via the LAMMPS calculator interface.

from atomistics.calculators import calc_molecular_dynamics_thermal_expansion_with_lammpslib

structure_md = structure_opt.copy().repeat(11)
result_dict = calc_molecular_dynamics_thermal_expansion_with_lammpslib(
    structure=structure_md,  # atomistic structure
    potential_dataframe=potential_dataframe,  # interatomic potential defined as pandas.DataFrame
    Tstart=15,  # temperature to for initial velocity distribution
    Tstop=1500,  # final temperature
    Tstep=5,  # temperature step
    Tdamp=0.1,  # temperature damping of the thermostat
    run=100,  # number of MD steps for each temperature
    thermo=100,  # print out from the thermostat
    timestep=0.001,  # time step for molecular dynamics
    Pstart=0.0,  # initial pressure
    Pstop=0.0,  # final pressure
    Pdamp=1.0,  # barostat damping
    seed=4928459,  # random seed
    dist="gaussian",  # Gaussian velocity distribution
)
temperature_md_lst, volume_md_lst = result_dict["temperatures"], result_dict["volumes"]

100%|██████████| 298/298 [06:26<00:00,  1.30s/it]

The calc_molecular_dynamics_thermal_expansion_with_lammpslib() function defines a loop over a vector of temperatures in 5K steps. For each step 100 molecular dynamics steps are executed before the temperature is again increased by 5K. For ~280 steps with the Morse Pair Potential this takes approximately 5 minutes on a single core. These simulations can be further accelerated by adding the cores parameter. The increase in computational cost is on the one hand related to the large number of force and energy calls and on the other hand to the size of the atomistic structure, as these simulations are typically executed with >5000 atoms rather than the 4 or 108 atoms in the other approximations. The volume for the individual temperatures is stored in the volume_md_lst list.

Summary

To visually compare the thermal expansion predicted by the three different approximations, the matplotlib is used to plot the volume over the temperature:

import matplotlib.pyplot as plt

plt.plot(
    temperature_md_lst,
    np.array(volume_md_lst) / len(structure_md) * len(structure_opt),
    label="Molecular Dynamics",
    color="C2",
)
plt.plot(temperatures_qh_qm, volume_qh_qm, label="Quasi-Harmonic (qm)", color="C3")
plt.plot(temperatures_qh_cl, volume_qh_cl, label="Quasi-Harmonic (classic)", color="C0")
plt.plot(temperatures_ev, volume_ev, label="Moruzzi Model", color="C1")
plt.axhline(structure_opt.get_volume(), linestyle="--", color="red")
plt.legend()
plt.xlabel("Temperature (K)")
plt.ylabel("Volume ($\AA^3$)")

Text(0, 0.5, 'Volume ($\\AA^3$)')

Both the Moruzzi, V. L. et al. and the quantum mechanical version of the quasi-harmonic approach start at a larger equilibrium volume as they include the zero point vibrations, resulting in an over-prediction of the volume expansion with increasing temperature. The equilibrium volume is indicated by the dashed red line. Finally, the quasi-harmonic approach with the classical harmonic oszillator agrees very well with the thermal expansion calculated from molecular dynamics for this example of using the Morse Pair Potential.