Anisotropic frictional response and corresponding heating in cyclotrimethylene-trinitramine molecular crystals are studied using molecular dynamics simulations. The nature of damage and temperature rise due to frictional forces is monitored along different sliding directions on the primary slip plane, (010), and on non-slip planes, (100) and (001). Correlations between the friction coefficient, deformation, and frictional heating are established. We find that the friction coefficients on slip planes are smaller than those on non-slip planes. In response to sliding on a slip plane, the crystal deforms easily via dislocation generation and shows less heating. On non-slip planes, due to the inability of the crystal to deform via dislocation generation, a large damage zone is formed just below the contact area, accompanied by the change in the molecular ring conformation from chair to boat/half-boat. This in turn leads to a large temperature rise below the contact area.

Frictional sliding of surfaces is a fundamental process that governs a wide variety of material properties.1–5 An archetypal example is frictional heating of energetic molecular crystals.6,7 While the frictional heating is believed to play an essential role in the initiation of chemical reactions and creation of hot spots in these materials, molecular mechanisms of frictional heating remain largely unknown.

Cyclotrimethylene trinitramine (RDX), a widely used high energy density material (HEDM), is sensitive to thermal and mechanical insults.8–13 A polymer binder is often used to desensitize RDX. Potential initiation mechanisms in RDX include anisotropic plasticity and fracture,14,15 resulting in the decomposition of molecules.16 Compressive and shear dynamics models have been developed, based on reactive10 and nonreactive17 force fields, to study anisotropic shock sensitivity of RDX.18 Both experiments and theoretical models indicate larger sensitivity and more chemical reactivity normal to (100) and (210) planes.

The anisotropic shock response in RDX has been explained by steric hindrance model.19–21 For shock sensitive directions, shock impact triggers deformation on slip planes which has higher steric hindrance for shear deformation. Shock loading in RDX on (100) shock plane triggers deformation on (110) ⟨110⟩ slip plane, which has higher steric hindrance compared to (120) shock plane which triggers deformation on (010) ⟨100⟩ slip plane of low steric hindrance.10 Hence, (100) and (210) requires large mechanical work for deformation which results into more chemical decomposition and large temperature rise in these systems. Similar anisotropic behavior in observed in Pentaerythritol tetranitrate (PETN) where [100] direction has less steric hindrance for shear deformation compared to [110] direction.22 

Cawkwell et al. suggested that viscous flow may cause intense heating within the shear bands and promote high rates of thermal molecular decomposition.23 They proposed that stacking faults enclosed by partial dislocation loops can give rise to plastic hardening and thermal molecular decomposition in RDX.24 

A long-standing problem in high energy density materials (HEDM) is how microstructures affect molecular decomposition processes under external stimuli.23,25 When shock is applied to HEDM, certain microstructures concentrate the impact energy spatially, creating localized hot spots. Under certain conditions, some of these hot spots may ignite and grow. For example, energetic crystals with defects such as voids, grain boundaries, and cracks are known to be more sensitive to shock ignition than a perfect crystal. Mechanisms underlying these hot-spot nucleation and growth have attracted a great deal of attention.26–29 Thus, hot spots are critical factors in the impact initiation of HEDM, and understanding, predicting, and controlling the behavior of hot spots is the key to optimizing the performance of HEDM and making them safer and more resistant to accidents.

Among the potential hot-spot forming defects, frictional heating and associated plastic deformation at interfaces have attracted much attention.1,2 When a delaminated interface is present in HEDM, friction at that interface produces intense heating. Understanding friction at the atomistic level is of great importance for the safe handling of HEDMs.30 

