Next Article in Journal
Exploring the Biological and Phytochemical Potential of Jordan’s Flora: A Review and Update of Eight Selected Genera from Mediterranean Region
Previous Article in Journal
Cationic Vitamin E-TPGS Mixed Micelles of Berberine to Neutralize Doxorubicin-Induced Cardiotoxicity via Amelioration of Mitochondrial Dysfunction and Impeding Apoptosis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals

1
School of Chemistry and Chemical Engineering, North University of China, Taiyuan 030051, China
2
Gansu Yinguang Chemical Industry Group, Baiyin 730900, China
*
Author to whom correspondence should be addressed.
Molecules 2024, 29(5), 1159; https://doi.org/10.3390/molecules29051159
Submission received: 18 February 2024 / Revised: 26 February 2024 / Accepted: 28 February 2024 / Published: 5 March 2024
(This article belongs to the Section Materials Chemistry)

Abstract

:
External electric fields are an effective tool to induce phase transformations. The crystallization of ionic crystals from solution is a common phase transformation. However, understanding of mechanisms is poor at the molecular level. In this work, we carried out an experimental and theoretical investigation of the external electric-field-induced crystallization of TKX-50 from saturated formic acid solution by finite-temperature string (FTS) with order parameters (OPs) as collective variables for ionic crystals. The minimum-free-energy path was sketched by the string method in collective variables. The results show that the K-means clustering algorithm based on Euclidean distance and density weights can be used for enhanced sampling of the OPs in external electric-field-induced crystallization of ionic crystal from solution, which improves the conventional FTS. The crystallization from solution is a process of surface-mediated nucleation. The external electric field can accelerate the evolution of the string and decrease the difference in the potential of mean forces between the crystal and the transition state. Due to the significant change in OPs induced by the external electric field in nucleation, the crystalline quality was enhanced, which explains the experimental results that the external electric field enhanced the density, detonation velocity, and detonation pressure of TKX-50. This work provides an effective way to explore the crystallization of ionic crystals from solution at the molecular level, and it is useful for improving the properties of ionic crystal explosives by using external electric fields.

1. Introduction

Recently, external electric fields have been an effective tool to induce variations in phase transformations [1,2], with the formation of a potential new crystal phase and morphology, so as to realize control of the production process and improvements in the structures and properties of materials [3,4,5,6,7]. As mentioned by Alexander et al. [8], electric fields can reduce the nucleation time [9,10,11], increase the product yield [12,13,14], enhance overall crystal quality [10,13,15,16], and control the location of nucleation [16], the product crystal size [10,17,18], the crystal orientation [19,20,21] and polymorphism [22,23,24,25], etc.
At the molecular level, understanding the effects of electric fields on the nucleation mechanism can potentially and essentially control and enhance the crystalline product [1,5,7]. To this end, a large number of investigations have been carried out on the effects of electric fields on nucleation. In the experimental studies, Mahmood et al. reported the electric-field-induced supramolecular phase transition [5]. The electric-field-induced phase transition in perovskite relaxor ferroelectric crystals [26] and Li doping in VO2 film [27] were studied. For the chiral magnet Cu2OSeO3, the metastable skyrmion lattice emerged under electric fields [28]. As for the theoretical investigations conducted using ab initio density functional theory, the effects of electric fields on the phase boundaries, crystal growth rates, nucleation rates, and interfacial free energies have been studied [29,30]. The electric-field-tuned topological phase transition was investigated by first-principles calculations [31,32,33], with a change in the dipole moment [34]. Moreover, resistive switching in phase change materials [35] has recently attracted considerable attention for its application in non-volatile memory devices, and the magnetic properties have been studied [36]. Despite the extensive investigations being carried out, it remains poorly understood how the molecules aggregate and organize into a new form under external electric fields [1,3,4,5,6]. The main reason is that it exceeds the range of instrument testing for small molecules in experiments, and in theoretical simulations, the time scale is many orders of magnitude lower than that of real occurrence [37,38]. In fact, the nucleation itself is a rare event [39,40], and it is very difficult to find a set of suitable collective variables that can describe the reaction coordinate [39], which is the challenge of nucleation simulation and understanding nucleation mechanisms at the molecular level [41].
As an important rare-event method, finite-temperature string (FTS) has received great interest in exploring the phase transformation mechanisms [42,43,44,45,46,47]. With the reparameterization of the images along the path of the phase transformation [48], the minimum-free-energy path (MFEP) [49] could be obtained with the string method in collective variables (SMCV) [49] to pursue the principal curve. Furthermore, as a function of the collective motions of the system [50], the order parameters (OPs) involving only the structural parameters [39] could be used to the collective variable in describing the reaction coordinate for the phase transformation, and a corresponding MFEP could be obtained with SMCV at the molecular level. The OPs have been successfully used to monitor nucleation [39,51,52]. Moreover, K-means clustering is an algorithm by which a dataset containing n-dimensional vectors is clustered and divided into k sub-clusters by the iterations [53]. As a sampling method, it has been widely used in machine learning [54,55]. Therefore, the theoretical investigation of the phase transformation mechanism by the finite-temperature string method with order parameters as collective variables and the K-means clustering algorithm has been feasible. Recently, to improve the convergence rate of FTS with OP as collective variables for molecular crystal, we used a K-means clustering algorithm based on the Euclidean distance and sample weights to smooth the initial FTS and the minimum-free-energy path connecting between β-CL-20 and ε-CL-20 was sketched, and the free-energy profile along the path was obtained [56]. The crystallization of ionic crystals from solution is a common phase transition phenomenon, and this technique has been used abundantly in industry, with practical applications in pharmaceutical, chemical, petrochemical, food industries, and biotechnology, making it a multitrillion-dollar business [22,57,58,59,60,61]. To our knowledge, however, no theoretical investigation into the crystallization of ionic crystal from solution by the FTS method with the K-means clustering algorithm and OPs as the collective variables under the external electric field is present.
1,1’-dihydroxy-5,5’-tetrazolium dihydroxyamine salt (TKX-50) is a novel high-energy-density compound and one of the most promising third-generation energetic materials. As one of the important properties of explosives, the crystal morphology of TKX-50 has been extensively studied from theoretical simulation [62] to experiment [63] and from single component crystal [64] to co-crystal [65]. Industrial grade TKX-50 has many shortcomings, such as irregular crystal morphology, small particle size, large aspect ratio, and low intramolecular oxygen content [66], which affects its application in mixed explosives and propellants. Recrystallization from solution is a key process in the production of TKX-50 and is a common method to improve crystallization quality and application performance.
In this work, an experimental investigation of the crystallization of TKX-50 from the supersaturated formic acid solution was carried out under a direct current (DC) external electric field. To explain the experimental results, a theoretical study of the external electric-field-induced crystallization of TKX-50 from solution was presented by using the finite-temperature string method with order parameters as collective variables for ionic crystals. The local OPs of the TKX-50 crystal and the corresponding forms in the formic acid solution were constructed by the structural variables involving the bond distances, bond orientations, and relative orientations. A K-means clustering algorithm based on the dimensionally weighted Euclidean distance and sample weights was used to construct the smooth initial FTS. An MFEP involving the TKX-50 crystallization from the solution was sketched by using the SMCV method, and then the potential of mean force (PMF) was calculated under the external electric field. The essence of the crystallization of TKX-50 from the supersaturated formic acid solution was revealed. This study extends the application scope of FTS in the investigation into the crystallization from solution for the explosive ionic crystals by the K-means clustering technique with OPs as the collective variables and provides an effective way to explore the crystallization from the solution for the explosive ionic crystals at the molecular level. It is useful for further optimizing the technological process to obtain the desired explosive ionic crystals under the external electric field in the experiment.

2. Theory

2.1. Order Parameters

