Modeling Ramp Compression

151 Pages·2011·25.1 MB·English

SANDIA REPORT
SAND2011-6845
Unlimited Release
Printed December 2011

Modeling Ramp Compression
Experiments using Large-Scale
Molecular Dynamics Simulation

!"#$%&'(&)*+,-.+%/&0(&1$22*34&5(&6$%3/&0+%$2*$%&7",,38,$%/&1"9*$3:&'(&53.;$8:$"./&<$8=&>(&<83.2/&
?33.3&@(&0+%3./&)*+,$.&?(&1$22..+%/&0383,=&!(&)3,-:32+%

Prepared by
Sandia National Laboratories
Albuquerque, New Mexico 87185 and Livermore, California 94550

Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation,
a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's
National Nuclear Security Administration under contract DE-AC04-94AL85000.

Approved for public release; further dissemination unlimited.
Issued by Sandia National Laboratories, operated for the United States Department of Energy
by Sandia Corporation.

NOTICE: This report was prepared as an account of work sponsored by an agency of the
United States Government. Neither the United States Government, nor any agency thereof,
nor any of their employees, nor any of their contractors, subcontractors, or their employees,
make any warranty, express or implied, or assume any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product, or process
disclosed, or represent that its use would not infringe privately owned rights. Reference herein
to any specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise, does not necessarily constitute or imply its endorsement,
recommendation, or favoring by the United States Government, any agency thereof, or any of
their contractors or subcontractors. The views and opinions expressed herein do not
necessarily state or reflect those of the United States Government, any agency thereof, or any
of their contractors.

Printed in the United States of America. This report has been reproduced directly from the best
available copy.

Available to DOE and DOE contractors from
U.S. Department of Energy
Office of Scientific and Technical Information
P.O. Box 62
Oak Ridge, TN 37831

Telephone: [phone]
Facsimile: [phone]
E-Mail: [email protected]
Online ordering: http://www.osti.gov/bridge

Available to the public from
U.S. Department of Commerce
National Technical Information Service
5285 Port Royal Rd.
Springfield, VA 22161

Telephone: [phone]
Facsimile: [phone]
E-Mail: [email protected]
Online order: http://www.ntis.gov/help/ordermethods.asp?loc=7-4-0#online

2
SAND2011–6845
UnlimitedRelease
PrintedDecember2011
Modeling Ramp Compression Experiments using
Large-Scale Molecular Dynamics Simulation
AidanP.Thompson,J.MatthewD.Lane,JonathanZimmerman,
MichaelP.Desjarlais,GaryS.Grest,ReeseE.Jones,ThomasR.Mattsson,JeremyA.Templeton
SandiaNationalLaboratories
P.O.Box5800
Albuquerque,NM 87185-1322
[email protected]
ToddDitmire,HernanQuevedo
CenterforHighEnergyDensityScienceTheUniversityofTexasatAustin
MichaelI.Baskes
UniversityofCalifornia,SanDiego
YogendraM.Gupta,J.MichaelWiney
InstituteforShockPhysicsWashingtonStateUniversity
Abstract
Moleculardynamicssimulation(MD)isaninvaluabletoolforstudyingproblemssensitivetoatom-
scale physics such as structural transitions, discontinuous interfaces, non-equilibrium dynamics, and
elastic-plasticdeformation.Inordertoapplythismethodtomodelingoframp-compressionexperiments,
several challenges must be overcome: accuracy of interatomic potentials, length- and time-scales, and
extractionofcontinuumquantities. Wehavecompleteda3yearLDRDprojectwiththegoalofdevelop-
ingmoleculardynamicssimulationcapabilitiesformodelingtheresponseofmaterialstorampcompres-
sion. Thetechniqueswehavedevelopedfallintothreecategories(i)moleculardynamicsmethods(ii)
interatomicpotentials(iii)calculationofcontinuumvariables. Highlightsincludethedevelopmentofan
accurateinteratomicpotentialdescribingshock-meltingofBeryllium,ascalingtechniqueformodeling
slow ramp compression experiments using fast ramp MD simulations, and a technique for extracting
plasticstrainfromMDsimulations. AllofthesemethodshavebeenimplementedinSandia’sLAMMPS
MDcode,ensuringtheirwidespreadavailabilitytodynamicmaterialsresearchatSandiaandelsewhere.
3
Contents
I ShockCompression 7
1 Dynamiccompressionofpolymers 8
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
DensityFunctionalTheorywithAM05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Classicalreactiveinteractionpotentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Classicalnon-reactiveinteractionpotentials . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Constructionofthesimulationcells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 Modelingshockresponse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.4 Resultsindensepolymer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
ShockHugoniot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Chemicalstructure/dissociation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.5 Resultsinpolymerfoam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
ShockHugoniot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
Shocktemperatureandvoidcollapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2 Dynamiccompressionofsiliconandgermanium 28
2.1 Introductionandmethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2 Resultsanddiscussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3 Laser-drivendynamicmeltingofsilicon 34
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2 Experimentalsetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.5 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4 Dynamiccompressionofaluminum 47
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 ComputationalMethods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
MDSimulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
ContinuumVariablesfromMDSimulations . . . . . . . . . . . . . . . . . . . . . . . . . . 48
NonlinearElasticContinuumCalculations . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4
Elasticshockwaveprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
[100]Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
[111]Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
[110]Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
TemperatureCalculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5 ConcludingRemarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
II RampCompression 76
5 DynamicCompressionofBeryllium 77
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 PhysicalProperties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.3 PhaseDiagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4 DynamicCompression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6 Time-PositionScalingrelationinramp-loading 88
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.2 Methodofcharacteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.3 Methodofscaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
Kinematicanddynamicsimilarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Invarianceinthewaveequation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Rampscalinginthenonlinearelasticresponseregime . . . . . . . . . . . . . . . . . . . . . 100
Rampscalingintheplasticresponseregime . . . . . . . . . . . . . . . . . . . . . . . . . . 100
III Methods 106
7 ContinuumProperties 107
7.1 EulerianAtomistic-ContinuumFormulationbyHardy . . . . . . . . . . . . . . . . . . . . . 107
BalanceLaws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
DensitiesandLocalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
EnergyandForceAssumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
ContinuumExpressions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
7.2 LagrangianAtomistic-ContinuumFormulation . . . . . . . . . . . . . . . . . . . . . . . . 111
BalanceLaws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
DensitiesandLocalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
EnergyandForceAssumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
ContinuumExpressions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5
7.3 EvaluationandComparisonofSpatialandMaterialFrameExpressions . . . . . . . . . . . 113
Stressforaconstrainedfinitetemperaturesystem . . . . . . . . . . . . . . . . . . . . . . . 113
Finitetemperaturedeformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Tensilestretchingofacenter-crackedbody . . . . . . . . . . . . . . . . . . . . . . . . . . 116
7.4 Cauchy-BornStressandFreeEnergyCalculation . . . . . . . . . . . . . . . . . . . . . . . 122
Quasi-harmonicCauchy-Bornmodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Quasi-harmonicmodelforpairwisepotentials . . . . . . . . . . . . . . . . . . . . . . . . . 125
Quasi-harmonicmodelforEAMpotentials . . . . . . . . . . . . . . . . . . . . . . . . . . 125
ExamplesoffreeenergyandCauchy-Bornstresscalculation . . . . . . . . . . . . . . . . . 127
7.5 PlasticDeformationGradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
IncrementalApproach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
TotalStrainApproach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
Examplecalculation: uniaxialstretching . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Appendix: Thermodynamicintegration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
8 NewMethodsinLAMMPS 139
8.1 TheLAMMPSSHOCKPackage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
Hugoniostat(fix nphug) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
MultiscaleSimulationTechnique(fix msst) . . . . . . . . . . . . . . . . . . . . . . . 142
DirectRamp-loading(fix wall/piston) . . . . . . . . . . . . . . . . . . . . . . . . . 145
Materialadditionandremoval(fix append atoms) . . . . . . . . . . . . . . . . . . . 146
8.2 ExtensionstoATCpackage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
9 Communications 154
9.1 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
9.2 Presentations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
10 Acknowledgements 156
6
Part I
Shock Compression
1 Dynamic compression of polymers
PrincipalAuthors: ThomasR.MattssonandJ.MatthewD.Lane
Thissectioncontainsapaperentitled“First-PrinciplesandClassicalMolecularDynamicsSimulationof
Shock Polymers ” that was published in Phys. Rev. B in 2010. It was authored by Thomas R. Mattsson, J.
MatthewD.Lane,Kyle. R.Cochrane,MichaelP.Desjarlais,AidanP.Thompson,FlintPierce,andGaryS.
Grest. [SAND2009-3918J]
Abstract
Organic polymers and nanocomposites are increasingly being subjected to extreme environments.
Molecular-scalemodelingofthesematerialsoffersinsightintofailuremechanismsandresponse. Den-
sity functional theory (DFT) molecular dynamics (MD) and classical MD simulations of the principal
shock Hugoniot are presented for two hydrocarbon polymers, polyethylene (PE) and poly(4-methyl-
1-pentene) (PMP). DFT results are in excellent agreement with experiment data, which is currently
available up to 80 GPa. Further, we predict the PE and PMP Hugoniots up to 350 GPa and 200 GPa,
respectively. Forcomparison, westudiedtworeactiveandtwonon-reactiveinteractionpotentials. For
the latter, the exp-6 interaction of Borodin et al. showed much better agreement with experiment than
OPLS. For the reactive force fields, ReaxFF displayed decidedly better agreement than AIREBO. For
shocks above 50 GPa, only the DFT results are of high fidelity, establishing DFT as a reliable method
for shocked macromolecular systems. We extend these results to include low-density polymer foams
using NEMD techniques. We find good quantitative agreement with both experiment and hydrocode
simulations. Further,wehavemeasuredlocaltemperaturestoinvestigatetheformationofhotspotsand
polymerdissociationnearfoamvoids.
1.1 Introduction
Over the last few years, first-principles simulations in combination with increasingly accurate shock ex-
periments at multi-Mbar pressure have yielded important insights into how matter behaves under extreme
conditions. While comprehensive advances have been made for many light elements, for example deu-
terium[2,3]andcarbon[4],progresshasbeenslowerforequallyimportant,albeitmorechallenging,mate-
rialslikemolecularcrystalsandpolymers[5].
Modeling a macro-molecular material requires trade-offs in system size and fidelity of the atomic inter-
action. Itisnotclearapriorihowtobeststrikethebalancebetweenthecompetingrequirements,sincethe
criticalvariabledeterminingsystemresponsemayinsomecasesbethesize-specificstructure/geometry,or
in other cases the bond reactivity and interaction fidelity. In order to resolve this important question, we
havesimulatedshockcompressionoftwodifferentpolymersusingfirst-principlesDensityFunctionalThe-
ory[6](DFT)andclassicalmoleculardynamics(MD)simulationsusingfourdifferentinteractionpotentials:
ReaxFF[7],AIREBO[8],OPLS[9],andBorodinetal.’sexp-6potentials[10].
Polymerfoamisacrucialpartofthedynamichohlraumplatformusedtogenerateintensex-rayradiation
atSandia’sZ-machine[11–13]anddesignsforInertialConfinedFusion(ICF)ofteninvolvefoams. Asex-
perimentaldesignsarebeingrefined–primarilyasaresultofaninterplaybetweenradiation-hydrodynamics
7
modeling,improveddiagnostics,andimproveddrivers–theimportanceofhavingaccuratematerialmodels
based on a fundamental understanding of the material’s behavior is becoming increasingly clear. Model-
ing both the qualitative and quantitative aspects of shock propagation in polymers and polymer foams is
thereforeofsignificantimportance.
Figure 1: Simulation snapshots of the initial foam structure (top) and of the propagating shock front (bot-
tom). Theshockshownisstrongenoughtodissociatetheatomsofthepolymer,whichcanbeseenvaporized
withinthevoidspaceaheadofthefront. Thesnapshotsareshadedtoindicatedepth. Blackspacesareareas
inwhichvoidmergersallowunobstructedviewsthroughtheentiresample.
For generality, we chose to study two polymers, polyethylene (PE) and poly(4-methyl-1-pentene) (PMP
or TPX). PE is an extensively used general-purpose plastic with the simplest possible linear alkane struc-
ture. Atactic PMP is commonly used in shock studies as a low-density polymer foam and has specialized
applicationsintargetmaterialsforinertiallyconfinedfusion(ICF)studies. PMPisabranchedalkanewitha
bulkysidechain. Thetwopolymerswereselectedtoberepresentativeoftwodifferentclassesofpolymers,
as PMP is amorphous at room temperature, while PE is semi-crystalline. We believe the results are likely
to be widely applicable to other macromolecular materials. Several experimental studies of shocked PE
havebeenreported[14–17]inadditiontothequantitativedataforbothPEandPMPfromtheLASLshock
handbook[18].
In expanding our work to foams, we note the inherent difficulty presented by the multiple length scales
responsible for the shock response in foams. In addition to the atomic length scale response which we
have shown DFT and MD can capture well, foams introduce the larger scale of the void structure, wall
thickness, and even longer range density inhomogeneity within samples. The likely need to model these
larger scales drove us to shift our modeling from DFT/MD to MD/Mesoscale. Thus, although DFT was
notablysuccessfulinmodelingdensepolymers,wedidnotapplythesetechniquestofoamsystems. Onthe
mesoscalewereportbrieflyoncontinuumhydrocoderesultsforcomparisonwithMD.Completecontinuum
analysis[19]andexperimentalresults[20]willbereportedelsewhere.
1.2 Methodology
WhileDFTisacomputationallycostlymethodinwhichitisnecessarytoreducethenumberofsimulation
atomstoseveralhundredatmost,itallowsforahigh-fidelitydescriptionofchemicalbondsandinter-atomic
8
repulsion. Reactiveforce-fieldslikeReaxFF[7]andAIREBO[8],whichhavebeenappliedtostudyshocks
inhydrocarbons[21–23],cannotaccuratelycapturetherangeofresponsescomparedtoDFT,howeverthey
allowforchemicalreactionsandsignificantlylargersystemsizes. Non-reactiveforce-fields,liketheOPLS
potential of Jorgensen et al. [9] and the exp-6 potential of Borodin et al. [10], are computationally much
more efficient and allow for even larger system sizes, but do not allow covalent bonds to break or form,
which can become important for strong shocks. The relative speed of the different methods is naturally of
interest,OPLSandexp-6arethefastestpotentials,AIREBOisapproximatelytwiceasslowwhileReaxFF
is30timesslowerthanOPLSandexp-6. AllclassicalMDsimulationswereruninLAMMPS[1,24]while
theDFT-MDsimulationswereperformedwithVASP5.1.40[25,26].
1.2.1 DensityFunctionalTheorywithAM05
DFT is a formally exact representation of the Schro¨dinger equation. However, in practice, the choice
of exchange-correlation functional determines the accuracy. We employed the recently developed multi-
purpose Armiento-Mattsson (AM05) functional [27], it is a functional with no empirically determined pa-
rameters. Itimprovesuponthelocaldensityapproximation(LDA)byreproducingtwomodelsystemswith
known solutions: the uniform electron gas and the surface jellium [27,28]. AM05 has demonstrated high
fidelity for many solids [28,29]; in reference [28], the performance of seven functionals was compared
for twenty representative semiconductors, simple metals, transition metals, alkali-halides, and oxides. For
AM05,thepredictionbias(meanaverageerror)inlatticeconstantissmallwhileLDAandPBEbothexhibit
significantbias. Onaverage,AM05isbetterthanchoosingbetweenLDAandPBEforeachsolidseparately
(TableIofreference[28])anddoesaswellasthedecidedlymorecomputationallydemandinghybridfunc-
tionals studied in Refs. [30] and [31]. Furthermore, AM05 also works well for hydrogen bonding in the
waterdimer[32]andforchemicalreactionenergiesforalargenumberofmolecularreactions[33].
Both polyethylene and poly(4-methyl-1-pentene) are materials where van der Waals forces are impor-
tant. AlthoughpreviousworkusingDFTbyByrdandRicedemonstratedadifficultyinmodelingenergetic
molecularsolidsatlowpressure,theyfoundanincreasingaccuracyastheexternalpressureincreases[34],
andthebehaviorunderstrongshocksisdominatedbythehighpressureresponse. Furthermore,thelackof
van der Waals attraction in AM05 [29] makes it decidedly different from most other exchange-correlation
functionals. Sincethefunctionaldisplaysamonotonicbehavioruponexpansionandcompression[29],itis
arguablysuitableforstudyingcompressioninvanderWaalssystems. TheapplicabilityofAM05toshocked
energeticmaterialswillbethesubjectoffuturework[35].
Finally,AM05wasrecentlythebestsuitedfunctionaltomodelquartz[36]upto1TPainthedevelopment
of a high-pressure shock impedance standard. Taken together, there are ample reasons not only to employ
AM05,buttoexpecthigh-fidelityresultsforshockcompression.
The DFT-MD simulations were performed with VASP 5.1.40 [25,26] using stringent convergence set-
tings [3,37,38]. The plane wave cutoff was above 800 eV in order to converge the stress-tensor [3,39].
Pulayerrorsareminimalduetothehighcutoffandtheappreciablechangesinvolumeforthedifferentcal-
culationsalongtheHugoniot[40]. Theionictimestepisbetween0.1and0.5fsdependingontemperature.
Steady-state simulations in the NVT ensemble used a Nose´-Hoover thermostat with a time-constant of 80
timesteps. Velocitieswerescaledtocontroltemperatureintheramped-temperaturesimulations. Partitionof
kinetic energy between hydrogen and carbon was verified by monitoring the temperature for each element
separately [41]. Complex k-point sampling with mean-value point (1,1,1) was used due to its high accu-
4 4 4
racy for disordered structures at high temperature. Electronic states were occupied according to Mermin’s
9
finite-temperatureformulationofDFT[42],afactorthatisparticularlyimportantinthewarm-densematter
regimewhereelectronsthermalexcitationsaresignificant.
1.2.2 Classicalreactiveinteractionpotentials
ReaxFF and AIREBO are reactive potentials which allow the possibility of dynamic bond formation and
breaking. ReaxFF uses bond-order and charge equilibration to model local chemical changes. It has been
used to simulate a wide variety of materials and processes, including molecular solids under shock and
detonation [23,43,44]. AIREBO is based on Brenner’s REBO potential augmented with explicit 6-12
dispersiontermsandhaspreviouslybeenusedtomodelshockpropagationinshortchainhydrocarbons[22,
45].
BoththeAIREBOandReaxFFcalculationswereperformedusingthestandardLAMMPSparallelimple-
mentations [1], which have been validated against the original serial codes. For AIREBO, a 10.2 A˚ cutoff
wasused. ForReaxFF,a10 A˚ cutoffwasusedandthechargeequilibrationconvergencetolerancewasset
to10 6. ThespecificformoftheReaxFFpotentialinLAMMPSisdescribedinChenowethetal.[7]. The
−
ReaxFFparametervaluesarethesameasthoseusedbyStrachanetal.[43].
1.2.3 Classicalnon-reactiveinteractionpotentials
The OPLS and exp-6 potentials have pre-assigned non-breakable bonds. In these potentials, interactions
within and between molecules are described by a set of atomic potentials which include van der Waals,
electrostatic,molecularbond,angle,andtorsioninteractionterms.
U = U +U +U +U (1)
tot nonbond bond ang tor
where U is the sum of the van der Waals and Coulomb potentials. The bond potential and angle
nonbond
potential are harmonic. The dihedral potentials for the OPLS and exp-6 force fields differ only slightly in
form.
Animportantdifferencebetweenthesepotentialsisintheformofthenonbondedinteractions. TheOPLS
nonbondedpotentialiscomposedofstandard12-6Lennard-Jones(LJ)andCoulombpotentials[9,46],while
theexp-6forcefieldofBorodinetal.[47,48]utilizesaBuckinghamexponential-6formforthenonbonded
pairpotential. Nonbondedinteractionsarecalculatedbetweenallatompairswithindifferentmoleculesand
betweendistantatomswithinthesamemolecule. Therespectivefunctionalformsare
12 6
σ σ q q
UOPLS = 4￿ ij ij +k i j (2)
nonbond ij r − r coul r
￿￿ ij￿ ￿ ij￿ ￿ ij
and
C
exp-6 ij
U = A exp( B r ) (3)
nonbond ij − ij ij − r6
ij
where ￿ sets the energy scale, σ sets the separation scale for the ij pair, q is partial charge. A is the
ij ij ij
strength of the potential’s repulsive component, Bi−j1 is the characteristic decay length, and Cij indicates
the strength and range of the attractive component. The atoms in the exp-6 potential are uncharged. Both
potentials were cut off at 12A˚ for the large simulation cells. An 8A˚ cutoff was used to equilibrate the
10