Transcript A Facility for Simulating the Dynamic Response of Materials

MD and Force Field Development for HMX
• Level 0 - Generic Force Field (Dreiding) calculations
– Density of States
– Pressure Loading
– Phase transitions
• Level 1 - Vibrationally accurate force field for DMN,
HMX & RDX
–
–
–
–
–
DFT (B3LYP) calculations on isolated monomers
QUEST calculations on condensed phase systems
FFOPT parameterization
Intra-, inter-molecular VET, phonon - phonon couplings
H-bond effects
Crystallographic Forms of HMX
E = 0.20 kg/cm2
E = 0.10 kg/cm2
r = 1.58
r = 1.78
r = 1.894
429 K - to
melting point
r = 1.839
E = 0.20 kg/cm2
Impact Energy
E = 0.75 kg/cm2
Stable @ 300K
377 - 429 K
Correlation of Density of States (from MD) with Sensitivity (h50
measurement)
b-HMX h50 = 0.33m
TATB h50 = 3.2m
a-HMX h50 = N/A
g-HMX h50 = 0.14m
d-HMX h50 = N/A
most sensitive
Chair Form - a, Boat Form - b,g,d
Short intramolecular
H:::O contacts,
responsible for
VET and
energy
localization in
N-N bond?
QM [B3LYP/6-31G**] E = 0.0
N-N
N-C
N-O
N-NO2(inv)
N-C-N-C
Experimental Crystal
Gas Phase DFT/6-31G**
1.354/1.373
1.392/1.401
1.437-1.472
1.442-1.478
1.205-1.233
1.224-1.231
8.0/20.3
7.5/18.6
18.1/43.3/101.7/117.1
19.4/43.1/99.3/113.5
QM [B3LYP/6-31G**] E = 2.218 Kcal/mol
N-N
N-C
N-O
N-NO2(inv)
N-C-N-C
Experimental Crystal
Gas Phase DFT/6-31G**
1.353/1.366
1.406/1.417
1.444-1.471
1.450/1.458
1.214-1.237
1.223/1.224
3.4/12.2
5.0/20.6
64.7/68.2/100.7/105.0
72.2/72.3/93.4/93.6
Comparison of Molecular & Cell Parameters for HMX b and a
forms
HMX-b
N-N
N-C
N-O
O-N-O
O-N-N
N-N-C
C-N-C
CH:::O-N-O
Experiment
1.354(l)/1.373(s)
1.437/1.448/1.455/1.472
1.205/1.209/1.221/1.233
125.9/126.7
115.9/116.2/117.4/117.9
115.1/117.4/118.2
122.4(l)/123.8(s)
2.192/2.337
Cell Params
a (A)
b (A)
c (A)
a
b
g
Volume (A**3)
s (gm/cc)
Expt.
6.54
11.05
8.70
90.0
124.3
90.0
519.4
1.894
HMX-a
Dreiding Exp-6
1.329(s)/1.332(l)
1.441/1.443
1.204/1.206
115.2/116.1
122.0/123.0
121.5/122.4
115.1(s)/116.4(l)
2.343/2.432
HMX-b
Dreiding Exp-6
6.603
11.094
8.94
90.00
125.55
90.00
532.799
1.846
% Error
-0.96
-0.40
-2.76
-1.01
-2.58
2.53
Experiment
1.353(s)/1.366(l)
1.444/1.448/1.451/1.471
1.214/1.223/1.236/1.237
124.8/124.9
116.4/117.2/118.0/118.8
115.9/119.2/119.3/119.6
119.6/123.7
2.117/2.229
Expt
15.14
23.89
5.91
90.0
90.0
90.0
2138.7
1.839
Dreiding Exp-6
1.31
1.441/1.444
1.2
115.6
122.2
119.3/119.6/119.8/121.1
120.1/121.9
2.245/2.272
HMX-a
Dreiding Exp-6
14.448
24.383
5.809
90.0
90.0
90.0
2046.46
1.922
% Error
4.57
-2.06
1.76
4.31
-4.51
Dreiding Frequencies Comparable to DFT Calculation
Comparison of Frequencies for HMX beta form
3500
Frequency (cm-1)
3000
Dreiding
B3LYP/6-31G**
2500
2000
1500
1000
500
0
0
10
20
30
40
50
Mode Number
60
70
80
Comparison of different van der Waals functions for H:::O
interaction
•Dreiding Exponential6 has the softest inner
wall
Energy comparison Exp-6 vs Lennard-Jones
Energy Exp-6
55
Energy LJ 12-6
Energy LJ 9-6
E (kcal/mol)
45
Energy COMPASS
•Dreiding LJ 12-6 is too
steep on the inner wall
35
25
15
5
-5 1.5
1.7
1.9
2.1
2.3
2.5
2.7
2.9
•COMPASS force field
(9-6) is a good
approximation to the
exp-6 form
R (A)
Conclusion:
•Accuracy - Exponential - 6
•Speed - COMPASS
•Dreiding 9-6 is not
adequate
Comparison of Crystallographic Cell Data with Experimental
Values from Cady & Ollinger for b-HMX
Cell Length vs Pressure Dreiding Exp-6
12
Cell Length (A)
11
10
A(expt)
A(exp-6)
9
B (expt)
B(exp-6)
8
C(expt)
C(exp-6)
7
6
0.00
1.00
2.00
3.00
4.00
5.00
6.00
Pressure (GPa)
7.00
8.00
9.00
10.00
Comparison of b-HMX Elastic Constants and Bulk Modulus with
Experimental Data (Joe Zaug, LANL)
HMX beta Elastic Constants & Bulk Modulus
30
Experiment
25
Exp-6
GPa
20
15
10
5
0
-5
C11 C12 C13 C15 C22 C23 C25 C33 C35 C44 C46 C55 C66
Constants
K
HMX Cold Compression Curves for the 4 crystal morphologies
HMX EOS Dreiding Exp-6
beta
alpha
gamma
delta
100
90
Pressure (GPa)
80
a
70
• d form is the most
compressible
60
• stable b form is
intermediate in
compressibility
50
40
d
30
b
20
10
0
0.5
0.6
• a form is the least
compressible followed
closely by the g form
0.7
0.8
V/V_0
0.9
1
Isothermal P-V curves for b-HMX
EOS HMX beta Dreidxing Exp-6 Temp
0K
300K
600K
900K
1200K
100
90
Pressure (GPa)
80
70
60
50
• minimization for 0K
• 20ps NPT
Molecular Dynamics
at elevated
temperatures
• evidence of melting
above 600K?
40
30
20
Melting behavior?
10
0
250
• P-V curves obtained
from
350
450
550
Volume (A**3)
650
750
850
Calculation of Shock Adiabat intersection with P-V isotherms
HMX beta EOS Dreiding Exp-6
0K
100
300K
90
600K
Tang's Adiabat
Pressure (GPa)
80
70
60
50
40
30
20
10
0
0.45
0.55
0.65
V/V_0 0.75
0.85
0.95
• Proof-ofprinciple
exercise to
calculate
temperatures
from
intersection
points
• Volume needs
to be converted
to engineering
units
b-HMX Cv from Phonon Dispersion Curves of Crystal
b -HMX Specific Heat
310
C_v (cal/mol-K)
260
111
222
210
333
160
110
60
10
100
600
1100
1600
Temperature (K)
2100
2600
• Series
converged at
222 directions
in Brillouin
zone
• 90% of
asymptotic
high T limit
reached at
1400K
Convergence of Gruneisen Parameter from 50ps (0.5fs step) MD
of 4x3x2 supercell of b-HMX
Gruneisen Param eter Convergence
0.42
Gruneisen
0.41
0.4
0.39
Gruneisen
0.38
0.37
0.36
10
20
30
40
Tim e Interval (ps)
• Gruneisen Parameter shows converged behavior by
40ps of MD
50
Level 1 - Vibrationally Accurate Force Field Development Dimethylnitramine
H
O
125.9
H
H
C
1.238
N
1.456
N
H
125.4
O
H
C
H
C
O
1.239
126.0
N
1.361
H
1.357
1.458
125.8
N
O
H
Parameter QM (B3LYP/6-31G**) Gas Phase Expt. (ED)
N-N
1.385
1.382
N-O
1.232
1.223
N-C
1.459
1.46
ONO
125.9
130.4
ONN
117
114.8
CNC
120.3
116.2
CNN
115.6
127.6
NCH(ip)
107.4
NC2 inversion
28.6
0.0
NO2 inversion
1.7
0.0
ONNC dihedral
16.9,164.9
0.0/180.0
H:::O
2.386
H
H
C
H
H
Crystal (X-ray)
1.331
1.237/1.242
1.455/1.449
124.0
117.8/118.2
124.6
117.6/117.7
0.0
0.0
0.0/180.0
2.586
Comparison of Vibrational Frequencies for DMN
Experiment B3LYP/6-31G** Scale factor
3033
3171
0.956
2948
3061
0.963
1462
1532
0.954
(1441)
1504
0.958
1304
1353
0.964
1248
1271
0.982
1023
1039
0.985
838
853
0.982
626
630
0.994
427
408
1.047
Mode
s'(CH3)
s(CH3)
ds'(CH3)
ds(CH3)
s(NO2)
(NN)
r||(CH3)
s(CNC)
d(NO2)
d(CNC)
B2
s'(CH3)
s(CH3)
a(NO2)
ds'(CH3)
ds(CH3)
a(CNC)
r||(CH3)
r(NO2)
d(NN)
3033
2948
1528
(1454)
1411
1292
1110
619
350
3159
3054
1634
1518
1448
1356
1061
630
325
0.960
0.965
0.935
0.958
0.974
0.953
1.046
0.983
1.077
3024
2913
1747
1620
1687
1267
1061
485
330
3011
2947
1600
1444
1495
1227
1035
521
376
22
1
-72
10
-84
65
75
98
-26
A2
a(CH3)
da(CH3)
r_|_(CH3)
(NO2)
(CH3)
2993
(1456)
(1144)
162
107
3120
1478
1151
168
114
0.959
0.985
0.994
0.964
0.939
3027
1659
1042
156
131
3012
1421
1021
66
129
-19
35
123
96
-22
B1
a(CH3)
da(CH3)
r_|_(CH3)
(NO2)
g(NN)
(CH3)
2993
1450
1050
762
225
(115)
3122
1487
1147
771
-117
162
0.959
0.975
0.915
0.988
3028
1665
1133
724
397
117
153.78
3013
1429
1091
758
327
99
54.1
-20
21
-41
4
-102
16
0.710
Rms. error
Dreiding
3031
3061
1649
1686
1717
1016
1249
764
539
427
Error
This work
3017
16
2949
-1
1463
-1
1536
-95
1327
-23
1251
-3
997
26
780
58
615
11
462
-35
Sym
A1
Geometric Parameters for DMN Crystal
Parameter
Crystal
This work % error
Dreiding
1.40
a
6.088
6.003
6.296
-6.77
b
6.266
6.690
6.919
-1.40
c
5.995
6.079
5.987
0.00
a
90.0
90.0
90.0
b
-2.12
114.28
116.7
117.17
g
0.00
90.0
90.0
90.0
-2.78
N-N
1.331
1.368
1.320
N-O 1.237/1.242 1.230/1.232
1.201/1.202
N-C 1.455/1.449 1.451/1.457
1.445/1.449
ONO
124.0
125.0
ONN 117.8/118.2
117.9
CNC
124.6
127.0
CNN 117.6/117.7 116.2/116.9
114.3
122.4/123.3
119.8
119.7/120.5
NC2 inversion
NO2 inversion
0.0
0.0
0.0
0.0
0.0
0.0
ONNC dihedral
0.0/180.0
0/180.0
0.0/180.0
H:::O
H:::H
2.586
2.226
2.572
2.390
2.709
2.277
%error
-3.42
-10.42
0.13
0.00
-2.53
0.00
0.83
TATB
Overview
•
•
TATB (1,3,5-triamino2,4,6-trinitrobenzene)
has planar structure.
This makes it easy to
pack and can have high
density. Experimental
density at STP is
1.9374 g/cc.
TATB crystal has low
symmetry: triclinic (P-1).
a (A)
Experiment
Exp6-Dreiding
error
b (A)
c (A)
density
(g/cc)
9.01 9.028 6.812 1.9374
8.986 8.976 6.883 1.8852
0.27% 0.58% 1.04% 2.69%
TATB
Force Field Result
Compare with
Experiment
Exp6-Dreiding Force
Field uses Morse
bond stretch and
Exp6 van der Waals
interaction. The 300K
isothermal curve fits
well with the
experimental data.
50
Olinger Experiment
Exp6-Dreiding
40
Pressure (GPa)
•
300K Isotherms
Low Pressure Region
7
6
5
4
30
3
2
1
20
0
0.41
0.43
0.45
0.47
0.49
0.51
10
0
0.3
0.35
0.4
0.45
Specific Volume (cc/g)
0.5
0.55
0.53
TATB
Isothermal Curves with Exp6-Dreiding Force Field
Exp6-Dreiding Isotherms
cold
300K
600K
900K
50
40
low Pressure Region
10
Pressure (GPa)
8
30
6
4
20
2
0
0.39
10
0.44
0.49
0
0.3
0.4
0.5
Specific Volume (cc/g)
0.6
0.54
0.59
0.64
TATB
Force Field Improvement
and Ab initio calculations
•
Lacking experimental data,
we use ab initio calculations
to improve and validate our
force field. Optimization of
the initial results is on-going.
Specific Heat at constant
pressure for gas phase TATB
are calculated at different
temperatures from
vibrational frequencies.
140
120
Specific Heat (cal/mol.K)
•
Variation of Specific Heat C p with Temperature
for gas phase TATB from QM
100
80
60
40
20
0
0
500
1000
1500
Temperature (K)
2000
2500
TATB
Dimer Binding Energy
Ab initio calculations of
Dimer Binding Energy
•
Dimer binding Energy of
two TATB molecules as a
function of the separation
distance is calculated to
explore H-bond potential.
At STP, the H2N-NO2
separation distance in
TATB crystal is 3.400A in a
direction, 3.421A in b
direction.
(3.4,-1.618)
separation
distance in
TATB
crystal
0
Energy (Kcal/mol)
•
0.5
-0.5
-1
-1.5
-2
-2.5
3.1
3.3
3.5
3.7
3.9
H2N-NO2 separation distance
4.1
4.3
TATB
Ab initio Calculations of Dihedral Angle Torsion Energy
Energy Variation with O-N-C-C dihedral angle
1

