BOPcat software package for the construction and testing of tight-binding models and bond-order potentials

TL;DR

BOPcat is a Python software package that streamlines the construction and testing of tight-binding and bond-order potential models, facilitating accurate and transferable atomistic simulations.

Contribution

The paper introduces BOPcat, a modular Python tool that automates the creation and validation of TB/BOP models using a flexible, user-friendly interface and parallel processing capabilities.

Findings

Successfully constructed a transferable magnetic BOP for Fe.

Demonstrated the software's ability to handle diverse property predictions.

Validated the model against unseen crystal structures.

Abstract

Atomistic models like tight-binding (TB), bond-order potentials (BOP) and classical potentials describe the interatomic interaction in terms of mathematical functions with parameters that need to be adjusted for a particular material. The procedures for constructing TB/BOP models differ from the ones for classical potentials. We developed the BOPcat software package as a modular python code for the construction and testing of TB/BOP parameterizations. It makes use of atomic energies, forces and stresses obtained by TB/BOP calculations with the BOPfox software package. It provides a graphical user interface and flexible control of raw reference data, of derived reference data like defect energies, of automated construction and testing protocols, and of parallel execution in queuing systems. We outline the concepts and usage of the BOPcat software and illustrate its key capabilities by exemplary constructing and testing of an analytic BOP for Fe. The parameterization protocol with a successively increasing set of reference data leads to a magnetic BOP that is transferable to a variety of properties of the ferromagnetic bcc groundstate and to crystal structures that were not part of the training set.

Figures

Tables

Table 1. Table 1: Parameters of the BOP model for Fe. The unit of b emb subscript 𝑏 emb b_{\mathrm{emb}} is Å − 2 superscript italic-Å 2 \AA^{-2} and c emb subscript 𝑐 emb c_{\mathrm{emb}} is Å italic-Å \AA .

	$a$ (eV)	$b$ (${\AA }^{-1}$)	c
$dd\sigma$	-24.9657	1.4762	0.9253
$dd\pi$	21.7965	1.4101	1.0621
$dd\delta$	-2.3536	0.7706	1.3217
$U_{\mathrm{rep}}$	1797.4946	3.2809	1.0067
$U_{\mathrm{emb}}$	-1.3225	1.3374	2.1572
$N_d$	6.8876
$I\left(eV\right)$	0.9994

Table 2. Table 2: Bulk properties of ferromagnetic bcc-Fe. The experimental values for the elastic constants, vacancy formation energy and vacancy migration energy are taken from Refs.  Rayne1961 ,   DeSchepper1983 and   Vehanen1982 , respectively. The DFT values for the vacancy migration energy and surface energies are taken from Refs.   Domain2001 and   Blonski2007 , respectively.

	BOP	DFT	experiment
V/atom (Å3)	11.48	11.46	11.70
$B$ (GPa)	171	176	168
$C_{11}$ (GPa)	265	257	243
$C_{12}$ (GPa)	125	154	138
$C_{44}$ (GPa)	87	85	122
$E_f^{\mathrm{vac}}$ (eV)	2.03	2.20	2.00
$E_{\mathrm{mig}}^{\mathrm{vac}}$ (eV)	1.33	0.65	0.55
$E_f^{100}$ (eV)	3.65	4.64
$E_f^{110}$ (eV)	3.13	3.64
$E_f^{111}$ (eV)	3.59	4.34
${\gamma }_{\left(100\right)}$ (J/m2)	1.44	2.47
${\gamma }_{\left(110\right)}$ (J/m2)	1.27	2.37
${\gamma }_{\left(111\right)}$ (J/m2)	2.04	2.58
${\gamma }_{\left(211\right)}$ (J/m2)	1.50	2.50

Equations

$$U_B=U_{\mathrm{bond}}+U_{\mathrm{prom}}+U_{\mathrm{ion}}+U_{\mathrm{es}}+U_{\mathrm{rep}}+U_{\mathrm{X}}$$

$$U_{\mathrm{bond}}=2\sum _{i\alpha }{\int }^{E_F}(E-E_{i\alpha })n_{i\alpha }(E)dE$$

$${\mu }_{i\alpha }^{(n)}$$

$${\chi }^2=\sum _p\frac{w_p}{N_p}\sum _s\frac{\widetilde{w_{ps}}}{N_s^{(p)}}\frac{X_{ps}^{\mathrm{model}}-X_{ps}^{\mathrm{ref}}}{{\bar{X}}_p^{\mathrm{ref}}}.$$

$$\widetilde{w_{ps}}=\frac{1}{{\gamma }_p{N}_{\mathrm{atoms}}^{(i)}}$$

$$U_{\mathrm{binding}}=U_{\mathrm{bond}}+U_{\mathrm{rep}}+U_{\mathrm{X}}+U_{\mathrm{emb}}.$$

$${\beta }_{ij}^{\sigma }=a\exp -b{R}_{ij}^c.$$

$$U_{\mathrm{rep}}=\frac{1}{2}\sum _{i,i\ne j}{a}_{\mathrm{rep}}\exp \left(-b_{\mathrm{rep}}{R}_{ij}^{c_{\mathrm{rep}}}\right).$$

