Example: Basic fitting

This example illustrates the optimization of potential parameter via the Python interface. It demonstrates both how to start optimizers that are integrated directly in atomicrex and how to use optimizers imported from other Python librariers, specifically scipy. Here, potential definition and structures are provided via xml files. Note that in the case of the internal optimizers the settings are read from the xml input file.

Location

examples/python_basic_fitting

Input files

  • fit_potential.py: 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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    import atomicrex
from scipy.optimize import minimize
import argparse

__doc__ = """ This script demonstrates the application of the python interface
for parameter optimization.  """

parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('mode', choices=['atomicrex', 'scipy'],
                    help='Use either the internal or a scipy optimizer.')
args = parser.parse_args()

# initialize job
job = atomicrex.Job()
job.parse_input_file('main.xml')
job.set_verbosity(2)

# fitting
job.prepare_fitting()
print('initial residual: %g' % job.calculate_residual())
if args.mode == 'atomicrex':
    job.perform_fitting()
elif args.mode == 'scipy':
    pass
    initial_params = job.get_potential_parameters()
    resopt = minimize(job.calculate_residual,
                      initial_params,
                      options={'disp': False,
                               'eps': 1e-4,
                               'gtol': 1e-10},
                      method='L-BFGS-B')
residual = job.calculate_residual()
print('final residual: %g' % job.calculate_residual())

# output properties
for name, structure in sorted(job.structures.items()):
    print('')
    print(' Structure: %s' % name)
    structure.compute_properties()
    structure.print_properties()

# parameters
job.print_potential_parameters()

  • 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
    <?xml version='1.0' encoding='iso-8859-1'?>
<job>
  <name>Al EAM potential</name>

  <fitting>
    <BFGS conv-threshold='1e-10' max-iter='200' />
  </fitting>
  
  <atom-types>
    <species>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: definition of potential form and parameters

     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
78
79
80
81
82
83
    <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>
  
  <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='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='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: definition of structures and properties

     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
    <group>

  <fcc-lattice id='fcc' relative-weight='5'>
    <atom-type>Al</atom-type>
    <lattice-parameter>4.05</lattice-parameter>
    <properties>
      <lattice-parameter relative-weight='2' fit='true' min='3' max='5' target='4.05' />
      <atomic-energy relative-weight='2' fit='true' target='-3.36' />
      <bulk-modulus fit='true' target='79' relative-weight='1'/>
      <C11 fit='true' target='114' relative-weight='1'/>
      <C12 fit='true' target='61.9' relative-weight='1'/>
      <C44 fit='true' target='31.6' relative-weight='1'/>
    </properties>
  </fcc-lattice>
  
  <bcc-lattice id='bcc'>
    <atom-type>Al</atom-type>
    <lattice-parameter>2.6</lattice-parameter>
    <properties>
      <lattice-parameter min='2.0' max='3.3' fit='true' target='4.0'/>
      <atomic-energy relative-weight='0.1' fit='true' target='-3.25' />
    </properties>
  </bcc-lattice>
  
  <hcp-lattice id='hcp'>
    <atom-type>Al</atom-type>
    <lattice-parameter>3.0</lattice-parameter>
    <ca-ratio>1.6</ca-ratio>
    <properties>
      <lattice-parameter min='2.0' max='3.3' />
      <atomic-energy relative-weight='0.1' fit='true' target='-3.33' />
    </properties>
  </hcp-lattice>
  
  <diamond-lattice id='diamond'>
    <atom-type>Al</atom-type>
    <lattice-parameter>5.0</lattice-parameter>
    <properties>
      <lattice-parameter min='4.0' max='6.0' />
      <atomic-energy relative-weight='0.01' fit='true' target='-2.36' />
    </properties>
  </diamond-lattice>

</group>

Output (files)

  • log.python_atomicrex: stdout from run using an internal optimizer from atomicrex

      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
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    initial residual: 1480.32
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.73
   Al[EAM].rho[user-function].twomu: 1.14
   Al[EAM].F[user-function].D: 8
