Unraveling the Key Role of the Benzyl Group in the Synthesis of CL-20 Precursor HBIW.

As one of the most important energetic material molecules, hexanitrohexaazaisowurtzitane (CL-20) can only be synthesized using an amine with a benzyl group. Moreover, the reaction mechanism remains unexplored and the special role of the benzyl group has not been revealed. To address these issues, we perform an extensive theoretical study to investigate the synthesis mechanism of CL-20 precursor HBIW by employing density functional theory. Our calculated results demonstrate that the benzyl group can reduce the energy of the intermediate and the energy barrier of the rate-determining step due to the π-π stack interaction between two benzene rings of the benzyl group. For the first time, we revealed that the reactions can produce 16 intermediates with different chiralities during the formation of the first two side five-membered rings and only two of which can further form the bottom six-membered ring and finally obtain the product HBIW. The steric hindrance effect of the benzyl group leads to the formation of a higher-energy intermediate first, thereby providing an opportunity to correct the wrong chirality. All of these factors make the diimine with the benzyl group the most suitable reactant for the synthesis of CL-20.

Introduction

So far, hexanitrohexaazaisowurtzitane (CL-20) is the best-performance energetic material that can be industrially produced. Although the acid catalysis conditions,[1,2] solvent conditions,[3] choice of raw materials,[4−7] nitration steps,[8−11] and production processes[12,13] have been improved recently, the main synthetic strategy of CL-20 remains unchanged in the past 30 years.

In 1990, Nielsen et al. reported the synthetic strategy of CL-20.[14] As shown in Figure , the synthetic strategy of CL-20 includes three steps: a cage structure HBIW is formed by benzylamine with glyoxal, then HBIW is debenzylated, and finally the product CL-20 is formed by nitration. The nitration reaction is almost a necessary reaction step for all energetic materials. However, the benzene ring of HBIW is directly attacked by the nitration reagent, and the nitration reaction site is affected by the presence of the benzyl group.[15,16] Therefore, to not affect the subsequent nitration reaction, CL-20 needs to be debenzylated, but the yield is low in this step, and the expensive heavy metal Pd is used, which greatly affects the synthetic cost of CL-20.

Figure 1

Synthetic strategy of CL-20.

To solve the problem of debenzylation, it can be carried out from two aspects. On the one hand, the debenzylation process can be optimized and more suitable catalysts can be found. The current debenzylation reactions mainly include hydrogenolysis debenzylation reaction (Pearlman’s catalyst),[17−19] oxidative debenzylation reaction,[20−22] nitrosolysis debenzylation reaction,[23−25] chloroformate debenzylation reaction,[26] and other debenzylation methods[27,28] among which the Pearlman’s catalyst is still considered the most commonly used method for debenzylation. Therefore, it remains a great challenge to find a suitable method to optimize the debenzylation process.

On the other hand, reactants other than benzylamine analogues can be used to directly synthesize CL-20,[29] so that the debenzylation step can be skipped, which is beneficial to improve the reaction efficiency and reduce the reaction cost. After many studies, in addition to the use of benzylamine analogues, allylamine and propargylamine can also be used as reactants to form the cage structure,[30,31] but the yield is not satisfactory. If the benzyl group is not involved during the reaction, the second step of debenzylation is not required. Therefore, it is necessary to theoretically explore the particularity of the benzyl group during the synthesis process and compare it with other groups to find a more suitable reaction method or reactive group.

As shown in Figure , the HBIW formation reaction by benzylamine with glyoxal includes two steps: diimine is formed by aldehyde amine condensation, and three diimine monomers are polymerized into a cage structure. Aldehyde amine condensation can also occurs when amines other than benzylamine are used, and the corresponding diimine substances can be extracted and distinguished by X-ray diffraction (XRD).[14] The cage structure can be successfully synthesized by N,N′-dibenzyl-1,2-ethylenediimine, but most of the diimine separated from other amines and glyoxal are stable substances.[14] Therefore, the formation of the cage structure is affected by the polymerization process of diimines. Pang et al. have calculated the reaction path of HBIW.[32] However, due to the large number of atoms in the molecule and the great flexibility of the benzyl group, the ONIOM (M062X/DGDZVP: PM6) method is used to simplify the benzyl group, and it is difficult to study the particularity of the benzyl group. In this study, the density functional theory (DFT) method is used to study the cage formation path of N,N′-dimethyl-1,2-ethylenediimine, N,N′-diallyl-1,2-ethylenediimine, and N,N′-dibenzyl-1,2-ethylenediimide; then, the energy barriers of the rate-determining step were compared, and the effect of chirality on the synthesis of HBIW is discovered for the first time.

