Citation
Development of Classical Interatomic Potentials for Applications in Corrosion and Phase Transitions

Material Information

Title:
Development of Classical Interatomic Potentials for Applications in Corrosion and Phase Transitions
Creator:
Noordhoek, Mark J
Place of Publication:
[Gainesville, Fla.]
Florida
Publisher:
University of Florida
Publication Date:
Language:
english
Physical Description:
1 online resource (5 p.)

Thesis/Dissertation Information

Degree:
Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Materials Science and Engineering
Committee Chair:
PHILLPOT,SIMON R
Committee Co-Chair:
DOUGLAS,ELLIOT PAUL
Committee Members:
SINNOTT,SUSAN B
MYERS,MICHELE V
CHEN,YOUPING
Graduation Date:
8/9/2014

Subjects

Subjects / Keywords:
Atomic interactions ( jstor )
Atoms ( jstor )
Corrosion ( jstor )
Databases ( jstor )
Eggshells ( jstor )
Energy ( jstor )
Hydrogen ( jstor )
Oxygen ( jstor )
Simulations ( jstor )
Zirconium ( jstor )
Materials Science and Engineering -- Dissertations, Academic -- UF
comb -- corrosion -- potentials -- zirconium
Genre:
bibliography ( marcgt )
theses ( marcgt )
government publication (state, provincial, terriorial, dependent) ( marcgt )
born-digital ( sobekcm )
Electronic Thesis or Dissertation
Materials Science and Engineering thesis, Ph.D.

Notes

Abstract:
Classical molecular dynamics (MD) is an established method used to simulate atomic scale phenomena in materials systems. The key component of MD is the description of the interaction energy between atoms. In this work, two independent formalisms of classical interatomic potentials are used to describe various aspects of zirconium-based and barium titanate-based systems. A third-generation Charge-Optimized Many-Body (COMB) potential for the zirconium-zirconium oxide-zirconium hydride (Zr-O-H) system is presented. These materials are widely found in the fuel cladding components of nuclear reactors. During energy production, Zr undergoes a number of mechanical and electrochemical transformations, which may lead to catastrophic failure if allowed to continue indefinitely. Thus, a better understanding of these processes is of critical importance. COMB allows various mechanisms of these transformations to be studied within the same simulation. The metallic Zr interactions in this work correctly describe many structural properties including stacking fault energies and surface energies. Mechanical testing of Zr shows that the slip mechanisms in MD match those in experiments. Deposition of O2 and H2O molecules onto various low-index Zr surfaces shows the resulting corrosion behavior is dependent on the surface structure and deposited species. Joint MD and experimental studies on the structural properties of the Ba1-xSrxTiO3 and Ba1-xCaxTiO3 systems are presented. These systems are ferroelectric for small Sr/Ca concentrations (< 30%), which is exploited in their use for many commercial devices such as capacitors. The mechanism of the phase transition between the ferroelectric tetragonal phase and paraelectric cubic phase is still an open debate. In this work, the shell model is used to describe interatomic interactions for each system. Direct comparison between the structural properties predicted in MD and those observed in experiment is provided. The model significantly overestimates the distortion within the tetragonal phase, which prevents a reliable determination of the transformation mechanism. Additional work to develop a new shell model description for BaTiO3 is presented. However, efforts to produce an improved potential have not been successful. Suggestions for improved potentials are provided. ( en )
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Thesis:
Thesis (Ph.D.)--University of Florida, 2014.
Local:
Adviser: PHILLPOT,SIMON R.
Local:
Co-adviser: DOUGLAS,ELLIOT PAUL.
Statement of Responsibility:
by Mark J Noordhoek.

Record Information

Source Institution:
UFRGP
Rights Management:
Applicable rights reserved.
Resource Identifier:
969976736 ( OCLC )
Classification:
LD1780 2014 ( lcc )

Downloads

This item has the following downloads:


Full Text

PAGE 1

GravitationalWavesfromtheSoundofaFirstOrderPhaseTransitionMarkHindmarsh,1,2 ,*StephanJ.Huber,1 , †KariRummukainen,2 , ‡andDavidJ.Weir2 ,§1DepartmentofPhysicsandAstronomy,UniversityofSussex,Falmer,BrightonBN19QH,UnitedKingdom2DepartmentofPhysicsandHelsinkiInstituteofPhysics,UniversityofHelsinki,PL64,FI-00014,Helsinki,Finland (Received16April2013;revisedmanuscriptreceived12October2013;published27January2014) Wereportonthefirstthree-dimensionalnumericalsimulationsoffirst-orderphasetransitionsintheearly Universetoincludethecosmicfluidaswellasthescalarfieldorderparameter.Wecalculatethe gravitationalwave(GW)spectrumresultingfromthenucleation,expansion,andcollisionofbubblesofthe low-temperaturephase,forphasetransitionstrengthsandbubblewallvelocitiescoveringmanycasesof interest.WefindthatthecompressionwavesinthefluidcontinuetobeasourceofGWslongafterthe bubbleshavemerged,aneweffectnottakenproperlyintoaccountinpreviousmodelingoftheGWsource. Forawiderangeofmodels,themainsourceoftheGWsproducedbyaphasetransitionis,therefore,the soundthebubblesmake.DOI: 10.1103/PhysRevLett.112.041301 PACSnumbers:98.80.Cq,04.30.Db,47.75.+f,95.30.LzInahotbigbang,therewerephasetransitionsinthe earlyUniverse [1,2] ,whichmaywellhavebeenoffirst order;onemajorconsequenceofsuchatransitionwouldbe thegenerationofgravitationalwaves [3 – 8] .Theelectroweaktransitioninthestandardmodelisknowntobea crossover [9 – 11] ,butitmaybefirstorderinminimal extensionsofthestandardmodel [12 – 17] .Itistherefore essentialtoproperlycharacterizetheexpectedpower spectrumfromfirst-orderphasetransitions. First-orderphasetransitionsproceedbythenucleation, growth,andmergerofbubblesofthelow-temperature phase [3,18 – 25] .Thecollisionofthebubblesisaviolent process,andboththescalarorderparameterandthefluidof lightparticlesgenerategravitationalwaves. Numericalstudieshavebeencarriedoutonthebehavior ofbubblesinsuchaphasetransitionusingspherically symmetric( 1 þ 1 )-dimensionalsimulations [23,24] .The calculationofthegravitationalwavespectrumhasbeen refinedintheinterveningyears,notablyusingthesemianalyticenvelopeapproximation [5,7,8,26,27] (butsee Ref. [28] foranalternativeapproach).Fullythree-dimensionalsimulationsofthescalarfieldonlyhavebeencarried out [29] ,qualitativelysupportingtheenvelopeapproximationandpointingoutimportantgravitationalwaveproductionfromthescalarfieldafterthebubblemerger. Inahotphasetransition,thefluidplaysanimportant role,firstasabrakeonthescalarfieldandsecondasa sourceofgravitationalwavesitself.Thefluidhasgenerally beenassumedtobeincompressibleandturbulent [30 – 33] . Animportantquestionforthegravitationalwavepower spectrumisthevalidityofthismodeling,whichgenerally borrowsfromtheKolmogorovtheoryofnonrelativistic drivenincompressibleturbulence. InthisLetter,wereportonthefirstfullythree-dimensional simulationofbubblenucleationinvolvingacoupled field-fluidsystem.Wemakeuseofthesesimulationsto calculatethepowerspectrumofgravitationalradiationfrom afirst-orderphasetransition,forarangeoftransition strengthsandbubblewallvelocitiesrelevantforanelectroweaktransitioninextensionsofthestandardmodel.Wefind thatthecompressionwavesinthefluid — soundwaves — continuetobeanimportantsourceofgravitationalwavesfor uptoaHubbletimeafterthebubblemergerhascompleted. ThisbooststhesignalbytheratiooftheHubbletimetothe transitiontime,whichcanbeordersofmagnitude. ThesystemdescribingthematterintheearlyUniverse consistsofarelativisticfluidcoupledtoascalarfield, whichacquiresaneffectivepotential V ð ;T Þ¼ 1 2 ð T2 T2 0Þ 2 1 3 T 3þ 1 4 4: (1) Therest-framepressure p andenergydensity are ¼ 3 aT4þ V ð ;T Þ T V T ;p ¼ aT4 V ð ;T Þ (2) with a ¼ð 2= 90 Þ g ,and g theeffectivenumberofrelativisticdegreesoffreedomcontributingtothepressureat temperature T .Thestress-energytensorforascalarfield andanidealrelativisticfluid Uis T¼ 1 2gð Þ2þ½ þ p UUþ gp (3) wherethemetricconventionis( þþþ ).Thescalarfield potentialisincludedinthedefinitionof p .Wesplit T¼ 0 (nonuniquely)intofieldandfluidpartswitha dissipativetermpermittingtransferofenergybetweenthe scalarfieldandthefluid ¼ U [22,23] .This simplifiedmodelcanbeimprovedbutisadequatefor parametrizingtheentropyproduction [24] . Giventheseexpressions,theequationsofmotioncanbe derived.Forthefield,wehave PRL 112, 041301(2014) PHYSICALREVIEWLETTERSweekending 31JANUARY20140031-9007 = 14 = 112(4) = 041301(5)041301-1©2014AmericanPhysicalSociety

PAGE 2

þ 2 V ¼ W ð _ þ Vii Þ (4) where W istherelativistic factorand Viisthefluid 3-velocity Ui¼ WVi.Forthefluidenergydensity E ¼ W , contracting ½ Tfluidwith Uyields _ E þ ið EViÞþ p ½ W:þ ið WViÞ V W ð _ þ Vii Þ ¼ W2ð _ þ Vii Þ2: (5) Theequationsofmotionforthefluidmomentumdensity Zi¼ W ð þ p Þ Uiread _ Ziþ jð ZiVjÞþ ip þ V i ¼ W ð _ þ Vjj Þ i : (6) Theprincipalobservableofinteresttousisthepower spectrumofgravitationalradiationresultingfrombubble collisions.Oneapproachistoproject Tijateverytimestep andthenmakeuseoftheGreen ’ sfunctiontocomputethe finalpowerspectrum [34,35] ;thisisquitecostlyin computertime.Instead,weusetheproceduredetailedin Ref. [36] .Weevolvetheequationofmotionforanauxiliary tensor uij, uij2uij¼ 16 G ð ijþ f ijÞ ; (7) where ij¼ i j and f ij¼ W2ð þ p Þ ViVj.Thephysicalmetricperturbationsarerecoveredinmomentumspace by hijð k Þ¼ ij;lmð ˆ k Þ ulmð t; k Þ ,where ij;lmð ˆ k Þ istheprojectorontotransverse,tracelesssymmetricrank2tensors. Wearemostinterestedinthemetricperturbationssourced bythefluid,asthefluidshearstressesgenerallydominate overthoseofthescalarfield,althoughitwillbeinstructive toalsoconsiderbothsourcestogether. Havingobtainedthemetricperturbations,wefindthat thepowerspectrumperlogarithmicfrequencyintervalis d GWð k Þ d ln k ¼ 1 32 GL3 k3ð 2 Þ3Z d j _ hlmð t; k Þj2: (8) Wesimulatethesystemonacubiclatticeof N3¼ 10243points,neglectingcosmicexpansionwhichisslowcomparedwiththetransitionrate.Thefluidisimplementedasa three-dimensionalrelativisticfluid [37] ,withdonorcell advection.Thescalarandtensorfieldsareevolvedusinga leapfrogalgorithmwithaminimalstencilforthespatial Laplacian.Principally,weusedlatticespacing x ¼ 1 T 1 candtimestep t ¼ 0 . 1 T 1 c,where Tcisthecriticaltemperatureforthephasetransition.Wehavecheckedthe latticespacingdependencebycarryingoutsinglebubble self-collisionsimulationsfor L3¼ 2563T 3 cat x ¼ 0 . 5 T 1 c, forwhichthevalueof GWat t ¼ 2000 T 1 cincreasedby 10%,whilethefinaltotalfluidkineticenergyincreasedby 7%.Simulatingwith t ¼ 0 . 2 T 1 cresultedinchangesof 0.3%and0.2%to GWandthekineticenergy,respectively. Startingfromasystemcompletelyinthesymmetric phase,wemodelthephasetransitionbynucleatingnew bubblesaccordingtotherateperunitvolume P ¼ P0exp ð ð t t0ÞÞ .Fromthisdistributionwegeneratea setofnucleationtimesandlocations(inasuitable untouchedregionofthebox)ateachofwhichweinsert astaticbubblewithaGaussianprofileforthescalarfield. Thebubbleexpandsandquicklyapproachesaninvariant scalingprofile [23] . Wefirststudiedasystemwith g ¼ 34 . 25 , ¼ 1 = 18 , ¼ 10 p = 72 , T0¼ Tc= 2 p ,and ¼ 10 = 648 ;thisallows comparisonwithprevious 1 þ 1 andsphericalstudiesofa coupledfield-fluidsystemwherethesameparameter choiceswereused [23] .Thetransitioninthiscaseis relativelyweak:intermsof T,theratiobetweenthelatent heatandthetotalthermalenergy,wehave TN¼ 0 . 012 at thenucleationtemperature TN¼ 0 . 86 Tc.Wealsoperformedsimulationswith ¼ 2 = 18 and ¼ 5 = 648 ,for which TN¼ 0 . 10 atthenucleationtemperature TN¼ 0 . 8 Tc,whichwerefertoasanintermediatestrength transition.Wenotethat TN 10 2isgenericforafirstorderelectroweaktransition,whereas TN 10 1would implysometuning [38] . Forthenucleationprocess,wetook ¼ 0 . 0125 Tc, P0¼ 0 . 01 ,and t0¼ tend¼ 2000 T 1 c.Thesimulationvolumeallowedthenucleationof100 – 300bubblessothatthe meanspacingbetweenbubbleswasoforder 100 T 1 c.The wallvelocityiscapturedcorrectly,butthefluidvelocity didnotquitereachthescalingprofilebeforecolliding. Typically,thepeakvelocitypriortocollisionis20% – 30% belowthescalingvalueforthedeflagrations. Fortheweaktransition,wechose ¼ 0 . 1 ,0.2,0.4, and0.6.Thefirstgivesadetonationwithwallspeed vw 0 . 71 ,andtheothersweakdeflagrationswith vw 0 . 44 ,0.24,and0.15,respectively.Theshockprofiles arefoundinFigs.2and3ofRef. [23] ;slicesofthetotal energydensityforoneofoursimulationsareshownin Fig. 1 .Theintermediatetransitionwassimulatedat FIG.1(coloronline).Slicesoffluidenergydensity E=T4 cat t ¼ 400 T 1 c, 800 T 1 c,and 1200 T 1 c,respectively,forthe ¼ 0 . 2 simulation.Theslicescorrespondroughlytotheendofthe nucleationphase,theendoftheinitialcoalescencephase,andthe endofthesimulation. PRL 112, 041301(2014) PHYSICALREVIEWLETTERSweekending 31JANUARY2014041301-2

PAGE 3

¼ 0 . 4 ,forwhichthewallspeedis vw 0 . 44 ,veryclose totheweaktransitionwith ¼ 0 . 2 . Figure 2 (top)showsthetimeevolutionoftwoquantities ¯ Uand ¯ Uf,definedsothat ð ¯ þ ¯ p Þ ¯ U2 ¼ 1 V Z d3x iiand ð ¯ þ ¯ p Þ ¯ U2 f¼ 1 V Z d3x f ii(9) where ¯ and ¯ p arethetime-dependent,volume-averaged rest-frameenergydensityandpressure,respectively. Thesquaresofthesequantitiesgiveanestimateofthe sizeoftheshearstressesofthefieldandthefluidrelativeto thebackgroundfluidenthalpydensity,whereas ¯ Uftends tothermsfluidvelocityfor ¯ Uf 1 .Weseethat ¯ Ugrows anddecayswiththetotalsurfaceareaofthebubblesof thenewphase,whilethemeanfluidvelocitygrowswiththe volumeofthebubblesandthenstaysconstantoncethe bubbleshavemerged.Wehavenoexplicitviscosity,and theslightdecreasingtrendin ¯ Uf,visibleforthe intermediatetransition,arisesfromthewell-knownnumericalviscosityofdonor-celladvection, num ¯ Uf x . Figure 2 (bottom)showstheGWenergydensityscaled bythefinalvalueof ð ¯ þ ¯ p Þ2¯ U4 fandtheaveragebubble sizeatcollision R¼ L=N1 = 3 b,where Nbisthenumberof bubblesinthesimulationvolume.Thescalingenables comparisontoamodeldiscussedaroundEq. (12) ,which predictsalineargrowthin GWatlatetimes,sourcedby persistentperturbationsinthefluid.TheGWenergydensity riseslinearlyafterthebubbleshavefullymergedwithsimilar slopes,whichsupportsthemodel.NotethattheGWsfrom detonations( ¼ 0 . 1 )behavesimilarlytothosefrom deflagrations. InFig. 3 weshowthetimedevelopmentoftheGW powerspectrumastheintermediatestrengthphasetransitionproceeds.Weseethatstronggrowthhappens between t ¼ 600 T 1 cand 1000 T 1 casthebubblesmerge (seeFig. 2 ).For t 1000 T 1 c,thereisevidenceofthe expected k 1powerspectrum,butitbecomeslessclearas theGWpowercontinuestogrow,sourcedbythepersistent fluidperturbations.Attheshortestlengthscales,weseea vw-dependentexponentialfalloff. Toestablishthenatureofthesefluidperturbations,we showinFig. 4 thetimedevelopmentofthelongitudinal (compressional)andtransverse(rotational)componentsof thefluidvelocitypowerspectrum.Atalltimes,itisclear thatmostofthefluidvelocityislongitudinal,indicating thattheperturbationsaremostlycompressionwaves. TurbulencegenerallydevelopsathighReynoldsnumber Reinthetransversecomponents,characterizedbyapowerlawbehaviorofthepowerspectrum.Giventhebubble separationscale R,wecanestimatethevalueofRe,due entirelytothenumericalviscosity,asRenum¼ ¯ UfR= num 102.Thereisnofirmevidenceofapowerlawat 0 500 1000 1500 2000 t ( Tc -1) 0 0.025 0.05 0.075 U Interm., =0.4 Weak, =0.1 Weak, =0.2 Weak, =0.4 Weak, =0.6 Uf U 0 500 1000 1500 2000 t ( Tc -1) 0 200 400 600 800 1000 1200gwR* -1 ( + p ) -2 Uf -4 ( G Tc -1) Interm., =0.4 Weak, =0.1 Weak, =0.2 Weak, =0.4 Weak, =0.6 FIG.2(coloronline).Top:timeseriesof ¯ Uand ¯ Uf(9) , showingtheprogressofthephasetransition;thecurvesfor ¯ Uand ¯ Ufareindividuallyidentifiedforthe “ intermediate ” case. Bottom:timeseriesof GWR 1 ½ð ¯ þ ¯ p Þ 2¯ U 4 ftend,showingthe evolutionofthegravitationalwaveenergydensityrelativetoan estimateofthesquareofthefinalfluidshearstresses. 1 1 . 0 1 0 . 0 k ( Tc) 1 × 10Š91 × 10Š81 × 10Š71 × 10Š61 × 10Š50.0001 0.001 0.01 0.1 1 10d GW/ d ln k ( G Tc 6) k-1 FIG.3(coloronline).Gravitationalwavepowerspectraduring thephasetransition,fortheintermediatestrengthtransition,from fluidonly(blacklines)andbothfluidandfield(graylines).From bottomtotop,thetimesare t ¼ 600 ,800,1000,1200,and 1400 T 1 c.Thereddashedlineindicatestheexpected k 1behavior. PRL 112, 041301(2014) PHYSICALREVIEWLETTERSweekending 31JANUARY2014041301-3

PAGE 4

high k ,butitisunclearwhetherReislargeenoughfor turbulencetodevelophere. WecannowformaclearerpictureofthefluidperturbationsandhowtheGWsaregenerated.First,wenotethat thefluidperturbationsareinitiallyintheformofa compressionwavesurroundingthegrowingbubble.The energyinthiswaveisproportionaltothevolumeofthe bubble R3andquicklyoutstripstheenergyinthescalar field,whichgrowsonlyas R2.Theenergyinthecompressionwavesremainsconstantafterthebubbleshave merged.Thisisduetolinearityandconservationofenergy: asthefluidvelocitiesaregenerallysmall,thereislittle transfertothetransversecomponents. Thebubblecollisiongeneratesgravitationalwaves,as predictedbytheenvelopeapproximation,andthereissome evidenceforthecharacteristic k 1spectrumbetween Randthehigh-frequencycutoff.ThegenerationofGWs continueslongafterthemergeriscompletedandthescalar fieldhasrelaxedtoitsnewequilibriumvalue.TheGWsare sourcedbythecompressionwavesinthefluid.Thissource ofgravitationalradiationfromaphasetransition — sound — hasnotbeenappreciatedbefore(exceptinRef. [4] ). Theresultingdensityofthegravitationalwavesisgiven fromtheunequaltimecorrelatoroftheshearstresstensor 2ð k;t1;t2Þ by [27,39] d GWð k Þ d ln k ¼ 2 Gk3 Ztdt1dt2cos ½ k ð t1 t2Þ 2ð k;t1;t2Þ : (10) Wemodelthesourceasturningonatthenucleationtime tNwithalifetime Ts(discussedbelow)andbeingafunctionof t1 t2betweenthosetimes,asisreasonableforstochastic soundwaves.Wesupposethecorrelatorispeakedat t1 t2¼ 0 withwidth xc=k ,where xcisadimensionless parameter.Thisresemblesthe “ top-hat ” correlatormodel ofRef. [27] ,exceptthatthesourceactsformuchlonger thanthedurationofthetransition 1.Weestimatethe amplitudeofthesourceas ½ð ¯ þ ¯ p Þ ¯ U2 f2anditslengthscale as R.Hence,for tN< ð t1;t2Þ
PAGE 5

fluidvelocities,exploretheeffectofdissipation,andlook forturbulentregimes.Parameterchoicesrecentlyidentified ashavingunstablebubblewalls [41] alsomeritinvestigation.Wehavenotstudiedthecasewherethewallsrun away,althoughhereweexpectthatthefluidisunimportant andtheenvelopeapproximationapplies. Inthecasesthatwedostudy,wefindthevelocity perturbationsareprincipallyacousticwavesandthatthe resultinggravitationalradiationdensityisparametrically largerthangivenintheenvelopeapproximationbytheratio ofthefluiddampingtime stothedurationofthephase transition 1.Weconcludethat,forawiderangeoffirstorderphasetransitionsofinterest,themainsourceofthe gravitationalwavebackgroundisthesoundtheymake. OursimulationsmadeuseoffacilitiesattheFinnish CentreforScientificComputingCSCandtheCOSMOS Consortiumsupercomputer(withintheDiRACFacility jointlyfundedbySTFCandtheLargeFacilitiesCapital FundofBIS).K.R.acknowledgessupportfromthe AcademyofFinlandproject1134018;M.H.andS.H. fromtheScienceandTechnologyFacilitiesCouncil(Grant No.ST/J000477/1). *m.b.hindmarsh@sussex.ac.uk†s.huber@sussex.ac.uk‡kari.rummukainen@helsinki.fi§david.weir@helsinki.fi [1]D.Kirzhnits,JETPLett. 15 ,529(1972). [2]D.KirzhnitsandA.D.Linde, Ann.Phys.(Paris) 101 ,195 (1976) . [3]E.Witten, Phys.Rev.D 30 ,272(1984) . [4]C.J.Hogan,Mon.Not.R.Astron.Soc. 218 ,629(1986). [5]A.Kosowsky,M.S.Turner,andR.Watkins, Phys.Rev.D 45 ,4514(1992) . [6]A.Kosowsky,M.S.Turner,andR.Watkins, Phys.Rev. Lett. 69 ,2026(1992) . [7]A.KosowskyandM.S.Turner, Phys.Rev.D 47 ,4372 (1993) . [8]M.Kamionkowski,A.Kosowsky,andM.S.Turner, Phys. Rev.D 49 ,2837(1994) . [9]K.Kajantie,M.Laine,K.Rummukainen,andM.E. Shaposhnikov, Phys.Rev.Lett. 77 ,2887(1996) . [10]M.LaineandK.Rummukainen, Phys.Rev.Lett. 80 ,5259 (1998) . [11]M.Laine,G.Nardini,andK.Rummukainen, J.Cosmol. Astropart.Phys.01( 2013 )011. [12]M.S.Carena,M.Quiros,andC.Wagner, Phys.Lett.B 380 , 81(1996) . [13]D.Delepine,J.Gerard,R.GonzalezFelipe,andJ.Weyers, Phys.Lett.B 386 ,183(1996) . [14]M.LaineandK.Rummukainen, Nucl.Phys. B535 ,423 (1998) . [15]C.Grojean,G.Servant,andJ.D.Wells, Phys.Rev.D 71 , 036001(2005) . [16]S.J.HuberandM.Schmidt, Nucl.Phys. B606 ,183(2001) . [17]S.J.Huber,T.Konstandin,T.Prokopec,andM.G.Schmidt, Nucl.Phys. B757 ,172(2006) . [18]P.J.Steinhardt, Phys.Rev.D 25 ,2074(1982) . [19]H.Kurki-Suonio, Nucl.Phys. B255 ,231(1985) .[20]K.KajantieandH.Kurki-Suonio, Phys.Rev.D 34 ,1719 (1986) . [21]K.Enqvist,J.Ignatius,K.Kajantie,andK.Rummukainen, Phys.Rev.D 45 ,3415(1992) . [22]J.Ignatius,K.Kajantie,H.Kurki-Suonio,andM.Laine, Phys.Rev.D 49 ,3854(1994) . [23]H.Kurki-SuonioandM.Laine, Phys.Rev.D 54 ,7163 (1996) . [24]H.Kurki-SuonioandM.Laine, Phys.Rev.Lett. 77 ,3951 (1996) . [25]J.R.Espinosa,T.Konstandin,J.M.No,andG.Servant, J.Cosmol.Astropart.Phys.06( 2010 )028. [26]S.J.HuberandT.Konstandin, J.Cosmol.Astropart.Phys. 09( 2008 )022. [27]C.Caprini,R.Durrer,T.Konstandin,andG.Servant, Phys. Rev.D 79 ,083519(2009) . [28]C.Caprini,R.Durrer,andG.Servant, Phys.Rev.D 77 , 124015(2008) . [29]H.L.ChildandJ.T.Giblin,Jr., J.Cosmol.Astropart.Phys. 10( 2012 )001. [30]A.Kosowsky,A.Mack,andT.Kahniashvili, Phys.Rev.D 66 ,024030(2002) . [31]G.Gogoberidze,T.Kahniashvili,andA.Kosowsky, Phys. Rev.D 76 ,083002(2007) . [32]C.CapriniandR.Durrer, Phys.Rev.D 74 ,063521(2006) . [33]C.Caprini,R.Durrer,andG.Servant, J.Cosmol.Astropart. Phys.12( 2009 )024. [34]S.KhlebnikovandI.Tkachev, Phys.Rev.D 56 ,653(1997) . [35]R.EastherandE.A.Lim, J.Cosmol.Astropart.Phys.04 ( 2006 )010. [36]J.Garcia-Bellido,D.G.Figueroa,andA.Sastre, Phys.Rev. D 77 ,043517(2008) . [37]J.WilsonandG.Matthews, RelativisticNumerical Hydrodynamics (CambridgeUniversityPress,Cambridge, England,2003). [38]S.J.HuberandT.Konstandin, J.Cosmol.Astropart.Phys. 05( 2008 )017.[39]D.G.Figueroa,M.Hindmarsh,andJ.Urrestilla, Phys.Rev. Lett. 110 ,101302(2013) . [40]P.B.Arnold,G.D.Moore,andL.G.Yaffe, J.HighEnergy Phys.11( 2000 )001. [41]A.MegevandandF.A.Membiela, arXiv:1311.2453. PRL 112, 041301(2014) PHYSICALREVIEWLETTERSweekending 31JANUARY2014041301-5



PAGE 1

DEVELOPMENT OF CLASSICAL INTERATOMIC POTENTIALS FOR APPLICATIONS IN CORROSION AND PHASE TRANSITIONS By MARK NOORDHOEK A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2014

