Limited-Sample Coarse-Grained Strategy and Its Applications to

Article pubs.acs.org/JPCC

Limited-Sample Coarse-Grained Strategy and Its Applications to Molecular Crystals: Elastic Property Prediction and Nanoindentation Simulations of 1,3,5-Trinitro-1,3,5-triazinane Jian Liu,*,† Qun Zeng,† Yalin Zhang,‡ and Chaoyang Zhang*,† †

Institute of Chemical Materials, China Academy of Engineering Physics (CAEP), P.O. Box 919-311, Mianyang, Sichuan 621900, China ‡ Institute of Computer Application, China Academy of Engineering Physics (CAEP), P.O. Box 919-1201, Mianyang, Sichuan 621900, China S Supporting Information *

ABSTRACT: Modeling plastic deformation of crystalline materials by all-atomistic methods remains a challenge, and large-scale methods, such as coarse-grained (CG) methods, are highly desirable. To overcome the difficulty in constructing CG potentials for stiff molecular crystals by conventional approaches, we propose a limited-sample coarse-grained (LSCG) strategy. We construct a CG potential of α-1,3,5-trinitro-1,3,5-triazinane (α-RDX), a widely used energetic material, and perform coarse-grained molecular dynamics (CGMD) simulations to validate the LSCG potential. We compare the calculated mechanical properties with other reported results. The results show that the LSCG method is effective when compared with the all-atomic methods and provides adequate insight into systems with larger scales. Therefore, through the LSCG method, the deformation mechanisms of α-RDX crystals under nanoindentation conditions are revealed by a series of CGMD simulations that resemble nanoindentation on their (100) surface, with nanoholes sited differently. Valuable results are obtained and understood. That is, the stress around the nanohole can trigger void collapses when the nanohole is located at a shallow position directly beneath the indentation surface. At a location deeper than 4 times the maximum impression depth, the stress around the hole is extremely weak to cause void collapse. Most of the dislocation loops are found to be parallel to the (001) plane, which is attributed to the low slip threshold of the (010) [100] slip system. This result shows that the LSCG strategy can deal with much larger systems and reveal mechanisms on the mesoscale.

1. INTRODUCTION Molecular crystals are an important class of substances that are widely applied in pharmaceuticals,1 energetics,2 electronics,3 and biomaterials,4 among others. Plastic deformation significantly influences the preparation and application of crystals, which can induce various physical and chemical changes in a short timescale.5,6 For example, the mechanical strength of pharmaceutical crystals can greatly be enhanced via plastic deformation during compaction into tablets.7 Plastic deformation and associated dislocation motion are responsible for the increased sensitivity of energetic materials to accidental initiation.8 Microstructure damage that evolves plastic deformation may result in photoelectric transfer hindrance in electronic materials, which eventually leads to final functional failure.9 These phenomena stress the significance of studying plastic deformation in molecular crystal materials. The underlying mechanism of plastic deformation remains unknown because it lacks experimental observations of elementary deformation processes at any strain rate.10 Meanwhile, molecular dynamics (MD) simulations are available and © 2016 American Chemical Society

complementary. MD simulations offer an advantage to obtain atomistic pictures of the plastic zone generated, including many details of dislocation generation, propagation, and vanishing. Thus, MD simulations have been applied extensively and have become a main tool for exploring the microscopic world.11−17 However, MD simulations for plastic behaviors of molecular crystals are more challenging than those for monoatomic, ionic, and metallic materials, as molecular crystals are often more structurally complicated. For example, additional issues, such as molecular conformations and orientations, exist.18,19 Therefore, the computation degrees of freedom are largely increased, and extending atomistic computer simulations to the real scale of deformation mechanisms, such as dislocation, microtwin, and stacking faults, is nearly impossible.20 Thus, modeling plastic deformation of molecular crystalline materials via the atomistic Received: April 27, 2016 Revised: June 24, 2016 Published: June 27, 2016 15198

DOI: 10.1021/acs.jpcc.6b04256 J. Phys. Chem. C 2016, 120, 15198−15208

Article

The Journal of Physical Chemistry C

Figure 1. AA-to-CG mapping scheme for α-RDX. One CG bead corresponds to one RDX molecule, and the centers of the CG beads are defined by the real centroid of the corresponding atomistic molecule.

Several interesting works have already applied CG approaches to molecular crystals. Izvekov et al. reported a multiscaled coarse-grain (MSCG) protocol based on FM theory and applied it to cyclotrimethylene trinitramine (RDX), a molecular crystal used as an explosive.27 By the CG approach, the lattice structure, elastic properties, and pressure−volume (P−V) curve were perfectly predicted. On the basis of Izvekov’s work, Brennan et al.28 developed a variant of DPD including chemical reactions and applied it to simulate shock compression and thermal decomposition of RDX. The CG potential parameters of these works were derived by matching the structural and thermodynamic properties of the CG model to those from atomistic calculations through a flexible force field developed by Smith and Bharadwaj.29 They used noncrystalline structures under high temperatures to avoid SDF issues raised by the ordered structure. Recently, Moore et al. followed up Izvekov’s work by employing a density dependence and a conservative force field form that conserves energy. They accurately predicted the density, structure, pressure−volume isotherm, bulk modulus, and elastic constants of the RDX crystal at ambient pressure, exhibiting the transferability of the force field to a liquid phase at melt conditions.30 In the present work, we develop a limited-sample coarsegrained (LSCG) protocol to acquire the interaction potential among CG particles by quantum chemical calculations. This research aims to find a new approach to construct the CG potentials of molecular crystals based directly on their ordered structures instead of AA simulations and corresponding potentials. This is a structure-matching approach25,26 because the interactive potential is derived on the basis of the structural nature of the given molecule. Our protocol directly uses the results of quantum chemistry calculation and is independent of AA potentials.22−26 We still choose RDX as the reference system. We employ CGMD simulations to calculate the elastic properties and gain insight into the plastic deformation mechanism under nanoindentation conditions. We discuss the influences of inner defects on the deformation behavior during indentation by presetting a nanohole at different sites in a perfect crystal model. The calculated properties, including lattice parameters, Young’s moduli (E), shear moduli (G), the P−V curve, and bulk modulus (B0) and its first derivatives (B0′ ), are comparable to the experimental and simulated results reported previously,27,31,49,52−54 verifying the feasibility of the LSCG potential of RDX. Furthermore, the nanoindentation simulations on the (100) surface of α-RDX indicate that when a nanohole is located at a depth less than twice the maximum