12






30
8.7
9.5
5.2
60
14.6
17.3
6.3
90
9.2
8.9
2.3
F0 is the O-N-C-C dihedral angle.
F1 is the left hand H-N-C-C dihedral angle.
F2 is the right hand H-N-C-C dihedral angle.
F3 is the non-planarity of C6 ring.
Energy (Kcal/mol)
0
1


9
6
3
0
0
30
60
O-N-C-C dihedral angle
90
TATB
Future Work
•
•
•
•
Shock Hugoniot from isotherms
Gruneisen Parameter
Further Force Field Improvement
Large scale calculations
Kel-F800
Overview
•
•
•
•
•
Kel-F 800 is a random copolymer of
chlorotrifluoroethylene and vinylidene
fluoride monomer units in a 3:1 ratio.
The presence of the vinylidene fluoride
disrupts the the crystallinity of the
chlorotrifluroethylene to form an essentially
amorphous polymer
Although amorphous, the polymer is very
dense due to the presence of the Cl and F
atoms
It is used in composites and as a binder for
many plastic-bonded explosive systems
First atomistic/molecular study of Kel-F 800
system.
Cl
F
H
F
C
C
C
C
F
F
H
F
3
Kel-F 800
1
Kel-F 800
The packing dilemma
•
Using 2 chains causes
alignment within unit cell giving
a crystalline type appearance.
•
Using more chains in the unit
cell overcomes this problem.
•
The Cerius2 Amorphous
builder initially builds to the
correct density but minimizes
to a much lower density than
given from experimental.
•
The MSC developed code for
Cohesive Energy Density
packs the molecules in such a
way as to maintain the correct
density.
24 monomers - 2 chains
24 monomers - 16 chains
20 Å
50 Å
The appearance of chain alignment is apparent when only 2 chains
are used however the relative complexity of the 16 chain case
should eliminate this problem.
Kel-F800
75
Validation
•
•
|CED| (cal/cm^3)
Due to the lack of experimental data for
the pure Kel-F 800 polymer system,
poly(chlorotrifluorethlene-co-vinylidene
fluoride) Some validation work was
done by calculating Cohesive Energy
Densities and Solubility parameters
using a MSC “in-house” developed
code.
Initial studies and choice of force field
were conducted on pure PCTFE,
poly(chlorotrifluoroethylene) for which
some experimental data is given.
The Dreiding-EXP6 force field
appears to be the force field of
choice. It is, however, somewhat
slower than the COMPASS force field.
Cohesive Energy Density
PCTFE
70
65
Upper limit of Experiment
60
Dreiding-EXP6
55
lower limit of Experiment
50
45
40
2
5
16
No of "Polymer" chains in cell
Cohesive Energy Density
Kel-F 800
100
|CED| (cal/cm^3)
•
COMPASS
80
60
Dreiding-EXP6
40
COMPASS
20
0
2
5
No of "Polymer"chains in cell
16
Kel-F 800
Kel-F800
Force Field Choice
•
•
•
Initial work was done using the
Dreiding force field. This uses
Lennard-Jones (LJ) 12-6 potential to
calculate the van der Waals
interactions.
This forcefield gives a very steep
inner wall slope for the pair potential
between 2 non-bonded atoms.
The Buckingham “EXP6” potential
gives a much a more gentle inner
wall slope, however is
computationally more demanding
and substantially slower.
The Compass force field from MSI
uses LJ 9-6 potential and is
supposedly optimized for polymer
simulations. It is also faster than the
EXP6 potential.
70
COMP
60
EXP6
LJ9-6
50
LJ12-6
Dreiding-EXP6
Pressure (GPa)
•
Cold Compression Curves
Force field comparisons
40
30
20
10
0
0.45
COMPASS
0.55
0.65
0.75
V/Vo
0.85
0.95
Kel-F800
Kel-F 800
Isothermal Compression Curves
Dreiding-EXP6 force field
Isothermal Compression
•
•
•
The Dreiding-EXP6 and
Compass force fields have
proven to be the best.
Compass has the advantage of
being computationally faster
than EXP6.
The disadvantage is that
Compass is not paramatized
sufficiently to handle the HE
materials.
Dreiding-EXP6, although slow,
will be able to handle the
inclusion of the HMX and TATB
molecules for a more complete
shock wave simulation on an
atomistic level.
70
0K
100K
60
200K
300K
50
Pressure (GPa)
•
40
30
20
Cold Compression
10
0
0.25
0.3
0.35
0.4
0.45
Volume (cm^3/g)
0.5
0.55
0.6
Kel-F800
Future work
•
•
•
•
Determining molecular weight dependence of chains
used in cell
• finding the compromise between accuracy
and speed.
Calculating the GRUNEISEN parameterand other
physical properties
• as a function of temperature and pressure
• longer Molecular Dynamics runs
Repeating for various polymer binders
• eg Estane
Huge simulations combining polymer binder and HE
materials.
HO CH2 O
4
O
O
C
CH2 C
4
O
O CH2 O C
4
ESTANE
H
N
CH2
N C O