md / qm / csm method for screening representative conformations of cb[7] host-guest systems

Abstract

The application discloses an MD / QM / CSM method for screening representative conformations of CB[7] host-guest system. The method comprises the following steps: firstly, topological structures and RESP charges of host molecules and guest molecules are calculated respectively; secondly, initial conformations of host-guest complexes in four different orientations are constructed by using Gaussview and Autodock molecular docking, and the search conformation space is expanded; thirdly, MD simulation is carried out on each initial conformation; fourthly, the advantage conformation is screened out through cluster analysis; fifthly, the complex free energy of the advantage conformation is calculated by using the QM / CSM method; and finally, the calculation value of the most stable advantage conformation is compared with experimental data to determine the representative conformation. The application can accurately screen the representative conformation of each CB[7] complex system, and obtain the complex free energy of the host-guest system. The calculated complex free energy has good linear correlation with the experimental value, R 2 = 0.98, and the change trend is consistent with the experimental value.

Description

Technical Field

[0001] This invention belongs to the field of theoretical calculation of supramolecular host-guest systems, specifically an MD / QM / CSM method for screening representative conformations of CB[7] host-guest systems. Background Technology

[0002] Neurotoxic drugs possess strong hallucinogenic and addictive properties, and their abuse seriously endangers individual and public health. Research has found that when some drugs are encapsulated by a host molecule, their side effects and toxicity in the body are significantly reduced (Organic & Biomolecular Chemistry, 2016, 14:7563-7569). Based on this finding, the application of host macromolecules as antidotes has been extensively studied. Currently, supramolecular compounds such as cyclodextrins, cucurbiturils, acyclic cucurbiturils derivatives, calixarenes, and columnar aromatics have shown good efficacy in addressing the side effects of neurotoxic drugs, pesticide poisoning, and the toxicity of illicit drugs, and have attracted widespread attention from researchers in the field of detoxification.

[0003] Cucurbituril[7] (CB[7]) is the most attention-grabbing host molecule among fourth-generation macrocyclic molecules. It has a hydrophobic cavity similar to that of cyclodextrin, with seven negatively charged carboxyl groups at the cavity opening. Therefore, the ideal structure of the guest molecule that interacts with cucurbituril[7] is a positively charged amphiphilic guest molecule, in which the positively charged part interacts with the negatively charged carboxyl groups at the cavity opening of cucurbituril by charge-dipole interaction, the hydrophobic part is complexed in the hydrophobic cavity of cucurbituril, and the hydrophilic part is exposed to the solvent to produce a solvent effect. Some guest molecules also form intermolecular hydrogen bonds with cucurbituril. The synergy of these non-covalent interactions makes cucurbituril have a high binding affinity with the guest molecule. Analyzing the structure and recombination mechanism of the CB[7] host-guest complex is of great significance for a deep understanding of the molecular recognition phenomenon that occurs in this system. With the development of computers and people's continuous exploration of computational methods, theoretical calculation plays an increasingly important role, providing reliable theoretical support for experimental data and providing a new perspective for the design of new supramolecular materials.

[0004] Currently, there are many theoretical calculation methods for predicting host-guest recombination ability. In recent years, Gilson et al. organized the SAMPL (Statistical Assessment of the Modeling of Proteins and Ligands) series of studies, which comprehensively evaluated the commonly used theoretical methods for calculating host-guest recombination free energy. In the summary of SAMPL research results, they pointed out that most methods can calculate the general trend of host-guest recombination free energy changes, but they still have accuracy problems. The theoretical calculation values ​​either have low linear correlation with experimental measurements or have large root mean square errors. At the same time, the universality and portability of the calculation methods still have large problems (J Comput Aided Mol Des, 2017, 31(1):1-19). Grimme's research group compared and evaluated the correlation between the two different methods, QM / CSM and MD, in calculating the binding affinity of CB[7] and 22 kinds of alkane guest molecules and experimental results. In the QM / CSM calculation, the conformation is classified and the structure is optimized. The DFT-D3 method was used to calculate the gas binding free energy, and the cosmo solvation model was used to calculate the solvation energy. In the MD calculation method, the attach-pull-release (APR) method was used to calculate the difference in free energy between the bound and unbound states. The results showed that the APR method had better correlation. In addition, both methods have their own limitations. The QM / CSM method has a large computational cost and uses fewer conformations, resulting in a large error. Although the MD method incorporates a wide range of sampled conformations and uses an explicit solvent model, it can more accurately handle the interaction between solvent molecules and solute molecules. However, the potential function used by the MD method is relatively coarse, and the error is large in handling the host-guest interaction (Journal of Physical Chemistry B, 2017, 121(49): 11144-11162).