To gain a molecular-level understanding of the way molecules aggregate and organize themselves into crystal structures, Santiso et al. presented a method to construct OPs that involve only the structural parameters [39]. They are a set of mathematical functions containing the intermolecular distance, bond orientation, relative orientation, internal configurations, etc., and these OPs can be used to characterize the translation, vibration, and rotation of the functional groups as well as the rearrangement of molecules. Among them, the bond distance OPs, bond orientation Ops, and relative orientation OPs, represented as φ α , i d , φ α , i b , and φ α , i r , have been used to describe the formation of crystals in solution [39]. They are defined as follows:
φ α , i d = j i 1 2 π σ α exp [ ( r ^ i j r ^ α ) 2 2 σ α 2 ]
φ α , i b = 1 2 π I 0 ( η r ^ α ) exp [ η r ^ α cos 2 ( ϕ r ^ ϕ r ^ α ) ]
φ α , i r = 1 2 π I 0 ( η q α ) exp [ η q α cos 2 ( ϕ q ϕ q α ) ]
where i denotes the ith molecule; rij is the distance between the centers-of-mass of two molecules; ϕ r ^ is defined as the bond orientation from the projection of r ^ i j onto the absolute orientation of the ith molecule; ϕ q is the relative orientation resulting from the rotation of the absolute orientation of the ith molecule onto jth molecule; and rα, ϕ r ^ α , and ϕ q α are the mean center-of-mass distance, mean bond orientation, and mean relative orientation corresponding to the α-peak, respectively. σ α , η r ^ α , and η q α are the standard deviation and concentration parameters, respectively, and I0 is the modified Bessel function of the second kind and order 0. Thus, the OPs of the bond distance combined with the bond orientation or relative orientation could be given by
φ α , i db = 1 2 π σ α 1 2 π I 0 ( η r ^ α ) j i exp [ ( r ^ i j r ^ α ) 2 2 σ α 2 ] exp [ η r ^ α cos 2 ( ϕ r ^ ϕ r ^ α ) ]
φ α , i dr = 1 2 π σ α 1 2 π I 0 ( η q α ) j i exp [ ( r ^ i j r ^ α ) 2 2 σ α 2 ] exp [ η q α cos 2 ( ϕ q ϕ q α ) ]
To simplify the OPs, two kinds of local OPs were often adopted [51]. One is the OP for a molecule i defined by summing over the different α peaks, and the other is defined by the averaged value within the divided cell (C) in the simulation box and summed over all the peaks, given by Equations (6) and (7), respectively:
θ i * = α φ i , α *
θ C * = 1 N C i C α φ i , α *
where * represents “d”, “b”, “r”, “db”, and “dr” in Equations (1)–(5), respectively.

2.2. Finite-Temperature String

Finite-temperature String is a technique to determine the minimum energy paths or MFEPs by evaluating locally the mean force and tensor, which accounts for the curvilinear nature of the collective variables [49,67].
Given a set of N collective variables used to describe the reaction of interest, θ(x) = (θ1(x), θ2(x), …, θN(x)), where xRn are the Cartesian coordinates, the function of the free energy depending on z = (z1, z2, …, zN) is defined as
F ( z ) = k B T ln ( Z 1 R n e β V ( x ) δ ( z 1 θ 1 ( x ) ) δ ( z N θ N ( x ) ) d x )
where V(x) is the potential energy, Z = R n e β V ( x ) d x , and β = 1/kBT, where kB and T are the Boltzmann constant and the temperature, respectively. Note that N is less than the dimensionality of the full system.
For a string denoted by z(α, t) at the evolutionary time t, where α ∈ [0,1] is used to parametrize the string, the evolution equation involving the free-energy gradient to make it converge to MFEP could be given by
z i ( α , t ) t = j , k = 1 N P i j ( α , t ) M j k ( z ( α , t ) ) F ( z ( α , t ) ) z k
where Pij(α, t) is the projector on the plane perpendicular to the path at z(α, t), Mjk(z(α, t)) is the tensor, and F ( z ( α , t ) ) z k (i.e., F ( z ) ) is the free-energy gradient (i.e., mean force). Given two minima of the free energy located at za and zb, Equation (9) can be solved according to the boundary conditions z(0, t) = za and z(1, t) = zb. As t→∞ the solution of (9) converges to an MFEP connecting za and zb, and simultaneously, the tangent vector is parallel to M ( z ) F ( z ) at every point z, i.e.,
0 = j , k = 1 N P i j ( α ) M j k ( z ( α ) ) F ( z ( α ) ) z k
where P i j ( α ) = δ i j t ^ i ( α ) t ^ j ( α ) , t ^ i ( α ) = z i / α | z i / α | , and Mjk(z) is defined as
M i j ( z ) = Z 1 e β F ( z ) R n k = 1 n θ i ( x ) x k θ j ( x ) x k e β V ( x ) δ ( z 1 θ 1 ( x ) ) δ ( z N θ N ( x ) ) d x
Based on their estimators under the ergodicity with the large enough values of k and T, Equations (10) and (11) could be solved to obtain [49]:
F ( z ) z j k T 0 T ( z j θ j ( x ( t ) ) ) d t
M ij ( z ) 1 T k = 1 n 0 T θ i ( x ( t ) ) x k θ j ( x ( t ) ) x k d t
Thus, Equation (8) could be solved, and a string could be evolved by the iteration of the mean force until convergence.

2.3. Minimum-Free-Energy Path from Finite-Temperature String

