Microstructures
All Journals
Search Log In
Home Articles Special Issues Volumes Webinars Videos
Microstructures
Search Submit
Home Articles Special Issues Volumes Webinars Videos
Download PDF
Research Article  |  Open Access  |  21 May 2024

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Views: 846 |  Downloads: 535 |  Cited:  2
Zhong-Hao Ye1
, 
Feng Guo1,*
  , ... 
Xiao-Chun Wang1
J Mater Inf 2024;4:5.
10.20517/jmi.2023.37 |  © The Author(s) 2024.
Author Information
Article Notes
Cite This Article

Abstract

In this work, we report the discovery of energy cocrystals using an efficient iterative workflow combining an evolutionary algorithm and a machine learning potential (MLP). The compound 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20) has attracted significant attention owing to its higher energy density than traditional energetic materials. However, the higher sensitivity has limited its applications. An important way to reduce its sensitivity involves cocrystal engineering with traditional explosives. Many cocrystal structures are expected to be composed of these two components. We developed an efficient iterative workflow to explore the phase space of CL-20 and 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (HMX) cocrystals using an evolutionary algorithm and an MLP. The algorithm was based on the Universal Structure Predictor: Evolutionary Xtallography (USPEX) software, and the MLP was the reactive force field with neural networks (ReaxFF-nn) model. A set of high-density cocrystal structures was produced through this workflow; these structures were further checked via first-principles geometry optimizations. After careful screening, we identified several high-density cocrystal structures with densities of up to 1.937 g/cm3 and HMX: CL-20 ratios of 1:1 and 1:2. The influence of hydrogen bonds on the formation of high-density cocrystals was also discussed, and a roughly linear relationship was found between energy and density.

Graphical Abstract

Keywords

Crystal search, machine learning, ReaxFF-nn, evolutionary algorithm
Download PDF 0 4

INTRODUCTION

Crystal structure prediction (CSP) techniques enable the design of new crystal structures with a given chemical composition. However, the high computational costs of density functional theory (DFT) calculations limit their applications in molecular systems, whose unit cells usually contain dozens or even hundreds of atoms. Recently developed machine learning potentials (MLPs), such as the Gaussian approximation potential (GAP)[1–6], moment tensor potential (MTP)[7–9], and DeepMD-kit package[10–13], can perform DFT-level calculations at the computational cost of classical force fields. CSP methods, such as Universal Structure Predictor: Evolutionary Xtallography (USPEX)[14,15] and Crystal structure AnaLYsis by Particle Swarm Optimization (CALPSO)[16], combining MLPs[17–23] will undoubtedly become the main approaches in the search for new materials. The work by Podryabinkin[24] on CSP combining an MLP and an evolutionary algorithm showed that an accuracy of 11 meV/atom could be achieved in the known structure predictions. It is clear that MLPs[25] can accelerate CSP by replacing most DFT calculations. Although attempts[26,27] have been made to use the classical force field as a local optimizer for structure minimization, the parameter set used by the force field cannot be dynamically adjusted according to the search results and, therefore, cannot guarantee the predicted structure to be the best one. In contrast, machine learning-based methods can provide predictions with quality close to that of first-principles calculations at the computational costs of classical techniques and thus can deal with large-scale[28] systems. In addition, the parameter set can be dynamically updated through a search–check–train–search workflow developed by our group. A similar workflow is used in the active learning procedure used to train deep-learning potential models, such as the Deep Potential GENerator (DP-GEN)[29] and Fast Learning of Atomistic Rare Events (FLARE)[30] software, which contains three Exploration–Labeling–Training iterative steps. In the Exploration step, a sample of trajectories is generated using methods such as molecular dynamics (MD) simulations and randomly sampling the crystal phase space. The search step in our workflow is similar to the Exploration of the crystal space group by the USPEX software. In the Labeling step, and correspondingly in our check step, the ab initio energies are calculated for the configurations generated in the previous step, and the configurations with a large variance are added to the training dataset. In the last step, the model is trained, and a new parameter set is generated. In this work, we combined MLP and evolutionary algorithm to identify new cocrystal structures of 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20)[31] and 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (HMX)[32], aiming to find an optimal balance between their energy and safety. CL-20 is known for its high detonation power; however, significantly reducing its sensitivity is a great challenge. Effectively addressing this problem involves combining CL-20 with other low-sensitivity materials, producing more attractive explosive cocrystal materials[33–37]. A highly promising CL-20:HMX (2:1) cocrystal was synthesized by Bolton et al.[38]; the cocrystal exhibits superior energetics to HMX while maintaining comparable impact sensitivity. Many cocrystal structures composed of CL-20 and HMX are expected to exist because of the large number of crystal structures the single components can form.

As thousands of crystal structures will need to be optimized in one search, using the DFT method as the local optimizer will require years to complete. The machine learning-based approach will speed up the search and provide similar accuracy, as the parameter set is dynamically updated after the search is corrected by DFT calculations. The MLP will "learn" the interactions between molecules and give more accurate predictions. A workflow combining the USPEX evolutionary algorithm and the ReaxFF-nn MLP was developed to explore the configuration space of CL-20 and HMX. Using this approach, we report the discovery of high-density cocrystal structures of HMX and CL-20 with 1:1 and 1:2 ratios. Furthermore, there are discussions on whether C–H...O hydrogen bond interactions[39] exist in crystal structures; hence, the role of C–H...O hydrogen bond interactions in the cocrystal structures is discussed by analyzing hydrogen bond energies.

METHODS

ReaxFF-nn machine learning potential