$$U_{\mathrm{X}}=-\frac{1}{4}\sum _{i}I{m}_i^2$$

$$U_{\mathrm{emb}}=-\sum _i\sqrt{\sum _{j,j\ne i}{a}_{\mathrm{emb}}\exp \left(-b_{\mathrm{emb}}(R_{ij}-c_{\mathrm{emb}}{)}^2\right)}$$

$$f_{\mathrm{cut}}=\frac{1}{2}\left[\cos \left(\pi \frac{R_{ij}-(r_{\mathrm{cut}}-d_{\mathrm{cut}})}{d_{\mathrm{cut}}}\right)+1\right]$$

Full text

I Introduction

Atomistic simulations are indispensable in materials science and require a robust description of the interatomic interaction. The application of highly accurate density functional theory (DFT) calculations is often limited by the computationally involved description of the electronic structure. A systematic coarse-graining of the electronic structure Drautz2006 ; Drautz2011 leads to the tight-binding (TB) bond model Sutton-88 and to analytic bond-order potentials (BOP) Hammerschmidt-09-IJMR ; Drautz2015 . These TB/BOP models can be evaluated at significantly reduced numerical effort and provide a transparent and intuitive model of the interatomic bond for the prediction of materials properties, see e.g. Refs. Mrovec2004 ; Mrovec2007 ; Mrovec2011 ; Seiser-11-2 ; Cak2014 ; Ford-14 ; Ford2015 ; Wang-19 .

The TB/BOP models employ adjustable parameters that need to be optimized for a particular material. This optimization is in principle comparable to the parameterization of classical potentials which can be performed with various existing software packages Brommer2007 ; Duff2015 ; Barrett2016 ; Stukowski2017 . The parameterization of TB/BOP models, however, requires sophisticated successive optimization steps, computationally-efficient handling of large data sets and an interface to a TB/BOP calculator Horsfield-96 ; Hammerschmidt2019 .

Here, we present BOPcat (Bond-Order Potential construction and testing), a software to parameterize TB/BOP models as implemented in the BOPfox software Hammerschmidt2019 . The parameters of the models are optimized to reproduce target properties like energies, forces, stresses, eigenvalues, defect formation energies, elastic constants, etc. With the interface of BOPfox to LAMMPS Plimpton-95 and ASE Larsen2017 , the list of target properties can in principle be extended to include dynamical properties. We illustrate the capability of the BOPcat software by constructing and testing an analytic BOP with collinear magnetism for Fe. Extensive tests show the good transferability of the BOP to properties which were not included in the parameterization, particularly to elastic constants, point defects, $\gamma$ surfaces, phonon spectra and deformation paths of the ferromagnetic bcc groundstate and to other crystal structures. The structures and target properties used here are taken from DFT calculations but could also include experimental data or other data sources.

We first provide a brief introduction of bond-order potentials in Sec. II. In Sec. III, the BOPcat program is outlined and the implemented modules are described. In Sec. IV, we discuss the construction of an analytic BOP for Fe and its testing. In the appendix we provide further examples of parameterization protocols for BOPcat.

II Bond-order potential formalism

Analytic BOPs provide a robust and transparent description of the interatomic interaction that is based on a coarse-graining of the electronic structure Drautz2006 ; Drautz2011 ; Drautz2015 from DFT to TB to BOP. In the following we compile the central equations and functions of the TB/BOP formalism that are parameterized in BOPcat.

In the TB/BOP models, the total binding energy is given by

[TABLE]

with the bond energy $U_{\rm{bond}}$ due to the formation of chemical bonds, the promotion energy $U_{\rm{prom}}$ from redistribution of electrons across orbitals upon bond formation, the onsite ionic energy $U_{\rm{ion}}$ to charge an atom, the intersite electrostatic energy $U_{\rm{es}}$, the exchange energy $U_{\rm{X}}$ due to magnetism and the repulsive energy $U_{\rm{rep}}$ that includes all further terms of the second-order expansion of DFT. The individual energy and force contributions are described in detail in Ref. Hammerschmidt2019 , their functional forms and according parameters for the exemplary construction of a magnetic BOP for Fe are given in Sec. IV.

The bond energy is obtained by integration of the electronic density of states (DOS) $n_{i\alpha}$,

[TABLE]

with atomic onsite levels $E_{i\alpha}$ for orbital $\alpha$ of atom $i$. In TB calculations, the DOS is obtained by diagonalization of the Hamiltonian $\hat{H}$ for a given structure. In analytic BOPs Drautz2006 ; Drautz2011 ; Drautz2015 , the DOS is determined analytically from the moments

[TABLE]

that are computed from the matrix elements $H_{i\alpha j\beta}$ of the Hamiltonian between pairs of atoms. The BOP DOS provides an approximation to the TB DOS at a computational effort for energy and force calculations that scales linearly with the number of atoms Teijeiro-16-2 . The quality of the BOP approximation can be improved systematically by including higher moments with a power-law scaling for the increase in computational effort Teijeiro-16-1 . A detailed introduction to TB/BOP calculations in BOPfox is given in Ref. Hammerschmidt2019 .

III Program Flow

