The force constant potential

From GPUMD
Jump to navigation Jump to search

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.

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 fcn.in, where n can be 2, 3, 4, 5, and 6. If highest_order is n, files from fc2.in to fcn.in need to be prepared. These files should be in the folder you specified in the driver potential file (see above).
  • In each fcn.in file, the first line should give the number of force constants to be listed. Starting from the second line, each line gives one force constant and its indices. The formats of one entry for the files from fc2.in to fc6.in are respectively:
    i j a b Phi_ij_ab
    i j k a b c Phi_ijk_abc
    i j k l a b c d Phi_ijkl_abcd
    i j k l m a b c d e Phi_ijklm_abcde
    i j k l m n a b c d e f Phi_ijklmn_abcdef
  • Important convention: the atomic indices take values from 0 to N-1; the direction indices take values from 0 to 2.
  • 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 [/math] in file fc2.in
    • [math]i \leq j \leq k [/math] in file fc3.in
    • [math]i \leq j \leq k \leq l [/math] in file fc4.in
    • [math]i \leq j \leq k \leq l \leq m [/math] in file fc5.in
    • [math]i \leq j \leq k \leq l \leq m \leq n[/math] in file fc6.in

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

How to prepare the potential files?

  • To prepare the above-mentioned potential files, we recommend to use the hiPhive package: https://hiphive.materialsmodeling.org/
  • The hiPhive team has implemented functions which can be used to output the files required for this potential.