The ReaxFF reactive force field, originally introduced by van Duin et al.[40], has been extensively used[41,42] in reactive MD simulations[43–51]. We developed[52] an MLP model labeled ReaxFF-nn based on ReaxFF with neural networks for bond orders in previous work[53]. This significantly improved various calculations, such as equations of states, reaction energies, and potential energy surfaces. The total energies for ReaxFF and ReaxFF-nn are respectively given in:

EsystemReaxFF=Ebond+Elp+Eover+Eunder+Eval+Epen+Ecoa+Etors+Econj+EHbond+EvdWaals+ECoulomb,

EsystemReaxFF−nn=Ebondnn+Eval+Etors+EHbond+EvdWaals+ECoulomb.

where Ebond, Elp, Eover, Eunder, Eval, Epen, Ecoa, Etors, Econj, EHbond, EvdWaals and ECoulomb correspond to the bond, lone pair, over coordinate, under coordinate, valence angle, penalty, three-body conjugation, torsion rotation barrier, four-body conjugation, hydrogen bond, van der Waals and Coulomb energy terms, respectively[40].

In the current form, as given in Equation (2), neural networks are employed to compute short-range interactions, while long-range interactions, such as Coulomb and van der Waals terms, still use the same formula as ReaxFF. Because most MLPs use prespecified static point charges, they focus on short-range interactions[54–56] and have difficulties dealing with long-range Coulomb interactions, which remain a major challenge for these potentials at present. An exception is represented by ReaxFF, and, similarly, ReaxFF-nn, which use charges dynamically calculated by the electronegativity equilibrium[57] method, and further related approaches are expected to be developed in the future[58].

The corrected bond order BOReaxFF and ReaxFF-nn is respectively calculated by:

BOReaxFF=BO′⋅f1(Δi′,Δj′)⋅f4(Δi′,BOij′)⋅f5(Δj′,BOij′),

where the f1, f4, and f5 terms are given by the ReaxFF formula in Ref.[40],

BOReaxFF−nn=BO′⋅fnn(Δi′,BOij′,Δj′),

fnn used by a hidden-layer neural network is calculated by:

fnn(x)=σ(wo⋅σ(wh⋅σ(wi⋅x+bi)+bh)+bo).

where x=(Δi′,BOij′,Δj′) represents the input vector, Δi′ is the sum of the BOij′ terms of atomi, and BOij′ is the uncorrected bond order, as defined in Ref.[40]; moreover, σ is the logistic activation function σ(x)=(1+exp⁡(−x))−1, while w and b are the weights and biases of the neural networks, respectively.

The bond energies employed by ReaxFF and ReaxFF-nn are respectively given in:

EbondijReaxFF=Deσ⋅BOijσ⋅exp⁡[pbe1(1−BOijσ)pbe2]−Deπ⋅BOijπ−Deππ⋅BOijππ,

EbondijReaxFF−nn=Deσ⋅fnne(BOijσ,BOijπ,BOijππ).

where Deσ, Deπ, Deππ, pbe1 and pbe2 are adjusting parameters, and BOijσ, BOijπ, and BOijππ are bond order components.

In this case, fnne represents a hidden-layer neural network that calculates the bond energy. The input vector for this network is the uncorrected bond order BOij=(BOijσ,BOijπ,BOijππ). The network output, multiplied by a scalable parameter Deσ, gives the bond energy of the pair of atoms atomi and atomj.

Training and testing the potential model

The quality of the dataset used for the MD interaction potential plays a crucial role in determining its accuracy. To ensure the complete sampling of the phase space of molecular structures such as CL-20,2,4,6-trinitrotoluene (TNT), 1,1-diamino-2,2-dinitroethylene (FOX-7), cyclotrimethylene trinitramine (RDX) and HMX, we manually collected a large amount of data using methods including ab initio MD simulations, stretching of specific chemical bonds, scanning covalent bond angles, and equations of state of a crystal structure. After DFT calculations, these data were stored in an atomic simulation environment (ASE)[59]Trajectory object, which can be interpreted as a time series containing many ASE Atoms objects moving in the configurational phase space. This data format is easy to prepare because the ASE package includes many DFT calculators, such as Vienna ab initio simulation package (VASP), Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA), and Quantum Espresso (QE). After performing calculations through these calculators, the Trajectory file can be generated automatically. One can also perform single-point energy calculations and gather the trajectory file after the calculations. The training function accepts a dataset variable, which is a Python dict object that contains all trajectory file names; moreover, the I-ReaxFF package can generate input information for the training function, such as collecting all bonds and all valence angles of the same type into variables that store this information. More details on the data preparation can be found in the documentation and examples of our I-ReaxFF package[53,60].

Figure 1A-C compares the lattice energies calculated by DFT and ReaxFF-nn as a function of the density for the β-HMX, ε-CL-20, and 1:1 HMX: CL-20 cocrystal. Scaling the box size in three directions changed the density from 1.4 to 2.4 g/cm3. Figure 1D and E shows the energies vs. the N–NO2 bond distance calculated by DFT and ReaxFF-nn. Figure 1F compares the lattice energies of a batch of 100 cocrystal structures of 1:1 HMX: CL-20 generated by the USPEX software. The results demonstrate a good agreement between our ReaxFF-nn potential model and the DFT calculations. Although active learning algorithms[29,61,62] to collect data automatically are available, the application to molecular systems requires further tests, and we also plan to develop an uncertainty-driven active learning algorithm based on the Z-matrix.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 1. Lattice energy as a function of density of (A) β-HMX, (B) ε-CL-20, and (C) 1:1 HMX: CL-20 cocrystal; Energy of molecular (D) HMX and (E) CL-20 as a function of N–NO2 bond distance; (F) Comparison of DFT- and ReaxFF-nn-calculated lattice energies of 1:1 cocrystal structures randomly generated by USPEX. HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; DFT: density functional theory; ReaxFF-nn: reactive force field with neural networks; USPEX: Universal Structure Predictor: Evolutionary Xtallography.