PAGE 2

© 2014 Mark Noordhoek

PAGE 3

To my parents

PAGE 4

4 ACKNOWLEDGMENTS who rescued my doctoral ambitions mere weeks before the qualifying exam deadline. His guidance as my advisor was truly an enjoyable experien ce due to his knowledge, kindness, and patience. I would like to thank my committee members Dr s . Susan Sinnott, Michele Manuel, Elliot Douglas and Youping Chen. Dr. Sinnott in particular is a special influence for he r classroom teaching, mentoring and our numerous scientific collaborations. My time with Dr s . Igor Levin and Victor Krayzman at NIST was a great experience as they taught me many lessons about what it takes to publish quality work and be a professional scientist. I would like to thank my fellow group members past and present for their shared knowledge, he lpful discussions and camaraderie . Drs. Tao Liang, Alex Chernat y nsk i y and Tzu Ray Shan deserve my praise for their mentoring in numerous instances. Also, I would like to thank m y fellow classmates Tsu Wu Chiang , Bowen Deng and Yangzhong Li who have been there for me throughout our many years in the group. I am forever thankful to my parents for raising me with love and providing me with the wonderful opportunities life. Th eir continued support has allowed me to pursue my professional interests past what I could have envisioned. This dissertation

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 ABSTRACT ................................ ................................ ................................ ................... 12 CHAPTER 1 GENERAL INTRODUCTION ................................ ................................ .................. 14 Projects Overview ................................ ................................ ................................ ... 14 Corrosion Overview ................................ ................................ ................................ 14 Corrosion of Zirconium in Nuclear Reactors ................................ ........................... 17 Crystallography of Zirconium Based Materials ................................ ....................... 18 Zirconium Metal ................................ ................................ ................................ 18 Zirconium Oxide ................................ ................................ ............................... 19 Zirconium Hydride ................................ ................................ ............................ 19 Fluorite Structure ................................ ................................ .............................. 19 Motivation for t his Work on Zr O H ................................ ................................ ......... 20 Overview of Ferroelectricity ................................ ................................ .................... 20 Phonons ................................ ................................ ................................ .................. 22 Perovskite Crystallography ................................ ................................ ..................... 25 Phases of Perovskite Materials ................................ ................................ ........ 26 Phase Transitions in Perovskites ................................ ................................ ..... 28 Motivation f or Studying Ferroelectrics With MD ................................ ...................... 29 2 SIMULATION METHODOLOGY ................................ ................................ ............. 35 Introduction to Computational Methods ................................ ................................ .. 35 Molecular Dynamics Overview ................................ ................................ ................ 35 Interatomic Potential ................................ ................................ ......................... 36 Periodic Boundary Conditions ................................ ................................ .......... 37 Time Integration Algorithms ................................ ................................ .............. 38 Thermodynamics Ensemble Conditions ................................ ........................... 39 Thermostats ................................ ................................ ................................ ..... 40 Barostats ................................ ................................ ................................ .......... 41 Shell Model Overview ................................ ................................ ............................. 42 Shell Mo del Formalism ................................ ................................ ........................... 44 Interatomic Interactions ................................ ................................ .................... 44 Intra atomic Interaction ................................ ................................ ..................... 46 Overview of the Charge Optimized Many Body (COMB) Potential ......................... 46 COMB Potential Formalism ................................ ................................ .................... 47

PAGE 6

6 Parameterization of Interatomic Potentials ................................ ............................. 52 Challenges of Optimization ................................ ................................ ............... 52 Training Database and Cost Function ................................ .............................. 53 Optimization Strategies ................................ ................................ .................... 55 Soft ware Used in t his Dissertation ................................ ................................ .......... 56 3 CHARGE OPTIMIZED MANY BODY (COMB) POTENTIALS FOR ZIRCONIUM .. 58 Overview of Slip in HCP Metals ................................ ................................ .............. 58 Param eterization of COMB Potentials f or Zr ................................ ........................... 60 COMB Static Property Results ................................ ................................ ................ 62 Bulk and Surface Properties ................................ ................................ ............. 62 Point Defects ................................ ................................ ................................ .... 63 Stacking Fault Energy ................................ ................................ ...................... 67 Mechanical Deformation of Zr ................................ ................................ ................. 68 Assessment of COMB Zr Potentials ................................ ................................ ........ 70 4 MECHANISMS OF ZR SURFACE CORROSION DETERMINED VIA MOLECULAR DYNAMICS SIMULATIONS ................................ ........................... 83 Overview of Potentials f or Zr O H ................................ ................................ ........... 83 Parameterization and Results of the Zr O H Potential ................................ ............ 85 Bulk Zr O System ................................ ................................ ............................. 86 Bulk Zr H System ................................ ................................ ............................. 87 Summary of Zr O H COMB Potential ................................ ............................... 88 Energetics of Oxygen and Hydrogen on Zr (0001) Surface ................................ .... 88 Overview of O 2 and H 2 O Deposition ................................ ................................ ....... 91 Deposition of O 2 ................................ ................................ ............................... 92 Deposition of H 2 O ................................ ................................ ............................. 94 Unified Vi ew of O 2 and H 2 O Deposition ................................ ............................ 95 5 ATOMISTIC STRUCTURE OF (BA,SR)TIO 3 : COMPARING MOLECULAR DYNAMICS SIMULATIONS WITH STRUCTURAL MEASUREMENTS ............... 117 Introduction to (Ba,Sr)TiO 3 Ferroelectrics ................................ ............................. 117 Overview of MD Simulations f or BaTiO 3 ................................ ............................... 117 Experimental and Simulatio n Procedure ................................ ............................... 119 BaTiO 3 Results ................................ ................................ ................................ ...... 121 Ba 1 x Sr x TiO 3 Results ................................ ................................ .............................. 125 Conclusions on Confronting MD Simulations w ith Experimental Results .............. 126 6 DEVELOPMENT OF SHELL MODEL POTENTIAL FOR (BA,CA)TIO 3 ................ 137 Approach to Potential Development ................................ ................................ ..... 137 Extension of Ba 1 x Sr x TiO 3 Potential to Ba 1 x Ca x TiO 3 ................................ ............. 137 Development of New Potential f or BaTiO 3 ................................ ............................ 141 Training Database ................................ ................................ .......................... 141 Meth odological Challenges ................................ ................................ ............ 143

PAGE 7

7 Static Properties ................................ ................................ ............................. 144 MD Results at Finite Temperature ................................ ................................ .. 145 Future Work on Shell Model Development ................................ ..................... 146 7 GENERAL CONCLUSIONS ................................ ................................ ................. 162 Summary of Project Outcomes ................................ ................................ ............. 162 Final Thoughts on Interatomic Potential Development ................................ ......... 163 APPENDIX A PARAMETERS DEVELOPED FOR COMB POTENTIAL ................................ ..... 166 B PARAMETERS DEVELOPED FOR SHELL MODEL ................................ ............ 170 LIST OF REFERENCES ................................ ................................ ............................. 175 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 185

PAGE 8

8 LIST OF TABLES Table page 3 1 Slip systems in zirconium ................................ ................................ ................... 72 3 2 Properties of zirconium metal ................................ ................................ ............. 73 3 3 Zr interstitial formation energies in DFT and COMB ................................ ........... 74 3 4 Divacancy formation and binding energies ................................ ......................... 75 3 5 Stacking fault energies of zirconium metal ................................ ......................... 76 4 1 Properties of cubic fluorite ZrO 2 ................................ ................................ .......... 97 4 2 Point defect formation energies in zirconia ................................ ......................... 98 4 3 Properties of zirconium hydride ................................ ................................ .......... 99 4 4 Properties of cubic fluorite ZrH 2 ................................ ................................ ........ 100 4 5 Binding energies of oxygen atoms o n Zr (0001) surface ................................ .. 101 4 6 Binding energies of hydr ogen atoms on Zr (0001) surface ............................... 102 5 1 Experimental and simulated lattice parameters for Ba 1 x Sr x TiO 3 ...................... 129 5 2 Agreement between the selected simulated and experimental structural characteristics in Ba 1 x Sr x TiO 3 ................................ ................................ .......... 130 6 1 Lattice parameters and internal coordinate positions for orthorhombic CaTiO 3 147 6 2 Experiment al and MD structural properties for the rhombohedral phase of BaTiO 3 . ................................ ................................ ................................ ............. 148 6 3 Experimental and MD structural properties for the or thorhombic phase of BaTiO 3 . ................................ ................................ ................................ ............. 149 6 4 Experimental and MD structural properties for the tetragonal phase of BaTiO 3 . ................................ ................................ ................................ ............. 150 6 5 Experimental and MD structural properties for the cubic phase of BaTiO 3 . ...... 151 6 6 Relative energy differences of BaTiO 3 ................................ .............................. 152 6 7 Energy profile for Ti displacements in cubic BaTiO 3 ................................ ......... 153 A 1 COMB p otential parameters for Zr ................................ ................................ .... 166

PAGE 9

9 A 2 Short range Zr Zr parameters for COMB Zr 1 and COMB Zr 2. ....................... 167 A 3 Potential parameters for Zr O and Zr H using COMB Zr 2 ............................... 168 A 4 COMB Zr 2 potential parameters for three body Legendre polynomials .......... 169 B 1 Two body parameters for the (Ba,Ca)TiO 3 potential ................................ ......... 170 B 2 Elemental parameters for the (Ba,Ca)TiO 3 potential ................................ ........ 171 B 3 Elemental parameters for BaTiO 3 in the new fitting attempt ............................. 172 B 4 Three body parameters for BaTiO 3 in the new fitting attempt ........................... 173 B 5 Two body parameters for BaTiO 3 in the new fitting attempt ............................. 174

PAGE 10

10 LIST OF FIGURES Figure page 1 1 A unit cell of the FCl 2 fluorite structure ................................ ............................... 31 1 2 Sche matic of a ferroelectric hysteresis loop of polarization as a function of applied electric field ................................ ................................ ............................ 32 1 3 A unit cell of the ABO 3 perovskite structure ................................ ........................ 33 1 4 A unit cell of tetragonal (I 4 /mcm) SrTiO 3 ................................ ............................. 34 2 1 Schematic of bond distances and ang le between three atoms ........................... 57 3 1 Common slip systems found in HCP metals ................................ ....................... 77 3 2 Intersti tial positions in a HCP lattice ................................ ................................ ... 78 3 3 Schematic of divacancy positions in a HCP lattice ................................ ............. 79 3 4 Stacking fault maps for basal (0001) in Zr ................................ .......................... 80 3 5 Stacking fault maps for prism in Zr ................................ .......................... 81 3 6 The structure of nanocrystalline Zr after a tensile test with 8.8% strain .............. 82 4 1 A schematic of the Zr (0001) surface ................................ ................................ 103 4 2 Snapshot of Zr (0001) surface viewed along the direction after 50 ps at 300 K ................................ ................................ ................................ ................ 104 4 3 Snapshots of Zr surfaces after O 2 deposition ................................ ................... 105 4 4 Dissociation of an O 2 molecule on Zr (0001) ................................ .................... 106 4 5 Oxygen profile in Zr as a function of depth after deposition of O 2 ..................... 107 4 6 Fraction of oxygen atoms on th e surface and within the first subsurface Zr layer after O 2 deposition ................................ ................................ ................... 108 4 7 Fraction of Zr atoms in the top two planes that have no nearby oxygen atom after O 2 deposition ................................ ................................ ............................ 109 4 8 Bond angle distribution in Zr (0001) after O 2 deposition ................................ ... 110 4 9 Snapshots of Zr surfaces after H 2 O deposition ................................ ................. 111 4 10 Dissociation of a H 2 O molecule on Zr (0001) ................................ .................... 112

PAGE 11

11 4 11 Fraction of atomic oxygen in the Zr structure after H 2 O deposition .................. 113 4 12 Fraction of atomi c hydrogen in the Zr structure after H 2 O deposition ............... 114 4 13 Fraction of atomic oxygen in the two outermost Zr planes after H 2 O deposition ................................ ................................ ................................ ......... 11 5 4 14 Fraction of atomic hydrogen in the two outermost Zr planes after H 2 O deposition ................................ ................................ ................................ ......... 116 5 1 Comparison of various experimental and simulated structural properties of BaTiO 3 at 300 K ................................ ................................ ................................ 131 5 2 The Ti PDD in tetragonal BaTiO 3 at 300 K ................................ ....................... 132 5 3 Local p z polarization component in MD as a function of time for BaTiO 3 at 350 K. ................................ ................................ ................................ ............... 133 5 4 Simulated Ti Ti Ti chain data in BaTiO 3 ................................ ........................ 134 5 5 Comparison of experimental and calculated (MD) data for Ba 1 x Sr x TiO 3 .......... 135 5 6 Bond length and Ti off centering data for simulated Ba 1 x Sr x TiO 3 ..................... 136 6 1 Experimental and calculated neutron PDFs for CaTiO 3 at 300 K ..................... 154 6 2 Comparison of various experimental and simulated structural properties of CaTiO 3 and Ba 0.8 Ca 0.2 TiO 3 at 300 K ................................ ................................ .. 155 6 3 Experimental an d calculated neutron PDFs for Ba 0.8 Ca 0.2 TiO 3 . ........................ 156 6 4 PDD for Ca in simulated Ba 0.8 Ca 0.2 TiO 3 and Sr in simulated Ba 0.8 Sr 0.2 TiO 3 ..... 157 6 5 The c/a ratio as function of Sr/Ca composition for BaTiO 3 , Ba 1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 ................................ ................................ ................................ ..... 158 6 6 Tetragonal < > cubic transition temperature for BaTiO 3 , Ba 0.8 Ca 0.2 TiO 3 and Ba 0.8 Sr 0.2 TiO 3 ................................ ................................ ................................ .... 159 6 7 Oxygen oxygen bond distance in BaTiO 3 as a function of temperature ........... 160 6 8 Average BaTiO 3 lattice parameters at 300 K ................................ .................... 161

PAGE 12

12 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DEVELOPMENT OF CLASSICAL INTERATOMIC POTENTIALS FOR APPLICATIONS IN CORROSION AND PHASE TRANSITIONS By Mark Noordhoek August 2014 Chair: Simon R. Phillpot Major: Materials Science and Engineering Classical molecular dynamics (MD) is an established method used to simulate atomic scale phenomena in materials systems. The key component of MD is the description of the interaction energy between atoms. In this work, two independent formalisms of classical interatomic potentials are used to describe various aspects of zirconium based and barium titanate based systems. A third generation Charge Optimized Many B ody (COMB) potential fo r the zirconium zirconium oxide zirconium hydride (Zr O H) system is presented . These materials are widely found in the fuel cladding components of nuclear reactors. During energy production , Zr undergoes a number of mechanical and electrochemical transfor mations, which may lead to catastrophic failure if allowed to continue indefinitely. Thus, a better understanding of these processes is of critical importance. COMB allows various mechanisms of these transformations to be studied within the same simulation . The metallic Zr interactions in this work correctly describe many structural properties including stacking fault energies and surface energies. Mechanical testing of Zr shows that the slip mechanisms in MD match those in experiments. Deposition of O 2 and H 2 O molecules onto various low index Zr surfaces shows the

PAGE 13

13 resulting corrosion behavior is dependent on the surface structure and deposited species. J oint MD and experimental studies on the structural properties of the Ba 1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 s ystem s are presented. These system s are ferroelectric for small Sr /Ca concentrations (< 30%), which is exploited in their use for many commercial devices such as capacitors. The mechanism of the phase transition between the ferroelectric tetragonal phase a nd paraelectric cubic phase is still an open debate. In this work, the shell model is used to describe interatomic interactions for each system. Direct comparison between the structural properties predicted in MD and those observed in experiment is provide d. The model significantly overestimates the distortion within the tetragonal phase, which prevents a reliable determination of the transformation mechanism. Additional work to develop a new shell model description for BaTiO 3 is presented. However, efforts to produce an improved potential have not been successful. Suggestions for improved potentials are provided.

PAGE 14

14 CHAPTER 1 GENERAL INTRODUCTION Projects Overview This dissertation contains information and results for two distinct projects. The first project covers various aspects of zirconium, zirconium oxide and zirconium hydride . One area of interest for these materials is in nuclear reactors, where zirconium based clad ding undergoes a number of mechanical and elect rochemical processes during the process of energy generation . Featured topics in this dissertation include the mechanical deformation of metallic Zr and corrosion of various Zr surfaces when exposed to gaseous O 2 and H 2 O. T his chapter provides the relevant background on corrosion and crystallography for these systems in order to better understand the subjects in later chapters . The second project covers structure property relationships in barium titanate based ferroelectrics. Barium titanate is widely used in a variety of devices such as capacitors and therm al switches . Its ferroelectric properties are determined by the time dependent atomic structure, which is the focus of this project. Thus, this chapter also provides background on ferroelectricity, crystallography and the motion of atoms in a crystal. Corrosion Overview Corrosion of metals is degradation caused by an electrochemical reaction with the environment. Common examples include the formation of rust on steel and the surface tarnish on silverware. Corrosion is a major economic problem , costing approximately 3% of gross domestic product in industrialized countries [1] . Thus, the

PAGE 15

15 need to fully understand corrosion mechanism s and reduce the ir adverse e ffects is of great societal import. The electrochemical process of corrosion involves two half cell reaction s that involve the transfer of charge and mass. The two half cell reactions are oxidation and reduction, which are the loss and gain of electrons, respectively. The ability to undergo an oxidation and reduction (redox) reaction leads to a potential energy gradient across the system . The cell potential ener gy (also called the electromotive force) of the system is (1 1) where E Cathode is the reduction half reaction and E Anode is the oxidation half reaction [2] . This work will focus on the corrosion of zirconium. The half cell reactions of zirconium interacting with hydrogen are (1 2) (1 3) where the superscript on E 0 implies standard conditions. The potential associated with the cathode reaction by definition is zero, as all reaction energies are given relative to hydrogen. Thus, t he total cell reaction (eq. 1 1) is +1.55 V. The ability of a reaction to spontaneously occur is determined by the Gibbs free energy (1 4 ) where n is number of electrons transferred, E Cell is the potential energy and F is [3] . A negative free energy corres ponds to a spontaneous reaction , while a positive free energy need s an external driving force t o induce the reaction .

PAGE 16

1 6 T he above discussion about half cell reactions only applies to a metal immersed in a solution with its own ions. This rarely happens for m ost corrosion processes . Combining the free energy discussion above with the definition of the Gibbs free energy, one can show that (1 5 ) w here R is the gas constant, T is the temperature and Q is the reaction quotient [3] . Equation 1 5 is commonly known as the Ner n st equation ; this equation is most useful in solving practical problem s in electrochemistry . Thermodynamics determines whether a reaction may occur, but the rate of reaction is highly dependent on the metal . The oxidation (or hydrogenation) of a metal surface results in the growth of an oxide film, which is often proportional to the time exposed to the corrosive environment. T he kinetics of reaction is based on the ability of the oxide to form a protect ive coating and the mechanism of oxide growth. The ability to form a protective oxide is quantified by the Pilling Bedworth ratio, defined as (1 6) where a value of unit y means that the oxide film covers the metal with no stress [2] . Commonly observed rate laws are cubic, parabolic and linear, which correspond to an oxide thickness that grows with time as t 1/3 , t 1/2 and t, respectively. Metal systems with X PB less than unity tend to have linear growth, while system s with X PB slightly greater than unity tend to have protective films that grow at a parabolic rate , t 1/2 [2] . Metals used in corrosive environments often have a cubic or parabolic growth rate. The mechanism of film growth is determined by the location of the oxide growth

PAGE 17

17 [2] . In some systems, metal cations and electrons diffuse to the outer surface where the new oxide is formed. Other systems have oxide films that grow at the metal/oxide interface, which requires diffusion of oxygen anions to the interface and diffusion of electrons to the outer oxide surface. This type of oxide growth generally experiences a parabolic to linear transition in the growth rate, due to cracks or pores forming at the interface. Corrosion of Zirconium in Nuclear Reactors Zirconium a lloys are commonly used for fuel claddings in nuclear reactors. These alloys were selected due to the low thermal neutron cross section of Zr, their good corrosion resistance and good mechanical properties [4] . Zircaloy 4 (typically Zr 1.5% Sn 0.2% Fe 0.1% Cr) is the current industry standard, while ZIRLO TM (Zr 1% Sn 1% Nb) is an emerging Westinghouse developed material. During the course of the fuel cycle, the fuel cladding undergoes a number of mechanical and chemical processe s. Current challenges that limit cladding performance include grid to rod fretting (GTRF), pellet cladding intera ctions (PCI) and corrosion [5] . Corrosion mechanisms of the Zr alloys are influenced by the chemistry of the matrix, microstruct ure of the surface, temperature and corrosion environment. Introduction of oxygen and/or hydrogen into hexagonal close packed (HCP) Zr decreases ductility, which causes a decrease in the frac ture toughness, impact strength and makes it more susceptible to brittle fracture [6,7] . The maximum solubility of oxygen in Zr is relatively high with experimental values ranging from 29 40% [8 11] . Oxygen occupies the octahedral interstitial position in the matrix, with long range order whose structure is highly dependent of the composition a nd heat treatment. An increase in the oxygen concentration results in the formation of a zirconia ( ZrO 2 ) layer. The presence of

PAGE 18

18 hydrogen results in a zirconium hydrogen solid solution or zirconium hydride precipitation [12] . The hydride morphology and relative orientation to the stress field are known to affect the degree of embrittlement [7] . The Pilling Bedworth ratio f or zirconium and its oxide is ~1.56 [4] . The oxide initially grows at a cubic rate with new oxide forming at the metal/oxide interface. The growth rate changes to linear at an oxide thickness of approximately 2 2.4 m [4] . The rate limiting step of oxide growth is the diffusion of electrons to the outer oxide surface because zirconia is an excellent electrical insulator. A complete understanding of the initial oxidation structure, transition from ZrO x to ZrO 2 and oxide grain formation/growth has been difficult to achieve within these different nuclear reactor environments. On e aim of this work is to begin to develop a model capable of describing the corrosion of Zr in conditions similar to those in nuclear reactor s . While a complete description of oxide growth is limited due to the exclusion of electrons and their transport in this model , individual aspects of the process may still be studied in detail. Having a comprehensive model framework capable of describing Zr O H in the same simulation is a current advancement in the field. Crystallography of Zirconium Based Materials Zirconium Metal The ground phase of Zr has an HCP structure with space group P6 3 /mmc . At 1136 K it transforms into body centered cubic (BCC; space group ) before melting at 2128 K [13] . At high pressure , the HCP structure transforms into a nother hexagonal phase with space group P6/mmm [14] .

PAGE 19

19 Zirconium Oxide At zero pressure, the ground state of zirconia is a mon oclinic phase ( space group P2 1 /c) called baddeleyite i n mineralogy. The cations are seven fold coordinated while the anions are 3 fold and 4 fold coordinated. Upon heating, it transforms to tetragonal ( space group P4 2 /nmc; T c ~ 1440 K) and cubic ( space gro up ; T c ~ 2640 K) phases [15] . In t he tetragonal and cubic phases , cations are 8 fold coordinated while anions are 4 fold coordinated. The cubic structure is isostructural to fluorite, CaF 2 . The high pressure phases are Ortho rhombic I ( space group Pbca) with 7 fold cation coordination and Ortho rhombic II ( space group Pnma) which has 9 fold cation coordinati on . Ortho rhombic II is isostructural to the cotunnite PbCl 2 structure. Zirconium Hydride Currently, four different zirconium hydride phases of various stoichiometries have been identified [16] . The metastable Zr 2 H ( space group P3m1) is coherent with the Zr matrix. With added hydrogen a metastable ZrH ( space group P4 2 /n) phase is observed, which is face cen tered tetragonal (c/a > 1). A face centered cubic phase is designated the ZrH 1.5 phase ( space group ). The phase has H atoms randomly occupying only 6 of the 8 tetrahedral sites. The ZrH 2 phase is face centered tetragonal (space group I4/mmm) , ZrH 2 , with c/a ratio less than 1. At standard conditions, ZrH 2 is the most energetically favorable phase. Fluorite Structure The fluorite structure (Figure 1 1 ) is in the cubic crystal system and has a space group of . The cations occupy Wyckoff position 4a while anions are at 8b. The structure may be thought as having simple cubic (SC) anion lattice within a cation face -

PAGE 20

20 centered cubic ( FCC ) network. The cations are 8 fold coordinated while the anions are 4 fold coordinated. Motivation f or t his W ork o n Zr O H The numerous complex processes that take place in the fuel cladding are typically difficult to separate experimentally. Simulation provides one vehicle to study the system as a whole, as well as to parse out individual unit mechanisms. This r equires a framework capable of handling each individual material system (metal, metal oxide, metal hydride), in addition to their integration. Molecular dynamics (MD) simulations of zirconium and other metals commonly use the embedded atom method (EAM) or its modified derivatives [17 19] . While these potentials provide a satisfactory characterization of many important properties in zirconium, their formalisms are incapable of describing nonmetallic environments and are thus not able to attack the issues of PCI or corrosion. The Charge Optimized Many Body (COMB) formalism is capable of providing a unified description of the interatomic interactions in systems with mixed bonding [20 22] . In this work, a COMB description of the Zr Z rO 2 ZrH 2 system is provided. A comparison of various calculated properties of z irconium metal, zirconium oxide and zirconium hydride to experimental and first principles results is made . Simulations of O 2 and H 2 O deposition on different Zr surfaces provide insight into corrosion behavior and form the grounds for further study. Overview of Ferroelectricity Ferroelectricity is the ability of a material to undergo a reversal of a spontaneous electric polarization in response to an applied electric field. The a toms in f erroelectric materials experience small displacements that lead to a reduction in symmetry [23] . In

PAGE 21

21 crystallography, there are 32 point groups that are compatible with three dimensional translation of atoms to fill space [24] . Twenty one of these point groups are non centrosymmetric. Twenty of these 21 point groups may be piezoelectric, which is the ability of generate electricity in response to a mechanical stress. Ten of these 20 piezoelectric point groups are pyroelectric, which produces polarization from a change in temperature. Pyroelectric and ferroelectric materials experience charge separation even in the absence of an electric field, and are thus called polar ma terials. F erroelectric materials are always piezoelectric and pyroelectric , but not all piezoelectric and pyroelectric materials are ferroelectric . The high symmetry structure without structural distortions is called paraelectric. The chemical bonding in f erroelectric materials is mostly ionic, which is characterized by electrostatic interactions between ions. These ions contain positively charged nuclei surrounded by negatively charged electron clouds. When an electric field is applied to a ferroelectric m aterial , a dipole moment is generated by the separation of the nuclei and electron clouds [25] . The electronic polarizability tensor , , is defined as (1 7 ) wh ere, is the dipole moment and is the electric field. In addition to electronic polarization, i onic polarization may contribute to the overall polarization mechanism . To illustrate ion ic polarization, c onsider a centrosymmetric structure such as rock salt (NaCl). In the absence of an electric field, the ions are located at thei r equilibrium lattice positions and the net dipole moment is zero. When an electric field is applied, the cations and anions are displace d in opposite directions. For crystal structures where the

PAGE 22

22 cation and anion displacements may be unequal, a net dipole moment with finite value is created [25] . The total polarization of a material is defined as the sum of all dipole moments (electronic , ionic , others ) per unit volume [25] . The polarization is correlated on length scales larger than the unit cell, which leads to the formation of domains in which the dipoles are aligned. These domains grow and shrink as the electric field is changed. As the electric fiel d is increased, the polarization reach es a maximum value , P s . The process of reversing the polarization in an electric field creates a hysteresis loop, which is shown schematically in Figure 1 2. Reducing the electric field to zero will not cause the polarization to be zero; instead, a remnant polarization (P r ) will remain. Continuing the electric field reversal leads to zero polarization; the value at which this occurs is called the coercive field (E c ). Phonons All atoms in a crystal vibrate around their equilibrium positions . The minimum energy of a crystal , calle d the zero point energy, occurs a t zero Kelvin [25] . Since atoms are bonded to one another their vibrations are coupled together. Phono ns are the quant a of lattice vibrational energy in which the atomic vibrations move in a wave like manner [26] . Vibrational modes where atoms move parallel to the wave propagation direction are longitudinal modes, while atomic motions perpendicular to the wave propagation are t ra ns verse waves. The number of vibrational mo des in a crystal is determined by the atomic structure . A unit cell of N atoms will have 3 acoustic modes and 3N 3 optical modes. Therefore, 3N is the maximum number of independent modes in a crystal [25] .