Figure 2

Proposed synthesis steps of HBIW.

Computational Details

All geometry optimizations were performed at the density functional theory (DFT) level using the hybrid B3LYP functional[33−35] with the basis set 6-31G(d).[36] Grimmie’s D3 dispersion correction[37] was employed in the calculations to appropriately evaluate the intermolecular interactions. The nature of the optimized stationary points was characterized by frequency calculations at the same level of theory, with all reactant, intermediate, and product configurations having no imaginary frequency, and all transition states (TSs) having one imaginary mode. Thermal corrections to the electronic energies of the optimized geometries were estimated at an experimental temperature of 298 K and 1 atm from the frequency calculations. The intrinsic reaction coordinate (IRC)[38] calculations from first-order saddle points were performed to locate the local minima for the reaction pathways, which are denoted by solid lines in the energy profile. Solvent effects were implicitly taken into account by means of the solvation model based on the density (SMD) method.[39] To obtain more accurate electronic energies as well as reaction barrier heights, we further carried out single-point calculations based on the optimized geometries using the B3LYP functional with larger basis sets 6-311+g(2d,p).[40] All DFT calculations were performed using Gaussian 16 program.[41] The molclus program[42] combined with the extended tight-binding (xTB) method[43,44] is used in the conformation search of the intermediate structure of the corresponding N,N′-dibenzyl-1,2-ethanediimine in the subsequent reaction. First, 800 initial structures were generated by turning the angle of the benzyl group based on the gentor module in the molclus program.[42] Then, the single-point energies were calculated and sorted by the low-precision GFN-xTB method, and the 10 structures with the lowest energy were obtained. Finally, these structures are calculated using high-precision quantum chemical methods to obtain the lowest-energy conformation. To study weak interactions, reduced density gradient (RDG) analysis is conducted based on Multiwfn software.[45]

Results and Discussion

Proposed Synthesis Mechanism of HBIW

All of the reactions are carried out under acidic conditions. Our calculation results reveal that the terminal N atom of diimine is easily combined with H+. The first reaction is the polymerization between diimine and protonated diimine. By scanning the energy profiles of this polymerization reaction at different binding sites, we found that the N of diimine is easily combined with the C of protonated diimine, and the above reaction site remains unchanged when the functional group of diimine is changed (more calculation details are shown in part 1 of the Supporting Information).

The HBIW cage structure is composed of two side five-membered rings and one bottom six-membered ring. The construction of the cage is affected by the ring formation pathway. As shown in path A in Figure , a cage formation method is that two side five-membered rings are formed first, then a bottom six-membered ring is formed. The first side five-membered ring is formed by the diimine with the first protonated diimine, and the second side five-membered ring is formed by the dimer with the second protonated diimine. Subsequently, the bottom six-membered ring is formed by first rotating the C–C bridge bond and then closing the ring. As shown in path B in Figure , another cage formation method is that a bottom six-membered ring is formed first, then two side five-membered rings are formed. The first bottom six-membered ring is formed by two protonated diimines, then the dimer is attacked from above by the diimine, forming two side five-membered rings. The cage formation in path B in Figure has been verified to be unable to react due to the high reaction energy barrier (more calculation details are shown in part 2 of the Supporting Information).

Figure 3

Proposed synthesis path A of HBIW.

Figure 4

Proposed synthesis path B of HBIW.

Reaction Pathway of the Cage Structure

Owing to a large number of atoms in the system and the great flexibility in the benzyl group, directly investigating the cage formation of N,N′-dibenzyl-,1,2-ethylenediimide remains a great challenge. Therefore, the cage formation reaction of N,N′-dimethyl-1,2-ethylenediimine is first studied to determine the rate-determining step, and then the special effect of the benzyl group is revealed through comparisons.