method is still a challenge, and large-scale methods are required. A feasible approach to reduce the overall degrees of freedom is by merging groups of chemically connected atoms into union atoms and to produce effective coarse-grained (CG) interaction potentials by averaging over microscopic details of the atomistic models. This is the so-called CG method, which allows computer simulations to be run on length and time scales that are 2−3 orders of magnitude larger compared to those in allatomistic (AA) simulations. The CG method also provides a bridge between microscopic and mesoscopic scales.21 Dissipative particle dynamics (DPD) and coarse-grained molecular dynamics (CGMD) are two conventional CG methods, and we are concerned with the latter one in the present work. CGMD approaches, including MARTINI,22 SDK,23 force matching (FM),24 iterative Boltzmann inversion (IBI),25 and relative entropy framework,26 have been developed in the past two decades. These approaches are widely used in soft condensed materials, such as polymers, surfactants, and biomaterials. However, they are rarely used in stiff molecular crystals such as energetic molecular crystals. The main issue of a CG simulation, as that in an AA simulation, is to construct effective interaction potentials. Constructing effective interaction potentials is the most difficult and complicated task in the entire simulation process. CG potentials are generally parameterized by reproducing structural properties, such as radial distribution function (RDF) and other distribution functions of bonds, angles, and torsions, and comparing the structural properties to those of AA simulations. Furthermore, CG potential parameters are developed by a series of iterative adjustments starting from guessed potentials.22−26 Two prerequisites should be provided before applying these approaches. First, the underlying AA interactions should be readily available, and second, all of the structural distribution functions (SDFs) should be consecutive and smooth. The CG potentials are based on the predetermined AA potentials; thus, the accuracy of CG potentials is dominated by that of the selected AA potentials. With regard to SDFs, they are often unavailable in the case of molecular crystals, as a consecutive and smooth SDF is available only in a disordered phase. Given the three-dimensional ordered structure nature, molecular crystals will inevitably present their SDF curves as periodical leapfrogging. For example, a classic RDF curve of a crystal material usually appears like shrill peaks, which are neither consecutive nor smooth. Application of CG approaches to crystal materials is still challenging for this reason. 15199


Article

The Journal of Physical Chemistry C impression depth (MID), local stress can trigger void collapse. When the location of the nanohole is deeper than 4 times the MID, the stress wave is extremely weak to cause void collapse. Most of the dislocation loops are parallel to the (001) plane, which is attributed to the low slip threshold of the (010)/[100] slip system. These results are consistent with the experimental and atomic simulation results of the slip characteristics of RDX crystals,10,65−67 validating that the LSCG method can effectively handle complex issues, such as plastic deformation of molecular crystals. The rest of this paper is organized as follows. Section 2 describes the principle of LSCG potential. Section 3 comprises the modeling process, simulation parameters, illustration of models, and software tools. The results and discussion are presented in Section 4, and conclusions are drawn in Section 5.

conditions.7 Thus, the conformational influence on total potential energy is negligible in the present work. The mutual orientations of a molecule pair show numerous possibilities (Figure 2). As displayed in the figure, each

2. PRINCIPLE OF LSCG The present LSCG protocol comprises five procedures: (1) mapping from AA to CG, (2) picking out pair samples, (3) scanning potential curves through quantum chemical calculations, (4) fitting the potential parameters, and (5) averaging on the potential parameters of all of the pair samples. These procedures are illustrated in detail as follows. Three polymorphs have been identified for RDX. The αpolymorph is the most stable under ambient conditions and crystallizes in the orthorhombic space group Pbca, with eight molecules per unit cell. The lattice parameters determined by neutron diffraction are a = 1.3182 nm, b = 1.1574 nm, and c = 1.0709 nm.31 A 21-to-1 mapping scheme is used in the present work, which means that each RDX molecule is represented by a single CG particle (Figure 1). Once the mapping scheme is established, the interaction potential between CG particles needs to be derived by averaging over all atomistic degrees of freedom contained in them. Constructing an analytical potential energy function that accurately represents the interparticle interactions is important in this work. This step influences the calculations of the structural parameters and mechanical properties. The potential function of an organic system is the sum of bonded and nonbonded terms.

Figure 2. Mutual orientations of a molecule pair.

Etot = E bonded + Enonbonded

molecule may rotate around any axis passing through its centroid. Thus, numerous possible values exist for their intermolecular potential energies. A CG potential has to average over numerous orientations, which can be expressed as x

E(rij) = lim

x →∞

1 ∑ f (rij, oijx) x 1

(3)

where oijx is the xth orientation of i and j or the xth molecular pair. According to the function in eq 3, the CG potential may be derived from averaging over a sufficient number of paired samples, which would be a formidable task. To simplify this task as far as possible, we attempt to develop an LS strategy. Suppose the pairs in an actual crystal are stable, they are regarded as pair samples. Then, on the basis of centroids of molecules (c. m.), the actual pairs are picked out from the experiment-measured structures of α-RDX, with a cutoff of 0.9 nm (Figure 3). Eventually, six independent pair types can be extracted from Figure 3 by removing the repetitive ones. The corresponding

(1)

