The thermal decomposition of condensed CL-20 was investigated using reactive force field molecular dynamics (ReaxFF MD) simulations of a super cell containing 128 CL-20 molecules at 800–3000 K. The VARxMD code previously developed by our group is used for detailed reaction analysis. Various intermediates and comprehensive reaction pathways in the thermal decomposition of CL-20 were obtained. Nitrogen oxides are the major initial decomposition products, generated in a sequence of NO₂, NO₃, NO, and N₂O. NO₂ is the most abundant primary product and is gradually consumed in subsequent secondary reactions to form other nitrogen oxides. NO₃ is the second most abundant intermediate in the early stages of CL-20 thermolysis. However, it is unstable and quickly decomposes at high temperatures, while other nitrogen oxides remain. At all temperatures, the unimolecular pathways of N―NO₂ bond cleavage and ring-opening C―N bond scission dominate the initial decomposition of condensed CL-20. The cleavage of the N―NO₂ bond is greatly enhanced at high temperatures, but scission of the C―N bond is not as favorable. A bimolecular pathway of oxygen-abstraction by NO₂ to generate NO₃ is observed in the initial decomposition steps of CL-20, which should be considered as one of the major pathways for CL-20 decomposition at low temperatures. After the initiation of CL-20 decomposition, fragments with different ring structures are formed from a series of bond-breaking reactions. Analysis of the ring structure evolution indicates that the pyrazine derivatives of fused tricycles and bicycles are early intermediates in the decomposition process, which further decompose to single ring pyrazine. Pyrazine is the most stable ring structure obtained in the simulations of CL-20 thermolysis, supporting the proposed existence of pyrazine in Py-GC/MS experiments. The single imidazole ring is unstable and decomposes quickly in the early stages of CL-20 thermolysis. Many C₄ and C₂ intermediates are observed after the initial fragmentation, but eventually convert into stable products. The distribution of the final products (N₂, H₂O, CO₂, and H₂) obtained in ReaxFF MD simulation of CL-20 thermolysis at 3000 K quantitatively agrees with the results of the CL-20 detonation experiment. The comprehensive understanding of CL-20 thermolysis obtained through this study suggests that ReaxFF MD simulation, combined with the reaction analysis capability of VARxMD, would be a promising method for obtaining deeper insight into the complex chemistry of energetic materials exposed to thermal stimuli.

Cu has been widely used as a substrate material for graphene growth. To understand the atomistic mechanism of growth, an efficient and accurate method for describing Cu-C interactions is necessary, which is the prerequisite of any possible large-scale molecular simulation studies. The semi-empirical density-functional tight-binding (DFTB) method has a solid basis from the density functional theory (DFT) and is believed to be a good tool for achieving a balance between efficiency and accuracy. However, existing DFTB parameters cannot provide a reasonable description of the Cu surface structure. At the same time, DFTB parameters for Cu-C interactions are not available. Therefore, it is highly desirable to develop a set of DFTB parameters that can describe the Cu-C system, especially for surface reactions. In this study, a parametrization for Cu-C systems within the self-consistent-charge DFTB (SCC-DFTB) framework is performed. One-center parameters, including on-site energy, Hubbard, and spin parameters, are obtained from DFT calculations on free atoms. Two-center parameters can be calculated based on atomic wavefunctions. The remaining repulsive potential is obtained as the best compromise to describe different kinds of systems. Test calculations on Cu surfaces and Cu-or C atom-adsorbed Cu surfaces indicate that the obtained parameters can generate reasonable geometric structures and energetics. Based on this parameter set, carbon dimerization on the Cu(111) surface has been investigated via molecular dynamics simulations. Since they are the feeding species for graphene growth, it is important to understand how carbon dimers are formed on the Cu surface. It is difficult to observe carbon dimerization in brute-force MD simulations even at high temperatures, because of the surface structure distortion. To study the dimerization mechanism, metadynamics simulations are performed. Our simulations suggest that carbon atoms will rotate around the bridging Cu atom after a bridging metal structure is formed, which eventually leads to the dimer formation. The free energy barrier for dimerization at 1300 K is about 0.9 eV. The results presented here provide useful insights for understanding graphene growth.