BOPcat is a collection of Python modules and tools for the construction and testing of TB/BOP models. The individual steps are specified by the user as a protocol in terms of a Python script. The overall construction and testing proceeds as follows:

• •

define and initialize input controls (CATControls)

• •

generate initial model (CATParam)

• •

read reference data (CATData)

• •

initialize calculator interface (CATCalc)

• •

build and run optimization kernel (CATKernel)

The individual tasks are modularized in the modules CATControls, CATParam, CATData, CATCalc and CATKernel that interact as illustrated in Fig. 1. Due to its modular structure, one can also use only a subset of the modules, e.g., for successive optimizations.

In the following, the individual modules are described in detail with snippets of the python protocol for a basic construction process of a magnetic BOP for Fe. Further example protocols are given in the appendix.

III.1 Input controls: CATControls

The variables required for running BOPcat are defined in the module CATControls. These include the list of chemical symbols for which the model is constructed, calculator settings, filters for the reference data, optimization variables, specifications of the model, etc.. The consistency of the input variables is checked within the module to catch inconsistent settings at an early stage. The CATControls object is then passed on to succeeding modules from which the relevant variables are assigned. As an example initialization of the CATControls object we define the functions of the TB/BOP model and the initial guess for the according parameters by providing an existing model from a file (models.bx).

from bopcat.catcontrols import CATcontrols ctr = CATControls() ctr.elements = [’Fe’] ctr.model_pathtomodels = ’models.bx’ ctr.calculator_nproc = 4 ctr.calculator_parallel = ’multiprocessing’ ctr.calculator_settings = {’scfsteps’:500} ctr.data_filename = ’Fe.fit’ ctr.data_free_atom_energies = {’Fe’:-0.689} ctr.data_system_parameters = {’spin’:[1,2], ’system_type’:[0]} ctr.opt_optimizer = ’leastsq’ ctr.opt_variables = [{’bond’:[’Fe’,’Fe’], ’rep1’:[True,True]}]

Alternatively, BOPcat can generate an initial guess from raw bond integrals Jenke2019 with additional user-specified information on the functions to be used for the individual energy contributions, the valence-type ($s$/$p$/$d$) and number of valence electrons of the elements, the cut-off radii for determining pairs of interacting atoms, etc. (see  B)

III.2 Reference data: CATData

The set of reference data of structures and properties is read from a text file by the CATData module. The entries for reference data have the following format:

data_type = 0 code = 7 basis_set = 0 pseudopotential = 30 strucname = bcc a1 = -1.3587 1.3587 1.3587 a2 = 1.3587 -1.3587 1.3587 a3 = 1.3587 1.3587 -1.3587 coord = cartesian Fe 0.0000 0.0000 0.0000 energy = -7.9561 forces 0.0000 0.0000 0.0000 stress = -0.08423 -0.0843 -0.0843 -0. -0. -0. system_type = 0 space_group = 229 calculation_type = 1 stoichiometry = Fe calculation_order = 10.033 weight = 1.0 spin = 1

In contrast to optimization procedures for classical potentials, the optimization of TB/BOP models with BOPcat can be carried out not only for atomic properties but also at the electronic-structure level by providing eigenvalues and corresponding $k$-points as reference data. The optimization can be further extended to other properties of TB/BOP calculations (e.g. band width, magnetic moment) by supplying the according keyword and the data type of the BOPfox software in the reference data module of BOPcat.

CATData stores the reference data in a set of Atoms objects of the ASE python framework for atomistic simulations Larsen2017 . Each of the ASE Atoms objects can have calculator-specific identifiers such as the settings of the DFT calculation (code, basis_set, pseudopotential) as well as structure-related identifiers such as calculation type (calculation_type) and calculation order (calculation_order), which are used as filters to use only a subset of the available reference data. A complete description of the implemented identifiers and filters is given in the manual of BOPcat. The corresponding reference data can include basic quantities like energies, forces, stresses as well as derived quantities such as defect energies. Extracting the reference structures and their corresponding properties is illustrated in the following:

from bopcat.catdata import CATData data = CATData(controls=controls) ref_atoms = data.get_ref_atoms( structures=[’bcc*’,’fcc*’,’hcp*’], quantities=[’energy’], sort_by=’energy’) ref_data = data.get_ref_data()

The user-specified reference data, i.e. structures and properties, are passed to the calculator and optimization kernel, respectively.

III.3 Initial model: CATParam

The construction of TB/BOP models in BOPcat relies on iterative optimization steps that require an initial guess that can be provided by the user or generated by the CATParam module. In the latter case the functions and parameters are initialized from a recently developed database of TB bond integrals from downfolded DFT eigenspectra Jenke2019 . In the following snippet, the initial model is read from the filename specified in the input controls.

from bopcat.catparam import CATParam param = CATParam(controls=controls) ini_model = param.models[0]

The CATParam module also provides meta-operations on sets of models such as identifying the optimum model for a given set of reference data.

III.4 Calculator: CATCalc