After data collection, the ReaxFF-nn model was trained using the I-ReaxFF software package, and the overall precision of this dataset was approximately 94.0%, as calculated by:

Accu=1.0−∑|EDFT−EReaxFF−nn|∑(EDFT−EZPE),

where EZPE represents the zero-point energy.

To further evaluate the reliability of the ReaxFF-nn potential model, we optimized the crystal structures of β-HMX (chair conformation[63]) and ε-CL-20 (the most stable crystal form of CL-20[27,38,64]). A comparison of the calculated lattice constants with the experimental data and the first-principles values is presented in Table 1. The analysis shows that the calculated values obtained using ReaxFF-nn agree well with the experimental and first-principles data.

Table 1

The comparison of lattice constants of ε-CL-20 and β-HMX from experiment, DFT calculated values, and GULP geometry optimization with potential ReaxFF-nn

a(Å)b(Å)c(Å)β(°)
ε-CL-20, Exp.[27,65]8.85212.55613.386106.82
ε-CL-20, Exp.[66,67]8.83712.55413.289106.87
ε-CL-20, Exp.[68]8.86312.59313.395106.92
ε-CL-20, DFT-D[69]8.87112.55713.477106.62
ε-CL-20, ReaxFF-nn8.98912.40113.428105.90
β-HMX, Exp.[70]6.52611.0377.364102.67
β-HMX, Exp.[71]6.53011.0307.350102.69
β-HMX, Exp.[72]6.54011.0507.371102.83
β-HMX, DFT-D[69]6.59310.8607.406102.84
β-HMX, ReaxFF-nn6.43410.5827.876103.40

CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; DFT: density functional theory; GULP: General Utility Lattice Program.

RESULTS AND DISCUSSION

Molecular CSPs were carried out using the USPEX package[14,15,73,74]. This package is based on the evolutionary algorithm developed by Oganov, Glass, Lyakhov, and Zhu. Initially, a group of 100 structures was randomly generated from the appropriate crystal space group. In subsequent generations, the total population of 100 structures was maintained after operations including heredity, randomization, permutation, rotation, and softmutation. The optimal structures of the previous generation were kept as keptbest[15] structures. Each structure underwent two local relaxation steps using the General Utility Lattice Program (GULP) package[75] and the ReaxFF-nn potential. The fitness of each structure was determined based on the corresponding energy calculated by GULP. Following the principles of biological evolution, structures with lower fitness were then eliminated. Figure 2 illustrates the main search workflow comprising four stages. First, the initial parameter set is generated by training the model against the existing dataset. Second, the USPEX algorithm is combined with the ReaxFF-nn MLP (structural relaxation) to search for cocrystals for specific needs, generating a set of crystal structures. Third, these crystal structures are checked by comparing the corresponding ReaxFF-nn- and DFT-calculated energies. Finally, structures with a very large difference between the calculated energies are added to the training dataset. These steps are repeated until a small difference between the ReaxFF-nn and DFT results or new cocrystal structures meeting the selected requirements have been found.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 2. Schematic illustration of four-step iterative workflow combining USPEX and ReaxFF-nn for searching the appropriate cocrystal structures and training the MLP. USPEX: Universal Structure Predictor: Evolutionary Xtallography; ReaxFF-nn: reactive force field with neural networks; MLP: machine learning potential.

Prediction of 1:2 HMX: CL-20 cocrystal configurations

After several searches, we obtained some cocrystal structures with lower enthalpy (potential energy) and higher density. However, to acquire densely packed cocrystals, we set a small pressure in the search process (p = 0.1 GPa). The obtained structures had to be optimized (relaxed) to zero pressure. The optimization results using the ReaxFF-nn and DFT (SIESTA[76]) method with the Perdew-Burke-Ernzerhof (PBE) functional and a double-ζ basis set at zero pressure are shown in Supplementary Table 1. Furthermore, an experimentally synthesized cocrystal structure was optimized by DFT (SIESTA) to serve as a reference, as shown in the first row of Supplementary Table 1.

Supplementary Table 1 reveals a good agreement between the cocrystal structures optimized by ReaxFF-nn and DFT (SIESTA), demonstrating the high performance of our MLP model. Combined with the data from the structures in Supplementary Table 1, the enthalpy and density of the cocrystal structures are shown as bar graphs in Figure 3. The enthalpy and density values calculated by SIESTA and ReaxFF-nn are shown in red and turquoise, respectively.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 3. (A) Enthalpy and (B) density vs. crystal index from USPEX of 1:2 HMX: CL-20 optimized by DFT (SIESTA) and MLP (ReaxFF-nn). USPEX: Universal Structure Predictor: Evolutionary Xtallography; HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; DFT: density functional theory; SIESTA: Spanish Initiative for Electronic Simulations with Thousands of Atoms; MLP: machine learning potential; ReaxFF-nn: reactive force field with neural networks.

In order to confirm the findings of the search workflow, we chose the five most significant cocrystal configurations with a 1:2 ratio and further optimized them using the van der Waals Density Functional (vdW-DF)[77] method and DFT-D3 method with the PBE functional and an energy cutoff of 600 eV. These DFT-D3 calculations were carried out using the VASP[78] software. The results are detailed in Table 2.

Table 2

The comparison of the enthalpy, density, and lattice constants from DFT (SIESTA), vdW-DF, and DFT-D3 (VASP) calculations of the top predicted 1:2 HMX: CL-20 cocrystal structures, respectively

