In a generalised form they are equal to Finnis-Sinclair, effective medium theory or glue potentials. Their total energy can be written as

$E = \sum_{ij}V(r_{ij}) + \sum_i F(\rho_i)$

with

$\rho_i = \sum_j \rho(r_{ij})$

The original functions for V, $\rho$ and F were derived from different theories, but they can be chosen freely.

Fitting is done using atomicrex https://atomicrex.org. In the fit process an objective or cost function is minimized. The objective function is defined as

$\chi^2 = \sum_i w_i r_i$

where $w_i$ is a weight and $r_i$ is a residual that describes the difference to target values. This residual can be defined in different ways, so it is not possible to simply compare the residual for different fitting processes or codes. A more in depth explanation and some examples can be found on https://atomicrex.org/overview.html#objective-function.

It is possible to assign different weights to certain structures or properties, depending on what should be investigated using the potential. Here every structure has the same weight, but the force vector with N*3 values is normalized to have the same total weight as the single value energy. Therefore it is divided by the number of atoms.

It is necessary to define a pair potential, an electronic density function and an embedding function.

For all of those it is possible to choose between different functional forms.

Classic pair potentials are physically motivated and have a very limited number of paramaters that are derived from a experimentally measured quantity.

Splines or polynomials offer more flexibility, but can lead to unphysical oscillations or overfitting. Compared with the machine learning potentials shown later the number of parameters is very low no matter which functions you choose and the problem is highly non linear.

In this case a generalized morse function is used for the pair interaction. It has the form

The parameters in the morse potential can be derived from phyiscal quantities, but in this case they are just educated guesses. For example $r_0$ is the equilibrium distance of a dimer. The nearest neighbor distance in fcc Cu is about 2.5 $\mathring A$ so it is taken as initial value.

In the case of analytic functions the initial parameter choices should not matter too much, since the functional form is constrained.

The electronic density and embedding function will be splines. Depending on the properties that are calculated other functional forms could give better results. The inital parameters require more testing and hand tuning than the parameters of analytic functions.

The electron density is chosen to be a spline function. The cutoff has to be defined. Derivatives left and right are optional, they default to 0. For the right cutoff this is fine, since the forces should smoothly go to 0. For the left this is not necessarily the best choice, since the function value should increase at very close distances. Very large absolute values will lead to osciallations and should be avoided.

For a spline function it is necessary to define node points. They can be equally spaced or sampled with higher density around turning points, f.e. the first neighbor distance.

Too few nodepoints mean low flexibilty, too many lead to overfitting. This requires some testing to find an optimal choice.

%% Cell type:code id:hidden-wildlife tags:

``` python

rho_nodes=np.linspace(0.5,7.0,7).round(2)

rho_nodes

```

%%%% Output: execute_result

array([0.5 , 1.58, 2.67, 3.75, 4.83, 5.92, 7. ])

%% Cell type:markdown id:binary-devil tags:

The nodes need initial values. The electron density should be proportional to $e^{-r}$, so this function is chosen to calculate them.

Additionally it is a good idea to define limits for the node points. This is optional for local minimizers, but the fit can quickly run away without limits. Global optimizers typically require them to constrain the sampled space.

A density can't be negative so the lower limit is set to 0. The upper limit is chosen to be 3 times the initial values. These choices aswell as the choice for $e^{-r}$ as initial values are somewhat arbitrary, but don't matter much. The electron density from single atoms does not directly influence the calculated energies and forces, instead the summed up density at some place is used in the embedding function, so the final numerical values are an interplay between electron density and embedding function. Since the latter will also be a spline function it can only be defined for a certain range of rho values as node points. Therefore it is better to limit the range of electron density values and define larger limits for the embedding function instead.

Finally the last node point at the cutoff range is set to 0 and fitting is disabled to prevent a discontinuous change of energy at the cutoff.

%% Cell type:code id:mexican-absence tags:

``` python

rho.parameters["node_7.0"].start_val=0

rho.parameters["node_7.0"].enabled=False

```

%% Cell type:markdown id:standard-relative tags:

$-\sqrt(\rho)$ can be used as initial guess for the embedding energy, which is taken from second moment approximation tight binding.

The node points have to be chosen in a range compatible to the electron density. This can be estimated by calculating it for a densely packed structure.

Alternatively atomicrex writes the maximum electron density of all structures to the output. This can be used as a hint for the node points for consequent fits.

Everything else is similar to the electron density.

Finally a few parameters need to be set that influence the fitting process.

The minimization can be done with different algorithms. Atomicrex itself implements the BFGS algorithm. Additionally the algorithms from the nlopt library can be used.

The job PotentialTest was saved and received the ID: 820

The job strain_0_9 was saved and received the ID: 821

The job strain_0_92 was saved and received the ID: 822

The job strain_0_94 was saved and received the ID: 823

The job strain_0_96 was saved and received the ID: 824

The job strain_0_98 was saved and received the ID: 825

The job strain_1_0 was saved and received the ID: 826

The job strain_1_02 was saved and received the ID: 827

The job strain_1_04 was saved and received the ID: 828

The job strain_1_06 was saved and received the ID: 829

The job strain_1_08 was saved and received the ID: 830

The job strain_1_1 was saved and received the ID: 831

job_id: 821 finished

job_id: 822 finished

job_id: 823 finished

job_id: 824 finished

job_id: 825 finished

job_id: 826 finished

job_id: 827 finished

job_id: 828 finished

job_id: 829 finished

job_id: 830 finished

job_id: 831 finished

%% Cell type:code id:vocal-heather tags:

``` python

murn.plot()

```

%%%% Output: display_data

%% Cell type:code id:attached-palestinian tags:

``` python

murn["output/equilibrium_bulk_modulus"]

```

%%%% Output: execute_result

140.0186998321462

%% Cell type:markdown id:injured-rainbow tags:

### Same cane be done for the 1000 structures dataset

The final parameters of the 100 structure fit can be used for the 1000 Structure fit. This can speed up the fitting process and often leads to better results, especially if the initially guessed values are far from the optimum.

In general it is a good idea to start with few structures and try around with different functions and initial parameters. This is much faster than using all structures from the beginning and gives good guesses for the initial values of the parameters. It also allows to use global optimization with millions of steps in short time spans.

### This can take long and writes 1000 seperate POSCAR files, TAKE CARE

Run these jobs with more cores or more time after the workshop. Also increase the number of iterations for better results

With this choice of functions and initial parameters starting directly from all structures gives the same residual. In a previous iteration of the potential it was about 7 times worse, so it is a good idea to test this.