University of Southern California Dissertations and Theses

Physics informed neural networks and electrostrictive cavitation in water

(USC Thesis Other)

Physics informed neural networks and electrostrictive cavitation in water

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Physics Informed Neural Networks and Electrostrictive Cavitation in Water
by
Shane Jackson
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Physics)
August 2020
Copyright 2020 Shane Jackson
Dedication
This dissertation is dedicated to my father Eric Jackson and my mother Deborah Jackson.
ii
Acknowledgements
If I were to be fully honest, there are countless people I should acknowledge here. So many
friends, family members and even strangers have contributed through small acts of kindness
which, perhapsunbeknownsttothem, havehelpedmenishtheworkthathaswentintothis
thesis. For brevity sake, I shall lie by omission and thank only those with obvious impact.
IwouldliketothankmyadvisorRajivKaliaforhisyearsofadvice,patience,enlightening
conversations and moral support. Without a doubt, he has been a dutiful and fantastic
advisor.
I would like to thank my committee members: Aiichiro Nakano, Priya Vashishta, James
BoedickerandMarkusWintererfortheirsupportandtime. InmytimeinCACCS,Professors
Nakano and Vashishta have periodically offered me useful advice, guidance and knowledge.
ProfessorBoedickerallowedmetodomyrstresearchatUSCforwhichIwillbeperpetually
thankful. Furthermore, he has occasionally performed as the most valuable player (M.V.P)
for my trivia team. Professor Winterer has gone to extraordinary lengths answering all of
my, perhaps inane, questions on X-ray absorption spectroscopy. He was generous enough to
allow me to join his students in experiments at Argonne National Laboratory. Furthermore
I'd like to specically thank the students of Markus Winterer: Viktor Mackert and Thomas
Winter for helping me understand the XAFS experiments and helping in analyzing the data.
iii
I would like to thank my friends Allen Blaylock, Sara McCarthy and Patrick Edwards
for always being available to discuss my research and offering useful insight despite the
conversations often being outside their area of expertise. Also, they all tend to be decent
drinking companions.
I would like to thank Debbie Jackson and John Saggese for their love and support.
Although my father is not here to see me nish my thesis, I doubt I would have started
down this path if not for him. Last but most denitely not least, I would like to thank my
ance, Ana Hollander (t.b.Jackson) for her patience and love. I love you and I look forward
to our wedding.
iv
Table of Contents
Dedication ii
Acknowledgements iii
List Of Tables vii
List Of Figures viii
Abstract xiii
Chapter 1: Introduction 1
Variation Principles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Chapter2: ElectrostrictiveCavitationinWaterInducedbyaSnO
2
Nanopar-
ticle 12
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
Computational Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Chapter 3: XAS study of colloidal SnO
2
27
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Electrochemical double layer . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
X-ray Absorption Fine Structure (XAFS) . . . . . . . . . . . . . . . . . . . . 31
Reverse Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Chapter 4: A Neural Network for Principle of Least Action 48
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
v
Chapter 5: Solving IVPs with Least Action Neural Networks 62
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Chapter 6: Conclusion 70
Reference List 73
Appendix A
Cavitation Supporting Information . . . . . . . . . . . . . . . . . . . . . . . . . . 82
Force Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
Water without Nanoparticle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Void volume measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Video . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Appendix B
NNMD S.I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Lennard-Jones Sytstem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Pre-training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Results for 108 and 256-atom systems . . . . . . . . . . . . . . . . . . . . . . . . 90
vi
List Of Tables
A.1 Buckingham: E
ij
= B
i
jexp(r
ij
=
ij
)C
ij
=r
6
ij
taken from [7] . . . . . . . . 82
A.2 Lennard-Jones: E
ij
= A
ij
=r
12
ij
C
ij
=r
6
ij
taken from [7] . . . . . . . . . . . . . 82
A.3 AnisotropicHydrogen-bondingterm: E
HB
ij
= ϵ
HB
[
5
(

HB
r
ij
)
12
6
(

HB
r
ij
)
10
]
cos
4