IDMethodEnthalpy (eV)ρ(g⋅cm−3)a(Å)b(Å)c(Å)α(°)β(°)γ(°)
C415PBE-200.5801.88411.02612.0528.25184.2890.7171.44
vdW-DF-199.7321.92110.90112.0148.20683.9491.1371.74
PBE-D3-648.6751.90211.09512.1468.12083.3991.0670.63
C420PBE-200.5321.87010.86612.1288.31185.0389.1272.60
vdW-DF-199.7251.92310.82912.0288.23183.9290.5671.95
PBE-D3-648.6991.87111.09012.1988.19883.6690.1270.89
D477PBE-200.3361.86810.85812.1848.32384.5588.8271.93
vdW-DF-199.4601.92010.87411.5968.65679.7782.8671.18
PBE-D3-648.5731.90510.95111.7598.61579.1982.9070.00
E400PBE-200.3731.91210.77112.2318.406100.0483.6171.59
vdW-DF-199.3611.95410.88211.9508.31399.5983.2571.88
PBE-D3-648.9351.93410.89112.2108.25899.2882.9671.03
F111PBE-200.3291.91411.23412.1368.26898.7182.1768.92
vdW-DF-199.2601.99111.23911.7518.05698.0081.7371.38
PBE-D3-648.9761.93711.24412.1458.09998.1482.1069.64

DFT: Density functional theory; SIESTA: Spanish Initiative for Electronic Simulations with Thousands of Atoms; vdW-DF: van der Waals Density Functional; VASP: Vienna ab initio simulation package; HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; PBE: Perdew-Burke-Ernzerhof.

Based on the information on the DFT-D3-optimized cocrystal structures from Table 2, the five cocrystal structures with lower energy and higher density and the structure experimentally synthesized by Bolton et al.[38] are shown in Figure 4. The optimized cocrystal structures in this picture show the characteristic alternating arrangements of monolayer HMX and bilayer CL-20, consistent with the cocrystal arrangement of the experimentally synthesized structure.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 4. (A) Experimentally synthesized cocrystal structure[38], and (B–F) 1:2 HMX: CL-20 cocrystal structures after DFT-D3 structure optimization at zero pressure. The (B–F) structures correspond to IDs C415, C420, D477, E400, and F111 in Table 2, respectively. HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; DFT: density functional theory.

Among the structures with the predicted 1:2 HMX: CL-20 ratio obtained from USPEX and MLP, the structure F111 shown in Figure 4F had the highest density of 1.914 g/cm3 after DFT optimization and 1.937 g/cm3 after DFT-D3 optimization. Furthermore, the F111 cocrystal also exhibited the lowest lattice energy.

Prediction of 1:1 HMX: CL-20 cocrystal configurations

A search of 1:1 HMX: CL-20 cocrystal structures was carried out using USPEX and ReaxFF-nn through the search–check–train–search iterative loop, as illustrated in Figure 2. A low pressure was applied to the candidate cocrystal cells; in the next step, the generated structures were optimized using the DFT (SIESTA) method with the PBE functional and a double-ζ basis set, as well as with the ReaxFF-nn potential to zero pressure. Supplementary Table 1 shows the structural information on the optimized cocrystals with lower enthalpy and higher density for the 1:1 ratio at zero pressure. The enthalpy and density of the cocrystal structures after optimization, obtained using ReaxFF-nn and DFT, are shown in Figure 5, which enables us to identify structures with lower enthalpy and higher density.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 5. Comparison of DFT (SIESTA)- and ReaxFF-nn-calculated (A) enthalpy and (B) density values vs. crystal index of cocrystal structures with 1:1 HMX: CL-20 ratio generated by USPEX. DFT: Density functional theory; SIESTA: Spanish Initiative for Electronic Simulations with Thousands of Atoms; ReaxFF-nn: reactive force field with neural networks; USPEX: Universal Structure Predictor: Evolutionary Xtallography.

Similarly, to confirm the findings of the search workflow, we chose the six most significant cocrystal configurations with a 1:1 ratio and further optimized them using the vdW-DF method and DFT-D3 (VASP) method with the PBE functional and an energy cutoff of 600 eV. The results are detailed in Table 3. We selected the better cocrystal structures calculated by DFT-D3, as shown in Figure 6. The optimized structures in Figure 6 show characteristic cocrystal structures with alternating arrangements of HMX and CL-20 monolayers. Theoretically, this is consistent with the characteristic arrangement of cocrystal structures.

Table 3

Comparison of enthalpy, density, and lattice constants from DFT (SIESTA), vdW-DF, and DFT-D3 (VASP) calculations of the top predicted cocrystal structures with 1:1 ratio, respectively

IDMethodEnthalpy (eV)ρ(g⋅cm−3)a(Å)b(Å)c(Å)α(°)β(°)γ(°)
D3PBE-127.4471.9238.1448.4919.34398.3190.8382.73
vdW-DF-127.4161.9338.1208.4579.36398.5190.7182.75
PBE-D3-411.6871.9378.1418.4559.31898.4991.1582.81
D558PBE-127.4991.9258.2058.6669.05598.1090.5183.40
vdW-DF-127.5211.9278.2178.6459.06198.0490.3883.15
PBE-D3-411.4231.8988.3078.6879.06398.3890.9783.29
D559PBE-127.4511.9258.1198.4209.44398.3090.6982.56
vdW-DF-127.4151.9338.1378.3799.44298.3990.5082.22
PBE-D3-411.4591.8958.2078.4629.45098.6691.2282.76
F377PBE-127.5451.9148.0519.6048.40381.0483.0589.45
vdW-DF-127.5911.9308.0719.6378.29380.9682.7289.63
PBE-D3-411.6821.9198.1809.2988.53280.5883.0188.90
G1766PBE-127.4981.9008.8548.3788.70984.6186.4491.46
vdW-DF-127.5321.9188.9188.3158.64784.3585.5791.32
PBE-D3-411.5471.9068.9538.3198.65484.1686.5791.41
H223PBE-127.3541.87911.6908.3646.70991.5796.9485.42
vdW-DF-127.2591.90811.5238.3876.73291.1998.9383.99
PBE-D3-411.6201.92111.4488.2826.82791.4099.5383.95