[0005] Chinese patent ZL2021103378002 discloses an improved MD / QM / CSM method for extracting representative conformations of the β-CD host-guest system. This method uses MD to simulate the recombination process of β-CD and guest molecules, uses cluster analysis to screen representative conformations, and finally compares the results calculated by the representative conformations with experimental values ​​to verify the rationality of the representative conformations. However, this method still has the following problems: (1) The initial conformation of the MD simulation can have a significant impact on the MD results. This method does not fully search and construct the initial conformation space of the MD simulation, which will affect the accuracy of the representative conformation screening; (2) This method first establishes the representative conformation at the MM level, and then uses QM / CSM calculation for verification. It does not fully search the representative conformation space, which will also affect the accuracy of the representative conformation screening; (3) When calculating the QM vacuum free energy, this method uses the PM3 method. The PM3 method has low calculation accuracy and will produce a large error when calculating the interaction between the host and guest. Summary of the Invention

[0006] The purpose of this invention is to provide an MD / QM / CSM method for screening representative conformations of CB[7] host-guest systems, and to synthesize supramolecular host-guest systems by computer simulation.

[0007] The technical solution to achieve the objective of this invention is as follows:

[0008] The MD / QM / CSM method for screening representative conformations of the CB[7] subject-object system includes the following steps:

[0009] Step 1: Apply the GAFF force field to the host molecule CB[7] to obtain the topological structure of CB[7] and calculate the RESP charge at the HF / 6-31G* level;

[0010] Step 2: Using a GAFF force field, the topological structure of the amino-containing guest molecule is obtained and the RESP charge is calculated at the HF / 6-31G* level.

[0011] Step 3: Using Gaussview, two initial conformations with different orientations are constructed for the host molecule after charge and topology calculations in Step 1 and the guest molecule after charge and topology calculations in Step 2. Two initial conformations with different orientations are also constructed using Autodock software, for a total of four initial conformations with different orientations.

[0012] Step 4: Perform MD simulations on the four initial conformations with different orientations to obtain four MD simulation trajectories. Then, adjust the threshold and perform cluster analysis and structural optimization on the complex conformations in each MD trajectory to screen out the dominant conformations of the complexes for each trajectory. At the same time, perform MD simulations on the host molecule in Step 1 and the guest molecule in Step 2 respectively, adjust the threshold, perform cluster analysis and structural optimization to screen out the dominant conformations of individual host molecules and individual guest molecules for subsequent QM / CSM calculations of individual host molecules and individual guest molecules.

[0013] Step 5: Perform QM / CSM calculations on the dominant conformation of the complex for each trajectory obtained in Step 4, and select the most stable dominant conformation of the complex.

[0014] Step 6: Compare the calculated results of the most stable dominant conformation of the complex obtained in Step 5 with the experimental values. Perform a linear correlation analysis between the calculated recombination free energy and the experimental values. If the relative stability of the host-guest system is consistent with the experimental results, and the calculated values ​​and experimental values ​​have a good linear correlation, then the most stable dominant conformation of the complex is the accurate representative conformation of the system. If the relative stability of the host-guest system is inconsistent with the experimental results, or the calculated values ​​and experimental values ​​do not have a good linear correlation, then return to Step 3.

[0015] In this invention, the structural formula of CB[7] is:

[0016]

[0017] In this invention, some of the guest molecules contain chiral structures. The guest molecules are selected from METH (containing chiral structures), FEN, PCP, COC, or KET (containing chiral structures), and the structural formulas of each guest molecule are as follows:

[0018]

[0019] Furthermore, in step 1, the GLYCAM04 force field and Amber99SB are used as force field parameters.

[0020] Furthermore, in step 3, the host molecule CB[7] and the guest molecule are in a 1:1 ratio to construct the initial conformation of the host-guest complex.