After the reference structures are read and the initial model is generated, the CATCalc module sets up the calculator for the computation of the specified properties for the reference structures by the initial and then optimized TB/BOP model. To this end, a list of BOPfox-ASE calculators is constructed which can be run in serial or in multiprocessing mode. The latter is optimized by a load balancing scheme that is based on the size of the structures. Additional properties for the model optimization can easily be included by extending this module accordingly. In the following snippet, an initial model is provided and an ASE calculator is initialized for each of the reference structures.

from bopcat.catcalc import CATCalc calc = CATCalc(controls=ctr, model=ini_model, atoms=ref_atoms)

Here, the input controls (ctr), the initial guess for the TB/BOP model (ini_model) and the reference structures (ref_atoms) are specified in the code snippets given in Sec. III.1, III.3 and III.2, respectively.

III.5 Construction and testing kernel: CATKernel

The reference data and associated calculators are then used to set up the CATKernel module. This module provides an interface to the objective function for the construction and testing of TB/BOP models. The default objective function for $N_{p}$ properties of $N_{s}$ structures is given by

[TABLE]

The user can choose other definitions of the objective function that are implemented in BOPcat or provide external implementations. The weights for the individual structures

[TABLE]

are specified by the user or determined by the module. The dimensionality factor $\gamma_{p}$ balances the relative weights of energies, forces and stress by taking values of $1,3,6$, respectively. The objective function can also be normalized by the variance of the properties as in Ref. Krishnapriyan2017 . For construction purposes, this module provides an interface between objective function and optimizer algorithms that are either readily available in the Python modules Scipy Scipy and NLopt Nlopt or provided by the user as external module.

In the following example we illustrate the initialization and running of the CATkernel object for the input controls (ctr), calculator (calc) and reference properties (ref_data) defined in Sec. III.1, III.4 and III.3.

from bopcat.catkernel import CATKernel kernel = CATKernel(controls=ctr, calc=calc, ref_data) kernel.run()

III.6 Parallel execution of BOPcat

The optimization kernel can be executed efficiently in parallel by distributing the required property calculations with Python subprocesses and the message passing interface. For this purpose, BOPcat also features a process management module to interact with the queuing system of compute clusters. This allows high-throughput optimizations of TB/BOP models with, e.g., different initial models or different sets of reference data. We illustrate this in the following snippet:

from bopcat.process_management
import Process_catkernels from bopcat.process_management
import Process_catkernel from bopcat.process_management import queue models = [ini_model.rattle(kernel.variables) for i in range(5)] subprocs = [] for i in range(len(kernels)): ckern = kernel.copy() ckern.calc.set_model(models[i]) proc = Process_catkernel(catkernel=ckern, queue=queue) subprocs.append(proc) proc = Process_catkernels(procs=subprocs) proc.run()

A list of random models (models) is first generated by applying noise to the parameters of the initial model (ini_model). A copy of the CATKernel object is then generated (ckern) and a model from the list is assigned. The kernel is then serialized by the Process_catkernel object which contains the details on how to execute the kernel. Finally, the individual subprocesses are wrapped in the main object Process_catkernels which submits them to a specified queue (queue) of the compute cluster.

III.7 Graphical user interface

The setup of a construction or testing process with BOPcat discussed so far is based on writing and editing Python scripts. Alternatively, BOPcat can be driven with a graphical user interface (GUI), see snapshots in Fig. 2.

Basic information of a construction such as the current value of the parameters and objective function can be visualized. The reference data can be inspected visually and corresponding attributes such as the weights can be set directly. For testing purposes, a given TB/BOP model can be evaluated and visualized for a given set of reference data. The GUI also comes with a parameterization builder which allows one to design and execute a parameterization protocol. In particular, one can add a number of optimization layers that consist of several optimization kernels with communication of intermediate TB/BOP models between subsequent layers. This provides the technical prerequisite for automatizing the construction of TB/BOP models by sophisticated protocols.

IV Construction and testing of simple Fe BOP

IV.1 Functional form

As an example application of BOPcat we construct a simple BOP model for Fe including magnetism. We note that this Fe BOP is constructed to demonstrate BOPcat and that capturing the reference data more precisely would require a more complex functional form. The atoms are assumed to be charge-neutral so that the binding energy given in Eq. 1 reduces to

[TABLE]

We choose a $d$-valent model including only two-center contributions, similar to previous TB/BOP models for Fe Mrovec2011 ; Ford2015 ; Madsen2011 , with an initial number of $d$-electrons of $N_{d}=6.8$. We use a BOP expansion up to the 9-th moment (Eq. 3) as a good compromise between performance and accuracy of the DOS.

The Hamiltonian matrix elements $H_{i\alpha j\beta}$ which enter the calculation of the moments are constructed from Slater-Koster bond integrals $\beta_{ij}^{\sigma}$ with the angular character $\sigma=dd\sigma,dd\pi,dd\delta$. The initial guess of the bond integral parameters is taken from projections of the DFT eigenspectrum on a TB minimal basis for the Fe-Fe dimer Jenke2019 . The distance dependence of the bond integrals is represented by a simple exponential function

[TABLE]

For this simple BOP model, the repulsive energy is given as a simple pairwise contribution.

[TABLE]

The magnetic contribution to the energy $U_{X}$ is evaluated using

[TABLE]