DFT: Density functional theory; SIESTA: Spanish Initiative for Electronic Simulations with Thousands of Atoms; vdW-DF: van der Waals Density Functional; VASP: Vienna ab initio simulation package; PBE: Perdew-Burke-Ernzerhof.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 6. Cocrystal structures with 1:1 HMX: CL-20 ratio after DFT-D3 (VASP) optimization. Structures (A–F) correspond to IDs D3, D558, D559, F377, G1766 and H223 in Table 3, respectively. HMX: 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane; CL-20: 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane; DFT: density functional theory; VASP: Vienna ab initio simulation package.

The results show the D558 and D559 structures had the highest density of 1.925 g/cm3 after DFT (SIESTA) optimization. However, the DFT-D3 calculations with the VASP software gave a different result: the D3 structure shown in Figure 6A exhibited the highest density of 1.937 g/cm3 and the lowest lattice energy. To our knowledge, there have been no previous theoretical or experimental reports of 1:1 HMX: CL-20 cocrystal structures, and this is the first report of 1:1 HMX: CL-20 cocrystal structures with a density higher than that of the β-HMX (d = 1.901 g/cm3)[70] crystal.

Hydrogen bond interactions

Hydrogen bond interactions play an important role in the design of molecular crystals; in practice, we found that increasing the C–H...O hydrogen bond contribution to the total energy can apparently raise the rate of high-density cocrystals. This leads to the conclusion[79] that the C–H...O interaction can be employed to design molecular crystals. Figure 7 shows the relationship between hydrogen bond energy (calculated by the ReaxFF-nn potential) and crystal density. The figure shows that the C–H...O hydrogen bond energy represents the main part of the total hydrogen bond energy; moreover, an approximately linear relationship is observed between the hydrogen bond energy and the density of the cocrystals, with a correlation factor η=0.67, which indicates that the C–H...O hydrogen bond energy is the key factor in the design of HMX/CL-20 cocrystals.

Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential

Figure 7. Relationship between hydrogen bond energy and density of predicted cocrystals. The hydrogen bond energy was calculated using the ReaxFF-nn potential. ReaxFF-nn: Reactive force field with neural networks.

CONCLUSIONS

In this work, we report an evolutionary algorithm and machine learning workflow denoted as search–check–train–search to explore new cocrystal structures of HMX and CL-20, which are expected to have higher energy than β-HMX and lower sensitivity than ε-CL-20. Using the present workflow, several high-density HMX/CL-20 energetic cocrystal structures with 1:1 and 1:2 HMX: CL-20 ratios were found and validated via DFT and DFT-D3 calculations. The most valuable finding was the D3 cocrystal structure, which had the highest density and lowest lattice energy among the 1:1 cocrystal structures. The density of D3 (1.937 g/cm3) was only slightly lower than that of 1:2 HMX: CL-20 (1.945 g/cm3) reported by Bolton et al.[38] at room temperature. Another high-density cocrystal with a 1:1 ratio was H223, which had a similar lattice energy and a density of 1.921 g/cm3. Among cocrystals with 1:2 HMX: CL-20 ratio, F111 had the lowest lattice energy and the highest density of 1.937 g/cm3. These high-density cocrystal structures can be accessed on the website [80]. In practice, we found that the C–H...O hydrogen bond interactions play an important role in the search for high-density cocrystals. Although these interactions are weak, with an equilibrium hydrogen bond length of 2.254 Å estimated through EOS calculations, they account for the vast majority of the total hydrogen bond energy of HMX/CL-20 cocrystals. An approximately linear relationship was found between the hydrogen bond energy and the density of cocrystal structures.

DECLARATIONS

Authors' contributions

Designed this project: Guo F

Complete all the calculations: Ye ZH

Provided technical and material support: Wen YS, Chai CG, Zhang ZR

Offered writing suggestions: Wang XC, Li HS, Zhang GQ, Cui SX

Availability of data and materials

The raw data (high-density cocrystals and parameter set for LAMMPS) are available from the website[80].

Financial support and sponsorship

This work was supported by the National Natural Science Foundation of China (NSAF. 12002325 and 11504153) and the Natural Science Foundation of Gansu Province (Grant No. I23JRRA1173). We greatly appreciate the assistance.

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2024.

Supplementary Materials

REFERENCES

1. Bartók AP, Csányi G. Gaussian approximation potentials: a brief tutorial introduction. Int J Quantum Chem 2015;115:1051-7.

2. Dragoni D, Daff TD, Csányi G, Marzari N. Achieving DFT accuracy with a machine-learning interatomic potential: thermomechanics and defects in bcc ferromagnetic iron. Phys Rev Mater 2018;2:013808.

3. Bartók AP, Payne MC, Kondor R, Csányi G. Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons. Phys Rev Lett 2010;104:136403.

4. Szlachta WJ, Bartók AP, Csányi G. Accuracy and transferability of Gaussian approximation potential models for tungsten. Phys Rev B 2014;90:104108.

5. Botu V, Ramprasad R. Learning scheme to predict atomic forces and accelerate materials simulations. Phys Rev B 2015;92:094306.

6. Li Z, Kermode JR, De Vita A. Molecular dynamics with on-the-fly machine learning of quantum-mechanical forces. Phys Rev Lett 2015;114:096405.

7. Novikov IS, Gubaev K, Podryabinkin EV, Shapeev AV. The MLIP package: moment tensor potentials with MPI and active learning. Mach Learn Sci Technol 2021;2:025002.