[0021] Furthermore, in step 4, the specific method for selecting the preferred conformation of each trajectory complex is as follows:

[0022] (1) The initial structure of the host-guest complex was dissolved in TIP3P water, and the buffer distance from the host-guest complex to the boundary of the periodic aqueous solution box was set to Long-range electrostatic interactions were handled using the PME method, with the cutoff distance for non-covalent interactions set as... MD simulations were performed under periodic boundary conditions, and the SHAKE algorithm was used to constrain the H-containing chemical bonds. First, the host molecule was constrained to perform the first step of structural optimization to eliminate any abnormal interactions that might exist between the solvent water molecules and the host molecule and its complex. Then, a second structural optimization was performed under unconstrained conditions. After optimization, the system was heated. The Langevin kinetics were used to heat the system from 400 ps to 300 K. After heating, the NPT ensemble was used to perform MD simulations at 1 atm and 300 K. The integration step size was set to 2 fs, and the energy and structure were saved every 2 ps. The total MD time was 4 ns, and a total of 2000 conformations were obtained from the MD trajectory.

[0023] (2) Cluster analysis was performed on the conformations generated by MD simulation. By adjusting the RMSD threshold, the conformations were divided into 6 to 8 categories. The structures of the top two largest clusters were selected and unconstrained structure optimization was performed at the MM level. After optimization, the energy of the optimized conformation and the size of the cluster in which the conformation is located were considered to select the dominant conformation. If the number of conformations in the two clusters differed greatly, the minimum point in the largest cluster was selected as the dominant conformation. If the number of conformations in the two clusters was close, the global minimum point in the two conformations was selected as the dominant conformation.

[0024] In this invention, a significant difference in the number of conformations between two clusters means that the proportion of conformations in one cluster reaches more than 30% of the total number of conformations. In this case, the global minimum point is selected as the dominant conformation from the conformations that meet this condition. Conversely, a similar number of conformations between two clusters means that the proportion of conformations in any of the clusters does not reach 30% of the total number of conformations. In this case, the minimum point is preferentially selected from the largest cluster as the dominant conformation.

[0025] To ensure the conformational quality of both the host and guest molecules, this invention also performed MD simulations on both molecules separately, with calculation parameters essentially consistent with those used in the MD simulations of the aforementioned complex. Since the guest molecule system is relatively small, the cutoff distance for non-covalent interactions was set to... Its conformational space is also relatively small, so the MD simulation time was set to 2 ns. At the same time, the conformations generated by the MD simulation of the host molecule and the guest molecule were clustered and analyzed, consistent with the clustering analysis method of the host-guest complex, to screen out the dominant conformations of the host molecule and the guest molecule.

[0026] Furthermore, in step 5, the QM / CSM calculation method is as follows:

[0027] (1) Remove water molecules from the preferred conformations of the host-guest complex, host molecule and guest molecule, optimize layer by layer using the ONIOM2 method, optimize the guest molecule using the B3LYP / 6-31G(d) method, optimize CB using the semi-empirical PM6 method[7], and calculate the vacuum recombination free energy using the QM / CSM method.

[0028] (2) Solvent effect was calculated using the PCM method. The preferred conformation for removing water molecules was calculated on the M062X / 6-31G(d,p) basis set to obtain the solvent effect.

[0029] (3) The total recombination free energy in the solution is obtained by adding the vacuum recombination free energy to the solvent effect.

[0030] Compared with the prior art, the advantages of this invention are:

[0031] (1) The computational model established in this invention has a huge impact on the final result. For some guest molecules that have poor matching with the CB[7] cavity, Autodock molecular docking is used to screen out relatively stable initial conformations from the complex conformation space. The initial conformation space is fully searched, providing accurate structural data for subsequent energy calculations.

[0032] (2) This invention selects the most stable composite conformation from different trajectories of the same subject-object system at a more precise QM level. The calculation results are then compared with experimental values ​​for verification, making the selection of representative conformation space more sufficient, reasonable and accurate.

[0033] (3) In the calculation of QM / CSM, the calculation model established in this invention adopts the ONIOM2 calculation method to perform hierarchical optimization of the host molecule and the guest molecule. The rigid CB[7] host molecule adopts the PM6 semi-empirical optimization, and the flexible guest molecule adopts the more accurate DFT method. The ONIOM2 calculation method achieves the unity of calculation accuracy and calculation speed, and speeds up the calculation while maintaining accurate calculation accuracy.

