# Structure generation¶

The construction of force constant potentials (FCPs) is crucially dependent on the availability of suitable training data. Once configurations are available the atomic forces can be readily computed, typically using methods such as density functional theory (DFT).

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.

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

To circumvent this problem hiPhive implements a modified rattle procedure, which combines randomly drawing displacements with a Monte Carlo 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.

## Source code¶

Rattled structures are generated in
tutorial/advanced/structure_generation/generate_rattle_structures.py
"""
Generate rattled structures using both the standard rattle approach and the
Monte Carlo rattle approach.

The parameters in this example are chosen to generate a similar
magnitude of the displacements using both methods.

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

from ase.build import bulk
from ase.io import write
from hiphive.structure_generation import (generate_rattled_structures,
generate_mc_rattled_structures)

# parameters
alat = 3.52         # lattice parameter
size = 6            # supercell size
n_configs = 25      # number of configurations to generate
rattle_std = 0.12
min_distance = 0.67 * alat

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

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_configs, rattle_std)
write('structures_rattle.extxyz', structures_rattle)

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

Structures from molecular dynamics (MD) simulations are generated in
tutorial/advanced/structure_generation/generate_md_structures.py
"""
Generate rattled structures using both the standard rattle approach and the
Monte Carlo rattle approach.

The parameters in this example are chosen to generate a similar
magnitude of the displacements using both methods.

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

from ase.build import bulk
from ase.io import write
from hiphive.structure_generation import (generate_rattled_structures,
generate_mc_rattled_structures)

# parameters
alat = 3.52         # lattice parameter
size = 6            # supercell size
n_configs = 25      # number of configurations to generate
rattle_std = 0.12
min_distance = 0.67 * alat

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

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_configs, rattle_std)
write('structures_rattle.extxyz', structures_rattle)

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

The force distributions are analyzed in
tutorial/advanced/structure_generation/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.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):
""" Get distributions of interatomic distances, and displacements.

Parameters
----------
structure_list : list of ASE Atoms objects
list of structures used for computing distributions
ref_pos : Numpy array (N, 3)
reference positions used for computing the displacements
calc : ASE Calculator object
calculator 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

temperature = 800
ref_pos = reference_structure.get_positions()

calc = EMT()

# generate distributions
rattle_distributions = get_distributions(rattle_structures, ref_pos, calc)
mc_rattle_distributions = get_distributions(
mc_rattle_structures, ref_pos, calc)
md_distributions = get_distributions(md_structures, ref_pos, calc)

# plot
fs = 14
lw = 2.0
fig = plt.figure(figsize=(12, 5))

units = OrderedDict(displacement='A', distance='A', force='eV/A')
for ax, key in zip([ax1, ax2, ax3], units.keys()):
ax.plot(*rattle_distributions[key], label='rattle', lw=lw)
ax.plot(*mc_rattle_distributions[key], label='Monte Carlo rattle', lw=lw)
ax.plot(*md_distributions[key], label='MD {}K'.format(temperature), lw=lw)

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