Because of broad potential applications in sensing, drug delivery, and molecular motors, two-dimensional (2D), flexible, responsive Janus materials have attracted considerable interest recently in many fields. Unfortunately, the molecular-level responsive deformation of these 2D Janus nanomaterials is still not clearly understood. Hence, investigating the influence factor and responsiveness of the deformation of the 2D flexible responsive Janus nanomaterials should be helpful to deepen our understanding of the deformation mechanism and may provide valuable information in the design and synthesis of novel functional 2D Janus nanomaterials. Therefore, a mesoscopic simulation method, dissipative particle dynamics simulation, based on coarse-grained models, is employed in this work to systematically investigate the effect of the chain length difference between grafted polymers within two compartments of each individual Janus nanosheet and the effect of solvent selectivity difference of these two compartments on the deformation of the polymer-grafted Janus nanosheet. Although the coarse-grained model within this simulation is relatively crude, it is still valid to provide a qualitative image of the deformation of the polymer-grafted Janus nanosheet. Furthermore, we find two basic principles: (1) with increasing length difference between grafted polymers on the two opposite surfaces, the nanosheet will bear an entropy-driven deformation with increasing curvature; (2) the solvent will preferentially wet the polymer layer with better compatibility, and such a swelling effect may also provide a driving force for the deformation process. Owing to the interplay of conformational entropy and mixing enthalpy, the equilibrium structures of the polymer-grafted Janus nanosheet result in several interesting structures, such as a tube-like structure with a hydrophobic outer surface, an envelope-like structure, and a bowl-like structure, with tuning of the chain length and solvent compatibility of grafted polymers. Additionally, an unusually tube-like structure with a hydrophobic outer surface has been observed for a relatively weak solvent selectivity, which may provide us a novel method to transfer materials into the incompatible environment and therefore has potential applications in many areas, such as controllable drug delivery and release, and industrial and medical detection. Our theoretical results first provide a fundamental insight into the controllable deformation of the flexible Janus nanosheet, which can then help in the design and synthesis of novel Janus nanodevices for potential applications in pharmaceuticals and biomedicine. Bearing the limited of the computational capabilities, our model Janus nanosheets are relatively small, which are not direct mappings from real system. We hope that a systematic simulation study on this topic would be possible soon with the rapid developments in computer technology and simulation methods, and this would provide an exhaustive and universal methodology to guide experimental studies and applications.

Selective molecular permeation through two-dimensional nanopores is of great importance for nanoporous graphene membranes. In this study, we investigate the selective permeation characteristics of gas molecules through a nitrogen-and hydrogen-modified graphene nanopore using molecular dynamics simulations. We reveal the mechanisms of selective molecular permeation from the aspects of molecular size and structure, pore configuration, and interactions between gas molecules and graphene. The results show that the permeances of different molecules are different, and the following order is observed in our study: H₂O > H₂S > CO₂ > N₂ > CH₄. Molecular permeance is related to the molecular size, mass, and molecular density on the graphene surface. The molecular permeation rate is inversely proportional to the molecular mass based on gas kinetic theory, while the molecular density on the graphene surface exerts a positive effect on molecular permeation. The permeance of H₂O molecules is the highest owing to their smallest diameter, while the permeance of CH₄ molecules is the lowest owing to their biggest diameter; in these cases, the molecular size is a dominating factor. For H₂S and CO₂ molecules, the diameters of H₂S molecules are larger than those of CO₂ molecules, but the interactions between H₂S molecules and graphene are stronger, resulting in a stronger permeation ability of H₂S molecules. Between CO₂ and N₂ molecules, CO₂ molecules show higher permeation rates owing to smaller diameters and stronger interactions with graphene. The graphene surface also shows nonuniform molecular density distribution owing to molecular permeation through graphene nanopores. Because of the doped nitrogen atoms, the CH₄ molecules prefer to permeate from the left and right sides of the graphene nanopore, while the other molecules prefer to permeate from the center of the nanopore owing to their small diameters. For the molecules that show stronger interactions with graphene, the molecular density on the graphene surface is higher; accordingly, the residence time on the graphene surface is longer and the experience time period during permeation is also longer. The mechanisms identified in this study can provide theoretical guidelines for the application of graphene-based membranes. In addition, the permeance of gas molecules in the graphene nanopore adopted in this study is on the order of 10^-3 mol·s^-1·m^-2·Pa^-1, and the selectivity of other molecules relative to CH₄ molecules is also high, showing that the membranes based on this type of nanopore can be employed in natural gas processing and other separation industries.