8. Shapeev AV. Moment tensor potentials: a class of systematically improvable interatomic potentials. Multiscale Model Simul 2016;14:1153-73.

9. Novikov IS, Suleimanov YV, Shapeev AV. Automated calculation of thermal rate coefficients using ring polymer molecular dynamics and machine-learning interatomic potentials with active learning. Phys Chem Chem Phys 2018;20:29503-12.

10. Zeng J, Zhang D, Lu D, et al. DeePMD-kit v2: a software package for deep potential models. J Chem Phys 2023;159:054801.

11. Li J, An Q. Quasiplastic deformation in shocked nanocrystalline boron carbide: grain boundary sliding and local amorphization. J Eur Ceram Soc 2023;43:208-16.

12. Pitike KC, Setyawan W. Accurate Fe-He machine learning potential for studying He effects in BCC-Fe. J Nucl Mater 2023;574:154183.

13. Zhai B, Wang HP. Accurate interatomic potential for the nucleation in liquid Ti-Al binary alloy developed by deep neural network learning method. Comput Mater Sci 2023;216:111843.

14. Lyakhov AO, Oganov AR, Stokes HT, Zhu Q. New developments in evolutionary structure prediction algorithm USPEX. Comput Phys Commun 2013;184:1172-82.

15. Zhu Q, Oganov AR, Glass CW, Stokes HT. Constrained evolutionary algorithm for structure prediction of molecular crystals: methodology and applications. Acta Crystal Sect B Struct Sci 2012;68:215-26.

16. Wang Y, Lv J, Zhu L, Ma Y. CALYPSO: a method for crystal structure prediction. Comput Phys Commun 2012;183:2063-70.

17. Yamashita T, Kanehira S, Sato N, et al. CrySPY: a crystal structure prediction tool accelerated by machine learning. Sci Technol Adv Mater Methods 2021;1:87-97.

18. Li C, Wang C, Sun M, et al. Correlated RNN framework to quickly generate molecules with desired properties for energetic materials in the low data regime. J Chem Inf Model 2022;62:4873-87.

19. Choudhary K, DeCost B, Chen C, et al. Recent advances and applications of deep learning methods in materials science. npj Comput Mater 2022;8:59.

20. Hong C, Choi JM, Jeong W, et al. Training machine-learning potentials for crystal structure prediction using disordered structures. Phys Rev B 2020;102:224104.

21. Chen WC, Schmidt JN, Yan D, Vohra YK, Chen CC. Machine learning and evolutionary prediction of superhard B-C-N compounds. npj Comput Mater 2021;7:114.

22. Chen WC, Vohra YK, Chen CC. Discovering superhard B-N-O compounds by iterative machine learning and evolutionary structure predictions. ACS Omega 2022;7:21035-42.

23. Kruglov IA, Yanilkin A, Oganov AR, Korotaev P. Phase diagram of uranium from ab initio calculations and machine learning. Phys Rev B 2019;100:174104.

24. Podryabinkin EV, Tikhonov EV, Shapeev AV, Oganov AR. Accelerating crystal structure prediction by machine-learning interatomic potentials with active learning. Phys Rev B 2019;99:064114.

25. Bereznikova LA, Propad YV, Kruglov IA. Nitrogen phase diagram at high P-T conditions by the T-USPEX method. J Phys Chem C 2023;127:5683-8.

26. Pakhnova M, Kruglov I, Yanilkin A, Oganov AR. Search for stable cocrystals of energetic materials using the evolutionary algorithm USPEX. Phys Chem Chem Phys 2020;22:16822-30.

27. Wang C, Ni Y, Zhang C, Xue X. Crystal structure prediction of 2,4,6,8,10,12-Hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20) by a tailor-made OPLS-AA force field. Cryst Growth Des 2021;21:3037-46.

28. Wespiser C, Mathieu D. Application of machine learning to the design of energetic materials: preliminary experience and comparison with alternative techniques. Prop Explos Pyrotech 2023;48:e202200264.

29. Zhang Y, Wang H, Chen W, et al. DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models. Comput Phys Commun 2020;253:107206.

30. Vandermause J, Torrisi SB, Batzner S, et al. On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events. npj Comput Mater 2020;6:20.

31. Venkata Viswanath J, Venugopal KJ, Srinivasa Rao NV, Venkataraman A. An overview on importance, synthetic strategies and studies of 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (HNIW). Def Technol 2016;12:401-18.

32. Taylor JW, Crookes RJ. Vapour pressure and enthalpy of sublimation of 1,3,5,7-tetranitro-1,3,5,7-tetra-azacyclo-octane (HMX). J Chem Soc Faraday Trans 1 1976;72:723-9.

33. Tariq QuN, Tariq MuN, Dong WS, Manzoor S, Arshad F, Zhang JG. Comparative studies of synthesis, performance, and applications of recently developed CL-20 based co-crystals. Cryst Growth Des 2023;23:6974-87.

34. Yang Z, Li H, Zhou X, Zhang C, et al. Characterization and properties of a novel energetic-energetic cocrystal explosive composed of HNIW and BTF. Cryst Growth Des 2012;12:5155-8.

35. Xu H, Duan X, Li H, Pei C. A novel high-energetic and good-sensitive cocrystal composed of CL-20 and TATB by a rapid solvent/non-solvent method. RSC Adv 2015;5:95764-70.

36. Liu N, Duan B, Lu X, et al. Preparation of CL-20/DNDAP cocrystals by a rapid and continuous spray drying method: an alternative to cocrystal formation. Cryst Eng Comm 2018;20:2060-7.

37. Wang Y, Yang Z, Li H, et al. A novel cocrystal explosive of HNIW with good comprehensive properties. Prop Explos Pyrotech 2014;39:590-6.