PAGE 23

23 The number of activated phonons in a crystal is dependent on temperature. As the temperature increase s , the internal energy of the system at constant volume also increase s proportional to the molar heat capacity, C v . The temperature at which all vibrational modes are excited is known as the Debye temperature, T D . Above T D , nearly all crystals have heat capacities close to , but not exactly, C v =3R, where R is the gas constan t [25] . Near the phase transition temperature, the heat capacity may exhibit a large increase which is caused by the anharmonic vibrati on of atoms [27] . The potential energy based on the atomic positions of the atoms may be approximated by a Taylor expansion around the equilibrium positions: (1 8 ) where l, l denote [28] . The 0 and is the displacement of atom l equili brium , the second term on the righ t hand side of eq. 1 8 vanishes since the force on the atoms is zero. I f higher order terms (greater than second order ) are neglected eq. 1 8 becomes the harmonic approx imation to the potential energy where (1 9 ) directions at equilibrium . (1 10) In a crystal, the atomic positions are repeated indefinitely along the three Cartesian coordinates. Bloch theorem provides a wave like solution to eq. 1 10 as:

PAGE 24

24 (1 11) Here, N is the number of lattice vectors and m is the mass of atom l [28] . The summation over k represents the allowable energy eigen states within the first Brillouin zone. Normal modes , Q k , are constant coefficients and polarization vectors are k . Substituting eq. 1 11 into 1 10 yields (1 12) where (1 13) is known as the dynamical matrix [28] . The allowable vibrational frequencies are those that solve the equation : (1 14) When an external stimulus is applied to a material (temperature gradient, stress, e tc .), atomic displacements away from their original positions may decrease the energy of the system, thereby creating a driving force f or the phase transition. Using the general one dimensional time depend ent wave solution of (1 15 ) along with eq. 1 10 results in : (1 16 ) Observation of eq. 1 16 shows that the wave amplitude grows exponentially when takes on imaginary values . A phonon mode where the frequency tends to imaginary

PAGE 25

25 values due a stimulus is called a soft phonon. The traverse optical phonons are responsible for displacive type phase transitions in ferroelectric materials, which will be discussed in a later section [29] . Perovskite C rystallography The focus of this work on ferroelectrics is a small subset of pero v skites with the chemical formula ABO 3 . Figure 1 3 shows the unit cell of the pero v skite structure. This structure has space group , which is in the cubic crystal system . The A site atoms are a t the corners of the cube (Wyckoff position 1a) and have 12 nearest oxygen neighbors. The B site atom is located in the center of the cube (Wyckoff position 1b) and has 6 nearest oxygen n eighbors that form an octahedron . Oxygen atoms lie at the face center s (Wyckoff position 3c) and have 4 A site and 2 B site neighbors. When the perovskite structure is replicated in three dimensions it forms a network of octahedr a with the oxygen atoms as the vertices. In order for the system to be mechanically stable, th e A site and B site atoms must fit within this network. The stability of perovskites can be estimated by a tolerance factor [23] (1 17 ) where r X is the atomic radius of atom located at the X site . A value of equals 1 represents an ideal perovskite ; however, most real materials deviate from this based on the r A /r B ratio . When >1 , the A site atoms are much larger than B site and a polar distortion occurs , such as that in barium titanate [23] . When < 1 , oxygen octahedra rotations and tilting will occur like those observed in strontium titanate and calc ium titanate. Generally, the A site cations are larger than the B site cations. The value of is

PAGE 26

26 a useful quantity because it can by used to find region of stability for many chemically unique perovskites [30] . Phases of Perovskite Materials This work will focus on barium titanate ( BaTiO 3 ) , strontium titanate ( SrTiO 3 ) , calcium titanate ( CaTiO 3 ) and their solid solutions. The high temperature paraelectric phase for each of these materials is the perovsk ite structure (space group ). The low temperature phases in each differ due to the stability criteria described in the previous section. Upon decreasing the temperature in BaTiO 3 , the cubic phase transforms successively into three ferroelectric phases of tetragonal (T , space group P4mm), orthorhombic (O , space group Amm2) and rhombohedral (R , space group R3m) symmetry [31] . The C T, T O and O R transitio n temperatures are 393 K, 278 K and 183 K, respectively [31] . Each of the three ferroelectric phases may be described as elongations of the cubic unit cell; distortions along [001], [ 011] and [111] occur for the T, O and R phases, respectively. The low temperature (T c = 10 5 K) structure of SrTiO 3 is tetragonal (space group I4/mcm) [32] . This tetragonal cell is shown in Figure 1 4; this repre sentation has four times the number of atoms as the high temperature cubic unit cell. The Ti pos itions occupy Wyckoff position 4 c (0, 0, 0), Sr is at Wyckoff position 4b (0, 0.5, 0.25), while the oxygen atoms occupy two different positions; t he first is Wy ckoff position 4a (0, 0, 0.25) and the second is 8h ( x+0.5, x, 0.25) where x=0.24. In this tetragonal phase, adjacent oxygen octahedra in the (100) plane undergo an antiferrodistortion (AFD) , which is characterized by rotations of the octahedra cage (Figu re 1 4) . This rotation is 2.1 o at 4.2 K [33] .

PAGE 27

27 The ground state of CaTiO 3 is orthorhombic (space group Pbnm), which transforms into tetragonal (space group I4/mcm) at 1423 K, with a transformation to cubic at 1523 K [34] . The orthorhombic phase Pbnm is sometimes equivalently labeled Pnma : w hen the c direction is the long axis the space group is Pbnm; if the b direction is the long axis, then it is Pnma [30] . The Pbnm phase has 20 atoms per unit cell and is related to the original perovskite cell by a 45 o rotation about <001> [30] . Using a cubic perovskite reference, o xygen octahedra in Pbnm experience tilting denoted a + b b in Glazer notation [35] . abb cubic [100] direction being a different magnitude than in [010] and [001]. A + indicates they have opposite tilt. Substituting the A site in the perovskite structure with a different atom species will form solid solutions in the BaTiO 3 SrTiO 3 CaTiO 3 system . B site substitutions , typically Zr and Hf, are possible, but not studied in this work. S olid solution s for BaTiO 3 SrTiO 3 and SrTiO 3 CaTiO 3 are completely miscible for all compositions. The a ddition of small amounts (< 3 0%) of Sr or Ca into BaTiO 3 at room temperature affects the ferroelectric properties of the tetragonal phase. Substituting Sr into BaTiO 3 reduces the c/a distortion and decreases the ferroelectric transit ion temperature [36] . In contrast, the inclusion of Ca decrease s the c/a distortion, but the transition temperature remains unchanged. The large difference in atom ic radii between Ba (1.35 Ã… ) and Ca (1.00 Ã… ) results in strained Ca O bonds and amplified Ti off centering caused by Ca displacements [37,38] . The smaller size difference between Ba (1.35 Ã… ) and Sr (1.18 Ã… ) results in the Sr O bonds being less strained than that of Ca O [38] . As a result, the Ti

PAGE 28

28 off centering in Ba 1 x Sr x TiO 3 is reduced compared to BaTiO 3 and this leads to a lower transition temperature. Phase Transitions in Perovskites Two widely used models to describe phase transitions in per ovskites are the displacive and order disorder models. A displacive transition is characterized by a soft mode phonon that tend s towards zero frequency as the temperature is lowered [29] . Near the transition temperature, T c , the traverse lattice vibrations, T , are no longer harmonic and the conditions for the transition are (2 18 ) where A is a coefficient dependent on temperature [29] . In BaTiO 3 , this instability occurs between Ti and the displacements o f the neighboring oxygen atoms. In the order disorder model, at any instant the B site atom occup ies one of e ight equivalent sites displaced along the [111] direction within the cubic phase [39,40] . Electron diffraction patterns in BaTiO 3 show diffuse sheets in the cubic, tetragonal and orthorhombic phases, but not in the rhombohedral phase [39] . The cubic phase has three equivalent [100 ] directi ons, tetragonal phas e has two, orthorhombic has one and rhombohedral has zero. Thus, the observed number of diffuse sheets is the s ame as the number of unique [100 ] directions in each phase, which corresponds to the available sites for Ti O Ti chain disp lacements. There is experimental evidence supporting both displacive and order disorder models in BaTiO 3 . A n x ray absorption fine structure (XAFS) study showed Ti displacements along [111] in the cubic phase [41] . The magnitude of these displacements did not change as the temperature was lowered to the tetragonal phase ;

PAGE 29

29 this is taken as evidence of order disorder behavior. Nuclear magnetic resonance (NMR) spectroscopy show ed Ti displacements along the [100] , with t he local tetragonal distortions having no long range order [42] ; this is taken as evidence of displacive behavior . The discre pancy could be caused by the different observation time of the two techniques; NMR has a resolution of ~ 10 8 s, while for XAFS the time resolution is ~10 15 s [43] . Since atomic vibration frequencies are close to 10 1 3 s , XAFS can be taken to correspond to the instantaneous structure, while the NMR probes a time averaged structure. Thus, it is widely accepted that the mechanism for phase transitions in BaTiO 3 likely lies between the displacive and order disorder models [ 43 45] . In contrast to B aTiO 3 , SrTiO 3 and CaTiO 3 do not have polar phonons that tend to imaginary values as the temperature is decreased. The phonon stabilization is thought to be due to the effect of quantum fluctuations at low temperatures ( < 5 K) [46] . Instead of the formation of ferroelectric phases, SrTiO 3 and CaTiO 3 undergo rotation and tilting of the oxygen octahedra . The low temperature phases of SrTiO 3 and CaTiO 3 are not ferroelectric. However, doping their ground state structures may result in the formation of a ferroelectric phase; therefore, these materials are known as incipient ferroelectrics. Dopant i nduced ferroelectricity in SrTiO 3 or CaTiO 3 can occur at temperatures where quantum effects can be neglected, which allo ws these materials to be studied using classical interatomic potentials. Motivation f or S tudying F erroelectrics W ith MD Molecular d ynam ics (MD) simulation is well suited to study phase transitions in ferroelectrics due to its ability to simulate materials at finite temperatures and access small timescales . Several MD studies on the phase transitions in BaTiO 3 have been performed over the years [47 51] . Howeve r, it is difficult to directly compare their results

PAGE 30

30 due to limited reporting of structure properties , such as lattice parameters or bulk properties. The goal of this work is to develop interatomic potential s for the BaTiO 3 SrTiO 3 CaTiO 3 system and their solid solutions. An important consideration in developing these potentials is correctly describing the distortions along the Ti O Ti O chains as these will greatly affect the resulting phase transitions. The focus of this work will be on the tetragonal cubic transition in BaTiO 3 , as this occurs within the range of technological importance (393 K). The effect of Sr/Ca substitution for Ba on the transition temperature and structural properties will be discussed.

PAGE 31

31 Figure 1 1. A unit cell of the FCl 2 f luorite structure . The blue spheres are the cations and the red spheres are anions.

PAGE 32

32 Figure 1 2. Schematic of a ferroelectric hysteresis loop of polarization as a function of applied electric field. The figure is reproduced with permission from Ref. [52] .

PAGE 33

33 Figure 1 3 . A unit cell of the ABO 3 p erovskite structure . Black spheres are A site cations, the blue sphere is a B site cation and the red spheres are oxygen.

PAGE 34

34 A B Figure 1 4. A unit cell of tetragonal ( I 4 /mcm ) SrTiO 3 . A) The unit cell viewed along the [110] direction. B) The unit cell viewed along the [001] direction . The sphere colors for Ti, O and Sr are blue, red and black , respectively.

PAGE 35

35 CHAPTER 2 SIMULATION METHODOLOGY Introduction to Computational Methods There are a number of computational meth odologies used to study materials. These methods may be broa dly classified as atomistic, mesoscale and continuum level modeling. The difference between them is the length and time scales accessible in each method. Typically, atomistic methods probe physical phenomena that occur on the same order of magnitude as n anometers or picoseconds, while for mesoscale modeling millimeters or seconds are the largest scale observed . Continuum level models, such as the finite element method (FEM), are capable of much longer length and time scales (meters and days). Various meth ods are used to bridge the gap in length and time scales, but the limits of modern computers prohibit the calculation of many phenomena across the various scales. In this dissertation , atomistic simulations based on first principles calculations and mole cular dynamics (MD) are used to study various material systems. Molecular dynamics in particular is extensively used throughout this dissertation. This chapter provides the foundations of MD and the methods required to develop a model for different materia ls. Molecular Dynamics Overview Molecular d ynamics (MD) is a simulation method that produces the trajectory of interacting atoms by solving equations of motion [53] . A complete description of atomic motion would requi re the calculation of all the many body interactions among all the elec trons and nuclei in the system, which would be a very time consuming calculation. Fortunately, the Born Oppenheimer approximation states that atomic motion can be

PAGE 36

36 described with only nu clei interactions . Electronic structure calculations exploit the Born Oppenheimer approximation by ascribing electron positions based on those of the nuclei. This is due to the large mass difference between the nuclei and electron clouds , such that electro n motion may be regarded as instantaneous in relation to that of the nuclei [53] . Molecular dynamics and other computational methods use this approximation to simplify the resulting calculations . the equation of motion in classical mechanics via: (2 1) where, is the force acting on atom i computed from its mass, , and acceleration . In principle, the motion of atoms is completely deterministic if given an initial set of positions and velocities. Interatomic Potential In density functional theory (DFT) methods, the energy of a system is described by a fun ctional of the electron den sity where e lectrons and their many body interactions are explicitly defined. Atomic motion in a b initio MD is achieved via calculation of the potential energy and use of classical mechanics. Using approximations from first principles, many analytic forms for the potential energy have been devised . However, the form of some analytical functions is simply assumed ad hoc . These functions are far less computationally time consuming in MD , as electrons and their intera ctions are excluded. Collectively these ana lytical forms are (2 2)

PAGE 37

37 wh ere, the force, , is the negative gradient of the energy function, V . The energy of the system is determined from the atomic configuration based on the positions of the nuclei, since this description averages the effective contributions from the electrons. The choice of potential form is critically important and is generally chosen based on the type of bonding found with in the system. One of the earliest interatomic potential s is the Lennard Jones potential, which accurately describes interactions in noble gases [54] . The f unctional form may be written as (2 3 ) where r is the distance between two atoms and separation of r 0 . More advanced potential forms have been developed based on the physics of particular types of atomi c bonding. Traditionally the pot ential form for covalent, ionic and metallic bonding have all been different from each other , but , as discussed in more detail in this dissertation, new advances in potential development have allowed these different interact ions to be described in one unified framework [ 20,55] . Periodic Boundary Conditions The initial atomic coordinates are given within a simulation box often called the supercell. The boundary conditions determine how this supercell is embedded in multidimensional space. Isolated boundary conditions (IBC) specify that atoms interact only with atoms listed within the original supe rcell. This type of boundary conditions is ideal for clusters or molecules. For liquids and solid bulk materials, periodic boundary conditions (PBC) are typically used in which the supercell is replicated in up to three dimensions to represent an infinite , periodic system. Atomic motions within the replicas

PAGE 38

38 are exactly the same as those in the original supercell and thus interactions are calculated for the original and replica supercells. To conserve the total number of atoms, when an atom moves outside the original supercell on one side it reenters from the opposite side. Problems arise when the supercell is too small and an atom interacts with itself in the adjacent replica image, which is unphysical. For example, atomic interactions described by th e Len nard Jones potential may be calculated out to infinite distances. At very large separations, the energy is essentially zero and their calculation only serves to increase the computational cost. Thus, a cut off value, r c , is specified in which atomic intera ctions longer than r c are ignored. Terminating interactions with a cut off leads to a discontinuity in the energy and forces at r ij = r c . Several functional forms have been developed to avoid these discontinuities by altering the energy function near the c ut off . Examples include the shifted force cut off and the Tersoff cut off [53,56] . The shifted force cut off modifies the energy function by a correction factor, while the Tersoff cut off uses trigonometric functions to ensure a smooth transition of the energy and forces. In order t o avoid having an atom interact with itself in the adjacent replica image, the size of the supercell must be at least twice as large as the cut off value r c . Time Integration Algorithms The heart of MD is the integration of atomic motion with time. The finite time difference b etween atomic positions is called the timestep, . The timestep must be very small in order to reduce accumulation of errors within the i nte gration algorithm . In particular, it is chosen so that dozens to hundreds of time points within a single atomic vibration are sampled. The vibration period for atoms is ~10 13 s, while c ommon timesteps in MD are on the order of femtoseconds (10 15 s) . Various integrators are

PAGE 39

39 used in MD such as the Verlet [57] , velocity Verlet, leap frog and others [53] . In this work, we use the Verlet and leap frog algorithms. The position of atoms are updated in one step within the Verlet algorithm by : (2 4) The leap frog algorithm is a two step process in which the half timestep velocities are generated first, followed by a full timestep update of the position: (2 5) (2 6) Although leap frog in faster than Verlet, it is n either numerically stable nor time reversible. The truncation errors are typically to . Thermodynamics Ensemble Conditions A thermodynamic state is characterized by state variables (also called state functions) that define the system. These state variables only depend on the value of a given property and not the process by which the value is determined [3] . Commonly used state variables include the number of atoms (N), volume (V), temperature (T), pressure energy (E) or entropy (S) . Their relationships may be derived from the laws of thermodynamics and values computed from experimental measurements such as the coefficient of thermal expansion or the heat capacity. The thermodynamic ensemble used in MD simulations depends on the state variables . Three im portant ensembles are microcanonical (NVE), canonical ( NVT) and [53] . The most

PAGE 40

40 common ensembles used in experiments are NPT (for solids and liq uids) and NVT (for gases), which will be the focus of this work. Thermostats Initialization of a n MD simulation usually involves assigning initial velocities to atoms based on the input temperature with directional velocities randomly chosen or based on a Maxwell Boltzmann distribution. For the random velocity method, a Maxwell Boltzmann distribution will be attained after a few timesteps. The instantaneous temperature at any time is readily calculated by the average kinetic energy via: (2 7 ) The total energy of the system at constant volume for a given time is (2 8 ) where V is the potential energy and E kinetic is the kinetic energy. As the simulation progresses, the total kinetic energy of the atoms in the system may change, thereby causing a change in temperature. These thermal fluctuations are caused by the interchange between kinetic and potential energy. To keep the desired temperature throughout the course of the simulation, a thermostat is applied to alter t he velocities of atoms. Many methods are available to control the temperature during the course of an MD simulation; these t hermostats include the Langevin [58 60] , Andersen [61] , Berendsen [62] and Nos é Hoover [63] algorithms . The thermostats used in this work are those of Berendsen and Nos é Hoover; the former will be used as an example in this section. The Berendsen thermostat rescales atomic velocities after each full step in the velocity V erlet algorithm by

PAGE 41

41 (2 9 ) where (2 10 ) is the rescaling factor. Here, is a constant used to control the temperature fluctuations; these fluctuations occur in the range of 0.5 2 ps. The target thermostat energy is (2 11 ) w here f is the total degrees of freedom in the system, k B is the Boltzmann constant and T ext is the ta rget external energy. The Berendsen thermostat does not conserve the total energy of the system , but does conserve total momentum. Barostat s In order to simulate isothermal isobaric (NPT) or isoenthalpic isobaric (NPH) conditions, it is imperative to allow fluctuations of the supercell size and /or shape. Barostats are applied to the supercell to maintain a desired pressure and/or stress, analogous to the thermostats controlling the temperature. In the Andersen scheme [61] , the supercell may be thought as a container of fluid compressed by pis tons in the three principle directions. The reduced coordinates of each particle is given by (2 12 ) where V is the volume. Thus, t he kinetic and potential energy may be rewritten as: (2 13 )

PAGE 42

42 (2 14 ) The isotropic expansion/contraction of the supercell is achieved via modification of the Lagrangian to : (2 15 ) Here, the first two terms on the right hand side (RHS) of eq. 2 15 represent the kinetic and potential energy of the particles in the supercell. The last two terms on the RHS of eq. 2 15 represent the kinetic and potential energy of the pistons where M and are constants . In order to allow anisotropic fluctuations to the supercell, more sophisticated approaches have been developed. Among these, the method of Parrinello and Rahman is widely used [64,65] . Shell Model Overview Materials with predominantly ionic bonding are well described by electrostatic interactions between the charged pa rticles expressed as (2 16 ) where U i is the potential en ergy of atom i, e is the charge and is the permittivity of free space . For a purely electrostatic interaction, the ions would collapse on themselves. This is the C oulombic catastrophe in which the electrostatic energy between two opposite charges tends towards negative infinity at small distances. To overcome this, a short range repulsion between oppositely charged ions must be added to the interaction energy. Commo n functional forms of repulsion are Buckingham and Morse terms, where the energy increases exponentially as ions approach one another [66,67] .

PAGE 43

43 Physically this represents the overlapping electron clouds with electrons being forbidden to have the same quantum state by the Pauli exclusion principle. The shell model [68] is a widely used method for describing polarizability in ionic systems, including ferroelectric materials. Electronic polarizability arises from the displacement of valence electrons within an atom, wh ile ionic polarizability occurs for displacements between ions . In the shell model, each ion is described as a positively charged core connected to a negatively charged shell via a spring. Physically this represents the displacement of the cloud of valence electrons with respect to the nucleus and core electrons . This displacement may occur due to the presence of an electric field or caused by random fluctuations during MD . Long range electrostatic interactions occur among cores and shells of different atom s, but not within the same atom. Short range repulsions are ascribed only between the shells of different atoms. Three body terms may also be included to improve the original shell model. In this work, the three body interactions are placed on the shells a nd limited to first nearest neighbors . The shell model has two implementations in MD codes : one in which the shells have mass (adiabatic shells) [69] and one where they do not (relaxed shells) [70] . In the The fictional mass assigned to the shell should account for the natural vibrational frequency of the atom. In this work, shells were assigned 10% of the total atomic mass. For the relaxed shells, the shell motion responds to that of the cores, with position such that the forces acting on the shell are zero [71] . This relaxation of the shell can be very time consuming ; thus the adiabatic shell model is preferred in this work.

PAGE 44

44 Shell Model Formalism In this work, the total potential form of the shell model i s : (2 17 ) The first four terms on the right hand side of eq. 2 17 are interatomic interactions, while U Core Shell is an intra atomic interaction. Interatomic Interactions The electrostatic energy is a simple r 1 contribution where : (2 18 ) Electrostatic interactions are calculated by the Smooth Particle Mesh Ewald (SPME) method [72] , which is a modification of the standard Ewald sum [53] . The Ewald sum makes two m odifications to the coulombic interactions of charged ions, which is an infinite sum in real space. The first is a finite sum in real space where a spherical Gaussian of opposite charge is centered on each ion. The second is a finite sum in reciprocal spac e over Gaussian charges with the same charge as the original ion. These two finite sums, in addition to a self energy correction term, allow the electrostatic energy to be calculated to a converged value. For the systems considered in this dissertation, t h e short range repulsion s have Buckingham [66] and Morse [67] terms : (2 19 ) (2 2 0 )

PAGE 45

45 For the Buckingham potential, the parameters A, and C are relate d to the ionic stiffness, radii and strength of the van der Waals force for the interaction between two ions , respectively. For the Morse term, E 0 is the well depth at equilibrium b ond length r 0 . Physically, the Buckingham potential represents the repulsion between electron clouds, while the Morse potential represents the anharmonic vibration of covalent bonds [66,67] . Materials with complete ly ionic bonding oftentimes have integer formal charges on each ion , although covalent systems may use this assumption as well . However, bonding in fe rroelectric materials has a substantial covalent component , which means that the use of non integer partial charges is appropriate . These partial charges, in addition to the Morse term, account for the covalent bonding in these materials. The use of Buckingham and Morse terms in conjunction allows the complex bonding in ferroelectric materials to be described better than the Buckingham potential alone. The angular terms are a superposition of unscreened and screened harmonic potentials [71] : (2 2 1 ) Here, the first three terms on the right hand side of eq. 2 2 1 only depend on the difference between the angle between three atoms, ijk , and a set value, 0 . The last term is scaled by an exponential factor that depends on the difference in bond lengths between r ij and r ik . As we shall see, t he screened harmonic term is useful in describing Ti O bonds in BaTiO 3 , where the ferroelectric phases each have differing Ti O bond distances.

PAGE 46

46 Intra atomic Interaction An intra electron cloud due to an electric field. In the shell model, the shells are connected to their cores via a spring (2 2 2 ) where w is the distance between the core and shell for the same atom [73] . The spring constant, k 2 , and shell charge, Y, are related to the electronic polarizability by : (2 2 3 ) Overview of the Charge Optimized Many Body (COMB) Potential The foundations of COMB are built upon two important concepts: bond order and self consistent charge equilibration. Abell [74] used chemical pseudopotential theory to show that the strength of a bond between atoms varies by the inverse square root of the number of neighbors such that (2 2 4 ) where Z is the total number of first nearest neighbors. Tersoff expanded this formalism for use in the silicon system . simulated using embedded atom method (EAM) potentials, while semiconductors (covalent bonding) were simulated using Abell Tersoff potentials. However, Brenner showed that these two formalisms are functionally equivalent [75] . The major drawback of traditional potentials is that electrostatic interactions were co nstrained to fixed point charges. This makes it impossible to simulat e heterogeneous interfaces (such as the metal /oxide interface) with their complex and changing bond

PAGE 47

47 environment at the interface. To overcome this, several methods were developed which al lowed dynamical charge transfer between atoms based on their local environment. Rappe and Goddard developed the charge equilibration (Qeq) approach [76] , which used ionization energ ies and electron affinities as inputs to calc ulate the chemical potential between two atom species. In this approach, a tomic c harges change so that the chemical potential becomes equal at all atomic sites . For charge equilibration, COMB uses the electronegativity equilibration (EE) principle based on the Qeq approach of Streitz and Mintmire [77] . Yasukawa extended the Abell Tersoff potential to include charges within the Qeq approach for the Si/SiO 2 system [78] . Yu et al. added correction terms to the Yasukawa potential, which became the first generation COMB potential [20,79] . Devine et al. [80] extended the COMB formalism by replacing the point charge mo del with Coulomb integrals over Slater 1s orbitals [81] . This formulation became the second generation COMB potential with many systems being parame terized for it [21,80,82] . The COMB potential is now in its third generation [22] . The newest inclusions are additions from the REBO potentials [83] , which are widely used to simulate hydrocarbons. This formalism also correctly describes molecular oxygen as the ground state, which the second generation COMB could not do. COMB Potential Formalism The terms of the COMB formalism that are relevant for t his work are presented in this section. A metal/oxide interface will be used to illustrate how the relevant COMB terms affect this system. The total energy of COMB is dependent on a series of elec trostatic and short range terms

PAGE 48

48 (2 25 ) where q is the charge, r is the distance between atoms, and ijk is the angle between a triad of atoms. The electrostatic energy portion is (2 26 ) where U Self is the energy required to form a charge on an isolated atom, U qq is the electron electron interaction, and U qZ is the electron nuclear interaction term. For non charged systems, such as pure metals, the electrostatic energy is set to zero. The self energy term , which is species dependent, is denoted as a fourth order expansion of the energy as a function of charge (2 27 ) where and J repre se nt the Mulliken electronegativity and chemical hardness, respectively. A correction term, E Field , changes the charge state of an atom when it is embedded in a molecule or ionic lattice : (2 28 ) Here, parameters P and P J are bond type specific and help modify the bulk properties. The electron electron and electron nuclear terms are calculated by integrating over the electron density of the system. The electron density of an atom, i , is based on S type Slater orbitals [81] : (2 29)

PAGE 49