To study friction and anisotropic plastic response across different slip planes and slip directions in α-RDX crystal,31 we perform non-equilibrium molecular dynamics (MD) simulations.32 Experimental and computational studies indicate that (010) is the primary and (011) and (021) are secondary slip planes in α-RDX and that [001] and [100] are the slip directions in the (010) plane.33 Figure 1(a) shows the RDX crystal structure and its primary slip plane. We study the effect of frictional sliding on damage and localized heating along different sliding directions on the primary slip plane, (010), as well as on non-slip planes, (100) and (001). We observe that slip systems are less susceptible to heat generation because they deform easily. Non-slip systems show larger increases in temperature than slip systems, in accordance with the steric hindrance model of heat generation and dissipation.

FIG. 1.

(a) Crystal structure of α-RDX showing its primary slip plane (magenta), where cyan, grey, red, and blue colors represent carbon, hydrogen, oxygen, and nitrogen atoms, respectively. (b) Schematic of the friction simulation. (c) and (d) Calculated kinetic friction coefficients for different planes and sliding directions as a function of time. The surfaces and sliding directions are denoted by parentheses and square brackets, respectively.

FIG. 1.

(a) Crystal structure of α-RDX showing its primary slip plane (magenta), where cyan, grey, red, and blue colors represent carbon, hydrogen, oxygen, and nitrogen atoms, respectively. (b) Schematic of the friction simulation. (c) and (d) Calculated kinetic friction coefficients for different planes and sliding directions as a function of time. The surfaces and sliding directions are denoted by parentheses and square brackets, respectively.

Close modal

Figure 1(b) shows the simulation setup, which consists of two blocks of RDX crystals. The upper block moves at a constant speed of 100 m/s across a surface of the lower block. The sliding direction is the x-axis, and we monitor frictional effects in the x-z plane where the z-axis is normal to the interface between the two blocks. Periodic boundary conditions are applied along the x and y axes. The dimensions of the lower and upper blocks are 200 Å × 120 Å × 100 and 100 Å × 30 Å × 100 Å, respectively, in the x, y, and z directions. The system contains a total of about 280,000 atoms. After creating the system, we minimize the energy using the conjugate gradient method and thermalize the system at a temperature of 300 K for 400 ps. During frictional sliding, a 10 Å-thick layer at the bottom of the lower block is kept fixed. The total distance traveled by the top block is 100 Å. The simulation time steps during equilibration and frictional sliding are 1 fs and 0.3 fs, respectively. The simulations are performed using LAMMPS package.34 

We use the interatomic potential developed by Smith and Bharadwaj (SB),35 which consists of bonding, non-bonding, and Coulombic interactions. Bonding interaction includes bond stretching and bending. Non-bonding interaction is descried by Lennard-Jones and Buckingham potentials with parameters taken from Wallis and Thompson36 and Bedrov et al.,37 respectively. The Coulombic interaction is represented by fixed atomic charges, and the long-range electrostatic interactions are calculated using the Particle-Particle Particle-Mesh (PPPM) method.38 The functional form of the SB potential is given in the supplementary material. The crystal density, coefficient of linear expansion, and elastic constants predicted by the SB potential agree well with experimental results. This interatomic potential has been used extensively to study shock response and mechanical properties of RDX.23,24,39

The friction coefficient (μ) is the ratio of the total normal force (FN) to the force (FS) in the sliding plane of the lower block of RDX [Fig. 1(b)]

FS=μFN.
(1)

The friction coefficient thus defined changes with time. Figures 1(c) and 1(d) show the calculated values of μ on different planes (denoted by parentheses) and sliding directions (denoted by square brackets) as a function of time. Initially, the values are high because of the bonding between the two blocks of RDX. After 80 ps, they reach constant values. We estimate μ by averaging between 80 and 100 ps, and these average values are presented in Table I. Overall, the estimated values of friction coefficient for non-slip planes are higher than those for slip planes.

TABLE I.

Time-averaged friction coefficient μ for slip and non-slip planes.

Surface/directionμ
Slip plane (010)/[001] 0.54 
Slip plane (010)/[100] 0.49 
Non-slip (100)/[010] 0.67 
Non-slip (100)/[001] 0.59 
Non-slip (001)/[100] 0.55 
Non-slip (001)/[010] 0.64 
Surface/directionμ
Slip plane (010)/[001] 0.54 
Slip plane (010)/[100] 0.49 
Non-slip (100)/[010] 0.67 
Non-slip (100)/[001] 0.59 
Non-slip (001)/[100] 0.55 
Non-slip (001)/[010] 0.64 

