Post on 10-Oct-2020
Representing High-Dimensional Potential-Energy Surfaces
by Artificial Neural Networks
Jörg Behler
Lehrstuhl für Theoretische ChemieRuhr-Universität Bochum
D-44780 Bochum, Germany
joerg.behler@theochem.ruhr-uni-bochum.de
Energy Landscapes, Chemnitz 2010
Introduction: Simulation of Complex Condensed Systems
⇒ Reliable potentials for large-scale simulations are needed!
Surface Science:
• crystal structure prediction• phase diagrams• properties of materials
Challenges:
• heterogeneous catalysis• corrosion• self-assembled monolayers
Materials Science:
• complicated bonding situations• unusual coordination• different types of bonding:covalent, metallic, vdW etc.
• large systems• long simulation times
Potentials for Atomistic Simulations
The reliability of the results obtained in theoretical simulations depends on the accuracy of the underlying potentials.
accu
racy
spee
d
atoms time
≈ 10 0
≈ 100 100 ps
≈ 1000 1 ns
≈ 10,000 10 ns
≈ 100,000 100 ns
ab initio (MP2, CC, CI)
classical force fields
Hartree Fock /density functional theory
semiempirical / tight binding
empirical potentials
Potentials for Atomistic Simulations
“List of wishes”
• should be very accurate
• should describe making and breaking of bonds
• should be applicable to all types of systems (metals and insulators, solids and molecules)
• should be fast to evaluate
• systematic improvement should be possible
• no “manual” work for construction
• should be transferable (“predictive”)
Artificial Neural Networks
y
x
x
x
1
2
3
Σ
dendritesoma
axon
axon hillock terminal button
W. McCulloch, and W. Pitts, Bull. Math. Biophys. 5, 115 (1943).
Signal Processing in “Bio-Neurons”
Artificial Neuron
Artificial Neural Networks
InputLayer
HiddenLayer
OutputLayer
T.B. Blank, S.D. Brown, A.W. Calhoun, and D.J. Doren, J. Chem. Phys. 103 (1995) 4129.S. Lorenz, A. Groß, and M. Scheffler, Chem. Phys. Lett. 395 (2004) 210.J. Behler, S. Lorenz, and K. Reuter, J. Chem. Phys. 127 (2007) 014705.
G1
G2
y21
y31
y11
Bias
E
w111
w231
w112
w212
w312
w011
w012
w031
⎟⎟⎠
⎞⎜⎜⎝
⎛⎟⎠
⎞⎜⎝
⎛+⋅+= ∑∑
==
2
1
110
13
1
21
201
2
iiijja
jja GwwfwwfE
Analytic Total Energy Expression:
Properties of a NN Potential:• very flexible functional form• very accurate• no knowledge of functional form required• properties of the system can be included
via mapping on symmetry functions• bad for extrapolation• many training points needed• analytic forces
Example:
Activation Functions
Activation functions enable the fitting of general nonlinear functions.
hyperbolic tangent sigmoid function
Examples:
Activation functions• converge for very small and very large arguments• have a nonlinear shape for intermediate values
Flexibility of the Activation Functions
⇒ very flexible but simple
( ) ( ) dbaxcxf ++= tanhBasic functional element of the NN:
Basic Example: Harmonic Oscillator
Fitting of the potential of a 1D harmonic oscillator for -3 < x < 3:
⇒ Just two functions provide a rather good approximation⇒ Fit is further improved by adding more functions
Basic Example: Dimer Potential
⇒ NN has a non-physical functional form, but can learn physics
r (a.u.)
E(a
rb. u
nits
)
More realistic example: Pair potential
Crystal Structure Prediction
Find the relevant structures
Determine the stability of many structures
• experiment• chemical intuition• molecular dynamics (e.g. Parrinello Rahman)• genetic algorithms• metadynamics• …
• fast empirical potentials for first evaluation ⇒ often unreliable• DFT ⇒ too expensive
Two Challenges:
⇒ We employ a Neural Network potential
Limitations of Standard Feed-Forward Neural Networks
Conceptual problems for a direct application of Neural Networks to condensed systems
• limited number of dimensions (up to 12)
• symmetry of the system is not included(exchange of atoms changes the energy output)
• energy possibly depends on rotation and translation
• fit is valid only for a given system size
⇒ A new Neural Network scheme is required to deal with high-dimensional systems
Is it possible to use Neural Network potentials to construct high-dimensional potential energy surfaces?
G1
G2
y21
y31
y11
Bias
E
w111
w231
w112
w212
w312
w011
w012
w031
The High-Pressure Crystal Structures of Silicon
cubic diamond β-tin
Review: A. Mujica, A. Rubio, A. Muñoz, and R.J. Needs, Rev. Mod. Phys. 75 (2003) 863.
Imma simple hexagonal
Cmcahcpfcc
11.7 GPa 13.2 GPa 15.4 GPa
38 GPa42 GPa79 GPa
⇒ Challenging model system to study pressure-induced phase transitions
DFT Calculations
Code: PWSCFCode: PWSCF
• plane waves• plane waves• plane waves• plane waves• plane waves• plane waves• plane waves• pseudo potentials• pseudo potentials• pseudo potentials• pseudo potentials• pseudo potentials• pseudo potentials• pseudo potentials• stress tensor for k-points• stress tensor for
Calculational Details:Calculational Details:Calculational Details:
• 20 Ry cutoff• 3x3x3 k-points
Timing: 4-10 h on a 2.8 GHz Opteron (serial run)
System Details:System Details:System Details:System Details:System Details:System Details:System Details:System Details:
• 64 Si atoms• 64 Si atoms• Pressure range 0 – 100 GPa, Temperature 0 – 3000 K• Pressure range 0 – 100 GPa, Temperature 0 – 3000 K• Pressure range 0 – 100 GPa, Temperature 0 – 3000 K
PWSCF:S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org/.
k-points• stress tensor for k-points
Calculational Details:
PWSCF:S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org/.
Code: PWSCF
• stress tensor for k-points
Calculational Details:
PWSCF:S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org/.
Code: PWSCF
• stress tensor for k-points
Calculational Details:
PWSCF:S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org/.
Code: PWSCF
• stress tensor for k-points
Calculational Details:
Timing: 4-10 h on a 2.8 GHz Opteron (serial run)
• 64 Si atoms
PWSCF:S. Baroni, A. Dal Corso, S. de Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org/.
Code: PWSCF
• stress tensor for k-points
Calculational Details:
Metadynamics: CPU Time Estimate
A metadynamics run takes about 100 metasteps, and each stepincludes a MD run of about 2 ps.
⇒ 200 000 energy and force evaluations per metadynamics run⇒ 1 600 000 CPU hours = 67 000 CPU days = 183 CPU years⇒ We cannot use DFT directly!
We need many metadynamics runs!
⇒ We need a significant speed-up in the potential evaluation, butwithout a significant loss in accuracy
Can we combine metadynamics with neural networks?
Neural Network for High-Dimensional Systems
Idea:
Represent the total energy as a sum of atomic energies:
∑=i
iEE
⇒ requires a new Neural Network architecture and suitablesymmetry functions
J. Behler, and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
Example: fcc structure (4 atom unit cell)
= + + +
Each Ei is given by an individual NN as a function of the chemicalenvironment
Neural Network for High-Dimensional Systems
AtomicEnergies
S1
S2
S3
E1
E2
E3
EΣ
Total Energy
AtomicSubnets
Ri={Xi,Yi,Zi}
R1
R2
R3
AtomicPositions
∑=i
iEE
G1
SymmetryFunctions
G2
G3
Input Invisible Output
Example: 3-atom system
Structuresidentical
J. Behler, and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).J. Behler, R. Martoňák, D. Donadio, and M. Parrinello, phys. stat. sol. (b) 245, 2618 (2008).
„Structuralfingerprint“
Neural Network Accuracy: Total Energies
-6925
-6920
-6915
-6910-6925 -6920 -6915 -6910
training settest set
-510
-509
-508
-507
-506
-505
-504
-503
-502
-501
-500-510 -508 -506 -504 -502 -500
training settest set
EDFT (Ry/64 atoms)
Training setTest set
Points:171441907
RMSE (meV/atom):56
MAD (meV/atom):45
Accuracy
cubicdiamond
EDFT (eV/64 atoms)
EN
N(eV)
EN
N(R
y)
Neural Network Accuracy: Total Energies
Energy error distribution:
Neural Network Accuracy: Crystal Structures of Silicon
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
60 70 80 90 100 110 120 130 140 150
volume per atom (Bohr3)
ener
gype
r ato
m(e
V)
sc DFTsc NNdiamond DFTdiamond NNβ-tin DFTβ-tin NNsh DFTsh NNhcp DFThcp NNfcc DFTfcc NNbcc DFTbcc NNbc8 DFTbc8 NNhex. diamond DFThex. diamond NN
Ideal Crystal Structures
Neural Network Accuracy: Radial Distribution Functions Liquid Silicon
0.0
0.5
1.0
1.5
2.0
2.5
3.0
1.5 2.5 3.5 4.5 5.5
DFTNNBazantLenoskyTersoff
r (Ǻ)
Liquid Si, 64 atom cell, 3000 K, a = 10.862 Ǻ
DFT data: T. Kühne, priv. comm.
Disordered Structures:
Neural Network: Computational Costs
64 atom silicon cell
DFT (scf): about 5-10 hours (k-points!)NN: ≈ 1 second (depending on symmetry functions)
Timing for the calculation of energy, forces and stress (single Opteron CPU) :
⇒ Speed up factor: 20000
⇒ Now 100 ps MD on a single CPU per day with about DFT accuracy
⇒ Time for a metadynamics run is reduced from 200 CPU years toabout 2 days on a standard PC.
Neural Network: Performance
0
200
400
600
800
1000
1200
0 100 200 300 400 500 600
number of atoms
t(ar
b. u
nits
)
Benchmark: MD Simulation
Linear Scaling with System Size: Parallelization:
⇒ Neural Network is about 20000 – 100000 times faster than DFT
number of cores
Metadynamics Simulations of Phase Transitions in Silicon
diamond β-Sn Imma sh Cmca hcp fcc
T = 300 K, p = 12 GPa
(a)
(b)
(c)
(d)
Enthalpy barrier ≈ 0.1 eV/atom
J. Behler, R. Martoňák, D. Donadio, and M. Parrinello, Phys. Rev. Lett. 100, 185501 (2008).
Metadynamics Simulations of Phase Transitions in Silicon
diamond β-Sn Imma sh Cmca hcp fcc
T = 300 K, p = 16 GPa
(a) (b)
(c)
Metadynamics Simulations of Phase Transitions in Silicon
diamond β-Sn Imma sh Cmca hcp fcc
T = 800 K, p = 45 GPa
(a) (b)
J. Behler, R. Martoňàk, D. Donadio, and M. Parrinello, phys. stat. sol. (b) 245, 2618 (2008).
Metadynamics Simulations: Check of the Potential
Recalculation of the final structure of each metastep with DFT:
diamond β-Sn Imma sh Cmca hcp fcc
T = 800 K, p = 45 GPa
-107,70
-107,60
-107,50
-107,40
-107,30
-107,20
-107,100 20 40 60 80 100
DFT finalNN final
metastep
E per atom(eV)
DFTNN
Metadynamics Simulations of Phase Transitions in Silicon
diamond β-Sn Imma sh Cmca hcp fccT = 300 K, p = 50 GPa
(a)
(b)
(d)
(c)
Metadynamics Simulations of Phase Transitions in Silicon
diamond β-Sn Imma sh Cmca hcp fccT = 300 K, p = 90 GPa
(a)
(b)
Construction of the Training Data Set
Self-consistency:NN Fit ready for
application
Start withrandom structures
DFT calculations
Neural Network Fit
Apply fit to generatestructures (e.g. MD)
Identify poorlyrepresented structures
Qualityok
Neural Network Potential for Sodium
Training Set: 16500 0.72Test Set: 1500 0.91
RMSE (meV/atom):Structures:
DFT Code: PWSC, PBE functional PWSCF: S. Baroni et al., http://www.pwscf.org/.
a0 (Å) 4.201 4.202 5.297 5.297
DFT NN
bcc fcc
DFT NN
H. Eshet, R.Z. Khaliullin, T.D. Kühne, J. Behler, and M. Parrinello, Phys. Rev. B, in press (2010).
B0 (GPa) 7.63 7.59 7.624 7.613
Neural Network Potential for Sodium
H. Eshet, R.Z. Khaliullin, T.D. Kühne, J. Behler, and M. Parrinello, Phys. Rev. B, in press (2010).
NPT simulation of liquid sodium at 600 K, 100 GPa (64 atom cell)
Summary
Outlook:• useful tool to speed up extended simulations
Neural Network Potentials:• now applicable to condensed systems (solids, large clusters, liquids)
• NN reproduces total energies very accurately (also in value!)• provide analytic derivatives (forces, stress)• fast and linear scaling with system size• can be constructed using any electronic structure method
Limitations of NN Potentials:• need many reference calculations• limited transferability, no extrapolation• limited number of chemical elements