38. Bolton O, Simke LR, Pagoria PF, Matzger AJ. High power explosive with good sensitivity: a 2:1 cocrystal of CL-20:HMX. Cryst Growth Des 2012;12:4311-4.

39. Corpinot MK, Bučar DK. A practical guide to the design of molecular crystals. Cryst Growth Des 2019;19:1426-53.

40. van Duin ACT, Dasgupta S, Lorant F, Goddard WA. ReaxFF: A reactive force field for hydrocarbons. J Phys Chem A 2001;105:9396-409.

41. Senftle TP, Hong S, Islam MM, et al. The ReaxFF reactive force-field: development, applications and future directions. npj Comput Mater 2016;2:15011.

42. Han Y, Jiang D, Zhang J, Li W, Gan Z, Gu J. Development, applications and challenges of ReaxFF reactive force field in molecular simulations. Front Chem Sci Eng 2016;10:16-38.

43. Chenoweth K, van Duin ACT, Goddard WA. ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. J Phys Chem A 2008;112:1040-53.

44. Guo F, Zhang H, Hu HQ, Cheng XL. Effects of hydrogen bonds on solid state TATB, RDX, and DATB under high pressures. Chin Phys B 2014;23:046501.

45. Liu J, Li X, Guo L, et al. Reaction analysis and visualization of ReaxFF molecular dynamics simulations. J Mol Graphics Modell 2014;53:13-22.

46. Zhou T, Song H, Liu Y, Huang F. Shock initiated thermal and chemical responses of HMX crystal from ReaxFF molecular dynamics simulation. Phys Chem Chem Phys 2014;16:13914-31.

47. Huang X, Guo F, Yao K, et al. Anisotropic hydrogen bond structures and orientation dependence of shock sensitivity in crystalline 1,3,5-tri-amino-2,4,6-tri-nitrobenzene (TATB). Phys Chem Chem Phys 2020;22:11956-66.

48. Daksha CM, Yeon J, Chowdhury SC, Gillespie JW Jr. Automated ReaxFF parametrization using machine learning. Comput Mater Sci 2021;187:110107.

49. Bu Y, Guo F, Li K, et al. High-temperature pyrolysis behavior and structural evolution mechanism of graphene oxide: a ReaxFF molecular dynamics simulation. Appl Surf Sci 2022;593:153451.

50. Jiang J, Wang HR, Zhao FQ, Xu SY, Ju XH. Decomposition mechanism of 1,3,5-trinitro-2,4,6-trinitroaminobenzene under thermal and shock stimuli using ReaxFF molecular dynamics simulations. Phys Chem Chem Phys 2023;25:3799-805.

51. Feng S, Guo F, Yuan C, et al. Effect of neutron irradiation on structure and decomposition of α-RDX: a ReaxFF molecular dynamics study. Comput Theor Chem 2023;1219:113965.

52. Xue LY, Guo F, Wen YS, et al. ReaxFF-MPNN machine learning potential: a combination of reactive force field and message passing neural networks. Phys Chem Chem Phys 2021;23:19457-64.

53. Guo F, Wen YS, Feng SQ, et al. Intelligent-ReaxFF: evaluating the reactive force field parameters with machine learning. Comput Mater Sci 2020;172:109393.

54. Friederich P, Häse F, Proppe J, Aspuru-Guzik A. Machine-learned potentials for next-generation matter simulations. Nat Mater 2021;20:750-61.

55. Schütt KT, Sauceda HE, Kindermans PJ, Tkatchenko A, Müller KR. SchNet - A deep learning architecture for molecules and materials. J Chem Phys 2018;148:241722.

56. Yoo P, Sakano M, Desai S, Islam MM, Liao P, Strachan A. Neural network reactive force field for C, H, N, and O systems. npj Comput Mater 2021;7:9.

57. Rappe AK, Goddard WA III. Charge equilibration for molecular dynamics simulations. J Phys Chem 1991;95:3358-63.

58. Islam MM, Kolesov G, Verstraelen T, Kaxiras E, van Duin ACT. eReaxFF: a pseudoclassical treatment of explicit electrons within reactive force field simulations. J Chem Theory Comput 2016;12:3463-72.

59. Larsen AH, Mortensen JJ, Blomqvist J, et al. The atomic simulation environment - a Python library for working with atoms. J Phys: Condens Matter 2017;29:273002.

60. Guo F. I-ReaxFF: intelligent-reactive force field. Available from: https://github.com/fenggo/I-ReaxFF. [Last accessed on 15 May 2024].

61. Sivaraman G, Krishnamoorthy AN, Baur M, et al. Machine-learned interatomic potentials by active learning: amorphous and liquid hafnium dioxide. npj Comput Mater 2020;6:104.

62. Shi B, Zhou Y, Fang D, et al. Estimating the performance of a material in its service space via Bayesian active learning: a case study of the damping capacity of Mg alloys. J Mater Inf 2022;2:8.

63. Landenberger KB, Matzger AJ. Cocrystals of 1,3,5,7-Tetranitro-1,3,5,7-tetrazacyclooctane (HMX). Cryst Growth Des 2012;12:3603-9.

64. Millar DIA, Maynard-Casely HE, Allan DR, et al. Crystal engineering of energetic materials: co-crystals of CL-20. Cryst Eng Comm 2012;14:3742-9.

65. Bidault X, Chaudhuri S. A flexible-molecule force field to model and study hexanitrohexaazaisowurtzitane (CL-20) - polymorphism under extreme conditions. RSC Adv 2019;9:39649-61.

66. Zhang XQ, Chen XR, Kaliamurthi S, Selvaraj G, Ji GF, Wei DQ. Initial decomposition of the co-crystal of CL-20/TNT: sensitivity decrease under shock loading. J Phys Chem C 2018;122:24270-8.

