The force constant potential
Contents
Special requirements
- When you use this potential, currently you need to add -DUSE_FCP in the makefile. We may remove this requirement in the future. In this case, you cannot use any other potential.
- Neighbor list is not used for this potential, but you still need to specify the
MN
andr_c
parameters in xyz.in. You can choose any reasonable values for them and just remember that they will not used for calculating the forces and related quantities.
Potential form
In the force constant potential, the potential energy is calculated as a Taylor expansion in terms of the atomic displacements [math]\{u_i^{a}\}[/math] relative to a set of reference (equilibrium) positions as
$$
U = U_2 + U_3 + U_4 + U_5 + U_6 + \cdots;
$$
$$
U_2 =
\frac{1}{2!}
\sum_{ij}
\sum_{ab}
\Phi_{ij}^{ab}
u_{i}^{a} u_{j}^{b};
$$
$$
U_3 =
\frac{1}{3!}
\sum_{ijk}
\sum_{abc}
\Phi_{ijk}^{abc}
u_{i}^{a} u_{j}^{b} u_{k}^{c};
$$
$$
U_4 =
\frac{1}{4!}
\sum_{ijkl}
\sum_{abcd}
\Phi_{ijkl}^{abcd}
u_{i}^{a} u_{j}^{b} u_{k}^{c} u_{l}^{d};
$$
$$
U_5 =
\frac{1}{5!}
\sum_{ijklm}
\sum_{abcde}
\Phi_{ijklm}^{abcde}
u_{i}^{a} u_{j}^{b} u_{k}^{c} u_{l}^{d} u_{m}^{e};
$$
$$
U_6 =
\frac{1}{6!}
\sum_{ijklmn}
\sum_{abcdef}
\Phi_{ijklmn}^{abcdef}
u_{i}^{a} u_{j}^{b} u_{k}^{c} u_{l}^{d} u_{m}^{e} u_{n}^{f}.
$$
Here, [math]\Phi_{ij}^{ab}[/math], [math]\Phi_{ijk}^{abc}[/math], [math]\Phi_{ijkl}^{abcd}[/math], [math]\Phi_{ijklm}^{abcde}[/math], and [math]\Phi_{ijklmn}^{abcdef}[/math] are the second-order, third-order, fourth-order, fifth-order, and sixth-order force constants. The indices [math]i[/math], [math]j[/math], [math]k[/math], [math]l[/math], [math]m[/math], and [math]n[/math] refer to the atoms and can take integer values from 0 to [math]N-1[/math], where [math]N[/math] is the number of atoms in the system. The indices [math]a[/math], [math]b[/math], [math]c[/math], [math]d[/math], [math]e[/math], and [math]f[/math] refer to the axes in the Cartesian coordinate system and can take integer values [math]0[/math], [math]1[/math], and [math]2[/math], which correspond to the [math]x[/math], [math]y[/math], and [math]z[/math] axes, respectively. In GPUMD
, we only consider force constants up to the sixth order. One can also consider higher order force constants, but that is usually not needed.
Parameters and units
Parameter | Units |
[math]\Phi_{ij}^{ab}[/math] | eV A^{-2} |
[math]\Phi_{ijk}^{abc}[/math] | eV A^{-3} |
[math]\Phi_{ijkl}^{abcd}[/math] | eV A^{-4} |
[math]\Phi_{ijklm}^{abcde}[/math] | eV A^{-5} |
[math]\Phi_{ijklmn}^{abcdef}[/math] | eV A^{-6} |
Potential file format
One needs to prepared quite a few files related to this potential, but they can be conveniently generated by using the hiPhive package[1,2], except for the driver potential file below (which is very easy to prepare). Therefore, if you use this package to generate the input files related to the force constant potential, you don't need to read the remaining descriptions except for the driver potential file.
The driver potential file
The potential file for this potential model reads
fcp number_of_atom_types highest_order path_to_force_constant_files
fcp
is the name of this potential and tells the code that we are using the force constant potential.number_of_atom_types
is the number of atom types defined in thexyz.in
file.highest_order
is the highest order of the force constants used in the potential. For example, whenhighest_order
is 4, second-order, third-order, and fourth-order forces constants will be used (and should be prepared).path_to_force_constant_files
is the path to the force constant files (see below). Important convention: Write something like/path/to/your/folder
instead of/path/to/your/folder/
. That is, there should be no/
after the folder name.
The force constant files
- The force constant data should be prepared in some files named as
clusters_order2.in clusters_order3.in clusters_order4.in clusters_order5.in clusters_order6.in fcs_order2.in fcs_order3.in fcs_order4.in fcs_order5.in fcs_order6.in
These files should be in the folder you specified in the driver potential file (see above). If you only consider force constants up to the 4th order, you don't need the files with numbers 5 and 6.
- In each
clusters_ordern.in
file, the first line should give the number of clusters to be listed. Starting from the next non-empty line, the formats of one entry for the files fromclusters_order2.in
toclusters_order6.in
are respectively:
i j idx_fcs i j k idx_fcs i j k l idx_fcs i j k l m idx_fcs i j k l m n idx_fcs
Here idx_fcs
is the index for a force constant matrix specified in the fcs_ordern.in
files (see below). Important convention: the permutation symmetry of the force constants with respect to their atomic indices is assumed. According to this symmetry, we require that [math]i \leq j \leq k \leq l \leq m \leq n[/math].
- In each
fcs_ordern.in
file, the first line should give the number of force constant matrices to be listed. Fromfcs_order2.in
tofcs_order6.in
, each force constant matrix is of size [math] 3 \times 3[/math], [math] 3 \times 3 \times 3 [/math], [math] 3 \times 3 \times 3 \times 3[/math], [math] 3 \times 3 \times 3 \times 3 \times 3[/math], [math] 3 \times 3 \times 3 \times 3 \times 3 \times 3[/math], respectively. Taking a second order force constant matrix as an example, it reads:
0 0 phi2_00 0 1 phi2_01 0 2 phi2_02 1 0 phi2_10 1 1 phi2_11 1 2 phi2_12 2 0 phi2_20 2 1 phi2_21 2 2 phi2_22
The formats for higher order force constant matrices are similar.
- Important convention: the atomic indices take values from
0
toN-1
; the direction indices take values from0
to2
.
The equilibrium position file
Because this potential is defined in terms of the atom displacements, one has to define the equilibrium (reference) positions of the atoms in the system. A file called r0.in
is used for this purpose. This file should be in the folder you specified in the driver potential file (see above). The format of this file is:
x_0 y_0 z_0 x_1 y_1 z_1 x_2 y_2 z_2 x_3 y_3 z_3 ...
That is, each line gives the position of one atom. The order of the atoms should be consistent with that in the xyz.in
file. The coordinates are in units of angstrom.
References
- [1] The
hiPhive
package: https://hiphive.materialsmodeling.org/ - [2] Fredrik Eriksson, Erik Fransson, and Paul Erhart, The Hiphive Package for the Extraction of High‐Order Force Constants by Machine Learning, Adv. Theory. Sim., 1800184 (2019). https://doi.org/10.1002/adts.201800184