We observe that different kinds of intermolecular and intramolecular defects are generated in RDX crystal. Intermolecular defects are identified by finding the number of nearest neighbors for each RDX molecule using Voronoi analysis.40 An RDX molecule has 14 nearest neighbors in a defect-free crystal. During frictional sliding, any RDX molecule whose coordination number is not 14 is considered a defect. Through this analysis, we can identify dislocations and other complex defects.

RDX molecules in ambient condition are in chair conformation. Intramolecular defects are generated when a chair conformation changes to a half-boat or boat conformation, even though the total number of nearest neighbors of an RDX molecule remains 14. Changes in ring conformation of RDX molecules are determined using ring puckering analysis.41 

Figures 2(a) and 2(b) show defect generation due to frictional sliding along the [001] direction in the (010) plane. As we move the upper block, partial dislocation loops are generated in the lower block at 10° of (001) plane just below the contact area. These partial dislocations loops grow on the (001) plane along [010] direction and dissipate energy supplied by the moving upper block. We observe slight pileup of RDX molecules but hardly any changes in ring conformations of RDX molecules or in the formation of a damage zone. Supplementary material video, S1.mp4, shows the deformation in the system during frictional sliding on (010) plane along [001] direction.39 

FIG. 2.

(a) and (b) Deformation in α-RDX crystal due to the frictional sliding on the (010) slip plane along the [001] direction. Red and brown spheres are RDX molecules whose coordination number is not 14, whereas the rings are RDX molecules whose conformation changes from chair to half-chair/boat. Red arrow indicates the sliding direction. Plane on which inter-molecular forms is at 10° of (001) which is shown in magenta in (a). (c) and (d) Time evolution of defects. The density of intermolecular defects (c) and % change of intramolecular defects (d) are plotted as a function of time during frictional sliding in the (010) plane along the [100] and [001] directions.

FIG. 2.

(a) and (b) Deformation in α-RDX crystal due to the frictional sliding on the (010) slip plane along the [001] direction. Red and brown spheres are RDX molecules whose coordination number is not 14, whereas the rings are RDX molecules whose conformation changes from chair to half-chair/boat. Red arrow indicates the sliding direction. Plane on which inter-molecular forms is at 10° of (001) which is shown in magenta in (a). (c) and (d) Time evolution of defects. The density of intermolecular defects (c) and % change of intramolecular defects (d) are plotted as a function of time during frictional sliding in the (010) plane along the [100] and [001] directions.

Close modal

When frictional sliding is along the [100] direction in the (010) plane, we observe changes in conformations of molecular rings and a small damage zone below the contact area between the two blocks but no dislocation lines.

Figures 2(c) and 2(d) show the densities of intermolecular and intramolecular defects for sliding along the [100] and [001] directions in the (010) slip plane. We find that changes in ring conformations for the [100] and [001] sliding directions are similar and correlated with the concentration of intramolecular defects in the crystal, and they both have similar intramolecular defect densities and damage zones. Lower damage in the RDX crystal on a slip plane is due to the facts: first, both (010) [001] and (010) [100] planes has low steric hindrance for shear deformation, and second, the friction coefficient is smallest along the [100] direction (μ = 0.49), whereas large deformation in the system is observed along the [001] direction (μ = 0.54) due to the generation of partial dislocation loops in the system. During frictional sliding on [001] direction on (010) plane, we observe large dislocation activity around 20 ps. Initially, partial dislocation loops are nucleated at the interface of the upper and lower blocks. These partial dislocation loops grow along the [01¯0] direction. Also, partial dislocation loops are nucleated from the bottom of the lower block; however, they grow along the [010] direction towards the upper block. These two sets of dislocation start annihilating each other around 20 ps and the dislocation density drops. This cycle of nucleation and annihilation of dislocation loops happens again between 40 and 60 ps and we observe multiple peaks of defect densities as shown in Fig. 2(c). This mechanism of cyclic nucleation and annihilation of partial dislocation loops is also observed in a larger system consisting of 3.16 million atoms as shown in Figs. S3 and S4, in supplementary material, movie S4.mp4, and is described in the supplementary material.