[0034] (4) The calculation model established in this invention combines the MD method with the QM / CSM method. It can not only calculate the recombination free energy in vacuum, but also examine the solvent effect, making the calculation of the recombination free energy of the host-guest system more complete. The calculation method of this model can calculate the trend of recombination free energy and provide guidance for experiments. Attached Figure Description

[0035] Figure 1 The flowchart is for the MD / QM / CSM method of the present invention for screening representative conformations of the CB[7] subject-guest system.

[0036] Figure 2Different initial structures for several subject-object systems: CB[7] / S-METH, CB[7] / R-METH, CB[7] / FEN, CB[7] / PCP, CB[7] / COC, CB[7] / S-KET and CB[7] / R-KET.

[0037] Figure 3 Cluster analysis diagram of the optimal MD trajectories of several subject-object systems, namely CB[7] / S-METH, CB[7] / R-METH, CB[7] / FEN, CB[7] / PCP, CB[7] / COC, CB[7] / S-KET and CB[7] / R-KET.

[0038] Figure 4 The graph shows the change of the centroid distance between the host and guest systems over time in MD simulations for several host-guest systems, namely CB[7] / S-METH, CB[7] / R-METH, CB[7] / FEN, CB[7] / PCP, CB[7] / COC, CB[7] / S-KET and CB[7] / R-KET.

[0039] Figure 5 The diagrams show representative structures of several host-guest systems, including CB[7] / S-METH, CB[7] / R-METH, CB[7] / FEN, CB[7] / PCP, CB[7] / COC, CB[7] / S-KET and CB[7] / R-KET. Among them, beige represents the CB[7] composite system RC. MM Green represents RC QM The dashed and solid lines represent RC. MM and RC QM Intermolecular hydrogen bonds, unit:

[0040] Figure 6 Linear correlation analysis of composite free energy and experimental values ​​for several host-guest systems CB[7] / S-METH, CB[7] / R-METH, CB[7] / FEN, CB[7] / PCP, CB[7] / COC, CB[7] / S-KET and CB[7] / R-KET using the MD / QM / CSM method, where green represents the CB[7] composite system. Blue is Red indicates Black represents the linear correlation line. Detailed Implementation

[0041] The present invention will be further described below with reference to specific embodiments and accompanying drawings.

[0042] The experimental data are shown in Table 1, and the reference is (Chembiochem, 2017, 18(16): 1583-1588).

[0043] Example 1

[0044] The verification process of the subject-object system is as follows: Figure 1 As shown:

[0045] Step 1: For CB[7], the topology is obtained by using the GAFF force field. The RESP charge is calculated at the HF / 6-31G* level, with the GLYCAM04 force field and Amber99SB as the force field parameters.

[0046] Step 2: Five amino-containing guest molecules—METH (chiral), FEN, PCP, COC, and KET (chiral)—were geometrically optimized on Gussian09 using the B3LYP / 6-31G(d) method and basis sets. Then, the topological structures were obtained in the GAFF force field using the antechamber, tleap, and parmchk2 programs in Amber software, and the RESP charges of the guest molecules were calculated at the HF / 6-31G* level.

[0047] Step 3: Using Gaussview and Autodock software, the initial structure of the CB[7] host-guest system, after charge and topology calculations, and the amino-containing guest molecule, was constructed at a 1:1 ratio. The constructed initial structure is shown below. Figure 2 As shown, IC A 1,j Built by Gaussview, IC A 2,j Built by Autodock software (A = S-METH, R-METH, FEN, PCP, COC, S-KET, R-KET).

[0048] Step 4: Perform MD simulations on the initial structures of each CB[7] host-guest system with different orientations using Amber software. Dissolve the initial structure of the host-guest complex in TIP3P water, and set the buffer distance between the host-guest complex and the boundary of the periodic aqueous solution box to be... Long-range electrostatic interactions were handled using the PME method, with the cutoff distance for non-covalent interactions set as... MD simulations were performed under periodic boundary conditions, with the SHAKE algorithm used to constrain H-containing chemical bonds. First, constraints were applied to the host molecule for initial structural optimization to eliminate potential abnormal interactions between the solvent water molecules, the host molecule, and their complexes. Then, a second structural optimization was performed under unconstrained conditions. After optimization, the system was heated. Langevin kinetics were used to heat the system from 400 ps to 300 K. Following heating, NPT ensemble simulations were performed at 1 atm and 300 K, with an integration step size of 2 fs. Energy and structure were saved every 2 ps, resulting in a total MD time of 4 ns. A total of 2000 conformations were obtained from the MD trajectory.