In Figure A, we show the energy profiles of the formation of the first side five-membered ring by N,N′-dimethyl-1,2-ethylenediimide with protonated N,N′-dimethyl-1,2-ethylenediimine. First, N1 of N,N′-dimethyl-1,2-ethylenediimide attacks C7 of protonated N,N′-dimethyl-1,2-ethylenediimine via the transition state TS1 with an energy barrier of 10.1 kcal/mol and leads to the formation of the intermediate IM1. To connect N5 to C2, IM1 first crosses a transition state TS2 with the energy barrier of 1.5 kcal/mol by rotating the dihedral angle N5–C6–C7–N1 to generate intermediate IM2. N5 and C2 are then combined via a transition state TS3 with the energy barrier of 7.0 kcal/mol and then leads to the formation of the intermediate IM3. IM3 is a relatively stable intermediate, and the first five-membered ring N1C2N5C6C7 in the cage structure is formed.

Figure 5

(A) Gibbs free energy profile for the first side five-membered ring formation reaction of N,N′-dimethyl-1,2-ethanediimine with protonated N,N′-dimethyl-1,2-ethylenediimine. (B) Gibbs free energy profile for the second side five-membered ring formation reaction of IM3 with protonated N,N′-dimethyl-1,2-ethylenediimine. (C) Molecular structures in the reaction pathways.

In Figure B, we show the energy profiles of the formation of the second side five-membered ring by IM3 with protonated N,N′-dimethyl-1,2-ethylenediimine in which the energy of IM3 was set to be zero. First, N4 of IM3 attacks C11 of protonated N,N′-dimethyl-1,2-ethylenediimine via the transition state TS4 with the energy barrier of 15.0 kcal/mol and leads to the formation of the intermediate IM4. Next, N9 and C3 of IM4 are bonded to close the ring via a transition state TS5 with the energy barrier of 7.8 kcal/mol, leading to the formation of the intermediate IM5. IM5 is a relatively stable intermediate in which the second five-membered ring N4C3N9C10C11 in the cage structure is formed. Since there are a large number of conformational transition reactions in the subsequent reaction steps after IM5 and the energy surface is relatively smooth, we thereby show the computed energy profiles of the subsequent reactions in Figures S3–S5, as well as a detailed discussion on these reactions. It can be seen that the formation process of the cage structure from three diimine monomers mainly includes three steps: formation of two side five-membered rings by polymerization of three diimine monomers, formation of the bottom six-membered ring by conformational transformation, and the proton extraction process.

The rate-determining step of the overall reaction is the formation of the second side five-membered ring, which is the reaction path from IM4 to IM5. Due to the reaction of free-motion monomers into multimers being an obvious entropy reduction process, it is also predictable that the transition state TS5 possesses the highest energy barrier in the overall reaction process. Therefore, it is necessary to compare the energy of TS5 of different diimines to reveal the substitute effects. In the following, the reaction pathways of N,N′-diallyl-1,2-ethylenediimine and N,N′-dibenzyl-1,2-ethylenediimide are described.

In Figures and 7, we show the energy profiles of the formations of the first and second side five-membered rings by N,N-diallyl-1,2-ethylenediimine and N,N′-dibenzyl-1,2-ethylenediimine, respectively. Combining Figures –7, we can see that the substituents do not affect the reaction paths of the formation of the first five-membered ring. However, both the allyl and benzyl groups can reduce the energy barrier of the overall reaction of this process. For the formation of the second five-membered ring, IM3A and IM3B need to rotate the ally and benzyl groups to facilitate the subsequent reactions, respectively. The energy barriers of these rotation processes are less than 2 kcal/mol (see Figures B and 7B). Importantly, from the methyl to allyl to benzyl groups, the energy barriers of the formation of the second five-membered ring decrease from 22.2 to 21.2 to 19.6 kcal/mol, respectively. Experimentally, one could not obtain the HBIW using N,N′-dimethyl-1,2-ethylenediimine as a reactant. When N,N′-diallyl-1,2-ethanediimine was used as a reactant, the yield was too low to be used in practice.[46] The yield can reach ∼20% when N,N′-dibenzyl-1,2-ethanediimine was used as a reactant, and this is currently the only method that can be used to synthesize CL-20. Our calculated results demonstrated that the yield increased gradually with the decrease of the rate-determining step energy barriers, which is consistent with the experimental expectation. Therefore, the lower rate-determining step energy barrier of N,N′-dibenzyl-1,2-ethanediimine is one of the key factors for the synthesis of the cage structure HBIW. However, the difference between the computed energy barriers is only ∼1 kcal/mol, which is not large enough to explain the significant difference in the experimentally observed yield. There should be other factors that still affect the reaction mechanism.