49 (2 30) wh ere, Z i is the core charge, is the Dirac delta function and is an orbital exponent that controls the radial decay of the electrostatic energy. The C oulombic repulsion is: (2 31) (2 32) The electron nuclear coul o mbic attraction is an interaction between the electrons and core charges such that: (2 33) (2 34) (2 35) In COMB, there are no nuclear nuclear coulombic repulsion terms. This repulsion is instead modeled by an exponentially decaying term, similar to the Buckingham term in the shell model. However, the short range interactions are modified by the local charge distribution in COMB, whereas charges are fixed in the shell mode l. The short range interaction is given by (2 36 ) w here (2 37 ) is the short range attraction term and

PAGE 50

50 (2 38 ) is the short range repulsion. The change dependence of the short range contributions change according to (2 39) (2 40) (2 41) w here D U and D L are parameters that reflect the change in atomic radius with charge. Likewise, Q U and Q L are the atomic charges that correspond to the limits of the valence shell. When the charges are zero, the potential is functionally equivalent to the Te rsoff potential. Equations 2 25 through 2 41 illustrate the variable charge distributions found in real systems. In the metal/oxide interface example, the oxide far from the interface has mostly electrostatic energy contributions. Covalent and metallic bo nding between nearest neighbors are taken into account through the short r ange attraction term in eq. 2 36 . In the metal far from the interface, the electrostatic energy is essentially zero and all interac tions are determined by eq. 2 36 . Atoms close to th e metal/oxide interface have mixed bonding with the charge distribution on atoms ranging from zero to +/ q max . T he bond order term, b ij , in eq. 2 36 modifies the short range attraction based on the local environment of the bond between atoms i and j

PAGE 51

51 (2 4 2 ) which is composed of terms based on bond angles and local coordination. The screening function is (2 4 3 ) which returns a value different from unity if r ij differs from r ik . This allows longer bonds to be weakened, thereby preferring shorter bonds. The angular function (2 4 4 ) is a sixth order polynomial dependent on the angle between a triad of atoms formed by the atom pairs ij and i k. The coordination term is given by: (2 45 ) Here i is the total number of neighbors of the central atom i, (2 46 ) where t he superscript represents the element identification number and N i l represents the total number of l th element neighbors of atom i . The function F c ( r ij ) is the Tersoff cut off function :

PAGE 52

52 (2 47 ) Th is cut off allows the bond order to transition fro m 1 for bond distances less than R min to 0 for bonds greater than R max . The last term in eq. 2 25 is the Legendre polynomials (LP) : (2 48 ) These Legendre polynomials are a series of functions that allows correct ion to the energies associated with bond angles. In practice, it is best to avoid these terms if possible as these often have multiple minima that may affect dynamic results. The LP contribution to the total energy should be very small. Parameterization of Inter atomic Potentials Challenges of O ptimization The use of interatomic potentials in MD simulations requires a determination of parameters based on its many functional forms . While there are times when a parameter may be directly related to an experimentally observable quantity, oftentimes an intuitive relationship between a parameter and specific physical quantity is missing. Simple analytic forms, such as Lennard Jones or Buckingham , have only two or three parameters to be determined. As the complexity of th e formalism increases, such as that for COMB, the number of parameters that must be determined rapidly increases to upwards of 50 . In order to navigate through this extremely complicated parameter space, it is often best to fix certain parameters so that o ptimization is constrained to a

PAGE 53

53 smaller space. itting is best done in a step by step process, in which parameters related to specific interactions are fitted and then fixed as the fitting procedure progresses. Inter atomic potential fitting is complex because there are often open ended bounds on the value of different parameters. For a set of any given parameters, the local (global) minima in the multi dimensional parameter space may be highly dep endent on the trainin g database and this database may require constant alteration throughout the fitting process. In nearly all cases, the global minima in unknown and may be inaccessible due to the optimization method or time constraints. Training Database and Cost Function The training database , or simply database, is a collection of structure files that contain information such atomic species type and their spatial positions. Material properties such as cohesive energy, elastic constants, surface energy etc. are supplied an d these become the target values for parameter optimization. These properties may be generated from experiment s or DFT calculations , but care must be taken to ensure that the database is self consistent, as contradictory information will lead to a poor pot ential. For example , the use of lattice parameters and internal coordinates from experiment , but elastic constants from DFT calculations may lead to an ill defined curvature in the energy function around the equilibrium atomic positions. Generating the ent ire database from DFT calculations alone is self consistent, but may not be desirable if the calculations significantly deviate from known experimental properties. Consider the construction of the database for an arbitrary unary system where the ground state is a trimer between atoms i j k as shown in Figure 2 1. To ensure that the ground state predicted by the potential is indeed a trimer, structures of varying

PAGE 54

54 coordination and spatial configurations should be added to the database, such as dimers and c ommon crystal structures like face centered cubic (FCC) and simple cubic (SC). For a given parameter set in an arbitrary potential formalism, the energy and other properties for each st ructure are calculated and recorded. If the desired properties are not achieved, such as the trimer not being the ground state, then the parameters will need to be altered. The difference between the calculated and target values for various properties may be summed into a cost function [84] : (2 49) (2 50) Here, the weight is a meas ure of the precision between target and calculated values and has units that ensure the total cost is dimensionless. A total cost of zero implies that all the properties in the database are perfectly described by the potential. In practice, the weights for different properties should be distributed such that importance is given to those that contribute most to the desired outcome of the potential. To help guide the parameter optimization towards the desired properties of the ground state, it is useful to supply structures that have slight deviations from the target structure. For the trimer example, the database should include a series of trimer structures that have different bond lengths and angles. If the target bond lengths are r ij and r ik (Figure 2 1) , then the database should include trimers with bond lengths r ij +/ r i k +/ and combinations thereof 0.1% of the target value. The same process may be used to accurately fit the bond

PAGE 55

55 angle DFT calculations are very useful in this regard as the energy differences between arbitrary trimers are readily obtained. However, addin g too many structures with slight deviations may over constrain the energy landscape in a potential formalism with limited flexibility. No classical potential is capable of accurately describing all the properties of a given material. Thus, the desired pro perties described by the potential should be identified at the beginning of the fitting process. The relevant structures needed to achieve this outcome are then added to the database. In this work, the database for COMB Zr O H included many Zr surface stru ctures with oxygen/hydrogen occupying different surface or subsurface sites. Meanwhile, the database for ferroelectric materials focused on the bulk structures and the relative energy differences between phases. Thus, there is no single procedure for creat ing a database; however, the number of properties included should be expansive enough to ensure that the final potential has no unreasonable predictions for alternate atomic environments. Optimization Strategies The simplex algorithm is downhill of optimization where only moves that minimize a cost function are accepted [85] . The simplex method is very numerically robust and highly efficient, but may have difficulty escaping shallow local minima. One major problem with the simplex method is the uniform step size, the percent age by which each parameter is changed, which may be unable to greatly change parameters tha t only slightly alter the cost function. Simulated annealing is another optimization technique, one that allows both uphill and downhill moves [86] . Simulated anneali ng conceptually is similar to the thermodynamic cooling of molten metal. As the temperature of a molten metal is cooled,

PAGE 56

56 atoms diffuse to find the global minima. If this process is slow enough, the global minima will be found upon cooling. However, if the temperature is rapidly dropped ( i.e., the system is quenched), then the system may remain in some higher energy local minimum . In terms of the number of calls to the cost function, simulated annealing requires many more than the simplex algorithm. However, its advantages include the ability to escap e local minima, which makes it a competitive approach. Software Used i n t his Dissertation Parameterization of the COMB Zr O H system used our in house code Potential Optimization Software for Materials (POSMat) [84] , the basic version of which was written by Dr. Dundar Yilmaz. The author added and modified several routines to better serve the purposes of the work described in this dissertation . For parameterization of ferroelectric oxides, POSMat and Gener al Utility Lattice Program ( GULP ) [87] were used in conjunction. Molecular dynamics simulations used the Large scale Atomic/Molecular Massively Par allel Simulator ( LAMMPS ) [88] and DL _ POLY [89] software packag es. First principles calculations were performed with the Vienna ab initio software package (VASP) [90] .

PAGE 57

57 Figure 2 1. Schematic of bond distances and angle between three atoms.

PAGE 58

58 CHAPTER 3 CHARGE OPTIMIZED MANY BODY ( COMB ) POTENTIAL S FOR ZIRCONIUM Overview of Slip in HCP Metals The hexagonal close packed (HCP) structure presents unique challenges in the description of interatomic interactions compared to face centered cubic (FCC) or body centered cubic (BCC) materials. While FCC and BCC only have one lattice parameter (a 0 ) that determines interatomic distance and planar spacing, HCP has two (a 0 and c 0 ). A c 0 /a 0 ratio (~ 1.633 ) corresponds to a n ideal close pack ing of hard spheres. The sequence of hard sphere packing between HCP and FCC differs in that HCP stacking may be represented by ABABAB repeating planes , while FCC has ABCABC repeating places . Deviations from the ideal c 0 /a 0 ratio in HCP metals change the a ngular distributions and interplanar distances. HCP metals with ideal or greater c 0 /a 0 ratios have basal {0001} planes as the most closely packed, with interplanar spacing of c 0 /2. When c 0 /a 0 is less than ideal, prism planes are more closely packed. Th is difference i n planar density play s a major role in the resulting mechanical behavior of HCP metals [91] . A dislocation is a one dimensional crystal defect where atoms do not lie in their regular crystallographic positions. The process of slip is the gliding motion of dislocations and occurs in specific crystallographic p lanes and directions based on the crystal structure. The combination of a slip plane and slip direction is called a slip system. Figure 3 1 shows the most common slip systems for HCP metals. The primary Mu ch of the work is this chapter has been published in M. J. Noordhoek, T. Liang, Z. Lu, T. R. Shan, S. B. Sinnott, and S. R. Phillpot, Journal of Nuclear Materials 441 , 274 (2013). The simulations on the mechanical properties of Zr discussed in the above pa per and in this chapter were carried out by Z. Lu.

PAGE 59

59 slip direction is in all HCP metals [92] . However, HCP metals with c 0 /a 0 equal to or greater than the ideal ratio undergo slip on the basal planes, while those with c 0 /a 0 less than the ideal value generally undergo slip on the prism planes. Exceptions to this trend in clude Mg, Co and Be where the c 0 /a 0 ratio is less than ideal, but slip occurs on the basal planes [91] . This contrasts to metals like Ti, Zr and Hf, where the partially filled d band electrons are thought to create directional bonding with covalent characte r [92] . Thus, the d band electrons affect the slip behavior in HCP metals. Also, t he orientation relationship between dislocations and crystal structure play s a major role in the resulting slip behavior of HCP materials. In particular, a dislocation may dissociate into partials in the presence of metastable stacking faults [92] . Thus, the difference in slip behavior for e energy landscape associated with the dissociation of dislocations into partials. The principle and secondary slip systems for Zr are listed in Table 3 1. A stacking fault is an irregularity in the stacking sequence between planes such as ABABCACA in HCP . To visualize a stacking fault, consider a crystal cleaved into two halves along a given crystallographic plane. The energy of the system will increase w hen the halves are displaced perpendicular to the cleaved direction . These plots of energy versus atomi c displacement are commonly called surfaces or stacking fault maps [92] . They play an important role in the development of interatomic potentials by providing information about the energy landscape . For single crystal materials, atoms are in their regular lattice positions and fill three dimensional space. A polycrystalline material contains many crystals with different size s and orientation s . The regions where two or more crystallites ( i.e., grains) meet are

PAGE 60

60 grain boundaries. Grain boundaries have different mechanical properti es than the perfect structure and can alter the strength of a material by impeding the motion of dislocations. The relationship between yield strength, , and grain diameter, d, is given by the Hall Petch equation [93] . Thus , smaller grain sizes generally have increase d yield strength. Howeve r, at very small grain sizes ( Petch relationship is observed where decreasing the grain size further will decrease the yield strength [93,94] . This occurs because grain boundary sliding becomes the dominant deformation mechanism for nanocrystalline materials. The goal of this work is to develop a COMB potential that is suitable for the simulation of nanocrystalline Zr under mechanical deformation . Emphasis during the fitting procedure is placed on correctly describing the c 0 /a 0 ratio as well as the energy landscape associated with slip processes. In addition, the potential must be acceptable for extension to the Zr O H system. In the following sections, parameterization and results of two different COMB potentials for Zr are given. Pa rameterization of COMB Potential s f or Zr This section contains a summary of t he COMB formalism related to metal fitting. A full description of COMB may be found in Chapter 2. For unary metal systems, the charge transfer between atoms may be safely neglect ed in MD. Thus , i n the description of Zr , atomic charges are set to zero, which leads to the formalism for the short range interactions to be equivalent to the Tersoff Brenner potential (equations 2 36 to 2 48 in Chapter 2). The short range interaction, U S hort ( r ij ), is summarized as (3 1) W here

PAGE 61

61 (3 2) and (3 3) The bond order term, b ij , is defined as (3 4) with (3 5) and (3 6) The cut off function is (3 7) Parameterization of the Zr Zr interactions is accomplished with our in house code Potential Optimization Software for Materials (POSMat) [84] . Target values are obtained from experiments and DFT calculations. The conditions for the DFT calculations in this work are the same as those used in Ref. [95] as described here. Vanderbilt U ltrasoft (US) pseudopotentials within the generalized gr adient approximation (GGA) were chosen [90,95 97] . The first Brillouin zone is sampled with a 3 3 4 grid; the energy cut -

PAGE 62

62 off is 225 eV. The calculated lattice parameters for this configuration are a 0 = 3.23 Ã… and c 0 = 5.18 Ã… , which results in c 0 /a 0 = 1.60 . The simplex [85] optimization technique is used to minimize a cost function, which is determined by differences between the target and calculated values. This cost function is calculated from weighted square absolute deviations of lattic e parameters, elastic constants and other structural properties. The Zr Zr parameters are fit to a database (13 structures total) consisting of bulk structures (5), stack ing faults (6 ) and surfaces (2). A full list of the properties and their target values is given in Table 3 2 . This dissertation contains two distinct parameter sets for Zr metal. These sets are dubbed COMB Zr 1 [98] and COMB Zr 2 ; the parameter values may be found in Appendix A . Although the fitting database use s data from experiment, DFT calculations from the literature and DFT cal culations from this work, proper care was taken to avoid potentially fatal inconsistencies in the database. This involve s placing the weights on the cost function such that all the necessary trends c an be obtained (e.g. phase order) without prohibitively d amaging the value of other properties. Emphasis is placed on the stacking fault energy trends, which are most pertinent to the mechanical behavior Zr . Point defects are excluded from the database and thus constitute predicted values. The resulting values are acceptable for the final use of these potentials. COMB Static Property Results Bulk and Surface Properties Table 3 2 compares the properties of zirconium determined from experiment, DFT calculations, COMB and an embedded atom method ( EAM ) potential . Th e EAM formalism contains two body interactions and have been shown to accurately describe many metallic systems [17] . The selected EAM potential is used as a direct comparison

PAGE 63

63 for the mechanical deformation simulations [99] . The experimental a axis lattice parameter is 3.23 0 Ã… at 4 K [100] ; values for COMB , at 0 K, are 3.226 Ã… for COMB Zr 1 and 3.229 Ã… for COMB Zr 2. The experimental c axis lattice parameter is 5.15 0 Ã…; COMB yields values of 5.159 Ã… for COMB Zr 1 and 5.155 Ã… for COMB Zr 2. The c 0 /a 0 ratio of 1.599 (COMB Zr 1) and 1.596 (COMB Zr 2) compares favorably with the experimental value of 1.594. These values are both much less than the ideal HCP c 0 /a 0 ratio of 1.633, which means prism planes are more closely spaced than basal {0001} planes. Thus, slip is more likely to occur on the prism plane upon mechanical deformation, which is observed experimentally [91] . Zirconium surface energies are an important property to consider when studying interfacial or corrosion behavior. The wide rang e of Zr alloy compositions and surface preparation techniques have a d ded to the difficulties in understanding Zr surface reactivity with corrosive species such as water, O 2 and H 2 [101] . Insights obtained from m odeling these systems with COMB or other methods may ultimately allow more precise control over these properties. The average experim ental surface energy is 2050 mJ/m 2 [102] , while DFT calculations are about 80% of that : 1600 mJ/m 2 and 1660 mJ/m 2 for the (0001) and surfaces, respectively [103] . Surface energies in COMB match DFT within 11%, but only COMB Zr 2 accurately predicts the basal (0001) surface energy to be lower than prism (Table 3 2 ) . Point Defects It is known that irradiation induced growth and creep plays an important role in material behavior when Zr alloys are used as nuclear fuel cladding applications [104] . There are eight interstitial configurati ons commonly referred to in the literature, as

PAGE 64

64 shown in Figure 3 2 [105,106] . The octahedral (O) and tetrahedral (T) sites lie between basal planes and diff er in their placement along the direction. Their basal counterparts (BO and BT) lie directly within the (0001) plane. The crowdion (C) interstitial is located between atoms in adjacent basal planes, while the basal crowdion (BC) lies between neighboring atoms in the <1000> direction within the basal plane. Split (S) and basal split (BS) interstitials correspond to two atoms sharing one lattice position, with atoms being displace d along the <0001> and <1000> directions, respectively. In total, four interstitial sites lie within the basal plane while four are located between basal planes. Results of first principles calculations showed a strong dependence of Zr interstitial energ i es on the conditions of the details of the computational approach [95,105,106] . In particular, early studies were limited to very small supercells, consisting of 4 4 3 repeat units containing 96 + 1 = 97 atoms [95,106] . In one study [106] , the calculations were performed with constant volume and shape of the cell, allowing the internal relaxation of atoms. In others [95,105] , the calculations were performed at constant volume per atom, while maintaining the shape of the supercell. The defect energies converge slowly with increasi ng system size for each calculation method [14 ]. In particular, it was found that the supe rcell should be at least 8 unit cells along <1000> for accurate results. Increasing the supercell size along <000 1 > exhibited much smaller relative differences than < 1 000> in the defect formation energy. The directional dependence of the supercell size on the defect energy is due long range strain effects , because bond distances with the interstitial within the (0001) plane are generally shorter than those out of plane .

PAGE 65

65 Since no point defect information was included in the fitting database, all of the COMB values presented here are predictions of the potential. The COMB calculations for interstitial formation energies were performed using constant supercell volume while allowing for internal relaxation. Two supercell sizes of 10 4 3 and 30 20 20 unit cells were used , the latter being inaccessible in DFT due to its large size . DFT calculations [105] suggest that interstitial sites within the basal plane have lower formation energies than their counterparts out of the basal plane. Calculations using large supercells showe d the basal octahedral (BO) position to be the most s table , followed closely by the basal split (BS) [10 5] . COMB Zr 1 predict s the BO position as the most stable site, but overall the interstitial energies are about 0.8 eV larger than DFT (Table 3 3 ) . COMB Zr 2 predictions for point defect energies are markedly worse than COMB Zr 1, with the relative ener gies between sites exhibiting the incorrect trends , as compared to DFT. Thus, COMB Zr 1 is preferred for studies related to irradiation damage in Zr. The full results are sho wn in Table 3 3 . Values l 3 indicate that the inters titial spontaneously moves to another interstitial site or forms a complex structure in which the local symmetry is lost. As a n example of the supercell size dependence on the defect formation energy, the BO position is chose n to compar e DFT and COMB Zr 1 . Using DFT, supercell sizes of 6 4 2 and 10 4 3 give values of 2.87 eV and 2.85 eV, respectively. For COMB Zr 1 , the values for 6 4 2 and 10 4 3 supercells are 3.74 eV and 3.62 eV, respectively. Thus, there is a much larger (6X) energy decrease in COMB Zr 1 , which may be attributed to the abrupt short range cut off function. As noted in Chapter 2, the maximum separation between interacting particles is given by the cut off distance, r c . In

PAGE 66

66 COMB Zr 1, the value is r c = 4.2 Ã…, which is much shorter than the strain field distance associated with an interstitial. In addition, the small supercell size of 6 4 2 repeat units leads to an unphysical strain interaction across the period boundary conditions. At a supercell size of 30 20 20, the BO energy in COMB converge s to a value of 3.53 eV. This signifies that the strain field around the interstitial positions is long range d . We have performed DFT calculations for vacancy and divacancy formation energies using VASP ( Chapter 2 ) . He re we used a larger supercell consisting of 5 6 3 repeat unit cells containing 180 atoms. The calculated lattice parameters for this configuration are 3.23 Ã… in the a direction and 5.18 Ã… in the c direction. Point defect formation energies were calculated using constant cell shape and volume, while internal relaxation of the atoms was allowed. Keeping the volume constant does not allow the strain from the point defect to alter the lattice parameters, which may occur for the relatively small supercells consi dered in DFT. This procedure resulted in a vacancy formation energy of 1.94 eV , which is similar to 1.86 eV obtained from a previous DFT study [95] . The vacancy fo rmation energies in COMB are 1.8 6 eV (COMB Zr 1) and 2.01 eV (COMB Zr 2) as shown in Table 3 2 . A schematic showing the only two possible structures of adjacent vacancies in an HCP lattice is given in Figure 3 3 . A distinction is made between vacancies next to each other within the basal (0001) plane (in plane) and those in adjacent basal planes (out of plane). The DFT results show that divacancies within the same basal plane have positive binding energies, while those in adjacent basal planes have negative binding energies. A negative binding energies means it is more energetically favorable for the vacancies to be well separated. For COMB and other empirical potentials, the

PAGE 67

67 divacancy binding energies are all positive and no significant differences for in plane or out of plane divacancies are observed. Thus, predictions for high temperature behavior where divacancies are likely to be generated, such as the melting process, will not likely be quantitatively reliable. The values predicted by COMB are larger than DFT, while LREP and AMEAM underestimate the DFT results. The se r esults are summarized in Tab le 3 4 . Stacking Fault Energy The stacking fault energy map for the (0 001) basal plane from COMB, DFT and EAM [107] is shown in Figure 3 4 . Only the basal I 2 stacking fault is included in the fitting database ; the rest of the energetics are predictions. The energies of the various translational configurations of the stacking faults were calculated by only allowing atomic relaxation in the direction perpendicular to the stacking fault plane. The indirect path , which results in the I 2 stacking fault, is leading to saddle points in the center and corners of the maps in Figure 3 4 . The values of the stable and unstable stacking faults on the basal plane from COMB, DFT and EAM are listed in Table 3 5 . The values for COMB are acceptable because the basal stacking fault energy is larger than those for the prism planes. In addition, the basal unstable stacking fault energy is much larger than the basal stable stacking fault , which prohibits (0001) basal slip from being the dominant slip system . The stacki ng fault energy map for the prism plane from COMB, DFT and EAM [107] is shown in Figure 3 5 . Here, the direct path is the prismatic dislocation. The stable and unstable stacking faults are included in the fitting database. The values of the stable and unstab le stacking faults on the prism plane from COMB,

PAGE 68

68 DFT, and EAM are listed in Table 3 5 . As noted previously, b oth the stable and unstable stacking fault energies in the prism plane are lower than those of the basal plane. This is crucial ly important as the energy difference between the prism stable and unstable stacking faults will influence the amount of stress required to cause slip. T he difference between the stable and unstable stacking fault for the prismatic dislocation from COMB Zr 1 (32 mJ/m 2 ) i s very similar to the DFT value (28 mJ/m 2 ) , while COMB Zr 2 (54 mJ/m 2 ) slightly overestimates this value . The Legrand [108] criteria for the type of slip observed is based on , where R<1 leads to basal slip, R>1 leads to prism slip, and R 1 lea ds to both types being activated. COMB Zr 1 has a value of R=1.61 while COMB Zr 2 is 1.79. This compare s to 1.44 for DFT [107] , 1.20 for an EAM potential [99] and 2.3 from Legrand who used a pseudopotential theory . Since the stacking fault energetics largely determine the mechanical response of a metal, the above analysis gives us confidence th at the COMB potential s for Zr will correctly reproduce at least the qualitative features of the mechanical response of Zr. As a check on this assertion, we show preliminary results for slip behavior in nanocrystalline Zr in the next section . Mechanical De formation of Zr In order to test the appropriateness of this potential for the study of the mechanisms of mechanical deformation, in this section we show preliminary results for a tensile test of a 3 D zirconium nanocrystal . These simulations were performe d by Z . Lu and published in a joint paper [98] . They are included here as they provide an important assessment of the mechanical performance of the COMB Zr potential. The simulations are performed in the LAMMPS software (Chapter 2) [88] using the COMB -

PAGE 69

69 Zr 1 potential. The nanocrystalline structure as shown in Figure 3 6 contains 16 randomly oriented grains with average grain size of 17 nm, consisting of total of approximately 1.8 million atoms. This structure is annealed at various temperatures up to 900K and then cooled to room temperature. This annealing process allows atoms within the grain boundar ies to diffuse to equilibrium positions. In order to distinguish FCC like atoms from normal HCP atoms, common neighbor analysis (CNA) i s used [109] . The structure i s visualized with AtomEye [110] . A tensile test is conducted at 300K with a strain rate of 4.0x10 9 s 1 using a timestep of 0.1 fs. An earlier study of Mg used a strain rate of 1.5x10 9 s 1 [111] , while a yet earlier study of FCC Al used a strain rate of 2.0x10 8 s 1 [112] . The high strain rates used in MD as compared to experiments ( typically less than ~ 10 3 s 1 ) are due to the very short timescales accessible by MD. Nevertheless , these high strain rates have been shown to describe the critical relevant physics in nanocrystalline metals [112] . First , even at these strain rates, the velocity of glide dislocations is less than the speed of sound in the material ; thus the system is not in a shock regime [94] . Second , Coble creep plays a major role in the deformation mechanism for small grain sizes . Coble creep is caused by atomic diffusion along the grain boundaries and has a strain rate dependence on the grain diameter given by [113] . Thus, the relative strain rate of a polycrystal with grain sizes 10 0 10 nm is 10 12 , which is on the same order as the ratio of experimental to MD strain rates . The common slip systems found in Zr are given in Table 3 1 and illustrated in Figure 3 1 . Prismatic slip and the pyramidal slip are frequen tly observed in experiments [92] , with prismatic slip being most easily activated.

PAGE 70

70 During the simulated tensile test, dislo cation activity was observed in 15 of th e 16 grains as shown in Figure 3 6 . Both prismatic and pyramidal dislocations were also observed. The prismatic is a full dislocation, while the pyramidal dislocation is a partial . T hese simulation results are co nsistent with the above experimental observations . T hus , it appears that the COMB potential can successfully describe the dislocation behaviors in zirconium . A full analysis of these dislocations is given by Lu et al. [114] . Assessment of COMB Zr P otentials Overall, the results for the COMB Zr 1 and COMB Zr 2 parameterizations are very similar. Both accurately describe the lattice parameters, elastic properties and stacking fault properties of metallic Zr. Most of these properties are within 20% of the target value, which is generally considered acceptable within the computation al field. The level of agreement for COMB as compared to other empirical potentials for HCP, FCC and BCC metals in the literature is quite satisfactory in terms of static properties as compared to experiments . In addition, the ability of COMB to simulate t he mechanical deformation of nanocrystalline Zr proves that the formalism is suitable for use in many other metallic systems. The major difference between COMB Zr 1 and COMB Zr 2 is for interstitial defect energies, which were not emphasized in this work. Their predictions for irradiation induced damage likely lie outside of experimentally observed phenomena; however, no empirical potential is capable of correctly describing all the properties of a given material.

PAGE 71

71 As part of continuing work, the COMB Zr 2 potential is extended to the Zr O H system and is the subject of the next chapter. COMB Zr 2 is chosen over COMB Zr 1 because the bond order scaling factor, in eq. 3 4, is unity and therefore does not adversely affect the Zr O and Zr H fitting process. T he resulting Zr O H potential makes significant progress towards understanding of complex phenomena in this system and shows how other materials may benefit from an atomistic description using COMB.

PAGE 72

72 Table 3 1 . Slip systems in zirconium [91] . Principle Secondary

PAGE 73

73 Table 3 2 . Properties of zirconium metal given by the COMB potentials compared with experiments, DFT and an EAM potential . A * indicates that the select property value was included in the fitting database. Properties Experiment [26,100,102] DFT [103] EAM [99] COMB Zr 1 [98] COMB Zr 2 [115] a 0 ( Ã…) 3.23 0 * 3.23 0 3.234 3.226 3.229 c 0 ( Ã…) 5.15 0 * 5.17 0 5.168 5.159 5.155 c 0 /a 0 1.594 1.60 1 1.598 1.599 1.596 E c (eV/atom) 6.32 * N/A 6.6 4 6.32 6.31 C 11 (GPa) 155 * 156 147 154 154 C 12 (GPa) 67 * 61 69 71 74 C 13 (GPa) 65 * 62 74 62 59 C 33 (GPa) 172 * 166 168 164 173 C 44 (GPa) 36 * 26 44 32 31 B (GPa) 97 94 100 96 96 G (GPa) 42 39 42 39 40 hcp) (eV/atom) N/A 0.040 * 0.054 0.074 0.073 hcp) (eV/atom)* N/A 0.078 * 0.103 0.105 0.152 Basal I 2 stacking fault (mJ/m 2 ) N/A 227 * 198 267 267 Prism stacking fault (mJ/m 2 ) N/A 197 * 175 215 193 Basal surface (mJ/m 2 ) 2050 1600 * 1770 1770 1550 Prism surface (mJ/m 2 ) 2050 1660 * 1930 1740 1553 E Vacancy (eV) N/A 1.86 2.31 1.86 2.01

