# 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).

## Source code¶

Rattled structures are generated in
examples/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 hiphive.fitting 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)
prepare_structures(rattled_structures, supercell, calc)

for structure in rattled_structures:
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/structure_generation/2_generate_md_structures.py
The force distributions are analyzed in
examples/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

T = 800
ref_pos = reference_structure.get_positions()

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')