Molecular dynamics simulations

The objective of this tutorial section is to demonstrate the usage of an FCP object in a molecular dynamics (MD) simulation. Such simulations can provide very rich information pertaining to the system including but not limited to free energies, the lattice thermal conductivity, or phonon lifetimes, fully accounting for anharmonic effects in the classical limit.

The integration of the equations of motion will be carried out using functionality provided by ASE while hiPhive is used to provide an interaction model in the form of an ASE calculator object.

Preparations

First a number of parameters are set that define

  • the size of the simulation cell (cell_size),

  • the number of MD steps (number_of_MD_steps),

  • the time step in fs (time_step),

  • the temperatures at which MD simulations will be carried out (temperatures),

  • the frequency at which configurations will be written to disk (dump_interval), and

  • the names of the output files (log_file and traj_file)

"""
Run molecular dynamics simulations using the fourth order
hiPhive force constant potential and the ASE MD module.

Runs in approximately 500 seconds on an Intel Core i5-4670K CPU.
"""

import os
import numpy as np
from ase.build import bulk
from hiphive import ForceConstantPotential
from hiphive.calculators import ForceConstantCalculator

from ase import units
from ase.io.trajectory import Trajectory
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.md import MDLogger


# parameters
cell_size = 6  # system size
number_of_MD_steps = 1000
time_step = 5  # in fs
dump_interval = 20
temperatures = [600, 1200]
log_file = 'md_runs/logs_T{}'
traj_file = 'md_runs/trajs_T{}.traj'
if not os.path.isdir(os.path.dirname(log_file)):
    os.mkdir(os.path.dirname(log_file))

Next we prepare a supercell.

atoms = bulk('Ni').repeat(cell_size)
reference_positions = atoms.get_positions()

Then we initialize an ASE calculator object using the FCP prepared previously. To this end, we construct the force constant matrices for the supercell to be used in the MD simulations using the ForceConstantPotential.get_force_constants function, which are subsequently used to initialize an ASE calculator via the ForceConstantCalculator class.

fcp = ForceConstantPotential.read('fcc-nickel.fcp')
fcs = fcp.get_force_constants(atoms)
calc = ForceConstantCalculator(fcs)

MD simulations

We are ready to carry out MD simulations using the functionality provided by ASE. The ForceConstantCalculator is attached to the supercell structure, an integrator that samples Langevin dynamics is initialized, and the output is prepared. Finally, the atomic velocities are initialized and the MD simulation is run

atoms.set_calculator(calc)
for temperature in temperatures:
    dyn = Langevin(atoms, time_step * units.fs, temperature * units.kB, 0.02)
    logger = MDLogger(dyn, atoms, log_file.format(temperature),
                      header=True, stress=False, peratom=True, mode='w')
    traj_writer = Trajectory(traj_file.format(temperature), 'w', atoms)
    dyn.attach(logger, interval=dump_interval)
    dyn.attach(traj_writer.write, interval=dump_interval)

    # run MD
    MaxwellBoltzmannDistribution(atoms, temperature * units.kB)
    dyn.run(number_of_MD_steps)

Once the MD simulations have concluded the mean-square displacements (MSDs) are computed.

for temperature in temperatures:
    traj_reader = Trajectory(traj_file.format(temperature), 'r')
    msd = []
    for atoms in [a for a in traj_reader][10:]:
        displacements = atoms.positions - reference_positions
        msd.append(np.mean(np.sum(displacements**2, axis=1)))
    print('T = {:4d}    MSD = {:.5f} A**2'.format(temperature, np.mean(msd)))

These values can now be compared with the results obtained in the harmonic approximation:

Molecular dynamics simulations based on FCP
-------------------------------------------
T =  600    msd = 0.02146 A**2
T = 1200    msd = 0.04700 A**2

Harmonic approximation based on FCP
-------------------------------------------
T =  600 K  MSD = 0.02550 A**2
T = 1200 K  MSD = 0.05049 A**2

This comparison indicates that the harmonic approximation systematically overestimates the MSDs. It must be noted though that the MSD analysis from phonopy includes quantum-mechanical effects, most notably zero-point motion, whereas the MD simulation is purely classical. More importantly with respect to the verification of the force constant model, a comprehensive study with tighter convergence parameters demonstrates that the fourth-order model constructed here closely reproduces the results from MD simulations that employ the original (EMT) potential directly, as shown in the figure below.

../_images/msds_nickel.svg

Mean-square displacement of FCC Ni as a function of temperature as obtained within the harmonic approximation as well as from molecular dynamics (MD) simulations based on the original potential (EMT) and a fourth-order Hamiltonian constructed from MC rattled configurations.

Source code

The complete source code is available in examples/tutorial/5_run_fourth_order_MD.py

"""
Run molecular dynamics simulations using the fourth order
hiPhive force constant potential and the ASE MD module.

Runs in approximately 500 seconds on an Intel Core i5-4670K CPU.
"""

import os
import numpy as np
from ase.build import bulk
from hiphive import ForceConstantPotential
from hiphive.calculators import ForceConstantCalculator

from ase import units
from ase.io.trajectory import Trajectory
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.langevin import Langevin
from ase.md import MDLogger


# parameters
cell_size = 6  # system size
number_of_MD_steps = 1000
time_step = 5  # in fs
dump_interval = 20
temperatures = [600, 1200]
log_file = 'md_runs/logs_T{}'
traj_file = 'md_runs/trajs_T{}.traj'
if not os.path.isdir(os.path.dirname(log_file)):
    os.mkdir(os.path.dirname(log_file))

# set up supercell
atoms = bulk('Ni').repeat(cell_size)
reference_positions = atoms.get_positions()

# get force constant calculator
fcp = ForceConstantPotential.read('fcc-nickel.fcp')
fcs = fcp.get_force_constants(atoms)
calc = ForceConstantCalculator(fcs)

# run molecular dynamics simulations
atoms.set_calculator(calc)
for temperature in temperatures:
    dyn = Langevin(atoms, time_step * units.fs, temperature * units.kB, 0.02)
    logger = MDLogger(dyn, atoms, log_file.format(temperature),
                      header=True, stress=False, peratom=True, mode='w')
    traj_writer = Trajectory(traj_file.format(temperature), 'w', atoms)
    dyn.attach(logger, interval=dump_interval)
    dyn.attach(traj_writer.write, interval=dump_interval)

    # run MD
    MaxwellBoltzmannDistribution(atoms, temperature * units.kB)
    dyn.run(number_of_MD_steps)

# compute mean-square displacement from MD trajectories
for temperature in temperatures:
    traj_reader = Trajectory(traj_file.format(temperature), 'r')
    msd = []
    for atoms in [a for a in traj_reader][10:]:
        displacements = atoms.positions - reference_positions
        msd.append(np.mean(np.sum(displacements**2, axis=1)))
    print('T = {:4d}    MSD = {:.5f} A**2'.format(temperature, np.mean(msd)))