[0049] Cluster analysis was performed on the conformations on each MD simulation trajectory, classifying them into 6–8 clusters by adjusting the RMSD threshold. The structures of the two largest clusters were selected for unconstrained structure optimization at the MM level. After optimization, the dominant complex conformation was selected by considering both the energy of the optimized conformation and the size of the cluster. If the number of conformations in the two clusters differed significantly, the minimum point from the largest cluster was chosen as the dominant complex conformation; if the number of conformations in the two clusters was similar, the global minimum point between the two conformations was selected as the dominant complex conformation. Water molecules were removed from the selected dominant complex conformation for subsequent QM / CSM calculations. The cluster analysis of the optimal MD trajectory for each system is shown below. Figure 3 As shown.

[0050] To ensure the conformational quality of both the host and guest molecules, this invention also performed MD simulations on both molecules separately, with calculation parameters essentially consistent with those used in the MD simulations of the aforementioned complex. Since the guest molecule system is relatively small, the cutoff distance for non-covalent interactions was set to... Its conformational space is also relatively small, so the MD simulation time was set to 2 ns. Simultaneously, cluster analysis was performed on the conformations generated by the MD simulations of the host and guest molecules, consistent with the cluster analysis method for complexes, to screen out the dominant conformations of the host and guest molecules. Water molecules were removed from the selected dominant conformations for subsequent QM / CSM calculations.

[0051] The distance between the centroids of the subject and object in a MD trajectory can measure the dynamic behavior of subject-object separation. In the optimal MD trajectory, the change of the subject-object centroid distance over time for each system is as follows: Figure 4As shown. The size of the centroid distance of the system is related to the size of the guest molecule; the smaller the guest molecule, the smaller the centroid distance. The fluctuation of the centroid distance is related to the degree of matching between the guest molecule and CB[7]. The more matched the guest molecule is with the CB[7] molecule, the smaller the centroid distance fluctuation and the more stable the system. The host-guest centroid distance can reflect the strength of nonbonding interaction. Generally speaking, a small centroid distance (but not too small) and small fluctuation indicate a strong host-guest nonbonding interaction. The order of centroid distance fluctuation is: FEN>R-KET>COC>S-METH>R-METH>S-KET>PCP. Around 100 ps, ​​the host-guest centroid distance changes significantly, indicating that the guest molecule has entered the CB[7] cavity to form a complex. After that, the host-guest centroid distance of S-METH, R-METH and CB[7] is in The small centroid distance and small fluctuation indicate that the non-bonded interaction is strong. However, the small centroid distance fluctuation of S-KET and CB[7] is due to the large volume of the composite part, which makes TΔS too positive, which is not conducive to the stability of the composite.

[0052] Step 5: The recombination free energy of the dominant conformations of the individual CB[7] host molecules, individual amino-containing guest molecules, and CB[7] complexes selected from the MD simulation and subjected to cluster analysis in Step 4 is calculated using the QM / CSM method. In calculating the vacuum recombination free energy, the ONIOM2 method is used for hierarchical optimization, the B3LYP / 6-31G(d) method is used for optimizing the guest molecules, and the semi-empirical PM6 method is used for optimizing CB[7]. The PCM method is used to calculate the solvent effect, and the dominant conformations are calculated on the MO62X / 6-31G(d,p) basis set to obtain the solvent effect. The vacuum recombination free energy and the solvent effect are added together to obtain the total recombination free energy in the solution. The most stable dominant conformation is selected. The recombination free energy of the most stable dominant conformation after calculation is compared with the experimental value for linear correlation analysis.

