The gap_fit program

In order to fit an interatomic potential to data with gap_fit, you need

Data
A definition of descriptors and kernels that serve as basis functions
Some parameters that control the least squares fit

We will go over each of these in turn, and reference the relevant gap_fit options.

Data

The input file to gap_fit is a series of atomic structures (they could be molecules, or periodic systems), with atomic numbers, cartesian positions of the atoms, and some associated data. The required file format is extended XYZ, which is similar to a concatenation of ordinary XYZ, except that the second line (after the first line that just contains the number of atoms) is not free format, but has to be a series of key=value pairs, with some mandatory keys and other restrictions. See the detailed definition [here].

The data associated with each structure is an energy, optionally forces, and optionally virial stress. All of these are specified in the extended XYZ file. Following the definition of the extended XYZ file, the energy needs to be in units of eV, the forces in units of eV/A and the virial stress in units of eV (this is the regular stress multiplied by the volume). Watch out for the definition of the sign of the stress! For example VASP uses the opposite sign convention to QUIP. The atomic positions need to be in Angstroms, and the periodic unit cell is specified by three cartesian lattice vectors.

An example structure in extended XYZ format:

1
config_type=isolated_atom gap_energy=-157.72725320 dft_virial="0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000" dft_energy=-158.54496821 nneightol=1.20000000 pbc="T T T" Lattice="20.00000000       0.00000000       0.00000000       0.00000000      20.00000000       0.00000000       0.00000000       0.00000000      20.00000000" Properties=species:S:1:pos:R:3:Z:I:1:map_shift:I:3:n_neighb:I:1:gap_force:R:3:dft_force:R:3
Si             10.00000000     10.00000000     10.00000000      14      -1      -1      -1       0      0.00000000      0.00000000      0.00000000      0.00000000      0.00000000      0.00000000```

The structure above has lots of data in it, not all of it is useful for fitting (e.g. gap_energy happens to be a prediction by some model, not relevant for making models, or nneightol which was inserted by some earlier manupulation of the file). Note how the value of the mandatory Properties key specifies the names and types of the columns in the rest of the configuration. S stands for ‘string’, R stands for ‘real’, I is for ‘integer’ and L is for logical. Columns species and pos (atom types and cartesian positions, respectively) are mandatory.

All structures in the extended XYZ format need a periodic lattice unit cell (and defined by the ``Lattice`` key), even if a particular direction (or any direction) is not supposed to be periodic. Just use a lattice vector which is bigger than twice the cutoff of any potential you plan to create or use.

To get good fits, it is extremely helpful of the energy and the forces are consistent, i.e. up to discretisation error (e.g. in k-space in case of a DFT code), the forces are the negative derivatives of the energy. E.g. if the data comes from density functional theory with a finite electronic temperature (as is almost always used to help convergence, even in the case of insulators, when there are defects), it is the free energy (i.e. including the electronic entropy) that is consistent with the Hellman-Feynman forces, rather than the potential energy. In ASE, this is obtained by specifying the force_consistent=true option to Atoms.get_potential_energy().

Fits are tremendously improved by including forces. For 3N atoms, there is only one energy per structure, but 3N force components (obtained at almost the same cost as the single energy), so therefore contributing a lot more data. Virial stresses are again really important, especially for obtaining elastic constants for periodic solids. This is because the forces are always zero for an isotropic change away from equilibrium, and give no information on the curvature of the potential energy surface.

The type of data included can vary from configuration to configuration, there is no restriction on having to include the same type of data for all configurations.

Data field names

Conceptually, the data the potential is fitted to are the energy and its derivatives (forces and stresses), but the names with which these data are referred to in the XYZ file are arbitrary (but have to be consistent in the whole XYZ file). Unless otherwise indicated, gap_fit looks for data with keys energy, force and virial, and this can be controlled using the command line options energy_parameter name, force_parameter_name, and virial_parameter_name, respectively. This facilitates specifying several types of energies and related quantities in the same XYZ, and selected the one you want to fit to on the command line. For example, energies computed with different DFT functionals can be stored in the same XYZ file, or even with a completely different method such as QMC.

The above example has the keys dft_energy for example, so when it is used in fitting, the appropriate key names need to be specified.

Isolated atom

Most of the time, it makes sense to fit an interatomic potential to the binding energy, i.e. the total energy minus the energy of isolated atoms. There are several mechanisms in gap_fit to aid this. The default behaviour is to look for a configuration among the input structures in the input XYZ file that contains only one atom, and has an energy specified for it. This energy will then be taken and used to calculate the binding energies of all configurations before the fit is made. We give it the name e0, and it is also stored in the XML file that specifies the potential. When predictions are made, the e0 value is added back to the prediction of the binding energy. If the input XYZ contains multiple atom types, each one needs to appear once by itself among the configurations.

Alternatively, if isolated atoms are not part of the input file, a value for e0 can be specified on the command line, in case of multiple elements as a series of numbers: e0={H:123.456:O:789.10...}. There is also an option to have an e0 value be taken as the average per-atom energy of the input data. This only makes sense if there is only one type of atom.

Note that if you create a dataset where the specified energies are already binding energies, you should still specify 0.0 as the isolated atom energy (either via a configuration, preferred, or on the command line), otherwise the interatomic potential might (depending on the descriptor) not formally be zero for isolated atoms.

SOAP hyperparameters

This section contains notes about choosing hyperparameters for the SOAP descriptor.

cutoff

Every finite range potential can be cast in the form of a sum over site energies or atomic energies, and the cut-off radius defines the range of this local term. The actual interaction range is of course twice the cut-off radius, because atoms up to this distance can potentially interact with one another via a many-body term centered on an atom in between them. When we approximate a quantum mechanical potential energy (which is not formally local) using a local site energy with cut-off radius, the error we necessarily incur can be characterised in theform of a force variance.

Command line example

Here is an annotated fitting example.

gap_fit atoms_filename=database.xyz # input data in extended XYZ format
  gap={                              # start of descriptor and kernel spec
  distance_Nb                       # first descriptor is interatomic distance based
  order=2                           # descriptor is 2-body (i.e. a pair potential)
  cutoff=5.0                        # distance cutoff in the kernel, in Angstrom
  n_sparse=15                       # number of representative points, M in Sec. II
  covariance_type=ard_se            # form of kernel: squared exponential (Gaussian)
  delta=2.0                         # scaling of kernel, per descriptor, here it is per atom pair, in eV
  theta_uniform=2.5                 # length scale in Gaussian kernel in Angstrom
  sparse_method=uniform             # choice of representative points, here a uniform grid up to the cutoff
  compact_clusters=T                # how cutoff is applied, here a spherical manner around each atom
  :                                 # separator between descriptors
  soap                              # second descriptor is a SOAP
  l_max=6 n_max=12                  # number of angular and radial basis functions for SOAP
  atom_sigma=0.5                    # Gaussian smearing width of atom density for SOAP, in Angstrom
  cutoff=5.0                        # distance cutoff in the kernel, in Angstrom
  radial_scaling=-0.5               # exponent of atom density scaling, power of distance
  cutoff_transition_width=1.0       # distance across which kernel is smoothly taken to zero, in Angstrom
  central_weight=1.0                # relative weight of central atom in atom density for SOAP
  n_sparse=8000                     # number of representative points, M in Sec. II
  delta=0.2                         # scaling of kernel, per descriptor, here for SOAP it is per atom, in eV
  covariance_type=dot_product       # form of kernel
  zeta=4                            # power kernel is raised to - together with dot_product gives a polynomial kernel
  sparse_method=cur_points          # choice of representative points, here CUR decomposition of descriptor matrix
  }                                 # end of descriptor and kernel spec
 default_sigma={0.002 0.2 0.2 0.0}  # default regularisation corresponding to energy, force, virial, hessian
 config_type_sigma={                # start of per configuration-group regularisation spec, using groups defined in the input data file
  isolated_atom:0.0001:0.01:1.0:0.0:
  rss_rnd:0.03:0.4:0.5:0.0:
  rss_005:0.02:0.3:0.4:0.0:
  rss_200:0.01:0.2:0.2:0.0:
  rss_3c:0.005:0.1:0.1:0.00:
  cryst_dist:0.0003:0.03:0.05:0.00:
  cryst_dist_hp:0.005:0.1:0.1:0.0:
  liq_P4:0.003:0.3:0.5:0.0:
  liq_network:0.003:0.3:0.5:0.0:
  2D:0.001:0.03:0.05:0.0:
  ribbons:0.01:0.5:0.2:0.0
  }                                 # end of per configuration-group regularisation spec
 energy_parameter_name=energy       # name of the key in the input data file corresponding to the total energy
 force_parameter_name=forces        # name of the key in the input data file corresponding to the forces
 virial_parameter_name=virial       # name of the key in the input data file corresponding to the virial stress
 sparse_jitter=1.0e-8               # extra diagonal regulariser
 do_copy_at_file=F                  # copy input data into potential XML file?
 sparse_separate_file=T             # write representative point data into a separate file not in the main potential XML
 gp_file=gap.xml                    # name of output potential XML file
 core_param_file=P_r6_innercut.xml  # name of XML file containing the baseline potential (QUIP format)
 core_ip_args={IP Glue}             # initialisation string to call baseline potential

Command line options

See gap_fit –help for description of options.

GAP options

This does not exist any more: quippy.gap_fit_parse_gap_str.

sparse_method options are:

RANDOM: default, chooses n_sparse random datapoints
PIVOT: based on the full covariance matrix finds the n_sparse “pivoting” points
CLUSTER: based on the full covariance matrix performs a k-medoid clustering into n_sparse clusters, returning the medoids
UNIFORM: makes a histogram of the data based on n_sparse and returns a data point from each bin
KMEANS: k-means clustering based on the data points
COVARIANCE: greedy data point selection based on the sparse covariance matrix, to minimise the GP variance of all datapoints
UNIQ: selects unique datapoints from the dataset
FUZZY: fuzzy k-means clustering
FILE: reads sparse points from a file
INDEX_FILE: reads indices of sparse points from a file
CUR_COVARIANCE: CUR, based on the full covariance matrix
CUR_POINTS: CUR, based on the datapoints

Running with MPI

In order to run gap_fit with MPI, you need to configure and compile it accordingly, see the QUIP Readme for more details.

Choosing the sparse points is not parallelised, therefore you have to provide the sparse points for each descriptor in a separate file. The main workflow is:

Calculate without MPI to get the sparseX output files.
- Use sparsify_only_no_fit=T to just create the sparseX files (less memory needed).
Convert output files to input files ($QUIP_ROOT/bin/gap_prepare_sparsex_input.py gp.xml).
- Optional: Rename them to something shorter (e.g. 1.input, 2.input etc.).
Change input method (sparse_method=FILE sparse_file=1.input).
- For more than one species use add_species=F and explicit input.
- Check the number of braces (gap={{$gap1}:{$gap2}:{$gap3}}).
Run with mpirun -np … (or srun … for Slurm).