PAGE 74

74 Table 3 3 . Zr interstitial formation energies in DFT and COMB . Units are in eV. DFT (GGA) [95] DFT (GGA) [105] COMB Zr 1 [98] COMB Zr 2 [115] Supercell dimensions 4 4 3 10 4 3 30 20 20 30 20 20 Octahedral (O) 2.84 2.93 4.12 3.01 Basal Octahedral (BO) 2.88 2.75 3.53 4.16 Crowdion (C) 3.08 N/A Unstable 3.21 Basal Crowdion (BC) 2.95 N/A Unstable 3.81 Split (S) 3.01 3.02 4.57 3.79 Basal Split (BS) N/A 2.87 3.60 4.25 Tetrahedral (T) Unstable N/A 3.78 4.15 Basal Tetrahedral (BT) 4.03 N/A 4.14 3.95

PAGE 75

75 Table 3 4 . Divacancy formation and binding energies in DFT, COMB and other empirical potentials . Units are in eV. DFT (This Work) LREP [116] AMEAM [117] COMB Zr 1 [98] COMB Zr 2 [115] in plane E Formation In 2v 3.78 3.30 3.18 4.91 4.03 E binding In 2v 0.07 0.20 0.22 0.12 0.15 out of plane E Formation Out 2v 3.95 3.31 3.18 4.91 4.03 E binding Out 2v 0.10 0.19 0.21 0.12 0.15

PAGE 76

76 Table 3 5 . Stacking fault energies of zirconium metal given by the COMB potentials compared with DFT and an EAM potential. Properties DFT [103] EAM [99] COMB Zr 1 [98] C OMB Zr 2 [115] Stable Basal stacking fault (mJ/m 2 ) 213 199 267 267 Unstable Basal stacking fault (mJ/m 2 ) 260 314 346 353 Sable Prism stacking fault (mJ/m 2 ) 166 14 5 215 193 Uns t able Prism stacking fault (mJ/m 2 ) 194 287 247 247

PAGE 77

77 Figure 3 1 . Common slip systems found in HCP metals. Reproduced with permission from Ref. [118] .

PAGE 78

78 Octahedral (O) Tetrahedral (T) Crowdion (C) Split (S) Basal Octahedral (BO) Basal Tetrahedral (BT) Basal Crowdion (BC) Basal Split (BS) Figure 3 2 . Interstitial positions in a HCP lattice. The black spheres represent the interstitial, while gray spheres are normal HCP atoms. Sites in the top row lie out of the basal plane, while sites in the bottom row lie within the basal plane.

PAGE 79

79 A B Fig ure 3 3 . Schematic of divacancy positions in a HCP lattice. A) An in plane divacancy . B) An out of plane divacancy. The atoms are black spheres, while vacancies are denoted by open circles.

PAGE 80

80 A B C D F igure 3 4 . Stacking fault maps for b asal (0001) in Zr. A) The stacking fault map from DFT [107] . B) The stacking fault map from EAM [107] . C) The stacking fault map f rom COMB Zr 1 [98] . D) The stacking fault map from COMB Zr 2 [115] .

PAGE 81

81 A B C D Figure 3 5 . Stacking fault maps for p rism in Zr. A) The stacking fault map from DFT [107] . B) The stacking fault map from EAM [107] . C) The stacki ng fault map from COMB Zr 1 [98] . D) The stacking fault map from COMB Zr 2 [115] .

PAGE 82

82 A B F igure 3 6 . The structure of nanocrystalline Zr after a tensile test with 8.8% strain. A) The c ommon neighbor analysis ( CNA ) map [109] after the tensile test . Red spheres are disordered atoms, while light green spheres are normal HCP atoms. B) The structure with n ormal HCP atoms removed. The orange spheres are atoms with coordination different from that of perfect HCP .

PAGE 83

83 CHAPTER 4 MECHANISMS OF Z R SURFACE CORROSION DETERMINED VIA MOLECULAR DYNAMICS SIMULATIONS Overview of Potentials f or Zr O H To date, interatomic potential formalisms for metal and oxide systems have usually been distinct and incompatible with one another due to the different bonding environments . Both metallic Zr and ZrO 2 have been studied us ing molecular dynamics, but a systematic approach capable of probing complex phenomena like corrosion behavior has been lacking. There are a number of classical potentials for ZrO 2 , which include both fixed and variable charge approaches [119 122] . The fixed charge potentials have been extensively used to describe oxygen diffusion in ZrO 2 based materials ; most notably in yttria stabilized zirconia [119,120,123,124] . Advanced v ariable charge potentials capable of describing zirconium/zirconium oxide interfaces have been developed, but none have been extended to encompass the zirconium h ydride system [121,122] . The challenges posed by the complex bonding environment in the Zr O H system necessitates some sacrifices in the flexibility (the scope of applications) of the cited potentials as well as the COMB potential described in this chapter. In particular, certain material properties are poorly described compared to experimental observations. Most notably , the correct phase order s of co mpeting crystal phases of ZrO 2 and ZrH 2 systems are difficult to capture. For ZrO 2 , a monoclinic phase is the experimental ground state at standard conditions. This monoclinic phase is structurally similar to the higher This chapter is based on M. J. Noordh oek, T. Liang, T. W. Chiang , S. B. Sinnott, and S. R. Phillpot , Journal of Nuclear Materials 452 , 285 (2014).

PAGE 84

84 temperature tetragonal and cubic pha ses. The difficulty of capturing the full phase behavior is due, in part, to the small relative energy differences between phases. DFT calculations predict the stabilities of the tetragonal and cubic phases relative to the monoclinic phase to be +3 7 meV/at om and +7 2 meV/atom, respectively [125] . Likewise for ZrH 2 , a ZrH 2 face centered tetragonal structure is the ground state. Experimentally, the fluorite ZrH 2 structure is not observed. However, DFT calculations place the relative energy difference of fluorite to ZrH 2 at only +7 meV/atom [126] . Proper selection of the short r ange cut offs (e.g. Zr Zr, Zr O, Zr H ) is crucial because different values are capable of changing the relative phase stability. Changing the cut off by mere tenths of an Angstrom results in specific bonds becoming included or excluded in the energy calcul ation. In this work, the Zr Zr , O O and H H cut offs were determined prior to the parameterization of the oxide and hydride , which constrains the flexibility during parameterization. This chapter focuses on the development and applications of a COMB poten tial for the Zr O H system. During the fitting procedure, emphasis was placed on the reaction energies of oxygen/hydrogen with the Zr (0001) surface and the enthalpy of formation for various oxide/hydride stoichiometries. Simulations of O 2 and H 2 O depositi on onto Zr surfaces correctly show chemisorption, which allows detailed atomistic analysis of corrosion behavior. The effect of surface structure on corrosion behavior is presented. We have not found it possible to parameterize a COMB potential for either ZrO 2 or ZrH 2 that captures all of the experimental phases while still reproducing the Zr structure. Experimentally, extensive electron diffraction studies have shown that the

PAGE 85

85 fluorite ZrO 2 phase may exist on various Zr surfaces upon oxygen deposition [127 129] . Therefore, we have focused on the fluorite cubic phase a s the ground state for ZrO 2 . Parameterization for the Zr O system was completed first, which constrain s the Zr response to charge transfer in the Zr H system. Thus, the fluorite cubic phase is also the ground state for ZrH 2 . Th is sacrifice in phase order a llow s a better description of Zr surfaces interacting with relatively small amounts of oxygen and hydrogen. Parameterization and Results of the Zr O H Potential Parameterization of the Zr O and Zr H interactions were accomplished with our in house code Potential Optimization Software for Materials (POSMat) [84] . Two optimization techniques, simulated annealing [86] and simplex [85] , were used in conjunction to minimize a least squares functio n determined by differences in the target and calculated values. The final parameter set may be found in Appendix A; this potential is also available in the official distribution of the LAMMPS code [88] . First principles, density functional theory (DFT) calculations in this wor k use d the Vienna ab initio software package (VASP) with the projected augmented wave (PAW) [130] method in the generalized gradient approximation (GGA) with non spin polarized electrons [90,131] . The k point mesh in the Monkhorst Pack scheme [13 2] was 12x12x12 and the energy cut off was 650 eV. The Zr Zr parameters used in the Zr O H potential are from the COMB Zr 2 set, which were fixed during Zr O and Zr H fitting. For Zr O and Zr H fitting, separate training databases were created from experimental and first principles calculations. The initial databases contained a wide range of stoichiometries, from a single point defect i n the Zr metal matrix to ZrO 2 and ZrH 2 . The database included various oxygen and hydrogen

PAGE 86

86 coverage on the Zr (0001) surface as well as all experimentally observed phases of ZrH x . After obtaining a reasonable candidate parameter set, each initial structure was quenched to an energy minimum by a steepest descents algorithm. If the resulting quenched structure had a lower energy than the initial structure, the quenched structure was added to the fitting database and given a target value higher than the initia l structure target value. This iterative method of continually adding structures to the training database helped guide the optimization procedure towards the desired values. In all, the final Zr O database contained 88 structures, while the Zr H database c ontained 93 structures. Bulk Zr O S ystem The bulk properties of cubic fluorite zirconia are l isted in Table 4 1 . With this parameterization, the monoclinic and tetragonal phases both transform into the fluorite structure upon energy minimization. The latt ice parameter for COMB (a 0 = 5.12 5 Ã… ) matches experiment (a 0 =5.086 Ã… ) within 1% error [15] . While the C 44 elastic constant is significantly larger than the DFT value, overall the bulk and shear modulus match the DFT values within 10%. Table 4 2 shows point defect properties from DFT and COMB. The point defect formation energies were calculated by ( 4 1) where E Defect is the total energy of the supercell containing a point defect and E Perfect is the energy of the fluorite supercell. Here, is the chemical potential of species x , which is the atomic energy of Zr an d O in their respective (metallic and molecular)

PAGE 87

87 ground states. Addition of the chemical potential in Eq. 4 1 is used for vacancies, while subtraction is used for interstitials. Our DFT predictions for vacancy formation energy are slightly different than t hose in Ref. [133] , which can be attributed to the different pseudopotentials used. For a given species type, the vacancy energy is much higher than that of the interstitial in DFT. COMB matches this trend, but vacancy formation energies are overestimated by a factor of approximately two. These point defects were included in the fitting database, but the emphasis was placed on having the fluorite structure as the ground state while ensuring positive defect formation energy. A similar sacrifice of the point defect behavior was made in the fitting of metallic Zr. Bulk Zr H S ystem Table 4 3 provides the enthalpies of formation for the various phases in COMB as compared to DFT and experiment. The enthalpies in COMB are much larger than in DFT, but the correct trend is observed with ZrH 2 being the most stable stoichiometry. The crystal structure of ZrH 2 is differentiated from th e fluorite ZrH 2 structure by the c 0 /a 0 ratio. The energy profile for ZrH 2 as a function of c 0 /a 0 ratio shows a local maxima at c 0 /a 0 =1 and global minima at c 0 /a 0 =0.878 [134] . Mechanical stability requirements for the cubic system are [135] : C 11 > 0; C 44 > 0; C 11 > |C 12 |; (C 11 + 2C 12 ) > 0; Thus, DFT predicts fluorite ZrH 2 to be unstable as indicated in Table 4 4 . For the parameterization of COMB, we fit to the lattice parameter, bulk modulus and ensured mechanical stability of the fluorite structure. The se results are achieved for COMB as shown in T able 4 4 . Upon energy minimization, ZrH 2 transforms into the fl uorite structure.

PAGE 88

88 Summary of Zr O H COMB P otential The current Zr metal potential (COMB Zr 2) correctly describes the HCP ground state we ll in terms of bulk, mechanical and surface properties making it well suited for the study of corrosion processes. The ground states for ZrO 2 and ZrH 2 are the cubic fluorite structures, in contrast to experimental observations of structures that are similar to, though distinct from, the fluorite structure. However, the enthalpies of formation for various stoichiometries are in line with first principle calculations, suggesting that the reactivity of zirconium with oxygen and hydrogen should be captured by the COMB potential: the next sections explore the interaction of oxygen and hydrogen with various zirconium surfaces. E nergetics of Oxygen and Hydroge n on Zr (0001) Surface As a test of the capabilities of this integrated set of Zr O H potentials and as a first effort to characterize the mechanism s of corrosion, we first perform static calculations of oxygen and hydrogen interacting with a Zr (0001) ba sal surface, for which DFT dat a is readily available. Tables 4 5 and 4 6 compare the binding energies of different oxygen and hydrogen monolayer (ML) coverage s on a Zr (0001) basal surface for DFT and COMB. Figure 4 1 shows the face centered cubic (FCC) an d hexagonal close packed (HCP) sites on the (0001) surface. The site occupancy nomenclature is taken from Refs. [136] and [137] . One monolayer correspo nds to oxygen/hydrogen coverage equal to the number of surface Zr atoms. A designation of Octa(1,2) corresponds to a layer of oxygen/hydrogen atoms in the octahedral sites situated between the 1 st and 2 nd basal planes of Zr. For oxygen surface coverage, t he FCC sites are more stable than HCP sites in DFT, a result that is reproduced with the COMB potential. The binding energy of 0.5 ML

PAGE 89

89 FCC is 9.0 eV/atom from DFT and 5.6 eV/atom in COMB. While the binding energies are much stronger in DFT, the relative e nergies between sites in COMB gives much mo re satisfactory results (Table 4 5 ). Binding energies for oxygen atoms within the bulk will give an indication whether it is energetically favorable for an atom to diffuse into the bulk. The oxygen binding energy in octahedral site in bulk Zr is 8.77 eV from DFT [136] . This value is slightly higher than the binding energy of a 0.5 ML surface coverage. For COMB, the binding energy of an oxygen atom in the octahedral site is 10.6 eV and 12.0 eV in the tetrahedral site. These values are about twice as large as the binding energy of 0.5 ML on the surface. Thus, it is very energetically favorable for oxygen to diffuse into the Zr bulk. For hydrogen surface coverage, DFT predicts that HCP sites are more stable than FCC sites, with the binding energy decreasing fr om 0.5 ML to 1.0 ML. These tren ds are reversed in COMB (Table 4 6 ). As in the case of oxygen coverage, the binding energies in COMB are higher than DFT and by approximately the same percentage. When a second layer of hydrogen atoms is added between the 1 st and 2 nd planes of Zr atoms, the binding energy is positive, which means the configuration is unstable. The binding energy of a hydrogen atom in bulk Zr is 0.46 eV in the tetrahedral site and 0.38 eV in the octahedral site from DFT [138] . COMB values are 5.3 eV for tetrahedral site and 2.1 eV for the octahedral sit e. Hence, the hydrogen energetics will favor diffusion into the Zr bulk akin to the oxygen case. The static calculations show COMB overestimates the binding of oxygen/hydrogen atoms in the bulk, and underestimates binding energy for surface atoms. These results suggest that although COMB may not reproduce the initial ( 0.5

PAGE 90

90 ML) O/H coverage, the surface morphology at larger coverage should still be captured with the amount of O/H within the first few surface layers of Zr, which will determine the initial oxide/hydride formation. To explore this further , we next simulate some of these scenarios at 300 K for up to 50 ps. Configurations of a single atom on the sur face, 1 ML in FCC surface sites and 1 ML in Octa(1,2) sites were considered. For configurations with only oxygen, the atoms remain in their original sites. In contrast, hydrogen readily diffuses into the bulk. Next, we try mixed configurations of different oxygen and hydrogen layers. The first structure is 1 ML of oxygen placed on FCC surface sites a nd 1 ML o f hydrogen in Octa(1,2) (Figure 4 2A ). The oxygen stays on the surface, while hydrogen diffuses into the bulk. When the species types are switched, hydrogen remains on the surface and ox ygen remains in Octa(1,2) (Figure 4 2B ). Thus, it appears tha t layers of oxygen act as barriers for hydrogen diffusion at the short timescales considered. Experimentally the diffusion of H through the oxide scale is enhanced by cracks and grain boundaries, with diffusion rates of ~t 1/4 to t 3/8 [139] . The simulated structure does not contain defects nor another readily accessible diffusion pathway for the timescales observed. Thus, the presence of a layer of oxygen prohibiting hydrogen diffusion in COMB is not in consistent with expe riment. We also performed nudge d elastic band (NEB) calculations [140] in DFT to quantify the energy barriers for oxygen/hydrogen penetrating from the surface into the first subsurface layer. The barrier is non existent for the tetrahedral site, as a single subsurface oxygen/hydrogen atom will spontaneously exit through the surface. For the octahedral site, the energy barriers are +1.89 eV and +0.74 eV f or oxygen and

PAGE 91

91 hydrogen, respectively. The Arrhenius equation for the rate at which an atom will diffuse to the subsurface is ( 4 2 ) where v 0 is the attempt frequency (~10 13 s 1 ), E AB is the energ y barrier, T is the temperature and k B is the Boltzmann constant [141] . Thus, the relative rate at which hydrogen diffuses is ~10 19 that of oxygen at 300 K. In addition, COMB overestimates the binding energy of a H atom in the bulk, leading to a lower energy barrier that allows H to readily diffuse. Similar routes were taken with the prism and surfaces at 300 K. For basal (0001), the Octa(1,2) site (z = c 0 /4) lies between the planes of Zr atoms. However, f or the prism structures, the Octa(1,2) site is inaccessible because it does not lie directly between the Zr planes for these surfaces. Therefore, crowdion sites , which lie between the prism planes, were chosen. In the case of hydrogen, these surfaces acted in a similar manner as in the basal case, with hydrogen ready diffusing. When oxygen is placed directly on top of a Zr atom it does not diffuse. However, the perfect surface structure is lost by rearrangement of the Zr atoms when oxygen is placed in the c rowdion site. Thus, the surface structure and placement of the oxygen will have a large effect on the resulting corrosion behavior. To better understand this, we also analyze the different surfaces after depositing O 2 and H 2 O. Overview of O 2 and H 2 O Depos ition O 2 or H 2 O atmospheres [4] . Generally, the prism planes corrode much faster than the basal plane [142,143] . However, detailed atomic

PAGE 92

92 level studies showing the structure at short timescales are scarce. Scanning tunneling microscopy (STM) of O 2 atmospheric conditions showed oxide clusters form on the surface with clustering being less dense on (0001) than [144] . At longer times, formation of voids and microcracks at the metal/oxide interface may be responsible for enhanced corrosion behavior [145] . The presence of hydrogen may help reduce the se defects by increasing Zr diffusion to the interface [146] . Due to the wide range of chemistries of Zr alloys in the liter ature, the exact nature of the surface morphology or surface chemistry may be corrosion atmosphere dependent. As a first step in elucidating the role of surface structure and chemistry on the oxidation behavior, we deposited O 2 and H 2 O molecules onto three zirconium surfaces: (0001), and . The supercell dimensions for the metal substrate are approximately 90 Ã… x 84 Ã… x 87 Ã… containing 22 , 500 Zr atoms. The volume is kept fixed and the set 300 K temperature is controlled with a Langevin thermostat f or a depth of ~35 Ã… from the surface. A thermostat does not control Zr atoms below this depth, which allows the kinetic energy from the deposited molecules to be dissipated from the structure. Impinging molecules are given energy of 1 eV, with 0.2 ps relax ation time between each deposited molecule. This energy corresponds to molecule velocities of 24.6 Ã… /ps for O 2 and 32.7 Ã… /ps for H 2 O. A total of 375 molecules were deposited, which corresponds to a coverage of 1 ML for O 2 on Zr (0001). Deposition of O 2 The following analysis is for oxygen atoms that chemisorb, which is ~90 95% under these conditions. Oxygen occupancy on the surface and between Zr planes was

PAGE 93

93 found to be highly dependent on the exposed structure. Figure 4 3 shows snapshots of these surface s after deposit ion has been completed while Figure 4 4 shows a close up view of an oxygen molecule dissociating after 10 ps. The charge transfer i s initially very localized (Figure 4 4A). Upon dissociation of O 2 , the atomic oxygen may diffuse across t he su rface (Figure 4 4B) or into the bulk. The interpla nar distance, Zr planar density and open space available for diffusion all play an important role in the formation of th e initial oxide structures. Figure 4 5 provides the oxygen profile as a function of de pth from the surface. The basal (0001) surface has significantly more oxygen on the surface (~38%), compared with the two prism planes (~ 23%). As the surfaces oxidizes, Zr atoms rearrange themselves as adatoms. This leads to small amounts of oxygen atoms chemisorbing at these reconstructed sites, which corresp onds to the +x axis tail of Figure 5. Figure 4 6 shows the total occupancy of oxygen on the surface and between the 1 st and 2 nd Zr planes. About 87% of oxygen lie within these top two layers for (0001 ), compared with only ~54% for and ~44% for . In all cases, the rest of the oxygen lies below the second Zr plane as shown in the depth profile of Figure 4 5. At 300 K, the bulk equilibrium interplanar spacing is 0.93 Ã… , 1.62 Ã… and 2.59 Ã… for the , and (0001) planes, respectively. Thus, the closed packed structure and large interplanar spacing of the (0001) plane can explain the limited oxygen diffusion as compared to the two prism planes. Correspondingly, we determine the number of Zr at oms that do not have oxygen nearest neighbors; that is, un oxidized atoms. Figure 4 7 illustrates the number of un oxidized Zr within the firs t two Zr planes. The value of ~ 24% for the (0001) surface layers is twice as large as that for the two prism surfa ces. Considering that 87% of

PAGE 94

94 oxygen is in the first two Zr planes (Figure 4 6), this corresponds to a significant clustering of oxygen atoms (Figure 4 3). The surface structure shows a uniform distribution of oxygen in each layer, penetrating to the f ourth Zr planes due to the close interplanar spacing and least dense Zr planar density. The surface structure shows a uniform distribution of oxygen within the top three Zr planes. Hence, the prism and oxide layers are approximately the same thickness. The overall structure of the oxide films on all three surfaces are very similar in terms of coordination and bo nd angles for all surfaces. Figure 4 8 provides the Zr O Zr and O Zr O bond angles for basal (0001). These angles show peaks around 1 00 110 ° , suggesting that oxygen has a tetrahedral bonded environment. Coordination analysis indicates that Zr atoms that have nearby oxygen prefer a coordination of four. Thus, these surfaces show the early formation of the oxide layer, which is the fluori te structure in our model. Deposition of H 2 O For H 2 O deposition, the rate of water molecule dissociation is dependen t on the surface structure. Figure 4 9 shows snapshots show the transfer of charge from Zr to H 2 O, atomic oxygen and atomic hydrogen of the se surfaces after deposition has been completed. Figure 4 10 shows a close up view of a water molecule dissociating after 10 ps. Atomic oxygen diffuses into the bulk while atomic hydrogen diffuses across the surface, which is characteristic for the (0001) surface. We differentiate atomic oxygen from bonded oxygen (H 2 O or OH molecules) by whether the oxygen has at least one positively charged hydrogen neighbor. Atomic hydrogen is negatively charged within the Zr H system in COMB since the electronegativity o f Zr is less than H. Figures 4 11 and 4 12 shows the fraction of atomic oxygen and hydrogen as a function of time for the

PAGE 95

95 different surfaces. For times less than ~24 ps, the prism surface has the least amount of atomic O/H in the structure. For longer times, the basal (0001) surface exhibits much more water and OH content compared to the two prism planes. Thus, a close packed Zr planar density is less likely to dissociate water at long times, while the two prism planes show different behavior at short times. Figures 4 13 and 4 14 provide the occupancy of atomic oxygen and hydrogen in the two outermost Zr planes. Atomic oxygen is more likely to penetrate into the third or fourth Zr layers t han hydrogen. Comparison of Figure 4 6 (O 2 deposition) to Figure 4 13 (H 2 O deposition) shows that the presence of hydrogen significantly increases the penetration of oxygen into the third and fourth layers of Zr for the prism planes. However, the presence of hydrogen has no effect on oxygen penetration in basal (0001). All the bonded oxygens are within top two layers for (0001), while for the prism planes ~65% of chemisorbed bonded oxygen is within the top two planes. The percentage of Zr atoms without a nearby oxygen or hydrogen is similar for H 2 O deposition as the O 2 case (Figure 4 7), meaning that the (0001) surface has clusters of oxygen atoms while the prism planes do not. Unified V iew of O 2 and H 2 O D eposition The corrosion behavior of the basal (0001) surface in zirconium is distinct from its prism surfaces for the conditions of this work. Oxygen and hydrogen are much more likely to remain on the surface of (0001) and experience an uneven spatial distribution d ue to clustering. These oxygen clusters represent highly localized oxide formation and may eventually coalesce into grains at greater oxygen overage similar to those that have been observed experimentally [144] . Additionally, our simulations predict that

PAGE 96

96 water dissociates much more rapidly on the surface than on the (0001) surface. Atomic oxygen and hydrogen diffuse much farther into the surface with a more even distribution than (0001), which may be the cause of the enhanced corrosion rate on [143] .

PAGE 97

97 Table 4 1 . Properties of cubic fluorite ZrO 2 given by the COMB potential compared with experimental data and DFT calculations. Property Experiment [15,147] DFT [125] COMB COMB/ DFT (% error) a 0 ( Ã…) 5.08 6 5.116 5.12 5 0 C 11 (GPa) 417 520 486 7 C 12 (GPa) 82 93 94 +1 C 44 (GPa) 47 61 98 +61 B (GPa) 194 235 228 3 G (GPa) 95 122 134 +10

PAGE 98

98 Table 4 2. Point defect formation energies in zirconia. Units are in eV. Ref. [148] is point defects in the monoclinic phase; Ref. [133] is in the cubic phase. Point Defect DFT (monoclinic phase) [148] DFT (cubic phase) [133] DFT (cubic phase) This Work COMB (cubic phase) O vacancy 6.15 6.58 6.15 11.2 0 O interstitial 1.31 N/A 3.46 0.37 Zr Vacancy 16.44 1 6.14 15.00 27.9 0 Zr Interstitial 3.61 N/A 1.78 4.17

PAGE 99

99 Table 4 3 . Properties of zirconium hydride given by the COMB potential compared with experiments and DFT calculations. Enthalpy unit) Experiment [149] DFT (This Work) COMB Fluorite ZrH 2 N/A 1.69 2.26 ZrH 2 1.69 1.71 N/A ZrH 1.5 1.46 1.2 4 1.60 ZrH N/A 0.85 1.41 Zr 2 H N/A 0.53 1.26

PAGE 100

100 Table 4 4 . Properties of cubic fluorite ZrH 2 given by the COMB potential compared with DFT calculations. Property DFT [134] COMB a 0 ( Ã…) 4.823 4.819 C 11 (GPa) 83 384 C 12 (GPa) 160 21 C 44 (GPa) 20 55 B (GPa) 133 142 G (GPa) 27 106

PAGE 101