The construction of a photo-controllable artificial molecular machine capable of realizing the light-driven motion on a molecular scale and of performing a specific function is a fascinating topic in supramolecular chemistry. The bistable switchable molecule, azobenzene (AZO), has been introduced into the supramolecular architecture as a key building block, owing to its efficient and reversible trans (E)-cis (Z) photoisomerization. The binding strength of the dibenzo[24]crown-8 (DB24C8) host and dialkylammonium-based rod-like guest consisting of an AZO moiety and the Z$\to $E photoisomerization process in an interlocked host-guest complex have been investigated by the density functional theory (DFT) calculations and the reactive molecular dynamics (RMD) simulations by considering both torsion and inversion paths. The strong host-guest binding strength provides a necessary premise to stabilize the complex during the E-Z photoisomerization of the AZO unit, which is a terminal stopper to control the directional motion of the guest. A stronger binding strength for the Z isomer can be induced by the stronger hydrogen-bonding interaction. The steric effect is introduced into the Z isomer to force the ring slipping exclusively over the cyclopentyl terminal (pseudostopper). The host-guest complexation has a slight effect on the conformation of the AZO functional subunit for the two isomers. The faster Z$\to $E photoisomerization process within the picosecond timescale is kinetically more favored than the dethreading of the ring through the pseudostopper subunit of the rod. After isomerization, a structure relaxation is observed for the crown ether ring within 500 ps. The flexible backbone of the crown ether ring is helpful in realizing steady and stable host-guest recognition during photoisomerization. Moreover, the orthogonality of the site-specific binding interaction is revealed by the similar binding energies obtained at similar hydrogen bonding recognition sites for various interlocked host-guest supramolecular systems although the constituents of the guests are different from each other. The introduction of two stereoisomers of the AZO subunit has little influence on the other conformations of guest subunits. These results are useful for the rational design of more sophisticated stimuli-controlled artificial molecular machines.

Ethylene carbonate (EC) liquid and its vapor-liquid interface were investigated using a combination of molecular dynamics (MD) simulation and vibrational IR, Raman and sum frequency generation (SFG) spectroscopies. The MD simulation was performed with a flexible and polarizable model of the EC molecule newly developed for the computation of vibrational spectra. The internal vibration of the model was described on the basis of the harmonic couplings of vibrational modes, including the anharmonicity and Fermi resonance coupling of C＝O stretching. The polarizable model was represented by the charge response kernel (CRK), which is based on ab initio molecular orbital calculations and can be readily applied to other systems. The flexible and polarizable model can also accurately reproduce the structural and thermodynamic properties of EC liquid. Meanwhile, a comprehensive set of vibrational spectra of EC liquid, including the IR and Raman spectra of the bulk liquid as well as the SFG spectra of the liquid interface, were experimentally measured and reported. The set of experimental vibrational spectra provided valuable information for validating the model, and the MD simulation using the model comprehensively elucidates the observed vibrational IR, Raman, and SFG spectra of EC liquid. Further MD analysis of the interface region revealed that EC molecules tend to orientate themselves with the C＝O bond parallel to the interface. The MD simulation explains the positive Im[$ \chi ^{(2)}$](ssp) band of the C＝O stretching region in the SFG spectrum in terms of the preferential orientation of EC molecules at the interface. This work also elucidates the distinct lineshapes of the C＝O stretching band in the IR, Raman, and SFG spectra. The lineshapes of the C＝O band are split by the Fermi resonance of the C＝O fundamental and the overtone of skeletal stretching. The Fermi resonance of C＝O stretching was fully analyzed using the empirical potential parameter shift analysis (EPSA) method. The apparently different lineshapes of the C＝O stretching band in the IR, Raman, and SFG spectra were attributed to the frequency shift of the C＝O fundamental in different solvation environments in the bulk liquid and at the interface. This work proposes a systematic procedure for investigating the interface structure and SFG spectra, including general modeling procedure based on ab initio calculations, validation of the model using available experimental data, and simultaneous analysis of molecular orientation and SFG spectra through MD trajectories. The proposed procedure provides microscopic information on the EC interface in this study, and can be further applied to investigate other interface systems, such as liquid-liquid and solid-liquid interfaces.