Figure 6

(A) Gibbs free energy profile for the first side five-membered ring formation reaction of N,N′-diallyl-1,2-ethanediimine with protonated N,N′-diallyl-1,2-ethanediimine. (B) Gibbs free energy profile for the second side five-membered ring formation reaction of IM3A with protonated N,N′-diallyl-1,2-ethanediimine. (C) Molecular structures in the reaction pathways.

Figure 7

(A) Gibbs free energy profile for the first side five-membered ring formation reaction of N,N′-dibenzyl-1,2-ethanediimine with protonated N,N′-dibenzyl-1,2-ethanediimine. (B) Gibbs free energy profile for the second side five-membered ring formation reaction of IM3B with protonated N,N′-dibenzyl-1,2-ethanediimine. (C) Molecular structures in the reaction pathways.

We noticed that for N,N′-dimethyl-1,2-ethylenediimine and N,N′-diallyl-1,2-ethanediimine, the energies of IM4 and IM4A are only lower than those of TS4 and TS4A, approximately 0.6 and 1.4 kcal/mol, respectively. These results demonstrated that IM4 and IM4A are unstable and can be easily returned to IM3 and IM3A, respectively, thereby impeding the subsequent reactions. For N,N′-dibenzyl-1,2-ethylenediimine, the energy of IM4B is lower than that of TS4B 2.1 kcal/mol, indicating that IM4B is more stable than IM4 and IM4A and can facilitate the subsequent reactions. The RDG analysis reveals that the stabilization of IM4B results from strong π–π stacking interaction between two benzene rings of benzyl groups (see Figure S6 and related discussion in the Supporting Information). Therefore, the other major role of benzyl is to lower the energy of IM4B, obtaining more stable reaction intermediates.

Reaction Pathway of the Cage Structure via the Cis-Reactant

Although trans-ethylenediamine is the most stable monomer, the energy of cis monomer is only about 1 kcal/mol higher than that of the trans monomer. Therefore, a small fraction of the cis-ethylenediamine monomer may also exist in the solution and participate in the reaction. Since the rate-determining step is the reaction to generate the second five-membered ring intermediate IM5, the energy barrier of this step has also been computed. As shown in Figure , the reaction of the cis-reactant with intermediate IM3 is a synergistic step, and the energy barriers (TS5′, TS5A′, TS5B′) are higher than those of the reactions of the trans-reactant with intermediate IM3. The reactions involving such cis-monomers may also contribute to the overall reactions but only in a limited way.

Figure 8

(A) From left to right, Gibbs free energy profile for the second side five-membered ring formation reaction of IM3 with protonated cis-N,N′-dimethyl-1,2-ethanediimine, cis-N,N′-diallyl-1,2-ethanediimine, and cis-N,N′-dibenzyl-1,2-ethanediimine. (B) Molecular structures in the reaction pathways.

Failure in Cage Formation Due to the Chiral Issue of Intermediate

Although the benzyl group can reduce the energy barrier of the rate-determining step, the reduced energy barrier is small and ∼2.6 kcal/mol when compared with that of the methyl group, which could not explain why no products can be obtained using N,N′-dimethyl-1,2-ethylenediimide. Moreover, it remains unexplored why the yield cannot be increased by increasing the temperature.[12] Therefore, there should be other factors affecting the reaction yield that have not been revealed.