-------------------------------------------------------
-------------------------------------------------------
Number of properties being fitted: 9
-------------------------------------------------------
  iter= 1 residual= 1480.32 grad= 20745.7
  iter= 2 residual= 203907 grad= 746759
  iter= 3 residual= 1159.03 grad= 13279.2
  iter= 4 residual= 1159.03 grad= 13279.2
  iter= 5 residual= 1119.64 grad= 6617.21
  iter= 6 residual= 1119.64 grad= 6617.21
  iter= 7 residual= 1103.83 grad= 5969.37
  iter= 8 residual= 1103.83 grad= 5969.37
  iter= 9 residual= 1099.92 grad= 14765.1
  iter= 10 residual= 1099.92 grad= 14765.1
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.78517
   Al[EAM].rho[user-function].twomu: 1.19806
   Al[EAM].F[user-function].D: 7.99568
-------------------------------------------------------
  iter= 11 residual= 1035.19 grad= 3109.96
  iter= 12 residual= 1035.19 grad= 3109.96
  iter= 13 residual= 1257 grad= 9397.81
  iter= 14 residual= 1041.8 grad= 3296.9
  iter= 15 residual= 1032.74 grad= 2872.63
  iter= 16 residual= 1035.88 grad= 5237.62
  iter= 17 residual= 1032.48 grad= 7413.92
  iter= 18 residual= 1032.48 grad= 7413.92
  iter= 19 residual= 1047.64 grad= 440.704
  iter= 20 residual= 1029.1 grad= 5825.19
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.8028
   Al[EAM].rho[user-function].twomu: 1.22477
   Al[EAM].F[user-function].D: 7.99466
-------------------------------------------------------
  iter= 21 residual= 1029.1 grad= 5825.19
  iter= 22 residual= 1040.72 grad= 6364.4
  iter= 23 residual= 1028.84 grad= 3684.53
  iter= 24 residual= 1028.84 grad= 3684.53
  iter= 25 residual= 1132.41 grad= 11573.9
  iter= 26 residual= 1031.78 grad= 4488.16
  iter= 27 residual= 1029.16 grad= 3640.6
  iter= 28 residual= 1028.89 grad= 3555.82
  iter= 29 residual= 1028.85 grad= 3811.61
  iter= 30 residual= 1028.85 grad= 3636.03
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.8047
   Al[EAM].rho[user-function].twomu: 1.22771
   Al[EAM].F[user-function].D: 7.99476
-------------------------------------------------------
  iter= 31 residual= 1028.85 grad= 4380.48
  iter= 32 residual= 1028.85 grad= 4005
  iter= 33 residual= 1028.85 grad= 3944.25
  iter= 34 residual= 1028.85 grad= 3978.34
  iter= 35 residual= 1028.85 grad= 3632.17
  iter= 36 residual= 1028.85 grad= 4053.33
  iter= 37 residual= 1028.85 grad= 4053.33
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.8047
   Al[EAM].rho[user-function].twomu: 1.22771
   Al[EAM].F[user-function].D: 7.99476
-------------------------------------------------------
Fitting process complete.
final residual: 1028.85
()
 Structure: bcc
total-energy: -2.86993 eV
atomic-energy: -2.86993 eV/atom (target: -3.25 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 180.567
total-volume: 8.788 A^3
atomic-volume: 8.788 A^3/atom
lattice-parameter: 2.6 A [2:3.3]
()
 Structure: diamond
total-energy: -5.21398 eV
atomic-energy: -2.60699 eV/atom (target: -2.36 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 76.256
total-volume: 31.25 A^3
atomic-volume: 15.625 A^3/atom
lattice-parameter: 5 A [4:6]
()
 Structure: fcc
total-energy: -3.58951 eV
atomic-energy: -3.58951 eV/atom (target: -3.36 eV/atom) Abs. weight: 0.15625; Tolerance: 0.01; Weighted residual: 82.3066
total-volume: 12.9217 A^3
atomic-volume: 12.9217 A^3/atom
bulk-modulus: 87.5623 GPa (target: 79 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 1.43191
C11: 100.857 GPa (target: 114 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 3.37394
C12: 80.9163 GPa (target: 61.9 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 7.06291
C44: 45.186 GPa (target: 31.6 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 3.60509
lattice-parameter: 3.72501 A (target: 4.05 A) Abs. weight: 0.15625; Tolerance: 0.005; Weighted residual: 660.134 (relaxed) [3:5]
()
 Structure: hcp
total-energy: -6.44751 eV
atomic-energy: -3.22375 eV/atom (target: -3.33 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 14.1104
total-volume: 37.4123 A^3
atomic-volume: 18.7061 A^3/atom
lattice-parameter: 3 A [2:3.3]
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.8047
   Al[EAM].rho[user-function].twomu: 1.22771
   Al[EAM].F[user-function].D: 7.99476
-------------------------------------------------------

  • log.python_scipy: stdout from run using optimize.minimize from scipy

     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
    initial residual: 1480.32
final residual: 1002.29
()
 Structure: bcc
total-energy: -2.77015 eV
atomic-energy: -2.77015 eV/atom (target: -3.25 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 287.814
total-volume: 8.788 A^3
atomic-volume: 8.788 A^3/atom
lattice-parameter: 2.6 A [2:3.3]
()
 Structure: diamond
total-energy: -5.16885 eV
atomic-energy: -2.58442 eV/atom (target: -2.36 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 62.957
total-volume: 31.25 A^3
atomic-volume: 15.625 A^3/atom
lattice-parameter: 5 A [4:6]
()
 Structure: fcc
total-energy: -3.60106 eV
atomic-energy: -3.60106 eV/atom (target: -3.36 eV/atom) Abs. weight: 0.15625; Tolerance: 0.01; Weighted residual: 90.7944
total-volume: 13.2414 A^3
atomic-volume: 13.2414 A^3/atom
bulk-modulus: 85.9308 GPa (target: 79 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 0.938209
C11: 98.8996 GPa (target: 114 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 4.45358
C12: 79.4466 GPa (target: 61.9 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 6.01336
C44: 44.4446 GPa (target: 31.6 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 3.22232
lattice-parameter: 3.75547 A (target: 4.05 A) Abs. weight: 0.15625; Tolerance: 0.005; Weighted residual: 542.169 (relaxed) [3:5]
()
 Structure: hcp
total-energy: -6.54788 eV
atomic-energy: -3.27394 eV/atom (target: -3.33 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 3.92864
total-volume: 37.4123 A^3
atomic-volume: 18.7061 A^3/atom
lattice-parameter: 3 A [2:3.3]
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.77894
   Al[EAM].rho[user-function].twomu: 1.21223
   Al[EAM].F[user-function].D: 7.9934
-------------------------------------------------------

  • log.no_python [for reference]: stdout from “pure” (no Python involved) execution of atomicrex

      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
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
    This is program version 0.1.6
Reading job file main.xml
-------------------------------------------------------
Parsing input file(s)
-------------------------------------------------------
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.73
   Al[EAM].rho[user-function].twomu: 1.14
   Al[EAM].F[user-function].D: 8
-------------------------------------------------------
-------------------------------------------------------
Number of properties being fitted: 9
-------------------------------------------------------
  iter= 1 residual= 1480.32 grad= 20791.4
  iter= 2 residual= 210872 grad= 778758
  iter= 3 residual= 1160.52 grad= 13390.1
  iter= 4 residual= 1160.52 grad= 13390.1
  iter= 5 residual= 1121.04 grad= 6634.79
  iter= 6 residual= 1121.04 grad= 6634.79
  iter= 7 residual= 1104.36 grad= 5922.8
  iter= 8 residual= 1104.36 grad= 5922.8
  iter= 9 residual= 1105.91 grad= 14901.1
  iter= 10 residual= 1085.19 grad= 6171.57
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.75779
   Al[EAM].rho[user-function].twomu: 1.18125
   Al[EAM].F[user-function].D: 7.99642
-------------------------------------------------------
  iter= 11 residual= 1085.19 grad= 6171.57
  iter= 12 residual= 1089.99 grad= 11320.4
  iter= 13 residual= 1061.54 grad= 5547.1
  iter= 14 residual= 1061.54 grad= 5547.1
  iter= 15 residual= 1686.51 grad= 8513.32
  iter= 16 residual= 1086.35 grad= 3507.56
  iter= 17 residual= 1058.83 grad= 5150.04
  iter= 18 residual= 1058.83 grad= 5150.04
  iter= 19 residual= 1116.08 grad= 7584.27
  iter= 20 residual= 1054.8 grad= 3933.06
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.79747
   Al[EAM].rho[user-function].twomu: 1.21346
   Al[EAM].F[user-function].D: 7.99374
-------------------------------------------------------
  iter= 21 residual= 1054.8 grad= 3933.06
  iter= 22 residual= 1278.59 grad= 2799.79
  iter= 23 residual= 1063.77 grad= 2936.99
  iter= 24 residual= 1055.87 grad= 3544.03
  iter= 25 residual= 1054.98 grad= 4030.37
  iter= 26 residual= 1054.83 grad= 3909.24
  iter= 27 residual= 1054.81 grad= 3651.84
  iter= 28 residual= 1054.8 grad= 3683.48
  iter= 29 residual= 1054.81 grad= 9896.22
  iter= 30 residual= 1054.81 grad= 9763.95
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.79747
   Al[EAM].rho[user-function].twomu: 1.21346
   Al[EAM].F[user-function].D: 7.99374
-------------------------------------------------------
  iter= 31 residual= 1054.81 grad= 9642.46
  iter= 32 residual= 1054.81 grad= 9248.07
  iter= 33 residual= 1054.81 grad= 9641.22
  iter= 34 residual= 1054.8 grad= 9445.2
  iter= 35 residual= 1054.8 grad= 9445.2
-------------------------------------------------------
Potential parameters being optimized (#dof=4):
   Al[EAM].V[user-function].A: 500
   Al[EAM].V[user-function].lambda: 2.79747
   Al[EAM].rho[user-function].twomu: 1.21346
   Al[EAM].F[user-function].D: 7.99374
-------------------------------------------------------
Fitting process complete.
-------------------------------------------------------
Computing structure properties
Structure 'fcc':
   total-energy: -3.6615 eV
   atomic-energy: -3.6615 eV/atom (target: -3.36 eV/atom) Abs. weight: 0.15625; Tolerance: 0.01; Weighted residual: 142.035
   total-volume: 12.9435 A^3
   atomic-volume: 12.9435 A^3/atom
   bulk-modulus: 88.415 GPa (target: 79 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 1.7313
   C11: 101.824 GPa (target: 114 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 2.8957
   C12: 81.712 GPa (target: 61.9 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 7.6663
   C44: 45.7691 GPa (target: 31.6 GPa) Abs. weight: 0.078125; Tolerance: 2; Weighted residual: 3.92113
   lattice-parameter: 3.72709 A (target: 4.05 A) Abs. weight: 0.15625; Tolerance: 0.005; Weighted residual: 651.68 (relaxed) [3:5]
Structure 'bcc':
   total-energy: -2.92726 eV
   atomic-energy: -2.92726 eV/atom (target: -3.25 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 130.203
   total-volume: 8.788 A^3
   atomic-volume: 8.788 A^3/atom
   lattice-parameter: 2.6 A [2:3.3]
Structure 'hcp':
   total-energy: -6.58809 eV
   atomic-energy: -3.29404 eV/atom (target: -3.33 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 1.61616
   total-volume: 37.4123 A^3
   atomic-volume: 18.7061 A^3/atom
   lattice-parameter: 3 A [2:3.3]
Structure 'diamond':
   total-energy: -5.32148 eV
   atomic-energy: -2.66074 eV/atom (target: -2.36 eV/atom) Abs. weight: 0.125; Tolerance: 0.01; Weighted residual: 113.056
   total-volume: 31.25 A^3
   atomic-volume: 15.625 A^3/atom
   lattice-parameter: 5 A [4:6]
-------------------------------------------------------
-------------------------------------------------------