Surfactant molecules, when dispersed in solution, have been shown to spontaneously form aggregates. Our previous studies on molecular dynamics (MD) calculations have shown that ionic sodium dodecyl sulfate molecules quickly aggregated even when the aggregation number is small. The aggregation rate, however, decreased for larger aggregation numbers. In addition, studies have shown that micelle formation was not completed even after a 100 ns-long MD run (Chem. Phys. Lett. 2016, 646, 36). Herein, we analyze the free energy change of micelle formation based on chemical species model combined with molecular dynamics calculations. First, the free energy landscape of the aggregation, ΔG_i+j^†, where two aggregates with sizes i and j associate to form the (i + j)-mer, was investigated using the free energy of micelle formation of the i-mer, G_i^†, which was obtained through MD calculations. The calculated ΔG_i+j^† was negative for all the aggregations where the sum of DS ions in the two aggregates was 60 or less. From the viewpoint of chemical equilibrium, aggregation to the stable micelle is desired. Further, the free energy profile along possible aggregation pathways was investigated, starting from small aggregates and ending with the complete thermodynamically stable micelles in solution. The free energy profiles, G(l, k), of the aggregates at l-th aggregation path and k-th state were evaluated by the formation free energy $\sum\limits_i {{n_i}\left( {l, k} \right)G_i^\dagger } $ and the free energy of mixing $\sum\limits_i {{n_i}(l, k){k_B}Tln({n_i}(l, k)/n(l, k))} $, where n_i(l, k) is the number of i-mer in the system at the l-th aggregation path and k-th state, with $n\left( {l, k} \right) = \sum\limits_i {{n_i}\left( {l, k} \right)} $. All the aggregation pathways were obtained from the initial state of 12 pentamers to the stable micelle with i = 60. All the calculated G(l, k) values monotonically decreased with increasing k. This indicates that there are no free energy barriers along the pathways. Hence, the slowdown is not due to the thermodynamic stability of the aggregates, but rather the kinetics that inhibit the association of the fragments. The time required for a collision between aggregates, one of the kinetic factors, was evaluated using the fast passage time, t_FPT. The calculated t_FPT was about 20 ns for the aggregates with N = 31. Therefore, if aggregation is a diffusion-controlled process, it should be completed within the 100 ns-simulation. However, aggregation does not occur due to the free energy barrier between the aggregates, that is, the repulsive force acting on them. This may be caused by electrostatic repulsions produced by the overlap of the electric double layers, which are formed by the negative charge of the hydrophilic groups and counter sodium ions on the surface of the aggregates.

A complex reaction, such as combustion, polymerization, and zeolite synthesis, involves a large number of elementary reactions and chemical species. Given a set of elementary reactions, the apparent reaction rates, population of chemical species, and energy distribution as functions of time can be derived using deterministic or stochastic kinetic models. However, for many complex reactions, the corresponding elementary reactions are unknown. Molecular dynamics (MD) simulation, which is based on forces calculated by using either quantum mechanical methods or pre-parameterized reactive force fields, offers a possibility to probe the reaction mechanism from the first principles. Unfortunately, most reactions take place on timescales far above that of molecular simulation, which is considered to be a well-known rare event problem. The molecules may undergo numerous collisions and follow many pathways to find a favorable route to react. Often, the simulation trajectory can be trapped in a local minimum separated from others by high free-energy barriers; thus, crossing these barriers requires prohibitively long simulation times. Due to this timescale limitation, simulations are often conducted on very small systems or at unrealistically high temperatures, which might hinder their validity. In order to model complex reactions under conditions comparable with those of the experiments, enhanced sampling techniques are required. The replica exchange molecular dynamics (REMD) is one of the most popular enhance sampling techniques. By running multiple replicas of a simulation system using one or several controlling variables and exchanging the replicas according to the Metropolis acceptance rule, the phase space can be explored more efficiently. However, most published work on the REMD method focuses on the conformational changes of biological molecules or simple reactions that can be described by a reaction coordinate. The optimized parameters of such simulations may not be suitable for simulations of complex reactions, in which the energy changes are much more dramatic than those associated with conformational changes and the hundreds elementary reactions through numerous pathways are unknown prior to the simulations. Therefore, it is necessary to investigate how to use the REMD method efficiently for the simulation of complex reactions. In this work, we examined the REMD method using temperature (T-REMD) and Hamiltonian (H-REMD) as the controlling variable respectively. In order to quantitatively validate the simulation results against direct simulations and analytic solutions, we performed the study based on a simple replacement reaction (A + BC = AB + C) with variable energy barrier heights and reaction energies described using the ReaxFF functional forms. The aim was to optimize the simulation parameters including number, sequence, and swap frequency of the replicas. The T-REMD method was found to be efficient for modeling exothermic reactions of modest reaction energy (< 3 kcal∙mol^-1) or activation energy (ca. < 20 kcal∙mol^-1). The efficiency was severely impaired for reactions with high activation and reaction energies. The analysis of the simulation trajectory revealed that the problem was intrinsic and could not be solved by adjusting the simulation parameters since the phase space sampled using T-REMD was localized in the region favored by high (artificial for speed-up) temperatures, which is different from the region favored by low (experimental) temperatures. This issue was aggravated in the case of endothermic reactions. On the other hand, the H-REMD run on a series of potential surfaces having different activation energies was demonstrated to be remarkably robust. Since the energy barrier only reduces the reaction rates, while the phase space controlled by the reaction energy differences remains unchanged at a fixed temperature, excellent results were obtained with fewer replicas by using H-REMD. It is evident that H-REMD is a more suitable method for the simulation of complex reactions.