[0053] Step 6: Perform stability ranking and linear correlation analysis between the calculated composite free energy of the most stable dominant conformation and the experimental values. If there is no good linear correlation between the two or the relative stability order of each system is inconsistent with the experimental values, return to Step 3, use Guassview software to readjust the initial orientation between the host and guest components or re-select conformations from Autodock, construct new MD initial conformations, and repeat the above steps until the composite free energy of the selected dominant conformation of the complex has a good linear correlation with the experimental values ​​and the relative stability order of each system is consistent with the experimental values. Figure 5 For the representative conformations selected from each system, RC MM For RC structures, unconstrained structural optimization is performed at the MM level after MD simulation. QM For RC MM The RC structure optimized by QM.

[0054] The calculated and experimental values ​​of the composite free energy of the selected dominant conformation showed a good linear correlation, and the trend of the calculated values ​​was consistent with that of the experimental values. The results are shown in Table 1 (the last two columns are the calculated and experimental values ​​of the composite free energy). Therefore, it is considered that the dominant conformation selected from MD simulation is reasonable.

[0055] Table 1. Composite free energy calculated by ONIOM(B3LYP / PM6) / PCM method versus experimental value. Unit: Kcal / mol

[0056]

[0057] Table 2. Solvent effects and contributions of each component to vacuum free energy calculated by ONIOM (B3LYP / PM6) and PCM methods, unit: kcal / mol

[0058]

[0059] The calculated results for each system (Table 1) were compared with the experimental values ​​to explore the correlation between the calculations and the experiments. For example... Figure 6 As shown, the results indicate that the calculated results for each system correlate well with the experimental values ​​(R0). 2 =0.98). Non-bonded interactions and hydrophobic interactions are the two main reasons affecting the binding trend of CB[7]. The range of non-bonded interactions is between -32.91 and -57.24 Kcal / mol, and the range of solvent effects is between 3.84 and 16.24 Kcal / mol. Non-bonded interactions have a greater impact on the binding affinity of each system. Generally speaking, there is a competitive relationship between non-bonded interactions and hydrophobic interactions: one is enhanced while the other is weakened. Guest molecules can be divided into two groups. One group includes METH and FEN, whose CB[7] complexes have high stability. The other group includes PCP, COC and KET, whose CB[7] complexes have low stability. The non-bonded-hydrophobic interaction competition is reflected in both groups. However, this competition is not obvious between the two groups. This is mainly because the latter group of molecules and CB[7] have poor shape matching. The hydrophobic benzene ring part is not fully compounded, which weakens the non-bonded interaction and is not conducive to the solvent effect. Therefore, rigid CB[7] host molecules can selectively compound different guest molecules better. Figure 6 The green line in and blue line It was found that the nonbonding interactions of each system gradually weakened, and the solvent effect also gradually weakened. This is mainly related to the structural differences of each system. Regarding the METH and FEN system, the charged amino group is located around the C=O portal of CB[7] after recombination, resulting in a strong Coulomb interaction and a strong nonbonding interaction. Regarding the PCP and COC system, the charged amino group is located far from the CB[7] portal after recombination, and cannot generate a strong Coulomb interaction, resulting in a weak nonbonding interaction. Regarding the KET system, chlorobenzene is recombinated in the cavity, and the degree of matching with CB[7] is low, resulting in a weak nonbonding interaction. The solvent effect is related to the depth of the benzene ring recombination in the CB[7] cavity. The more tightly the benzene ring is recombinated, the less the hydrophobic benzene ring part is exposed in the solution, and the more its hydrophilic group is fully exposed in the solution, resulting in a better solvent effect. In the METH and FEN system, the benzene ring is tightly recombinated, and there is less hydrophobic part exposed in the solution, resulting in a better solvent effect. In the PCP, COC and KET systems, the benzene ring complex is not as deep as that of METH and FEN, resulting in more hydrophobic parts exposed in the solution and a poorer solvent effect.

[0060] Based on the original MD / QM / CSM method, the present invention expands the search space of the subject-guest composite conformation by using the Autodock program and optimizes the screening method of representative conformations, making the calculation results more reliable; and in the QM calculation, a hierarchical calculation is adopted, which improves the calculation speed without losing the calculation accuracy. The comparison between the experimental values ​​and the calculation results shows that the results of the CB[7] subject-guest system calculated by the MD / QM / CSM method are reliable and can provide theoretical guidance for the experiment.

Claims