with $m_{i}=N_{i}^{\uparrow}-N_{i}^{\downarrow}$ the magnetic moment of atom $i$ and $I$ the Stoner exchange integral that is initially set to 0.80 eV similar to Ref. Mrovec2011 . In extension to Eq. 1 we use an additional empirical embedding contribution in Eq. 6

[TABLE]

which may be understood as providing contributions due to the missing $s$ electrons and non-linear exchange correlation.

The bond integrals, repulsive and embedding functions are multiplied with a cut-off function

[TABLE]

in the range of $[r_{\mathrm{cut}}-d_{\mathrm{cut}}$, $r_{\mathrm{cut}}]$ in order to restrict the range of the interatomic interaction. For the Fe BOP constructed here we used $r_{\mathrm{cut}}=3.8~{}\AA$, $d_{\mathrm{cut}}=0.5~{}\AA$ for the bond integrals and $r_{\mathrm{cut}}=6.0~{}\AA$, $d_{\mathrm{cut}}=1.0~{}\AA$ for the repulsive and embedding energy terms.

IV.2 Reference data

The reference data for constructing and testing the Fe BOP is obtained by DFT calculations using the VASP software Kresse-96-1 ; Kresse-96-2 with a high-throughput environment Hammerschmidt-13 . We used the PBE exchange-correlation functional Perdew-96 and PAW pseudopotentials Bloechl-94 with $p$ semicore states. A planewave cut-off energy of 450 eV and Monkhorst-Pack $k$-point meshes Monkhorst-76 with a density of 6 /Å*-1* were sufficient to converge the total energies to 1 meV/atom. The reference data for constructing the Fe BOP comprised the energy-volume curves of the bulk structures bcc, fcc, hcp, and the topologically close-packed (TCP) phases A15, $\sigma$, $\chi$, $\mu$, C14, C15 and C36 that are relevant for Fe compounds Ladines2015 .

For testing the Fe BOP, several additional properties were determined that are related to the ferromagnetic bcc groundstate structure. In particular, the elastic constants at the equilibrium volume were determined by fitting the energies as function of the relevant strains Golesorkhtabar2013 . The formation energies of point defects were calculated with the supercell approach with fixed volume corresponding to the bulk equilibrium volume. A $6\times 6\times 6$ supercell was found to be sufficient to converge the formation energies to an uncertainty of 0.1 eV. The energy barriers for vacancy migration were calculated with the nudged elastic band method. The phonon spectra are computed with the Phonopy software Togo2015 . The transformation paths from bcc to other crystal structures were determined according to Ref. Paidar1999 .

IV.3 Construction

The parameters of the Fe BOP with the functions and reference data given above were optimized with the following BOPcat protocol:

The magnitudes of the repulsive and embedding terms $a_{\mathrm{rep}}$ and $a_{\mathrm{emb}}$ are adjusted to reproduce the energies of hydrostatically deformed bcc, fcc and hcp with fixed values for the exponents taken from Ref. Madsen2011 . 2. 2.

The prefactors of the bond integrals are adjusted including the energies of the TCP phases. 3. 3.

The exponents of the repulsive functions are optimized including the randomly deformed structures in the reference set. 4. 4.

The exponents of the bond integrals are optimized by increasing the size of the reference data. 5. 5.

The other parameters are optimized further while increasing the number of reference structures.

The resulting parameters of the model are compiled in Tab. 1.

The initial and optimized bond integrals, and the empirical potentials are plotted in Fig. 3. There is no substantial change in the bond integrals except that the bond integrals become longer-ranged after optimization. At shorter distance the bond integrals become weaker which can be rationalized by the screening influence of the neighboring atoms in the bulk structure. The effective number of $d$ electrons and the Stoner integral increased during optimization.

The optimized pairwise repulsive term $U_{\mathrm{rep}}$ shown in Fig. 4 is repulsive for all distances. The embedding term $U_{\mathrm{emb}}$, in contrast, is negative for all distances.

IV.4 Testing

IV.4.1 Properties related to bcc-Fe

An important basic test of the Fe BOP is the performance for properties that are closely related to the ferromagnetic bcc groundstate structure at 0 K. In Tab. 2 we compare elastic constants as well as point defect and surface properties as predicted by the Fe BOP to available DFT and experimental values. The overall good agreement demonstrates the robust transferability of the Fe BOP to properties not included in the training set.

The predicted elastic constants $C_{11}$, $C_{12}$ and $C_{44}$ are in good agreement with DFT although the specific deformed structures were not included in the reference set for constructing the BOP. The prediction for the vacancy formation energy $E_{f}^{\mathrm{vac}}$ is also in good agreement with DFT while the vacancy migration energy $E_{\mathrm{mig}}^{\mathrm{vac}}$ is overestimated as in previous TB/BOP models for Fe Madsen2011 ; Mrovec2011 . The formation energies of the vacancy and self-interstitial atoms calculated by the Fe BOP exhibit the correct energetic ordering although they were not included in the parameterization. The absolute values are slightly underestimated as compared to DFT, similar to previous models Madsen2011 ; Mrovec2011 . The relative stability of the low-index surfaces are also satisfactorily reproduced by the present Fe BOP. However, similar to the previous BOP model of MrovecMrovec2011 , the energy difference of the other surfaces relative to $(110)$ is overestimated. The deviations for point defects and surfaces are attributed to the missing $s$ electrons in the model and the lack of screening in the orthogonal TB model.