The development of efficient catalysts for the hydrogenation of CO₂ to formic acid (FA) or formate has attracted significant interest as it can address the increasingly severe energy crisis and environmental problems. One of the most efficient methods to transform CO₂ to FA is catalytic homogeneous hydrogenation using noble metal catalysts based on Ir, Ru, and Rh. In our previous work, we demonstrated that the activity of CO₂ hydrogenation via direct addition of hydride to CO₂ on Ir(Ⅲ) and Ru(Ⅱ) complexes was determined by the nature of the metal-hydride bond. These complexes could react with the highly stable CO₂ molecule because they contain the same distinct metal-hydride bond formed from the mixing of the sd² hybrid orbital of metal with the 1s orbital of H, and evidently, this property can be influenced by the trans ligand. Since boryl ligands exhibit a strong trans influence, we proposed that introducing such ligands may enhance the activity of the Ru―H bond by weakening it as a result of the trans influence. In this work, we designed two potential catalysts, namely, Ru-PNP-HBcat and Ru-PNP-HBpin, which were based on the Ru(PNP)(CO)H₂ (PNP = 2, 6-bis(dialkylphosphinomethyl)pyridine) complex, and computationally investigated their reactivity toward CO₂ hydrogenation. Bcat and Bpin (cat = catecholate, pin = pinacolate) are among the most popular boryl ligands in transition metal boryl complexes and have been widely applied in catalytic reactions. Our optimization results revealed that the complexes modified by boryl ligands possessed a longer Ru―H bond. Similarly, natural bond orbital (NBO) charge analysis indicated `s\�vl VE`s\�vl Vr of the hydride in Ru-PNP-HBcat and Ru-PNP-HBpin was higher as compared to that in Ru-PNP-H₂. NBO analysis of the nature of Ru―H bond indicated that these complexes also followed the law of the bonding of Ru―H bond proved in the previous works (Bull. Chem. Soc. Jpn. 2011, 84 (10), 1039; Bull. Chem. Soc. Jpn. 2016, 89 (8), 905), and the d orbital contribution of the Ru atom in Ru-PNP-HBcat and Ru-PNP-HBpin was smaller than that in Ru-PNP-H₂. Consequently, the Ru-PNP-HBcat and Ru-PNP-HBpin complexes were more active than Ru-PNP-H₂ for the direct hydride addition to CO₂ because of the lower activation energy barrier, i.e., from 29.3 kJ∙mol^-1 down to 24.7 and 23.4 kJ∙mol^-1, respectively. In order to further verify our proposed catalyst-design strategy for CO₂ hydrogenation, the free energy barriers of the complete pathway for the hydrogenation of CO₂ to formate catalyzed by complexes Ru-PNP-H₂, Ru-PNP-HBcat, and Ru-PNP-HBpin were calculated to be 76.2, 67.8, and 54.4 kJ∙mol^-1, respectively, indicating the highest activity of Ru-PNP-HBpin. Thus, the reactivity of Ru catalysts for CO₂ hydrogenation could be tailored by the strong trans influence of the boryl ligands and the nature of the Ru―H bond.

Dye-sensitized solar cells (DSSCs) are one of the most promising renewable energy technologies. Charge transfer and charge transport are pivotal processes in DSSCs, which govern solar energy capture and conversion. These processes can be probed using modern electronic structure methods. Because of the heterogeneity and complexity of the local environment of a chromophore in DSSCs (such as solvatochromism and chromophore aggregation), a part of the solvation environment should be treated explicitly during the calculation. However, because of the high computational cost and unfavorable scaling with the number of electrons of high-level quantum mechanical methods, approaches to explicitly treat the local environment need careful consideration. Two problems must be tackled to reduce computational cost. First, the number of configurations representing the solvent distribution should be limited as much as possible. Second, the size of the explicit region should be kept relatively small. The purpose of this study is to develop efficient computational approaches to select representative configurations and to limit the explicit solvent region to reduce the computational cost for later (higher-level) quantum mechanical calculations. For this purpose, an ensemble of solvent configurations around a 1-methyl-8-oxyquinolinium betaine (QB) dye molecule was generated using Monte Carlo simulations and molecular mechanics force fields. Then, a fitness function was developed using data from inexpensive electronic structure calculations to reduce the number of configurations. Specific solvent molecules were also selected for explicit treatment based on a distance criterion, and those not selected were treated as background charges. The configurations and solvent molecules selected proved to be good representatives of the entire ensemble; thus, expensive electronic structure calculations need to be performed only on this subset of the system, which significantly reduces the computational cost.