1. A MD / QM / CSM method for screening representative conformations of CB[7] host-guest systems, characterized in that, Comprising the following steps: Step 1, the host molecule CB[7] is adopted GAFF force field, and the topological structure of CB[7] and the RESP charge calculated at the HF / 6-31G* level are obtained; Step 2, the amino-containing guest molecule is adopted GAFF force field, and the topological structure of the guest molecule and the RESP charge calculated at the HF / 6-31G* level are obtained; Step 3, the host molecule after charge and topological calculation in step 1 and the guest molecule after charge and topological calculation in step 2 are adopted Gaussview to build two different orientations of initial conformation, and Autodock software is adopted to build two different orientations of initial conformation, and four different orientations of initial conformation are built; Step 4, four different orientations of initial conformation are subjected to MD simulation, four MD simulation trajectories are obtained, and the threshold is adjusted, and the cluster analysis and structure optimization of the complex conformation in each MD trajectory are carried out, so as to screen out the dominant conformation of each trajectory complex; At the same time, the host molecule in step 1 and the guest molecule in step 2 are subjected to MD simulation, the threshold is adjusted, the cluster analysis and structure optimization are carried out, and the dominant conformation of the single host molecule and the single guest molecule is screened out, which is used for subsequent QM / CSM calculation of the single host molecule and the single guest molecule; Step 5, the dominant conformation of each trajectory complex obtained in step 4 is subjected to QM / CSM calculation, and the most stable complex dominant conformation is screened out, and the QM / CSM calculation method is as follows: (1) remove the water molecules in the dominant conformation of the host-guest complex, the host molecule and the guest molecule, and adopt ONIOM2 method for hierarchical optimization, use B3LYP / 6-31G(d) method to optimize the guest molecule, use semi-empirical PM6 method to optimize CB[7], and adopt QM / CSM method to calculate the vacuum complex free energy; (2) using PCM method to calculate the solvent effect, removing the water molecules in the dominant conformation, and calculating at M062X / 6-31G(d,p) basis set to obtain the solvent effect; (3) adding the vacuum complex free energy and the solvent effect to obtain the total complex free energy in the solution; Step 6, the calculation results of the most stable complex dominant conformation obtained in step 5 are compared with the experimental values, the calculated complex free energy is linearly correlated with the experimental values, if the relative stability of the host-guest system is consistent with the experimental results, and the calculated value has good linear correlation with the experimental value, then the most stable complex dominant conformation is the accurate representative conformation of the system; If the relative stability of the host-guest system is inconsistent with the experimental results, or the calculated value has no good linear correlation with the experimental value, then return to step 3.

2. The MD / QM / CSM method of claim 1, wherein, The guest molecule is selected from METH, FEN, PCP, COC or KET, and the structural formula of each guest molecule is:

3. The MD / QM / CSM method of claim 1, wherein, In step 1, GLYCAM04 force field and Amber99SB are used as force field parameters.

4. The MD / QM / CSM method of claim 1, wherein, In step 3, the host molecule CB[7] and the guest molecule are constructed into the initial conformation of the host-guest complex at a ratio of 1:

1.

5. The MD / QM / CSM method of claim 1, wherein, In step 4, the specific method for screening out the dominant conformation of each trajectory complex is as follows: (1) The initial structure of the host-guest complex was dissolved in TIP3P water, and the buffer distance of the host-guest complex to the periodic water solution box boundary was set to The long-range electrostatic interaction was treated using the PME method, and the cutoff distance of the non-covalent bond interaction was set to The MD simulation was performed under periodic boundary conditions, and the SHAKE algorithm was used to constrain the H-containing chemical bonds. First, the host molecule was subjected to a constraint force, and the system was subjected to a first-step structure optimization to eliminate possible abnormal interactions between the solvent water molecules and the host molecules and their complexes. Then, the system was subjected to a second structure optimization under unconstrained conditions. After optimization, the system was warmed up, and the Langevin dynamics was used to warm up the system to 300 K for 400 ps. After warming up, the NPT ensemble was used to perform MD simulation at 1 atm and 300 K. The integration step was set to 2 fs, and the energy and structure were saved every 2 ps. The total MD time was 4 ns, and a total of 2000 conformations were obtained in the MD trajectory. (2) Cluster analysis is performed on the conformations generated by MD simulation, the conformations are divided into 6-8 classes by adjusting the RMSD threshold, the structures of the first two largest clusters are screened, and unconstrained structure optimization is performed at the MM level. After optimization, the energy of the optimized conformation and the size of the cluster where the conformation is located are considered comprehensively, and the dominant conformation is selected. If the number of conformations in the two clusters is quite different, the minimum point in the largest cluster is selected as the dominant conformation. If the number of conformations in the two clusters is similar, the global minimum point in the two classes of conformations is selected as the dominant conformation.