As a test towards finite-temperature applications we determine the phonon bandstructure for ferromagnetic bcc Fe shown in Fig. 5.

The prediction of the Fe BOP is in good overall agreement with DFT and experimental results. The close match near the $\Gamma$ point is expected from the good agreement of the elastic constants. The experimental data were obtained at 300 K which can explain the deviation of the experimental from the DFT/BOP data for the $T2$ branch along $T-N$.

The transferability of the model in the case of large deformations of the bcc groundstate structure is tested using deformation paths that connect the high symmetry structures bcc, fcc, hcp, sc and bct. The energy profile versus the deformation parameter for the tetragonal, hexagonal, orthogonal and trigonal deformation paths is shown in Fig. 6.

The energy profiles for all transformation paths are predicted very well by the Fe BOP model. The energies at the high symmetry points are in good agreement with DFT except for the sc structure which is slightly underestimated by BOP.

IV.4.2 Transferability to other crystal structures

In order to verify the transferability of the Fe BOP model to other crystal structures, we determined the equilibrium binding energy, volume and bulk modulus for a broader set of crystal structures as shown in Fig. 7.

The properties of the bcc, fcc, hcp, A15, $\sigma$, $\chi$, $\mu$, C14, C15 and C36 structures that were included in the parameterization are reproduced with excellent agreement. This shows that the chosen functional form of the Fe BOP is sufficiently flexible to adapt to this set of reference data. The Fe BOP model also shows robust predictions of equilibrium binding energy, volume and bulk modulus across the set of tested structures. The open crystal structures (e.g. A4) show comparably larger errors which can be expected as the reference set used in the parameterization covers mostly close-packed local atomic environments Jenke2018 while local atomic environments of open structures were not part of the training set.

V Conclusions

We present the BOPcat software for the construction and testing of TB/BOP models as implemented in the BOPfox code. TB/BOP models are parameterized to reproduce reference data from DFT calculations including energies, forces, and stresses as well as derived properties like defect formation energies. The modular framework of BOPcat allows one to implement complex parameterization protocols by flexible python scripts. BOPcat provides a graphical user interface and a highly-parallelized optimization kernel for fast and efficient handling of large data sets.

We illustrate the key features of the BOPcat software by constructing and testing a simple $d$-valent BOP model for Fe including magnetism. The resulting BOP predicts a variety of properties of the groundstate structure with good accuracy and shows good quantitative transferability to other crystal structures.

Acknowledgements

The authors are grateful to Aparna Subramanyam, Malte Schröder, Alberto Ferrari, Jan Jenke, Yury Lysogorskiy, Miroslav Cak, and Matous Mrovec for discussions and feedback using BOPcat. We acknowledge financial support by the German Research Foundation (DFG) through research grant HA 6047/4-1 and project C1 of the collaborative research centre SFB/TR 103. Part of the work was carried out in the framework of the International Max-Planck Research School SurMat.

Appendix

We provide several examples to illustrate the script required to direct an optimization procedure. The examples show only the primary capabilities of BOPcat. For a complete list of possible functionalities of the software, the user is referred to the user manual.

Appendix A Constructing Fe-Nb EAM

Aside from parameterizing TB models and BOPs, the software is also capable of constructing empirical potentials. We illustrate this in the following script which can be used to parameterize an EAM for the Fe-Nb binary. The model consists of pair potentials described by polynomials, and square root embedding functions. The optimization proceeds in four steps, each parameterizing different interaction pairs (1) Fe-Fe, (2) Nb-Nb, (3) Fe-Nb and (4) all. In the first three steps, the reference data consists only of energies, while the last step include structures with energy, forces and stress.

from bopcat.catcontrols import CATControls from bopcat.catdata import CATData from bopcat.catparam import CATParam from bopcat.catcalc import CATCalc from bopcat.catkernel import CATKernel import bopcat.functions as Funcs

#################### Input controls #################### cat_controls = CATControls()

required, list of elements

cat_controls.set(elements = [’Fe’,’Nb’])

Model

cut_offs

cat_controls.set(model_cutoff =
{’r2cut’: {’Fe-Fe’:6.0, ’Nb-Nb’:6.0, ’Fe-Nb’:6.0} ,’d2cut’: {’Fe-Fe’:0.5, ’Nb-Nb’:0.5, ’Fe-Nb’:0.5}})

functional forms

cat_controls.set(model_functions =
{’rep1’:{’Fe-Fe’:Funcs.polynom4() ,’Nb-Nb’:Funcs.polynom4() ,’Fe-Nb’:Funcs.polynom4()} ,’rep2’:{’Fe-Fe’:Funcs.sqrt_polynom3 ,’Nb-Nb’:Funcs.sqrt_polynom3 ,’Fe-Nb’:Funcs.sqrt_polynom3}})

Calculator

dictionary of keyword-values used by calculator

cat_controls.set(calculator_settings =
{’version’:’eam’})