For non-slip (100) and (001) planes, we observe large damage zones just below the contact area of the two blocks, as shown in Figs. 3(a) and 3(b). The size of the damage zone increases as the upper block moves across the surface of the lower block. By the time a steady sliding state is reached, 4–5 layers of RDX molecules are significantly deformed. Deformation occurs mainly due to changes in the ring conformation of the RDX molecules from the chair to boat or half-boat conformations and also due to changes in the coordination numbers of RDX molecules. Just below these 4–5 layers, we observe a layer of intermolecular defects [Fig. 3(b)] and material pileup in front of the upper block. It can be seen from Fig. 3(b) that all molecules in the pileup region have boat conformation. Figures 3(c) and 3(d) show defect densities on the (100) and (001) planes. Supplementary material video, S2.mp4, shows the formation of a damage zone during frictional sliding on (100) plane along [010] direction.39 

FIG. 3.

(a) and (b) Deformation in α-RDX crystal due to frictional sliding on (100) non-slip plane along [010] sliding direction. Red and brown spheres are RDX molecules whose coordination number is not 14, whereas the rings are RDX molecules whose conformation is changed from chair to half-chair/boat. Red arrow indicates the sliding direction. Plane on which inter-molecular forms is (100) for (100)[010] frictional sliding which is shown in magenta in (a). (c) and (d) Time evolution of defects. The density of intermolecular defects (c) and % change of intramolecular defects (d) are plotted as a function of time during frictional sliding in the (100) and (001) planes.

FIG. 3.

(a) and (b) Deformation in α-RDX crystal due to frictional sliding on (100) non-slip plane along [010] sliding direction. Red and brown spheres are RDX molecules whose coordination number is not 14, whereas the rings are RDX molecules whose conformation is changed from chair to half-chair/boat. Red arrow indicates the sliding direction. Plane on which inter-molecular forms is (100) for (100)[010] frictional sliding which is shown in magenta in (a). (c) and (d) Time evolution of defects. The density of intermolecular defects (c) and % change of intramolecular defects (d) are plotted as a function of time during frictional sliding in the (100) and (001) planes.

Close modal

Sliding on the non-slip (001) plane along the [010] direction is a special case. Even though the (001) [010] system has a high friction coefficient (μ = 0.64), frictional sliding on (001) plane along [010] direction activates deformation by slip in the material on (010) slip plane along [001] direction. The mechanism of deformation for this case is explained in supplementary material (Figs. S1 and S2).39 

We also examine how the local temperature is correlated with the friction coefficient and damage in the system. Figure 4 shows how the volume fraction of the system with local temperature greater than 400 K increases with time during frictional sliding. Supplementary material video, S3.mp4, shows the region of the system with temperature greater than 400 K for slip and non-slip planes.39 

FIG. 4.

Volume % of the system where the local temperature is greater than 400 K.

FIG. 4.

Volume % of the system where the local temperature is greater than 400 K.

Close modal

These temperature rises for each frictional sliding direction depend on its friction coefficient, nature of deformation, and its steric hindrance for shear deformation. Let us first examine the primary slip plane, (010). The friction coefficient is high for the [001] direction and low for the [100] direction. But the temperature rise is similar in both cases. Note that dislocations are generated in response to frictional sliding only in the [001] direction and that the friction coefficient is higher along this direction than in other directions of the (010) plane. These dislocation lines help dissipate energy supplied to the system during frictional sliding. Therefore, the increase in temperature is insignificant even though (010) [001] has a high friction coefficient. Also, both (010) [001] and (010) [100] have low steric hindrance. Hence, it requires less mechanical energy for deformation and low temperature rise.