In Figure A, during the formation of the first side five-membered ring, since the protonated diimine has a planar structure, the protonated diimine can be flipped 180°, and C7 attacks the N on the other diimine from different directions, resulting in different chiralities of C7. Similarly, in Figure B, the involvement of the third diimine during the formation of the second side five-membered ring can also make C11 possess different chiralities. Finally, when the second side five-membered ring was formed, four atoms (C2, C3, C7, and C11) of the formed intermediate may possess different chiralities. C2 and C3 atoms are the bridge C atoms connecting the two side five-membered rings, whereas C7 and C11 atoms are the ring C atoms on the side five-membered ring. Whether the amino group −NR is on the same side or on the different side is mainly affected by the chiralities of C2 and C3. Whether the amino group −NR is inside or outside the two five-membered rings is mainly affected by the chirality of C7 and C11. Depending on the chirality of these four atoms, 16 possible intermediates can be obtained, as shown in Figure . If we do not consider the differences between the two five-membered rings, there are only eight different structures, i.e., SSSS and RRSS are the same structure. Importantly, the formation of the bottom six-membered ring is affected by the chirality of these four atoms. Only when the chirality of C2, C3, C7, and C11 conforms to SSSS or RRSS type, the bottom six-membered ring can be formed. Since the direction of attack is nonselective, it is impossible to avoid the generation of the wrong chiral structure. Figure reveals that, for the intermediate with SSRR or RRRR chirality, the hydrogens are located between the two five-membered rings, and the amine groups −NR are located above and below the five-membered ring. Therefore, the process of generating these two intermediates is likely to encounter minimal steric hindrance effects and possesses the lowest-energy barrier, which can be used as a typical representative for the wrong reaction path.

Figure 9

(A) Comparison of the attack direction of diimine during the formation of the first five-membered ring. (B) Comparison of the attack direction of diimine during the formation of the second five-membered ring.

Figure 10

Sixteen chiral structures during the cage formation process.

To unveil the effect of the formation of wrong chiral intermediates, the reaction path of the formation of the intermediate with SSRR chirality was taken as an example. In Figure A, we show the energy profile of the formation of the first side five-membered ring by N,N′-dimethyl-1,2-ethylenediimide with protonated N,N′-dimethyl-1,2-ethylenediimine. Comparing Figure A with Figure A, we can find that the formation of the first side five-membered ring with right chirality needs three elementary reaction steps, whereas the formation process with wrong chirality only needs two elementary reaction steps and possesses a lower energy barrier. The overall energy barrier of the formation of the second side five-membered ring with wrong SSRR chirality is 17.1 kcal/mol and lower than that with right chirality ∼5.1 kcal/mol, as shown in Figures B and 11B. These results demonstrate that the reactions prefer to form intermediates with wrong chirality and prevent further formation of the bottom six-membered ring when N,N′-dimethyl-1,2-ethylenediimide was used as the reactant. This should also be the major reason why no product can be obtained using N,N′-dimethyl-1,2-ethylenediimide as the reactant.

Figure 11

(A) Gibbs free energy profile for the first side five-membered ring formation reaction of N,N′-dimethyl-1,2-ethanediimine with protonated N,N′-dimethyl-1,2-ethanediimine (the chiral type is SSRR). (B) Gibbs free energy profile for the second side five-membered ring formation reaction of IM14 with protonated N,N′-dimethyl-1,2-ethanediimine (the chiral type is SSRR). (C) Molecular structures in the reaction pathways.

Since the reaction paths of different diimines are similar, only the reaction paths near the rate-determining step are described in the following. In Figure , we show the energy profiles of the formation of the second side five-membered ring with the chiral type of SSRR by IM14A with protonated N,N′-diallyl-1,2-ethylenediimine and IM14B with protonated N,N′-dibenzyl-1,2-ethylenediimine, respectively. Comparing with the formation processes with right chirality, we found that the overall energy barriers of the formation of the second side five-membered ring with wrong chirality reduced by ∼1.1 kcal/mol for both N,N′-diallyl-1,2-ethylenediimine and N,N′-dibenzyl-1,2-ethylenediimine. The small difference between the energy barriers of reactions with right and wrong chiralities indicates that both the intermediates with right and wrong chiralities can be formed simultaneously.

Figure 12

(A) Gibbs free energy profile for the second side five-membered ring formation reaction of IM14A with protonated N,N′-diallyl-1,2-ethanediimine (the chiral type is SSRR). (B) Gibbs free energy profile for the second side five-membered ring formation reaction of IM14B with protonated N,N′-dibenzyl-1,2-ethanediimine (the chiral type is SSRR). (C) Molecular structures in the reaction pathways.

Conformation Transformation of Intermediates