101 Table 4 5. Binding energies (eV/atom) of oxygen atoms on Zr (0001) surface. Structure DFT t otal energy [136] COMB t otal energy DFT relative energy COMB relative energy 0.5 ML FCC 9.0 5.6 0.0 0.0 0.5 ML HCP 8.8 5.3 +0.2 +0.2 1.0 ML FCC 8.7 5.0 +0.3 +0.6 1.0 ML HCP 8.5 4.9 +0.5 +0.7 1.0 ML Octa(1,2) 8.5 5.5 +0.5 +0.1 1. 0 ML FCC + 1.0 ML Octa(1,2) 8.4 3.4 +0.6 +2.2 1. 0 ML FCC + 1.0 ML Octa(2,3) 8.6 5.0 +0.4 +0.6

PAGE 102

102 Table 4 6. Binding energies (eV/atom) of hydrogen atoms on Zr (0001) surface. Structure DFT t otal energy [137] COM B total energy DFT relative energy COMB relative energy 0.5 ML FCC 1.00 0.65 +0.08 0.00 0.5 ML HCP 1.06 0.51 +0.02 +0.14 1.0 ML FCC 1.03 0.28 +0.05 +0.37 1.0 ML HCP 1.08 0.28 0.00 +0.37 1.0 ML Octa(1,2) N/A +1.61 N/A +2.26 1. 0 ML FCC + 1.0 ML Octa(1,2) 0.72 +0.35 +0.36 +1.00 1. 0 ML FCC + 1.0 ML Tetra I(1,2) 0.86 +0.46 +0.22 +1.11 1. 0 ML HCP + 1.0 ML Tetra I(1,2) 0.67 +0.45 +0.41 +1.10

PAGE 103

103 A B Figure 4 1. A schematic of the Zr (0001) surface . A) O xygen or hydrogen ( green spheres) is located at the FCC sites . B ) Oxygen or hydrogen (green spheres) is located at the HCP sites. The bl ue spheres are surface Zr atoms, while gold spheres are Zr subsurface atoms.

PAGE 104

104 A B Figure 4 2. Snapshot of Zr (0001) surface viewed along the direction after 50 ps at 300 K. A ) Initial sites are surface FCC for oxygen and subsu rface Octa(1,2) for hydrogen. B ) Initial sites are surface FCC for hydrogen and subsurface Octa(1,2) for oxygen. Atoms are color coded by charge.

PAGE 105

105 A B C Figure 4 3. Snapshots of Zr surfaces after O 2 deposition. A) The (0001) surface. B ) The surface. C ) The surface . Atoms are color coded by charge. A less uniform spatial distribution of oxygen atoms is seen on the (0001) surface as compared to and .

PAGE 106

106 A B Figure 4 4. Dissociation of an O 2 molecule on Zr (0001). A ) The O 2 molecule is intact on Zr (0001) . B ) After 10 ps the O 2 molecule has dissociated into atomic oxygen. Atoms are color coded by charge.

PAGE 107

107 Figure 4 5. Oxygen profile in Zr as a function of depth after deposition of O 2 . The dashed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. The original Zr surface location is zero on the x axis.

PAGE 108

108 Figure 4 6. Fraction of oxygen atoms on the surface and within the first subsurface Zr layer after O 2 deposition . The dashed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t= 0 ps is when the first molecule interacts with the surface.

PAGE 109

109 Figure 4 7. Fraction of Zr atoms in the top two planes that have no nearby oxygen atom after O 2 depos ition . The dashed b lack, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t=0 ps is when the first molecule interacts with the surface.

PAGE 110

110 Figure 4 8. Bond angle distribution in Zr (0001) after O 2 deposition. The O Zr O (dotted black line) and Zr O Zr (solid red line) bond angles are shown . The solid blue spikes represent O Zr O angles in the fluorite structure, which also has a single peak at 109 o for Zr O Zr angles. The simulated structure has O Zr O peaks at ~ 70 o , ~ 100 o , and ~ 170 o . O Zr O angles less than 60 o represent oxygen molecules that have not dissociated.

PAGE 111

111 A B C Figure 4 9. Snapshots of Zr surfaces after H 2 O deposition. A) The (0001) surface. B ) The surface. C ) The surface . Atoms are color coded by charge. The surface s contain varying amounts of H 2 O, OH, atomic O and atomic H.

PAGE 112

112 A B Figure 4 10. Dissociation of a H 2 O molecule on Zr (0001) . A) The H 2 O molecule is intact on Zr (0001). B ) After 10 ps the H 2 O molecule has dissociated into atomic oxygen and hydrogen. The oxygen diffuses into the bulk while hydrogen stays on the surface. Atoms are color coded by charge.

PAGE 113

113 Figure 4 11. Fraction of atomic oxygen in the Zr structure after H 2 O deposition . The dash ed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t=0 ps is when the first molecule interacts with the surface.

PAGE 114

114 Figure 4 12. Fraction of atomic hydrogen in the Zr structure after H 2 O deposition . The dashed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t=0 ps is when the first molecule interacts with the surface.

PAGE 115

115 Figure 4 13. Fraction of atomic oxygen in the two outermost Zr planes after H 2 O deposition . The dashed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t=0 ps is when the first molecule interacts with the surface.

PAGE 116

116 Figure 4 14. Fraction of atomic hydrogen in the two outermost Zr planes after H 2 O deposition . The dashed black, dotted green and solid red lines represent the (0001), and surfaces, respectively. Time t=0 ps is when the first molecule interacts with the surface.

PAGE 117

117 CHAPTER 5 ATOMISTIC STRUCTURE OF (B A ,S R )T I O 3 : COMPARING MOLECULAR DYNAMICS SIMULATIONS WITH STRUCTURAL MEASUREMENTS Introduction to ( Ba ,Sr) TiO 3 Ferroelectrics BaTiO 3 based systems find use in many commercial devices including ceramic capacitors and therm al switches . Upon cooling, BaTiO 3 , which crystallizes with a perovskite like structure, undergoes a series of polymorphic transiti ons from the paraelectric cubic (C) to ferroelectric tetragonal (T) to orthorhombic (O) to rhombohedral (R) phases. The key structural feature of these transitions is an off centering of Ti atoms within the oxygen octahedra. to exhibit a mixed displacive and order disorder character [42,43,45] , though its exact nature still remains under debat e. Overview of MD Simulations f or BaTiO 3 The most common methods for simulating BaTiO 3 are the classical shell model [47,150] , effective Hamiltonians [48,50,151] and density functional theory (DFT) [152,153] . The formalism of an effective Hamiltonian is a Taylor expansion based on the energy landscape of small perturbations [154] . These functions are based on the relevant degrees of freedom of the system , usually denoted by strains and lattice Wannier functions [154,155] . Both the shell model and effecti ve Hamiltonians require parameterization, which typically uses data from first principle calculations. Previous molecular d ynamics (MD) simulations of BaTiO 3 [48,49,150] came to conflicting conclusions on th e transition mechanism. For example, one study [48] This chapter is based on M. J. Noordhoek, V. Krayzman, A. Chernatynskiy, S. R. Phillpot, and I. Levin, Applied Physics Letters 103 , 022909 (2013).

PAGE 118

118 inferred an order disorder character from the saturation of the soft mode frequency near the Curie temperature . By contrast, another report [49] found Ti atoms vibrate around the centers of the cubic cells, implying a displacive type transition , with the Ti off centering attributed to the pref erential oxygen motion along the Ti O bonds. The actual mechanism remains elusive because comprehensive quantitative comparison s of the structural parameters from simulations with those from the appropriate experimental data are nearly nonexistent . Moreove r, s ystematic comparison of various MD simulations [ 48 51,150,151] is difficult because of the limited structural information provided in the literature ; for example, in many cases, not even the lattice parameters of the simulated models were reported. The ability of MD simulations with classical potentials to reproduce simultaneously all the fundamental structural features in BaTiO 3 and i ts solid solutions, such as (Ba,Sr)TiO 3 , a prerequisite for predicting the lattice dynamics and functional properties, has not yet been established. Such fundamental features include (i) dissimilar patterns of local and average Ti off centering within the [TiO 6 ] octahedra, (ii) extended one dimensional correlations of displacements that yield { 100 } sheets of diffuse scattering, (iii) modality/anisotropy of probability density distributions (PDD) of the Ti and O atoms and, in the solid solutions , (iv) local coordination of distinct species mixed on the same crystallographic sites and the effects of this mixture on features (i) (iii). All these structural characteristics are encoded in various diffraction and spectroscopic signals, which can in turn be calcula ted directl y from the atomic coordinates. Thus, quantitative comparison of the simulated and measured signals is a critical step in validating structural pr edictions from MD simulations.

PAGE 119

119 Reverse Monte Carlo (RMC) is a method that uses experimental data as inputs, such as the pair distribution function (PDF), Bragg diffraction and extended X ray absorption f ine structure (EXAFS) spectra, and calculates a best fit for the structural properties [156] . Recent advances in structural refinements using a RMC algorithm enable explicit experimental determination of the relevant displacement characteristics [156,157] , which can be used for an even more comprehens ive validation of simulations. RMC refinements, which are frequently under constrained, would also benefit from this comparison because the MD simulations could provide sensible starting atomic configurations and physics based constraints. In this chapter , we explore integration of MD simulations and structural measurements by confronting the simulation results for Ba 1 x Sr x TiO 3 (0 x 0.5) with five types of experimental data . In addition to a direct overlay of the simulated and measured signals, we also compare key simulated structural parameters to those derived from the RMC refinements . The simulated structures reproduce many aspects of the experimental results remarkably well , while the agreement for other parameters is less satisfactory. Further analysis of the simulations provides insights into (i) a possible mechanism of dielectric relaxation and (ii) chemically resolved lo cal distortions in ( Ba , Sr ) TiO 3 solid solutions. Experimental and Simulation Procedure The experimental data used in this study included ( i ) total pair distribution function s (PDF) determined from neutron scattering, ( ii ) neutron Bragg profile s , ( iii ) loc al Ti off centering determined from the pre edge peak in the Ti K edge X ray absorption spectra (XAS) , ( iv ) diffuse scattering patterns in electron diffraction and ( v ) Sr K edge extended X ray absorption f ine structure (EXAFS) spectra. All the experimental data

PAGE 120

120 were collected on the same ceramic samples by Levin et al. [158] . Direct comparison between MD and the experimental data i s provided in the next sections. These experiments provide detailed structural information such as atomic positions and thermal displacements, which are used as inputs for RMC refinement. W e perform structural refinements for tetragonal BaTiO 3 and Ba 0.8 Sr 0 .2 TiO 3 using simulated atomic configurations of 16×16×16 unit cells (20,480 atoms) in the RMCProfile software [156] . The data used for RMCProfile fitting includes the five types of experimental data stated above, as well as Ti off centering and the macroscopic value of polarization [37,157] . MD s imulations are performed usin g the DL_POLY software package [89] . Interatomic interactions are described using a classical shell model for the potential developed for Ba 1 x Sr x TiO 3 [150] with parameters adopted from Sepliarsky et al . [73] . The value s of the parameters are shown in Appendix B. In this work, the adiabatic shell model is used in the simulations (Chapter 2). A fraction (10%) of the atomic mass is assigned to the shells so that the core and shell motions can be numerical ly integrated usin g the leapfrog method [ 17 ]. Core shell interactions are described by an anharmonic spring according to V( w )=c 1 w 2 /2+c 2 w 4 /24, where w is the core shell separation. Short range shell shell interactions between shells are described by Buckingham terms V( r ) = A exp( r / ) C / r 6 . Long range Coulombic interactions are calculated using the Smoothed Particle Mesh Ewald method [72] . All simulations are performed for atomic configurations of 16×16×16 unit cells under periodic b oundar y conditions. A Berendsen thermostat and barostat are applied, allowing anisotropic variat ions of the lattice parameters [62] . A timestep of 0.4 fs is used. Structures are equilibrated for at least 50 ps before recording the configurations

PAGE 121

121 for analys e s. For all time dependent properties, configuration snapshots are recorded every 20 fs over a 20 ps period. The configurations for Ba 1 x Sr x TiO 3 ( x =0.2, 0.5) are generated by replacing an appropriate fraction of randomly selected Ba atoms with Sr. The pres ent potentials underestimate the tr ansition temperatures by 60 K. The phonon dispersion curves calculated for BaTiO 3 are similar t o those reported previously [150] . BaTiO 3 Results Table 5 1 compares calculated and experimental values of the lattice parameters for tetragonal BaTiO 3 at 300 K . The potential overestimates the c parameter , thus yielding a significantly enhanced c / a ratio (Table 5 1 ). The PDF calculated for the simulated tetragonal BaTiO 3 ( 300 K ) reproduces the experime ntal signal satisfactorily (Figure 5 1 A ). In this plot, the calculated PDF corresponds to the MD configuration with the lattice parameters scaled to their experimental va lues in order to facilitate comparison between the two datasets. Similar to the experiment, the 1 st Ti O peak in the calculated PDF (inset in Figure 5 1A ) exhibits a well defined doublet, which is commonly regarded as a signature of the local rhombohedral like distortion of [TiO 6 ] octahedra . The effective local Ti O bond distances were obtained by fitting the 1 st Ti O peak in the simulated and experimental PDFs using a parameterized monoclinic ( Cm ) model that allows for rotation of the local polarization ve ctor in the (110) plane. Evidently, both simulations (Ti O distances: 1.831 Å , 1.918 Å (×2), 2.088 Å (×2), 2.229 Å ) and experiments (1.831 Å , 1.922 Å (×2), 2.077 Å (×2), 2.241 Å ) yield similar rhombohedral like distortions with 3 short and 3 long Ti O distances. Relative peak intensities in the neutron Bragg profile of BaTiO 3 (300 K) are reproduced well (Figure 5 1B ) , after again accounting for the differences between the simul ated and ex perimental lattice parameters . The peak intensities are determined by

PAGE 122

122 the local variations of atomic position in the global structure. The MD simulations underestimate the magnitudes of atomic thermal displacements, as reflected in the stronger intensities of the calculated Bragg peaks. Thus, the repulsive component of the potential is slightly too large, leading to the smaller thermal displacements. Simulations performed at temperatures closer to the T C transition (e.g., 380 K) returned similar results. The { 100 } diffuse scattering sheets in electron diffraction, which are regarded as a signature feature of the BaTiO 3 structure [39] , are also reproduced (Figure 5 1C ). This diffuse scattering originates from extended displace ment correlations (positive Ti Ti and negative Ti O) along the 100 Ti O Ti chains. The simulated correlation paramete rs and correlation lengths (Figure 5 1D ) agree with those determined using the RMC refinements. A shorter correlation length encountered in the experiments is likely an artifact of under constrained RMC fits. Local off centering of Ti atoms, Ti , can be determined from the intensity of the Ti K pre edge peak in XAS [41] . For BaTiO 3 , both the total magnitude of this off centering and the ratio of its mean squared c and a axis components, = Ti 2 c / Ti 2 a , ha ve been measured as a function of temperature [41] . T he total value of Ti for the simulated structure (0.26 Ã… ) is close to the experimental value s of 0.23(1) Ã…. However, the simulated ratio =5 is considerably larger than 2.0 determi ned from single crystal XAS [41] . The RMC refinements yield 1.8. These values signify that the off centering of the Ti atom in MD is significantly overestimated along the c axis as compared to the a axis, in contr ast to the experimental data. This is caused in part by the exaggerated tetragonal distortion in MD.

PAGE 123

123 The simulated directions of the local Ti displacements, relative to the center of gravity of the Ba positions, reveal clear preferences for the (110) pla nes (i.e., a directional probability in the (001) plane exhibits well defined peaks at ± 45 o to the a axes). A similar trend is observed in the RMC refined structures. Consistent with the experimental results, the directional probability within the (110) pl anes displays peaks away from the c axis; however, in simulations these peaks occur at ± 7.5 o rather than at ±20 o determined from the RMC refinements. Thus, the simulations reproduce the qualitative features of the Ti displacements but significantly underestimate the a axis displacement components, as also reflected in the overestimated value of and the incorrect sh ape of the Ti PDD discussed below. The oxygen PDDs in the simulated BaTiO 3 structures are strongly anisotropic and flattened along the Ti O bonds in qualitative agreement with the RMC refinements, though the simulations exaggerates the flattening. This ani sotropy contrasts with the one proposed by Matsko et al . in which oxygen PDDs are flatted perpendicular to the Ti O bond [49] . The simulated Ti PDD in the ( xy ) plane (Figure 5 2A ) of the tetragonal phase is clearly unimodal. In contrast, the corresponding experimental T i PDD exhibits a broad top (Figure 5 2A ), which is consistent with a bimodal distribution ex pected for the 4 site model [39] . Thus, splitting of the Ti sites is not reproduced in the simulations. The simulated and experimental Ti PDDs agree on unimod al profiles along the c axis (Figure 5 2B ) . Previous MD studies [151] inferred the site hopping me chanism for Ti atoms from the occasional short lived (40 fs) reversals of the p z component of local polarization in . We observe similarly brief, though less frequent,

PAGE 124

124 dynamic reversals of p z at 350 K as shown in Figure 5 3 . A conclusive discrimination between dynamic site hopping and usual thermal displacements, which may also lead to occasional reversal of polarization, is possible only if the atoms reside at the excited sites significantly longer than th e characteristic phonon times. In BaTiO 3 , this discrimination is difficult, if even possible, because of the short relaxation time (50 fs) implied by the bare relaxation mode frequency of 100 cm 1 reported fro m spectroscopic measurements [151] . Thus far, the site hopping of Ti atoms remains unconfirmed. Considering the existence of extended displacement correlations a long the 100 s of gravity. The behaviors of the Ti Ti Ti and O O O sub chains of the same Ti O Ti O chain are considered separately. Figure 5 4 A demonstrates that the centers of the a axis Ti sub chains hop between two well defined positions separated by 0.08 Ã… with a characteristic frequency of 0.5 THz. A similar behavior is observed for the O sub chains, which move in anti phase with respect to their Ti counterparts, yielding a dy namic reversal of polarizat ion along the chain direction. In contrast, the c axis chains exhibit no hopping. Unlike a single peak Ti PDD, a distribution of the a axis chain center positions ( Figure 5 4 B ) exhibits a bimodal distribution consistent with a 4 site model (splitting along 100 T directions). In the cubic phase simulated at 440 K, such bimodal distributions are observed for the chain centers along all three cub ic axes (i.e., 6 split sites). The 6 site chain center and 8 site atomic models of dynam ic disorder involve distinct characteristic timescales with different implications for the high fre quency dielectric relaxation [151] . Namely, the 6 site chain center model corresponds to a displacive type transition, while the 8 site model corresponds to an order disorder transition. Further improved

PAGE 125

125 simulations, which reproduce the multiple site Ti PDD identified by the RMC refinements, are needed to differentiate between these two mechanisms. Thus, the MD simulations for BaTiO 3 replicate nearly all of the measured data: some quantitativ ely, some semi quantitatively. The principal deficiency of these simulations is in the strongly overestimated c axis anisotropy of the Ti off centering, which is manifested also in the incorrect unimodal appearance of the Ti PDD contrary to the 4 site Ti distribution obta ined from the RMC refinements. Below we perform a similar comparison of simulations and structural data for the Ba 1 x Sr x TiO 3 solid solutions followed by analyses of the local Ba and Sr coordination and the effects of the Ba/Sr disorder on the local Ti off centering. Ba 1 x Sr x TiO 3 Results The lattice parameters for the Ba 1 x Sr x TiO 3 compos itions are summarized in Table 5 1 . Similar to BaTiO 3 , the simulations for the tetragonal x =0.2 overestimate the distortion. For x =0.5, the structure is cubic. For x =0.2, while most of the features of the PDF are reproduced satisfactorily (Figure 5 5 A ), the simulated PDF exhibits a significantly narrower 1 st Ti O peak with a less pronounc ed peak splitting (inset in Figure 5 5 A ). That is, the simulations underestimate the local distortion for the [TiO 6 ] octahedra in the tetragonal solid solutions. However, for cubic x =0.5, the agreement is excellent (Figure 5 5 B and inset). The discrepancies between the simulated and experimental Bragg intensities are similar to those in BaTiO 3 for both x =0.2 and x =0.5 and are caused primarily by the relatively small thermal parameters in the simulations. The total Ti off centering in the simulated structures decreases from Ti =0.21 Ã… for x =0.2 to Ti =0.15 Ã… for x =0.5, in good agreement with the available experimental

PAGE 126

126 v alues (0.20 Ã… for x =0.2 and 0.15 Ã… for x =0.5) [159] . For x =0.2, the anisotropy ratio for the Ti off centering =1.8 also agrees well with =1.6 determ ined from the RMC refinements. Likewise, the Sr EXAFS for the simulated configuration [160] reproduces t he experimental EXAFS data (Figures 5 5 C, D ) surprisingly well, especially for x =0.5, for which the agreement is nearly perfect. For Ba 1 x Sr x TiO 3 (x=0.2, 0.5), the diffuse sca ttering sheets are retained [37] , consistent with experiments. The displacement correlations (Ti Ti and Ti O) t hat underlie this diffuse scattering are similar to those in BaTiO 3 . Likewise, the PDDs of the Ti chain centers feature bimodal distributions though with a smaller peak splitting than observed in BaTiO 3 . The average positions of the Sr atoms in the simul ated structure of Ba 0.8 Sr 0.2 TiO 3 are displaced along the c axis relative to the Ba atoms (in the same direction) by 0.08 Ã…. The distributions of the Sr O distances are broad and shifted to smaller distances relative to th ose for the Ba O distances (Figure 5 6 A ), as also observed in the RMC refinement s (not shown). The relaxation of the [SrO 12 ] coordination appears to be approximately isotropic. Both the MD simulations and RMC refinements reveal a systematic decrease in the Ti off centering with the increasing local Sr/Ba ratio around the Ti atoms (Fig ure 5 6 B ), which explains a decrease in the total Ti off centering with increasing x . This trend correlates with the smaller volumes of the oxygen octahedra surrounded preferentially by the Sr atoms, as observed both in the simu lations and RMC refinements [37] . Conclusion s on Confronting MD Simulations w ith Experimental Results In summary, we have presented a systematic compa rison of the results of MD simulations in (Ba,Sr)TiO 3 to a comprehensive suite of experimental structural data.

PAGE 127

127 Various simulated and measured signals were compared in addition to specific structural parameters (Table 5 2 ). For several signals, the agre em ent was excellent (e.g., Fig ure s 5 1, 5 5 B , 5 5 D ). The simulated models reproduced many key features of the local structures, although the agreement for several important parameters is still imperfect. The simulations combined with the structural refinemen ts yielded new insights into the local structures of the (Ba,Sr)TiO 3 solid solutions. In addition, the simulations identified anti correlated hopping of the Ti Ti Ti and O O O chains as a potential mechanism for the terahertz dielectric relaxation in t hese materials. Our results highlight the importance of detailed validations of simulated structures before the resulting materi al properties can be analyzed. For example, a dynamic response of polarization in BaTiO 3 cannot be attributed to the Ti site hop ping until the existence of the split Ti sites in the simulated structure is established. The detailed analyses performed here constitute an extremely rigorous test of the interatomic potential. The relatively simple potential used in the present simulatio ns reproduced a wide range of structural properties rather well. Future development of even more materials specific interatomic potentials for ferroelectric materials should use training database s that combine the experimental methods used in this study with results of density functional theory calculations , which can provide high accuracy structural information such as forces, stre sses and energies for each atom. This parameterization strategy has been used with success for other oxide syste ms [161,162] . Additional fidelity might also be expected by the combination of such advanced training sets wi th more advanced potentials [163] in which the ionic charges and bonding strength dynamically vary in respo nse to differences in the local

PAGE 128

128 environment. This strategy is followed in the next chapter, albeit with very limited success .

PAGE 129

12 9 Table 5 1. Experimental and simulated lattice parameters for Ba 1 x Sr x TiO 3 . The experimental values for all compositions are quoted for 300 K. Numbers in parenthesis refer to a single standard deviation. The MD values correspond to 300 K for x =0 and to 220 K for x =0.2 and x =0.5. x Experiments MD a c c / a a c c / a 0 3.991(1) 4.034(1) 1.011(1) 3.987(1) 4.073(1) 1.022(1) 0.2 3.974(1) 3.999(1) 1.006(1) 3.974(1) 4.028(1) 1.014(1) 0.5 3.949(1) 3.954(1)

PAGE 130

130 Table 5 2. Agreement between the selected simulated and experimental structural characteristics in Ba 1 x Sr x TiO 3 (E excellent, S satisfactory, P poor). Parameter x =0 x =0.2 x =0.5 Distortion of [TiO 6 ] E S/P E d Ti E E E a = Ti 2 c / Ti 2 a P E E Ti PDD P S N/A Ti Ti/Ti O correlations S S N/A Sr O/Ba O distances N/A S E

PAGE 131

131 A B C D Figure 5 1. Comparison of various experimental and simulated structural properties of BaTiO 3 at 300 K . A ) Experimental (red) and calculated (blue) neutron PDFs . B ) Experimental (red) and calculated (blue) Bragg profiles . C ) Experimental electron diffraction pattern a long the [401] zone axis. The inset in the lower left quadrant displays the diffuse scattering calculated for the MD simulated structures. D ) C orrelations (filed symbols simulated, open symbols experimental) in Ti Ti and O O for the Ti O displacements along the Ti O Ti chains that give rise to t he diffuse scattering .

PAGE 132

132 A B Figure 5 2. The Ti PDD in tetragonal BaTiO 3 at 300 K. A) T he Ti PDD for simulated (blue) and RMC refined (red) models along the [110] direction. B ) The Ti PDD for simulated (blue) and RMC refined (red) models along the [001] direction .

PAGE 133

133 Figure 5 3. Local p z polarization component in MD as a funct ion of time for BaTiO 3 at 350 K.

PAGE 134

134 A B Figure 5 4 . Simulated Ti Ti Ti chain data in BaTiO 3 . A ) The t ime dependent x coordinate for the center of gravity of the a axis Ti Ti Ti chain. B ) The PDD (along the a axis) for the Ti Ti Ti chain center.

PAGE 135

135 A B C D Figure 5 5 . Comparison of experimental and calculated (MD) data for Ba 1 x Sr x TiO 3 . A) Experimental (red) and calcu lated (blue) neutron PDFs for Ba 0.8 Sr 0.2 TiO 3 . B). Experimental (red) and calculated (blue) neutron PDFs for Ba 0.5 Sr 0.5 TiO 3 . C) Experimental (red) and calculated (blue) Sr EXAFS for Ba 0.8 Sr 0.2 TiO 3 . D) Experimental (red) and calculated (blue) Sr EXAFS for Ba 0.5 Sr 0.5 TiO 3 .

PAGE 136

136 A B Figure 5 6 . Bond length and Ti off centering data for simulated Ba 1 x Sr x TiO 3 . A ) The Sr O (red) and Ba O (blue) distance distributions in Ba 0.8 Sr 0.2 TiO 3 (solid line) and Ba 0.5 Sr 0.5 TiO 3 (dashed line). B ) Dependence of the local Ti off centering on the local Sr/Ba ratio for Ba 0.8 Sr 0.2 TiO 3 and Ba 0.5 Sr 0.5 TiO 3 .

PAGE 137

137 CHAPTER 6 DEVELOPMENT OF SHELL MODEL POTENTIAL FOR ( BA,CA ) TI O 3 Approach to Potential Development This chapter covers the development of a shell model po tential for ( Ba ,Ca) TiO 3 using two different approaches. The first is to extend the Ba 1 x Sr x TiO 3 potential by Tinte et al. [150] to encompass Ba 1 x Ca x TiO 3 . The original parameters are taken from Ref. [73] and may also be found in Appendix B . This approach is advantage ous because the parameterization of BaTiO 3 is complete and its structural deficiencies have been well quantified as shown in Chapter 5 . The second method involves a completely new parameterization for BaTiO 3 in an attempt to find better agreement with expe riment. Once completed, Ca interactions will be added and fine tuned to match experimental properties. Extension of Ba 1 x Sr x TiO 3 Potential to Ba 1 x Ca x TiO 3 The results in Chapter 5 show that the shell model is capable of describing many important structural aspects of Ba 1 x Sr x TiO 3 ferroelectrics. In an attempt to extend the previous potential [150] to Ba 1 x Ca x TiO 3 , parameterization for Ca based interactions is performed. The original parameters are kept fixed during the fitting process to ensure full transferability between the BaTiO 3 , SrTiO 3 and CaTiO 3 sys tems. The Ca O two body interactions and Ca shell constant are fit to the orthorhombic and cubic phases of Ca TiO 3 . Since Ti O and O O interactions are fixed, the available flexibility to describe Ca TiO 3 is very limited. If values for the Ca spring constant are too small, the shells are easily displaced from the cores and the material is unstable in MD simulations. In addition, we found that a Buckingham term alone for Ca O interactions is insufficient to accurately describe CaTiO 3 . Thus, a Morse term is added to help capture the

