# The force constant potential

## 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 and r_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 $\{u_i^{a}\}$ 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, $\Phi_{ij}^{ab}$, $\Phi_{ijk}^{abc}$, $\Phi_{ijkl}^{abcd}$, $\Phi_{ijklm}^{abcde}$, and $\Phi_{ijklmn}^{abcdef}$ are the second-order, third-order, fourth-order, fifth-order, and sixth-order force constants. The indices $i$, $j$, $k$, $l$, $m$, and $n$ refer to the atoms and can take integer values from 0 to $N-1$, where $N$ is the number of atoms in the system. The indices $a$, $b$, $c$, $d$, $e$, and $f$ refer to the axes in the Cartesian coordinate system and can take integer values $0$, $1$, and $2$, which correspond to the $x$, $y$, and $z$ 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 $\Phi_{ij}^{ab}$ eV A-2 $\Phi_{ijk}^{abc}$ eV A-3 $\Phi_{ijkl}^{abcd}$ eV A-4 $\Phi_{ijklm}^{abcde}$ eV A-5 $\Phi_{ijklmn}^{abcdef}$ 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 the xyz.in file.
• highest_order is the highest order of the force constants used in the potential. For example, when highest_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 from clusters_order2.in to clusters_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 $i \leq j \leq k \leq l \leq m \leq n$.

• In each fcs_ordern.in file, the first line should give the number of force constant matrices to be listed. From fcs_order2.in to fcs_order6.in, each force constant matrix is of size $3 \times 3$, $3 \times 3 \times 3$, $3 \times 3 \times 3 \times 3$, $3 \times 3 \times 3 \times 3 \times 3$, $3 \times 3 \times 3 \times 3 \times 3 \times 3$, 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 to N-1; the direction indices take values from 0 to 2.

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