where Etot is the total potential energy, Ebond is the bonded term, and Enonbond is the nonbonded term. Given that an entire molecule is merged into one bead, the intramolecular interactions, such as bonds, angles, dihedrals, and out-ofplane interactions, are frozen; only nonbonded interactions contribute to the total potential energy. Therefore, the potential function between CG beads can be derived from atomistic intermolecular interactions. The intermolecular potential energy between two organic molecules is defined by distance and orientation and conformation. It is represented as Etot = Enonbond = f (rij, oij , ci , cj) ≈ f (rij, oij)

(2)

where i and j stand for molecules i and j, respectively, rij is the distance between centroids of i and j, oij is the mutual orientation factor, and ci and cj are conformation factors of i and j. Previous experimental and theoretical studies have showed that α-RDX is stable below 3.77 GPa and 466 K,32,33 and its conformational fluctuations are undetectable under ambient

Figure 3. Intermolecular pairs picked out from the α-RDX crystal. Atomistic molecules are displayed in the line style, and their c. m.’s, overlapped with those of the CG beads, are displayed as balls. The selected pairs are interpreted as interactions between the central molecule (in red) and its neighbors (in green). 15200


Article

The Journal of Physical Chemistry C statistical weight, wi, of each pair type is provided by the RDF curve of α-RDX on the basis of c. m. (Figure 4), which is listed

Figure 5. Potential energy curves scanned at the B3LYP/6-311G** level. Pi stands for the ith molecular pair. Figure 4. RDF curve of α-RDX on the basis of c. m. (cutoff = 0.9 nm).

6

r0 =

∑ wri 0i

(7)

i=1

in the last row of Table 1. The potential energy curve of each pair type is scanned at the B3LYP/6-311G** level (Figure 5), which has been proven by previous reports to be reliable in calculating interaction potentials.29,34−38 These potential energy curves are parameterized according to the Morse equation39 E = D0[e−2α(r − r0) − 2e−α(r − r0)]

3. MODELING AND SIMULATION DETAILS Two groups of simulations were performed. One group was for calculation of mechanical properties to verify the feasibility of the current LSCG potential. The other was for resembling nanoindentation on the perfect and defected crystals, to reveal the related plastic deformation mechanism. A 5 × 5 × 5 supercell of α-RDX, with 1000 CG beads, was established. Then, a normal pressure and temperature (NPT) (P = 1 atm and T = 298 K) MD simulation was run for 50 ps. Therefore, the final equilibrated cell was used as the initial one for elastic constants and P−V curve calculations. Young’s (E) and shear (G) moduli were derived from the stress−strain relationship versus a series of uniaxial compression simulations along the three axes. The compression was up to 1%, with increments of 0.1%. Each compressed cell was equilibrated by a normal volume and temperature (NVT) simulation (T = 298 K) for 50 ps. The P−V curve of RDX at ambient temperature was obtained by a series of NPT MD simulations over the pressure range from 0.2 to 5 GPa, with increments of 0.2 GPa. At each pressure point, the cell was relaxed at 298 K for 200 ps. Subsequently, the bulk modulus was obtained by fitting the P− V curve to the third order Birch−Murnaghan equation43

(4)

where r is the c. m. distance between two RDX molecules, r0 is the equilibrium distance, D0 is the well depth, and α controls the “width” of the potential. Given the rigorous physical meaning, Morse potential has been applied to CG simulation since 200340 and has achieved good consistency of various properties, including density, pressure, isothermal compressibility, and viscosity.41,42 Thus, Morse potential was employed in constructing the current CG protocol for RDX to simulate its mechanical properties. The parameterized potentials and weight factors of the pairs are listed in Table 1. Finally, the LSCG potential parameters are calculated to be D0 = 6.919, α = 1.338, and r0 = 6.466, using eqs 5, 6, and 7, respectively. 6

D0 =

∏ D0i w

i

(5)

i=1 6

α=

∑ wiαi

(6)

i=1

Table 1. Potential Parameters Fitted According to the Morse Equationa

a

Pi

P1

P2

P3

P4

P5

P6

D0i αi r0i wi

13.1896 1.3333 4.2603 16.6756

6.3697 1.5763 6.3385 15.6505

5.1472 1.1638 6.4700 15.1296

6.2604 1.3604 6.8725 13.4781

8.9108 1.3526 7.0889 24.4317

3.5455 1.2222 7.6922 14.6988

Pi is the type name of pair i; D0i, αi, and r0i are the potential parameters; and wi is the weight factor. 15201


Article

The Journal of Physical Chemistry C −7/3 ⎡ ⎛ V ⎞−5/3⎤ 3B0 ⎢⎛ V ⎞ ⎥ −⎜ ⎟ p= ⎜ ⎟ ⎥⎦ 2 ⎢⎣⎝ V0 ⎠ ⎝ V0 ⎠

⎧ ⎫ ⎡⎛ ⎞−2/3 ⎤⎪ ⎪ V 3 ⎥ ⎨1 + (B0′ − 4)⎢⎜ ⎟ ⎬ −1 ⎢⎣⎝ V0 ⎠ ⎥⎦⎪ ⎪ 4 ⎩ ⎭

(8)

where B0, B0′, p, V, and V0 are the bulk modulus, pressure derivative of B0, hydrostatic pressure, volume at current pressure, and volume of the equilibrated cell, respectively. With respect to the nanoindentation simulation, a system that contains a cuboid simulation box filled with RDX CG beads as a substrate and a 90° conical diamond indenter was built (Figure 6). The total beads in the system amount to 6 912

Figure 7. A schematic plot of the four void models. Void 1 (VR1) is 20 nm below the indentation center, void 4 (VR4) is 40 nm below the indentation center, void 2 (VR2) is located at the site at which void 1 moves along the [010] orientation by 20 nm, and void 3 (VR3) is located at the site at which void 1 moves along the [001] orientation by 20 nm. Every hole is spherical, with a diameter equal to 8 nm.

