Login In
2018 Volume 34 Issue 10
2018, 34(10): 1095-1096
doi: 10.3866/PKU.WHXB201803291
Abstract:
2018, 34(10): 1097-1105
doi: 10.3866/PKU.WHXB201712131
Abstract:
The development of efficient catalysts for the hydrogenation of CO2 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 CO2 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 CO2 hydrogenation via direct addition of hydride to CO2 on Ir(Ⅲ) and Ru(Ⅱ) complexes was determined by the nature of the metal-hydride bond. These complexes could react with the highly stable CO2 molecule because they contain the same distinct metal-hydride bond formed from the mixing of the sd2 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)H2 (PNP = 2, 6-bis(dialkylphosphinomethyl)pyridine) complex, and computationally investigated their reactivity toward CO2 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 that the nucleophilic character of the hydride in Ru-PNP-HBcat and Ru-PNP-HBpin was higher as compared to that in Ru-PNP-H2. 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-H2. Consequently, the Ru-PNP-HBcat and Ru-PNP-HBpin complexes were more active than Ru-PNP-H2 for the direct hydride addition to CO2 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 CO2 hydrogenation, the free energy barriers of the complete pathway for the hydrogenation of CO2 to formate catalyzed by complexes Ru-PNP-H2, 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 CO2 hydrogenation could be tailored by the strong trans influence of the boryl ligands and the nature of the Ru―H bond.
The development of efficient catalysts for the hydrogenation of CO2 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 CO2 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 CO2 hydrogenation via direct addition of hydride to CO2 on Ir(Ⅲ) and Ru(Ⅱ) complexes was determined by the nature of the metal-hydride bond. These complexes could react with the highly stable CO2 molecule because they contain the same distinct metal-hydride bond formed from the mixing of the sd2 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)H2 (PNP = 2, 6-bis(dialkylphosphinomethyl)pyridine) complex, and computationally investigated their reactivity toward CO2 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 that the nucleophilic character of the hydride in Ru-PNP-HBcat and Ru-PNP-HBpin was higher as compared to that in Ru-PNP-H2. 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-H2. Consequently, the Ru-PNP-HBcat and Ru-PNP-HBpin complexes were more active than Ru-PNP-H2 for the direct hydride addition to CO2 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 CO2 hydrogenation, the free energy barriers of the complete pathway for the hydrogenation of CO2 to formate catalyzed by complexes Ru-PNP-H2, 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 CO2 hydrogenation could be tailored by the strong trans influence of the boryl ligands and the nature of the Ru―H bond.
2018, 34(10): 1106-1115
doi: 10.3866/PKU.WHXB201701083
Abstract:
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.
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.
2018, 34(10): 1116-1123
doi: 10.3866/PKU.WHXB201801151
Abstract:
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.
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.
2018, 34(10): 1124-1135
doi: 10.3866/PKU.WHXB201801291
Abstract:
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[\begin{document}$ \chi ^{(2)}$\end{document}
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[
2018, 34(10): 1136-1143
doi: 10.3866/PKU.WHXB201801301
Abstract:
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: H2O > H2S > CO2 > N2 > CH4. 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 H2O molecules is the highest owing to their smallest diameter, while the permeance of CH4 molecules is the lowest owing to their biggest diameter; in these cases, the molecular size is a dominating factor. For H2S and CO2 molecules, the diameters of H2S molecules are larger than those of CO2 molecules, but the interactions between H2S molecules and graphene are stronger, resulting in a stronger permeation ability of H2S molecules. Between CO2 and N2 molecules, CO2 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 CH4 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 CH4 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.
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: H2O > H2S > CO2 > N2 > CH4. 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 H2O molecules is the highest owing to their smallest diameter, while the permeance of CH4 molecules is the lowest owing to their biggest diameter; in these cases, the molecular size is a dominating factor. For H2S and CO2 molecules, the diameters of H2S molecules are larger than those of CO2 molecules, but the interactions between H2S molecules and graphene are stronger, resulting in a stronger permeation ability of H2S molecules. Between CO2 and N2 molecules, CO2 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 CH4 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 CH4 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.
2018, 34(10): 1144-1150
doi: 10.3866/PKU.WHXB201802122
Abstract:
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.
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.
2018, 34(10): 1151-1162
doi: 10.3866/PKU.WHXB201802261
Abstract:
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 NO2, NO3, NO, and N2O. NO2 is the most abundant primary product and is gradually consumed in subsequent secondary reactions to form other nitrogen oxides. NO3 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―NO2 bond cleavage and ring-opening C―N bond scission dominate the initial decomposition of condensed CL-20. The cleavage of the N―NO2 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 NO2 to generate NO3 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 C4 and C2 intermediates are observed after the initial fragmentation, but eventually convert into stable products. The distribution of the final products (N2, H2O, CO2, and H2) 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.
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 NO2, NO3, NO, and N2O. NO2 is the most abundant primary product and is gradually consumed in subsequent secondary reactions to form other nitrogen oxides. NO3 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―NO2 bond cleavage and ring-opening C―N bond scission dominate the initial decomposition of condensed CL-20. The cleavage of the N―NO2 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 NO2 to generate NO3 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 C4 and C2 intermediates are observed after the initial fragmentation, but eventually convert into stable products. The distribution of the final products (N2, H2O, CO2, and H2) 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.
2018, 34(10): 1163-1170
doi: 10.3866/PKU.WHXB201802271
Abstract:
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, ΔGi+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, Gi†, which was obtained through MD calculations. The calculated ΔGi+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\begin{document}$\sum\limits_i {{n_i}\left( {l, k} \right)G_i^\dagger } $\end{document} \begin{document}$\sum\limits_i {{n_i}(l, k){k_B}Tln({n_i}(l, k)/n(l, k))} $\end{document} \begin{document}$n\left( {l, k} \right) = \sum\limits_i {{n_i}\left( {l, k} \right)} $\end{document}
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, ΔGi+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, Gi†, which was obtained through MD calculations. The calculated ΔGi+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
2018, 34(10): 1171-1178
doi: 10.3866/PKU.WHXB201803024
Abstract:
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\begin{document}$\to $\end{document} \begin{document}$\to $\end{document}
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
2018, 34(10): 1179-1188
doi: 10.3866/PKU.WHXB201803161
Abstract:
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.
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.