PAGE 138

138 anharmonic vibrations . The critical value of note during the fitting process is the gamma point frequency in the cubic phase. Fine tuning to a specific gamma point frequency value helped guide the parameterization to the desired results. The parameters developed for Ca may be foun d in Appendix B. Table 6 1 shows the structural properties of orthorhombic (Pbnm) Ca TiO 3 at 0 K. T he pair distribution function (PDF) at 300 K is shown in Figure 6 1. The MD results satisfactorily match experiment with all the relevant peaks appearing at the appropriate locations. As noted above, the Ca TiO 3 and Ba 1 x Ca x TiO 3 properties are very sensitive to the gamma point frequency in the cubic phase. The Ca TiO 3 cubic phase has a imaginary gamma point fre quency of 142i cm 1 , which i s very close to the value of 141 i cm 1 from DFT [164] . A composition of Ba 0.8 Ca 0.2 TiO 3 is used to study the ferroelectric properties and provide a direct comparison to those in Ba 0.8 Sr 0.2 TiO 3 [159] . Figure 6 2 shows the PDF and Ca EXAFS for Ba 0.8 Ca 0.2 TiO 3 . Th e bond distances in MD are slightly shifted to shorter distances as compared to experiment. This is made clearer by showing only t he first O Ba and O Ca peaks as seen in Figure 6 3 [165] . The O Ba bonds lengths are quite similar for MD and experiment, but the O Ca peak in MD is shifted towards much shorter distances. Thus, the O Ca bond distance in Ba 0.8 Ca 0.2 TiO 3 for MD is similar to the O Ca bond distance in pure Ca TiO 3 (Figure 6 3) . This is caused by the location of Ca in the tetragonal unit cell. Reitveld refinements show that Ca is positioned at (0, 0, 0.029) in experiment [165] and (0, 0, 0.043) is MD for Ba 0.8 Ca 0.2 TiO 3 . This overestimation of the Ca position by nearly 50% leads to a relaxed O Ca bond, which contra sts with experiment. Figure 6 4 shows the probability density distributions ( PDDs ) for Ca and Sr position within the tetragonal unit cell for Ba 0.8 Ca 0.2 TiO 3 and

PAGE 139

139 Ba 0.8 Sr 0.2 TiO 3 , respectivel y. The Ca position along the z direction is much larger than for the Sr case, which contrasts experimental observation that their positions are approximately the same [159] . Figure 6 5 shows the c/a ratio as a function of composition for BaTiO 3 , Ba1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 . T he tetragonal distortion was captured with the addition of small amounts of Ca to BaTiO 3 (< 30%). We saw in the previous chapter that the addition of Sr to BaTiO 3 reduces the c/a ratio, both experimentally and in MD . Experime ntally, the addition of Ca to BaTiO 3 does not decrease the c/a ratio, unlik e the effect of Sr. The MD results show that the c/a rat io for Ba 1 x Ca x TiO 3 decreases with increasing Ca composition . The decreased c/a ratio in MD is related to the O Ca bond lengths discussed above . In MD, O Ca bonds are relaxed, which allows the Ti to b e shifted closer to the unit cell center thereby reducing the c/a ratio. Thus, additional repulsion between O Ca is likely needed in the potential, but this modification is not possible du e to the fitting database limitations. However, the c/a r atio for Ba 1 x Ca x TiO 3 is still greater than Ba 1 x Sr x TiO 3 for MD in accordance with experiment. Another feature of Ba 1 x Ca x TiO 3 in experiments is the nearly constant tetragonal to cubic transition temperature (T c ) with increa sed Ca concentration, which again is related to the O Ca bond length discussion [37] . Figure 6 6 shows T c as a function of concentration for BaTiO 3 , Ba 1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 . The MD results systematically underestimate T c for each system as compared to experiment. Obtaining accurate T c in MD has proven problematic for other works as well [47,166] . The transition temperatures for the solid solutions are restric ted from being greatly modified, because Ti O interactions are fixed. As shown in the next section, the energy versus

PAGE 140

140 displacement of Ti within the unit cell will affect the resulting T c . In addition to the underestimated absolute value of T c , the trend of T c decreases with increasing Ca composition for Ba 1 x Ca x TiO 3 in MD. This contrasts with experiment, where T c stays constant [37] . Thus, the c/a ratio and T c is greater for Ba 1 x Ca x TiO 3 than Ba 1 x Sr x TiO 3 , but the trends of these values do not match experiment. The Gibbs free energy of mixing signifies whether the addition of Ca/Sr int o BaTiO 3 is thermodynamically possible . A negative Gibbs free energy indicates that the structure is stable , while positive values indicate an unstable chemical configuration . Experimentally, th e maximal solubility of Ca in BaTiO 3 is ~30%. As a check on the potential, we calculated Gibbs free energy using [3] : (6 1) The enthalpy is defined as (6 2) while the entropy given by the ideal solution theory a s: (6 3) Here, k B is the Boltzmann constant and x is the composition. The subscripts BT, CT and BCT represent BaTiO 3 , CaTiO 3 and Ba 1 x Ca x TiO 3 , respectively. The calculations performed at 300 and 450 K show that the Gibbs free energy is negative for all Ca compositions. Thus, our model shows that BaTiO 3 and CaTiO 3 have complete solubility in one another, which contrasts with the experimental observations. From the aforementioned data, the properties of Ba 1 x Ca x TiO 3 are quantitatively too similar to that of Ba 1 x Sr x TiO 3 to be used for a full study. For each system, the c/a ratio and T c decrease with increasing Ca/Sr composition. In addition, each is soluble

PAGE 141

141 across the entire composition range. Given these qualitative failings, we attempt a completely new parameterization for BaTiO 3 ; this is the subject of the next section. Development of New Potential f or BaTiO 3 In this section, an attempt is made to develop a new shell model potential for BaTiO 3 , in an effort to later include Ba 1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 . The numerous challenges in completing this task are outlined in this section. These include creating a suitable training database, determination of en ergy differences between phases and computational limitations. As shown below, considerable effort in altering the database is required to achieve a reasonable potential. The target inputs for fitting are deduced from DFT calculations and from intuition of the key physics of the system. The resulting potential produced acceptable static properties for BaTiO 3 , but failed provide reliable data during dynamics. Training Database The creation of a training database is often straightforward. Various str uctures are included that allow calculation of properties such as cohesive energy, phase order, surface energy and elastic properties . Data is taken from experiments, first principle calculations or a mixture of both. In this work, multiple training databases were de veloped because the fitting process is very sensitive to the spatial location of atoms within the unit cell. By definition, the ferroelectric nature of these materials is due to the symmetry breaking distortions in a unit cell. For BaTiO 3 , these distortion s are very small, with bond lengths varying by a little as a tenth of an Angstrom between phases. Thus, the database had to be modified by optimizing the internal coordinates of the atoms.

PAGE 142

142 Two different databases were created for the fit ting process. The first is from experiment al data given in Ref. [31] , while the second is created by optimizing the internal coordinates of the atoms. For the latter method, several assumptions involving Ti O and O O bond lengths were made based on the experimental structural data as shown in Figure 6 7 [31] . The assumptions are: The shortest O2 O2 bond length in the orthorhombic phase is equal to the shortest O O bond length in the rhombohedral phase. The shortest O1 O2 bond length in the orthorhombic phase is equal to the short est O1 O2 bond length in the tetragonal phase. The second shortest Ti O bond lengths in the orthorhombic and tetragonal phases are equal. The longest O1 O2 bond length in the tetragonal phase is equal to the longest O O bond length in the rhombohedral phas e. S imulated annealing is used to minimize a least squares function based on the above criteria. These assumptions appear to be reasonable, as the parameter set described below in based on fitting to this optimized database. The final internal coordinates of t he atoms are given in Tables 6 2 through 6 5 . Once the database is created, the importance of the various properties in the fitting database must be prioritized . For instance, since the main interest in this project is the structural properties of the tetragonal and cubic phases, the energy difference between these two must be correct in order for the phase transition to occur near the experimental transition temperature. The energy difference between the tetragonal and cubic phases is 5.6 meV in DFT [167] . Meanwhile, the Tinte et al. classical potential has a tetragonal to cubic energy difference of 48.4 meV [150] . The large overestimation of the energy difference is likely related to the larger c/a distortion in MD. However, from thermodynamic considerations, the larger energy difference in MD may be needed to

PAGE 143

143 account for energy due to configurational entropy at finite temperatures. T he product of Boltzmann constant and experimental transition temperature (393 K) is k b T c = 33.8 meV. Thus, the overestimated energy difference in MD may be a necessary feature to describe phase transitions. In this work, various scalar factors of the DFT energy difference were used in the fitting. The results below are for a factor of eight times the energy difference in DFT (i.e. 8 X 5.6 = 44.8 meV) , which is close to the value used by Tinte et al . The final training database contained the ground state structures as well as many others that acted as reference points in the energy profile. For instance, the Ti atom position in the target tetragonal ph ase is (0.5, 0.5, 0.5137). Therefore, a structure with the Ti positioned at (0.5, 0.5, 0.5137 +/ database and weighted to have a gr eater energy than the reference structure. This helped guide the optimization process towards the ground state for each of the phases. Methodological Challenges The BaTiO 3 system present s numerous challenges for an accurate description in both MD and firs t principle calculations. As noted above, the energy differences between phases are likely underestimated in DFT. In addition, the optimized database did not take into account O Ba bond lengths or total volume considerations. Thus, some trial and error app roximations are needed in the fitting process. Another significant challenge is the ability to quickly and accurately find the ground state energy for each structure throughout the fitting process. This is necessary in order to accurately define the energ y differences between phases. Energy minimization requires moving the atoms such that the forces on the cores and shells are

PAGE 144

144 zero. Unfortunately, this can be a very time consuming process. Since parameter optimization usually requires tens of thousands of calls to the optimization subroutine, taking more than a few seconds per iteration makes the fitting process prohibitively slow . Thus, in this work the core position is kept fixed while the shells were allowed to r elax, which usually occurred on reasonable timescales. Static Properties The final parameter set is obtained by fitting to the optimized database using simulated annealing optimization. These values may be found in Appendix B. Tables 6 2 to 6 5 shows the lattice parameters and internal coordinate s for each of the four phases as compared to experiment and the optimized database. Of critical importance is that the c/a ratio of the tetragonal phase (Table 6 4) in MD (c/a = 1.012) satisfactorily matches experimental values (c/a = 1.011) [31] . This is a significant improvement over the Tinte et al. potential where c/a = 1.022 [150] . As noted above, the energy difference between phases is one indicator of the phase transition temperature in MD. Table 6 6 shows these energy differences in DFT and MD. The results of this work are much closer to DFT than the Tinte et al. potential [150] . In addition to the energy difference between phases, the energy profile of Ti displacement within the cubic phase is critical. In DFT calculations, whe n cubic BaTiO 3 is expanded, the off centering of Ti atoms will lower the total energy of the system. This is related to the formation of low temperature ferroelectric phases, as SrTiO 3 and CaTiO 3 do not have this feature [168] . Displacements of Ti along <100>, <110> and <111> provide information about the energy differences between cubic and the three ferroelectric phases.

PAGE 145

145 Table 6 7 shows the energy prof ile for Ti displacements in cubic BaTiO 3 for the equilibrium 0 K lattice parameter and an expanded lattice parameter. For the equilibrium lattice parameter case, the total energy increases as Ti is moved along <100>, <110> and <111>, respectively. For Tint e et al. [150] , the energy increase is approximately a factor of eight, while in this work energy the increases by less than a factor of two as compare d to DFT. For the expanded lattice parameter case, the total energy decreases as Ti is moved along <100>, <110> and <111>, respectively. For Tinte et al. [150] , the energy decrease is approximately a factor of twelve, while in this work energy the increases by about half as compared to DFT. While the current parameterization has static properties that match DFT closer than Tinte et al. [150] , it appears that this may not be beneficial. As we shall see in the next section, this underestimation of the various energy contributions is likely the culprit for the undesired MD results at finite temperatures. Considerable effort was extended to find a better parameterization, but none proved successful. MD Results at Finite Temperature The simulations are performed for atomic configurations of 8x8x8 unit cel ls under periodic boundary conditions. A Nos é Hoover [63] thermostat and barostat are applied, allowing anisotropic variations of t he lattice parameters. The time step is 0.4 fs. Configurational snapshots are recorded eve ry 20 fs. Figure 6 8 shows the lattice parameter for BaTiO 3 as a function of time at 300 K. The plot shows that the structure is oscillating with time. Typical approaches to eliminate this behavior include various equilibration schemes, changing the time i ntegration scheme ( V erlet vs. leap frog) and altering the thermostat or barostat settings. Unfortunately, no methods have been successful in showing a well defined phase at any non zero temperature.

PAGE 146

146 The symptoms in Figure 6 8 appear to be caused by the shallow energy profiles of atomic displacement. When the shape of the supercell is restricted to change in only the principle directions, a well defined tetragonal phase remains a favorable starting condition. However, this is n ot an equilibrium structure as there is external pressure on the supercell. In addition, the tetragonality is slowly decreasing and eventually will reach unity for long enough times (> 50 ps). This is unlike the Tinte et al. potential [150] , where the structure reaches equilibrium usually within the first 10 ps. Future Work on Shell Model Development The results given in this chapter show that the abilit y of the shell model to capture the correct structural properties of ABO 3 perovskites is quite limited. We are able to correctly describe the phase order and obtain reasonable structural properties, but the desired phase transitions are not achieved. This inability to capture the phase transitions is likely caused by the shallow energy profiles of atomic displacements within each phase. The small differences in bond lengths and bond angles between phases makes correcting these energy differences very diffic ult to achieve. Whether this is even possible in MD is an open question for future work.

PAGE 147

147 Table 6 1. Lattice parameters and internal coordinate positions for orthorhombic CaTiO 3 from experiment ( extrapolated to 0 K ) [169] , DFT [164] and MD . Property Experiment DFT Thi s Work a ( Å) 5.3565 5.2898 5.3663 b ( Å) 5.4316 5.4122 5.4358 c ( Å) 7.6103 7.5374 7.6213 Ca x 0.9921 0.9892 0.9935 y 0.0410 0.0480 0.0349 z ¼ ¼ ¼ Ti x 0 0 0 y ½ ½ ½ z 0 0 0 O(1) x 0.0784 0.0838 0.0853 y 0.4839 0.4776 0.470 3 z ¼ ¼ ¼ O(2) x 0.7102 0.7063 0.713 2 y 0.2878 0.2927 0.285 1 z 0.0404 0.0441 0.0442

PAGE 148

148 Table 6 2 . Experimental and MD structural properties for the rhombohedral phase of BaTiO 3 . Experiment [31] Experiment [170] Database Target This Work a ( Ã…) 4.003 4.004 4.003 4.003 (degrees) 89.84 89.87 89.84 89.92 Ba (0, 0, 0) Ti (0.5+ x Ti , 0.5+ x Ti , 0.5+ x Ti ) x Ti 0.0128 0.0111 0.0128 0.0127 O (0.5+ x O , 0.5+ x O , 0.5+ z O ) x O 0.0109 0.0110 0.0109 0.0093 z O 0.0193 0.0180 0.01 93 0.0173

PAGE 149

149 Table 6 3 . Experimental and MD structural properties for the orthorhombic phase of BaTiO 3 . Experiment [ 31] Experiment [171] Database Target This Work a ( Ã…) 3.9873 3.9900 3.9 873 3.969 9 b ( Ã…) 5.6750 5.6690 5.6750 5.671 8 c ( Ã…) 5.6900 5.6820 5.6900 5.684 2 Ba (0, 0, 0) Ti (0.5, 0, 0.5+ z Ti ) z Ti 0.0169 0.0100 0.0126 0.0140 O1 (0, 0, 0.5+ z O1 ) z O1 0.0100 0.0100 0.0136 0.0109 O2 (0.5, 0.25+ y O2 , 0.25+ z O2 ) y O2 0.0061 0.0030 0.0043 0.0048 z O2 0.0143 0.0130 0.0157 0.0145

PAGE 150

150 Table 6 4 . Experimental and MD structural properties for the tetragonal phase of BaTiO 3 . Experiment [31] Experiment [172] Database Target This Work a ( Ã…) 3.9910 3.9945 3.9910 3.9779 c ( Ã…) 4.0352 4.0335 4.0352 4.0265 c/a 1.0111 1.0098 1.0111 1.0122 Ba (0, 0, 0) Ti (0.5, 0.5, 0.5+ z Ti ) z Ti 0.0224 0.0135 0.0137 0.0155 O1 (0.5, 0.5, 0.5+ z O1 ) z O1 0.0224 0.0250 0.0240 0.0216 O2 (0.5, 0, 0.5+ z O2 ) z O2 0.0105 0.0150 0.0146 0.0113

PAGE 151

151 Table 6 5 . Experimental and MD structural properties for the cubic phase of BaTiO 3 . Experiment [173] Database Target This Work a ( Ã…) 4.000 4.000 3.9848 Ba (0, 0, 0) Ti (0.5, 0.5, 0.5) O (0.5, 0.5, 0.0)

PAGE 152

152 Table 6 6. Relative energy differences of BaTiO 3 from DFT and MD calculations. Energy differences are in units of eV/BaTiO 3 . DFT [167] MD [150] MD ( This Work ) Tetragonal Cubic 0.0056 0.0484 0.0094 Orthorhombic Cubic 0.0076 0.0569 0.0151 Rhombohedral Cubic 0.0081 0.0611 0.0187

PAGE 153

153 Table 6 7. Energy profile for Ti displacements in cubic BaTiO 3 from DFT and MD calculations. The position of the Ti atom is indicated in parenthesis. Energy differences are in units of meV. DFT (This work) MD [150] MD ( This Work ) Equilibrium a 0 <100> (0.51, 0.50, 0.50) 0.6 4.4 1.0 <110> (0.51, 0.51, 0.50) 1.2 8.2 2.2 <111> (0.51, 0.51, 0.51) 1.9 11.7 3.4 a 0 * 1.0123 <100> (0.51, 0.50, 0.50) 0.8 8.7 0.4 <110> (0.51, 0.51, 0.50) 1.4 16.9 0.8 <111> (0.51, 0.51, 0.51) 2.0 24.5 1.0

PAGE 154

154 Figure 6 1. Experimental (black) [165] and calculated (red) neutron PDFs for CaTiO 3 at 300 K.

PAGE 155

155 A B Figure 6 2. Comparison of various experimental and simulated structural properties of CaTiO 3 and Ba 0.8 Ca 0.2 TiO 3 at 300 K. A) Experimental (black) and simulated (red) neutron PDFs for Ba 0.8 Ca 0.2 TiO 3 . B) Ca EXAFS for experimental CaTiO 3 (blue), experimental Ba 0.8 Ca 0.2 TiO 3 (red) and simulated Ba 0.8 Ca 0.2 TiO 3 (green) .

PAGE 156

156 Figure 6 3. Experimental (black) [165] and calculated (red) neutron PDFs for Ba 0.8 Ca 0.2 TiO 3 .

PAGE 157

157 Figure 6 4 . PDD for Ca in simulated Ba 0.8 Ca 0.2 TiO 3 (blue) and Sr in simulated Ba 0.8 Sr 0.2 TiO 3 (red) .

PAGE 158

158 Figure 6 5 . The c/a ratio as function of Sr/Ca composition for BaTiO 3 (green), Ba 1 x Sr x TiO 3 (blue) and Ba 1 x Ca x TiO 3 ( red). For this comparison, th e temperatures are 300 K, 220 K and 280 K for BaTiO 3 , Ba 1 x Sr x TiO 3 and Ba 1 x Ca x TiO 3 , respectively. Experiment al values are denoted by triangles and calculated results denoted by squares .

PAGE 159

159 Figure 6 6 . Tetragonal < > cubic transition temperature for BaTiO 3 (green), Ba 0.8 Ca 0.2 TiO 3 (red) and Ba 0.8 Sr 0.2 TiO 3 (blue). Experimental values are denoted by triangles and calculated results denoted by squares.

PAGE 160

160 Figure 6 7. Oxygen oxygen bond distance in BaTiO 3 as a function of temperature. Reprinted with permission from G. H. Kwei, A. C. Lawson, S. J. L. Billinge, and S. W. Cheong, Journal of Physical Chemistry 97 , 2368 (1993). Copyright 1993 American Chemical Society.

PAGE 161

161 Figure 6 8 . Average BaTiO 3 lattice paramete rs at 300 K. The blue, red and green lines represent the x, y and z directions, respectively.

PAGE 162

162 CHAPTER 7 GENERAL CONCLUSIONS Summary of Project Outcomes After introductory material in Chapters 1 and 2, Chapters 3 and 4 presented an interatomic potential for the zirconium oxygen hydrogen ( Zr O H ) system using the Charge Optimized Many B ody (COMB) potential. The first step for creating a model for this system is an accurate description for Zr. It was shown in Chapter 3 that many material properties satisfactorily match experimental observations inc luding t he ground state hexagonal close packed (HCP) structure . Perhaps most importantly is the correct reproduction of the slip behavior during mechanical deformation. In Chapter 4, Zr interactions with oxygen and hydrogen are developed. The predicted gro und state structures do not match experimental observations, but the resulting potentials are adequate for the study of Zr surface corrosion. The results of O 2 and H 2 O deposition on various Zr surfaces show that the surface structure greatly affects the co rrosion behavior. The use and development of s hell m odel interatomic potentials for BaTiO 3 based ferroelectrics is presented in Chapters 5 and 6. Chapter 5 shows how MD and experimental data may be directly compared in order to validate MD results. The po tential used in Chapter 5 overestimated the distortion in the tetragonal phase of BaTiO 3 , which lead to inconclusive results about the phase transition mechanism between the tetragonal and cubic phases. The work in Chapter 6 attempts to correct this flaw, but the many challenges posed prevents a satisfactory result.

PAGE 163

163 Final T houghts on Interatomic Potential Development This dissertation covers the development of classical interatomic potentials for use in numerous applications. These applications required dif ferent interatomic potential formalisms to describe their interaction, which illustrates the difficulty of capturing the wide range of bonding environments found in materials. While a perfect description of a ll of the properties of a material system using the current models is a laudable intention, it is naïve to expect this to be possible. Even within relatively simple systems, such as metals, the energy landscapes are too complex to be fully describe d using current methods. Thus, the limitations of each model must be kept in mind when beginning the process of potential development. One of the major challenges for classical potentials is the ability to capture the energy profile for small atomic displa cements. Examples studied in this work include surfaces for metals and symmetry breaking distortions in ferroelectric ABO 3 perovskites. For each, the energy landscape found in static calculations has a direct correlation to the resulting dynamical proper ties. For the case of ferroelectrics (Chapters 5 and 6), the correct value of the energy difference between phases remains an open question as DFT calculations likely underestimate these energies. Thus, additional work correlating energy differences betwee n phases and the resulting transition temperature should be performed. One of the fundamental concepts of COMB and other empirical potentials is the concept of bond order. The bond order is a measure of bond strength and is usually applied to the attracti ve portion of the energy functions (see eq. 2 36 in Chapter 2). In [174] , a modification to the

PAGE 164

164 repulsive energy term is proposed. However, this term has historically been discarded due to being seen as an unnecessary complication. The results in this dissertation suggest that the repulsive energy between atoms does not follow an exponential increase, which COMB and many other potentials use. For example, the calculation of stacking faults in metals invol ves the shifting of atomic planes in relation to one another. This creates very short bond distances with local coordination numbers larger than the equilibrium value. Modifications to the repulsive energy components may be required for a more accurate des cription of these surfaces. A limiting factor in MD is the interaction distance imposed by the cut off function. As seen in this work, different formalisms use different methods to account for interatomic interactions. In the s hell m odel, two body intera ctions (electrostatic and non electrostatic) are usually limited to 10 Ã…. For COMB, non electrostatic two body interactions are generally limited to the first nearest neighbors, while electrostatic interactions have a n interaction distance up to 11 Ã… . The s hell m odel and COMB both limit three body interactions to the first nearest neighbors. The fact that many three body potentials in the literature have successfully described numerous systems (e.g., Si, SiO 2 , C, Zr ) despite these limitations proves that th ese assumptions are quite reasonable, at least for some systems. The three body energy contributions at longer distances may not be significant, but the energy profile within a given distance of an atom is often unknown. Thus, a full description of the three dimensional energy landscape is beneficial in understanding the resulting dynamic behavior. This would allow the implementation of additional functions that could correct for these small energy differences.

PAGE 165

165 In conclus ion, there is still much work to be done on the description of atomic interactions in materials. As computational speeds continue to increase, it may be worth the effort to include more interaction terms at the expense of increases in computational cost. I n particular, screening functions that alter the energy landscape as a function of atomic distance and coordination should be considered.

PAGE 166

166 APPENDIX A PARAMETERS DEVELOPED FOR COMB POTENTIAL Table A 1 . COMB p otential parameters for Zr. The numb er in parenthesis refers to the equation number in Chapter 2 . Parameter COMB Zr 1 [98] COMB Zr 2 [115] (e q. 2 42 ) 0.444536 1.0 m (eq. 2 43) 1.0 1.0 (eq. 2 27 ) N/A 3.30277 J (eq. 2 27 ) N/A 2.95741 K (eq. 2 27 ) N/A 0.13315 L (eq. 2 27 ) N/A 0.03161 (eq. 2 30 ) N/A 1.83154 Z (eq. 2 29 ) N/A 0.38710 Q U (eq. 2 41 ) N/A 4.0 Q L (eq. 2 41 ) N/A 2.0 D U (eq. 2 41 ) N/A 0.1040 2 D L (eq. 2 41 ) N/A 1.9998 8

PAGE 167

167 Table A 2 . Short range Zr Zr parameters for COMB Zr 1 and COMB Zr 2. The number in parenthesis refers to the equation number in Chapter 2 . Parameter COMB Zr 1 Zr Zr [98] COMB Zr 2 Zr Zr [115] A (eq. 2 38) 1122.16000 504.680745 B 1 (eq. 2 37) 711.89600 113.606465 B 2 (eq. 2 37) 171.35700 N/A B 3 (eq. 2 37) 170.50400 N/A (eq. 2 38) 1.67129 1.82140750 1 (eq. 2 37) 1.85834 1.06228155 2 (eq. 2 37) 1.29329 N/A 3 (eq. 2 37) 1.22281 N/A b 6 (eq. 2 4 4 ) 0.025674 0.57646500 b 5 (eq. 2 4 4 ) 0.025066 0.64032200 b 4 (eq. 2 4 4 ) 0.010837 0.00424300 b 3 (eq. 2 4 4 ) 0.002500 0.05968300 b 2 (eq. 2 4 4 ) 0.009815 0.43707600 b 1 (eq. 2 4 4 ) 0.002271 0.05652700 b 0 (eq. 2 4 4 ) 0.001545 0.00720200 (eq. 2 43) 2.27174 2.07148800 c 0 (eq. 2 4 5 ) N/A 0.00587989 c 1 (eq. 2 4 5 ) N/A 1.01073431 c 2 (eq. 2 4 5 ) N/A 0.00619448 c 3 (eq. 2 4 5 ) N/A 0.91015157 R Min (eq. 2 47) 4.10 3.72 R Max (eq. 2 47) 4.20 3.92

PAGE 168