related plastic behaviors were analyzed using a dislocation extraction algorithm (DXA).46 All simulations were carried out using a Nose−Hoover thermostat, with a time step of 10 fs and a cutoff distance for nonbonded interactions of 2 nm. The simulations were performed with the large-scale atomic/ molecular massively parallel simulator code,47 and the local atomic structures were visualized with the open visualization tool code.48

Figure 6. Crystallographic orientations of RDX samples and indentation loading direction.

000, and their sizes are sufficiently large to avoid image interactions and to accurately reflect plastic deformation. The indenter consists of 157 054 beads that are arranged as a diamond anvil structure, and it remains rigid during simulations. Furthermore, all beads in the indenter move collectively as a rigid body in the simulations. As illustrated in Figure 6, indentation was performed along the y axis of the substrate or the (100) face of the RDX crystal. The interaction between the indenter and substrate was described by a purely repulsive potential, which was a conventional treatment11,44,45 to simplify the interactions between indenter and sample in indentation simulations. Before the loading indentation, the substrate was relaxed by an NPT MD simulation (0 GPa and 298 K) for 100 ps. In subsequent simulations, the five bottom layers of the substrate were frozen to realize the suppression on the substrate during the indentation. Overall, each nanoindentation simulation comprised three phases. The first phase was uploading. The indenter was pushed into the substrate along the y direction and the loading force was increased from 0 to 275 nN, with increments of 5.5 nN. After each uploading, the system was relaxed at 298 K by an NVT simulation for 100 ps. Then, the loading with 275 nN was held for 1 ns. The final phase was unloading to 0, with increments of 11 nN, in which the indenter was gradually moved out of the substrate. To investigate the influence of void defect on the plastic behaviors during the nanoindentation process, four substrates were built, each with a nanohole sited differently relative to the location and impression orientation of the indenter tip (Figure 7). Therefore, similar simulations were conducted, and the

4. RESULTS AND DISCUSSION 4.1. Lattice Structure, Elastic Properties, and P−V Curve. The lattice parameters were calculated after the initial NVT−NPT dynamics process reached equilibrium at 298 K and 1 atm; thereafter, they were compared with those in the experimental work by Choi and Prince,31 the AAMD simulation by Sorescu et al.,49 and the MSCG simulation by Izvekov et al.27 As shown in Figure 8, the sites of CG beads are geometrically in accordance with the experimental atomistic structure. As listed in Table 2, the unit volume, V, calculated by LSCG potential, was overestimated by 4.8%. The calculated lattice parameters, a and b, are consistent with the reported experimental and simulation values, and c is overestimated by 8%. Generally, the accuracy of LSCG is close to that of MSCG and worse than that of AAMD. Similar cases of structural deviation or accuracy decrease were also reported by other CG works.40−42,50,51 This result indicates that structural simplification in CG models gives rise to low accuracy because of missing detailed structural information. Young’s moduli, shear moduli, and the bulk modulus of αRDX, calculated via the current LSCG method, were compared with those in earlier reports (Table 3). Experimental52,53 and simultated27,54 data were used. The elastic properties calculated by the LSCG method are generally consistent with those from experimental measures and simulations, and the calculated values of Young’s moduli are better than those of shear moduli. Moreover, G12 underestimates shear moduli. 15202


Article


Figure 9. Comparison of the p−V curve derived from experiment and that derived from the present simulation.

during the past decades over the issue of plastic deformation of RDX by simulations, such as shear bands in response to shock,55,56 stress response and decomposition under nanoindentation, and so forth. In this article, we are interested in the deformation mechanisms of α-RDX crystals under nanoindentation conditions. Nanoindentation examination is now a highly effective tool to explore plastic deformation behaviors on a nanoscale.57−59 With the aid of other detection methods, several processes, such as dislocation, shear, and phase transformation, can be examined.60 In the recent years, nanoindentation has been adopted to study the stress response and decomposition of RDX.10,60,61 Armstrong and Elban61 investigated the relation between dislocation and fracture of RDX crystals, and they revealed that the orthorhombic crystal structure offers intrinsic resistance to dislocation motion. Chen et al.60 studied the nanoindentation behavior of RDX crystals using reactive-forcefield MD simulations and found that dramatic heating and decomposition of RDX molecules occur from the substrate, and the maximum displacement of RDX molecules occurs on the (210) plane. Ramos et al.10 investigated elastic−plastic transitions along the (210), (021), and (001) planes of single RDX crystals. They found that indentation on the (210) and (001) planes causes deformation pileup, featuring an association of the zone axes of slip planes. The current nanoindentation simulations aim to investigate the plastic deformation mechanism of RDX crystals, influenced by nanoholes at various depths and orientations. The mechanism differences among the perfect RDX (PR) crystal and four defect crystals will be discussed in terms of load−depth curves, impression shapes, and dislocation characteristics. During the uploading stage, as illustrated in Figure 10, the load−depth curves of five crystals increased differently. Evident sudden displacements or “pop-ins” were observed on the load− depth force curves of VR1 and VR4, whereas the curves are relatively smooth for the remaining models. As shown in Figure 7, the voids in VR1 and VR4 were directly beneath the indentation center, which should be the reason for the obvious pop-ins. Before the first displacement burst, the slope of the load−depth force curve for VR4 was larger than that for VR1, indicating that the material was softened with a decrease in void depth. For example, the location of the void in VR4 was deeper than that in VR1. Our results are consistent with recent nanoindentation simulations on metal materials that stated that at the elastic stage the slope of the load−depth curves increases with an increase in the depth of the void.62−64

Figure 8. Crystal structure of α-RDX (green spheres) calculated by the LSCG potential at 298 K, which is compared with the experimental structure (shown as stick bonds) and its centroid (red spheres).