(1) Initial trajectory.
The first step for constructing FTS is to generate an initial trajectory connecting two basins. It is very difficult for the polymorphic transformation of the molecular crystals due to the size effects of the system and the dimensionality of the OP space. In most cases, either it cannot be found, or an intermediate form or a form that is not interested is present, or it leads to an irreversible or inconvergent initial string [68,69]. This can also be caused by finite size effects, i.e., if the critical nuclei are of similar size or larger than the simulation box. If system size is not an issue, one way to improve the starting string is to run an ensemble of molecular dynamics (MD) simulations starting from the point corresponding to the highest potential of mean force and neighboring points, with velocities drawn from a Maxwell–Boltzmann distribution. We can then choose the trajectories that connect two crystal-form basins and restart the string method using images from it as the initial replicas, and the corresponding initial OPs θ i , m * ( 0 ) and θ C , m * ( 0 ) can be calculated using Equations (6) and (7) by an MD simulation.
(2) K-means clustering to obtain a smooth initial string.
In order to obtain a smooth initial string, a K-means clustering algorithm can be adopted. It consists of four steps:
Step 1. Restricting sampling. In each independent space around each given representative point, the samples can be taken according to the following harmonic functions:
ψ ( θ C d , b ) = k d 2 [ θ C , m d θ C , m d ( 0 ) ] 2 + k b 2 [ θ C , m b θ C , m b ( 0 ) ] 2
ψ ( θ C d , r ) = k d 2 [ θ C , m d θ C , m d ( 0 ) ] 2 + k r 2 [ θ C , m r θ C , m r ( 0 ) ] 2
where k d , k b , and k r are the spring coefficients. In this way, each simulation was localized in the independent space determined by each of the representative replicas.
Step 2. K-means clustering. Assumed that there is always a clustering center with the maximum weight in each independent space, and the OPs corresponding to the clustering center with the maximum weight are different in the different replicas. Obviously, such a center could be regarded as a locally important sampling point, which, to some extent, shows the evolution direction of the representative points along the string. Therefore, the application of the clustering algorithm can not only find the clustering center with the largest weight to solve the problem of distinguishing highly similar OPs in the adjacent replicas but also accelerate the convergence rate.
For the clustering algorithm, the dimensionally weighted Euclidean distance was used to measure the distance between sample points; after calculating the density and weight of all samples, the point with the highest weight can be selected in each independent sampling space. The weight wi is given by
ω i = ρ i × s i × 1 a i
where ρi, si, and ai represent the density of sample point i (i.e., the number of points i within a certain range), the distance between different clusters or datasets, and the distance between different sample points in the same cluster, respectively. They can be calculated by Equations (17) and (18) as follows:
ρ i = j = 1 n f [ d ω ( x i , x j ) 1 C n 2 i = 1 n j = i + 1 n d ω ( x i , x j ) ]
s i = { min ( d ω ( x i , x j ) ) , ρ j < ρ i max ( d ω ( x i , x j ) ) , ρ j ρ i
a i = 2 ρ i ( ρ i 1 ) i = 1 ρ i j = i + 1 ρ i d ω ( x i , x j )
where d ω ( x i , x j ) is the Euclidean distance, and f ( x ) = { 1 , x < 0 0 , x 0 . These adjacent clustering centers with the maximum weights can be connected to form a temporary string.
Step 3. Searching for successive strings by smooth function. In order to obtain a successive string in the collective coordinate space, the smooth function was defined as follows:
S * = m = 2 L 1 exp [ θ C , m * θ C , a v e r a g e ( m 1 , m + 1 ) * ]
where θ C , a v e r g e ( m 1 , m + 1 ) = ( θ C , m 1 + θ C , m + 1 ) / 2 . The smoothness of the temporary string was calculated by the Monte Carlo method. Firstly, the smoothing score of the current temporary string was calculated. Then, in the same space, another clustering center, where the weight difference from the first clustering center is small, was selected, and the smoothing score was recalculated. According to the new smoothing score, whether to accept the new change was determined according to the metropolis criterion. Finally, a string with the lowest smoothing score was obtained.
Step 4. Calculating  z C , m * . Based on the optimization string, OPs of the mth replica were calculated as the initial OPs to evolve string, represented by z C , m * . For comparison, the values of z C , m * were also directly calculated without K-means clustering, i.e., average-based sampling.
(3) Determining the MFEP using the SMCV method.
The MFEPs can be determined with the SMCV method by the following steps:
Step 1. Estimating  z F ( z C , m * )  and  M ( z C , m * ) . Based on the initial string with or without clustering, the free-energy gradient z F ( z C , m * ) and metric tensor M ( z C , m * ) can be estimated by using Equations (12) and (13), respectively.
Step 2. Calculating a new target OPs  z C , m ( new ) * . Following Maragliano [70], z C , m ( new ) * was calculated by the string evolution equation,
z C , m ( new ) * = z C , m * Δ τ M ( z C , m * ) F ( z C , m * )
Step 3. Reparameterizing. Interpolate a curve z(α) through z C , m ( new ) * by the b-spline fitting method and compute the new targets OPs (represented as z C , m + 1 * ) by interpolating to make the arc length at a constant between the consecutive replicas.
Step 4. Obtaining MFEP under convergence. Compute the difference between z C , m + 1 * and z C , m * . If it is small enough, the string is converged, and the MFEP is obtained by the converged target OPs; otherwise, go back to Step 1, and the iterative process is continued until it converges.
Step 5. Calculating local OPs and PMF. Based on the final converged string, the local OPs and PMF were calculated.

3. Results and Discussion

3.1. Crystallization of TKX-50 by DC External Electric Field

The experiment on the crystallization of TKX-50 from the supersaturated formic acid solution was carried out under the DC external electric field. The experimental setup is in Figure 1, and the molecular structure of TKX-50 and crystal data are shown in Figure S1 and Tables S1–S10.
Although the space group of the TKX-50 crystal was not affected by the external electric field, and the crystalline structure of TKX-50 has a monoclinic crystalline structure with space group P21/c, the crystal data and structural parameters are changed, obviously. In particular, the volume decreased, and the density increased by 0.1 g/cm3. For the explosive, the higher the density, the higher the detonation velocity and pressure become. Therefore, the introduction of the external electric field into the explosive system is beneficial for improving the performance of explosives. Moreover, from Tables S3 and S8, it can be seen that the bond length of C–N on the TKX-50 tetrazole ring is shortened by an average of 0.04 Å under the external electric field, suggesting that the external electric field enhances the stability of TKX-50.
The X-ray powder diffraction is shown in Figure 2. Without the external electric field, TKX-50 exhibits the characteristic diffraction peaks at 15.03, 15.46, 26.12, 27.00, and 31.15. Under the external electric field, the diffraction peaks undergo changes, appearing at positions 15.01, 16.27, 26.54, 27.05, and 30.89. In particular, compared with the diffraction peaks without the electric field, the non-diffraction characteristic peaks under the electric field are fewer, indicating that the crystalline quality and purity of TKX-50 are enhanced under the external electric field.

3.2. MD Simulation on Crystallization of TKX-50 without External Electric Field

3.2.1. Peaks in Pair Distribution Function

The peaks in the pair distribution functions for the reference structure of the TKX-50 crystal and the corresponding average peak locations and concentration parameters at 300 K with a cutoff of 10.0 Å are given in Table 1. The peak locations are very close to the experimental results. The smaller standard deviation and larger concentration parameter display the localization characteristics of the average centroid distance, bond orientation, and relative orientation distribution.
For the homogeneous nucleation of crystals of the ionic liquid [dmim+][Cl] from its supercooled liquid phase, two [dmim+][Cl] ion pairs show the variables relevant for the construction of OPs, which are based on the [dmim+] cations. The vector r joins the center of mass of the two cations. The vector normal to the plane of the imidazolium ring gives the absolute orientation q of each of the cations [71]. Similarly, according to the construction method of OPs by Santiso et al. [39], the vector r joins the center of mass of two tetrazolium anion rings, and the absolute orientation q was given by two anions. The bond orientation ϕ r ^ and the relative orientation ϕ q were also given by two anions; see “Section 4.2”. Thus, the local OPs of the TKX-50 crystal and the corresponding forms in the formic acid solution were constructed by the structural variables involving the bond distances, bond orientations, and relative orientations (see Figure 3). After obtaining the maximum-likelihood estimators, they were calculated (see Figure 4). Although three OPs of the TKX-50 crystal phase and the corresponding forms in the formic acid solution overlap slightly, the peaks exhibit significantly different values. Therefore, any of the OPs can serve as a good metric to detect the orders of the crystalline and the corresponding forms in the formic acid solution. In this work, we chose the bond orientation OPs and relative orientation OPs during all our simulations.

3.2.2. Convergence of FTS and K-Means Clustering

According to the trajectories from the simulation of 3 × 5 ns for 21 independent spaces, the average sampling and K-means clustering of OPs for bond orientation and relative orientation were performed, as shown in Figure 5. It was found that for ϕ C d b , the ln(Convergence) values are in the range of −5.5~−2.5 after about 20 iterations, indicating that there has been convergence. However, whether for the ln(Convergence) values from the average sampling or the K-means clustering, the changes are still uncertain and remain oscillating after 150 iterations. In this case, the simulation may cover adjacent spaces along the transition path, leading to an uneven initial string connecting adjacent points. For the non-smooth string, there are always some independent sampling spaces that are not perpendicular to the current string, leading to confusion in the evolution direction of the string.
As for ϕ C d r , compared to the ln(Convergence) values from the average sampling, those from the K-means clustering have decreased after 60 iterations. This shows that the application of the K-means clustering algorithm can not only quickly find the cluster center with the highest weight and solve the differentiation of highly similar OPs in adjacent replicas but also accelerate the convergence speed of strings. Generally, the K-means clustering algorithm combining dimensionally weighted Euclidean distance with density weight has a strong outlier interference ability, which can effectively avoid local optima and quickly converge to global optima to a certain extent. These are in agreement with our recent investigation [72].

3.2.3. Minimum-Free-Energy Path

PMF is the line integral of the restraint force along the path, and it shows the MFEP as the maximum-likelihood path in the space of collective variables along FTS [67,70].
The convergent FTS path from the TKX-50 crystal (isolated from the formic acid phase) to its supersaturated formic acid solution was obtained using the SMCV method. Figure 6I,II show the PMF curves as the function of the arclength along the converged FTS paths corresponding to the ϕ C d b as the collective variables by K-means clustering and without K-means clustering, respectively. (III) and (IV) correspond to the FTS path obtained from ϕ C d r without K-means clustering and with K-means clustering. The difference in PMF between the supersaturated formic acid solution and the transition state is small, with the value less than 4.0 kJ∙mol−1 in all the cases, while the difference in PMF between the TKX-50 crystal and the transition state is large, close to 50.0 kJ∙mol−1 from ϕ C d b in two cases and ϕ C d r as the collective variables by K-means clustering, and about 40.0 kJ∙mol−1 from ϕ C d r without K-means clustering. These results show that, at 300 K, the solubility of TKX-50 in formic acid is not high, while its saturated solution is prone to crystallization. This is in agreement with our experimental result (0.5 g TKX-50 is dissolved in 45 mL formic acid in this work). Note that due to PMF corresponding to the dimensions (3 × 3 × 3), instead of one-dimensional reaction coordinates, the values of PMF are higher than the real free-energy barrier. Moreover, the differences in PMF between the transition state and TKX-50 are close to each other, and those between the transition state and the supersaturated formic acid solution are also close to each other. This indicates that the difference is not related to the type of OPs or whether it is K-means clustering, as is in agreement with our previous investigation [56].
In order to gain further mechanistic insights into the crystallization of TKX-50 from the supersaturated formic acid solution, the local OPs of TKX-50 were shown for the key snapshots at the points of the PFM curve along the trajectory going from crystal to solution (see Figure 7). The molecules that can be regarded as existing in the crystal were labeled in blue with the values of ϕ C d b > ~140 or ϕ C d r > ~160, and the green molecules were regarded as TKX-50 in solution with the values of ϕ C d b < ~50 or ϕ C d r < ~40, whereas the wireframe molecules in red formed the interface between crystal and solution with the values in the range of 50~140 or 40~160 for s ϕ C d b and ϕ C d r , respectively.
From Figure 7, it can be seen that along the (a) → (b) → … → (l), the type of molecule gradually changes from the “all blue” initial conformation, the formation of molecular clusters of which “Scattered green areas are surrounded by blue areas” (interface induction), to gradually increasing “sporadic green” molecular clusters (i.e., local initiation, multinuclear asynchronous generation, and increase), and then to the “blue being engulfed by green” conformation. In most cases, OPs show a central region with the TKX-50 structure surrounded by molecules with intermediate values of OPs corresponding to the surface of the crystal. This indicates that similar to the polymorphic transformation [56,72], the crystallization from solution is also an interface-induced and locally initiated development process.

3.3. MD Simulation of Crystallization of TKX-50 under the External Electric Fields

The average peak locations and concentration parameters under the external electric fields with the strength of 51.40 × 108 V/m are shown in Table 2. Compared with the values without the external electric fields, both the peak values and the concentration parameters of the pair distribution functions for TKX-50 crystal change significantly, especially for them under the external electric field along the direction of the c-axis. As expected, due to two hydroxyamine groups having a positive charge and the dihydroxy-5,5’-tetrazolium group having a negative charge, they will undergo significant displacement under the external electric field, resulting in a significant change in the peak values and the concentration parameters, as shown in Figure 1c. Most of the directions of the positive and negative charges for hydroxylamine and dihydroxy-5,5’-tetrazolium groups are exactly along the direction of the c-axis of the crystal (see Figure 1h), so the changes are more significant along the c-axis. Moreover, in most cases, the values of ϕ q have decreased while those of ϕ r ^ α increased, showing that the π∙∙∙π stacking is more significant under the external electric field.
Figure 5b gives the convergence of ϕ C d r under an external electric field with a strength of 51.40 × 108 V/m by K-means clustering. Compared with those without the external electric field (see (IV)), the values of the ln(Convergence) of the relative orientation OPs are decreased greatly, and there is convergence after 60 iterations in all cases. This shows that the external electric field can accelerate the evolution of the string. In 2017, Jha et al. found that external electric fields can obviously increase the nucleation rate [59]. Koizumi et al. also observed the same phenomenon [73]. These results confirmed that, shown by the accelerated evolution of the string, the external electric field can accelerate the nucleation rate.
For the ϕ C d r as the collective variables, similar to the cases without an electric field, the difference in PMF between the supersaturated formic acid solution and the transition state is also small, with the value less than 3.5 kJ∙mol−1 under the external electric field. The difference in PMF between the TKX-50 crystal and the transition state is less than 40.0 kJ∙mol−1 under the external electric fields with a strength of 51.40 × 108 V/m along the a-, b-, and c-axes. These values are decreased in comparison with those in no field. Moreover, the fields parallel to the b- and c-axes affect the difference in PMF between the TKX-50 crystal and the transition state less than those parallel to the a-axis. The crystallization of TKX-50 from the formic acid solution is mainly closely related to the intermolecular interaction. From Figure 1f–h, along the a-axis, the strong H-bonding interactions are formed. They are strengthened by the increased dipole moments more significantly under the external electric fields, leading to a large value of the PMF between the TKX-50 crystal and the transition state along the a-axis. The value of the PMF between them is the smallest along the b-axis. From Figure 1f–h, the intermolecular π∙∙∙π stackings are formed along the b-axis. Due to the weaker π∙∙∙π stacking than H-bonding interaction, the more notable change will occur under the external electric field along the b-axis, leading to the low energy of the transition state and thus the smaller difference in PMF. Therefore, applying an external electric field along the b-axis direction has more practical value for achieving the crystallization of TKX-50 from solution.
In order to further verify the influence of the external electric field along the b-axis direction on the crystallization of TKX-50 from the supersaturated formic acid solution, the external electric fields were added with the field strengths of 0.514 × 108, 5.14 × 108, and 10.28 × 108~102.80 × 108 V/m with a step of 10.28 × 108 V/m. The values of the difference in PMF between the TKX-50 crystal and the transition state were calculated to be 51.3, 48.2, 37.6, 42.5, 40.1, 38.3, 36.6, 35.3, 39.7, 42.5, 46.8, and 32.5 kJ/mol. The corresponding 3D surface plots of the free-energy landscape obtained from θ C d and θ C b as the collective variables are shown in Figure 8. The influence of the external electric field on PMF has no regular pattern in values. For the former, a good linear result (R2 = 0.9882) was found, showing that the external electric fields (F) have an obvious influence on the polymorphic transformation from o-TNT to m-TNT. For the latter, this relationship was not discovered from the above values of the difference in PMF between the TKX-50 crystal and the transition state. As mentioned above, we only found that the values of the difference in PMF are decreased in comparison with those in no field. It is well-known that the movement of ions in an external electric field is complex, and in the process of the crystallization from solution, this movement pattern may be more complex, leading to irregular changes in the difference in PMF. However, one conclusion can be confirmed that the external electric fields can reduce the difference in PMF and thus accelerate the nucleation rate, as is also confirmed by the experimental result [59].

3.4. Detonation Performance

According to the densities of TKX-50 from the experiments (1.920 g/cm3, see Table S1 without field vs. 1.933 g/cm3, see Table S6 under external electric field), the calculated detonation velocity (VD) and detonation pressure (PD) were 9604 m/s and 42.49 GPa in no field, and 9650 m/s and 43.07 GPa with an external electric field, respectively. The results of VD and PD without field are consistent with those in Reference [74]. Compared with the values without the external electric fields, the values of VD and PD under the external electric field were increased by 0.48% and 1.37%, respectively. This indicates that the external electric field can enhance the detonation performance of TKX-50.
In summary, according to the experimental results, the structural parameters of the TKX-50 crystal were changed under the external electric field, leading to increased density, detonation velocity, and detonation pressure. From the MD simulation, the K-means clustering algorithm based on Euclidean distance and density weights can be used for enhanced sampling of the OPs in external electric-field-induced crystallization of ionic crystals from solution, which can improve the conventional FTS. The crystallization from solution is a process of surface-mediated nucleation. The external electric fields can improve the conventional FTS and decrease the difference in PMF between the crystal and the transition state. This is attributed to the significant changes in the OPs under the external electric field, which leads to the enhanced crystalline quality and purity of TKX-50 observed from the experiment. This suggests that the process of the crystallization of TKX-50 from solution can be controlled by changing the external electric fields to enhance the crystalline quality and purity.

4. Experiment and MD Simulation Details

4.1. Experiment

An electrochemical workstation was used to recrystallize TKX-50 under the external electric field. Dissolve TKX-50 (0.5 g, purity: 99.2%) in formic acid (45 mL, concentration: 88%) and prepare it as a supersaturated solution. The distance between the working electrode and the auxiliary electrode is 0.5 cm. The working electrode is a square platinum plate electrode with a side length of 1.0 cm, and a constant voltage of 2V was used.
The TKX-50 crystals with the dimensions of 0.16 mm × 0.09 mm × 0.08 mm obtained without the external electric field and those with the dimensions of 0.12 mm × 0.07 mm × 0.05 mm obtained under the external electric field were selected. They were placed on a Bruker D8 VENTURE single-crystal diffractometer. MoKα Ray (λ = 0.71073 Å) monochromated by graphite monochromator was as the X-ray source, with ω-scanning method within a certain scanning range θ.

4.2. MD Simulation Details

Based on our experimental structure (see Table S1), the TKX-50 crystal structure was constructed with the size of 6 × 6 × 6 (216 molecules). Based on the molecular structure and strategy of the OP construction for molecular crystal [39], a set of OPs for TKX-50 was built with the pair distribution function (see Figure 2).
The crystallization of TKX-50 from solution involves mainly the disruption of the cohesive forces in the crystal lattice. Due to the insignificant –NH3 rotation, the internal degrees of freedom of conformational transformation can be ignored. Based on molecular symmetry and point molecule representation [69], The vector r joins the center of mass of the two tetrazolium rings. The vector normal to the plane of the tetrazolium ring gives the absolute orientation q of each of the tetrazolium anions. The bond orientation ϕ r ^ is given by the angle formed by the vectors r and q1, and the relative orientation ϕ q is given by the angles between the vectors q1 and q2, as shown in Figure 2.
In order to construct the above OPs, a short (500 ps) MD simulation was first carried out for the structure of TKX-50 from the experimental crystal and its supersaturated formic acid solution at 300 K and 1 bar. After equilibration, the structures were annealed to zero temperature at constant volume and relaxed by energy minimization. Then, an MD simulation (2.5 ns) was carried out at 300 K and 1 bar. The trajectory was used to obtain maximum-likelihood estimators of parameters in Equations (6) and (7) in the case of the simulated box with the 3 × 3 × 3 grid.
A 1.0 ns MD simulation was carried out on the system with the 216 TKX-50 molecules at 300 K and 1 bar. The long-range electrostatics were handled using the particle mesh Ewald method [75] with a cutoff of 12 Å, and a Langevin thermostat with an oscillation frequency of 25 ps−1 was introduced for the temperature control. The Nosé–Hoover–Langevin piston with a damping time of 50 fs was used to control the pressure. The hydrogen bond lengths were constrained with the LINCS algorithm [76]. Periodic boundary conditions were applied in all directions. After equilibrium, annealing treatment was carried out at a constant volume to zero temperature, and energy minimization relaxation was carried out. Then, the temperature was raised to 300 K, and a 2.5 ns MD simulation was carried out at a constant volume with a time step of 0.5 fs. Based on the obtained trajectory, the parameters of Equations (1)–(7) were calculated using the maximum-likelihood estimation method, and the initial order parameters were determined.
The external electric fields with a strength of 51.40 × 108 V/m were added along the positive direction of the a-, b-, and c-axes. For the c-axis, the field strengths are 0.00, 0.514 × 108, 5.14 × 108, and 10.28 × 108~102.80 × 108 V/m with a step of 10.28 × 108 V/m were also applied. The MD simulation of the structural optimization was carried out by Materials Studio 5.0 software package with COMPASS force field.
For the MD simulation of FTS evolution, the ensemble remains at NPT (300 K and 1 bar). For each initial configuration, we choose the configuration that is closest to the new target OP from the previous step. To maintain the stability of the simulation, the center of the OP constraint is moved from the corresponding value in the previous step to a new value in the previous period of time. Average the constraint force for each step in the later period, and the initial velocity is randomly assigned from the Maxwell–Boltzmann distribution of the temperature at which it is located. The CHARMM22 force field was used for the melting of benzene, and the calculated results were consistent with experimental values [69]. Therefore, the CHARMM22 force field was selected. All the calculations were carried out by a modified version of the NAMD 2.14 software package [77] that included implementations of the OPs and the SMCV. All the simulations were performed in the NPT ensemble with P = 1 bar and T = 300 K.

4.3. Detonation Performance Calculations

The detonation velocity (VD) and detonation pressure (PD) can be evaluated by Kamlet approximation, as shown by Equation (22) [78]:
V D = 1.01 ( N 2 M ¯ Q D ) 1 4 ( 1 + 1.30 d ) P D = 1.558 ( N 2 M ¯ Q D ) 1 2 d 2
where N is the moles of gaseous detonation products per gram of explosive, M ¯ is the average molecular weight of gaseous products, and QD and d are the heat of detonation reaction and density of explosives, respectively.

5. Conclusions

In this work, an experimental investigation on the crystallization of TKX-50 from the supersaturated formic acid solution was carried out under the DC external electric field. To explain the experimental results, a theoretical study was presented by FTS with the local OPs as collective variables for ionic crystals. A K-means clustering algorithm based on the dimensionally weighted Euclidean distance and sample weights was used to smooth the FTS. The MFEP involving the crystallization of TKX-50 from solution was sketched by the SMCV method. The following conclusions were drawn from this study:
(1) From the experiment, the external electric field enhanced the density, detonation velocity, and detonation pressure of TKX-50.
(2) The K-means clustering algorithm based on Euclidean distance and density weights can be used for enhanced sampling of the OPs in the external electric-field-induced crystallization of the TKX-50 ionic crystal from solution, which can improve the conventional FTS.
(3) The crystallization of TKX-50 from the formic acid solution under the external electric field could be regarded as a process of surface-mediated nucleation.
(4) Due to the significant changes in OPs, the external electric fields accelerate the evolution of the string and decrease the difference in PMF between the crystal and the transition state. The process of the crystallization of TKX-50 from solution can be controlled by changing the external electric fields with enhanced crystalline quality and purity.
This work provides an effective way to explore the crystallization of ionic crystals from solution at the molecular level, and it is useful for realizing the control of the production process and the improvement in the properties of ionic crystal explosives by using external electric fields.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules29051159/s1, Figure S1: Molecular structure of TKX-50; Table S1: Crystal data and structural parameters of TKX-50 without the external electric field; Table S2: Atomic coordinates and equivalent temperature factor of TKX-50 without the external electric field; Table S3: Selected bond lengths of the TKX-50 without the external electric field; Table S4: Selected bond angles of the TKX-50 without the external electric field; Table S5: Selected dihedral angles of the TKX-50 without the external electric field; Table S6: Crystal data and structural parameters of TKX-50 under the external electric field; Table S7: Atomic coordinates and equivalent temperature factor of TKX-50 under the external electric field; Table S8: Selected bond lengths of the TKX-50 under the external electric field; Table S9: Selected bond angles of the TKX-50 under the external electric field; Table S10 Selected dihedral angles of the TKX-50 under the external electric field.

Author Contributions

F.R.: conceptualization, project management, simulation and calculation using FTS with OP as collective variables, data curation, and writing—original draft preparation. X.W. (Xiaolei Wang) and X.W. (Xiaojun Wang): crystallization experiment and data curation. Q.Z., L.C. and Z.Z.: simulation and figures. All authors have read and agreed to the published version of the manuscript.

Funding

The authors are grateful for the financial support from the Shanxi Province Natural Science Foundation of China (No. 201801D121067 and No. 202203021221111).

Institutional Review Board Statement

We allow the journal to review all the data, and we confirm the validity of the results. There are no financial relationships. This work was not published previously, and it has not been submitted to more than one journal. The submission is also not split up into several parts. No data have been fabricated or manipulated.

Informed Consent Statement

Not applicable. This study did not involve humans.

Data Availability Statement

The data related to this research can be accessed upon reasonable request via email.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Li, H.; Ren, D.; Cheng, X. The theoretical investigation of the β-crystobalite structure under the effect of electric field. Comp. Mater. Sci. 2015, 96, 306–311. [Google Scholar] [CrossRef]
  2. Simura, R.; Nakamura, K.; Uda, S. Change of melting temperature of non-doped and Mg-doped lithium niobate underan external electric field. J. Cryst. Growth 2008, 310, 3873–3877. [Google Scholar] [CrossRef]
  3. Roland, C.M.; Fragiadakis, D.; Bogoslovov, R.; Urban, S.; Dabrowski, R.; Tykarska, M.; Osiecka, N.; Czub, J. Volumetric, dielectric, calorimetric and X-ray studies of smectogenic 10PBO8 at atmospheric and elevated pressures. Liq. Cryst. Rev. 2012, 39, 993–1001. [Google Scholar] [CrossRef]
  4. Mukherjee, P.K.; Rahman, M. Electric-field induced isotropic to smectic-Cphase transition. J. Mol. Liq. 2014, 196, 204–207. [Google Scholar] [CrossRef]
  5. Mahmood, A.; Chan, M.S.Y.; Saleemi, A.S.; Guo, J.; Lee, S.L. Synergic effect: Temperature-assisted electric-field-induced supramolecular phase transitions at the liquid/solid Interface. Langmuir 2019, 35, 8031–8037. [Google Scholar] [CrossRef] [PubMed]
  6. Lu, N.; Zhang, P.; Zhang, Q.; Qiao, R.; He, Q.; Li, H.B.; Wang, Y.; Guo, J.; Zhang, D.; Duan, Z.; et al. Electric-field control of tri-statephase transformation with a selective dual-ion switch. Nature 2017, 546, 124–128. [Google Scholar] [CrossRef]
  7. Jeong, J.; Aetukuri, N.; Graf, T.; Schladt, T.D.; Samant, M.G.; Parkin, S.S.P. Suppression of metal-insulator transition in VO2 by electric field-induced oxygen vacancy formation. Science 2013, 339, 1402–1405. [Google Scholar] [CrossRef]
  8. Alexander, L.F.; Radacsi, N. Application of electric fields for controlling crystallization. CrystEngComm 2019, 21, 5014–5031. [Google Scholar] [CrossRef]
  9. Radacsi, N. Process Intensification in Crystallization: Submicron Particle Generation Using Alternative Energy Forms. Ph.D. Thesis, University of Debrecen, Hungary, Germany, 2012. Available online: https://www.academia.edu/10747704/Process_Intensification_in_Crystallization_Submicron_Particle_Generation_using_Alternative_Energy_Forms (accessed on 1 October 2012). [CrossRef]
  10. Flores-Hernández, E.; Stojanoff, V.; Arreguín-Espinosa, R.; Moreno, A.; Sánchez-Puig, N.J. An electrically assisted device for protein crystallization in a vapor-diffusion setup. J. Appl. Crystallogr. 2013, 46, 832–834. [Google Scholar] [CrossRef]
  11. Walter, T.K.; Ferreira, C.F.D.G.; Iulek, J.; Benelli, E.M. Use of Protein Thin Film Organized by External Electric Field as a Template for Protein Crystallization. ACS Omega 2018, 3, 8683–8690. [Google Scholar] [CrossRef]
  12. Martínez-Caballero, S.; Cuéllar-Cruz, M.; Demitri, N.; Polentarutti, M.; Rodríguez-Romero, A.; Rodríguez-Romero, A. Glucose lsomerase Polymorphs Obtained Using an Ad Hoc Protein Crystallization Temperature Device and a Growth Cell Applying an Electric Field. Cryst. Growth Des. 2016, 16, 1679–1686. [Google Scholar] [CrossRef]
  13. Rodríguez-Romero, A.; Esturau-Escofet, N.; Pareja-Rivera, C.; Moreno, A. Crystal growth of high-quality protein crystals under the presence of an alternant electric field in pulse-wave mode, and a strong magnetic field with radio frequency pulses characterized by x-ray diffraction. Crystals 2017, 7, 179. [Google Scholar] [CrossRef]
  14. Sazaki, G.; Moreno, A.; Nakajima, K. Novel coupling effects of the magnetic and electric fields on protein crystallization. J. Cryst. Growth 2004, 262, 499–502. [Google Scholar] [CrossRef]
  15. Koizumi, H.; Uda, S.; Fujiwara, K.; Okada, J.; Nozawa, J. Effect of an External Electric Field on the Kinetics of Dislocation-Free Growth of Tetragonal Hen Egg White Lysozyme Crystals. Crystals 2017, 7, 170. [Google Scholar] [CrossRef]
  16. Li, F.; Lakerveld, R. Influence of Alternating Electric Fields on Protein Crystallization in Microfluidic Devices with Patterned Electrodes in a Parallel-Plate Configuration. Cryst. Growth Des. 2017, 17, 3062–3070. [Google Scholar] [CrossRef]
  17. Rubin, E.; Owen, C.; Stojanoff, V. Crystallization under an External Electric Field: A Case Study of Glucose lsomerase. Crystals 2017, 7, 206. [Google Scholar] [CrossRef]
  18. Taleb, M.; Didierjean, C.; Jelsch, C.; Mangeot, J.P.; Capelle, B.; Aubry, A. Crystallization of proteins under an external electric field. J. Cryst. Growth 1999, 200, 575–582. [Google Scholar] [CrossRef]
  19. Nanev, C.N.; Penkova, A. Nucleation of lysozyme crystals under external electric andultrasonic fields. J. Cryst. Growth 2001, 232, 285–293. [Google Scholar] [CrossRef]
  20. Nanev, C.N.; Penkova, A. Nucleation and growth of lysozyme crystals under external electric field. Colloids. Surf. A 2002, 209, 139–145. [Google Scholar] [CrossRef]
  21. Mirkin, N.; Frontana-Uribe, B.A.; Rodríguez-Romero, A.; Hernández-Santoyo, A.; Moreno, A. The influence of an internal electric field upon protein crystallization using the gel-acupuncture method. Acta Crystallogr. Sect. D Biol. Crystallogr. 2003, 59, 1533–1538. [Google Scholar] [CrossRef]
  22. Profio, G.D.; Reijonen, M.T.; Caliandro, R.; Guagliardi, A.; Curcio, E.; Drioli, E. Insights into the polymorphism of glycine: Membrane crystallization in an electric field. Phys. Chem. Chem. Phys. 2013, 15, 9271. [Google Scholar] [CrossRef]
  23. Moreno, A.; Sazaki, G. The use of a new ad hoc growth cell with parallel electrodes for the nucleation control of lysozyme. J. Cryst. Growth 2004, 264, 438–444. [Google Scholar] [CrossRef]
  24. Hammadi, Z.; Astier, J.P.; Morin, R.; Veesler, S. Protein Crystallization Induced by a Localized Voltage. Cryst. Growth. Des. 2007, 7, 1472–1475. [Google Scholar] [CrossRef]
  25. Nieto-Mendoza, E.; Frontana-Uribe, B.A.; Sazaki, G.; Moreno, A. Investigations on electromigration phenomena for protein crystallization using crystal growth cells with multiple electrodes: Effect of the potential control. J. Cryst. Growth 2005, 275, 1437–1446. [Google Scholar] [CrossRef]
  26. Chen, C.-J.; Zhu, W.-B.; Chao, J.-H.; Shang, A.; Lee, Y.G.; Liu, R.-J.; Yin, S.S.; Mark, D.; Hoffman, R.C. Study of thermal and spatial dependent electric field-induced phase transition in relaxor ferroelectric crystals using Raman spectroscopy. J. Alloys Compd. 2019, 804, 35–41. [Google Scholar] [CrossRef]
  27. Chen, Y.-L.; Wang, Z.-W.; Chen, S.; Ren, H.; Li, B.-W.; Yan, W.-S.; Zhang, G.-B.; Jiang, J.; Zou, C.-W. Electric-field Control of Li-Doping Induced Phase Transition in VO2 Film with Crystal Facet-Dependence. Nano. Energy 2018, 51, 300–307. [Google Scholar] [CrossRef]
  28. Okamura, Y.; Kagawa, F.; Seki, S.; Tokura, Y. Transition to and from the skyrmion lattice phase by electric fields in a magnetoelectric compound. Nat. Commun. 2016, 7, 12669. [Google Scholar] [CrossRef] [PubMed]
  29. Zaragoza, A.; Espinosa, J.; Ramos, R.; Cobos, J.A.; Aragones, J.L.; Vega, C.; Sanz, E.; Ramirez, J.; Valeriani, C. Phase boundaries, nucleation rates and speed of crystal growth of the water-to-ice transition under an electric field: A simulation study. J. Phy-Condens. Mat. 2018, 30, 174002. [Google Scholar] [CrossRef]
  30. Deng, J.-K.; Chang, Z.-Y.; Zhao, T.; Ding, X.-D.; Sun, J.; Liu, J.Z. Electric Field Induced Reversible Phase Transition in Li Doped Phosphorene: Shape Memory Effect and Superelasticity. J. Am. Chem. Soc. 2016, 138, 4772–4778. [Google Scholar] [CrossRef]
  31. Sawahata, H.; Yamaguchi, N.Y.; Fumiyuki, I. Electric-field-induced Z2 topological phase transition in strained single bilayer Bi(111). Appl. Phys. Express 2019, 12, 075009. [Google Scholar] [CrossRef]
  32. Sawahata, H.; Yamaguchi, N.Y.; Kotaka, H.; Fumiyuki, I. First-principles study of electric-field-induced topological phase transition in one-bilayer Bi(111). Jpn. J. Appl. Phys. 2018, 57, 030309. [Google Scholar] [CrossRef]
  33. Collins, J.L.; Tadich, A.; Wu, W.; Gomes, L.C.; Rodrigues, J.N.B.; Liu, C.; Hellerstedt, J.; Ryu, H.; Tang, S.; Mo, S.-K.; et al. Electric-field-tuned topological phase transition in ultrathin Na3Bi. Nature 2018, 564, 390–394. [Google Scholar] [CrossRef] [PubMed]
  34. Matvija, P.; Rozbořil, F.; Sobotík, P.; Ošťádal, I.; Pieczyrak, B.; Jurczyszyn, L.; Kocán, P. Electric-field-controlled phase transition in a 2D molecular layer. Sci. Rep. 2017, 7, 7357. [Google Scholar] [CrossRef] [PubMed]
  35. Jeong, D.S.; Thomas, R.; Katiyar, R.S.; Scott, J.F.; Kohlstedt, H.; Petraru, A.; Hwang, C.S. Emerging memories: Resistive switching mechanisms and current status. Rep. Prog. Phys. 2012, 28, 076502. [Google Scholar] [CrossRef] [PubMed]
  36. Mallah, T.; Cavallini, M. Surfaces, thin films and patterning of spin crossover compounds. C. R. Chim. 2018, 21, 1270–1286. [Google Scholar] [CrossRef]
  37. Dickson, A.; Warmflash, A.; Dinner, A.R. Nonequilibrium umbrella sampling in spaces of many order parameters. J. Chem. Phys. 2009, 130, 074104. [Google Scholar] [CrossRef] [PubMed]
  38. Faradjian, A.K.; Elber, R. Computing time scales from reaction coordinates by milestoning. J. Chem. Phys. 2004, 120, 10880–10889. [Google Scholar] [CrossRef]
  39. Santiso, E.E.; Trout, B.L. A general set of order parameters for molecular crystals. J. Chem. Phys. 2011, 134, 064109. [Google Scholar] [CrossRef] [PubMed]
  40. Carter, E.; Ciccotti, G.; Hynes, J.; Kapral, R. Constrained reaction coordinate dynamics for the simulation of rare events. Chem. Phys. Lett. 1989, 156, 472–477. [Google Scholar] [CrossRef]
  41. Anwar, J.; Zahn, D. Uncovering molecular processes in crystal nucleation and growth by using molecular simulation. Angew. Chem. Int. Ed. 2011, 50, 1996–2013. [Google Scholar] [CrossRef]
  42. Zinovjev, K.; Tuñón, I. Adaptive finite temperature string method in collective variables. J. Phys. Chem. A 2017, 121, 9764–9772. [Google Scholar] [CrossRef]
  43. Dickson, B.M.; Huang, H.; Post, C.B. Unrestrained computation of free energy along a path. J. Phys. Chem. B 2012, 116, 11046–11055. [Google Scholar] [CrossRef]
  44. Díaz, L.G.; Ensing, B. Path finding on high-dimensional free energy landscapes. Phys. Rev. Lett. 2012, 109, 020601. [Google Scholar] [CrossRef]
  45. Maragliano, L.; Roux, B.; Vanden-Eijnden, E. Comparison between mean forces and swarms-of-trajectories string methods. J. Chem. Theory Comput. 2014, 10, 524–533. [Google Scholar] [CrossRef]
  46. Song, H.-D.; Zhu, F.-Q. Finite temperature string method with umbrella sampling: Application on a side chain flipping in Mhp1 transporter. J. Phys. Chem. B 2017, 121, 3376–3386. [Google Scholar] [CrossRef]
  47. Cao, L.R.; Lv, C.; Yang, W. Hidden conformation events in DNA Base extrusions: A generalized-ensemble path optimization and equilibrium simulation study. J. Chem. Theory. Comput. 2013, 9, 3756–3768. [Google Scholar] [CrossRef]
  48. Weinan, E.; Ren, W.; Vanden-Eijnden, E. String method for the study of rare events. Phys. Rev. B Condens. Matter Mater. Phys. 2002, 66, 052301. [Google Scholar]
  49. Vanden-Eijnden, E.; Venturoli, M. Revisiting the finite temperature string method for thecalculation of reaction tubes and free energies. J. Chem. Phys. 2009, 130, 194103. [Google Scholar] [CrossRef] [PubMed]
  50. Bellucci, M.A.; Trout, B.L. Bezier curve string method for the study of rare events in complex chemical systems. J. Chem. Phys. 2014, 141, 074110. [Google Scholar] [CrossRef] [PubMed]
  51. Shah, M.; Santiso, E.E.; Trout, B.L. Computer simulations of homogeneous nucleation of benzene from the melt. J. Phys. Chem. B 2011, 115, 10400–10412. [Google Scholar] [CrossRef]
  52. Gobbo, G.; Bellucci, M.A.; Tribello, G.A.; Ciccotti, G.; Trout, B.L. Nucleation of molecular crystals driven by relative information entropy. J. Chem. Theory Comput. 2018, 14, 959–972. [Google Scholar] [CrossRef]
  53. Mac, Q.J. Some methods for classification and analysis of multivariate observations. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability; Berkeley Symposium on Mathematical Statistics and Probability: Berkeley, CA, USA, 1967. [Google Scholar]
  54. Abo-Elnaga, Y.; Nasr, S. K-means cluster interactive algorithm-based evolutionary approach for solving bilevel multi-objective programming problems. Alex. Eng. J. 2021, 61, 811–827. [Google Scholar] [CrossRef]
  55. Rong, H.; Ramirez-Serrano, A.; Guan, L.; Gao, Y. lmage object extraction based on semantic detection and improved K-Means agorithm. IEEE Access 2020, 8, 171129–171139. [Google Scholar] [CrossRef]
  56. Ren, F.-D.; Liu, Y.-Z.; Ding, K.-W.; Chang, L.-L.; Cao, D.-L.; Liu, S.-B. Finite temperature string by K-means clustering sampling with order parameter as collective variables for molecular crystal: Application to polymorphic transformation between β-CL-20 and ε-CL-20. Phys. Chem. Chem. Phys. 2024, 26, 3500–3515. [Google Scholar] [CrossRef] [PubMed]
  57. Gao, Z.; Rohani, S.; Gong, J.; Wang, J. Recent Developments in the Crystallization Process: Toward the Pharmaceutical Industry. Engineering 2017, 3, 343–353. [Google Scholar] [CrossRef]
  58. Aber, J.E.; Arnold, S.; Garetz, B.A.; Myerson, A.S. Strong dc Electric Field Applied to Supersaturated Aqueous Glycine Solution Induces Nucleation of the γ Polymorph. Phys. Rev. Lett. 2005, 94, 145503. [Google Scholar] [CrossRef] [PubMed]
  59. Jha, P.K.; Sadot, M.; Vino, S.A.; Jury, V.; Curet-Ploquin, S.; Rouaud, O.; Havet, M.; Le-Bail, A. A review on effect of DC voltage on crystallization process in food systems. Innov. Food. Sci. Emerg. 2017, 42, 204–219. [Google Scholar] [CrossRef]
  60. Panda, S. The Internet of Things: Breakthroughs in Research and Practice; IGI Global: Pennsylvania, PA, USA, 2017. [Google Scholar]
  61. Myerson, A.S.; Ginde, R. Industrial Crystallization of Melts; CRC Press: Boca Raton, FL, USA, 2004; pp. 183–240. [Google Scholar]
  62. Chen, F.; Zhou, T.; Li, J.; Wang, X.; Cao, D.; Wang, J.; Yang, Z. Crystal morphology of dihydroxylammonium 5,5′-bistetrazole-1,1′-diolate (TKX-50) under solvents system with different polarity using molecular dynamics. Comp. Mater. Sci. 2019, 168, 48–57. [Google Scholar] [CrossRef]
  63. Xiong, S.-L.; Chen, S.-S.; Jin, S.-H.; Li, L.-J. Additives effects on crystal morphology of dihydroxylammonium 5,5’-bistetrazole-1,1’-diolate by molecular dynamics simulations. J. Energ. Mater. 2016, 34, 384–394. [Google Scholar] [CrossRef]
  64. Xu, X.; Chen, D.; Li, H.; Xu, R.; Zhao, H. Crystal Morphology Modification of 5, 5′-Bisthiazole-1, 1′-dioxyhydroxyammonium Salt. ChemistrySelect 2020, 5, 1919–1924. [Google Scholar] [CrossRef]
  65. Xiao, L.; Guo, S.; Su, H.; Gou, B.; Liu, Q.; Hao, G.; Hu, Y.; Wang, X.; Jiang, W. Preparation and characteristics of a novel PETN/TKX-50 co-crystal by a solvent/non-solvent method. RSC Adv. 2019, 9, 9204–9210. [Google Scholar] [CrossRef] [PubMed]
  66. Wang, X.-J.; Su, Q.; Chen, S.-S. Synthesis of cumulative nitrogen rich compound of dihydroxylammonium 5,5-bistetrazole -1,1-diolate (TKX-50). Initiators Pyrotech. 2014, 3, 38–41. [Google Scholar]
  67. Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String method in collective variables: Minimum free energypaths and isocommittor surfaces. J. Chem. Phys. 2006, 125, 024106. [Google Scholar] [CrossRef] [PubMed]
  68. Beckham, G.T.; Peters, B.; Starbuck, C.; Variankaval, N.; Trout, B.L. Surface-mediated nucleation in the solid-state polymorph transformation of terephthalic acid. J. Am. Chem. Soc. 2007, 129, 4714–4723. [Google Scholar] [CrossRef] [PubMed]
  69. Santiso, E.E.; Trout, B.L. A general method for molecular modeling of nucleation from the melt. J. Chem. Phys. 2015, 143, 174109. [Google Scholar] [CrossRef]
  70. Maragliano, L.; Vanden-Eijnden, E. On-the-fly string method for minimum free energy paths calculation. Chem. Phys. Lett. 2007, 446, 182–190. [Google Scholar] [CrossRef]
  71. He, X.; Shen, Y.; Hung, F.R.; Santiso, E.E. Molecular simulation of homogeneous nucleation of crystals of an ionic liquid from the melt. J. Chem. Phys. 2015, 143, 124506. [Google Scholar] [CrossRef]
  72. Beckham, G.T.; Peters, B.; Trout, B.L. Evidence for a size dependent nucleation mechanism in solid state polymorph transformations. J. Phys. Chem. B 2008, 112, 7460–7466. [Google Scholar] [CrossRef]
  73. Koizumi, H.; Fujiwara, K.; Uda, S. Control of Nucleation Rate for Tetragonal Hen-Egg White Lysozyme Crystals by Application of an Electric Field withVariable Frequencies. Cryst. Growth Des. 2009, 9, 2420–2424. [Google Scholar] [CrossRef]
  74. Zhu, Z.-S.; Jiang, Z.-M.; Wang, P.-C.; Lu, M.; Shu, Q.; Yu, X.-H. Synthesis and Properties of Dihydroxylammonium 5,5’-Bistetrazole-1,1’-diolate. Chin. J. Energ. Mater. 2014, 3, 332–336. [Google Scholar]
  75. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. [Google Scholar] [CrossRef]
  76. Hess, B.; Bekker, H.; Berendsen, H.J.C.; Fraaije, J. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463–1472. [Google Scholar] [CrossRef]
  77. Phillips, J.C.; Braun, R.; Wang, W.; Gumbart, J.C.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R.D.; Kale, L.V.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. [Google Scholar] [CrossRef] [PubMed]
  78. Kamlet, M.J.; Jacobs, S.J. Chemistry of detonations. I. A simple method for calculating detonation properties of C-H-N-O explosives. J. Chem. Phys. 1968, 48, 23–35. [Google Scholar] [CrossRef]
Figure 1. Experimental setup of DC external electric field and crystallization of TKX-50. (a): Experimental setup; (b): Schematic diagram of anion motion in an external electric field; (c):Molecular structure of TKX-50 obtained from experiment; (d): Molecular stacking in crystals; (e): Perspective of molecular stacking from the BC plane direction of the crystal (fh): The directions of the external electric fields along the three axes. Note: In (fh), a, b, and c mean the a-axis, b-axis, and c-axis, respectively.
Figure 1. Experimental setup of DC external electric field and crystallization of TKX-50. (a): Experimental setup; (b): Schematic diagram of anion motion in an external electric field; (c):Molecular structure of TKX-50 obtained from experiment; (d): Molecular stacking in crystals; (e): Perspective of molecular stacking from the BC plane direction of the crystal (fh): The directions of the external electric fields along the three axes. Note: In (fh), a, b, and c mean the a-axis, b-axis, and c-axis, respectively.
Molecules 29 01159 g001
Figure 2. XRD powder diffraction of TKX-50. (a): XRD under external electric field; (b): XRD without external electric field.
Figure 2. XRD powder diffraction of TKX-50. (a): XRD under external electric field; (b): XRD without external electric field.
Molecules 29 01159 g002
Figure 3. Construction of OPs for TKX-50. Red, blue, gray, and white represent O, N, C, and H atoms, respectively. The vector r joins the center of mass of the two molecules. The direction of the axis across the center of the ring can be used as an approximate measure of the absolute orientation (q1 or q2). The bond orientation ϕ r ^ is defined as the projection of r ^ i j onto q1, and the relative orientation ϕ q shows the rotation of q1 onto q2.
Figure 3. Construction of OPs for TKX-50. Red, blue, gray, and white represent O, N, C, and H atoms, respectively. The vector r joins the center of mass of the two molecules. The direction of the axis across the center of the ring can be used as an approximate measure of the absolute orientation (q1 or q2). The bond orientation ϕ r ^ is defined as the projection of r ^ i j onto q1, and the relative orientation ϕ q shows the rotation of q1 onto q2.
Molecules 29 01159 g003
Figure 4. The distribution of (a) distance, (b) bond orientation, and (c) relative orientation OPs for TKX-50 crystal and the corresponding forms in the formic acid solution. These distributions were obtained by considering 2.5 ns MD simulations in the NPT ensemble from the averaged values within the divided cells with the 3 × 3 × 3 grid for each of the systems with 216 TKX-50 molecules.
Figure 4. The distribution of (a) distance, (b) bond orientation, and (c) relative orientation OPs for TKX-50 crystal and the corresponding forms in the formic acid solution. These distributions were obtained by considering 2.5 ns MD simulations in the NPT ensemble from the averaged values within the divided cells with the 3 × 3 × 3 grid for each of the systems with 216 TKX-50 molecules.
Molecules 29 01159 g004
Figure 5. Convergence of collective variables during the evolution of the string. (I) and (II) correspond to ϕ C d b with K-means and without K-means clustering, respectively; (III) and (IV) correspond to ϕ C d r without K-means clustering and with K-means clustering, respectively. (V), (VI), and (VII) correspond to ϕ C d r with K-means clustering for the convergence of collective variables under an external electric field with the strength of 51.40 × 108 V/m along the a-, b-, and c-axes, respectively. (a) Convergence of collective variables without external electric field; (b) convergence of collective variables under external electric field.
Figure 5. Convergence of collective variables during the evolution of the string. (I) and (II) correspond to ϕ C d b with K-means and without K-means clustering, respectively; (III) and (IV) correspond to ϕ C d r without K-means clustering and with K-means clustering, respectively. (V), (VI), and (VII) correspond to ϕ C d r with K-means clustering for the convergence of collective variables under an external electric field with the strength of 51.40 × 108 V/m along the a-, b-, and c-axes, respectively. (a) Convergence of collective variables without external electric field; (b) convergence of collective variables under external electric field.
Molecules 29 01159 g005
Figure 6. PMF as a function of arclength along the FTS path for the converged string obtained from the SMCV. The initial point at arclength zero is the supersaturated formic acid solution of TKX-50, and the endpoint is TKX-50 crystal. (I,II) correspond to the FTS path obtained from ϕ C d b with K-means and without K-means clustering; (III,IV) correspond to the FTS path obtained from ϕ C d r without K-means clustering and with K-means clustering; (VVII) correspond to the FTS path obtained from ϕ C d r with K-means clustering for the convergence of collective variables under external electric field with the strength of 51.40 × 108 V/m along the a-, b-, and c-axes, respectively. In every image, B means the transition state, and A and C mean the state similar to the solution and the crystal, respectively.
Figure 6. PMF as a function of arclength along the FTS path for the converged string obtained from the SMCV. The initial point at arclength zero is the supersaturated formic acid solution of TKX-50, and the endpoint is TKX-50 crystal. (I,II) correspond to the FTS path obtained from ϕ C d b with K-means and without K-means clustering; (III,IV) correspond to the FTS path obtained from ϕ C d r without K-means clustering and with K-means clustering; (VVII) correspond to the FTS path obtained from ϕ C d r with K-means clustering for the convergence of collective variables under external electric field with the strength of 51.40 × 108 V/m along the a-, b-, and c-axes, respectively. In every image, B means the transition state, and A and C mean the state similar to the solution and the crystal, respectively.
Molecules 29 01159 g006aMolecules 29 01159 g006b
Figure 7. Changes in local order parameters on the FTS path shown by the key snapshots. Along the (al), the type of molecule gradually changes from the “all blue” initial conformation, the formation of molecular clusters of which “Scattered green areas are surrounded by blue areas”, to gradually increasing “sporadic green” molecular clusters, and then to the “blue being engulfed by green” conformation. (al) mean the several key snapshots from crystal to supersaturated formic acid solution of TKX-50 (the molecules of formic acid were deleted).
Figure 7. Changes in local order parameters on the FTS path shown by the key snapshots. Along the (al), the type of molecule gradually changes from the “all blue” initial conformation, the formation of molecular clusters of which “Scattered green areas are surrounded by blue areas”, to gradually increasing “sporadic green” molecular clusters, and then to the “blue being engulfed by green” conformation. (al) mean the several key snapshots from crystal to supersaturated formic acid solution of TKX-50 (the molecules of formic acid were deleted).
Molecules 29 01159 g007
Figure 8. Three-dimensional surface plots of the free-energy landscape obtained from θ C d and θ C b as the collective variables for the crystallization of TKX-50 from the supersaturated formic acid solution under the external electric fields along the b-axis: (a) 10.28 × 108 V/m; (b) 20.56 × 108 V/m; (c) 30.84 × 108 V/m; (d) 41.12 × 108 V/m; (e) 51.40 × 108 V/m; (f) 61.68 × 108 V/m; (g) 71.96 × 108 V/m; (h) 82.24 × 108 V/m.
Figure 8. Three-dimensional surface plots of the free-energy landscape obtained from θ C d and θ C b as the collective variables for the crystallization of TKX-50 from the supersaturated formic acid solution under the external electric fields along the b-axis: (a) 10.28 × 108 V/m; (b) 20.56 × 108 V/m; (c) 30.84 × 108 V/m; (d) 41.12 × 108 V/m; (e) 51.40 × 108 V/m; (f) 61.68 × 108 V/m; (g) 71.96 × 108 V/m; (h) 82.24 × 108 V/m.
Molecules 29 01159 g008
Table 1. Average peak locations and concentration parameters for TKX-50 crystal at 300 K a.
Table 1. Average peak locations and concentration parameters for TKX-50 crystal at 300 K a.
r (Å)1/σ2−1) ϕ r ^ α (°) η r ^ α ϕ q (°) η q α
5.53 (5.46)28.5611.53 (12.68)31.2575.31 (69.68)11.26
7.56 (7.62)5.1969.82 (67.53)28.6243.26 (44.52)23.18
8.69 (8.43)16.4078.93 (83.26)10.6314.83 (16.83)17.53
a The values in parentheses are the experimental results.
Table 2. Average peak locations and concentration parameters for the TKX-50 crystal at 300 K under the external electric field with the strength of 51.40 × 108 V/m.
Table 2. Average peak locations and concentration parameters for the TKX-50 crystal at 300 K under the external electric field with the strength of 51.40 × 108 V/m.
r (Å)1/σ2−1) ϕ r ^ α (°) η r ^ α ϕ q (°) η q α
5.62 a (5.50) b
6.03 c
18.69 (26.50)
15.83
42.82 (61.29)
68.66
32.83 (18.39)
30.26
25.69 (27.57)
19.68
11.53 (12.63)
9.10
7.56 (7.16)
6.98
29.02 (10.31)
6.28
69.82 (66.53)
70.15
25.72 (26.17)
36.01
32.55 (11.53)
28.62
20.12 (16.93)
22.26
8.69 (8.53)
8.55
15.62 (23.19)
21.51
78.18 (82.31)
80.66
15.70 (18.63)
19.25
16.76 (15.92)
17.39
23.17 (9.13)
16.52
a The values of the average peak locations and concentration parameters corresponding to the positive direction of the a-axis. b The values of the average peak locations and concentration parameters (in parentheses) corresponding to the positive direction of the b-axis. c The values of the average peak locations and concentration parameters (in italics) corresponding to the positive direction of the c-axis.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ren, F.; Wang, X.; Zhang, Q.; Wang, X.; Chang, L.; Zhang, Z. Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals. Molecules 2024, 29, 1159. https://doi.org/10.3390/molecules29051159

AMA Style

Ren F, Wang X, Zhang Q, Wang X, Chang L, Zhang Z. Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals. Molecules. 2024; 29(5):1159. https://doi.org/10.3390/molecules29051159

Chicago/Turabian Style

Ren, Fude, Xiaolei Wang, Qing Zhang, Xiaojun Wang, Lingling Chang, and Zhiteng Zhang. 2024. "Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals" Molecules 29, no. 5: 1159. https://doi.org/10.3390/molecules29051159

APA Style

Ren, F., Wang, X., Zhang, Q., Wang, X., Chang, L., & Zhang, Z. (2024). Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals. Molecules, 29(5), 1159. https://doi.org/10.3390/molecules29051159

Article Metrics

Back to TopTop