168 Table A 3 . Potential parameters for Zr O and Zr H using COMB Zr 2 [115] . The number in parenthesis refers to the equation number in Chapter 2 . Parameters Zr O O Zr Zr H H Zr A (eq. 2 38) 4490.86780 4490.86780 269.781927 269.781927 B 1 (eq. 2 37) 127.347859 127.347859 154.972908 154.972908 (eq. 2 38) 3.59864133 3.59864133 2.02956767 2.02956767 1 (eq. 2 37) 1.34724491 1.34724491 1.51066566 1.51066566 b 6 (eq. 2 4 4 ) 0.56392624 27.2539415 0.34295146 0.43794103 b 5 (eq. 2 4 4 ) 0.90346893 22.9602048 0.66632834 0.95821166 b 4 (eq. 2 4 4 ) 0.61730129 12.4593517 0.10497560 0.10365350 b 3 (eq. 2 4 4 ) 1.18807908 17.5358603 0.41125598 0.70017004 b 2 (eq. 2 4 4 ) 0.70784180 5.19218715 0.06595734 0.29997625 b 1 (eq. 2 4 4 ) 0.25470056 1.65341425 0.22614598 0.16514480 b 0 (eq. 2 4 4 ) 0.13937666 1.29877933 0.20794053 0.10114535 (eq. 2 43) 2.46031875 4.69413585 3.07395511 0.00044979 c 0 (eq. 2 4 5 ) 0.02532794 0.12254375 0.00262925 0.16244931 c 1 (eq. 2 4 5 ) 0.55513242 4.63977437 0.37656389 4.35168311 c 2 (eq. 2 4 5 ) 2.26883417 0.28390107 4.99975453 0.00297161 c 3 (eq. 2 4 5 ) 0.55513242 4.63977437 0.37656389 4.35168311 R Min (eq. 2 47) 2.7 2.7 2.6 2.6 R Max (eq. 2 47) 3.1 3.1 3.0 3.0 P X (eq. 2 28 ) 2.10649777 3.05359197 0.20374809 1.00588094 P J (eq. 2 28 ) 4.59187166 4.37899001 0.34413698 0.09799903

PAGE 169

169 Table A 4 . COMB Zr 2 p otential parameters for three body Legendre polynomials [115] . The number in parenthesis refers to the equation number in Chapter 2 . Zr O Zr O Zr O Zr H Zr H Zr H LP 1 (eq. 2 48 ) N/A N/A N/A N/A LP 2 (eq. 2 48 ) 0.123860E+00 0.450477E+00 0.248917E 01 0.956763E 02 LP 3 (eq. 2 48 ) 0.141204E+00 0.144971E+00 0.515734E 01 0.503766E 01 LP 4 (eq. 2 48 ) 0.279613E+00 0.136639E+00 0.169995E 01 0.431441E 01 LP 5 (eq. 2 48 ) 0.803593E 01 0.449541E+00 0.460514E 01 0.514659E 01 LP 6 (eq. 2 48 ) 0.126040E+00 0.344180E 04 0.503008E 01 0.213845E 01

PAGE 170

170 APPENDIX B PARAMETERS DEVELOPED FOR SHELL MODEL Table B 1 . Two body p arameters for the (Ba,Ca)TiO 3 potential that is compatible with those from Ref. [150] . These parameters may be found in Ref. [73] . The number in parenthesis refers to the equation number in Chapter 2. Parameter Ca O Ba O Sr O Ti O O O A (eq. 2 19) 568.71 1061.30 167.32 3769.93 4740.0 19) 0.4805 0.3740 0.4949 0.2589 0.2686 C (eq. 2 19) N/A N/A N/A N/A 160.0 E 0 (eq. 2 20) 6.851 N/A N/A N/A N/A k (eq. 2 20) 2.443 N/A N/A N/A N/A r 0 (eq. 2 20) 1.778 N/A N/A N/A N/A

PAGE 171

171 Table B 2 . Elemental parameters for the (Ba,Ca)TiO 3 potential that is compatible with those from Ref. [150] . These parameters may be found in Ref. [73] . The number in parenthesis refers to the equation number in Chapter 2. Parameter Ca Ba Sr Ti O k 2 (eq. 2 22) 670 251.8 623.3 322.0 31.0 k 4 (eq. 2 22) N/A N/A N/A 500.0 4000.0 q Core (eq. 2 18 ) 5.62 5.62 5.62 4.76 0.91 q Shell (eq. 2 18 ) 3.76 3.76 3.76 1.58 2.59

PAGE 172

172 Table B 3 . Elemental parameters for BaTiO 3 in the new fitting attempt. The number in parenthesis refers to the equation number in Chapter 2. Parameter Ba Ti O k 2 (eq. 2 22) 603.043 674.497 123.581 k 4 (eq. 2 22) N/A 1318855.423 321946.801 q Core (eq. 2 18) 1.60 6.02 0.77 q Shell (eq. 2 18) 0.38 2.26 2.43

PAGE 173

173 Table B 4 . Three body parameters for BaTiO 3 in the new fitting attempt. The number in parenthesis refers to the equation number in Chapter 2. The cut offs between all interactions is 3.1 Parameter O Ti O O Ba O k 5 (eq. 2 21) 8 5277 1 855 . 3 1 (eq. 2 21) 0 .0676 1 43 . 46 2 (eq. 2 21) 0.3968 0.3321 0 (eq. 2 21) 90 .0 60 .0

PAGE 174

174 Table B 5 . Two body parameters for BaTiO 3 in the new fitting attempt. The number in parenthesis refers to the equation number in Chapter 2. Parameter Ba O Ti O O O A 1 (eq. 2 19) 465.4287 390.35956 235.5834 1 (eq. 2 19) 0.5117 0.3429 0.6724 A 2 (eq. 2 19) 597.8622 517.9391 29.5667 2 (eq. 2 19) 0.5558 0.347 7 0.5222 A 3 (eq. 2 19) 539.0492 1292.5525 30.0926 3 (eq. 2 19) 0.6009 0.287 1 0.6566 E 01 (eq. 2 20) 9.8970 3.6787 2.1720 k 1 (eq. 2 20) 3.9944 3.8582 2.0149 r 0 1 (eq. 2 20) 2.2383 1.5153 1.7541 E 02 (eq. 2 20) 6.1333 2.7727 4.6288 k 2 (eq. 2 20) 4.2600 6.0959 1.9688 r 0 2 (eq. 2 20) 2.3370 1.1723 2.3521 E 03 (eq. 2 20) 7.8671 1.6671 0.5945 k 3 (eq. 2 20) 1.8011 4.6117 1.1018 r 0 3 (eq. 2 20) 2.6327 1.4415 2.0593

PAGE 175

175 LIST OF REFERENCES [1] E. McCafferty, Introduction to Corrosion Science (Springer, New York, NY, 2010). [2] R. Tilley, Understanding Solids: The Science of Materials (John Wiley & Sons, West Sussex, United Kingdom, 2004). [3] R. DeHoff, Thermodynamics in Ma terials Science (CRC Press, Boca Raton, FL, 2006). [4] B. Cox, Journal of Nuclear Materials 336 , 331 (2005). [5] K. Edsinger, C. R. Stanek, and B. D. Wirth, Jom 63 , 53 (2011). [6] A. W. Cronenberg and M. S. Elgenk, Journal of Nuclear Materials 78 , 390 (197 8). [7] C. E. Ells, Journal of Nuclear Materials 28 , 129 (1968). [8] S. Yamaguch i , Journal of the Physical Society of Japan 24 , 855 (1968). [9] M. Hirabaya shi , S. Yamaguch i , T. Arai, H. Asano, and S. Hashimot o , Physica Status Solidi A Applied Research 23 , 331 (1974). [10] T. Arai and M. Hirabayashi, Journal of the Less Common Metals 44 , 291 (1976). [11] T. Tsuji and M. Amaya, Journal of Nuclear Materials 223 , 33 (1995). [12] J.S. Bradbroo k , N. Ridley, and G. W. Lorimer, Journal of Nuclear Materials 42 , 142 (1972). [13] D. R. Lide, Zirconium (CRC Press, New York, 2007 2008), Vol. 4. [14] J. C. Jamieson, Science 140 , 72 (1963). [15] C. J. Howard, R. J. Hill, and B. E. Reichert, Acta Crystallographica Section B Structural Science 44 , 116 (1988). [16] Z. Zhao, J . P. Morniroli, A. Legris, A. Ambard, Y. Khin, L. Legras, and M. Blat Yrieix, Journal of Microscopy 232 , 410 (2008). [17] M. S. Daw and M. I. Baskes, Physical Review B 29 , 6443 (1984). [18] M. I. Baskes, Physical Review B 46 , 2727 (1992). [19] R. A. Johnso n and D. J. Oh, Journal of Materials Research 4 , 1195 (1989). [20] J. G. Yu, S. B. Sinnott, and S. R. Phillpot, Physical Review B 75 , 13, 085311 (2007).

PAGE 176

176 [21] T. R. Shan, B. D. Devine, J. M. Hawkins, A. Asthagiri, S. R. Phillpot, and S. B. Sinnott, Physical Review B 82 , 9, 235302 (2010). [22] T. Liang, B. Devine, S. R. Phillpot, and S. B. Sinnott, Journal of Physical Chemistry A 116 , 7976 (2012). [23] K. M. Rabe, C. H. Ahn, and J. M. Triscone, Physics of Ferroelectrics (Springer, Berlin, 2007), Vol. 105. [24] M. De Graef and M. McHenry, Structure of Materials (Cambridge University Press, New York, NY, 2007). [25] S. O. Kasap, Principles of Electronic Materials and Devices (McGraw Hill, New York, NY, 2006). [26] C. Kittel, Introduction to Solid State Physic s (Wiley, New York, 1986). [27] R. A. Cowley, Advances in Physics 12 , 421 (1963). [28] W. Jones and N. H. March, Theoretical Solid State Physics (John Wiley & Sons, Bristol, England, 1973). [29] W. Cochran, Physical Review Letters 3 , 412 (1959). [30] B. Steele, A. D. Burns, A. Chernatynskiy, R. W. Grimes, and S. R. Phillpot, Journal of Materials Science 45 , 168 (2010). [31] G. H. Kwei, A. C. Lawson, S. J. L. Billinge, and S. W. Cheong, Journal of Physical Chemistry 97 , 2368 (1993). [32] F. W. Lytle, Jo urnal of Applied Physics 35 , 2212 (1964). [33] H. Unoki and T. Sakudo, Journal of the Physical Society of Japan 23 , 546 (1967). [34] S. A. T. Redfern, Journal of Physics Condensed Matter 8 , 8267 (1996). [35] A. M. Glazer, Acta Crystallographica Section B S tructural Science B 28 , 3384 (1972). [36] C. Menoret, J. M. Kiat, B. Dkhil, M. Dunlop, H. Dammak, and O. Hernandez, Physical Review B 65 , 224104 (2002). [37] I. Levin, V. Krayzman, and J. C. Woicik, Applied Physics Letters 102 , 162906 (2013). [38] R. D. Sh annon, Acta Crystallographica Section A 32 , 751 (1976). [39] R. Comes, M. Lambert, and A. Guinier, Solid State Communications 6 , 715 (1968).

PAGE 177

177 [40] A. S. Chaves, F. C. S. Barreto, R. A. Nogueira, and B. Zeks, Physical Review B 13 , 207 (1976). [41] B. Ravel, E. A. Stern, R. I. Vedrinskii, and V. Kraizman, Ferroelectrics 206 , 407 (1998). [42] B. Zalar, V. V. Laguta, and R. Blinc, Physical Review Letters 90 , 4, 037601 (2003). [43] E. A. Stern, Physical Review Letters 93 , 3, 037601 (2004). [44] B. Zalar, A. Lebar , J. Seliger, R. Blinc, V. V. Laguta, and M. Itoh, Physical Review B 71 , 064107 (2005). [45] R. Pirc and R. Blinc, Physical Review B 70 , 8, 134107 (2004). [46] M. Itoh, R. Wang, Y. Inaguma, T. Yamaguchi, Y. J. Shan, and T. Nakamura, Physical Review Letters 82 , 3540 (1999). [47] S. Tinte, M. G. Stachiotti, M. Sepliarsky, R. L. Migoni, and C. O. Rodriguez, Ferroelectrics 237 , 345 (2000). [48] L. Xie, Y. L. Li, R. Yu, and J. Zhu, Journal of Applied Physics 109 , 6, 054101 (2011). [49] N. L. Matsko, E. G. Maksimov, and S. V. Lepeshkin, Journal of Experimental and Theoretical Physics 115 , 309 (2012). [50] I. Ponomareva, L. Bellaiche, T. Ostapchuk, J. Hlinka, and J. Petzelt, Physical Review B 77 , 4, 012102 (2008). [51] G. Geneste, Journal of Physics Condensed Matter 23 , 16, 125901 (2011). [52] M. Bibes, Nature Materials 11 , 354 (2012). [53] M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford University Press, Oxford, 1987). [54] J. E. Jones, Proceedings of the Royal Society of London Series a Containing Papers of a Mathematical and Physical Character 106 , 463 (1924). [55] A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, Journal of Physical Chemistry A 105 , 9396 (2001) . [56] J. Tersoff, Physical Review Letters 56 , 632 (1986). [57] L. Verlet, Physical Review 159 , 98 (1967). [58] S. A. Adelman and J. D. Doll, Journal of Chemical Physics 61 , 4242 (1974).

PAGE 178

178 [59] J. D. Doll, L. E. Myers, and S. A. Adelman, Journal of Chemical Physics 63 , 4908 (1975). [60] S. A. Adelman and J. D. Doll, Journal of Chemical Physics 64 , 2375 (1976). [61] H. C. Andersen, Journal of Chemical Physics 72 , 2384 (1980). [62] H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Ha ak, Journal of Chemical Physics 81 , 3684 (1984). [63] S. Nose, Molecular Physics 52 , 255 (1984). [64] M. Parrinello and A. Rahman, Physical Review Letters 45 , 1196 (1980). [65] M. Parrinello and A. Rahman, Journal of Chemical Physics 76 , 2662 (1982). [66] R. A. Buckingham, Proceedings of the Royal Society of London Series a Mathematical and Physical Sciences 168 , 264 (1938). [67] P. M. Morse, Physical Review 34 , 57 (1929). [68] B. G. Dick and A. W. Overhauser, Physical Review 112 , 90 (1958). [69] P. J. Mitchell and D. Fincham, Journal of Physics Condensed Matter 5 , 1031 (1993). [70] P. J. D. Lindan and M. J. Gillan, Journal of Physics Condensed Matter 5 , 1019 (1993). [71] I. T. Todorov and W. Smith, (STFC Daresbury Laboratory, Cheshire, England, U nited Kingdom, 2012). [72] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, Journal of Chemical Physics 103 , 8577 (1995). [73] M. Sepliarsky, A. Asthagiri, S. R. Phillpot, M. G. Stachiotti, and R. L. Migoni, Current Opinion in Solid State & Materials Science 9 , 107 (2005). [74] G. C. Abell, Physical Review B 31 , 6184 (1985). [75] D. W. Brenner, Physical Review Letters 63 , 1022 (1989). [76] A. K. Rappe and W. A. Goddard, Journal of Physical Chemistry 95 , 3358 (1991). [77] F. H. Streitz and J. W. Mintmire, Physical Review B 50 , 11996 (1994). [78] A. Yasukawa, J SME International Journal Series A Mechanics and Material Engineering 39 , 313 (1996).

PAGE 179

179 [79] J. G. Yu, S. R. Phillpot, and S. B. Sinnott, Physical Review B 75 , 4, 233203 (2007 ). [80] B. Devine, T. R. Shan, Y. T. Cheng, A. J. H. McGaughey, M. Lee, S. R. Phillpot, and S. B. Sinnott, Physical Review B 84 , 17, 125308 (2011). [81] J. C. Slater, Physical Review 36 , 0057 (1930). [82] T. R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and S. R. Phillpot, Physical Review B 81 , 12, 125328 (2010). [83] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, Journal of Physics Condensed Matter 14 , 783 (2 002). [84] J. A. Martinez, D. E. Yilmaz, T. Liang, S. B. Sinnott, and S. R. Phillpot, Current Opinion in Solid State & Materials Science 17 , 263 (2013). [85] J. A. Nelder and R. Mead, Computer Journal 7 , 308 (1965). [86] W. L. Goffe, G. D. Ferrier, and J. Rogers, Journal of Econometrics 60 , 65 (1994). [87] J. D. Gale and A. L. Rohl, Molecular Simulation 29 , 291 (2003). [88] S. Plimpton, Journal of Computational Physics 117 , 1 (1995). [89] I. T. Todorov, W. Smith, K. Trachenko, and M. T. Dove, Journal of Materials Chemistry 16 , 1911 (2006). [90] G. Kresse and J. Furthmuller, Physical Review B 54 , 11169 (1996). [91] E. Tenckhoff, Deformation Mechanisms, Texture, and Anisotropy in Zirconium an d Zircaloy (ASTM 966, Philadelphia, PA, 1988). [92] D. J. Bacon and V. Vitek, Metallurgical and Materials Transactions A Physical Metallurgy and Materials Science 33 , 721 (2002). [93] M. A. Meyers, A. Mishra, and D. J. Benson, Progress in Materials Science 51 , 427 (2006). [94] D. Wolf, V. Yamakov, S. R. Phillpot, A. Mukherjee, and H. Gleiter, Acta Materialia 53 , 1 (2005). [95] C. Domain and A. Legris, Philosophical Magazine 85 , 569 (2005). [96] D. Vanderbilt, Physical Review B 41 , 7892 (1990). [97] J. P. Pe rdew and Y. Wang, Physical Review B 45 , 13244 (1992).

PAGE 180

180 [98] M. J. Noordhoek, T. Liang, Z. Lu, T. R. Shan, S. B. Sinnott, and S. R. Phillpot, Journal of Nuclear Materials 441 , 274 (2013). [99] M. I. Mendelev and G. J. Ackland, Philosophical Magazine Letters 87 , 349 (2007). [100] G. Simmons and H. Wang, Single Crystal Elastic Constants: A Handbook (MIT Press, Cambridge, MA, 1971). [101] N. Stojilovic, E. T. Bender, and R. D. Ramsier, Progress in Surface Science 78 , 101 (2005). [102] W. R. Tyson and W. A. Mille r, Surface Science 62 , 267 (1977). [103] Y. Udagawa, M. Yamaguchi, H. Abe, N. Sekimura, and T. Fuketa, Acta Materialia 58 , 3927 (2010). [104] A. V. Nikulina, V. A. Markelov, M. M. Peregud, V. N. Voevodin, V. L. Panchenko, and G. P. Kobylyansky, Journal of Nuclear Materials 238 , 205 (1996). [105] G. D. Samolyuk, S. I. Golubov, Y. N. Osetsky, and R. E. Stoller, Philosophical Magazine Letters 93 , 93 (2013). [106] F. Willaime, Journal of Nuclear Materials 323 , 205 (2003). [107] A. Poty, J. M. Raulot, H. Xu, J. Bai, C. Schuman, J. S. Lecomte, M. J. Philippe, and C. Esling, Journal of Applied Physics 110 , 15, 014905 (2011). [108] B. Legrand, Philosophical Magazine B Physics of Condensed Matter Statistical Mechanics Electronic Optical and Magnetic Properties 49 , 171 (1984). [109] H. Tsuzuki, P. S. Branicio, and J. P. Rino, Computer Physics Communications 177 , 518 (2007). [110] J. Li, Modelling and Simulation in Materials Science and Engineering 11 , 173 (2003). [111] D. H. Kim, F. Ebrahimi, M. V. Manuel, J. S. Tule nko, and S. R. Phillpot, Materials Science and Engineering A Structural Materials Properties Microstructure and Processing 528 , 5411 (2011). [112] V. Yamakov, D. Wolf, S. R. Phillpot, and H. Gleiter, Acta Materialia 50 , 5005 (2002). [113] R. L. Coble, Jour nal of Applied Physics 34 , 1679 (1963). [114] Z. Lu, M. J. Noordhoek, A. Chernatynskiy, S. B. Sinnott, and S. R. Phillpot (in preparation).

PAGE 181

181 [115] M. J. Noordhoek, T. Liang, T. W. Chiang, S. B. Sinnott, and S. R. Phillpot, J ournal of Nuclear Materials 452 , 285 (2014). [116] Y. Dai, J. H. Li, and B. X. Liu, Journal of Physics Condensed Matter 21 , 10, 385402 (2009). [117] W. Y. Hu, B. W. Zhang, B. Y. Huang, F. Gao, and D. J. Bacon, Journal of Physics Condensed Matter 13 , 1193 (2001). [118] D. H. Kim, M. V. Ma nuel, F. Ebrahimi, J. S. Tulenko, and S. R. Phillpot, Acta Materialia 58 , 6217 (2010). [119] P. K. Schelling, S. R. Phillpot, and D. Wolf, Journal of the American Ceramic Society 84 , 1609 (2001). [120] M. S. Khan, M. S. Islam, and D. R. Bates, Journal of M aterials Chemistry 8 , 2299 (1998). [121] X. W. Zhou, H. N. G. Wadley, J. S. Filhol, and M. N. Neurock, Physical Review B 69 , 20, 035402 (2004). [122] A. C. T. van Duin, B. V. Merinov, S. S. Jang, and W. A. Goddard, Journal of Physical Chemistry A 112 , 3133 (2008). [123] X. M. Bai, Y. F. Zhang, and M. R. Tonks, Physical Chemistry Chemical Physics 15 , 19438 (2013). [124] R. L. Gonzalez Romero, J. J. Melendez, D. Gomez Garcia, F. L. Cumbrera, and A. Dominguez Rodriguez, Solid State Ionics 219 , 1 (2012). [125] G. Fadda, L. Colombo, and G. Zanzotto, Physical Review B 79 , 13, 214102 (2009). [126] C. Domain, R. Besson, and A. Legris, Acta Materialia 50 , 3513, Pii s1359 6454(02)00173 8 (2002). [127] R. A. Ploc, Journal of Nuclear Materials 110 , 59 (1982). [128] R. A. Ploc, Journal of Nuclear Materials 113 , 75 (1983). [129] R. A. Ploc, Journal of Nuclear Materials 115 , 110 (1983). [130] P. E. Blochl, Physical Review B 50 , 17953 (1994). [131] G. Kresse and D. Joubert, Physical Review B 59 , 1758 (1999). [132] H. J. Monkhorst and J. D. Pack, Physical Review B 13 , 5188 (1976).

PAGE 182

182 [133] O. I. Malyi, P. Wu, V. V. Kulish, K. W. Bai, and Z. Chen, Solid State Ionics 212 , 117 (2012). [134] P. Zhang, B. T. Wang, C. H. He, and P. Zhang, Computational Materials Science 50 , 3297 (2 011). [135] J. F. Nye, Physical Properties of Crystals (Oxford University Press, 1985). [136] F. H. Wang, S. Y. Liu, J. X. Shang, Y. S. Zhou, Z. Y. Li, and J. L. Yang, Surface Science 602 , 2212 (2008). [137] P. Zhang, S. X. Wang, J. Zhao, C. H. He, and Y. L. Zhao, Journal of Applied Physics 113 , 013706 (2013). [138] P. A. Burr, S. T. Murphy, S. C. Lumley, M. R. Wenman, and R. W. Grimes, Corrosion Science 69 , 1 (2013). [139] M. Grosse, E. Lehmann, M. Steinbruec k, G. Kuehne, and J. Stuckert, Journal of Nuclear Materials 385 , 339 (2009). [140] G. Henkelman, B. P. Uberuaga, and H. Jonsson, Journal of Chemical Physics 113 , 9901 (2000). [141] A. Kushima and B. Yildiz, Journal of Materials Chemistry 20 , 4809 (2010). [ 142] G. Bakradze, L. P. H. Jeurgens, and E. J. Mittemeijer, Surface and Interface Analysis 42 , 588 (2010). [143] H. G. Kim, T. H. Kim, and Y. H. Jeong, Journal of Nuclear Materials 306 , 44 (2002). [144] G. Bakradze, L. P. H. Jeurgens, and E. J. Mittemeijer , Surface Science 606 , 846 (2012). [145] N. Ni et al. , Acta Materialia 60 , 7132 (2012). [146] C. Anghel, G. Hultquist, and M. Limback, Journal of Nuclear Materials 340 , 271 (2005). [147] H. M. Kandil, J. D. Greiner, and J. F. Smith, Journal of the American Ceramic Society 67 , 341 (1984). [148] J. X. Zheng, G. Ceder, T. Maxisch, W. K. Chim, and W. K. Choi, Physical Review B 75 , 104112 (2007). [149] E. Zuzek, J. P. Abriata, A. San Martin, an d F. D. Manchester, Bulletin of Alloy Phase Diagrams 11 , 385 (1990).

PAGE 183

183 [150] S. Tinte, M. G. Stachiotti, S. R. Phillpot, M. Sepliarsky, D. Wolf, and R. L. Migoni, Journal of Physics Condensed Matter 16 , 3495 (2004). [151] J. Hlinka, T. Ostapchuk, D. Nuzhnyy, J. Petzelt, P. Kuzel, C. Kadlec, P. Vanek, I. Ponomareva, and L. Bellaiche, Physical Review Letters 101 , 4, 167402 (2008). [152] Q. Zhang, T. Cagin, and W. A. Goddard, III, Proceedings of the National Academy of Sciences of the United States of America 10 3 , 14695 (2006). [153] W. Zhong, D. Vanderbilt, and K. M. Rabe, Physical Review B 52 , 6301 (1995). [154] S. Tinte, J. Iniguez, K. M. Rabe, and D. Vanderbilt, Physical Review B 67 , 064106 (2003). [155] K. M. Rabe and U. V. Waghmare, Physical Review B 52 , 13 236 (1995). [156] M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin, and Q. Hui, Journal of Physics Condensed Matter 19 , 16, 335218 (2007). [157] V. Krayzman and I. Levin, Journal of Applied Crystallography 45 , 106 (2012). [158] I. Levin, V. Krayzman, an d J. C. Woicik, Physical Review B 89 , 024106 (2014). [159] M. J. Noordhoek, V. Krayzman, A. Chernatynskiy, S. R. Phillpot, and I. Levin, Applied Physics Letters 103 , 022909 (2013). [160] V. Krayzman, I. Levin, J. C. Woicik, T. Proffen, T. A. Vanderah, and M. G. Tucker, Journal of Applied Crystallography 42 , 867 (2009). [161] P. Tangney and S. Scandolo, Journal of Chemical Physics 117 , 8898 (2002). [162] P. Tangney and S. Scandolo, Journal of Chemical Physics 119 , 9673 (2003). [163] T. Liang et al. , Annual R eview of Materials Research, Vol 43 43 , 109 (2013). [164] E. Cockayne and B. P. Burton, Physical Review B 62 , 3735 (2000). [165] I. Levin, private communication ( 2014). [166] M. Sepliarsky, Z. Wu, A. Asthagiri, and R. E. Cohen, Ferroelectrics 301 , 55 (2004). [167] A. I. Lebedev, Physics of the Solid State 51 , 362 (2009). [168] Z. X. Chen, Y. Chen, and Y. S. Jiang, Journal of Physical Chemistry B 106 , 9986 (2002). [169] X. Liu and R. C. Liebermann, Physics and Chemistry of Minerals 20 , 171 (1993).

PAGE 184

184 [170] W. Schildkamp and K. Fischer, Zeitschrift Fur Kristallographie 155 , 217 (1981). [171] G. Shirane and A. Takeda, Journal of the Physical Society of Japan 7 , 1 (1952). [172] J. Harada, T. Pedersen, and Z. Barnea, Acta Crystallographica Section A Crystal Phy sics Diffraction Theoretical and General Crystallography A 26 , 336 (1970). [173] P. Pruzan, D. Gourdain, J. C. Chervin, B. Canny, B. Couzinet, and M. Hanfland, Solid State Communications 123 , 21 (2002). [174] J. Tersoff, Physical Review B 37 , 6991 (1988).

PAGE 185

185 BIOGRAPHICAL SKETCH Mark Noordhoek was raised in McBain, MI. He earned a m aterials s cience and e ngineering (MSE) from the University of Michigan in 2009. That August , Mark joined the University of Florida (UF) for his graduate studies in MSE . In April of 2010 he joined the C omputational Materials Science Focus G roup under the guidance of Dr. Simon Phillpot. In 2012, Mark had the opportunity to be a guest researcher for six months at the National Institute of Sta ndards and Technology (NIST) in Gaithersburg, MD working closely with Drs. Igor Levin and Victor Krayzman. In August of 2014 he received his Ph.D. in MSE at UF.