Exploring the parameter landscape

In general the objective function \(\chi^2\) is a very complex function of the potential parameters \(p_i\) as well as the weight factors \(w_i\). The resulting space \((\{p_i\};\{w_i\})\) can be excruciatingly difficult to navigate. In this context some basic analysis of the parameter can yield useful insight.

Consider specifically the very simple embedded atom method (EAM) type potential for aluminum that was already considered in several examples: [1], [2], [3]. The embedded atom method (EAM) potential format used here is a special case of the Tersoff potential with a purely repulsive pair potential

\[V(r) = A \exp(-\lambda r),\]

the embedding function

\[F(\rho) = -D \exp(-\rho),\]

and the electron density

\[\rho(r) = \exp(-2\mu r).\]

The four parameters are fitted to 108 atom face-centered cubic (FCC) configurations at the equilibrium lattice parameter as well as 10% compression and expansion. Also included in the fit is the Al bulk modulus and the FCC lattice parameter.

The previous examples started from the initial parameter set

\[(A,\lambda,\mu,D) = (500, 2.73, 1.14, 8)\]

When using local optimizers, in particular L-BFGS (internal L-BFGS minimizer, several local SciPy minimizers), the final parameter set tends to end up rather close to this initial set. Especially the \(A\) parameter appears to be rather “immobile”. Considering the functional form it is apparent that \(A\) and \(\lambda\) are coupled and because of the exponential function for our purpose of parameter optimization are related in unfavorable fashion.

To quantify their correlation one can for example conduct a scan of the \((A,\lambda)\) parameter plane. This can be readily achieved using atomicrex’s Python interface. Assuming that input files with potential and structure definitions have been prepared (main.xml, potential.xml, structures.xml; see below) the job can be set up as follows:

1
2
3
4
5
6
7
from __future__ import print_function
import atomicrex
import numpy as np

job = atomicrex.Job()
job.parse_input_file('main.xml')
job.prepare_fitting()

After storing the initial parameter vector:

params0 = job.get_potential_parameters()

one can sweep over a range of parameters:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
del_A = 0.02 * params0[0]
del_Lambda = 0.02 * params0[1]
n = 15
data = []
for dA in np.arange(-n * del_A, n * del_A, del_A):
  row = []
  data.append(row)
  params[0] = params0[0] + dA
  for dLambda in np.arange(-n * del_Lambda, n * del_Lambda, del_Lambda):
    params[1] = params[0] + dLambda
    row.append(job.calculate_residual(params))

The resulting matrix can be visualized using e.g., matplotlib. This clearly illustrates the extreme imbalance between the \(A\) and \(\lambda\) parameters as small (relative) changes in \(\lambda\) have a much more dramatic effect than in \(A\). This renders the minimization of the objective cumbersome, even in this very simple example.

A more systematic way to quantify parameter correlation is achieved by computing the second derivative (Hessian) of the objective function \(\Delta \chi^2 = \partial^2 \chi^2/\partial p_i\partial p_j\):

1
2
3
4
5
hessian = job.calculate_hessian()
for row in hessian:
  for d in row:
    print('{:10g}'.format(d), end='')
  print('')

This reveals a very large spread of values in particular on the diagonal and the extreme ratio between the off-diagonal values already indicates this to be a poorly conditioned problem. This becomes more apparent by analyzing the ratio of eigenvalues:

eig, vecs = np.linalg.eig(hessian)
print(np.max(np.abs(eig)) / np.min(np.abs(eig)))

This analysis thus illustrates the origin of the numerical diffculties encountered when optimizing this particular functional form. Yet, combinations such as \(A\) and \(\lambda\) are frequent in commonly employed potential forms.

The present analysis suggests that an analysis of the eigenvalue spectrum of the Hessian of the objective function can be a useful tool for testing the behavior of the functional form, especially when encountering difficulties in exploring the (local) parameter space.

The following example is revisited in the section on Monte Carlo sampling of the objective function, which turns out to be an interesting approach to navigating “narrow” parameter spaces such as the one discussed above.

Location of files

The example presented above employs input files from this example, which are also included here for convenience.

Input files

  • main.xml: main input file

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
    <?xml version="1.0" encoding="iso-8859-1"?>