Comparing Figures –7 and Figures and 12, we can find that the reaction paths of different reactants are similar. However, there is an important difference in conformation of the formed intermediate with the second side five-membered ring. Intermediates IM5B and IM16B generated by N,N′-dibenzyl-1,2-ethanediimine possessed relatively high energy due to the steric hindrance caused by the large volume of the benzyl group, which will convert to other conformations with lower energy. By performing the conformation search of IM5B with molclus program[42] combined with the xtb method,[43] we found that the conformation IB5B′ (see Figure S7) possesses the lowest energy. The computed energy profiles along the constructed linearly interpolated internal coordinate (LIIC) pathways between IM5B and IM5B′ (see Figure S7) reveal that the energy barrier of conformation changes from IB5B to IM5B′ is not larger than 8.5 kcal/mol. IM5B′ is the required structure for the next step to form the bottom six-membered ring, whereas it is impossible to form the bottom six-membered ring directly from IM5B. For N,N′-dimethyl-1,2-ethanediimine and N,N′-diallyl-1,2-ethanediimine, the formed IM5, IM5A, IM16, and IM16A possess a similar conformation to that of IM5B′. IM5 and IM16 can be used to proceed directly to the subsequent reaction without conformational transformation. Therefore, owing to the steric hindrance of the benzyl group of N,N′-dibenzyl-1,2-ethanediimine, a high-energy intermediate conformation is always formed first in the reaction. This high-energy conformation can be gradually converted into low-energy conformations through the rotation of the benzyl group or side five-membered ring (more calculation details are shown in the Supporting Information part 5). Since the conformation transformation does not happen easily, high-energy conformation has a certain lifetime. In the path of wrong chirality of N,N′-dibenzyl-1,2-ethanediimine (see Figure B), only 9.5 kcal/mol is needed during the reverse reaction from IM16B to IM15B, which contributes to the decomposition of the wrong chiral structure. Compared with other diimines, the reverse reaction energy barrier is lower. The reactant produced by the decomposition of the wrong chiral structure can participate in the reaction again, possibly following the right reaction path and forming the intermediate with the right chirality. Compared with other groups, intermediates with benzyl groups are more capable of correcting wrong chiral structures.

Conclusions

In summary, the synthesis mechanism of CL-20 precursor HBIW using diimine with three different functional groups has been investigated with the DFT method, and the particularity of the benzyl group in the cage formation reaction has been revealed. First, owing to the π–π stacking interaction between two benzene rings of benzyl groups, intermediate IM4B produced by N,N′-dibenzyl-1,2-ethanediimine possesses lower energy and can facilitate the subsequent reaction to form the second side five-membered ring. Second, comparing the diimine with methyl and allyl groups, the diimine with benzyl group possesses the lowest-energy barrier in the rate-determining step. Therefore, the cage structure of HBIW is easier to be generated by N,N′-dibenzyl-1,2-ethanediimine at room temperature. Third, depending on the attack direction between two reactants, four C atoms of the formed intermediate with the second side five-membered ring possess different chiralities. Among the generated 16 intermediates, only the ones with SSSS and RRSS chirality can further form the bottom six-membered ring. For N,N′-dimethyl-1,2-ethanediimine, the energy barrier of the formation of intermediates with wrong chirality is lower than that with right chirality, ∼5.1 kcal/mol, thereby facilitating the formation of intermediates with wrong chirality and preventing the formation of the final cage structure. These results well explain the experimental observation that no cage structure can be obtained using N,N′-dimethyl-1,2-ethanediimine as the reactant. For N,N′-diallyl-1,2-ethylenediimine and N,N′-dibenzyl-1,2-ethylenediimine, the energy barrier of the formation of the intermediate with wrong chirality is smaller than that with right chirality, ∼1.1 kcal/mol. Therefore, both the intermediates with right and wrong chiralities can be obtained. Finally, owing to the steric hindrance effect of the benzyl group, some high-energy intermediates are first produced in the formation of the intermediate with the second side five-membered ring, and the transition to a low-energy conformation is slower. The reverse reaction of these high-energy structures is easier, and the wrong chiral structure may be decomposed into reactants and then enter the correct reaction path. When a complex spatial structure is formed, more attention needs to be paid to the influence of chirality issues. Therefore, our calculation results reveal the key roles of the benzyl group in the synthesis of CL-20 and can well explain the experimentally observed substitute effects. Finally, it should be noted that it is difficult to find all of the possible reaction paths only by static DFT calculations. Moreover, the solvent acetonitrile may also affect the reactions via noncovalent interactions, which may not be properly described by the polarizable continuum model (PCM) solvation model. To recover these two overlooked effects, further ab initio molecular dynamics and QM/MM calculations are required.