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

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)