The angle in the hydrogen-bonding interaction is the angle between the
oxygen atom accepting the hydrogen bond and the oxygen atom donating the
hydrogen bond. Taken from [40]. . . . . . . . . . . . . . . . . . . . . . . . . 82
A.4 Atom charges taken from [7] and [40] M and D refer to the dummy atom and
Drud charge respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
B.1 Lennard Jones interaction parameters and lattice constant. . . . . . . . . . . 89
vii
List Of Figures
1.1 Example of periodic boundary conditions. The simulation cell's (green) par-
ticle (red and blue) moves (top box) from one boundary to the opposite (bot-
tom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1 Initialmagnitudeandstructureofelectriceld. (a)Theaverage xcomponent
oftheelectriceldversustimefordifferentappliedelectriceldstrengths. (b)
The time it takes the electric eld to reach 95% of the saturation value. (c)
The magnitude of the saturation value for the electric eld versus applied
electric eld. (d) The x dependence of the x component of the electric eld at
different time steps correspond to the initial application of the eld, the rise
and the saturation for an applied eld of 0.125 V/

A. . . . . . . . . . . . . . 15
2.2 (a) Snapshot of the system just after cavitation. (b)

E vs t: Blue markers
are the data used to calculate the exponential relaxation coefficient . The
temporal data spans an order of magnitude (60 to 635 ps). The red line is the
tted line of best t. (c) The blue line is the plotted stretched exponential
function with = 0:424 and = 134:6 ps. . . . . . . . . . . . . . . . . . . . 16
2.3 Evolution of cavity with an applied 0.42 V/A electric eld. Snapshots at T
= 700 (a), 800 (b), 900 (c), 1000 (d), 1050 (e), 1100 (f) ps. Part (g) shows
the cavity volume as a function of time with red Xs over the data points
corresponding to parts a-f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Cavity volume analysis. The cavity volume as a function of time for different
applied electric elds. In many of the curves, two distinct regions of near
linear growth appear. The rst corresponds to the spherical growth phase
while the second corresponds to the cylindrical growth phase. . . . . . . . . 22
2.5 The x component of the electric eld at 10 ps (a, b, c) 40 ps(d, e, f) and
100 ps(g, h, i) for an applied eld of 0.125 V/

A. (a,d,g) Heat map of E
x
the
averagedoverthexaxis. (b,e,h)PlotofE
x
vsxpositionaveragedovery and
z axis. The red lines show the locations of the edge of the nanoparticle. (c,
f, i) Snapshot of the cavity and nanoparticle at timestep under consideration.
The electric eld spikes in the region around the cavity while decaying as
distance is increased from the cavity. . . . . . . . . . . . . . . . . . . . . . . 23
viii
3.1 Schematic diagram of a Electrode-electrolyte interface according to GCS the-
ory. The positive charge builds on the surface of the working electrode or, in
our case, the nanoparticle (Green) while the Stern layer is formed by negative
ions in the solution (Blue). Between the IHP and OHP is the solvated atom
layer. This is followed by the diffuse atom layer of the EDL and nally the
bulk solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 An incoming beam of X-rays passes through a material of thickness d. The
intensity is then measured by the detector. . . . . . . . . . . . . . . . . . . . 31
3.3 SchematicofX-rayabsorptioncoefficientversusincidentphotonenergyforan
arbitrarymaterial. Thetrendofdecreasing withincreasingEisshown. The
edges are labeled (K,L) and the oscillations above the edge which comprise
XAFS calculations are visible. Figure reproduced from Rehr et. al [77]. . . 32
3.4 The experimental apparatus at APS. The three ionization chambers are la-
beled as I
0
;I
1
;I
2
. The
ouresence detector is labeled as I
f
. The sample is
placed between I
0
and I
1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.5 Atomic conguration of SnO
2
after RMC analysis of Rietveld renement of
XRDdatafortheperiodicslabmodel(a),theclustermodel(b)andLAMMPS
simulation (c). Unfortunately, after renement, all three models produced
different results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.6 EXAFSanalysisbyRMCusingperiodicboundaryconditionmodelink-space
(a,c, e) and real space (b,d,f). The EXAFS data (red) is compared with the
EXAFS analysis from the RMC (blue). Their difference (green) is shown
shifted by -10 (a,c,e) and -5 (b,d,f) for clarity. Each row corresponds to a dif-
ferent model: the slab model (a,b), the cluster model (c,d) and the LAMMPS
generated particle (e,f). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.7 Partialpairdistributionfunctions(pPDFs)resultingfromRMCanalysisusing
periodic boundary model for Sn-O (a) and Sn-Sn (b). Included pPDFs from
initial coordinates (dotted line). . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.8 EXAFSanalysisbyRMCusingperiodicboundaryconditionmodelink-space
(a) and real space (b). The EXAFS data (red) is compared with the EXAFS
analysis from the RMC (blue). Their difference (green) is shown shifted by
-10 (a) and -5 (b) for clarity. . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.9 Partialpairdistributionfunctions(pPDFs)resultingfromRMCanalysisusing
periodic boundary model for Sn-O (a) and Sn-Sn (b). Included pPDFs from
initial coordinates (dotted line). . . . . . . . . . . . . . . . . . . . . . . . . . 43
ix
3.10 EXAFSdatainfourierspace(a)andrealspace(b)forthecolloidandpowder.
XANES spectra for the colloid and powder (c). . . . . . . . . . . . . . . . . 44
3.11 Comparison of partial distribution functions from cluster models for dry pow-
der and aqueous colloid. The systematically larger amplitude in the case of
the colloid is an artifact of the applied normalization and cancels out in the
RMC analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.12 Comparison of the RMC ts for the different models used for the pellet (a)
and the colloid (b). The shift in the MD model is due to the difference in
lattice constant. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.13 Comparison of Sn-O pPDFs from RMC ts for the different models for the
dry powder (a) and the colloid (b). The systematically higher amplitude for
theMDmodelsin(b)isduetotheappliednormalization(density). Itcancels
out in RMC analysis and is a graphical problem. . . . . . . . . . . . . . . . . 46
3.14 Comparison of the Sn-Sn pPDFs from RMC ts for the different models for
pellets (a) and colloid (b). Like the previous gure, the systematically higher
amplitude in the MD model is due to the applied normalization (density). . . 46
4.1 Schematic of the neural network. The system conguration for all spatial
degreesoffreedomiscalculatedbyinsertinganinstantoftimetintotheinput
layer. A single hidden layer (blue) is used and the output layer (red) gives
theatomiccoordinateswhoserstandsecondderivativesprovidevelocityand
acceleration, respectively. The blue and red lines connecting layers represent
the weights and each node also has a bias. . . . . . . . . . . . . . . . . . . . 53
4.2 Four selected NN (solid blue) and MD (solid red) trajectories show that the
NN can reproduce non-trivial paths. These are for a time period of 1 ps . . . 55
4.3 Comparisonbetweentheneuralnetwork(NN)andMDsimulation. (a)Kinetic
(blue), potential (green) and total (red) energy per particle as a function of
time for MD (square dots) and NN (solid lines). Radial distribution function
and velocity autocorrelation function for MD (blue) and NN (red) are shown
in (b) and (c), respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.4 (a) shows how the three components of the loss function decrease with the
number of epochs. The boundary condition term (red) and energy conser-
vation term (blue) fall more rapidly than the OM action term (green). (b)
Showsthemeansquareerror(MSE)fortrajectory(red)andmomentum(blue)
which was calculated by comparing the MD trajectory to the neural network
predicted trajectory. (c) shows differences between MD and NN trajectories
(red) and velocities (blue) as a function of time. Here = 2ps. . . . . . . . . 56
x
4.5 Scaling of NN for 108-particle system. (a), (b) and (c) are the run time for
every 1000 training epoch as a function of the number of units in the hidden
layer, the size of the time grid, and the number of atoms in the system,
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.1 (a)Theboundaryconditionterm(orange),totalenergyterm(green),Onsager-
Machlup term (red) and total loss (blue) as a function of training epoch.
Note that these values include the corresponding
i
terms. After the initial
few hundred epochs, the loss is dominated by the OM action term. (b) The
meansquareerrorintheposition(blue)momentum(orange)andacceleration
(green). Theslopestendtowardszeroimplyingthatfurtherepochswouldnot
produce signicant improvements in the trajectory. . . . . . . . . . . . . . . 64
5.2 (a) The kinetic (blue), potential (green) and total energies (red) as a function
of time for conventional molecualar dynamics (dashed lines) and our NNMD
(solid lines) after training. (b) The MSE in acceleration (green), momentum
(orange)andposition(Q)comparingtheNNtotheconventionalMDsimulation. 65
5.3 The radial distribution function (a) and velocity autocorrelation function (b)
for regular MD (red) and NNMD after 500,000 epochs (blue). . . . . . . . . 65
5.4 Particletrajectoriesfor4selectedatomsofbothMD(red-dashed)andNNMD
(blue solid). Both trajectories agree at t=0; however, some particle trajecto-
ries diverge signicantly towards the end of the time span. . . . . . . . . . . 66
A.1 Tin dioxide nanoparticle in wate (a)r and pure water (b) with no applied
electric eld. Transient cavities on the order of 100

A
3
and smaller form. . . 84
A.2 Cavity volume analysis. (a) Raw void volume versus time. (b) Cavity volume
versus time. Calculated using voxel chains of 15 or larger. (c) Spherical and
cylindrical growth rates as a function of applied electric eld. Some of the
simulations runs were not run for sufficient lengths of time to see cylindrical
growth. Comparing, for example, the 0.042 curve from part (b), the red star
corresponds to the rst linear region from 600 ps to 750ps while the green
star corresponds to the slope in the region between 760 ps and 850 ps. . . 86
A.3 (a)

E vs T: Blue markers are the data used to calculate the exponential
relaxation coefficient . The temporal data spans an order of magnitude
(60 to 635 ps). The green (blue) line is plotted stretched exponential with
= 0:324 ( = 0:424) and = 90:16 ps ( = 134:6 ps). The blue line is
plotted stretched exponenital (b) log(log(

E) vs logT
onset
. (c)log(log(

E)
vs log(T
onset
min(T
onset
)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
xi
B.1 NN and MD trajectories for four randomly selected atoms in a system of
256 LJ atoms. Atomic trajectories computed by the NN (blue) are nearly
coincident with MD trajectories (red). . . . . . . . . . . . . . . . . . . . . . 91
B.2 a) NN and MD results for potential, kinetic and total energies per atom as a
function of time for a system of 256 atoms. (b) Radial distribution function
and (c) velocity auto-correlation function show the near-perfect agreement
between the NN and MD for a 256-atom system. . . . . . . . . . . . . . . . . 91
B.3 NN and MD trajectories for four randomly selected atoms in a system of
108 LJ atoms. Atomic trajectories computed by the NN (blue) are nearly
coincident with MD trajectories (red). . . . . . . . . . . . . . . . . . . . . . 92
B.4 (a) Shows that the NN results for potential, kinetic and total energy per
particle as a function of time are nearly coincident with the ground truth MD
results. Figures S2 (b) and (c) compare NN and MD results for the radial
distribution function and velocity auto-correlation function for the 108-atom
system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
B.5 NN results for a 256-atom LJ system. (a) loss-function components vs the
number of epochs. (b) Mean-squared error between MD and NN for trajec-
tories (red) and velocities (blue) as a function of epochs. (c) Mean-squared
error between MD and NN trajectories (red) and momenta (blue) vs time. . 93
B.6 Results for the 108-atom system. (a) loss function versus number of training
epochs. (b) mean-squared error between MD and NN results for positions
(red) and velocities (blue). (c) The red curve is the mean-squared error be-
tween MD and NN trajectories (Q) as a function of time. The blue curve is
the error between MD and NN velocities (P). . . . . . . . . . . . . . . . . . 93
xii
Abstract
Cavitation phenomenon in dielectric
uids has been a recent topic of interest in theory and
experiment. Westudyadielectric
uid-nanoparticlesystemsubjectedtoanexternalelectric
eld using molecular-dynamics simulations. Electric elds ranging from 0.042 to 0.25 V/

A
are applied to a water and tin dioxide system. Cavitation is observed in simulations with
both SPC/E water and the hydrogen bonding polarizable model. The cavitation onset time
displays a stretched exponential relaxation response with respect to the applied electric eld
withanexponent = 0:4230:08. Thisisinaccordancewiththeexacttheoreticalvaluefor
systemswithlong-rangedforces. Cavitygrowthratesaredividedintotwophases,aspherical
growth phase and cylindrical. Both are reported on as a function of the applied electric
eld. The structure of the electric eld is analyzed both spatially and temporally. Due
to their relatively high surface area to volume ratio, nanoparticles are potential candidates
to investigate the solid-liquid interface. X-ray Absorption Spectroscopy experiments were
performedonpowderedandcolloidalSnO
2
nanoparticlesusingtheB-12beamlineatArgonne
National Laboratory. The EXAFS data was analyzed using reverse Monte Carlo (RMC)
simulations with starting congurations generated from Rietveld renement results of X-ray
diffraction data.
xiii
The principle of least action is the cornerstone of classical mechanics, special and general
theory of relativity, quantum mechanics, and thermodynamics. A description of the training
of neural networks to learn the principle of least action through minimizing the Onsager-
Machlup action is provided. The network, with slight alterations to the loss function, is
used to solve both initial value problems (IVP) and boundary value problems (BVP). The
resultant phase-space trajectories for the BVP network are in excellent agreement with the
trajectoriesobtainedbyconventionalmoleculardynamics. Thepotentialandkineticenergies
as a function of time as well as radial distribution and velocity autocorrelation functions
computedbyboth neural networksare inexcellentagreementwiththecorrespondingresults
from the ground-truth MD simulation. From the rst epoch, the neural network approach
produces as full trajectory which increases in accuracy as training progresses. The time step
is a factor of twenty larger than the time step in MD simulation for the BVP and 5 times
largerfortheIVP.Theneural-networkapproachcanhelpscientistsdiscovernewphenomena
hidden in data, impute missing data in experiments and simulations and, most importantly,
allow atomistic simulations to reach long timescales.
xiv
Chapter 1
Introduction
The underlying subject of this thesis is the computational study of the movement of
atoms. Countless physical phenomenon can be understood through the investigation of
atomic dynamics and statistics. A perfectly physically faithful simulation is not possible.
Simply using computers implicitily requires limitations on the precision of any calculation
due to the bit-wise representation of numbers. Furthermore, non-classical simulations are
computationally expensive. For example, Ab initio quantum mechanical methods, such
as density functional theory (DFT) calculations are limited to hundreds or thousands of
atoms [1]. Fortuitously, many phenomenon have been successfully explored using classical
dynamics simulations.
We will primarily concern ourselves with classical simulations of N particles in d dimen-
sional boxes with atomic systems differentiated by the rules governing their interactions and
theconditionsplacedontheboundaryoftheirboxes. Theinteractionstaketheformofforce
elds which, when used in conjunction with Newton's second law, allows the calculation of
particle dynamics. The atomic box problem will be solved as an initial value problem (IVP),
boundary value problem (BVP) and with statistical methods (Monte Carlo Methods). The
rst two methods result in atomic trajectories for the the system. In an IVP the initial posi-
tions and velocities are specied and integrated repeatedly to calculate trajectories whereas
a BVP is dened by two congurations of the system separated in time which is solved
1
when the connecting trajectory is identied. The Monte Carlo method does not produce
time sequential states; however, it does produce statistically related states and is useful for
statistical calculations.
Chapters 2 and 3 involve the study of tin dioxide nanoparticles and water. When a polar
liquid, such as water, is subject to a strong non-homogeneous electric eld, cavities form.
These cavities, unlike gas lled bubbles, can be devoid of any atoms. They are sustained by
the negative pressure in the liquid induced by the electric eld. The theoretical background
describing this phenomon was developed by Schneider et al [84,85]. In the work described
in Chapter 2, molecular dynamics (MD) simulations were performed on SnO
2
-H
2
O systems.
The formation of cavity voids were observed when the system was subject to a homogeneous
external electric eld. The presence of the nanoparticle enabled the formation of the voids.
The time between turning on the electric eld and void formation was found to follow a
stretched exponential response (SER).
Thereisconsiderableresearchinterestinsolid-liquidinterfaces. Theformationofelectric
doublelayersisofparticularconcerntoelectrochemists. Consequently,X-rayabsorptionne
structure (XAFS) spectroscopy measurements were performed at Argonne National Labora-
tory (ANL) on SnO
2
pellets and tin dioxide-water colloidal systems. The data was analyzed
using reverse Monte Carlo (RMC) simulations starting from congurations generated from
Rietveld renement results of X-ray diffraction data. Furthermore, MD simulations of a
tin dioxide nanoparticle in both water and vacuum were used to generate starting congu-
rations for Rietveld renement. The resulting RMC analysis was then compared with the
experimental data and analysis.
2
Articial intelligence and machine learning techniques have found applications through-
out the physics and material science. Machine learning is enabling the prediction of material
properties using the data from rst principle calculations [92]. In MD, Neural networks have
beentrainedtoactasforceelds[57]. Recurrentneuralnetworkshaveevenbeenusedtoan-
alyze the trajectories for ligand-protein systems [11]. We have trained a feed-forward neural
network to learn the trajectories for 3-D Lennard Jones atomic systems. Our method relies
on minimizing the Onsager-Machlup (OM) action. The neural network can successfully cal-
culate the trajectories between two points in conguration space, solving the BVP (Chapter
4). In addition, the same network was trained to solve the IVP, akin to conventional MD,
with a altered loss function (Chapter 5).
To aid the reader, a review of pertinent background information follows. A brief discus-
sion of the calculus of variations will be provided. The Euler-Lagrange (EL) equations will
be used to derive Newton's equations of motion from both classical and Onsager-Machlup
action. From there the basics of MD will be discussed.
Variation Principles
Action principles are ubiquitous throughout physics. Richard Feynman formulated non-
relativistic quantum mechanics through action principles as his thesis [23]. Here we focus
on the formulation of Newton's equations of motion from the principle of stationary action.
Hamilton's principle of stationary action states:
3
The motion of the system from time t
i
to t
f
is such that the line integral (also
known as the action integral):
I =
∫
t
f
t
i
Ldt =
∫
t
f
t
i
TVdt (1.1)
has a stationary value for the actual path of the motion [28].
Before we derive the equations of motion, we will rst dene what it means for a functional
to be stationary. Assuming standard dot notation, _ y :=
dy
dt
, let there be some function of the
trajectory with at least two existing time derivatives, f = f(y; _ y; y;t) with y = y(t). We can
then dene the functional:
J =
∫
t
f
t
0
f(y; _ y; y;t)dt
We want the value of J to remain constant with respect to variations of the trajectory, y,
over the integrated time interval. We shall augment y(t) by adding an arbitrary function
(t) with the condition that both (t) and _ (t) vanish at the boundary points t
0
;t
f
:
y(t;) = y(t)+(t) (1.2)
whereisanarbitraryparameter. WethendeneJ asstationarywhenthederivativeunder
the integral vanishes:
dJ
d
=
∫
t
f
t
0
@f
@y
@y
@
+
@f
@_ y
@_ y
@
+
@f
@ y
@ y
@
+
@f
@x
@t
@
dt = 0
4
Of course this simplies, noting that the derivative of t is zero. We can integrate the
second term third term by parts to yield:
∫
@f
@_ y
@_ y
@
dt =
∫
@f
@_ y
@
2
y
@@t
dt =
@f
@_ y
@y
@

t
f
t
0

∫
@y
@
d
dt
@f
@_ y
dt (1.3)
∫
@f
@ y
@ y
@
dt =
@f
@ y
@
2
y
@@t

t
f
t
0

@f
@ y
@y
@

t
f
t
0
+
∫
@y
@
d
2
dt
2
@f
@ y
dt (1.4)
Recall that
@y
@
= and ;
′
vanish at the boundaries by denition. Now we can set
dJ=d

=0
:
dJ
d

0
=
∫ (
@f
@y

d
dt
@f
@_ y
+
d
2
dt
2
@f
@ y
)
(t)dt = 0 (1.5)
@f
@y
=
d
dt
@f
@_ y

d
2
dt
2
@f
@ y
(1.6)
Equation 1.5 only holds for arbitrary if the term inside the parenthesis is zero for all t. It
is straightforward to extend the above calculations to the case where y is a d dimensional
quantity; however, to avoid spoiling the reader, it shall be left as an exercise. Equation
1.6 is known as the Euler-Lagrange (EL) equation. If we let f equal the Lagrangian, L =
∑
i
m
2
_ q
2
V(q),weimmediatelyrecoverNewton'sequationsofmotionpluggingintoEquation
1.6:
m q =
@V
@Q
=F (1.7)
The necessary and sufficient condition for the E.L equations to be satsied is that the action
is stationary. Classical action cannot be a maximum for a physical trajectory; however,
5
trajectories can be local minima or even saddle points [30,63]. There exists an action which,
as we will show, produces Newton's equations of motion only at a global minimum. The
Onsager-Machlup (OM) action is given by:
S
OM
=
∫
t
f
t
0
(
M
ij
d
2
q
j
+@
q
i
V
)
2
dt =
∫
t
f
t
0
Q
k
Q
k
dt (1.8)
Q
k
= M
jk
d
2
q
j
+@
q
k
V (1.9)
where M
ij
= m
i
;(0) for i = j;(i̸= j) , m
i
is the mass of the particle corresponding to the
ith degree of freedom, q
j
is the coordinate of the jth degree of freedom, d refers to a time
derivative, @
X
=
@
@X
and summation is implied by repeated indices [63]. From inspection
of Equations and it is clear that S
OM
has a global minimum of zero which occurs when Q
k
vanishes for all t. Incidentally, when Q
k
(t) = 0 Equation reduces to Equation 1.7. In the
case of solving BVPs through optimization of trajectories (which will be addressed in detail
inChapter4)theglobalminimumofOMactionallowsitsuseasametricforsolutionquality
as opposed to classical action where the numerical value is not, in general, useful.
Molecular Dynamics
Data: q
i
(t = 0);v
i
(t = 0);i2 [0;N1]
Result: q
i
(t = t
f
);v
i
(t = t
f
)
initialization: assign initial value8q
i
; _ q
i
;
while t < t
f
do
a
i
Calc Acceleration(q
i
) ;
q
i
;v
i
Time Integration(q
i
;v
i
;a
i
);
q
i
Apply Boundary Condtions(q
i
);
end
Algorithm 1: Basic Molecular Dynamics Algorithm
6
Figure 1.1: Example of periodic boundary conditions. The simulation cell's (green) particle
(red and blue) moves (top box) from one boundary to the opposite (bottom).
The most powerful tool available for studying the atomic IVP is molecular dynamics.
MD has provided advancements in material science [41], chemistry [12], biology [42] and
physics[55]. ThebasicMDalgorithm(Algorithm1)requiresaninitiallistofatomicpositions
and velocities, an integrator to advance to the next time point, a force eld describing
the atomic interactions and a simulation box with specied boundary conditions. As a
computational approach, it requires the discritization of time to represent the trajectories
onacomputer. Aconstanttimestepof∆tmustbeselectedtorunanMDsimulation. While
the time step depends on the algorithm and the system under study, usually a time step on
the order of 1 fs is used as this satises the requirement of being smaller than the shortest
vibrational period in the system.
Atomicpotentialsorforceeldscanbedevelopedinvariousways. Theparametersfound
in models such as Lennard Jones or Buckingham potentials can be selected as a "rst-guess"
7
potential for a system. A cycle of simulating, comparing the results to experimental data
andtuningisonewaytogenerateaforceeldforaparticularmaterial. Someforceeldsare
developed using the results of ab initio quantum mechanics calculations instead. A forceeld
for tin dioxide and water was created using a combination of density funtional theory and
LJ tuning [7]. Recently, neural networks been trained as force elds for metals [57] and
water [99] in lieu of conventional models.
In a MD algorithm, the force calculation step is the most computationally expensive
portionoftheprogram. Thisisduetothecalcluationsoftwooreventhreebodyinteractions
which, in general, scales as N
2
or N
3
respectively. Consequently, multiple methods have
been developed to allow for the simulations of larger systems. The rst, is the application
of periodic boundary conditions (PBC) (See Figure 1.1). The simulation cell is replicated to
innity. Whenaparticlemovestotheedgeofthesimulationbox,itreappearsontheopposite
side. This allows even modest simulation of only a few thousand atoms to approximate bulk
systems. Several potential issues with force elds immediately arise with PBC. There are
an innite number of pair-wise interactions which now need to be calculated. The usual
solution is the implementation of minimum image convention (MIC) in which each particle
i only interacts with the particle corresponding to min(r
im
ij
;r
real
ij
). Short-range force elds,
such as the LJ potential, can be assigned a cutoff radius, r
c
, after which the force is assumed
to be zero. In simulations with both PBC and MIC, the cutoff radius should be less than
thehalf the shortest boxdimension. In order toavoiddiscontinuitiesinthe potentialenergy,
potentials can be shifted to have their zero values occur exactly as r
c
or have r
c
values
large enough that the value is sufficiently close to zero. There exist long range forces, such
8
as those arising from the Coulombic potential which cannot accurately be assigned a nite
cutoff radius.
The Coulombic potential can be estimated using the Ewald sum techinque. The inter-
action energy between ions in a simulation cell with periodic boundary conditions is given
by:
V
ion
=
1
2
∑
⃗ n
(
∑
i
∑
j
z
i
z
j
j⃗ r
ij
+⃗ nj
1
)
+
1
2
(
∑
i
∑
j≠=i
z
i
z
j
j⃗ r
ij
j
1
)
(1.10)
where z
i
is the net charge, ⃗ n is the lattice vector corresponding to different image cells and
⃗ r
ij
is the separation between the two ions in question. This sum is conditionally convergent
which means that the answer depends on the order in which the terms are added. We shall
add the terms radially outwards from the simulation cell i.ej⃗ nj = L,j⃗ nj =
p
2L, etc. In
the limit of large ⃗ n we start to approximate a sphere surrounded by a medium. The general
approach for the method is as follows. Center a Gaussian charge distribution around every
point particle with opposite total charge. Calculate the interaction of this new short-range
potential. Add in a second, canceling charge distribution with the same total charge as the
initial ion but the shape of the Gaussian charge density. Take the fourier transform of the
cancel charge distribution, sum and transform back to real space. Finally add the term to
the total effect of the ionic potential. The form of the two added potentials is arbitrary;
however, Gaussian's are notoriously pleasant to work with.
Modern MD software most commonly uses the velocity verlet algorithm to perform the
timeintegrationandupdatepositionsandvelocities[97]. VelocityVerletiscarriedoutinthe
followingfashion. Firstcalculatetheacceleration

⃗ q
t
fortimetusingthecurrentconguration
9
⃗ q
t
. Nextcalculatethehalf-stepofthevelocities
_
⃗ q
t+1=2
=
_
⃗ q
t
+
∆t
2

⃗ q
t
. Third,updatethepositions
to the next timestep. Recalculate the accelerations and nally update the velocities with a
second half step
_
⃗ q
t+1
=
_
⃗ q
t
+
∆t
2
+(

⃗ q
t
+

⃗ q
t+1
). This can be simplied into two lines producing:
⃗ q
t+1
=⃗ q
t
+
_
⃗ q
t
∆t+
1
2

⃗ q
t
∆t
2
(1.11)
_
⃗ q
t+1
=
_
⃗ q
t
+
1
2
(

⃗ q
t
+

⃗ q
t+1
)∆t (1.12)
This implicitly assumes that the accelerations are only a function of the positions at time
t+∆t which turns out to be an acceptable approximation for many systems. The Velocity
Verlet algorithm has local errorO(∆t
4
);O(∆t
2
) for the position and velocity respectively.
In simulating solids, the ideal crystal lattice is an obvious starting point for MD simu-
lation. For liquids, there are multiple options. Depending on the density desired, randomly
lling the simulation cell with molecules can lead to an explosive start if two atoms or
molecules are placed too close together. It is more simple to place the atoms or molecules in
a lattice with the correct density and allow the system to evolve from there. Initial velocities
are assigned from a Boltzman distribution corresponding to the ideal temperature. Enforc-
ing zero net linear momentum can then be done by sampling N1 velocities and using the
degrees of freedom of the nal particle to ensure total momentum vanishes.
Wediscussedperiodicboundaryconditionsabovebutthereexistothertypesofboundary
conditions. When a particle reaches the edge of the box, there are a 4 possibilities: the
particle can bounce off the edge, the particle can leave the simulation box and be removed
from the system (xed boundary conditions), the simulation box can expand to include
10
the new space the particle occupies (shrink wrapped boundary conditions) or nally the
aforementioned PBC.
11
Chapter 2
Electrostrictive Cavitation in Water Induced by a
SnO
2
Nanoparticle
Introduction
The
1
formation and use of cavities within
uids is a rich area of research and practical appli-
cations. Nanobubbles have been used for water treatment, aiding fermentation and targeted
delivery of pharmaceuticals [2,94]. The water jets formed from collapsing nanobubbles can
cause structural damage and erosion to important components and, more productively, aid
in manufacturing materials [62]. Bubble collapse induced sonoluminescene has been ob-
served [102]. Typically cavities arise in the form of gaseous bubbles; however, in this paper
the formation of cavities due to electrostrictive forces will be discussed.
Cavity formation due to negative pressure in water has been the subject of experimental
and theoretical work [5,85]. According to classic nucleation theory, the energy required to
generate a spherical void of radius R in a
uid with surface tension and pressure P, is
equal to [7]:
(R) =
4
3
R
3
P +4R
2
: (2.1)
1
Chapter originally published in ACS Omega [39]
12
The term for the cavity's internal vapor pressure has been excluded as the cavities discussed
in this paper have a negligible vapor pressure compared to the internal
uid pressure. Equa-
tion 2.1 implies, for negative pressures, the existence of a critical radius above which it will
be energetically favorable for the cavity to continue to grow. Negative pressure can be in-
duced in water (or any dielectric
uid) through the application of a spatially inhomogeneous
external applied electric eld. This induced negative pressure in a a dielectric
uid is given
by:
P =
1
2
ϵϵ
0
E
2
(2.2)
where ϵ
0
is the vacuum permittivity, ϵ is the dielectric constant of the liquid media, E is the
applied electric eld and is a constant which ranges between 1.3 and 1.5 for most polar
liquids [84]. Negative pressure in water has been measured through the Berthelot method
withpressuresrangingfrom-3.4to-12MPaacrosstemperaturesfrom275to300K;however,
theaccuracyoftheseresultshasbeenquestioned[16]. Negativepressuresexceeding-100MPa
wereachievedusingwaterinclusionsinquartz[108]; however,inthestudiesdonebyGreenet
al.[31],eachinclusionwasonlystudiedafewtimesorevenonlyonce. Abetterunderstanding
of a cavity formation statistics was provided by Azouzi et al. [5] in which a single inclusion
was studied repeatedly. Direct observation of cavities through Rayleigh scattering has been
proposed but to date no experiments have been carried out [86]. The Shneider, Pekker
and Fridman theory of void formation has been experimentally veried [86]. Water subject
to pulses with rise times of 600 ps to 3 ns were found to generate nanocavities through
electrostrictive forces [80]. Cavitation did not occur with longer rise times as the liquid had
13
sufficient time to respond to the forces without rupturing. The theoretical explanation of
cavity formation in the presence of a strong, inhomogeneous electric eld is found in the
work of Shneider et al [85].
The previous theoretical work cited above has shown the genesis of cavities when a
spatially inhomogeneous external eld is applied; however, through molecular dynamic sim-
ulations, we have witnessed the formation of cavities in water subject to a spatially constant
applied electric eld. Inhomogeneities in the local electric eld occur due to the presence
of the tin dioxide nanoparticle. This phenomenon in turn leads to void cavitation forma-
tion. Simulations of a 1.5 nm SnO
2
nanoparticle immersed in water were performed using
LAMMPS [71,88,93]. Electric elds ranging from 0.042 to 0.25 V/

A were applied along the
x axis. While the simulated electric eld strengths are quite high, water has been exper-
imentally subjected to electric elds between 0.1 and 0.5 V/

A using electrode tips on the
order of hundreds of nanometers [79].
Results and Discussion
The local electric eld was calculated at the location of each atom by dividing the Coulomb
force felt by the atom by the atomic charge. After the initial application of the electric eld
(occurs over 1 fs), the average x component of the effective electric eld rapidly increased
until reaching a saturation value (Figure 2.1 (a)). We have dened the saturation time, t
s
to
be the time that required for the effective electric eld to reach 95% of its maximum value.
The average electric eld remained stable with oscillations of at most 5% for t
p
< t < t
onset
.
The changes in the electric eld near the onset of cavitation will be discussed below. The
14
0 10 20 30
0
0.5
1
1.5
2
0.250 V/A
0.200 V/A
0.160 V/A
0.125 V/A
0.109 V/A
0.096 V/A
0.073 V/A
0.063 V/A
0.052 V/A
(a)
0 0.1 0.2 0.3
1
1.5
2
(b)
0 0.1 0.2 0.3
0
5
10
15
20
(c)
50 100 150
0
0.5
1
1.5
0 ps
1 ps
3 ps
5 ps
10 ps
(d)
Figure 2.1: Initial magnitude and structure of electric eld. (a) The average x component
of the electric eld versus time for different applied electric eld strengths. (b) The time
it takes the electric eld to reach 95% of the saturation value. (c) The magnitude of the
saturationvaluefortheelectriceldversusappliedelectriceld. (d)Thexdependenceofthe
x component of the electric eld at different time steps correspond to the initial application
of the eld, the rise and the saturation for an applied eld of 0.125 V/

A.
saturationvaluedecreased(Figure2.1(b))andthepulsetimeincreasedexponentially(Figure
2.1 (c)) as the applied electric eld decreased. The time resolution in these calculations is
one ps. For E
app
> 0:125 V/

A, a smaller time scale is needed for accurate measurements.
The saturation time increases from between 3 and 4 ps for a 0.25 V/

A applied eld to 20 ps
for the 0.0524 V/

A eld. Analysis of the logE vs logt
p
of the data in Figure 2.1 yielded a
15
(a)
0 2 4 6
-2
-1
0
1
= 0.424 0.09
(b)
0 500 1000
0
0.05
0.1
0.15
0.2
= 0.424 0.09
(c)
Figure 2.2: (a) Snapshot of the system just after cavitation. (b)

E vs t: Blue markers are
the data used to calculate the exponential relaxation coefficient . The temporal data spans
an order of magnitude (60 to 635 ps). The red line is the tted line of best t. (c) The blue
line is the plotted stretched exponential function with = 0:424 and = 134:6 ps.
16
slope of -0.97 which indicates an exponential dependence. The total electric eld increased
linearlly with the applied electric eld; however, there was a scaling factor of approximately
3.6! This increase in total electric eld strength is attributed to the polarization of the water
molecules and resulting dipole moment of the periodic simulation boxes [15]. The effective
electric eld for a constant applied electric eld is given by:
E
eff
=E
app
+E
dp
(2.3)
where E
dp
is the electric eld due to the dipole moment of the box:
E
dp
=
p
4ϵ
0
∑
;r
2cos
2
sin
2

r
3
(2.4)
where p is the total dipole of the box, the sum is over all dipoles within the box, r is the
distance from an individual dipole to the position where E is being measured, and is
the angle between the dipole and the distance vector between the dipole and measurement
point [22]. If we letp = E
app
where is the effective polarizability of the box, the effective
electric eld becomes:
E
eff
=E
app
(
1+

4ϵ
0
∑
;r
2cos
2
sin
2

r
3
)
: (2.5)
Noting that our electric eld scales as 3.6 times the applied electric eld, the second term in
Eq 2.5 must evaluate to 2.6. For steady state electric elds, the dipole component has been
found to be larger than the applied electric eld itself which is in agreementwith our results.
This would suggest that the dipole relaxation of the water molecules is responsible for the
17
electriceldpulsesseeninFigure2.1(a). Figure2.1(d)showsthexcomponentoftheelectric
eld at t = 0, 1, 5 and 10 ps for E
app
= 0:125 V/

A. From these results it can be seen that,
during this initial pulse, only the magnitude of the electric eld changes while the structure
of the electric eld does not change. In the absence of a nanoparticle, the pulse occurs on
similar timescales (between 3 and 4 ps for 0.25 V/

A); however, cavitation is not observed.
This reaffirms that the nanoparticle is required for cavitation in this system. Furthermore,
this suggests that alteration of the local electric eld is what induces the cavitation.
Figure 2.2 (a) shows a snapshot of a simulation shortly after cavity genesis. The time of
cavitation onset (T
onset
) was found to relate to the applied electric eld through a stretched
exponential relaxation response:

E = exp
(
(t
onset
=)

)
(2.6)
(2.7)
where

E = E=E
max
is the electric eld divided by the maximum applied electric eld
(E
max
=0.25 V/

A). The value for was 0:424 0:09 (Figure 2.2). The value for were
found by plotting log(log(

E)) vs logT
onset
and nding the weighted line of best t (Figure
2.2b,c). ThettingweightswerebasedontheerrorassociatedwitheachT
onset
measurement
(1ps). These values are very close to the "magic" values found in stretched exponential
response theory.
Stretched exponential relaxation has been used to explain relaxation in glasses and di-
electrics[69,70,72]. Acrossdecadesofexperimentsinawiderangeofsystems, afew"magic"
values for have arisen. One of the early explanations was the trapping model, in which
18
(a) (b)
(c) (d)
(e) (f)
400 600 800 1000
0
0.5
1
1.5
2
10
4
(g)
Figure 2.3: Evolution of cavity with an applied 0.42 V/A electric eld. Snapshots at T =
700 (a), 800 (b), 900 (c), 1000 (d), 1050 (e), 1100 (f) ps. Part (g) shows the cavity volume
as a function of time with red Xs over the data points corresponding to parts a-f.
19
entropic excitations diffuse freely through phase space until they reach sinks in which they
becomestuck[70]. Thestretchedexponentialnaturethenarisesasaresultoftheexcitations
near the sink being diminished quickly while with farther away excitations taking increasing
lengths of time due to their distance from the sink. This model predicts a value for of:
=
d
d+2
(2.8)
where d is the dimensionality of the system. Phillips, using the wealth of experimental data,
proposed an axiomatic model with the axiom:
=
d

d

+2
=
d
d+2=f
(2.9)
where d

= fd and f is the number of short range forces available for relaxation divided
by the number of long ranged forces. A topological justication of this model can be found
in [56]. For a system like the one studied in this paper with one short ranged interaction per
atom type (Lennard Jones/Buckingham) and one long ranged force (the coulomb force) one
would expect f = 1=2 yielding = 3=7. Further discussion on the interpretation of can
be found in the Appendix A.
Figure 2.3 shows the evolution of the cavity with respect to time for an applied eld of
0.042 V/

A. Snapshots shown are from 500 (a), 650 (b), 750 (c), 850 (d), 1000 (e) and 1100
(f) ps. The cavityvolumeas a function of time is shownfor the duration of the simulation in
Figure 2.3 g. Between 500 and 650 ps the cavity begins to form in front of the nanoparticle.
Between approximately 650 and 730 ps the cavity forms and expands roughly spherically at
20
a rate of approximately 24.7

A
3
/ps. The cavity then expands cylindrically at a rate of 92.78

A
3
/ps for 150 ps. After the initial sphere phase, the cavity grows in a cigar like fashion:
First a small tendril of the cavity extends in the direction of the electric eld, away from
the midpoint of the nanoparticle and then it expands outwards until the segment reaches
the same cyclindrical radius as the rest of the cavity. The cavity stabilizes for a period of 90
ps until continuing to grow at a reduced rate of 29

A
3
/ps. Between 1050 and 1100 ps, the
outermost (from the nanoparticle) segment broke off and collapsed. Yet, in that time span
the overall cavity volume continued to expand. During this new regime, the volume increase
comes from an expansion in the radial direction.
Anexplanationforthecylindricalgrowthofthevoidcavityisfoundintheliterature[84].
The equation for square magnitude of the electric eld near a spherical void is given by:
E
2
(r;) = E
2
app
(
1+(3cos
2
+1)
(
ϵ1
2ϵ+1
)
2
R
6
r
6
(5cos
2
1)
R
3
r
3
ϵ1
2ϵ+1
)
(2.10)
with donating the angle between r and the external applied electric eld E
app
. The volu-
metric force near the pore ( R < r < 2.5 R) then becomes:
F
3
2
ϵ
0
ϵE
2
0
(ϵ1)R
3
(2ϵ+1)r
4
(
5cos
2
1)
2(3cos
2
+1)
R
3
(ϵ1)
r
3
(2ϵ+1)
)
:
Notice that when = 0, which is inline with the external applied electric eld the electric
eld is minimized and the force is positive (pointing radially outwards). As increases,
the electric eld in this region increases and the force becomes less positive until becoming
21
negative between =8 and =4. This would extend explain the cigar shape found in the
simulations.
0 500 1000
0
0.5
1
1.5
2
10
4
0.250 V/A
0.200 V/A
0.160 V/A
0.125 V/A
0.109 V/A
0.096 V/A
0.073 V/A
0.063 V/A
0.052 V/A
0.042 V/A
Figure 2.4: Cavity volume analysis. The cavity volume as a function of time for different
applied electric elds. In many of the curves, two distinct regions of near linear growth
appear. The rst corresponds to the spherical growth phase while the second corresponds
to the cylindrical growth phase.
There were two apparent growth phases in cavitation genesis: an initial spherical phase
and a longer lasting cylindrical growth phase. From inspection of Figure 2.4, two distinct
growth slopes emerge. Generally, both phases grew at increased rates as the electric eld
strengthwasincreased. Thevolumegrowthappearstogrowlinearly;however,thedataistoo
noisy to determine this denitively. A description of how the void volumes were calculated
can be found in Appendix A.
The x component of the electric eld was calculated by dividing the coulomb force by
atomic charge for each hydrogen and TIP4P dummy atom of water and each tin and oxygen
atom corresponding to tin dioxide (Figure 2.5). The x component of the electric eld was
averaged and compressed into the y-z plane (Figure 2.5 a,d,g) and along the x axis (Figure
2.5 b,e,f). The snapshots were taken from the E
app
= 0.125 V/

A run. The nanoparticle
22
(a)
50 100 150
1.45
1.5
1.55
1.6
(b) (c)
0

93
y [A]
93

0
z [A]
0.8
1
1.2
1.4
1.6
1.8
(d)
50 100 150
1.45
1.5
1.55
1.6
(e) (f)
0

93
y [A]
93

0
z [A]
0.8
1
1.2
1.4
1.6
1.8
(g)
50 100 150
1.45
1.5
1.55
1.6
(h) (i)
Figure 2.5: The x component of the electric eld at 10 ps (a, b, c) 40 ps(d, e, f) and 100
ps(g, h, i) for an applied eld of 0.125 V/

A. (a,d,g) Heat map of E
x
the averaged over the
x axis. (b, e, h) Plot of E
x
vs x position averaged over y and z axis. The red lines show the
locations of the edge of the nanoparticle. (c, f, i) Snapshot of the cavity and nanoparticle at
timestep under consideration. The electric eld spikes in the region around the cavity while
decaying as distance is increased from the cavity.
disrupts the otherwise constant electric eld which produces the inhomogonuity required for
cavity formation. Once the bubble has formed (t=100ps), the average electric eld in the
water spikes heavily in the plane around the bubble in agreement with Eq 2.10 while seeing
a slighly reduced electric eld near the leftmost edge. The internal electric eld strength in
the absence of an externally applied electric eld has been found to be between 1.5 and 2.5
23
V/

A [21] with only modest alterations due to a comparatively weaker E
APP
. This is in line
with the electric eld results in Figure 2.5.
Summary
We have observed cavitation in MD simulations of a SnO
2
nanoparticle immersed in water
modeled by SPC/E and the hydrogen bonding polarizable force elds. The results from
these simulations describe a general phenomenon. Cavity growth rates indicate spherical
and cylindrical growth phases. These growth phases are analyzed as a function of the ap-
plied electric eld and the simulation data are consistent with the theoretical calculation
for a cigar-shaped cavity near the nanoparticle. A more precise anaylsis of the cavity shape
and evolution will be the subject of future work. The cavitation onset time as a function
of the applied electric eld displays a stretched exponential relaxation response with a uni-
versal value for the exponent, = d

=(d

+2) = 3=7, for systems with long-ranged forces.
Identifying the long timescale value of provides evidence for the most likely mechanism
behind cavitation, i. e. the thermally activated cavities diffuse into the main cavity by the
nanoparticle. These results can be veried by creating a solution of tin dioxide nanoparti-
cles suspended in water and creating a sufficiently strong electric eld through the use of a
fast pump laser and observing cavitation with the free electron laser facility at SLAC Na-
tional Accelerator Laboratory [104]. The Linac Coherent Light Source (LCLS) at SLAC has
been used to successfully study fragile proteins [44] and to investigate femtosecond timescale
phenomenon [87]. The same techniques could be used to probe our system. Experimental
verication of these simulations will provide an interesting direction for further research.
24
Computational Methods
Cavitation was observed with both SPC water and polarizable water force elds; however,
the simulations discussed within this paper were all done with a combination of the SnO2
forceeld developed by Bandura et al [7] with the Hydrogen Bonding Polarizable (HBP)
water force eld [40]. Further details on the force elds can be found in Appendix A.
Bulk crystal SnO
2
was simulated to generate nano nanoparticles. The equilibrated
nanoparticle was then placed into a prepared water box. Finally, the combined system
was used to run the electric eld simulations described below. All simulations were done
using LAMMPS . The default velocity-Verlet integrator was used with a time step of 1 fs.
Periodic boundary conditions were used in every simulation.
The199,970atomcombinedwaterandnanoparticlesystemwasusedatthestartingpoint
for the simulations described in this paper. Simulations ran from 100 ps to 1.2 ns. Uniform
electric elds ranging from 0.042 to 0.25 V/

A were applied along the x direction. Electric
eld strengths above 0.25 V/

A were ignored because force elds which can deal with bond
breaking like ReaxFF are required as water will disassociate at such high eld strengths.
The Ewald summation method was used with tin-foil boundary conditions to calculate the
Coulomb force with a maximum relative error in the forces of 1.010
5
.
Acknowledgments: This work was supported by the U.S. Department of Energy, Of-
ce of Science, Basic Energy Sciences, Materials Science and Engineering Division, Grant
DE-SC0018195. This work was performed, in part, at the Center for Integrated Nanotech-
nologies, an Office of Science User Facility operated for the U.S. Department of Energy
25
(DOE) Office of Science. Sandia National Laboratories is a multi-mission laboratory man-
aged and operated by National Technology and Engineering Solutions of Sandia, LLC., a
wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOEs National Nu-
clear Security Administration under contract DE-NA-0003525. The views expressed in the
article do not necessarily represent the views of the U.S. DOE or the United States Gov-
ernment. We would like to thank M. Winterer for his aid in developing an experiment
proposal.
26
Chapter 3
XAS study of colloidal SnO
2
Introduction
Colloids play an important role in soils or biological cells and are used in many products
of chemical industry such as polymers, paper, varnishes, paints, pastes, for surface nishing,
in ceramic processing, production and processing of nanoparticles, and in food industry [46].
The solid-liquid interface plays a major role in photocatalytic water splitting [6], electro-
catalysis, solid-electrolyte batteries and fuel cells [89]. Furthermore, the structure of the
electrochemical double layer has a profound in
uence on reaction kinetics [59]. Using X-
ray absorption ne structure (XAFS)techniques at Argonne National Laboratory (ANL) we
studied both powdered and colloidal SnO
2
nanoparticles. We recorded nite differences be-
tween the two systems; however, the difference in signal was small to the point of making
analysis of the liquid-solid interface difficult. The experimental data was supplemented with
a LAMMPS simulation of a 6 nm diameter SnO
2
nanoparticle in water and vacuum.
Background
Electrochemical double layer
One of the most important features of the liquid-solid interface is the electrochemical double
layer (EDL). Since its inception, the EDL has been a topic of great interest in the eld of
27
Figure 3.1: Schematic diagram of a Electrode-electrolyte interface according to GCS theory.
The positive charge builds on the surface of the working electrode or, in our case, the
nanoparticle (Green) while the Stern layer is formed by negative ions in the solution (Blue).
Between the IHP and OHP is the solvated atom layer. This is followed by the diffuse atom
layer of the EDL and nally the bulk solution.
electrochemistry. A complete, experimentally veried molecular level understanding of the
EDL would represent a major breakthrough in the elds of batteries, solar panels and fuel
cells [90,98,107]. Despite a great deal of experimental effort, the structure of the potential
within the EDL has remained elusive.
Much of the early investigations into the EDL structure involve electrocapillarity tech-
niques [29]. In recent years a wide range of experimental techniques have been used to
analyze different components of the EDL. The physical structure of the EDL has been di-
rectly investigated through the use of nuclear magnetic resonance (NMR) spectroscopy to
calculate ion numbers in the EDL of super capacitors [32]. Others have used scanning probe
microscopy (SPM) to obtain atomic resolutions of the adsorption of ions into the electrode
interface[26]. Brownetal. developedatechniquetodeterminetheabsolutesurfacepotential
as a function of specic cations, the location of the shear plane and the capacitance of the
Stern layer from XPS data [14]. The direct analysis of the surface potential was a result of
28
XPS on a stream of a nano particles suspended in a liquid microjet. Surface X-ray adsorp-
tion spectroscopy (sXAS) has been used to look at water solvent molecules at the water-gold
electrode surface [96]. The water molecules were found to have an ice like structure near the
surface when a no potential or a positive potential was applied. When a negative potential
was applied, approximately 50% of the surface molecules laid
at on the electrode surface
while the rest of the molecules behaved as bulk water.
ManymodelsoftheEDLhavebeenproposedincludingtheoriginalHelmholtsmodel,the
Gouy-Chapman (GC) model, and theGouy-Chapman-Stern (GCS)model[67]. The original
Helmholtz model describes the EDL like a parallel plate capacitor with charge density:
=
ϵϵ
0
d
V (3.1)
and differential capacitance:
@
@V
=
ϵϵ
0
d
(3.2)
which,whilesimple, almostuniversallyprovidespredictionsincompatablewithexperimental
results [8].
The Gouy-Chapman model describes a diffuse EDL with ionic concentration diminishing
rapidlywithdistancefromtheinterfaceuntilthebulkconcentrationisreached[95]. Consider
an innite electrode extending in thex;y;z directions with a face at the xy plane
surrounded by an innite electrolyte. If the electrolyte is split into layers dz thick, there
exists a layer, dz
0
such that the concentration is arbitrarily close to the bulk concentration.
29
At equilibrium, each layer will be at equilibrium with its neighbor; however, they will have
a different concentration because the electric potential and therefore the chemical potential
change with distance from the electrode. The concentration of ionic species i in the layer
with electric potential ϕ relative to the bulk system is:
n
i
= n
0
i
exp
(
z
i
eϕ
k
b
T
)
(3.3)
where n
0
i
is the bulk concentration, z
i
is the Z number for the ionic charge, e is the charge of
an electron and T is the temperature. The charge density can then be solved for and used
to calculate the electric potential by solving the Poisson-Boltzmann equation:
(
dϕ
dx
)
2
=
2k
b
T
ϵϵ
0
∑
i
n
0
i
(
exp
(
z
i
eϕ
k
b
T
)
1
)
: (3.4)
whichinturncanbefurthersimpliedbyassumingasymmetricelectrolytewithtwospecies:
dϕ
dx
=
(
8k
b
Tn
0
ϵϵ
0
)
1=2
sinh
(
zeϕ
2k
b
T
)
where n
0
is the numerical concentration of both types of ions and z is the magnitude of the
charge. This differential equation can be solved analytically; however, the typical take away
is that when the potential on the electrode is small (ϕ 50=z mV at room temperature ) it
decreases exponentially with distance from the electrode.
The GC model was a signicant improvement upon the original Helmholtz model; how-
ever, it treats ions as point particles which allows for unrealistically high concentrations of
ions near the interface.
30
Figure 3.2: An incoming beam of X-rays passes through a material of thickness d. The
intensity is then measured by the detector.
The GCS model solves the issue of ionic concentration near the electrode by introducing
a layer of electrically bonded ions connected to the interface bound by what is known as
the inner Helmholtz plane (IHP). The IHP is followed by a layer of solvent covered ions
over which the potential drops linearly. This layer is capped with the outer Helmholtz plane
(OHP).PasttheOHPisthefamiliardiffuselayerofGCtheorywithanexponentialdecaying
potential. A diagram of the GCS model is shown in Figure 3.1.
X-ray Absorption Fine Structure (XAFS)
The scattering and absorption of X-rays are used to study atomic structures. The ordered
structure of crystals are determined using X-ray diffraction. X-ray absorption ne structure
(XAFS)spectroscopyusetheoscillationsintheX-rayabsorptioncoefficientjustabovetheX-
rayabsorptionedgetodeterminetheatomicelectronic, structuralandvibrationalproperties
of non-crystalline materials.
In a typical X-ray absorption spectroscopy experiment (Figure 3.2) a sample material of
thickness d is placed in front of an X-ray beam of photon energy E and intensity I
0
. The
transmitted beam signal, I = I
0
e
x
, is then compared with the initial intensity to calculate
31
Figure 3.3: Schematic of X-ray absorption coefficient versus incident photon energy for an
arbitrary material. The trend of decreasing with increasing E is shown. The edges are
labeled (K,L) and the oscillations above the edge which comprise XAFS calculations are
visible. Figure reproduced from Rehr et. al [77].
the so-called X-ray absorption coefficient, (E) =
d
dx
ln(I). The photon energy is then
varied to measure the energy dependence of . A represenation of the typical (E) response
created by Rehr et. al is shown in Figure 3.3 [77]. Across different materials, has three
common characteristics. The rst is a tendency to decrease as energy increases. The other
two are the appearances of sharp increases known as edges which are then followed by a
period of oscillations. Edges correspond to the excitiations of core electrons to unoccupied
states with energies larger than the Fermi energy. The oscillations contain the structural
information about the material. The XAFS spectrum is dened as:
(E) = [(E)
0
(E)]=∆
0
(3.5)
where
0
(E) is the smoothly varying background absorption and ∆
0
is a normalization
factor which is due to the increase in total atomic background absorption from the edge
being investigated [77]. The normalization factor is usually estimated by looking at the
increase in at the edge.
32
The post-edge XAS response can be broken into two regions. The rst 30 eV correspond
to X-ray absorption near-edge structure (XANES) in which scattering processes and local
atomic resonances in X-ray absorption are the primary contributions to . The second
region of interest is the oscillations which occur after 30 eV from the absorption edge known
as extended X-ray absorption ne-structure (EXAFS). While the modern XAFS equation
relating the structure and is more complex, the essential concepts can be understood by
looking at the approximate form proposed by Sayers et al which related the XAFS signal
from the EXAFS region and the structure:
(k) =
∑
R
S
2
0
N
R
jf(k)j
kR
2
sin(2kR+2
c
+)e
2R=(k)
e
2
2
k
2
where R is the interatomic distance, N
R
; is the coordination number, is the temperature
dependent root mean square
uctuation in bond length, f(k) =jf(k)je
i(k)
is the backscat-
tering amplitude, (k) is the pase of the backscattering amplitude,
c
is the central-atom
partialwavephaseshiftofthenalstate,(k)istheenergydependentXAFSmeanfreepath,
S
0
is the amplitude factor and k is the photo electron wave vector k =
p
~!E
c
[77,78].
Here E
c
is the energy of corresponding edge.
Reverse Monte Carlo
ThereverseMonteCarlo(RMC)methodusesatomiccongurationstotexperimentaldata.
RMC combined with XAFS measurements has been used to determine the structure of a
wide variety of systems [47,105]. Unlike molecular dynamics and traditional Monte Carlo
simulations which require appropriately selected force elds to accurately simulate specic
33
Figure 3.4: The experimental apparatus at APS. The three ionization chambers are labeled
as I
0
;I
1
;I
2
. The
ouresence detector is labeled as I
f
. The sample is placed between I
0
and
I
1
.
materials, the RMC method minimizes the difference in measured structural terms between
experimental data and the simulated structure. The Metropolis MC algorithm is used to
minimize:

2
=
∑
i;j
(S
exp
S
sim
)
2
=
2
(3.6)
where S
exp
is the experimentally calculated term (usually the structure factor, scattering
intensity, pair distribution function, or EXAFS spectra), S
sim
is the term calculated from
the simulation, and is the accuracy of the experimental measurement [43].
Methodology
We performed X-ray absorption spectroscopy (XAS) experiments at Argonne National Lab-
oratory (ANL) using the advanced photon source (APS) (Figure 3.4). The main purpose of
34
the XAS experiments was to study the structure of SnO
2
and the liquid-solid interface of
colloidal SnO
2
nanoparticles. To that end, both SnO
2
pellets and SnO
2
colloid dispersed in
ultrapure water were analyzed. In the colloid samples, the pH was adjusted through hydro-
bromic acid (HBr) and rubidium hydroxide (RbOH). The synthesis of the powdered SnO
2
nanoclusters is out of the scope of this thesis but a detailed description can be found in the
Masters thesis of Thomas Winter (University of Duisburg, Duisburg Germany) [103].
The experiments were carried out using beamline 12-BM-B. The X-ray beam was ltered
using both Si(111) and Si(311) crystals which allowed for the selection of 4.5-30 keV and 10-
40 keV photons respectively. A double mirror setup was used to focus the beam to a 0.5 mm
x 1.0 mm area. A Canberra 13-element Ge detector was used as a
uorescence detector. In
addition, three ionization chambers were used to measure the transmission signal. For each
measurement, the same gas at atmospheric pressure was used in all three chambers. Argon
gas was used in the Sn-K edge measurements while nitrogen gas was used for the Br-K and
Rb-K edge measurements. The pellets were taped to the sample holders provided by ANL.
The solutions of SnO
2
in water were both placed in Kapton capillaries which were held by
a PTFE holder and lled inside liquid cells with Kapton windows mounted 45 degrees with
respect to the beam. Both capillaries and windows were made of Kapton HN foil with a
thickness of 25 m. Further details of the construction and design of the cells can be found
in the Appendix of Thomas Winters thesis.
A variable energy step size was used when scanning the energy range around the edges.
One second was spent on each energy value. The three regions were the pre-edge, edge and
post edge. From -200 to -20 eV (pre-edge) a step size of 5 eV was used. This was decreased
to 0.7954 (0.5) eV from -20 to +50 eV (edge) for the Sn-K (Br-K, Rb-K) edge. The step
35
size was increased again to 1.0 eV until 150 eV (edge) from the absorption edge. In the
nal segment, the step size was variable in energy units but constant in wave vector units
at 0:05

A
1
and ranged from 16 to 18

A
1
.
The number of measurements performed on each sample was determined on location
based on the quality of the measurements. Typically, pellets only required one measurement
whereas the colloids required up to 12 measurements for an adequate signal-to-noise ratio.
The xafsX software developed by Markus Winterer was used to calculate the absorption
coefficient from the raw data.
Results
X-ray diffraction data of the nanocrystalline SnO
2
powder were analyzed by Rietveld re-
nement. The corresponding crystallographic data are used to generate the initial atomic
congurations for RMC analysis. The Rietveld renement provided starting positions and
a structural factor which was then used as an input for RMC. The nal system positions
generated by RMC were then used to generate simulated EXAFS data to compare with the
experimental data. In addition, a 6 nm nanoparticle in vacuum and in water was simulated
in LAMMPS and used as the input for RMC. The SnO
2
and water force elds are described
in Chapter 2. Figure 3.2 shows the structure after the RMC analysis of the EXAFS data
for the periodic model, the cluster model and the LAMMPS nanocluster. The Metropolis
factor was lowered by a factor of 2 relative to the slab model for the cluster and LAMMPS
MD systems for the RMC analysis to improve convergence.
36
Figure 3.6 shows the EXAFs analysis of XAFS data for the SnO
2
pellet in the periodic
model (a,b) the cluster model (c,d) and the LAMMPs 6 nm particle in vacuum (e,f). Al-
though the results appear to be similar, there are noticeable differences at the largest peak
at 2.5

A and the peak at approximately 6

A. The average mean square error was calculated
as:
MSE =
1
N
∑
i
√
(y
data
i
y
model
i
)
2
; (3.7)
where y
model
i
corresponds to the model prediction at point i in phase/real space and y
data
i
corresponds to the data at point i. The results for the Fourier space and coordinate space
calculations respectively were 0.511 (slab model), 0.558 (cluster model), 0.617 (MD particle)

A
1
and0.399,0.448,0.491

Arespectively. Incomparingtheerrors,theslabmodelproduced
anoticeablybetterresult. Thismightbeduetothelargesizeoftheparticlewithbulkeffects
dominating surface effects in the data.
We also looked at the partial pair distribution functions (pPDFs). Figure 3.7 shows
the pPDFs resulting from the nal positions provided by the RMC analysis. The initial
positions from the Rietveld renement and initial LAMMPS conguration are also shown
(black dotted line). The slab and cluster model produce similar pPDfs however the cluster
model produces smoother oscillations especially in the regions between large peaks. The
LAMMPS model demonstrates signicant differences.
As mentioned above, the colloidal nanoparticle data was analyzed as well. Figure 3.8
shows the EXAFS analysis for the colloidal SnO
2
water system using the periodic model
(a,b), colloid model (c,d) and a LAMMPS simulation (e,f). The average mean square error
37
(a)
(b)
(c)
Figure3.5: AtomiccongurationofSnO
2
afterRMCanalysisofRietveldrenementofXRD
data for the periodic slab model (a), the cluster model (b) and LAMMPS simulation (c).
Unfortunately, after renement, all three models produced different results.
38
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(a)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(b)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(c)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(d)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(e)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(f)
Figure 3.6: EXAFS analysis by RMC using periodic boundary condition model in k-space
(a,c, e)and realspace (b,d,f). TheEXAFSdata(red)iscomparedwiththeEXAFSanalysis
from the RMC (blue). Their difference (green) is shown shifted by -10 (a,c,e) and -5 (b,d,f)
for clarity. Each row corresponds to a different model: the slab model (a,b), the cluster
model (c,d) and the LAMMPS generated particle (e,f).
39
(a) (b)
(c) (d)
(e) (f)
Figure 3.7: Partial pair distribution functions (pPDFs) resulting from RMC analysis using
periodic boundary model for Sn-O (a) and Sn-Sn (b). Included pPDFs from initial coordi-
nates (dotted line).
40
fortheFourierspaceandcoordinatespacecalculationswas0.640(slabmodel),0.700(cluster
model), 0.657(MDparticle)

A
1
and0.427, 0.490, 0.461

Arespectively. Similartothepellet
analysis,thebestttingmodelwastheslabmodel; however,therethedifferencesbetweenall
three were still not signicant. The corresponding pPDFs are shown in Figure 3.9. Relative
to the slab model, the colloid and MD systems are sharper with larger peaks. While the
two Rietveld models maintain a similar shape, the MD RMC rened system's pPDFs show
structural differences.
Comparison
The purpose of these experiments was to investigate the liquid-solid interface. The XANES
spectra and EXAFS analysis for the colloid and powder samples are shown in Figures 3.10.
There are small variations between the powder and colloidal EXAFS curves in both k space
andrealspace. ThedrypowdershowsamarginallylargerpeakintheXANEScurverelative
to the aqueous colloid. Comparing the pPDFs (Figure 3.11 we see a systematic reduction
in the peaks of both the SnO and Sn-Sn for the dry powder versus aqueous colloid. This
could be an artifact of analysis or a consequence of the improved relaxation at the particle's
surface due to the interaction with water. Unfortunately, the data is quite similar which
makesfurtherinvestigationintothesurface-liquidboundarydifficult. Goingforward,smaller
particles could be used to increase the surface area to volume ratio for future experiments.
We now directly compare the results from the three different models. The RMC analysis
in Figure 3.12 shows all three models converging to the same curve; however there is a
shift in the colloid systems Figure 3.12 (b) due to the different lattice constant in the MD
system. In comparing the pPDFs of the different systems, Figures 3.13, 3.14, there are
41
0 2 4 6 8 10 12 14 16
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(a)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(b)
0 2 4 6 8 10 12 14 16
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(c)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(d)
0 2 4 6 8 10 12 14 16
k (Å
−1
)
−10
−5
0
5
10
χ * k
3

data
rmc
diff
(e)
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
r (Å)
−5
0
5
10
15
(χ * k
3
)
data
rmc
diff
(f)
Figure 3.8: EXAFS analysis by RMC using periodic boundary condition model in k-space
(a) and real space (b). The EXAFS data (red) is compared with the EXAFS analysis from
the RMC (blue). Their difference (green) is shown shifted by -10 (a) and -5 (b) for clarity.
42
(a) (b)
(c) (d)
(e) (f)
Figure 3.9: Partial pair distribution functions (pPDFs) resulting from RMC analysis using
periodic boundary model for Sn-O (a) and Sn-Sn (b). Included pPDFs from initial coordi-
nates (dotted line).
43
(a) (b)
(c)
Figure 3.10: EXAFS data in fourier space (a) and real space (b) for the colloid and powder.
XANES spectra for the colloid and powder (c).
44
(a)
(b)
Figure3.11: Comparisonofpartialdistributionfunctionsfromclustermodelsfordrypowder
and aqueous colloid. The systematically larger amplitude in the case of the colloid is an
artifact of the applied normalization and cancels out in the RMC analysis
(a) (b)
Figure 3.12: Comparison of the RMC ts for the different models used for the pellet (a) and
the colloid (b). The shift in the MD model is due to the difference in lattice constant.
45
(a) (b)
Figure 3.13: Comparison of Sn-O pPDFs from RMC ts for the different models for the dry
powder (a) and the colloid (b). The systematically higher amplitude for the MD models in
(b) is due to the applied normalization (density). It cancels out in RMC analysis and is a
graphical problem.
(a) (b)
Figure 3.14: Comparison of the Sn-Sn pPDFs from RMC ts for the different models for
pellets (a) and colloid (b). Like the previous gure, the systematically higher amplitude in
the MD model is due to the applied normalization (density).
46
notable differences between the MD and experimental systems. The systematically larger
amplitude is an artifact of the applied normalization and cancels out in the RMC analysis.
It is therefore only a graphical problem.
Conclusion
The primary objective of this research was to gain a better understanding of the liquid-solid
interface through studying tin dioxide in water. Although high quality of data was collected
forbothsolidandcolloidalSnO
2
systems,thesimilaritybetweentheresultsmakesdiscerning
newinformationabouttheSnO
2
-waterinteractiondifficult. Infutureexperiments, thisissue
might be mitigated by using smaller nanoparticles as the higher surface to volume ratio
should produce data with more dramatic differences. Furthermore, more work can be done
on the molecular dynamics side. In the MD simulation, bond breaking/formation was not
allowed. The inclusion of bond-breaking through, for example, a REAXFF could improve
the accuracy of the surface-liquid interface representation. In addition, varying the pH of
the simulated water might further improve the simulation [58].
47
Chapter 4
A Neural Network for Principle of Least Action
Introduction
Least action
1
is the foundational principle of physics. It pervades classical mechanics [51],
special and general theory of relativity [52], quantum mechanics [23], electrodynamics [24],
and thermodynamics [27]. An articial neural network (ANN) is a network or circuit of
articial neurons or nodes that have collective computational properties [35]. In this paper,
we describe how an ANN learns the principle of least action by minimizing an action:
A =
∫
t
1
t
0
Ldt (4.1)
and provides the time evolution or dynamics of a many-body system between initial and
nal congurations at time t
0
and t
1
, respectively. In classical mechanics, the L in the action
integralistheLagrangian,itisthedifferencebetweentheinstantaneouskineticandpotential
energiesofasystemandtheclassicalactioncanbeviewedasthedifferencebetweenthetime-
averaged kinetic and potential energies. Action A is the key quantity in the loss function of
our neural network. The minimization of the loss function provides the trajectory q(t) and
its time derivative _ q(t) which collectively represent positions and velocities, respectively, of
particles in a system [51].
1
Chapter currently under review pending publishing
48
The principle of least of action embodies a boundary-value problem in which A is min-
imized with respect to q(t) subject to constraints that the initial and nal congurations,
q(t
0
) and q(t
1
), are xed [28]. Another formulation arising from the principle of least ac-
tion is Euler-Lagrange equations, which constitute an initial-value problem. These ordinary
differential equations are at the core of molecular dynamics (MD), the preeminent dynamic
simulation approach in physics, chemistry, materials science, and biology [73]. The essential
input to MD is interatomic potential energy from which forces are calculated and New-
tons equations or Hamiltons equations of motion are integrated over discretized time with
a nite-difference scheme [91]. The output is phase-space trajectoriesfq(t); _ q(t)g from
which structural, mechanical, thermodynamic, and dynamical properties of the system can
be computed and compared with experimental measurements.
While the neural-network approach to the principle of least action is new, the use of
machine learning (ML) in MD simulations is quite common these days especially in the
development of MD force elds [10]. ANN trained by data from quantum mechanical cal-
culations are rapidly replacing the painstakingly slow process of force-eld development.
ANN-based force elds are now widely used in MD simulation studies of material synthesis
andcharacterization. MLisalsomakingahugeimpactondataanalyticsrelatedtoMDsim-
ulations. Predictive models based on kernel ridge regression, support vector machines, and
random forests have been developed to predict material properties such as band gaps [75],
elastic constants [100], dielectric constants [76], and thermoelectric properties [38,101]. ML
combined with classical or quantum MD simulations have been used in a variety of other
ways in condensed matter physics, chemistry and materials science. For example, Bayesian
49
optimization methods together with MD have been used to discover optimal layered materi-
als for targeted material properties such as electronic band structure and thermal-transport
coefficients [9]. Generative models based on variational autoencoder (VAE), adversarial
networks, and normalizing
ows have been used to solve inverse problems of materials de-
sign from material properties. Two leading exemplars are the VAE models for the design
of drug molecules with targeted properties [54] and the study of structural transformation
pathways to design novel heterostructures coupling semiconducting (2H) and metallic (1T)
phases in strained two-dimensional MoWSe2 [74]. Another notable example of ML in MD is
the application of Restricted Boltzmann Machine (RBM) and its enhancement called Lim-
ited Boltzmann Machine (LBM) on a quantum annealer to model the synthesis of a MoS2
monolayer by chemical vapor deposition [34].
In this paper, we have taken the rst step in demonstrating that a feed-forward neural
network (NN) can learn the principle of least action and provide atomic trajectories that are
in excellent agreement with trajectories obtained by the conventional MD approach based
on a time-stepping algorithm for sequential integration of Newtons equations of motion. We
have chosen to train our NN on the Onsager-Machlup (OM) action:
S
OM
=
∫
T
0
[
N
∑
i=0

m
i
q
i
(t)+
@V(Q(t))
@⃗ q
i

2
]
dt (4.2)
where m
i
is the mass and

⃗ q
i
(t) is the second-order time derivative (i.e. acceleration) of the
position vector ⃗ q
i
(t) of ith particle. In Eq. 4.2, Q(t) collectively represents the coordinates
of all the particles in the system at time t and V(Q(t)) is the potential energy of the system.
Thisactionwasderivedforirreversiblestochasticthermodynamicprocesses[66]. Thegenesis
50
oftheOMactionwasOnsagersseminalpapers inwhichhedevelopedreciprocalrelationsfor
transport coefficients and proposed that the probability for
uctuation histories is exponen-
tially small with the exponent proportional to the time integral of dissipation function of the
history [64,65]. Later on, Onsager and Machlup showed that the most probable history is
given by the principle of least dissipation in which the dissipation function takes the form of
aLagrangianinclassicalmechanicsandtheminimizationofEq. 4.2wouldproducethemost
probable trajectory [66]. Here we show how a neural network minimizes S
OM
to produce
the entire trajectory of a classical many-body system in one fell swoop. Our NN for the
OM action can be easily implemented in any statistical ensemble and it is straightforward
to include constraints, invariances, and conservation laws in the neural network.
An intriguing aspect of using a neural-network approach to derive the dynamics of a
system from an action such as Eq. 4.2 is the general possibility of splicing the trajectory
into time segments and performing parallel calculations in time. This might be a way to
address the long-standing problem in atomistic simulations that they cannot reach long
timescales except for a few particular cases where the temporal locality of natural processes
andtransition-statetheoryallowthereformulationofthesequentiallong-timedynamicsasa
parallelsearchforlowactivation-barriertransitionevents. However, ingeneral, ithasbeena
challenge to simulate events involving slow atomic processes or slow conformational changes
such as protein folding [82] or dynamics of recrystallization in phase-change materials [48].
Space-time parallelization through domain decomposition combined with highly efficient
algorithms and massively parallel machines would enable large-scale atomistic simulations
of complex materials to reach long timescales [61].
51
We minimize the loss function in our NN model with the Nestrov-and-Adam optimizer
(NADAM) [20]. During training, the derivatives of the loss function with respect to pa-
rameters are calculated numerically using Googles JAX library. After training the network,
we compute atomic trajectories and compare them with the ground-truth MD simulations
performed using the same set of initial conditions. A full description of MD simulations for
the liquid Argon system is given in Appendix B.
Method
We train a neural network to learn the OM action (Eq. 4.2) for a many-body system of
Lennard-Jones (LJ) atoms [4]. The OM action has become a popular optimization approach
to solving atomistic boundary-value problems because it is guaranteed to have a global
minimum [106]. However, selecting a starting trajectory for the OM action is a non-trivial
task. One approach is to start with a straight-line trajectory between the initial and nal
congurations, express the trajectory in terms of a Fourier series [68], and vary the Fourier
coefficients to minimize the action in Eq. 4.2 under the constraint of energy conservation.
Unsurprisingly, the nal trajectory is found to depend heavily upon the initial values of
Fourier coefficients [68]. Lee et al. attempted to address this problem by introducing in
the cost function an additional term that makes the average kinetic energy consistent with
the temperature of the system [53]. Later on, they used conformational space annealing to
obtain many low OM action trajectories [53].
Fig. 4.1 shows our feedforward neural network which consists of an input node, a hid-
den layer with n units, and an output layer with d N units, where N is the number of
52
.
Hidden layer Output layer
Input layer
t
q
1, x
(t)
q
1, y
(t)
q
1, z
(t)
q
N, x
(t)
q
N, y
(t)
q
N, z
(t)
.
.
.
.
.
Figure4.1: Schematicoftheneuralnetwork. Thesystemcongurationforallspatialdegrees
offreedomiscalculatedbyinsertinganinstantoftimetintotheinputlayer. Asinglehidden
layer (blue) is used and the output layer (red) gives the atomic coordinates whose rst and
second derivatives provide velocity and acceleration, respectively. The blue and red lines
connecting layers represent the weights and each node also has a bias.
53
atoms. The neural network outputs the trajectories of all the atoms in the system, i.e.
cartesian components q
i;x
; q
i;y
;q
i;z
for i = 1 to N. We train the neural network to compute
conguration-space trajectories of Lennard-Jones systems consisting of N = 108, 256 and
500 atoms in three dimensions. Computational details of neural-network training and LJ
systems are given in Appendix B.
We evaluate the OM action numerically by discretizing the time integral in Eq. 4.2
into N
t
grid points. The time increment is ∆t = T=Nt, where T = 1ps and N
t
= 25 (see
Appendix B for results with different values of N
t
). The loss function L for the neural
network includes not only the OM action but also constraints due to energy conservation
and boundary conditions:
L(w;b;) =
1
S
dis
OM
+
2
C
pos
+
3
C
vel
+
4
T
∑
t=0
[E(Q(t
;
_
Q(t))E
0
]
2
(4.3)
C
pos
= (Q(0)Q
0
)
2
+(Q(T)Q
T
)
2
(4.4)
C
vel
= (
_
Q(0)V
0
)
2
+(
_
Q(T)V
T
)
2
(4.5)
S
dis
OM
=
T
∑
t=0

Q+
dV(Q)
dQ

2
∆t (4.6)
E(Q;
_
Q) =
N
∑
i=1
1
2
m
i
_
⃗ q
2
i
+
1
2
N
∑
i;j̸=i
4ϵ
(

12
j⃗ q
i
⃗ q
j
j
1
2

6
j⃗ q
i
⃗ q
j
j
6
)
(4.7)
wherefw; bg represents all the weights and biases in the neural network. For simplicity,
we use Q(t) to denote Q(fw; bg; t). Q
0
, Q
T
, V
0
, V
T
collectively represent initial and nal
coordinates and velocities of all the atoms in the system and s are the hyper-parameters.
Note that the last term in Eq. 4.4 is the energy-conservation constraint, where E
0
is the
initial total energy of the system. In Eq. 4.7, E
(
Q(t);
_
Q(t)
)
is the total energy with " and
54
(a) (b)
(c) (d)
Figure 4.2: Four selected NN (solid blue) and MD (solid red) trajectories show that the NN
can reproduce non-trivial paths. These are for a time period of 1 ps
as parameters of the LJ potential. The parameters are chosen with care to ensure that
all the boundary conditions are well satised by the neural network. We nd that setting

1
=
4
= (max
t
(jV (Q(t))j))
1
ensures that the neural network sets the initial conditions
correctly before minimizing the action.
Results
Our neural network learns to compute atomic trajectories of two- and three-dimensional
LJ systems from the OM least action principle. Here we present results for a 3D system
containing500atoms, whichwastrainedonaneuralnetworkwith125neuronsinthehidden
layer. Neural network results for 108 and 256-atom LJ systems are presented in SM. Results
55
(a) (b) (c)
Figure 4.3: Comparison between the neural network (NN) and MD simulation. (a) Kinetic
(blue), potential (green) and total (red) energy per particle as a function of time for MD
(square dots) and NN (solid lines). Radial distribution function and velocity autocorrelation
function for MD (blue) and NN (red) are shown in (b) and (c), respectively.
(a) (b) (c)
Figure 4.4: (a) shows how the three components of the loss function decrease with the
number of epochs. The boundary condition term (red) and energy conservation term (blue)
fall more rapidly than the OM action term (green). (b) Shows the mean square error (MSE)
for trajectory (red) and momentum (blue) which was calculated by comparing the MD
trajectory to the neural network predicted trajectory. (c) shows differences between MD
and NN trajectories (red) and velocities (blue) as a function of time. Here = 2ps.
(a) (b) (c)
Figure 4.5: Scaling of NN for 108-particle system. (a), (b) and (c) are the run time for every
1000 training epoch as a function of the number of units in the hidden layer, the size of the
time grid, and the number of atoms in the system, respectively.
56
in this section are generated on Nt = 25 time increments which correspond to a time step
( 40 fs) that is 20 times the MD timestep (2 fs).
Figure 4.2 shows a comparison of four randomly chosen single-particle trajectories in the
system over a time period of 1 ps. These trajectories cover a time scale beyond the straight-
line ballistic motion and possess sharp changes in direction. It is remarkable that the neural
network (NN) is able to reproduce the dynamics at the individual particle level. From this
level of matching the MD results for particle trajectories, it is not surprising that the NN
results for all other quantities are in excellent agreement with the ground truth MD results.
InFigure 4.3, wecomparethe results fromthe NN andMD for some standardsimulation
quantities. Since we use the microcanonical ensemble, the total energy should be constant.
This condition is part of the loss function, and thus we expect it to hold quite well. Indeed,
Figure 4.3a shows that the total energy is constant, and that for the kinetic and potential
terms the NN matches the MD values. The structure of the liquid is characterized by the
radial distribution function [45] g(r) shown in Fig. 4.3b. The NN matches the MD result.
Slightly below r = 1, g(r) = 0 because the particles cannot overlap. The NN trajectory
maintains this feature as well. The peaks in g(r) corresponding to the minimum energy
separation for nearest and second nearest neighbors indicate that the NN is able to obtain
thecorrectliquidstructure. Tocharacterizetheparticledynamics,wecalculatedthevelocity
autocorrelation function [3] (VAF) shown in Fig. 4.3c. The match between NN and MD is
good, and the correlation time tc, when the VAF = 0 for the rst time, is about 0.16 t
c
.
Beyond t
c
the velocities and, in general, the particle dynamics are not correlated with the
initial dynamics of the particles. This quantity is important in calculating time averages, as
57
a time series of data with points separated by at least t
c
are not correlated and contribute
independently to the average quantity
Figure4.4showshowvariouscontributionstothelossfunctiondecreasewiththenumber
of epochs. Since the system is in the microcanonical ensemble, the constraint of energy
conservation is part of the total loss function. All three terms in the loss function decrease
rapidly in the rst few thousand epochs. The contribution from the OM action is the largest
and determines the minimal value of the total loss. The loss function corresponding to
boundary-conditionconstraintsdecreasesrapidlyandreachesbelow10
6
. Itisnotsurprising
that this loss is so low because these conditions are only on two points in the trajectory.
Contrast this with the loss for the energy-conservation constraint, which goes below 10
5
after 200,000 epochs since the constraint must be satised for the whole trajectory.
We also show in Fig. 4.4 the mean-squared error (MSE) of the position (Q) and momen-
tum (P) with respect to the corresponding MD variables. These errors are computed from
trajectories generated by the neural network and the ground-truth MD simulations. The
general shape of the MSE is similar for the two quantities although they differ by orders of
magnitude. We offer two explanations for the rising error per time derivative. First, refer-
ringto equations 4.4 - 4.7, the time derivativesare described with fewerparameters than the
positions as the biases outside the activation function are killed by the derivative. Second,
all quantities in the loss function depend explicitly on position, two depend on the velocity
andonlyonedependsontheacceleration. Sincewehaveanitenumberofparameters, they
can only approximate the function up to some error. These errors magnify as we take higher
derivatives. The mean square error for the positions is very small, which as we saw in Fig.
4.2, results in excellent matching of the neural network and MD trajectories. The behavior
58
of the momentum as a function of epoch is similar to the positions, but with a larger MSE
at convergence.
We also studied how the NN scales with the number of units in the hidden layer and
discretizationofthetimedomainonIntel(R)Xeon(R)CPUE5-2640machinewith2.60GHz
clock speed. Surprisingly, with the help of Google JAX, Fig.4.5(a) indicates that increasing
the number of neurons in the hidden layer doesnt affect the run time appreciably. The
runtimeremainsaround210secondsfor1,000epochsoftraining. However,changingthegrid
size of the time domain increases the computing time signicantly, as shown in Fig. 4.5(b).
Figure 4.5(c) shows how the run time scales with the system size. The behavior appears
quadratic because the run time is dominated by the force calculation. Its straightforward
to reduce the scaling of the run time with the system size to O(N) by using the linked-list
method [4].
Conclusion
We have shown that a feedforward neural network can correctly generate atomic trajecto-
ries of a Lennard-Jones system by minimizing the Onsager-Machlup action. Besides the
phase-space trajectories, the potential and kinetic energies, radial distribution and velocity
autocorrelation functions are in excellent agreement with the corresponding results from the
ground-truth MD simulation. Our neural network approach is generally applicable to sys-
tems with more complex interatomic interactions than the Lennard-Jones potential. The
neural-network approach to principle of least action has a few advantages over the conven-
tional MD method:
59
The neural network provides the entire trajectory in one fell swoop over a time period
of 1 ps.
The time step in the neural-network approach to action is a factor of 20 larger than
the time step in conventional MD simulation.
Evenwithamuchlargertimestep,theenergyconservationovertheentiretimedomain
is as good as in the conventional MD simulation.
Goingforward,theneural-networkapproachpresentedherehasthepotentialtobeagame
changer in more than one way. For example, Hills et al. have developed an algorithm based
on the least action principle to predict the dynamics of physical systems using observed
data [33]. Our approach combined with theirs can model observational data for different
types of physical systems. This can be very useful in predicting the future dynamics of
physical systems outside the range of measurements and quite conceivably the two neural
networks could help scientists discover new phenomena hidden in data.
Another important area in which our neural-network approach can be impactful is impu-
tation of missing data. Suppose we have disconnected trajectories of a physical system and
we want to nd the missing piece between the end point of the rst and the starting point
of the second trajectory. Using the algorithm of Hills et al., one can nd the Lagrangian of
the underlying system. Subsequently, one can use our neural network to interpolate the end
point of the rst and the starting point of the second trajectory.
Lastly and most importantly, the neural-network approach to principle of least action
has the potential to parallelize the time axis in atomistic simulations in a way similar to the
60
divide-and-conquerschemesusedinspatialdecompositionoflarge-scaleMDsimulations[83].
These schemes have made it possible to simulate multimillion-to-multibillion atom systems
[60]. In the time parallelization approach, one would split a long-time trajectory, say 1ns
long, into shorter time segments on the order of 1 ps, then use the Hamiltonian Monte
Carlo approach [17] accelerated by the Boltzmann generator to generate a Markov chain of
1,000 intermediate congurations between t = 0 and t = 1 ns. Subsequently, one could use
the least action neural network to run 1,000 simulations concurrently on 1,000 processors
and generate a continuous trajectory over a time period of 1 ns. Thus, the neural-network
approach would provide substantial gains in parallel efficiency.
61
Chapter 5
Solving IVPs with Least Action Neural Networks
Introduction
Evenwithmoderncomputingcapabilities, traditionalmoleculardynamics(MD)simulations
are limited to hundreds of nanoseconds for large systems. This limitation is due, in part, to
the sequential nature of MD calculations. Each timestep must be calculated in order. This
is compounded by many systems requiring a timestep on the order of 1 fs. As discussed
in the previous chapter, action derived methods have been used to run simulations with an
effective time-step orders of magnitude higher than allowable in typical molecular dynamics
simulations. In Chapter 4 we showed how an Onsager-Machlup (OM) action trained neural
network has been used to nd the trajectories between two points in conguration space. In
this chapter, we demonstrate that the same neural network can solve the IVP with minor
alterations of the loss function.
Method
Wetrainedtheneuralnetworktolearnthespatialtrajectoriesfora108bodyLennard-Jones
atomic system. The LJ parameters for liquid argon were implemented. Our feedforward
network is comprised of an input node, an n unit hidden layer and an output layer with dxN
units where N is the number of atoms and d is the dimensionality of the system (Figure 4.1.
62
A scalar value for time is used as the input and the outputs are the corresponding cartesian
coordinates (and time derivatives) of each atom at that time.
The OM action is numerically evaluated by splitting the time grid into N
t
gridpoints.
The time difference here corresponds to the usual timestep used in conventional MD. The
loss function for the network is given by:
L(w;b;) =
1
S
dis
OM
+
2
(Q(0)Q
0
)
2
+
2
(
_
Q(0)V
0
)
2
+
3
T
∑
t=0
[E(Q(t
;
_
Q(t))E
0
]
2
(5.1)
S
dis
OM
=
T
∑
t=0

Q+
dV(Q)
dQ

2
∆t (5.2)
wherefw; bg represents all the weights and biases in the neural network. As in the previous
chapter, we use Q(t) to denote Q(fw; bg; t). The initial conditions for the system are
collectively labeled Q
0
;V
0
respectively. The 's are hyperparameters which allow us to
increase the probability that the system learns the correct boundary conditions. The term
E is given by Equation 4.7 and E
0
is the total energy in the system. We found that a value
ofjmax
t
(V(Q(t))j for
1
and
3
produced trajectories with accurate initial conditions. Note
that the potential energy term within the factors was treated as a constant with respect
to loss gradient calculations.
The loss of the system is minimized using the Nestrov-and-Adam optimizer (NADAM)
[20]. The derivatives of the loss function were calculated using Google's JAX library [13].
The network predicted trajectories were compared with conventional molecular dynamics
starting from identical initial conditions. The full description of the conventional MD can
be found in Appendix B.
63
0 1 2 3 4 5
Epoch
1e5
10
−7
10
−6
10
−5
10
−4
10
−3
Loss
Loss
BC
H
LA
(a)
0 100000 200000 300000 400000 500000
Epoch
10
−4
10
−3
10
−2
10
−1
10
0
10
1
10
2
MSE
Q
P
A
(b)
Figure 5.1: (a) The boundary condition term (orange), total energy term (green), Onsager-
Machlup term (red) and total loss (blue) as a function of training epoch. Note that these
values include the corresponding
i
terms. After the initial few hundred epochs, the loss
is dominated by the OM action term. (b) The mean square error in the position (blue)
momentum (orange) and acceleration (green). The slopes tend towards zero implying that
further epochs would not produce signicant improvements in the trajectory.
Results
We looked at the quality of the solution as both a function of epoch, (Figure 5.1) and as a
function of time (Figure 5.2). Much like the BVP case, the leading term in the loss function
is the OM-action term. After an initial increase which is primarily due to rapid
uctuations
of the 1=jPE
max
j component during the initial hundred epochs, the loss experiences a rapid
decay in all of its components followed by an increasingly slow continued decline. The loss is
dominated by the OM-action term. The energy conservation term reaches a minimum and
then proceeds to rise back up; however, after the initial descent it remains multiple orders
of magnitude lower than the OM term. The term corresponding to the initial conditions
oscillates about around 9*10
8
for the last 10000 epochs. Figure 5.1 b shows the MSE in
position (blue), velocity (orange) and acceleration (green) as a function of epoch. All three
terms continue to improve with epochs however, after 300,000 epochs the rate of reduction
becomes much lower.
64
0.00 0.05 0.10 0.15 0.20 0.25
Time (τ)
−5
−4
−3
−2
−1
0
1
Energy (ε )
MD: Etot(t)
MD: Epot(t)
MD: Ekin(t)
NN: Etot(t)
NN: Epot(t)
NN: Ekin(t)
(a)
0.00 0.05 0.10 0.15 0.20 0.25
Time (τ)
10
−8
10
−6
10
−4
10
−2
10
0
10
2
MSE
Q
P
A
(b)
Figure 5.2: (a) The kinetic (blue), potential (green) and total energies (red) as a function
of time for conventional molecualar dynamics (dashed lines) and our NNMD (solid lines)
after training. (b) The MSE in acceleration (green), momentum (orange) and position (Q)
comparing the NN to the conventional MD simulation.
0.0 0.5 1.0 1.5 2.0 2.5
r (σ)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
g(r)
MD g(r)
NN g(r)
(a)
0.00 0.05 0.10 0.15 0.20 0.25
time (τ)
−0.2
0.0
0.2
0.4
0.6
0.8
1.0
VAF
MD VAF
NN VAF
(b)
Figure 5.3: The radial distribution function (a) and velocity autocorrelation function (b) for
regular MD (red) and NNMD after 500,000 epochs (blue).
65
x(σ)
4.18
4.21
4.24
4.27
y(σ)
3.26
3.29
3.32
3.35
z(σ)
1.91
1.94
1.97
2.01
x(σ)
3.11
3.15
3.19
3.23
y(σ)
2.56
2.62
2.67
2.72
z(σ)
2.08
2.15
2.22
2.30
x(σ)
4.08
4.10
4.12
4.14
y(σ)
1.21
1.28
1.36
1.44
z(σ)
3.37
3.39
3.41
3.43
x(σ)
3.31
3.33
3.34
3.36
y(σ)
4.47
4.50
4.53
4.57
z(σ)
4.92
4.94
4.95
4.96
Figure 5.4: Particle trajectories for 4 selected atoms of both MD (red-dashed) and NNMD
(blue solid). Both trajectories agree at t=0; however, some particle trajectories diverge
signicantly towards the end of the time span.
66
Figure 5.2 shows how the system performs as a function of time. Across the entire tra-
jectory, energy is conservation is observed. Comparing the conventional MD and NNMD
potential and kinetic energies, we see there is strong agreement throughout the entire times-
pan. The mean square error of the systems trajectory and velocity rise slowly with time
but start below 10
8
. The acceleration's MSE, much like the case of the BVP, starts with
a high value although its rate of increase with respect to time is small relative to its total
magnitude. The MSE increases for all three quantities with time which can be explained by
compounding differences in the trajectories. Small difference in position and velocity will
propagate into larger differences farther down the time line. As can be seen in Figure 5.4,
despite the energy agreement individual trajectories vary signicantly from the correct solu-
tionas time is increased. Arepresentativesampling of higherand loweragreeing trajectories
were selected for the gure. Despite these errors, the network was able to reproduce the
radial distribution function (Figure 5.3 a) and the velocity autocorrelation function (Figure
5.3 b) with a high degree of delity.
Conclusion
The endgame of this research is time-parallel molecular dynamics and access to larger
timescales. While this work represents a rst step towards these goal, there are several
issues which have yet to be addressed.
How does one address spatial decomposition? In conventional MD, spatial decomposi-
tions is handled by either verlet tables or linked lists both of which update periodically [4].
Although the OM action method permits a lager time step; the frequency of rebuilding
67
neighborlistsisdependentonthesystemunderstudyandwouldlikelyneedtobedoneafter
thesameamountofsimulationtimeregardlessofmethodused. Intheinitialhundredepochs
of training, particle trajectories can change signicantly from epoch to epoch requiring the
frequent rebuilding of neighborlists. As training progresses one would expect the number of
neighbor lists to be rebuilt to be reduced. Alternatively, one can restrict the timespan the
network is trained on based such that atoms do not move too far from their initial positions
removing the need to perform spatial decomposition more than once.
In order for NNMD to become a useful tool in solving BVP or IVPs, the time to train
the network must be reduced. The current implementation of the NNMD code has not been
optimized for speed, so it is likely that revisions of the code will produce some efficiency
improvements. From a theoretical perspective, training time can be reduced by reducing
the number of epochs required to reach an acceptable accuracy or by reducing the time per
epoch. While code efficiency improvements address, we should implement parallel training.
Additionally, other action methods [53,68] have begun their training with fewer time points
and gradually increased the number of points until the desired accuracy was reached. Algo-
rithms to train feed-forward neural networks in parallel using both data parallelization and
model parallelization have been developed [81] and should be incorporated into our system.
How can we improve the accuracy of the results? Although systemwide parameters such
as the correlation function, energy and velocity autocorrelation function agree, the particle
trajectoriesareclearlyquitedifferentevenaftertrainingappearstosaturate. Currentlyeach
degreeoffreedomshares2nparameters. Differentnetworkstructureswithmoreindependent
terms per degree of freedom could improve the results. As noted above and in the previous
chapter, the initial network parameters can have a signicant effect on the nal trajectory.
68
The NADAM alogrithm usually produces good results. In our work in both the BVP and
IVP systems, we found that in every case tested it trained faster in terms of loss per epoch
than Nestrov assisted gradient decent (NAG) however, is a local search which comes with
certain limitations. Using a global search algorithm, such as a combined genetic algorithm
and gradient decent [49], could improve the accuracy per epoch as well as the nal satura-
tion accuracy at the cost of computation time. the NADAM training algorithm is a local
optimizer.
69
Chapter 6
Conclusion
Thecomputationalstudyofatomicmovementremainsanessentialtoolinbiology, chem-
istry,materialscienceandphysics. Simulationsenhanceexperimentsbyshowinginformation
which would be difficult or prohibitively expensive to observe directly through experiment
and provides direction for new experiments by suggesting likely productive systems to in-
vestigate further. In response to the global pandemic of COVID19, penta
ops of computing
resources have been devoted to simulating the virus and related proteins [18,36,37]. Large
scale quantum mechanical simulations of the virus and the top 1000 potential inhibitors
are being performed to aid in drug development. Novel potential cure molecules have been
identied through reinforcement learning combined with these simulations. MD simulations
are allowing investigation into why certain individuals develop more severe symptoms.
The molecular dynamics simulations of SnO
2
nanoparticles described in Chapter 2 are
not as urgent as a response to COVID19 however, nanocavities have a range of applications
fromindustrialcleaningtomedicaldrugdelivery. Inrunningthesesimulationswediscovered
the formation of void-cavities in polar liquids subject to homogeneous external elds. The
simulation has become the impetus for future experimental work. The time to cavity forma-
tion follows the stretched exponential response. Characterizing the time dependence as such
supports the theory that the cavities form from the diffusion of local density
uctuations
into a sink region (the void). The rapid collapse of the cavity upon electric eld removal
70
could be used in breaking apart larger nanoclusters. The appearance of the voids with both
polar and non-polar force elds substantiates this discovery. The timescales of the cavity
formation required the use of extremely large external elds which, due to the absence of
bond breaking in the force eld, leaves room for better simulation methodology in future
simulations.
In Chapter 3, I deviated from the overlaying theme of computational investigation and
performed physical experiments at Argonne National Laboratory. The purpose was to inves-
tigate the solid-liquid interface using tin dioxide. High quality data was collected, analyzed
using multiple Reiveld renement models and compared to MD simulation. There were not
signicant differences between the solid and colloidal SnO
2
nanoparticle systems data which
made detailed probing of the solid-liquid interface difficult. The data suggests that smaller
particles might be required to properly investigate this interaction using XAS.
Chapter 4 details the development of the action-based neural network solution to the
atomic boundary value problem. Using our new method, we achieved a timestep 40 times
larger than that of conventional molecular dynamics. There are several promising direc-
tions to continue this research which were unable to be pursued by the author due to time
constraints as this project was started late in his PhD career. Taking cues from other
action-methods, implementing a global search algorithm could ameliorate the dependence
of the nal trajectory on the initialization procedure. Furthermore, to better estimate the
relative success of our method, we should compare the nal OM action of our system with
that of other action-methods on identical systems. The problem solved by our network is
over dened as we supply the initial positions and velocities to the network at both the
initial and nal time point. We performed limited testing in specifying only the positions
71
and system energy; however, we were able to generate a low OM trajectory (albeit not the
trajectory generated by MD). Reducing the requirement of providing initial/nal velocities
would make our network applicable to a wider range of systems. Currently, our network has
been unable to reproduce trajectories longer than 2 ps for any of the tested system sizes and
it is unclear as to what is causing this limitation. We have not yet demonstrated spatial
decomposition which is required for any practically relevant molecular dynamics simulation.
Chapter 5 shows that the same network can be applied to solve the IVP. The accuracy of
the results leaves something to be desired; however, it remains an exciting rst step towards
the development of time-parallel molecular dynamics.
72
Reference List
[1] Jolyon Aarons, Misbah Sarwar, David Thompsett, and Chris-Kriton Skylaris. Per-
spective: Methods for large-scale density functional calculations on metallic systems.
The Journal of chemical physics, 145(22):220901, 2016.
[2] AshutoshAgarwal, WunJernNg, andYuLiu. Principleandapplicationsofmicrobub-
ble and nanobubble technology for water treatment. Chemosphere, 84(9):1175{1180,
2011.
[3] BJAlderandTEWainwright. Decayofthevelocityautocorrelationfunction. Physical
review A, 1(1):18, 1970.
[4] Michael P Allen and Dominic J Tildesley. Computer simulation of liquids. Oxford
university press, 2017.
[5] Mouna El Mekki Azouzi, Claire Ramboz, Jean-Fran cois Lenain, and Fr ed eric Caupin.
A coherent picture of water at extreme negative pressure. Nature Physics, 9(1):38,
2013.
[6] DowonBae, BrianSeger, PeterCKVesborg, OleHansen, andIbChorkendorff. Strate-
giesforstablewatersplittingviaprotectedphotoelectrodes. Chemical Society Reviews,
46(7):1933{1954, 2017.
[7] AV Bandura, JO Sofo, and JD Kubicki. Derivation of force eld parameters for sno2-
h2osurfacesystemsfromplane-wavedensityfunctionaltheorycalculations. The Jour-
nal of Physical Chemistry B, 110(16):8386{8397, 2006.
[8] Allen J Bard, Larry R Faulkner, Johna Leddy, and Cynthia G Zoski. Electrochemical
methods: fundamentals and applications, volume 2. wiley New York, 1980.
[9] Lindsay Bassman, Pankaj Rajak, Rajiv K Kalia, Aiichiro Nakano, Fei Sha, Jifeng
Sun, David J Singh, Muratahan Aykol, Patrick Huck, Kristin Persson, et al. Active
learning for accelerated design of layered materials. npj Computational Materials,
4(1):1{9, 2018.
[10] J org Behler. Perspective: Machine learning potentials for atomistic simulations. The
Journal of chemical physics, 145(17):170901, 2016.
[11] Vladimir P Berishvili, Valentin O Perkin, Andrew E Voronkov, Eugene V Radchenko,
Riyaz Syed, Chittireddy Venkata Ramana Reddy, Viness Pillay, Pradeep Kumar,
Yahya E Choonara, Ahmed Kamal, et al. Time-domain analysis of molecular dy-
namics trajectories using deep neural networks: Application to activity ranking of
tankyrase inhibitors. Journal of chemical information and modeling, 59(8):3519{3532,
2019.
73
[12] Oleg Borodin and Grant D Smith. Quantum chemistry and molecular dynamics simu-
lation study of dimethyl carbonate: ethylene carbonate electrolytes doped with lipf6.
The Journal of Physical Chemistry B, 113(6):1763{1776, 2009.
[13] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary,
Dougal Maclaurin, and Skye Wanderman-Milne. JAX: composable transformations of
Python+NumPy programs, 2018.
[14] MatthewABrown,ZareenAbbas,ArminKleibert,RichardGGreen,AlokGoel,Sylvio
May, and Todd M Squires. Determination of surface potential and electrical double-
layer structure at the aqueous electrolyte-nanoparticle interface. Physical Review X,
6(1):011007, 2016.
[15] CarlCalemanandDavidVanDerSpoel. Picosecondmeltingoficebyaninfraredlaser
pulse: A simulation study. Angewandte Chemie International Edition, 47(8):1417{
1420, 2008.
[16] Fr ed eric Caupin and Eric Herbert. Cavitation in water: a review. Comptes Rendus
Physique, 7(9-10):1000{1017, 2006.
[17] Khue-DungDang,MatiasQuiroz,RobertKohn,Minh-NgocTran,andMattiasVillani.
Hamiltonian Monte Carlo with energy conserving subsampling. MIT Press, 2019.
[18] Cesar de la Fuente. Effects of genetic variants of human ace2 protein on sars-cov-2
spike protein binding, May 2020.
[19] Vesselin Dimitrov and Takayuki Komatsu. Classication of simple oxides: a polariz-
ability approach. Journal of Solid State Chemistry, 163(1):100{112, 2002.
[20] Timothy Dozat. Incorporating nesterov momentum into adam. 2016.
[21] Niall J English and JMD MacElroy. Hydrogen bonding and molecular mobility in
liquid water in external electromagnetic elds. The Journal of chemical physics,
119(22):11806{11813, 2003.
[22] Niall J English and Conor J Waldron. Perspectives on external electric elds in molec-
ular simulation: progress, prospects and challenges. Physical Chemistry Chemical
Physics, 17(19):12407{12440, 2015.
[23] RichardPFeynman. Theprincipleofleastactioninquantummechanics. InFeynman's
ThesisA New Approach To Quantum Theory, pages 1{69. World Scientic, 2005.
[24] Richard P Feynman, Albert R Hibbs, and Daniel F Styer. Quantum mechanics and
path integrals. Courier Corporation, 2010.
[25] PW Fowler, JH Harding, and NC Pyper. The polarizabilities and dispersion coeffi-
cients for ions in the solid group iv oxides. Journal of Physics: Condensed Matter,
6(48):10593, 1994.
74
[26] Yongchun Fu and Alexander V Rudnev. Scanning probe microscopy of an elec-
trode/ionic liquid interface. Current Opinion in Electrochemistry, 1(1):59{65, 2017.
[27] Vladimir Garc a-Morales, Julio Pellicer, and Jos e A Manzanares. Thermodynamics
based on the principle of least abbreviated action: Entropy production in a network
of coupled oscillators. Annals of Physics, 323(8):1844{1858, 2008.
[28] Herbert Goldstein, Charles Poole, and John Safko. Classical mechanics, 2002.
[29] David C Grahame. The electrical double layer and the theory of electrocapillarity.
Chemical reviews, 41(3):441{501, 1947.
[30] CGGrayandEdwinFTaylor. Whenactionisnotleast. American Journal of Physics,
75(5):434{458, 2007.
[31] JL Green, DJ Durben, GH Wolf, and CA Angell. Water and solutions at negative
pressure: Raman spectroscopic study to-80 megapascals. Science, 249(4969):649{652,
1990.
[32] John M Griffin, Alexander C Forse, Wan-Yu Tsai, Pierre-Louis Taberna, Patrice Si-
mon, and Clare P Grey. In situ nmr and electrochemical quartz crystal microbalance
techniquesrevealthestructureoftheelectricaldoublelayerinsupercapacitors. Nature
materials, 14(8):812, 2015.
[33] Daniel JA Hills, Adrian M Gr utter, and Jonathan J Hudson. An algorithm for discov-
ering lagrangians automatically from data. PeerJ Computer Science, 1:e31, 2015.
[34] Sungwook Hong, Ken-ichi Nomura, Aravind Krishnamoorthy, Pankaj Rajak, Chun-
yang Sheng, Rajiv K Kalia, Aiichiro Nakano, and Priya Vashishta. Defect healing in
layeredmaterials: Amachinelearning-assistedcharacterizationofmos2crystalphases.
The journal of physical chemistry letters, 10(11):2739{2744, 2019.
[35] JJ Hopeld. Human memory, error correction codes and spin glasses. Proc. Natl.
Acad. Sci., 79:2554, 1982.
[36] Stratos Davlos Innoplexus. 19 hpc covid-19 novel molecule generation with reinforced
learning, May 2020.
[37] Stephan Irle. Quantum mechanics-based renement of sars-cov-2 inhibitors from clas-
sical docking, May 2020.
[38] YumaIwasaki,IchiroTakeuchi,ValentinStanev,AaronGiladKusne,MasahikoIshida,
Akihiro Kirihara, Kazuki Ihara, Ryohto Sawada, Koichi Terashima, Hiroko Someya,
et al. Machine-learning guided discovery of a new thermoelectric material. Scientic
reports, 9(1):1{7, 2019.
[39] Shane Jackson, Aiichiro Nakano, Priya Vashishta, and Rajiv K Kalia. Electrostrictive
cavitation in water induced by a sno2 nanoparticle. ACS omega, 2019.
75
[40] Hao Jiang, Othonas A Moultos, Ioannis G Economou, and Athanassios Z Pana-
giotopoulos. Hydrogen-bonding polarizable intermolecular potential model for water.
The Journal of Physical Chemistry B, 120(48):12358{12370, 2016.
[41] Jin-Wu Jiang, Jian-Sheng Wang, and Baowen Li. Youngs modulus of graphene: a
molecular dynamics study. Physical Review B, 80(11):113405, 2009.
[42] Martin Karplus and Gregory A Petsko. Molecular dynamics simulations in biology.
Nature, 347(6294):631{639, 1990.
[43] DA Keen and RL McGreevy. Structural modelling of glasses using reverse monte carlo
simulation. Nature, 344(6265):423{425, 1990.
[44] Jan Kern, Roberto Alonso-Mori, Rosalie Tran, Johan Hattne, Richard J Gildea,
Nathaniel Echols, Carina Gl ockner, Julia Hellmich, Hartawan Laksmono, Raymond G
Sierra, et al. Simultaneous femtosecond x-ray spectroscopy and diffraction of photo-
system ii at room temperature. Science, 340(6131):491{495, 2013.
[45] John G Kirkwood and Elizabeth Monroe Boggs. The radial distribution function in
liquids. The Journal of Chemical Physics, 10(6):394{402, 1942.
[46] J Klingler. Rj hunter: Introduction to modern colloid science. Berichte der Bunsen-
gesellschaft fur physikalische Chemie, 99(3):591{592, 1995.
[47] Alexander Kompch, Ayaskanta Sahu, Christian Notthoff, Florian Ott, David J Nor-
ris, and Markus Winterer. Localization of ag dopant atoms in cdse nanocrystals by
reverse monte carlo analysis of exafs spectra. The Journal of Physical Chemistry C,
119(32):18762{18772, 2015.
[48] Christophe Krzeminski, Quentin Brulin, V Cuny, Emmanuel Lecat, Evelyne Lampin,
and Fabrizio Cleri. Molecular dynamics simulation of the recrystallization of amor-
phous si layers: Comprehensive study of the dependence of the recrystallization veloc-
ity on the interatomic potential. Journal of applied physics, 101(12):123506, 2007.
[49] Somesh Kumar, Manu Pratap Singh, Rajkumar Goel, and Rajesh Lavania. Hybrid
evolutionary techniques in feed forward neural network with distributed error for clas-
sication of handwritten hindi swars. Connection Science, 25(4):197{215, 2013.
[50] Guillaume Lamoureux and Benot Roux. Modeling induced polarization with classical
drude oscillators: Theory and molecular dynamics simulation algorithm. The Journal
of chemical physics, 119(6):3025{3039, 2003.
[51] Lev Davidovich Landau and Evgenii Mikhailovich Lifshitz. Course of theoretical
physics, volume 1. Elsevier, 2013.
[52] Lev Davidovich Landau and Evgenii Mikhailovich Lifshitz. Course of theoretical
physics, volume 2. Elsevier, 2013.
76
[53] In-Ho Lee, Sukky Jun, Hanchul Kim, Seung-Yeon Kim, and Jooyoung Lee. Exploring
dynamic pathways by action-derived molecular dynamics. International journal of
nanotechnology, 3(2-3):334{352, 2006.
[54] Jaechang Lim, Seongok Ryu, Jin Woo Kim, and Woo Youn Kim. Molecular genera-
tive model based on conditional variational autoencoder for de novo molecular design.
Journal of cheminformatics, 10(1):1{9, 2018.
[55] Jun Liu, Yangyang Gao, Dapeng Cao, Liqun Zhang, and Zhanhu Guo. Nanoparti-
cle dispersion and aggregation in polymer nanocomposites: insights from molecular
dynamics simulation. Langmuir, 27(12):7926{7933, 2011.
[56] JRMacdonaldandJCPhillips. Topologicalderivationofshapeexponentsforstretched
exponential relaxation. The Journal of chemical physics, 122(7):074510, 2005.
[57] M arioRGMarques,JakobWolff,ConradSteigemann,andMiguelALMarques. Neural
networkforceeldsforsimplemetalsandsemiconductors: constructionandapplication
to the calculation of phonons and melting temperatures. Physical Chemistry Chemical
Physics, 21(12):6506{6516, 2019.
[58] John Mongan, David A Case, and J Andrew McCammon. Constant ph molecular
dynamics in generalized born implicit solvent. Journal of computational chemistry,
25(16):2038{2048, 2004.
[59] P Mulvaney, V Swayambunathan, F Grieser, and D Meisel. Effect of. zeta. potential
on electron transfer to colloidal iron oxides. Langmuir, 6(3):555{559, 1990.
[60] A Nakano, RK Kalia, K-i Nomura, A Sharma, P Vashishta, and F Shimojo. Act v.
duin, iwa goddard, r. biswas, and d. srivastava. Comput. Mater. Sci, 38:642, 2007.
[61] Aiichiro Nakano. A space{time-ensemble parallel nudged elastic band algorithm for
molecular kinetics simulation. Computer Physics Communications, 178(4):280{289,
2008.
[62] KNomura,RKKalia,ANakano,PVashishta,andACTvanDuin. Mechanochemistry
of shock-induced nanobubble collapse near silica in water. Applied Physics Letters,
101(7):073108, 2012.
[63] Roberto Olender and Ron Elber. Calculation of classical trajectories with a very
large time step: Formalism and numerical examples. The Journal of chemical physics,
105(20):9299{9315, 1996.
[64] Lars Onsager. Reciprocal relations in irreversible processes. i. Physical review,
37(4):405, 1931.
[65] Lars Onsager. Reciprocal relations in irreversible processes. ii. Physical review,
38(12):2265, 1931.
[66] Lars Onsager and Stefan Machlup. Fluctuations and irreversible processes. Physical
Review, 91(6):1505, 1953.
77
[67] Roger Parsons. The electrical double layer: recent experimental and theoretical devel-
opments. Chemical Reviews, 90(5):813{826, 1990.
[68] Daniele Passerone and Michele Parrinello. Action-derived molecular dynamics in the
study of rare events. Physical review letters, 87(10):108302, 2001.
[69] JC Phillips. Stretched exponential relaxation in molecular and electronic glasses. Re-
ports on Progress in Physics, 59(9):1133, 1996.
[70] JCPhillips. Axiomatictheoriesofidealstretchedexponentialrelaxation(ser). Journal
of non-crystalline solids, 352(42-49):4490{4494, 2006.
[71] Steve Plimpton. Fast parallel algorithms for short-range molecular dynamics. Journal
of computational physics, 117(1):1{19, 1995.
[72] Marcel Potuzak, Roger C Welch, and John C Mauro. Topological origin of stretched
exponentialrelaxationinglass. The Journal of chemical physics,135(21):214502,2011.
[73] AneesurRahman. Correlationsinthemotionofatomsinliquidargon. Physical review,
136(2A):A405, 1964.
[74] Pankaj Rajak, Aravind Krishnamoorthy, Aiichiro Nakano, Priya Vashishta, and Rajiv
Kalia. Structural phase transitions in a mowse 2 monolayer: Molecular dynamics
simulations and variational autoencoder analysis. Physical Review B, 100(1):014108,
2019.
[75] Arunkumar Chitteth Rajan, Avanish Mishra, Swanti Satsangi, Rishabh Vaish, Hi-
roshi Mizuseki, Kwang-Ryeol Lee, and Abhishek K Singh. Machine-learning-assisted
accurate band gap predictions of functionalized mxene. Chemistry of Materials,
30(12):4031{4038, 2018.
[76] Rampi Ramprasad, Rohit Batra, Ghanshyam Pilania, Arun Mannodi-Kanakkithodi,
and Chiho Kim. Machine learning in materials informatics: recent applications and
prospects. npj Computational Materials, 3(1):1{13, 2017.
[77] John J Rehr and Robert C Albers. Theoretical approaches to x-ray absorption ne
structure. Reviews of modern physics, 72(3):621, 2000.
[78] Dale E. Sayers, Edward A. Stern, and Farrel W. Lytle. New technique for investigat-
ing noncrystalline structures: Fourier analysis of the extended x-ray|absorption ne
structure. Phys. Rev. Lett., 27:1204{1207, Nov 1971.
[79] Dawn L Scovell, Tim D Pinkerton, Valentin K Medvedev, and Eric M Stuve. Phase
transitions in vapor-deposited water under the in
uence of high surface electric elds.
Surface science, 457(3):365{376, 2000.
[80] Yohan Seepersad, Alexander Fridman, and Danil Dobrynin. Anode initiated impulse
breakdown in water: the dependence on pulse rise time for nanosecond and sub-
nanosecond pulses and initiation mechanism based on electrostriction. Journal of
Physics D: Applied Physics, 48(42):424012, 2015.
78
[81] Christopher J Shallue, Jaehoon Lee, Joseph Antognini, Jascha Sohl-Dickstein, Roy
Frostig, and George E Dahl. Measuring the effects of data parallelism on neural net-
work training. arXiv preprint arXiv:1811.03600, 2018.
[82] DEShaw,PMaragakis,KLindorff-Larsen,SPiana,RODror,MPEastwood,JABank,
JM Jumper, JK Salmon, and Y Shan. Molecular dynamics simulations of lipids. M2
Proton Channel from In
uenza A: Example of Structural Sensitivity to Environment,
301:1796, 2013.
[83] Fuyuki Shimojo, Rajiv K Kalia, Aiichiro Nakano, and Priya Vashishta. Embedded
divide-and-conquer algorithm on hierarchical real-space grids: parallel molecular dy-
namicssimulationbasedonlinear-scalingdensityfunctionaltheory. Computer Physics
Communications, 167(3):151{164, 2005.
[84] Mikhail N Shneider and Mikhail Pekker. Pre-breakdown processes in a dielectric
uid
in inhomogeneous pulsed electric elds. Journal of Applied Physics, 117(22):224902,
2015.
[85] MN Shneider and M Pekker. Cavitation in dielectric
uid in inhomogeneous pulsed
electric eld. Journal of Applied Physics, 114(21):214906, 2013.
[86] MN Shneider and M Pekker. Rayleigh scattering on the cavitation region emerging in
liquids. Optics letters, 41(6):1090{1093, 2016.
[87] Raymond G Sierra, Hartawan Laksmono, Jan Kern, Rosalie Tran, Johan Hattne,
RobertoAlonso-Mori,BenediktLassalle-Kaiser,CarinaGl ockner,JuliaHellmich,Don-
ald W Schafer, et al. Nano
ow electrospinning serial femtosecond crystallography.
Acta Crystallographica Section D: Biological Crystallography, 68(11):1584{1587, 2012.
[88] Timothy W Sirk, Stan Moore, and Eugene F Brown. Characteristics of thermal con-
ductivity in classical water models. The Journal of chemical physics, 138(6):064505,
2013.
[89] Vojislav R Stamenkovic, Dusan Strmcnik, Pietro P Lopes, and Nenad M Markovic.
Energy and fuels from electrochemical interfaces. Nature materials, 16(1):57{69, 2017.
[90] BrianCHSteeleandAngelikaHeinzel. Materialsforfuel-celltechnologies. InMaterials
For Sustainable Energy: A Collection of Peer-Reviewed Research and Review Articles
from Nature Publishing Group, pages 224{231. World Scientic, 2011.
[91] WilliamCSwope, HansCAndersen, PeterHBerens, andKentRWilson. Acomputer
simulation method for the calculation of equilibrium constants for the formation of
physical clusters of molecules: Application to small water clusters. The Journal of
Chemical Physics, 76(1):637{649, 1982.
[92] KeisukeTakahashiandYuzuruTanaka. Materialsynthesisanddesignfromrstprinci-
ple calculations and machine learning. Computational materials science, 112:364{367,
2016.
79
[93] Aidan P Thompson, Steven J Plimpton, and William Mattson. General formulation
of pressure and stress tensor for arbitrary many-body interaction potentials under
periodic boundary conditions. The Journal of chemical physics, 131(15):154107, 2009.
[94] Vladimir P Torchilin. Recent advances with liposomes as pharmaceutical carriers.
Nature reviews Drug discovery, 4(2):145, 2005.
[95] REG Van Hal, Jan CT Eijkel, and Piet Bergveld. A general model to describe the
electrostatic potential at electrolyte oxide interfaces. Advances in colloid and interface
science, 69(1-3):31{62, 1996.
[96] Juan-Jesus Velasco-Velez, Cheng Hao Wu, Tod A Pascal, Liwen F Wan, Jinghua Guo,
David Prendergast, and Miquel Salmeron. The structure of interfacial water on gold
electrodes studied by x-ray absorption spectroscopy. Science, page 1259437, 2014.
[97] Loup Verlet. Computer" experiments" on classical
uids. i. thermodynamical proper-
ties of lennard-jones molecules. Physical review, 159(1):98, 1967.
[98] Michael G Walter, Emily L Warren, James R McKone, Shannon W Boettcher, Qixi
Mi, Elizabeth A Santori, and Nathan S Lewis. Solar water splitting cells. Chemical
reviews, 110(11):6446{6473, 2010.
[99] Hao Wang and Weitao Yang. Force eld for water based on neural network. The
journal of physical chemistry letters, 9(12):3232{3240, 2018.
[100] Juan Wang, Xiaoyu Yang, Zhi Zeng, Xiaoli Zhang, Xushan Zhao, and Zongguo Wang.
New methods for prediction of elastic constants based on density functional theory
combinedwithmachinelearning. Computational Materials Science,138:135{148,2017.
[101] Tian Wang, Cheng Zhang, Hichem Snoussi, and Gang Zhang. Machine learning
approaches for thermoelectric materials research. Advanced Functional Materials,
30(5):1906041, 2020.
[102] K. R. Weninger, C. G. Camara, and S. J. Putterman. Observation of bubble dynamics
within luminescent cavitation clouds: Sonoluminescence at the nano-scale. Phys. Rev.
E, 63:016310, Dec 2000.
[103] Thomas Winter. X-ray absorption spectroscopy of colloidal tin dioxide nanoparticles.
Master's thesis, University of Duisburg-Essen, 2019.
[104] Markus Winterer.
[105] Markus Winterer. Reverse monte carlo analysis of extended x-ray absorption ne
structure spectra of monoclinic and amorphous zirconia. Journal of Applied Physics,
88(10):5635{5644, 2000.
[106] Kunio Yasue. The role of the onsager{machlup lagrangian in the theory of stationary
diffusion process. Journal of Mathematical Physics, 20(9):1861{1864, 1979.
80
[107] Yifan Ye, Cheng Hao Wu, Liang Zhang, Yi-Sheng Liu, Per-Anders Glans-Suzuki,
and Jinghua Guo. Using soft x-ray absorption spectroscopy to characterize elec-
trode/electrolyte interfaces in-situ and operando. Journal of Electron Spectroscopy
and Related Phenomena, 221:2{9, 2017.
[108] Q Zheng, DJ Durben, GH Wolf, and CA Angell. Liquids at large negative pressures:
water at the homogeneous nucleation limit. Science, 254(5033):829{832, 1991.
81
Appendix A
Cavitation Supporting Information
Force Fields
SnO2 Force Field
B
ij
[kJ/mol]
ij
[

A] C
ij
[kJ/mol A
6
]
Sn-Sn 90547100 0.205 608
Sn-O 10364150 0.169 1544
O-O 465467 0.273 4699
O
water
-O
water
279737 3.51485 3672.79
Table A.1: Buckingham: E
ij
= B
i
jexp(r
ij
=
ij
)C
ij
=r
6
ij
taken from [7]
A
ij
[kJ/mol

A
12
] C
ij
[kJ/mol

A
6
]
Sn-O
water
600451 1167
O-O
water
1978673 3502
(H)O-O
water
2561695 3502
O
water
-O
water
2634130 2617
Table A.2: Lennard-Jones: E
ij
= A
ij
=r
12
ij
C
ij
=r
6
ij
taken from [7]
ϵ
HB
[kJ/mol]
HB
[

A]
H-H 1.346 2.75
Table A.3: Anisotropic Hydrogen-bonding term: E
HB
ij
=
ϵ
HB
[
5
(

HB
r
ij
)
12
6
(

HB
r
ij
)
10
]
cos
4
The angle in the hydrogen-bonding interaction
is the angle between the oxygen atom accepting the hydrogen bond and the oxygen atom
donating the hydrogen bond. Taken from [40].
The force eld developed by Bandura et. al was used to describe the SnO
2
-SnO
2
interac-
tionsaswellastheSnO
2
-H
2
OinteractionsinLAMMPS[7](TablesA.1,A.2,A.3). TheSnO
2
parameters were designed to work with a modied SPC/E water force eld with harmonic
angle bending instead of rigid molecules. The intra and inter forces between SnO
2
and H
2
O
are described by a combination of electrostatic and short-range forces. The Buckingham
potential describes the short range interactions for the SnO
2
while the Lennard-Jones model
82
Atom q (e) [

A]
Sn 2.050 -
O -1.025 -
O
water
2.0 0.7114
H 0.597 0.4592
M -1.194 0.5483
D -2.0 0.7114
Table A.4: Atom charges taken from [7] and [40] M and D refer to the dummy atom and
Drud charge respectively.
is used for the SnO
2
-H
2
O interactions. The coefficients of both forces are given in Tables
A.1, A.2, A.3.
Theatomicchargevalues(TableA.4)wereextractedfromDFTcalculationsusingtheVi-
enna Ab initio Simulation Package (VASP). The coefficients for the SnO2 short range Buck-
ingham interactions were calculated using the Generalized Utility Lattice Program (GULP)
with the exception of the disperssion coefficeint, C
ij
for the Sn-O interaction. The Sn-O dis-
persion coefficient was theoretically determined by using data from ref [19] and the method
of Fowler et al [25].
H2O Force Field
The Hydrogen-Bonding Polarizable model of water was used to describe the H
2
O-H
2
O in-
teractions. Specically, the code provided in the supplemental material of [40] was used to
implement the force eld in LAMMPS. The HBP model uses the rigid molecule structure
of TIP4P with an O-H-O bond angle of 104.52
◦
and an O-H bond length of 0.9572

A. The
dummy atom with the charge equal and opposite of the two hydrogen atoms was placed
at the bisector of the HOH angle. To capture the polarizability of water, a single Drude
oscillator is attached to the oxygen atom.
Inlieuofpointcharges, themodelusesGaussianchargescenteredonthepositionofeach
atom to represent the electrostatic interactions. The charge density on atom i is given by:

i
=
q
i
(2
2
i
)
3=2
exp
(
jrr
i
j
2
2
2
i
)
The values of q
i
and
i
for the 4 different species in each water molecule are shown in Table
A.4. The anisotropic hydrogen-bonding term is used to account for the short-ranged effects
causedbyhydrogenbondingwhiletheVanderWaalsinteractionsbetweenoxygenatomsare
modeled through the Buckingham potential. The charge values for the oxygen and Drude
oscillator were found to have little impact on the model when tested in [40]. The only
electrostatic interataction between the water molecules and the nanoparticle comes from the
hydrogen and dummy atoms.
A polarizable model of water was selected in part due to the ability to reproduce the
dielectric constant of liquid water. The HBP model has a dielectric constant of 764
compared to an experimetnal value of 78 at 298.15 K at 1 bar [40] which is a signicant
improvement over TIP4P water which predicts a dielectric value of 58.
83
The HBP O-O parital correlatoin function over estimates the rst peak and under esti-
mates the rst minimum implying that the model predicts a more structured form of water
than found in experiment [40].
The ab initio work done by Bandura et. al shows that water has a tendency to disso-
ciatively adsorb on the surface of SnO2; however, our simulations with HBP did not allow
bond breaking or formation. To verify the validity of the HBP water with the Bandura
SnO2 eld, water adsorption simulations were performed. Water droplets of radii 7 and 10
Angstrom were placed on the (001) surface of a 46.004x47.86x47.86 Angstrom block of SnO2
at room temperature (297.3 K). In each simulation the water droplet was allowed to expand
over the surface over 300 ps; however, the water stabilized on the surface within 150 ps in
both cases. The Ab Initio results of Bandura et al. say that the associative adsorption is
unlikely on the 001 surface; however mixed adsorption has an energy of -117.8 eV. The HBP
force eld produced energies of -163 kJ/mol and -173 kJ/mol for rst and second layer water
molecules respectively. The HBP force eld produces adsorption energies of correct order of
magnitudewithrespecttotheabinitioresults. WewouldnotexpectfulldelityastheHBP
model uses different atomic charges than the SPC/E water model and it does not allow for
dissociation as mentioned above.
Water without Nanoparticle
(a) (b)
FigureA.1: Tindioxidenanoparticleinwate(a)randpurewater(b)withnoappliedelectric
eld. Transient cavities on the order of 100

A
3
and smaller form.
Cavitation was observed in every electric eld simulation with applied electric elds,
E
app
, ranging from 0.042 V/

A to 0.25 V/

A. Water only simulations were also performed for
high electric eld values E
app
= 0:25;0:125; however, no cavitation was observed within 200
ps. Before the electric eld has been applied, when the system is still in equilibrium, small
cavities arise and collapse do to thermal
uctuations. These transient cavities are on the
order of 100

A
3
and occur both in the presence and absence of a nanoparticle (Figure A.1).
Once the electric eld is applied the water atoms rapidly polarize with a relaxation time
dependent on the strength of the electric eld. At some time, t
onset
, cavities begin to congeal
84
withinonenanoparticleradius(r
np
= 1:5nm)ofthenanoparticle. Thecavitybeginstogrow
either through the absorption of the thermally induced cavities or the pressure differential
caused by the electric eld in the void and in the water.
Simulation Details
Creation of the tin dioxide nanoparticle began by rst simulating bulk SnO2 crystal in a
84x95x95

A
3
box with 62400 atoms. The bulk structure was heated to 1000 K over 25 ps
with a Berendsen thermostat in the NVE ensemble. The Berendsen thermostat wasdisabled
and the system continued to run for another 150 ps. From this simulation, a spherical
nanoparticle with a 15 angstrom radius (1110 atoms) was extracted. Charge neutrality and
stoichometrywereensuredbyremovingtheappropriatenumberofSnandOatomsbasedon
their radial distance from the center of the nanoparticle. The nanoparticle was then heated
to 1400 K over 25s with a Langevin thermostat to adequately sample phase space before
being cooled and then cooled to 300 K.
Separately,a139.8x93.5x93.5

A
3
waterboxofHBPwaterwaspreparedto300K,1atm
and a density of 0.997 g/cm
3
. A 16.5

A radius sphere of water was removed from the center
of the box and replaced with the SnO
2
nanoparticle. The combined SnO
2
-H
2
O system was
run for 100 ps while the atoms (H, O, Sn, M) were attached to a Langevin thermostat at 300
K and the Drude oscillators were simultaneously connected to a 200 K Langevin thermostat.
The total system was comprised of 199970 atoms.
The combined water and nanoparticle system was used as the starting point for the
simulations below. Simulations ran from 100 ps to 1.2 ns. Uniform electric elds ranging
from 0.042 to 0.25 V/

A were applied along the x direction on the simulation timescale using
the"xeeld"commandinLAMMPS.Thisaddsaneffectiveelectriceldbydirectlyadding
a force f = qE to each atom.
To prevent Joule heating, a Langevin thermostat was applied and kept at 300 K for the
atoms and 200 K for the Drude oscillator particles. Unfortunately, modern experiments able
toapplysufficientelectriceldsareunabletosimulataneouslydealwiththeresultantheating
making the application of the thermostats somewhat unrealistic [22]. Solving the equations
ofmotionfortheDrudeparticlesintheself-consistenteld(SCF)regimerequiresminimizing
the energy of the Drude particles at each timestep which is computationally expensive [50].
In keeping with standard practice, the Drude particle was assigned a small mass to allow for
faster relaxation and thermalized to a lower temperature than the bulk system. Lamoureux
et al showed that in keeping the Drude temperature low, the thermodynamic and dynamic
properties of the system closely approximated those given by SCF [50]. The lower Drude
particle temperature serves to minimize the kinetic energy of the Drude particles while still
allowing them the ability to respond to changes in the electric eld or movement of their
corresponding oxygen atoms.
The system's pressure was calculated in LAMMPS. The force due to the external electric
eld was taken account in the LAMMPS calculated virial. Pressure distributions were ob-
served using the compute/stress/atom feature and binning the atoms into rectangular slices
along the x direction. In the longer simulations (those on the order of a nanosecond), the
total system pressure show large
uctuations on the order of MPa. At the simulation onset
however, thetotalsystempressurevalueswerewithinanorderofmagnitudeofthepredicted
85
negative pressure by Equation (2) in the main text. For example, using = 1:4, ϵ = 79
the predicted negative pressure for an applied eld was 0.125 V/

A is -700 MPa while the
observed system pressure oscillated around -150 MPa for the rst 50 ps.
Void volume measurement
0 500 1000
2
2.5
3
3.5
4
10
4
0.250 V/A
0.200 V/A
0.160 V/A
0.125 V/A
0.109 V/A
0.096 V/A
0.073 V/A
0.063 V/A
0.052 V/A
0.042 V/A
(a)
0 500 1000
0
0.5
1
1.5
2
10
4
0.250 V/A
0.200 V/A
0.160 V/A
0.125 V/A
0.109 V/A
0.096 V/A
0.073 V/A
0.063 V/A
0.052 V/A
0.042 V/A
(b)
(c)
Figure A.2: Cavity volume analysis. (a) Raw void volume versus time. (b) Cavity volume
versus time. Calculated using voxel chains of 15 or larger. (c) Spherical and cylindrical
growth rates as a function of applied electric eld. Some of the simulations runs were not
run for sufficient lengths of time to see cylindrical growth. Comparing, for example, the
0.042 curve from part (b), the red star corresponds to the rst linear region from 600 ps to
750ps while the green star corresponds to the slope in the region between 760 ps and 850
ps.
The cavity volume was calculated using a home made python script. The simulation box
is divided into a grid of cubic voxels. Atoms are assigned to a voxel based on their position
coordinates. The void volume is calculated by summing the number of empty voxels. In
86
order to more accurately measure the rate of growth of the main cavity, only voxel chains
above a certain threshold, n
thresh
value were considered. The chains of contiguous empty
voxels were identied through depth rst search. To analyze the data from the simulations,
avoxellengthof3

Awasused. TherawvoxelvoiddataisshowninFigureA.2a(n
thresh
= 1)
while the ltered data is shown (n
thresh
= 15) in Figure A.2 b.
From the ltered data, the two regions of growth, spherical and cylindrical, can be more
easilyidentied; however, duetothe shortsimulationlengthsforsomeoftheappliedelectric
eld runs, there are missing data points. Generally, the rate of growth in both the phases
tended to increase as electric eld was increased (A.2 c). The volume growth may scale
linearly with electric eld strength; however, currently the data is too noisy to determine
this denitively.
Discussion
When subtracting out the initial time, a value of = 1=3 arises from our data see (A.3).
Analysis of many sets of dielectric data produced a value for equal to 1/3 for low fre-
quency electric elds [70]. In a strong electric eld, the excitations (alignments of the polar
molecules) do not nd sinks but instead relax locally to the applied external eld. In a
parallel plate capacitor, there is only one effective dimension, that between the two plates
which yields a value for f of 1/3. Plugging into Eq 7 from the main paper in turns gives
= 1=3 as reported which is in agreement with our alternatively calculated value.
Identifying the long timescale value of will provide evidence for the mechanism behind
cavitation. If, the value is = 3=7 then the most likely mechanism, is that the exitations,
i.e the thermally activated cavities, diffuse into the main cavity by the nanoparticle. The
next step would be then identifying the mechanism in which the thermal cavities interact.
if instead, the long term value ends up being 1/3, then the mechanism will likely be tied
directly to long term relaxation of the water dipoles in the electric eld. Given the stability
of the electric eld after the initial pulse, the author nd this unlikely; however this needs to
be looked into further. In order to verify the value for , longer simulations will need to be
run. Typically, is calculated bylooking atdata correspondingto 100 [69]. The calculated
from our simulations is on the order of 100 ps. This would require simulations lasting 10
ns; however, cavitation may not occur over such lengthy timescales [80].
Video
The attached video is the tin dioxide simulation (red and yellow) with an applied electric
eld of 0.629 V/

A from 100 to 300 ps. Initially, one can see the formation of local density

uctuations. These eventually congeal into the cavity. The void (teal) appears over the
course of 100s of picoseconds.
87
(a)
(b) (c)
FigureA.3: (a)

E vsT:Bluemarkersarethedatausedtocalculatetheexponentialrelaxation
coefficient . The temporal data spans an order of magnitude (60 to 635 ps). The green
(blue) line is plotted stretched exponential with = 0:324 ( = 0:424) and = 90:16 ps
( = 134:6 ps). The blue line is plotted stretched exponenital (b) log(log(

E) vs logT
onset
.
(c)log(log(

E) vs log(T
onset
min(T
onset
))
88
Appendix B
NNMD S.I
Lennard-Jones Sytstem
Our neural network (NN) for the principle of least action is tested against the ground truth
moleculardynamics(MD)simulationsofasysteminwhichatomsinteractviaLennard-Jones
(LJ) potential. In reduced units, the pair-wise LJ interaction is given by:
V
ij
= 4
(
1
jq
i
q
j
j
12

1
jq
i
q
j
j
6
)
(B.1)
The interaction parameters in the LJ potential are given in table S1. Molecular dynamics
simulation was initiated from the face-centered cubic (FCC) structure. The reduced number
density of the system is 0.787 and the reduced time unit = 2:068 ps. Periodic boundary
conditions were imposed, and equations of motion were integrated with the velocity-Verlet
algorithm using a timestep of 2 fs. The system was equilibrated for 2 ns at a reduced
temperature of 0.695. We use the coordinates and velocities of atoms in the equilibrated
state as the initial conguration for the NN simulation. To obtain a nal state for our
boundary-value problem, we run another MD simulation for 0.5 using the same starting
state.
Pre-training
The pre-training strategy we use is as follows: we rst randomize the parametersfw; bg,
then with the knowledge of the initial and nal state of the system, we pre-train the network
using the following loss function:
L
pre
(fw; bg) = (B.2)
where ⃗ q
0
and ⃗ q
T
collectively represent the initial and nal positions of atoms, and⃗ v
0
and
⃗ v
T
are initial and nal velocities of atoms in the system, respectively. Note that the rst
two terms are constraints on initial and nal congurations of the system, while the last
Simulation Parameter Value
ϵ=k
b
(K) 125.7
(nm) 0.3345
m (amu) 39.948
Lattice Constant (nm) 0.575
Table B.1: Lennard Jones interaction parameters and lattice constant.
89
two terms provide an estimate of trajectories by making particles move at different but
constant velocities for the rst and second half of the trajectories. Although this pretraining
method has led to interesting results as presented in the main text, there is still a chance
that particles may get too close or even cross their trajectories after pre-training and give
rise to a tremendous loss at the beginning of the main training iteration. One way to avoid
this during the pre-training is to put large weights on the minimum inter-particle distance
if it becomes smaller than a threshold value.
Results for 108 and 256-atom systems
Inadditiontothe500-atomsystem, westudiedsystemswith108and256atoms. Inthenext
two gures, we compare NN results against MD simulations for a LJ system of 256 atoms.
Figures B.1 shows trajectories generated by the NN for four randomly chosen atoms. The
NN trajectories are nearly coincident with the ground-truth MD trajectories. The averaged
mean-square error between NN and MD trajectories is on the order of 10-5 (LJ units).
Figures B.5 and B.6 shows components of the loss-function and mean-square errors be-
tween MD and NN trajectories for the 256- and 108-atom systems, respectively. Figure
B.5(a) and B.6(a) show how different parts of the loss function decreases with the number
of epochs. The energy-conservation constraint dominates at the beginning of the training,
but it becomes very close to the action term toward the end of the training. The boundary
loss reduces quickly at an early stage of the training. Figure B.5(b) and B.6 (b) show the
mean-square error between NN and the ground truth MD trajectories; and Fig. B.5 (c) and
B.6(c) show the mean-square error in the time domain over which the NN is trained.
90
Figure B.1: NN and MD trajectories for four randomly selected atoms in a system of 256
LJ atoms. Atomic trajectories computed by the NN (blue) are nearly coincident with MD
trajectories (red).
Figure B.2: a) NN and MD results for potential, kinetic and total energies per atom as a
function of time for a system of 256 atoms. (b) Radial distribution function and (c) velocity
auto-correlation function show the near-perfect agreement between the NN and MD for a
256-atom system.
91
Figure B.3: NN and MD trajectories for four randomly selected atoms in a system of 108
LJ atoms. Atomic trajectories computed by the NN (blue) are nearly coincident with MD
trajectories (red).
Figure B.4: (a) Shows that the NN results for potential, kinetic and total energy per particle
as a function of time are nearly coincident with the ground truth MD results. Figures S2
(b) and (c) compare NN and MD results for the radial distribution function and velocity
auto-correlation function for the 108-atom system.
92
Figure B.5: NN results for a 256-atom LJ system. (a) loss-function components vs the
number of epochs. (b) Mean-squared error between MD and NN for trajectories (red) and
velocities (blue) as a function of epochs. (c) Mean-squared error between MD and NN
trajectories (red) and momenta (blue) vs time.
Figure B.6: Results for the 108-atom system. (a) loss function versus number of training
epochs. (b)mean-squarederrorbetweenMDandNNresultsforpositions(red)andvelocities
(blue). (c) The red curve is the mean-squared error between MD and NN trajectories (Q)
as a function of time. The blue curve is the error between MD and NN velocities (P).
93

Abstract (if available)

Abstract Cavitation phenomenon in dielectric fluids has been a recent topic of interest in theory and experiment. We study a dielectric fluid-nanoparticle system subjected to an external electric field using molecular-dynamics simulations. Electric fields ranging from 0.042 to 0.25 V/ ̊A are applied to a water and tin dioxide system. Cavitation is observed in simulations with both SPC/E water and the hydrogen bonding polarizable model. The cavitation onset time displays a stretched exponential relaxation response with respect to the applied electric field with an exponent β = 0.423 ± 0.08. This is in accordance with the exact theoretical value for systems with long-ranged forces. Cavity growth rates are divided into two phases, a spherical growth phase and cylindrical. Both are reported on as a function of the applied electric field. The structure of the electric field is analyzed both spatially and temporally. Due to their relatively high surface area to volume ratio, nanoparticles are potential candidates to investigate the solid-liquid interface. X-ray Absorption Spectroscopy experiments were performed on powdered and colloidal SnO₂ nanoparticles using the B-12 beamline at Argonne National Laboratory. The EXAFS data was analyzed using reverse Monte Carlo (RMC) simulations with starting configurations generated from Rietveld refinement results of X-ray diffraction data. ❧ The principle of least action is the cornerstone of classical mechanics, special and general theory of relativity, quantum mechanics, and thermodynamics. A description of the training of neural networks to learn the principle of least action through minimizing the Onsager-Machlup action is provided. The network, with slight alterations to the loss function, is used to solve both initial value problems (IVP) and boundary value problems (BVP). The resultant phase-space trajectories for the BVP network are in excellent agreement with the trajectories obtained by conventional molecular dynamics. The potential and kinetic energies as a function of time as well as radial distribution and velocity autocorrelation functions computed by both neural networks are in excellent agreement with the corresponding results from the “ground-truth” MD simulation. From the first epoch, the neural network approach produces as full trajectory which increases in accuracy as training progresses. The time step is a factor of twenty larger than the time step in MD simulation for the BVP and 5 times larger for the IVP. The neural-network approach can help scientists discover new phenomena hidden in data, impute missing data in experiments and simulations and, most importantly, allow atomistic simulations to reach long timescales.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Simulation and machine learning at exascale

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Heat-initiated oxidation of aluminum nanoparticles

PDF

Designing data-effective machine learning pipeline in application to physics and material science

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Hypervelocity impact damage in alumina

PDF

Compiler and runtime support for hybrid arithmetic and logic processing of neural networks

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites

PDF

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

PDF

Exploring properties of silicon-carbide nanotubes and their composites with polymers

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

Asset Metadata

Creator Jackson, Shane William (author)

Core Title Physics informed neural networks and electrostrictive cavitation in water

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 06/22/2020

Defense Date 05/15/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cavitation in water,dielectric fluids,feed-forward,molecular dynamics,neural networks,OAI-PMH Harvest,Physics,tin dioxide

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kalia, Rajiv (committee chair), Boedicker, James (committee member), Nakano, Aiichiro (committee member), Vashishta, Priya (committee member), Winterer, Markus (committee member)

Creator Email shebaiscool@gmail.com,swjackso@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-321146

Unique identifier UC11665695

Identifier etd-JacksonSha-8614.pdf (filename),usctheses-c89-321146 (legacy record id)

Legacy Identifier etd-JacksonSha-8614.pdf

Dmrecord 321146

Document Type Dissertation

Rights Jackson, Shane William

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA