Structure generation

The construction of force constant potentials (FCPs) crucially depends on the availability of suitable training and validation data. Once configurations are available the atomic forces can be readily computed, typically using methods such as density functional theory (DFT). The latter is demonstrated in the tutorial on reference calculations.

It is, however, essential that the configurations in the training and validation sets represent a distribution of displacements (and ultimately forces) that resemble configurations that are observed under dynamic conditions. While it is of course possible to run (ab-initio) molecular dynamics (MD) simulations to generate such configurations, it is computationally desirable to limit the number of such calculations.

Standard rattle (rattle)

It is straightforward to generate “rattled” configurations by randomly drawing displacements from a normal distribution. This is achieved e.g., by the rattle() function of the ASE Atoms class. In the present example it is shown that while this approach yields displacement distributions that resemble MD simulations the distribution of forces, however, exhibits a long tail with some atoms experiencing extremely large forces. This behavior can be traced to some interatomic distances becoming very small whence the atoms sample the very steep repulsive branch of the interaction potential. This is for example apparent from the pair distribution function, which is shown in the middle panel of the figure below.

Monte Carlo modified rattle procedure (MC-rattle)

To circumvent this problem hiphive provides a modified rattle procedure, which combines randomly drawing displacements with a Monte Carlo (MC) trial step that penalizes displacements that lead to very small interatomic distances. The resulting displacement and force distributions mimic the results from MD simulations and thus provide sensible input configurations for reference calculations without the need to carry out MD simulations using the reference method. This method does, however, require some material dependent fine tuning.

Superposition of normal modes (phonon-rattle)

Physically sound displacement patterns can in principle be constructed from a knowledge of the harmonic force constants. Superposing the respective normal modes with randomized phase factors allows one to readily generate structures representative of different temperatures without introducing additional parameters. This approach, referred to here as phonon-rattle, will also automatically account for differences in the thermal displacements between different species and crystallographic sites.

While this approach requires a set of harmonic force constants, a rough approximation is usually sufficient as a starting point as it can be iteratively improved as needed. One would thus start with a small number of rattled structures to construct a first approximation of the harmonic force constants. Using the latter one can then generate new structures by superposition of normal modes and then use these structures to construct a new model. This process can be continued in iterative fashion as needed.

The phonon-rattle approach yields distributions in rather close agreement with MD simulations and does not require tuning of hyperparameters (as in the case of MC-rattle).

Distribution of displacements, distances, and forces for structures obtained by the standard rattle and Monte Carlo rattle and phonon-rattle procedures as well as MD simulations using the reference potential.

Source code

Rattled structures are generated in
examples/advanced_topics/structure_generation/1_generate_rattle_structures.py

"""
Generate displaced structures using
* Standard rattle, gaussian displacements
* Monte Carlo rattle, gaussian displacements w penatly for short atomic dists
* Phonon rattle, construct a rough guess of fc2 and populate phonon modes with
  thermal energy corresponding to a temperature

The hyper parameters for the different methods are chosen such that the
magnitude of the displacements will be roughly the same

This script may take a few minutes to run.
"""

from ase.build import bulk
from ase.io import write
from ase.calculators.emt import EMT
from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from trainstation import Optimizer
from hiphive.utilities import prepare_structures
from hiphive.structure_generation import (generate_rattled_structures,
                                          generate_mc_rattled_structures,
                                          generate_phonon_rattled_structures)


# parameters
alat = 3.52         # lattice parameter
size = 6            # supercell size
n_structures = 25      # number of configurations to generate

# rattle parameters
T = 800
rattle_std = 0.12
min_distance = 0.67 * alat

# reference structure
prim = bulk('Ni', a=alat)
calc = EMT()

supercell = prim.repeat(size)
reference_positions = supercell.get_positions()
write('reference_structure.xyz', supercell)

# standard rattle
print('standard rattle')
structures_rattle = generate_rattled_structures(
    supercell, n_structures, rattle_std)
write('structures_rattle.extxyz', structures_rattle)

# Monte Carlo rattle
print('Monte Carlo rattle')
structures_mc_rattle = generate_mc_rattled_structures(
    supercell, n_structures, 0.25*rattle_std, min_distance, n_iter=20)
write('structures_mc_rattle.extxyz', structures_mc_rattle)


# Phonon rattle

# initial model
cs = ClusterSpace(supercell, [5.0])
sc = StructureContainer(cs)
rattled_structures = generate_rattled_structures(supercell, 1, 0.05)
rattled_structures = prepare_structures(rattled_structures, supercell, calc)

for structure in rattled_structures:
    sc.add_structure(structure)
opt = Optimizer(sc.get_fit_data(), train_size=1.0)
opt.train()
fcp = ForceConstantPotential(cs, opt.parameters)

# generate phonon rattled structures
fc2 = fcp.get_force_constants(supercell).get_fc_array(order=2, format='ase')
structures_phonon_rattle = generate_phonon_rattled_structures(
    supercell, fc2, n_structures, T)
write('structures_phonon_rattle_T{}.extxyz'.format(T), structures_phonon_rattle)  # noqa

Structures from molecular dynamics (MD) simulations are generated in
examples/advanced_topics/structure_generation/2_generate_md_structures.py

"""
Run molecular dynamics simulations and dump snapshots
"""
from ase.calculators.emt import EMT
from ase.build import bulk

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


# parameters
cell_size = 6  # system size
number_of_MD_steps = 500
time_step = 5  # in fs
dump_interval = 50
temperature = 800
traj_file = 'md_trajectory_T{}.traj'

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

# setup molecular dynamics simulations
atoms.set_calculator(calc)
dyn = Langevin(atoms, time_step * units.fs, temperature * units.kB, 0.02)
traj_writer = Trajectory(traj_file.format(temperature), 'w', atoms)
MaxwellBoltzmannDistribution(atoms, temperature * units.kB)

# run equilibration
dyn.run(number_of_MD_steps)

# run sample
dyn.attach(traj_writer.write, interval=dump_interval)
dyn.run(number_of_MD_steps)

The force distributions are analyzed in
examples/advanced_topics/structure_generation/3_analyze_force_distribution.py

"""
Generate and plot distributions of displacements, distances, and forces.
Forces are calculated using the EMT calculator.
"""
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
from ase.io import read
from ase.calculators.emt import EMT

bins = {'displacement': np.linspace(0.0, 0.7, 80),
        'distance': np.linspace(1.0, 4.5, 150),
        'force': np.linspace(0.0, 8.0, 50)}


def get_histogram_data(data, bins=100):
    counts, bins = np.histogram(data, bins=bins, density=True)
    bin_centers = [(bins[i+1]+bins[i])/2.0 for i in range(len(bins)-1)]
    return bin_centers, counts


def get_distributions(structure_list, ref_pos, calc):
    """ Gets distributions of interatomic distances and displacements.

    Parameters
    ----------
    structure_list : list(ase.Atoms)
        list of structures used for computing distributions
    ref_pos : numpy.ndarray
        reference positions used for computing the displacements (`Nx3` array)
    calc : ASE calculator object
        `calculator
        <https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html>`_
        used for computing forces
    """
    distances, displacements, forces = [], [], []
    for atoms in structure_list:
        distances.extend(atoms.get_all_distances(mic=True).flatten())
        displacements.extend(np.linalg.norm(atoms.positions-ref_pos, axis=1))
        atoms.set_calculator(calc)
        forces.extend(np.linalg.norm(atoms.get_forces(), axis=1))
    distributions = {}
    distributions['distance'] = get_histogram_data(distances, bins['distance'])
    distributions['displacement'] = get_histogram_data(
        displacements, bins['displacement'])
    distributions['force'] = get_histogram_data(forces, bins['force'])
    return distributions


# read atoms
T = 800
reference_structure = read('reference_structure.xyz')
ref_pos = reference_structure.get_positions()

structures_rattle = read('structures_rattle.extxyz@:')
structures_mc = read('structures_mc_rattle.extxyz@:')
structures_phonon = read('structures_phonon_rattle_T{}.extxyz@:'.format(T))
structures_md = read('md_trajectory_T{}.traj@:'.format(T))

calc = EMT()

# generate distributions
distributions_rattle = get_distributions(structures_rattle, ref_pos, calc)
distributions_mc = get_distributions(structures_mc, ref_pos, calc)
distributions_phonon = get_distributions(structures_phonon, ref_pos, calc)
distributions_md = get_distributions(structures_md, ref_pos, calc)

# plot
fs = 14
lw = 2.0
fig = plt.figure(figsize=(15, 5))
ax1 = fig.add_subplot(1, 3, 1)
ax2 = fig.add_subplot(1, 3, 2)
ax3 = fig.add_subplot(1, 3, 3)

units = OrderedDict(displacement='A', distance='A', force='eV/A')
for ax, key in zip([ax1, ax2, ax3], units.keys()):
    ax.plot(*distributions_rattle[key], lw=lw, label='Rattle')
    ax.plot(*distributions_mc[key], lw=lw, label='Monte Carlo rattle')
    ax.plot(*distributions_phonon[key], lw=lw, label='Phonon rattle {}K'
            .format(T))
    ax.plot(*distributions_md[key], lw=lw, label='MD {}K'.format(T))

    ax.set_xlabel('{} ({})'.format(key.title(), units[key]), fontsize=fs)
    ax.set_xlim([np.min(bins[key]), np.max(bins[key])])
    ax.set_ylim(bottom=0.0)
    ax.tick_params(labelsize=fs)
    ax.legend(fontsize=fs)

ax1.set_ylabel('Distribution', fontsize=fs)
plt.tight_layout()
plt.savefig('structure_generation_distributions.svg')