<job>
  <name>EAM potential</name>
  
  <verbosity>medium</verbosity>

  <fitting>
    <BFGS conv-threshold="1e-10" max-iter="200" />
  </fitting>

  <atom-types>
    <species mass='26.981' atomic-number='13'>Al</species>
  </atom-types>

  <potentials>
    <xi:include href="potential.xml"
		xmlns:xi="http://www.w3.org/2003/XInclude" />
  </potentials>
  
  <structures>
    <xi:include href="structures.xml"
		xmlns:xi="http://www.w3.org/2003/XInclude" />
  </structures>

</job>

  • potential.xml: initial parameter set (included in main input file via XML Inclusions)

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
    <eam id='Al potential' species-a='*' species-b='*'>
  <mapping>
    <pair-interaction species-a='*' species-b='*' function='V' />
    <electron-density species-a='*' species-b='*' function='rho' />
    <embedding-energy species='*' function='F' />
  </mapping>
  <export-eam-file>Al.eam.fs</export-eam-file>
  <functions>
    <user-function id='V'>
      <input-var>r</input-var>
      <expression>
        A * exp(-lambda * r)
      </expression>
      <derivative>
        -lambda * A * exp(-lambda * r)
      </derivative>
      <param name='A'>500</param>
      <param name='lambda'>2.73</param>
      <fit-dof>
        <A/>
        <lambda/>
      </fit-dof>
      <screening>
        <user-function id='V_screening'>
          <cutoff>6.5</cutoff>
          <input-var>r</input-var>
          <expression>
            1 - 1/(1 + ((r - cutoff) / h)^4)
          </expression>
          <derivative>
            4 * h^4 * (r - cutoff)^3 / ((h^4 + (r - cutoff)^4)^2)
          </derivative>
          <param name='h'>3</param>
        </user-function>
      </screening>
    </user-function>
    <user-function id='rho'>
      <input-var>r</input-var>
      <expression>
        exp(-twomu * r)
      </expression>
      <derivative>
        -twomu * exp(-twomu * r)
      </derivative>
      <param name='twomu'>1.14</param>
      <fit-dof>
        <twomu/>
      </fit-dof>      
      <screening>
        <user-function id='rho_screening'>
          <cutoff>6.5</cutoff>
          <input-var>r</input-var>
          <expression>
            1 - 1/(1 + ((r - cutoff) / h)^4)
          </expression>
          <derivative>
            4 * h^4 * (r - cutoff)^3 / ((h^4 + (r - cutoff)^4)^2)
          </derivative>
          <param name='h'>3</param>
        </user-function>
      </screening>
    </user-function>
    <user-function id='F'>
      <input-var>rho</input-var>
      <expression>
        -D * sqrt(rho)
      </expression>
      <derivative>
        -D / (2 * sqrt(rho))
      </derivative>
      <param name='D'>8</param>
      <fit-dof>
        <D/>
      </fit-dof>            
    </user-function>
  </functions>
</eam>

  • structures.xml: input structures (included in main input file via XML Inclusions)

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
    <group>

  <fcc-lattice id='FCC' relative-weight='3'>
    <atom-type>Al</atom-type>
    <lattice-parameter>4.0</lattice-parameter>
    <properties>
      <lattice-parameter fit='true' target='4.032' relative-weight='3'/>
      <atomic-energy fit='true' target='-3.36' relative-weight='3'/>
      <bulk-modulus target='80.9' relative-weight='1'/>
      <C11/>
      <C12/>
      <C44/>
    </properties>
  </fcc-lattice>
  
  <user-structure id='structure_1'>
    <poscar-file>POSCAR_res_POSCAR_0.9.traj_forces</poscar-file>
    <properties>
      <atomic-forces fit='true' output-all='true'/>
    </properties>
  </user-structure>
  
  <user-structure id='structure_2'>
    <poscar-file>POSCAR_res_POSCAR_1.0.traj_forces</poscar-file>
    <properties>
      <atomic-forces fit='true' output-all='true'/>
    </properties>
  </user-structure>
  
  <user-structure id='structure_3'>
    <poscar-file>POSCAR_res_POSCAR_1.1.traj_forces</poscar-file>
    <properties>
      <atomic-forces fit='true' output-all='true'/>
    </properties>
  </user-structure>
  
</group>