Table 2. Lattice Parameters of α-RDX at 298 K and 1 atm, Compared with Experimentally Determined Parameters and Those from Earlier Simulations a

Exp. AAMDb MSCGc LSCGd a

a, nm

b, nm

c, nm

V, nm3

ΔV/V, %

1.318 1.340 1.396 1.343

1.157 1.180 1.209 1.097

1.071 1.073 1.013 1.163

1.633 1.696 1.709 1.713

2.8 4.6 4.8

From ref 31. bFrom ref 48. cFrom ref 27. dThis work.

Table 3. Comparison of Elastic Properties (Unit in GPa) of α-RDX at 298 K and 1 atm to Those from Experimental Determinations and Simulations Exp. E1 E2 E3 G23 G31 G12 B0 B0′ a

a

20.9 16.0a 15.6a 5.2a 4.1a 6.9a 12.1c 8.6c

AAMD b

19.3 16.2b 18.5b 7.6b 4.9b 3.0b 13.0d 9.2d

MSCG b

19.7 19.4b 21.9b 7.2b 7.1b 4.1b

LSCGe 18.3 14.7 15.0 5.4 4.6 4.3 11.49 11.09

From ref 52. bFrom ref 27. cFrom ref 53. dFrom ref 54. eThis work.

We also examined the p−V curve under low-pressure conditions. Figure 9 shows that the curve derived from LSCG is highly consistent with that determined experimentally.53 Overall, the lattice parameters, elastic properties, and p− V curves derived from LSCG simulations are comparable to those from other methods. These results proved that the current potential can reliably simulate the mechanical properties of α-RDX at the CG level, which is helpful to investigate the physical phenomena on a mesoscale. 4.2. Nanoindentation Simulations. As mentioned above, the proposed LSCG potential can describe physical phenomena on the mesoscale, which is helpful for discussing plastic deformation behaviors. A great deal of work has been reported 15203


Article


during the dislocation identification process, as displayed in Figure 11. Figure 11 shows that the impression shape of VR2 is evidently different from that of the remaining cases, as a sinkhole lies on the surface and along the [010] orientation. VR2 is speculated to be related to slip behaviors of the RDX crystal and will be discussed later. The aforementioned DXA46 was employed to analyze the formation and propagation of dislocations during indentation. Burgers vectors were determined through this procedure, and a group of loop representations of the dislocation network could be obtained. The evolution of the dislocation network in PR is displayed in Figure 12. Figure 12 shows the dynamic process in the dislocation network during the three stages of indentation. At the load-holding stage, the dislocation network was propagating with the increment of the impression load. At the load-holding stage, the dislocation network was stably maintained, and after unloading, the dislocation network shrunk. Overall, the dislocation network was parallel to the (001) plane, which indicates that the propagation of dislocations occurs in the (001) plane. The initial dislocation was nucleated right below the indentation center when the shear stress approaches the theoretical slip threshold.44 A number of experimental works have been performed to correlate deformation features with possible slip systems in RDX crystals. Those works indicate that the (010) plane is the primary slip plane and that (011) and (021) are additional potential slip planes.10,66,67 The [100] orientation was suggested to be a cross-slip direction and suggests that (010) [100] was the most activated slip system. The current (001) indentation provides a well-suited condition for the (010) [100] slip. The particles in the (010) plane slipped along the [100] orientation as a whole, and no relative motion between particles existed in this plane. Consequently, the propagation of dislocation in the (010) plane was not feasible, and the main part of dislocation networks was parallel to the (001) plane. The slip potential barriers of (010) and (001) under the [100] drag force were calculated from the current LSCG potential (Figure 13), which clearly illustrates a lower slip threshold of the (010) plane. Four defected models (Figure 7) were established to study the dislocation behavior during indentation when a void defect was preset in the perfect model. The dislocation loops after unloading are displayed in Figure 14, including the four defected models and the perfect one. The dislocation

Figure 10. Load−depth force curves of RDX crystals derived from nanoindentation simulations.

At the stage of load holding, they exhibited a decreasing order of load−depth: VR1 > (VR2, VR3) > (VR4, PR). Compared with the absence of a void, the presence of a void led to an increment in MID, and the magnitude of increment was controlled by the depth and site of the void. For a void depth of 20 nm, the increment was obvious, as shown by cases VR1, VR2, and VR3 in Figure 7. Meanwhile, for a void depth of 40 nm, nearly 4 times the MID, the increment was unperceivable. This result indicates that the influence of a void defect on MID can be neglected if the void depth is larger than 4 times the MID. Among the models with a void at the same depth, the largest increment of MID occurred when the void was situated directly beneath the indentation center. Figure 7 shows that the voids of VR1 are located 20 nm under the indentation center but those of VR2 and VR3 shifted by 20 nm along the x and z axes, respectively. When unloading, the load−depth order is initially maintained, and then, the curves of PR, VR2, VR3, and VR4 cross one another. Finally, their depths reach about 7 nm; VR1 always achieves the largest depth and reaches 7.8 nm in the end. This result indicates that VR1 produces the largest plastic deformation. The impression shapes of the five models can be depicted by their free surfaces, which were obtained via DXA analysis46,47

Figure 11. Impression shapes after unloading. 15204


Article


Figure 12. Snapshots of the dislocation network in PR. Projected in the (001) and (010) planes.