Reference Data

database filters

cat_controls.set(data_parameters =
{’xc_functional’:’PBE’,’encut’:450})

filename

cat_controls.set(data_filename = ’FeNb.fit’)

optimizer

cat_controls.set(opt_optimizer_options =
{’max_nfev’:100,’diff_step’:0.01}) cat_controls.set(opt_optimizer = ’least_squares’)

others

cat_controls.set(verbose = 2)

#################### Initialization ####################

generate reference data

cat_data = CATData(controls=cat_controls)

generate initial model

cat_param = CATParam(controls=cat_controls)

set up calculator

cat_calc =
CATCalc(controls=cat_controls ,model=cat_param.models[-1])

################### Optimization Fe ####################

optimize Fe

fetch reference structures and properties

structures can be identified by structure name or

by ’stoich/sgrp/sys-type/calc-type/calc-order’

ref_atoms =
cat_data.get_ref_atoms(structures=[’Fe*//0//*’] ,quantities=[’energy’]) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms)

optimization variables

var = [{"bond":["Fe","Fe"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]}]

set up optimization kernel

optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls)

run optimization

optfunc.optimize()

update models

new_model = optfunc.optimized_model cat_param.models.append(new_model)

################### Optimization Nb ####################

start from previous model

cat_calc.set_model(new_model) ref_atoms =
cat_data.get_ref_atoms(structures=[’Nb*//0//*’] ,quantities=[’energy’]) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms) var = [{"bond":["Nb","Nb"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]}] optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls) optfunc.optimize() new_model = optfunc.optimized_model cat_param.models.append(new_model)

################### Optimization Fe-Nb ###################

start from previous model

cat_calc.set_model(new_model) ref_atoms =
cat_data.get_ref_atoms(structures=[’FeNb//0//*’] ,quantities=[’energy’]) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms) var = [{"bond":["Fe","Nb"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]}] optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls) optfunc.optimize() new_model = optfunc.optimized_model cat_param.models.append(new_model)

#################### Optimization all ####################

start from previous model

cat_calc.set_model(new_model) ref_atoms =
cat_data.get_ref_atoms(structures=[’FeNb//0//*’] ,quantities=[’energy’,’forces’ ,’stress’]) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms) var = [{"bond":["Fe","Fe"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]} ,{"bond":["Nb","Nb"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]} ,{"bond":["Fe","Nb"],’rep1’:[True,True,True,True] ,’rep2’:[True,True,True]}] optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls) optfunc.optimize() new_model = optfunc.optimized_model cat_param.models.append(new_model)

write out new model

new_model.write(filename=’models_optimized.bx’)

Appendix B Construction of BOP from raw TB bond integrals

A BOP or TB model can also be constructed directly from the numerical data for the bond integrals using BOPcat. One simply needs to supply the function to fit the data. In the following example, a BOP is constructed for Fe. The electronic part of the potential is parameterized with respect to the DFT eigenvalues while the other terms are optimized by fitting to the energies of the bcc structures.

from bopcat.catcontrols import CATControls from bopcat.catdata import CATData from bopcat.catparam import CATParam from bopcat.catcalc import CATCalc from bopcat.catkernel import CATKernel import bopcat.functions as Funcs from bopcat.eigs import get_relevant_bands

#################### Input controls #################### cat_controls = CATControls()

required, list of elements

cat_controls.set(elements = [’Fe’])

Model

cut_offs

cat_controls.set(model_cutoff =
{’rcut’ : {’Fe-Fe’:4.0},’r2cut’: {’Fe-Fe’:6.0} ,’dcut’ : {’Fe-Fe’:0.5},’d2cut’: {’Fe-Fe’:0.5}})

TB basis

cat_controls.set(model_orthogonal = True) cat_controls.set(model_valences={’Fe’:’d’}) cat_controls.set(model_betabasis=’tz0’) cat_controls.set(model_betatype=’loewdin’)

filenames to bond integrals and onsites

cat_controls.set(model_pathtobetas=’../betas’) cat_controls.set(model_pathtoonsites=’../onsites’)

functional forms

cat_controls.set(model_functions =
{’ddsigma’:{’Fe-Fe’:Funcs.exponential()} ,’ddpi’:{’Fe-Fe’:Funcs.exponential()} ,’dddelta’:{’Fe-Fe’:Funcs.exponential()} ,’rep1’:{’Fe-Fe’:Funcs.exponential()} ,’rep2’:{’Fe-Fe’:Funcs.sqrt_gaussian()}})

Calculator

dictionary of keyword-values used by calculator

cat_controls.set(calculator_settings =
{’bandwidth’:’findeminemax’ ,’terminator’:’constantabn’,’bopkernel’:’jackson’ ,’nexpmoments’:100,’moments’:9,’version’:’bop’})

Reference Data

database filters

cat_controls.set(data_parameters =
{’xc_functional’:’PBE’,’encut’:450})

filename

cat_controls.set(data_filename = ’Fe.fit’)

optimizer

cat_controls.set(opt_optimizer_options =
{’max_nfev’:100,’diff_step’:0.01}) cat_controls.set(opt_optimizer = ’least_squares’)

others

cat_controls.set(verbose = 2)

#################### Initialization ####################

generate reference data

cat_data = CATData(controls=cat_controls)

generate initial model

cat_param = CATParam(controls=cat_controls)

set up calculator

cat_calc =
CATCalc(controls=cat_controls ,model=cat_param.models[-1])

################## Optimization elec ################### ref_atoms =
cat_data.get_ref_atoms(structures=[’Fe/229/0/0/*’] ,quantities=[’eigenvalues’]) ref_data = cat_data.get_ref_data() orb_char = cat_data.orb_char cat_calc.set_atoms(ref_atoms)

select only relevant bands in dft eigs

ref_data, orb_char =
get_relevant_bands(cat_calc,ref_data,orb_char)

optimize bond integrals

var = [{"bond":["Fe","Fe"],"ddsigma":[True,True] ,"ddpi":[True,True], "dddelta":[True,True] ,’atom’:’Fe’,’onsitelevels’:[True]}]

set up optimization

optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls ,weights=orb_char)

run optimization

optfunc.optimize()

update models

new_model = optfunc.optimized_model cat_param.models.append(new_model)

################## Optimization rep #################### ref_atoms =
cat_data.get_ref_atoms(structures=[’Fe/229/0/1/*’] ,quantities=[’energy’ ,’vacancy_energy’] ) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms)

start from previous model

cat_calc.set_model(new_model)

optimization variables

var = [{"bond":["Fe","Fe"],’rep1’:[True,True] ,’rep2’:[True,True]}]

set up optimization

optfunc = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log.cat’ ,controls=cat_controls)

run optimization

optfunc.optimize()

update models

new_model = optfunc.get_optimized_model() cat_param.models.append(new_model)

write out new model

new_model.write(filename=’models_optimized.bx’)

Appendix C Multiprocessing and usage on clusters

In order to speed up the optimization process, BOPcat can also be executed on clusters. The following script illustrates this by optimizing a BOP for Re with respect to different sets of reference data. A kernel for each case is generated and serialized by the Process_catkernel object. These subprocesses are then wraped in a main process using the Process_catkernels object which facilitates the submission to the queing system. A unique queuing system can be easily set up by extending the queue object.

from bopcat.catcontrols import CATControls from bopcat.catdata import CATData from bopcat.catparam import CATParam from bopcat.catcalc import CATCalc from bopcat.catkernel import CATKernel from process_management import Process_catkernel,vulcan from process_management import Process_catkernels

#################### Input controls #################### cat_controls.set(elements = [’Re’])

model name

cat_controls.set(model =
"Cak-2014") cat_controls.set(model_pathtomodels = ’models.bx’) cat_controls.set(calculator_settings = {’scfsteps’:100}) cat_controls.set(data_parameters =
{’xc_functional’:’PBE’,’encut’:450, ’spin’:1})

reference data

cat_controls.set(data_filename = ’Re.fit’)

free atom energies

cat_controls.set(data_free_atom_energies = {’Re’:-0.5})

optimizer

cat_controls.set(opt_optimizer_options = {’maxiter’:100}) cat_controls.set(opt_optimizer = ’nelder-mead’)

others

cat_controls.set(calculator_nproc = 16) cat_controls.set(calculator_parallel = ’multiprocessing’) cat_controls.set(verbose = 2)

#################### Initialization ####################

generate reference data

cat_data = CATData(controls=cat_controls)

generate initial model

cat_param = CATParam(controls=cat_controls)

set up calculator

cat_calc =
CATCalc(controls=cat_controls ,model=cat_param.models[-1])

####################### Kernels ######################### var = [{"bond":["Re","Re"],’rep1’:[True,False,False,True ,False,False,False,False],’rep2’:[True,True,True ,True,True,False,False,False,False]},{"atom":"Re" ,’valenceelectrons’:[True]}]

use different reference structures for each kernel

strucs_list = [[’bcc*’,’hcp*’],[’fcc*’,’hcp*’] ,[’A15*’,’hcp*’],[’C14*’,’hcp*’]] kernels = [] for i in range(len(strucs_list)): ref_atoms =
cat_data.get_ref_atoms(structures=strucs_list[i] ,quantities=[’energy’]) ref_data = cat_data.get_ref_data() cat_calc.set_atoms(ref_atoms) kern = CATKernel(calc=cat_calc,ref_data=ref_data ,variables=var,log=’log_%d.cat’%i ,controls=cat_controls) kernels.append(kern)

####################### Queues ########################## queue = vulcan(cores=16,qname=’parallel16.q’ ,pe=’mpi16’,runtime=606024*14) subprocs = [] for kern in kernels: subprocs.append(Process_catkernel(catkernel=kern ,queue=queue,directives=[’source ~/.bashrc’]))

proc = Process_catkernels(procs=subprocs) proc.run() res = proc.get_results(wait=True) for i in range(len(proc._procs)): kern = proc._procs[i].get_kernel()._catkernel opt_model - kern.get_optimized_model() cat_param.models.append(opt_model)

ave_model = cat_param.average_modelsbx() ave_model.write(filename=’models_optimized.bx’) proc.clean()