6. The MD / QM / CSM method of claim 5, wherein, The specific method for screening the dominant conformation of the single host molecule is as follows: (1) The initial structure of the host molecule was dissolved in TIP3P water, and the buffer distance of the host-guest complex to the periodic water solution box boundary was set to The PME method was used to process long-range electrostatic interactions, and the cutoff distance of non-covalent bond interaction was set to MD simulation was performed under periodic boundary conditions, and SHAKE algorithm was used to constrain H-containing chemical bonds. First, the host molecule was subjected to a constraint force, and the system was optimized in the first step to eliminate possible abnormal interactions between solvent water molecules and the host molecule and its complex. Then, the structure was optimized again under unconstrained conditions. After optimization, the system was warmed up, and Langevin dynamics was used to warm up the system to 300 K for 400 ps. After warming up, NPT ensemble was used to perform MD simulation at 1 atm and 300 K. The integration step was set to 2 fs, and the energy and structure were saved every 2 ps. The total MD time was 4 ns, and a total of 2000 conformations were obtained in the MD trajectory. (2) Cluster analysis is performed on the conformations generated by MD simulation, the conformations are divided into 6-8 classes by adjusting the RMSD threshold, the structures of the first two largest clusters are screened, and unconstrained structure optimization is performed at the MM level. After optimization, the energy of the optimized conformation and the size of the cluster where the conformation is located are considered comprehensively, and the dominant conformation is selected. If the number of conformations in the two clusters is quite different, the minimum point in the largest cluster is selected as the dominant conformation. If the number of conformations in the two clusters is similar, the global minimum point in the two classes of conformations is selected as the dominant conformation.

7. The MD / QM / CSM method of claim 5, wherein, The specific method for screening the dominant conformation of the single guest molecule is as follows: (1) The initial structure of the guest molecule was dissolved in TIP3P water, and the buffer distance of the host-guest complex to the periodic water solution box boundary was set to The long-range electrostatic interactions were treated using the PME method, and the cutoff distance for non-bonded interactions was set to The MD simulation was performed under periodic boundary conditions and the SHAKE algorithm was used to constrain the H-containing chemical bonds. First, the host molecule was subjected to a constraint force, and the system was optimized in the first step to eliminate any abnormal interactions that may exist between the solvent water molecules and the host molecule and its complex. Then, a second structure optimization was performed under unconstrained conditions. After optimization, the system was warmed up using Langevin dynamics to 300 K over 400 ps. After warming up, the MD simulation was performed under the NPT ensemble at 1 atm and 300 K, with an integration step size of 2 fs. The energy and structure were saved every 2 ps, and the total MD time was 2 ns. A total of 2000 conformations were obtained in the MD trajectory. (2) Cluster analysis is performed on the conformations generated by MD simulation, the conformations are divided into 6-8 classes by adjusting the RMSD threshold, the structures of the first two largest clusters are screened, and unconstrained structure optimization is performed at the MM level. After optimization, the energy of the optimized conformation and the size of the cluster where the conformation is located are considered comprehensively, and the dominant conformation is selected. If the number of conformations in the two clusters is quite different, the minimum point in the largest cluster is selected as the dominant conformation. If the number of conformations in the two clusters is similar, the global minimum point in the two classes of conformations is selected as the dominant conformation.

8. The MD / QM / CSM method according to any one of claims 5 to 7, characterized in that, The number of conformations in the two clusters is quite different, which means that the proportion of the number of conformations in the cluster reaches more than 30% of the whole, and the global minimum point in the conformation that meets the condition is selected as the dominant conformation. The number of conformations in the two clusters is similar, which means that the proportion of the number of conformations in all clusters does not reach 30% of the whole, and the minimum point in the largest cluster is selected as the dominant conformation.

Images