For non-slip planes, where friction coefficients are higher and the planes do not deform easily via slip, we observe massive heating due to the formation of large damage zones and also due to changes in ring conformation. The increase in temperature correlates well with the density of intramolecular defects. Detailed analysis correlating the work done on the system due to frictional sliding and rise in temperature are shown in the supplementary material (Figs. S5 and S6).39 

As shown in Figs. 2(d) and 3(d) and Figs. S7(a) and S7(b), defect density generated during frictional sliding on (010)[100] slip plane is comparable to that of non-slip plane (100) [001]; however, (010) [100] slip plane shows less heating. An et al. has shown that the strain energy density required for shear deformation during shock loading on (100) is 1.67, which is 25% higher than that on (120) and its associated slip plane (010) [100].10 Similarly, we observe that larger work done is required during frictional sliding on (100) [001] non-slip plane as compared to that on the (010)[100] slip plane as shown in Fig. S8 and Table S1. This happens because steric hindrance for shear deformation on (010) [100] is lower compared to that on the non-slip plane (100).10 Directions with higher steric hindrance have large overlap between molecules. Hence, they require large mechanical work for deformation which eventually results in high temperature along those directions.

To locate the region of the highest temperature during frictional heating in these systems, contour plots of temperature on slip and non-slip systems are calculated as a function of displacement of the upper block as shown in Fig. 5 and Fig. S9, respectively.39 It can be seen from these plots that the location of the highest temperature is always at the interface between the upper and lower blocks, but its location at the interface changes as the upper block moves. Initially, the maximum temperature rise is in front of the upper block as shown in Figs. 5(a) and 5(b), and the maximum temperature is around 350 K near the leading edge of the sliding block. However, later its position shifts to the center of the upper block and the maximum temperature becomes around 500 K, as shown in Fig. 5(d). We observe that heat is generated at the leading edge of the sliding block as well as at the interface of the upper and lower blocks. The generated heat is confined at the interface and dissipates at the leading and tailing edges of the upper block, which may explain the formation of the maximum temperature rise region at the center of the interface.

FIG. 5.

Contour plots of temperature for (010) [100] frictional heating several values of displacement, d, of the upper block. Red arrows indicate the sliding direction.

FIG. 5.

Contour plots of temperature for (010) [100] frictional heating several values of displacement, d, of the upper block. Red arrows indicate the sliding direction.

Close modal

In summary, our simulation shows that deformations in the slip planes are caused by the motion of dislocations, while changes in molecular conformations give rise to deformations in the non-slip planes. Temperature rise during frictional sliding on (010) slip planes is lower compared to the non-slip planes, which agrees with the steric hindrance model for shear deformation. The simulation also shows that the friction coefficient tends to be lower on slip planes. This atomistic understanding of anisotropic frictional responses may provide insight toward the ultimate goal of predicting and controlling the behavior of hot spots.

See supplementary material for the description of Smith and Bharadwaj Potential, computation of the work done on the systems during frictional sliding, effect of steric hindrance during frictional sliding, and temperature contour plot on non-slip plane (001) [100] during frictional sliding.

S1.mp4 and S2.mp4 show frictional sliding on slip and non-slip plans, respectively. S3.mp4 shows temperature rise in the systems during frictional sliding, and S4.mp4 shows mechanism of dislocation nucleation and annihilation on (010) [001] slip plan during frictional sliding.

This work was supported by the Air Force Office of Scientific Research Grant No. FA9550-16-1-0042. Simulations were performed at the Center for High Performance Computing of the University of Southern California. We thank Professor Dana Dlott for valuable advice and encouragement.