The dislocation network in VR2 shifted toward the hole and connected with it. In VR3, the shift of the dislocation loops was not obvious, and the dislocation network was nonadjacent to the hole. In VR4, the two types of dislocations were independent from each other. Overall, the defected models differed from the perfect one in terms of load characteristics and structure evolution under indentation. The differences were induced by the preset nanoholes, and the details of these deviations depended on the specific location of the nanohole. Theoretically, the stress wave induced by indentation can propagate toward the interior of the substrate, and the nanohole in VR1 is at a shallow position under the indentation surface. Thus, the strength of the stress wave around the nanohole in VR1 is stronger than that in the others. Stronger stress around a hole may trigger void collapse, which increases impression depth. Continued collapse in the load-holding stage may increase indenter displacement and lead to a broader loading platform. The hole in VR2 or VR3 was not right beneath the acting point that was the main propagation direction of the stress wave; thus, stress around the hole was not sufficiently strong to trigger void collapse. The nanohole in VR4 was located at a deeper position, nearly 4 times the MID, at which the strength of the stress wave became too weak to cause collapse but triggered few type II dislocation loops. Surface sinking in VR2 was also attributed to the (010) [100] slip system, where the nanohole and indent center were located in a plane perpendicular to (010). Thus, the hole was involved

Figure 13. Slip potential barriers of (010) and (001) under the [100] drag force.

evolutions of the four defected models can be observed in Figures s1−s4 of the Supporting Information. As displayed in Figure 14, the dislocation network can be divided into two types of dislocation loops. Type I occurs close to the impression surface, and type II occurs around the nanohole, which occurs only in VR1 and VR4. In VR1, the two types of dislocation loops were mixed with each other, and evident plastic deformation could be observed around the hole. 15205


Article


Figure 14. Snapshots of the dislocation loops after unloading.

that the present LSCG protocol is generally effective in calculating the structural characteristics and elastic properties of α-RDX. The proposed protocol can also simulate certain deformation behaviors. Using the proposed protocol, the deformation mechanism of α-RDX crystals was investigated by a series of large-scale nanoindentation simulations on the (100) surface. The influence of interior defects was discussed by presetting nanoholes in the simulation model at different depths and orientations. The stress around the nanohole may trigger void collapse when the nanohole is located at a shallow position and is directly beneath the indentation surface. At a location deeper than 4 times the MID, the stress around the hole was extremely weak to cause collapse. Most of the dislocation loops were parallel to the (001) plane, which was attributed to the low slip threshold of the (010) [100] slip system. Surface sinking in the indentation surface was also caused by the (010) [100] slip system when the nanohole was located at the orientation of dislocation movement. The hole provided an interior space for dislocation climb; thus, the upper particles moved with the evolution of dislocation along the

in the dislocation movement. Given that the nanohole in VR2 provided an interior space for the dislocation climb, the particles around it moved with the evolution of dislocation along the orientation of the stress wave. Thus, a surface sinkhole occurred.

5. CONCLUSIONS We proposed an LSCG protocol and applied it to α-RDX. The protocol was performed by averaging on the interaction energy of actual molecular pairs from experimentally measured structures. The potential energy curves were scanned at the B3LYP/6-311G** level, and the potential parameters were fit according to the Morse equation. The final potential parameters were derived to be D0 = 6.919, α = 1.338, and r0 = 6.466. The calculated structural characteristic and elastic properties, such as lattice parameters, Young’s moduli (E), shear (G) moduli, p−V curve, and bulk modulus B0 and its first derivatives B0′ , were highly consistent with the published results of earlier experimental and theoretical works. This consistency proved 15206


Article


emission in aluminum using multiscale molecular dynamics simulation and continuum modeling. J. Mech. Phys. Solids 2014, 65, 35−53. (14) Tsuru, T.; Shibutani, Y. Anisotropic effects in elastic and incipient plastic deformation under (001), (110), and (111) nanoindentation of Al and Cu. Phys. Rev. B 2007, 75, No. 035415. (15) Wang, W.; Zhong, Y.; Lu, K.; Lu, L.; McDowell, D. L.; Zhu, T. Size effects and strength fluctuation in nanoscale plasticity. Acta Mater. 2012, 60, 3302−3309. (16) Pal, A.; Picu, R. C. Rotational defects in cyclotrimethylene trinitramine (RDX) crystals. J. Chem. Phys. 2014, 140, No. 044512. (17) Zamiri, A. R.; De, S. Deformation distribution maps of β-HMX molecular crystals. J. Phys. D 2010, 43, No. 035404. (18) Baer, M. R.; Hall, C. A. Isentropic loading experiments of a plastic bonded explosive and constituents. J. Appl. Phys. 2007, 101, No. 034906. (19) Dick, J. J.; Hooks, D. E.; Menikoff, R.; Martinez, A. R. Elastic− plastic wave profiles in cyclotetramethylene tetranitramine crystals. J. Appl. Phys. 2004, 96, 374. (20) Zhou, Y.; Strachan, A. Thermal conduction in molecular materials using coarse grain dynamics: role of mass diffusion and quantum corrections for molecular dynamics simulations. J. Chem. Phys. 2009, 131, No. 234113. (21) Qian, H.; Carbone, P.; Chen, X.; Karimi-Varzaneh, H. A.; Liew, C. C.; Müller-Plathe, F. Temperature-transferable coarse-grained potentials for ethylbenzene polystyrene and their mixtures. Macromolecules 2008, 41, 9919−9929. (22) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (23) Shinoda, W.; DeVane, R.; Kleinb, M. L. Coarse-grained molecular modeling of non-ionic surfactant self-assembly. Soft Matter 2008, 4, 2454−2462. (24) Izvekov, S.; Voth, G. Multiscale coarse graining of liquid-state systems. J. Chem. Phys. 2005, 123, No. 134105. (25) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. J. Comput. Chem. 2003, 24, 1624−1636. (26) Chaimovich, A.; Shell, M. S. Coarse-graining errors and numerical optimization using a relative entropy framework. J. Chem. Phys. 2011, 134, No. 094112. (27) Izvekov, S.; Chung, P. W.; Rice, B. M. Particle-based multiscale coarse graining with density-dependent potentials: application to molecular crystals (hexahydro-1,3,5-trinitro-s-triazine). J. Chem. Phys. 2011, 135, No. 044112. (28) Brennan, J. K.; Lísal, M.; Moore, J. D.; Izvekov, S.; Schweigert, I. V.; Larentzos, J. P. Coarse-grain model simulations of nonequilibrium dynamics in heterogeneous materials. J. Phys. Chem. Lett. 2014, 5, 2144−2149. (29) Smith, G. D.; Bharadwaj, R. K. Quantum chemistry based force field for simulations of HMX. J. Phys. Chem. B 1999, 103, 3570−3575. (30) Moore, J. D.; Barnes, B. C.; Izvekov, S.; Lísal, M.; Sellers, M. S.; Taylor, D. E.; Brennan, J. K. A coarse-grain force field for RDX: density dependent and energy conserving. J. Chem. Phys. 2016, 144, No. 104501. (31) Choi, C. S.; Prince, E. The crystal structure of cyclotrimethylenetrinitramine. Acta Crystallogr. 1972, B28, 2857−2862. (32) Baer, B. J.; Oxley, J.; Nicol, M. The phase diagram of RDX (hexahydro-1,3,5-trinitro-s-triazine) under hydrostatic pressure. High Pressure Res. 1990, 2, 99−108. (33) Miller, P. J.; Block, S.; Piermarini, G. J. Effects of pressure on the thermal-decomposition kinetics, chemical-reactivity and phase-behavior of RDX. Combust. Flame 1991, 83, 174−184. (34) Vasilios, E. R.; Vasilios, S. M. Force field development for poly (dimethylsilylenemethylene) with the aid of ab initio calculations. J. Phys. Chem. B 2006, 110, 14929−14938. (35) Lin, F.; Wang, R. Systematic derivation of AMBER force field parameters applicable to zinc-containing systems. J. Chem. Theory Comput. 2010, 6, 1852−1870.

