Force constant potential construction
This section is intended to provide an overview of the basic procedure for constructing a force constant potential (FCP) using structures generated previously. Specifically, we will obtain a fourth-order FCP for FCC Ni to be utilized and analyzed in the subsequent sections.
The construction of an FCP comprises the following key steps
set up a cluster space using the
ClusterSpaceclassexpress a set of reference structures in the cluster space generated previously using the
StructureContainerclasstrain the parameters associated with each orbit in the cluster space using an optimizer class (e.g.,
Optimizer,EnsembleOptimizer, orCrossValidationEstimator)set up an FCP from the final parameters using the
ForceConstantPotentialclass
General preparations
After importing necessary functions and classes from ASE and hiPhive, we read the structures produced in the previous section.
"""
Construct a ForceConstantPotential from training data generated previously.
Runs in approximately 100 seconds on an Intel Core i5-4670K CPU.
"""
from ase.io import read
from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from hiphive.utilities import prepare_structures
from trainstation import Optimizer
# read structures containing displacements and forces
prim = read('prim.extxyz')
atoms_ideal = read('supercell_ideal.extxyz')
rattled_structures = read('supercells_rattled.extxyz', index=':')
Preparation of cluster space
In order to be able to build an FCP, it is first necessary to create a
ClusterSpace object based on a prototype
structure. Here, we simply employ the first structure from our reference data
set. Furthermore, one must specify the cutoffs up to the desired order. Here,
the cutoffs are set to 5, 4, and 4 Å for pairs, triplets, and quadruplets,
respectively.
cutoffs = [5.0, 4.0, 4.0]
cs = ClusterSpace(prim, cutoffs)
print(cs)
cs.print_orbits()
As with many other hiPhive objects, it is possible to print core
information in a tabular format by simply calling the print() function
with the instance of interest as input argument. For example print(cs)
results in the following output:
=================== Cluster Space ====================
Spacegroup : Fm-3m (225)
symprec : 1e-05
Length scale : 0.1
Cutoffs : {2: 5.0, 3: 4.0, 4: 4.0}
Cell :
[[ 0. 1.76 -1.76]
[ 1.76 1.76 0. ]
[ 1.76 0. -1.76]]
Basis : [[ 0. 0. 0.]]
Numbers : [28]
Total number of orbits : 20
Total number of parameters : 19
Total number of clusters : 18
------------------------------------------------------
order | num-orbits | num-params | num-clusters
------------------------------------------------------
2 | 1 | 2 | 3
3 | 1 | 7 | 12
4 | 1 | 10 | 3
======================================================
A special function (ClusterSpace.print_orbits) is available for printing a list of the
orbits associated with the cluster space:
===================================== List of Orbits =====================================
index | order | elements | radius | prototype | clusters | parameters
------------------------------------------------------------------------------------------
0 | 2 | Ni Ni | 0.0000 | (0, 0) | 1 | 1
1 | 2 | Ni Ni | 1.2445 | (0, 1) | 6 | 3
2 | 2 | Ni Ni | 2.4890 | (0, 2) | 6 | 3
3 | 2 | Ni Ni | 2.1556 | (0, 3) | 12 | 4
4 | 2 | Ni Ni | 1.7600 | (0, 12) | 3 | 2
5 | 3 | Ni Ni Ni | 1.1062 | (0, 0, 1) | 12 | 5
6 | 3 | Ni Ni Ni | 1.5644 | (0, 0, 12) | 6 | 3
7 | 3 | Ni Ni Ni | 1.4370 | (0, 1, 5) | 8 | 7
8 | 3 | Ni Ni Ni | 1.6279 | (0, 1, 13) | 12 | 7
9 | 4 | Ni Ni Ni Ni | 0.0000 | (0, 0, 0, 0) | 1 | 2
10 | 4 | Ni Ni Ni Ni | 0.9334 | (0, 0, 0, 1) | 12 | 9
11 | 4 | Ni Ni Ni Ni | 1.3200 | (0, 0, 0, 12) | 6 | 5
12 | 4 | Ni Ni Ni Ni | 1.2445 | (0, 0, 1, 1) | 6 | 9
13 | 4 | Ni Ni Ni Ni | 1.3621 | (0, 0, 1, 5) | 24 | 30
14 | 4 | Ni Ni Ni Ni | 1.4239 | (0, 0, 1, 13) | 12 | 17
15 | 4 | Ni Ni Ni Ni | 1.6044 | (0, 0, 1, 14) | 24 | 28
16 | 4 | Ni Ni Ni Ni | 1.7600 | (0, 0, 12, 12) | 3 | 6
17 | 4 | Ni Ni Ni Ni | 1.5242 | (0, 1, 5, 17) | 2 | 6
18 | 4 | Ni Ni Ni Ni | 1.6291 | (0, 1, 5, 39) | 12 | 24
19 | 4 | Ni Ni Ni Ni | 1.7600 | (0, 1, 13, 14) | 3 | 10
==========================================================================================
Compilation of structure container
Next one needs to compile a StructureContainer that will be used for training the parameters
associated with each orbit. The latter is initialized by providing a
ClusterSpace object, after which several
structures are added to the container.
structures = prepare_structures(rattled_structures, atoms_ideal)
sc = StructureContainer(cs)
for structure in structures:
sc.add_structure(structure)
print(sc)
We can print a concise summary of the StructureContainer via print(sc), which yields:
=============== Structure Container ================
Total number of structures : 5
----------------------------------------------------
index | num-atoms | avg-disp | avg-force | max-force
----------------------------------------------------
0 | 256 | 0.1294 | 1.6524 | 3.7407
1 | 256 | 0.1395 | 1.7285 | 4.6880
2 | 256 | 0.1296 | 1.5673 | 3.6984
3 | 256 | 0.1275 | 1.5965 | 5.0317
4 | 256 | 0.1243 | 1.6245 | 3.4472
====================================================
Training of parameters
The StructureContainer object
created in the previous section contains all the information required
for constructing an FCP. The next step is thus to train the parameters
using the target data. More precisely, the goal is to achieve the best
possible agreement with a set of training data, which represents a
subset of all the data stored in the StructureContainer. In practice, this is a two step process
that involves (1) the initiation of an Optimizer object with the list of target properties
produced by the StructureContainer.get_fit_data method as input argument,
and (2) calling the train
function of the optimizer object.
opt = Optimizer(sc.get_fit_data())
opt.train()
print(opt)
Note 1: By default the optimizer will employ a least-squares optimization
algorithm. Other algorithms can be easily selected via the fit_method
keyword of the optimizer object.
Note 2: More elaborate optimizers are available via the
EnsembleOptimizer and
CrossValidationEstimator
class.
A concise overview of state of the optimizer after training is obtained via
print(opt):
===================== Optimizer ======================
fit_method : least-squares
number_of_target-values : 3840
number_of_parameters : 119
rmse_train : 0.01469551
rmse_test : 0.0161488
train_size : 2880
test_size : 960
======================================================
Various data can also be accessed directly as attributes of the
Optimizer object. For example one can
access the root-mean-squared errors (RMSE) over the training and testing sets as
follows:
print('RMSE train : {:.4f}'.format(opt.rmse_train))
print('RMSE test : {:.4f}'.format(opt.rmse_test))
Set up force constant potential
Finally, we join the final parameters from the optimization step with the
orbits of the cluster space into a ForceConstantPotential that can be written to file for later
use.
fcp = ForceConstantPotential(cs, opt.parameters)
fcp.write('fcc-nickel.fcp')
print(fcp)
Source code
The complete source code is available in
examples/tutorial/2_construct_fcp.py
"""
Construct a ForceConstantPotential from training data generated previously.
Runs in approximately 100 seconds on an Intel Core i5-4670K CPU.
"""
from ase.io import read
from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
from hiphive.utilities import prepare_structures
from trainstation import Optimizer
# read structures containing displacements and forces
prim = read('prim.extxyz')
atoms_ideal = read('supercell_ideal.extxyz')
rattled_structures = read('supercells_rattled.extxyz', index=':')
# set up cluster space
cutoffs = [5.0, 4.0, 4.0]
cs = ClusterSpace(prim, cutoffs)
print(cs)
cs.print_orbits()
# ... and structure container
structures = prepare_structures(rattled_structures, atoms_ideal)
sc = StructureContainer(cs)
for structure in structures:
sc.add_structure(structure)
print(sc)
# train model
opt = Optimizer(sc.get_fit_data())
opt.train()
print(opt)
# construct force constant potential
fcp = ForceConstantPotential(cs, opt.parameters)
fcp.write('fcc-nickel.fcp')
print(fcp)