67. Zhang XQ, Yuan JN, Selvaraj G, Ji GF, Chen XR, Wei DQ. Towards the low-sensitive and high-energetic co-crystal explosive CL-20/TNT: from intermolecular interactions to structures and properties. Phys Chem Chem Phys 2018;20:17253-61.

68. Bolotina NB, Hardie MJ, Speer RL Jr, Pinkerton AA. Energetic materials: variable-temperature crystal structures of γ- and ε-HNIW polymorphs. J Appl Cryst 2004;37:808-14.

69. Cawkwell MJ, Zecevic M, Luscher DJ, Ramos KJ. Dependence of the elastic stiffness tensors of PETN, α-RDX, γ-RDX, ε-RDX, ε-CL-20, DAAF, FOX-7, and β-HMX on hydrostatic compression. Prop Explos Pyrotech 2022;47:e202100281.

70. Deschamps JR, Frisch M, Parrish D. Thermal expansion of HMX. J Chem Crystal 2011;41:966-70.

71. Eiland PF, Pepinsky R. The crystal structure of cyclotetramethylene tetranitramine. Z Krist Cryst Mater 1954;106:273-98.

72. Cady HH, Larson AC, Cromer DT. The crystal structure of α-HMX and a refinement of the structure of β-HMX. Acta Cryst 1963;16:617-23.

73. Oganov AR, Glass CW. Crystal structure prediction using ab initio evolutionary techniques: principles and applications. J Chem Phys 2006;124:244704.

74. Oganov AR, Lyakhov AO, Valle M. How evolutionary crystal structure prediction works - and why. Acc Chem Res 2011;44:227-37.

75. Gale JD, Raiteri P, van Duin ACT. A reactive force field for aqueous-calcium carbonate systems. Phys Chem Chem Phys 2011;13:16666-79.

76. Soler JM, Artacho E, Gale JD, et al. The SIESTA method for ab initio order-N materials simulation. J Phys Condens Matter 2002;14:2745.

77. Dion M, Rydberg H, Schröder E, Langreth DC, Lundqvist BI. Van der Waals density functional for general geometries. Phys Rev Lett 2004;92:246401.

78. Kresse G, Furthmüller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 1996;54:11169-86.

79. Gavezzotti A, Presti LL. Building blocks of crystal engineering: a large-database study of the intermolecular approach between C-H donor groups and O, N, Cl, or F acceptors in organic crystals. Cryst Growth Des 2016;16:2952-62.

80. Guo F. ReaxFF-nn for lammps. Available from: https://github.com/fenggo/ReaxFF-nn_for_lammps. [Last accessed on 15 May 2024].

Cite This Article

Research Article
Open Access
Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential
Zhong-Hao Ye, ... Xiao-Chun Wang

How to Cite

Ye, Z. H.; Guo, F.; Chai, C. G.; Wen, Y. S.; Zhang, Z. R.; Li, H. S.; Cui, S. X.; Zhang, G. Q.; Wang, X. C. Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential. J. Mater. Inf. 2024, 4, 5. http://dx.doi.org/10.20517/jmi.2023.37

Download Citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click on download.

Export Citation File:

Type of Import

Tips on Downloading Citation

This feature enables you to download the bibliographic information (also called citation data, header data, or metadata) for the articles on our site.

Citation Manager File Format

Use the radio buttons to choose how to format the bibliographic data you're harvesting. Several citation manager formats are available, including EndNote and BibTex.

Type of Import

If you have citation management software installed on your computer your Web browser should be able to import metadata directly into your reference database.

Direct Import: When the Direct Import option is selected (the default state), a dialogue box will give you the option to Save or Open the downloaded citation data. Choosing Open will either launch your citation manager or give you a choice of applications with which to use the metadata. The Save option saves the file locally for later use.

Indirect Import: When the Indirect Import option is selected, the metadata is displayed and may be copied and pasted as needed.

About This Article

Special Issue

This article belongs to the Special Issue Application of Machine Learning to Crystal Structure Prediction
© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, sharing, adaptation, distribution and reproduction in any medium or format, for any purpose, even commercially, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Data & Comments

Data

Views
846
Downloads
535
Citations
2
Comments
0
4

Comments

Comments must be written in English. Spam, offensive content, impersonation, and private information will not be permitted. If any comment is reported and identified as inappropriate content by OAE staff, the comment will be removed without notice. If you have any queries or need any help, please contact us at support@oaepublish.com.

0
Download PDF
Download XML 6 downloads
Cite This Article 8 clicks
Export Citation 4 clicks
Like This Article 4 likes
Share This Article
Scan the QR code for reading!
See Updates
Contents
Figures
Related
Microstructures
ISSN 2770-2995 (Online)
Navigation
Follow Us
Navigation

Committee on Publication Ethics

https://members.publicationethics.org/members/microstructures

Portico

All published articles are preserved here permanently:

https://www.portico.org/publishers/oae/

Committee on Publication Ethics

https://members.publicationethics.org/members/microstructures

Portico

All published articles are preserved here permanently:

https://www.portico.org/publishers/oae/
Discover Content
Follow OAE
Twitter
Facebook
LinkedIn
YouTube
BiLiBiLi
WeChat
Rss
© 2016-2025 OAE Publishing Inc., except certain content provided by third parties

Your privacy, your choice

We use essential cookies to make sure the site can function. We also use optional cookies for advertising, personalisation of content, usage analysis, and social media.
By accepting optional cookies, you consent to the processing of your personal data - including transfers to third parties. Some third parties are outside of the European Economic Area, with varying standards of data protection. See our privacy policy for more information on the use of your personal data.