orientation of the stress wave. Therefore, our LSCG method can be used to build CG models and to deal with mechanical response issues on a mesoscale.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b04256. Snapshots showing the dislocation evolutions of the four defected models (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. Tel: 86-816-2488489 (J.L.). *E-mail: [email protected]. Tel: 86-816-2493506 (C.Y.Z.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful for the financial support from the National Natural Science Foundation of China (11572296 and 11502244) and Science & Technology Innovation Foundation of ICM (KJCX-201403).

■

REFERENCES

(1) Datta, S.; David, J. W. Crystal structures of drugs: advances in determination, prediction and engineering. Nat. Rev. Drug Discovery 2004, 3, 42−57. (2) Lipp, M. J.; Evans, W. J.; Baer, B. J.; Yoo, C. S. High-energydensity extended CO solid. Nat. Mater. 2005, 4, 211−215. (3) Horiuchi, S.; Yoshinori Tokura, Y. From our readers: materials science in secondary education: non-MRSEC initiatives. Nat. Mater. 2005, 7, 357. (4) Sanchez, C.; Arribart, H.; Guille, M. M. G. Biomimetism and bioinspiration as tools for the design of innovative materials and systems. Nat. Mater. 2005, 4, 277−288. (5) Zamiri, A. R.; De, S. Modeling the anisotropic deformation response of b-HMX molecular crystals. Propellants, Explos., Pyrotech. 2011, 36, 247−251. (6) Zhou, T.; Lou, J.; Song, H.; Huang, F. Anisotropic shock sensitivity in a single crystal δ-cyclotetramethylene tetranitramine: a reactive molecular dynamics study. Phys. Chem. Chem. Phys. 2015, 17, 7924−7935. (7) Munday, L. B.; Mitchell, R. L.; Knap, J.; Chung, P. W. Role of molecule flexibility on the nucleation of dislocations in molecular crystals. Appl. Phys. Lett. 2013, 103, No. 151911. (8) Dick, J. J.; Ritchie, J. P. Molecular mechanics modeling of shear and the crystal orientation dependence of the elastic precursor shock strength in pentaerythritol tetranitrate. J. Appl. Phys. 1994, 76, 2726. (9) Sherwood, J. N. Fifty years as a crystal gazer: life as an imperfectionist. Cryst. Growth Des. 2004, 4, 863−877. (10) Ramos, K. J.; Hooks, D. E.; David, F. B. Direct observation of plasticity and quantitative hardness measurements in single crystal cyclotrimethylene trinitramine by nanoindentation. Philos. Mag. 2009, 89, 2381−2402. (11) Wang, B.; Gao, Y.; Urbassek, H. M. Microstructure and magnetic disorder induced by nanoindentation in single-crystalline Fe. Phys. Rev. B 2014, 89, No. 104105. (12) Dudeck, K. J.; Marques, L. A.; Knights, A. P.; Gwilliam, R. M.; Botton, G. A. Sub-angstrom experimental validation of molecular dynamics for predictive modeling of extended defect structures in Si. Phys. Rev. Lett. 2013, 110, No. 166102. (13) Yamakov, V. I.; Warner, D. H.; Zamora, R. J.; Saether, E.; Curtind, W. A.; Glaessgen, E. H. Investigation of crack tip dislocation 15207


Article

The Journal of Physical Chemistry C (36) Laura, C. B.; Hernán, E. L.; Carlos, G. N.; Silvia, A. B. Density functional theory calculations of the molecular force field of L-ascorbic acid, vitamin C. J. Phys. Chem. A 2010, 114, 4997−5004. (37) Orlando, A.; William, L. J. Exploring solvent effects upon the menshutkin reaction using a polarizable force field. J. Phys. Chem. B 2010, 114, 8425−8430. (38) Neyertz, S.; Mathieu, D.; Khanniche, S.; Brown, D. An empirically optimized classical force-field for molecular simulations of 2,4,6-trinitrotoluene (TNT) and 2,4-dinitrotoluene (DNT). J. Phys. Chem. A 2012, 116, 8374−8381. (39) Morse, P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev. 1929, 34, 57−64. (40) Liew, C. C.; Mikami, M. A coarse-grained model for particle dynamics simulations of complex fluids. Chem. Phys. Lett. 2003, 368, 346−351. (41) Chiu, S. W.; Scott, H. L.; Jakobsson, E. A coarse-grained model based on Morse potential for water and n-alkanes. J. Chem. Theory Comput. 2010, 6, 851−863. (42) Zhang, N.; Zhang, P.; Kang, W.; Bluestein, D.; Deng, Y. Parameterizing the Morse potential for coarse-grained modeling of blood plasma. J. Comput. Phys. 2014, 257, 726−736. (43) Birch, F. Finite elastic strain of cubic crystals. Phys. Rev. 1947, 71, 809−824. (44) Ukwatta, A.; Achuthan, A. A molecular dynamics (MD) simulation study to investigate the role of existing dislocations on the incipient plasticity under nanoindentation. Comput. Mater. Sci. 2014, 91, 329−338. (45) Chen, Y.; Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P. Molecular dynamics nanoindentation simulation of an energetic material. Appl. Phys. Lett. 2008, 93, No. 171908. (46) Stukowski, A.; Bulatov, V. V.; Arsenlis, A. Automated identification and indexing of dislocations in crystal interfaces. Modell. Simul. Mater. Sci. Eng. 2012, 20, No. 085007. (47) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1−19. (48) Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO - the open visualization tool. Modell. Simul. Mater. Sci. Eng. 2010, 18, No. 015012. (49) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Intermolecular potential for the hexahydro-1,3,5-trinitro-1,3,5-s-triazine crystal (RDX): a crystal packing, Monte Carlo, and molecular dynamics Study. J. Phys. Chem. B 1997, 101, 798−808. (50) Rosch, T. W.; Brennan, J. K.; Izvekov, S.; Jan Andzelm, W. Exploring the ability of a multiscale coarse-grained potential to describe the stress-strain response of glassy polystyrene. Phys. Rev. E 2013, 87, No. 042606. (51) Sliozberg, Y. R.; Mrozek, A.; Schieber, J. D.; Kröger, M.; Lenhart, J. L.; Andzelm, J. W. Effect of polymer solvent on the mechanical properties of entangled polymer gels: coarse-grained molecular simulation. Polymer 2013, 54, 2555−2564. (52) Haussühl, S. Elastic and thermoelastic properties of selected organic crystals. Z. Kristallogr. 2001, 216, 339. (53) Olinger, B.; Roof, B.; Cady, H. H. Symposium international sur le comportement des milieux denses sous hautes pressions dynamiques; Commissariat a l’Energie Atomique Centre d’Etudes de Vajours: Paris, France, 1978; pp 3−8. (54) Munday, L. B.; Chung, P. W.; Rice, B. M.; Solares, S. D. Simulations of high-pressure phases in RDX. J. Phys. Chem. B 2011, 115, 4378−4386. (55) Cawkwell, M. J.; Sewell, T. D. Shock-induced shear bands in an energetic molecular crystal: application of shock-front absorbing boundary conditions to molecular dynamics simulations. Phys. Rev. B 2008, 78, No. 014107. (56) An, Q.; Liu, Y.; Zybin, S. V.; Kim, H.; Goddard, W. A. Anisotropic shock sensitivity of cyclotrimethylene trinitramine (RDX) from compress-and-shear reactive dynamics. J. Phys. Chem. C 2012, 116, 10198−10206.

(57) Zhang, L.; Ohmura, T. Plasticity initiation and evolution during nanoindentation of an iron−3% silicon crystal. Phys. Rev. Lett. 2014, 112, No. 145504. (58) Smith, G. S.; Tadmor, E. B.; Kaxiras, E. Multiscale simulation of loading and electrical resistance in silicon nanoindentation. Phys. Rev. Lett. 2000, 84, 1260−1263. (59) Ramamurty, U.; Jang, J. Nanoindentation for probing the mechanical behavior of molecular crystalsa review of the technique and how to use it. CrystEngComm 2014, 16, 12−23. (60) Chen, Y.; Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P. Molecular dynamics nanoindentation simulation of an energetic material. Appl. Phys. Lett. 2008, 93, No. 171908. (61) Armstrong, R. W.; Elban, W. L. Cracking at hardness microindentations in RDX explosive and MgO single crystals. Mater. Sci. Eng. A 1989, 111, 35−43. (62) Pejman, R.; Salehi, S. H.; Salehi, M. Numerical study of interfacial crack growth effects on nanoindentation mechanical properties in presence of pre-existing defect. Thin Solid Films 2013, 545, 408−413. (63) Tan, C.; Jeng, Y. Computer simulations of nanoindentation on Cu (111) with a void. Int. J. Solids Struct. 2009, 46, 1884−1889. (64) Zhu, P.; Hu, Y.; Wang, H. Atomistic simulations of the effect of a void on nanoindentation response of nickel. Sci. China: Phys., Mech. Astron. 2010, 53, 1716−1719. (65) Mathew, N.; Picu, R. C. Slip asymmetry in the molecular crystal cyclotrimethylenetrinitramine. Chem. Phys. Lett. 2013, 582, 78−81. (66) Munday, L. B.; Solares, S. D.; Chung, P. W. Generalized stacking-fault energy surfaces in the molecular-crystal alpha-RDX. Philos. Mag. 2012, 92, 3036−3050. (67) Mcdermott, I. T.; Phakey, P. P. A method of correlating dislocations and etch pits: application to cyclotrimethylene trinitramine. J. Appl. Crystallogr. 1971, 4, 479−481.

15208


Limited-Sample Coarse-Grained Strategy and Its Applications to

Recommend Documents