1.
Y. F.
Mo
,
K. T.
Turner
, and
I.
Szlufarska
,
Nature
457
(
7233
),
1116
(
2009
).
2.
S.
Vinod
,
C. S.
Tiwary
,
L. D.
Machado
,
S.
Ozden
,
J.
Cho
,
P.
Shaw
,
R.
Vajtai
,
D. S.
Galvao
, and
P. M.
Ajayan
,
Nano Lett.
16
(
2
),
1127
(
2016
).
3.
T.
Onodera
,
Y.
Morita
,
A.
Suzuki
,
M.
Koyama
,
H.
Tsuboi
,
N.
Hatakeyama
,
A.
Endou
,
H.
Takaba
,
M.
Kubo
,
F.
Dassenoy
,
C.
Minfray
,
L.
Joly-Pottuz
,
J.
Martin
, and
A.
Miyamoto
,
J. Phys. Chem. B
113
(
52
),
16526
(
2009
).
4.
E. H.
Cook
,
M. J.
Buehler
, and
Z. S.
Spakovszky
,
J. Mech. Phys. Solids
61
(
2
),
652
(
2013
).
5.
G. S.
Smith
,
N. A.
Modine
,
U. V.
Waghmare
, and
E.
Kaxiras
,
J. Comput.-Aided Mater. Des.
5
(
1
),
61
(
1998
).
6.
W. P.
King
,
S.
Saxena
,
B. A.
Nelson
,
B. L.
Weeks
, and
R.
Pitchimani
,
Nano Lett.
6
(
9
),
2145
(
2006
).
7.
S.
You
,
M. W.
Chen
,
D. D.
Dlott
, and
K. S.
Suslick
,
Nat. Commun.
6
,
6581
(
2015
).
8.
K.
Nomura
,
R. K.
Kalia
,
A.
Nakano
,
P.
Vashishta
,
A. C. T.
van Duin
, and
W. A.
Goddard
,
Phys. Rev. Lett.
99
(
14
),
148303
(
2007
).
9.
D.
Bedrov
,
J. B.
Hooper
,
G. D.
Smith
, and
T. D.
Sewell
,
J. Chem. Phys.
131
(
3
),
034712
(
2009
).
10.
Q.
An
,
Y.
Liu
,
S. V.
Zybin
,
H.
Kim
, and
W. A.
Goddard
,
J. Phys. Chem. C
116
(
18
),
10198
(
2012
).
11.
Y.
Li
,
R. K.
Kalia
,
A.
Nakano
,
K.
Nomura
, and
P.
Vashishta
,
App. Phys. Lett.
105
(
20
),
204103
(
2014
).
12.
A.
Strachan
,
E. M.
Kober
,
A. C. T.
Van Duin
,
J.
Oxgaard
, and
W. A.
Goddard
,
J. Chem. Phys.
122
(
5
),
054502
(
2005
).
13.
L.
Zhang
,
S. V.
Zybin
,
A. C. T.
Van Duin
, and
W. A.
Goddard
,
J. Energ. Mater.
28
,
92
(
2010
).
14.
J.
Sharma
,
R. W.
Armstrong
,
W. L.
Elban
,
C. S.
Coffey
, and
H. W.
Sandusky
,
Appl. Phys. Lett.
78
(
4
),
457
(
2001
).
15.
Y. C.
Chen
,
K.
Nomura
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
,
Appl. Phys. Lett.
93
(
17
),
171908
(
2008
).
16.
S. M.
Walley
,
J. E.
Field
, and
M. W.
Greenaway
,
Mater. Sci. Technol.
22
(
4
),
402
(
2013
).
17.
J. D.
Clayton
and
R.
Becker
,
J. Appl. Phys.
111
(
6
),
063512
(
2012
).
18.
R. A.
Yetter
,
F. L.
Dryer
,
M. T.
Allen
, and
J. L.
Gatto
,
J. Propul. Power
11
(
4
),
683
(
1995
).
19.
X. B.
Jiang
,
S.
Guo
,
M.
Yao
, and
J. H.
Peng
,
Combust., Explos. Shock Waves
50
(
1
),
118
(
2014
).
20.
M. C.
Gwak
,
T. Y.
Jung
, and
J. J. I.
Yoh
,
J. Mech. Sci. Technol.
23
(
7
),
1779
(
2009
).
21.
J. K. A.
Amuzu
,
B. J.
Briscoe
, and
M. M.
Chaudhri
,
J. Phys. D: Appl. Phys.
9
(
1
),
133
(
1976
).
22.
J. J.
Dick
,
R. N.
Mulford
,
W. J.
Spencer
,
D. R.
Pettit
,
E.
Garcia
, and
D. C.
Shaw
,
J. Appl. Phys.
70
(
7
),
3572
(
1991
).
23.
M. J.
Cawkwell
,
T. D.
Sewell
,
L. Q.
Zheng
, and
D. L.
Thompson
,
Phys. Rev. B
78
(
1
),
014107
(
2008
).
24.
M. J.
Cawkwell
,
K. J.
Ramos
,
D. E.
Hooks
, and
T. D.
Sewell
,
J. Appl. Phys.
107
(
6
),
063512
(
2010
).
25.
A.
Prakash
,
A. V.
McCormick
, and
M. R.
Zachariah
,
Nano Lett.
5
(
7
),
1357
(
2005
).
26.
A.
Tokmakoff
,
M. D.
Fayer
, and
D. D.
Dlott
,
J. Phys. Chem.
97
(
9
),
1901
(
1993
).
27.
J. W.
Mintmire
,
D. H.
Robertson
, and
C. T.
White
,
Phys. Rev. B
49
(
21
),
14859
(
1994
).
28.
B. L.
Holian
,
T. C.
Germann
,
J. B.
Maillet
, and
C. T.
White
,
Phys. Rev. Lett.
89
(
28
),
285501
(
2002
).
29.
K.
Kadau
,
T. C.
Germann
,
P. S.
Lomdahl
, and
B. L.
Holian
,
Science
296
(
5573
),
1681
(
2002
).
30.
S. P.
Sivapirakasam
,
M.
Surianarayanan
, and
G.
Swaminathan
,
J. Loss Prev. Process Ind.
22
(
2
),
254
(
2009
).
31.
C. S.
Choi
and
E.
Prince
,
Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater.
28
(
9
),
2857
(
1972
).
32.
C. J.
Mundy
,
S.
Balasubramanian
,
K.
Bagchi
,
M. E.
Tuckerman
,
G. J.
Martyna
, and
M. L.
Klein
,
Reviews in Computational Chemistry
(
John Wiley & Sons, Inc.
,
2007
), pp.
291
.
33.
N.
Mathew
and
R. C.
Picu
,
Chem. Phys. Lett.
582
,
78
(
2013
).
34.
S.
Plimpton
,
J. Comput. Phys.
117
(
1
),
1
(
1995
).
35.
G. D.
Smith
and
R. K.
Bharadwaj
,
J. Phys. Chem. B
103
(
18
),
3570
(
1999
).
36.
E. P.
Wallis
and
D. L.
Thompson
,
J. Chem. Phys.
99
(
4
),
2661
(
1993
).
37.
D.
Bedrov
,
A.
Chakravarthy
,
G. D.
Smith
,
T. D.
Sewell
,
R.
Menikoff
, and
J.
Zaug
,
J. Comput.-Aided Mater. Des.
8
,
77
(
2001
).
38.
S.
Plimpton
,
R.
Pollock
, and
M.
Stevens
, in
Proceedings of the Eighth SIAM Conference on Parallel Processing for Scientific Computing
(
1997
).
39.
K.
Joshi
,
M.
Losada
, and
S.
Chaudhuri
,
J. Phys. Chem. A
120
(
4
),
477
(
2016
).
40.
C. H.
Rycroft
,
Chaos
19
(
4
),
041111
(
2009
).
41.
D.
Cremer
and
J. A.
Pople
,
J. Am. Chem. Soc.
97
(
6
),
1354
(
1975
).

Supplementary Material