A GPU-based fast molecular docking method

Abstract

The present invention provides a GPU-based rapid molecular docking method, the purpose of which is to further improve the docking accuracy and speed of the QuickVina2-GPU method. In improving the docking accuracy, the Gypsum-DL method is used to optimize the input small molecule structure; the cross-docking method is used to select the optimal receptor structure; and the Coach-D method is used to optimize the docking pocket of the receptor. In improving the docking speed, the global grid cache is used to accelerate the intermolecular energy calculation of QuickVina2-GPU. In the preparation of the small molecule library, the present invention makes up for the shortcomings of Openbabel in handling small molecule charges, selects receptors with good docking accuracy and can effectively restore the true conformation, reduces the number of data transmissions between the CPU and the GPU, greatly improves the speed of molecular docking, and reduces communication overhead. The present invention can further meet the infinite pursuit of docking speed and accuracy under large-scale drug virtual screening.

Description

Technical Field

[0001] The present invention relates to a GPU-based rapid molecular docking method, belonging to the technical field of software and hardware collaborative accelerated drug design. Background Art

[0002] Molecular docking is a computer simulation method used to predict the binding mode between small molecule ligands and macromolecular receptors (usually proteins or nucleic acids). It samples and scores different receptor-ligand complex conformations to identify the binding state with the lowest interaction energy and the most stable binding state. This technique can be used to validate experimentally derived receptor-ligand conformations and for high-throughput virtual screening of compound libraries to rapidly identify potential first-choice ligands. Molecular docking is also widely used in enzyme engineering, antibody engineering, and protein structure prediction. QuickVina2 is an upgraded version of the well-known protein-small molecule docking software AutoDock Vina. Compared to AutoDock Vina, it significantly improves computational speed and result reliability. Specifically, QuickVina2 uses heuristic algorithms to avoid unnecessary local searches, significantly reducing computation time. Compared to its predecessor, QuickVina1, QuickVina2 ensures necessary local searches by determining whether random docking results require local optimization. As a result, QuickVina2 can more reliably identify the docking conformation with the global lowest energy. Although QuickVina2 has improved its computing speed through algorithm optimization, the CPU computing bottleneck makes it difficult to meet the speed requirements of modern drug discovery.

[0003] Therefore, we designed a GPU-based rapid molecular docking method to further enhance the docking accuracy and speed of QuickVina2-GPU. Through three strategies to improve docking accuracy, we can obtain highly accurate ligand-receptor binding conformations. By incorporating a global grid cache, we can rapidly dock ligands and receptors for high-throughput virtual screening, reducing the cost of drug development. Summary of the Invention

[0004] In view of the problems existing in the prior art, the present invention provides a GPU-based rapid molecular docking method.

[0005] In order to achieve the above object, the technical solutions adopted by the present invention are as follows:

[0006] The Gypsum-DL method is used to optimize the input small molecule structure; the cross-docking method is used to select the optimal receptor structure; the Coach-D method is used to optimize the docking pocket of the receptor; and the global grid cache is used to accelerate the intermolecular energy calculation of QuickVina2-GPU.

[0007] As a further improvement of the present invention, the Gypsum-DL method comprises the following steps:

[0008] S1. Prepare the small molecule file to be converted;

[0009] S2. Write the configuration file according to the .json configuration format requirements;

[0010] S3, generate and optimize multiple conformations of small molecules;

[0011] S4. Energy evaluation and ranking of conformations from the same small molecule are performed, and the one with the lowest energy is output as an .sdf file with a three-dimensional structure;

[0012] S5. Use a script to convert the output .sdf file with the three-dimensional structure into the .pdbqt format required for docking.

[0013] As a further improvement of the present invention, in S2, the .json configuration format specifically includes: the path of the small molecule file to be converted; the working mode, including one of multi-process, serial, and high-performance computing; the output folder path; the maximum number of output conformations; the maximum pH, and the minimum pH.

[0014] As a further improvement of the present invention, S3 specifically includes the following steps:

[0015] S3.1. Desalt small molecules to eliminate molecular fragments irrelevant to docking;

[0016] S3.2, hydrogenate small molecules and adjust pH;

[0017] S3.3. Enumerate tautomers of small molecules;

[0018] S3.4. Enumerate small molecule chirality and double bonds;

[0019] S3.5. Generate and optimize 3D conformations of small molecules.

[0020] As a further improvement of the present invention, the S4 specifically includes the following steps:

[0021] S4.1. Use the Chem.AddHs() method in rdkit to implicitly add hydrogen to the molecule.

[0022] S4.2. Apply the MMFF force field to the small molecule using the AllChem.MMFFGetMoleculeForceField() method in rdkit.Chem.

[0023] S4.3. Use the ff.CalcEnergy() method to calculate the energy of the small molecule and determine the corresponding small molecule for output based on the minimum value found.

[0024] As a further improvement of the present invention, the cross-docking method specifically comprises the following steps:

[0025] S6. Align and separate the prepared protein-ligand structures;

[0026] S7. Each ligand is docked with different protein ligands, and the error between the top1score conformation and the ligand conformation under the protein crystal structure is calculated, expressed as the root mean square error, recorded as RMSD. If RMSD < 2, it indicates that the docking successfully restores the true conformation;

[0027] S8. The docking acceptor was determined by the number of RMSD < 2 and the average RMSD in each protein crystal docking result.

[0028] As a further improvement of the present invention, the Coach-D method specifically comprises the following steps:

[0029] S9. Upload the receptor and optional ligand to be predicted on the Coach-D server website and wait for the server to return the results;

[0030] S10. After obtaining the results, find the predicted binding residues of the highest-ranked pocket;

[0031] S11. Calculate the geometric center of the binding residue as the center point of the docking pocket, and calculate the distance from this center point to the farthest residue in the XYZ direction as the size of the pocket to contain all the binding residues.

[0032] As a further improvement of the present invention, the method of accelerating the intermolecular energy calculation of QuickVina2-GPU by using the global grid cache specifically includes the following steps:

[0033] S12, read receptor, ligand files and operation parameter files;

[0034] S13. Create and initialize the OpenCL platform;

[0035] S14, generating an intramolecular energy table;

[0036] S15, the GPU-side kernel 1 generates a global grid cache and stores it in the global memory area;

[0037] S16, randomize the initial conformation and perform iterative search;

[0038] S17. After the CPU further optimizes, scores, and sorts the conformations, the next ligand repeats the operation of S16.

[0039] As a further improvement of the present invention, the S15 specifically includes the following steps:

[0040] S15.1. Divide the docking pocket into multiple squares. The size of each square after division is After obtaining the quantized coordinates of each grid point after division, the grid point thread is started in GPU core 1. Each thread calculates the force between the atom and the receptor at its corresponding grid point.

[0041] S15.2. Convert the quantized coordinates of the grid points into their real coordinates in the receptor according to the provided docking pocket file;

[0042] S15.3. Divide a cube with the grid point as the center. The default size of the cube is

[0043] S15.4: If the cube can contain the atoms in the receptor, and the distance between the receptor atoms and the lattice atoms is less than Then calculate the interaction force between the atom at this lattice point and these included receptor atoms. 17 interaction forces are calculated for this lattice point for subsequent use.

[0044] S15.5: After all threads complete the above operations, 17 global grid caches are obtained according to the atomic type. The global grid cache is stored in the global memory and is shared by all GPU threads.

[0045] As a further improvement of the present invention, the S16 specifically includes the following steps:

[0046] S16.1, perturbing random molecular conformations;

[0047] S16.2. For a randomly perturbed object, if the heuristic condition is not met, continue perturbing; if the condition is met, perform further local search;

[0048] S16.3, locally optimize the conformation and determine whether to update the conformation based on the acceptance criterion;

[0049] S16.4. Determine whether the number of cycles of S16.1, S16.2, and S16.3 reaches the search step number;

[0050] S16.5. If the number of steps is met, the optimal configuration of each thread in GPU core 2 is transferred to the CPU-side memory.

[0051] Beneficial effects:

[0052] 1. In the preparation of small molecule libraries, Gypsum-DL was used to optimize receptors, which made up for the shortcomings of Openbabel in handling small molecule charges;

[0053] 2. In the selection of receptors, a cross-docking strategy was used to select receptors with good docking accuracy and that could effectively restore the true conformation;

[0054] 3. In terms of pocket selection, Coach-D is used to predict and optimize pockets, further selecting the best docking pocket to improve docking accuracy;

[0055] 4. The intermolecular energy calculation of QuickVina2-GPU is accelerated by using the global grid cache, which reduces the number of communications between the device and the host, and greatly improves the speed of molecular docking. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 This is the overall architecture diagram of the present invention;

[0057] Figure 2 This is the flow chart of S1;

[0058] Figure 3 This is a flowchart of the QuickVina2-GPU method;

[0059] Figure 4 is a flow chart of S15;

[0060] Figure 5 Flowchart for generating the optimal configuration for each thread of core 2 of the GPU in the present invention. DETAILED DESCRIPTION

[0061] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail with reference to the following specific embodiments and accompanying drawings, but this does not limit the scope of the present invention.

[0062] It should be noted that the GPU model used in the following embodiments is NVIDIA GeForce RTX 30606GB, and the CPU model is Intel(R) Core(TM) i7-12700H CPU, but the present invention is not limited to this.

[0063] refer to Figure 1 , which is the overall architecture diagram of the present invention.

[0064] refer to Figure 2 , a flow chart of S2 of a GPU-based rapid molecular docking method of the present invention, comprising the following steps: optimizing the input small molecule structure using the Gypsum-DL method, including:

[0065] Step 1: Prepare the small molecule file to be converted.

[0066] Table 1 shows some examples of the .smi format of small molecule files:

[0067] Table 1

[0068]

[0069] Table 2 shows some examples of the small molecule file .sdf format:

[0070] Table 2

[0071]

[0072]

[0073] Step 2: Write the .json file,

[0074] Table 3 shows an example of a .json file:

[0075] Table 3

[0076]

[0077] The parameters in Table 3 include:

[0078] source: the path of the small molecule to be converted;

[0079] seprate_output_files: separation output, one output for each small molecule;

[0080] job_manager: working mode;

[0081] output_folder: output folder;

[0082] add_pdb_output: add pdb output;

[0083] max_variants_per_compound: The maximum number of conformations output for each small molecule;

[0084] thoroughness: thoroughness. The larger the setting, the more conformations will be enumerated during the enumeration process. num_processors: the number of processors.

[0085] min_ph: minimum pH value;

[0086] max_ph: maximum pH value;

[0087] Step 3, generating and optimizing multiple conformations of the small molecule, includes:

[0088] Step 3.1: Desalt the small molecules to exclude molecular fragments not relevant to docking;

[0089] Step 3.2: Hydrogenate the small molecule so that the small molecule is within the user-specified pH range;

[0090] Step 3.3: Enumerate the tautomers of the small molecule;

[0091] Step 3.4: Enumerate small molecule chirality and double bonds;

[0092] Step 3.5: Generate and optimize the 3D conformation of the small molecule;

[0093] Step 4: Energy evaluation and ranking of conformations from the same small molecule, and taking the output with the lowest energy. This step 4 includes:

[0094] Step 4.1: Use Chem.AddHs() in rdkit to implicitly add hydrogen to the molecule;

[0095] Step 4.2: Use AllChem.MMFFGetMoleculeForceField() in rdkit.Chem to apply the MMFF force field to the small molecule.

[0096] Step 4.3: Use the ff.CalcEnergy() method to calculate the energy of the small molecule and store the result in a dictionary. Use Python's built-in min function to find the minimum value of the values ​​field in the dictionary. Use the minimum value found to determine the corresponding small molecule for output.

[0097] Step 5: Use a script to convert the output .sdf file to .pdbqt format.

[0098] Table 4 and Table 5 are examples of converted pdbqt files:

[0099] Table 4

[0100]

[0101]

[0102] Table 5

[0103]

[0104] The cross-docking method was used to select suitable receptors, including:

[0105] Step 6: The eight protein-ligand structures of 4ITH, 4ITH, 6C3E, 6C4D, 5TX5, 6NW2, 6HHO, and 6R5F were aligned and separated from their ligands. The small molecule ligands were then not processed to maintain their alignment.

[0106] Step 7: Each ligand is then docked with a different protein ligand, and the top 1 scoring conformation is selected and analyzed for RMSD to evaluate their docking quality. Generally speaking, an RMSD < 2 indicates that the docking has successfully restored the true conformation.

[0107] Step 8: Determine the docked protein crystals by the number of RMSD < 2 in each protein crystal docking result and the average RMSD. The results are shown in Table 6. For RIPK1, we selected 5TX5 as the molecular docking target.

[0108] Table 6

[0109]

[0110] The docking pocket of the receptor was selected by Coach-D prediction, including:

[0111] Step 9: Upload the receptor 5TX5 and ligand to be predicted on the Coach-D website and wait for the server’s results;

[0112] Step 10: Get the results and find the predicted binding residues of the highest-ranked pockets;

[0113] Step 11: Calculate the geometric center of the binding residue as the center of the docking pocket; calculate the distance from this center point to the farthest residue in the XYZ direction as the size of the pocket to contain all the binding residues.

[0114] The small molecule libraries were converted to pdbqt and then sent to QuickVina2-GPU for molecular docking;

[0115] Figure 3 The flowchart for this step includes:

[0116] Step 12: Read the receptor, ligand files, and run parameter files, including the following files:

[0117] Tables 7 and 8 show some of the receptor files for 5TX5:

[0118] Table 7

[0119]

[0120]

[0121] Table 8

[0122]

[0123] Table 9 shows the docking pocket and other parameter files:

[0124] Table 9

[0125]

[0126] The parameters in Table 9 include:

[0127] receptor: the path where the receptor file of 5TX5 is located;

[0128] ligand_directory: the directory where all ligand files of 5TX5 are located;

[0129] output_directory: the path where the output files of 5TX5 are located;

[0130] center_x, center_y, center_z: the three-dimensional coordinates of the center site of the docking pocket; size_x, size_y, size_z: the size of the docking pocket;

[0131] thread: number of threads, the default value is 8000;

[0132] Step 13: Create and initialize the OpenCL platform, including the following:

[0133] Query the platform and use the clGetPaltformIDs() function in OpenCL to obtain the number of available OpenCL platforms. Then, call this function again based on the obtained number of platforms to initialize an available platform.

[0134] Query the device and use the clGetDeviceIDs() function in OpenCL to obtain the number of devices in the platform. Then call the function again based on the obtained number of devices to initialize an available device.

[0135] Create a context. The context contains memory, programs, kernel objects, etc. Use the clCreateContext() function to create a context for the specified device, which is used for communication between the device and the program, and between devices.

[0136] Create a command queue. The command queue can operate multiple objects in the context and use the cl_queue_properties() property value to create a command queue.

[0137] Create a program object, use the clCreateProgramWithSource() function to create a program object from the source code, and use the clBuildProgram() function to build the program object for the specified device;

[0138] Create kernel objects and use the clCreateKerne() function to create two kernel objects;

[0139] Set kernel parameters using the clSetKernelArg() function to enable execution of the kernel.

[0140] Execute the kernel. The clEnqueueNDRangeKernel() function can queue the kernel to be executed on the device using the command queue.

[0141] Step 14: Generate an intramolecular energy table for subsequent calculation of intramolecular energy. Generate a random number lookup table based on the random number seed for random perturbation of the molecular conformation. Transfer these two tables to the GPU.

[0142] Step 15: Kernel 1 on the GPU side generates a global grid cache and stores it in the global memory area.

[0143] Figure 4 This is a flowchart of step 15, which includes:

[0144] Step 15.1: Divide the docking pocket into multiple small grids. The size of the small grids after division is After obtaining the quantized coordinates of the vertices of each small grid (hereinafter referred to as grid points), GPU core 1 also starts several grid point threads, each of which calculates the forces acting on the atoms at its corresponding grid point. For each grid point in each thread, the following operations are performed:

[0145] Step 15.2: Convert the quantized coordinates of the grid points into their real coordinates in the receptor according to the provided docking pocket file.

[0146] Tables 10 and 11 show the quantized coordinates and true coordinates of some grid points in the case of docking pockets.

[0147] Table 10

[0148]

[0149] Table 11

[0150]

[0151] Step 15.3: Divide a cube with the grid point as the center. The default size of the cube is

[0152]

[0153] Step 15.4: If the cube can contain the atoms in the receptor, and the distance between the receptor atoms and the lattice atoms is less than Then calculate the interaction force between the atom at the lattice point and the included receptor atoms. Since there may be 17 different types of atoms at the lattice point, 17 interaction forces can be calculated for the lattice point for subsequent use.

[0154] Tables 12 and 13 show the values ​​of 17 interaction forces obtained for atoms at some lattice points:

[0155] Table 12

[0156]

[0157]

[0158] Table 13

[0159]

[0160] Step 15.5: After all threads complete the above operations, 17 global grid caches are obtained according to the atomic type. The global grid cache is stored in global memory and is shared by all GPU threads.

[0161] Step 16: This step is serial for small molecule libraries, but parallel for single molecules. The CPU randomizes the conformation of a single ligand into many different molecular conformations and transmits them to Kernel 2 on the GPU. Kernel 2 on the GPU uses the random molecular conformations generated by the CPU to perform an iterative Monte Carlo local search to obtain the top-scoring conformations. The score obtained for a conformation is positively correlated with the total energy of the conformation; that is, the lower the total energy of the conformation, the lower the score (the lower the better). The total energy of the conformation includes intermolecular energy and intramolecular energy. The intermolecular energy is calculated from the global grid cache obtained by Kernel 1 on GPU, and the intramolecular energy is calculated from the intramolecular energy table. Kernel 2 on the GPU transmits the final conformation back to the CPU.

[0162] The molecular conformation is represented by the formula:

[0163] C: molecular conformation,

[0164] x, y, z: the position of the molecular conformation in the search space,

[0165] a, b, c, d: quaternions of molecular conformation,

[0166] The torsion angle of the rotatable bond,

[0167] N rot : The number of rotatable bonds.

[0168] Figure 5 The flowchart for generating the best conformation for each thread in core 2 of the GPU includes the following steps:

[0169] Step 16.1: Perturb the random molecular conformation by adding probability-based perturbations to the following three sets of data: perturb the three-dimensional coordinates of the molecular conformation in the search space with a uniformly distributed random number; perturb the quaternion of the molecular conformation with a normally distributed random number; and replace a specific dimension of the torsion angle of a rotatable bond in the molecular conformation with a uniformly distributed random number.

[0170] Step 16.2: For the randomly perturbed conformation, the thread uses the conformations and heuristic conditions stored by other threads to detect. If the heuristic conditions are not met, it returns to continue to randomly perturb the conformation. If the conditions are met, it continues to perform further local search. This can reduce unnecessary local searches and improve search efficiency.

[0171] Step 16.3: Perform BFGS local optimization on the conformation and use the acceptance criterion to determine whether to accept the updated conformation.

[0172] Step 16.4: Determine whether the number of cycles of steps 16.1, 16.2, and 16.3 reaches the number of search steps.

[0173] Step 16.5: If the set number of search steps is met, the optimal configuration obtained by each thread in GPU core 2 is transferred to the memory on the CPU side.

[0174] Step 17: After further optimization, scoring, and sorting of the conformations on the CPU side, multiple random conformations of the next ligand are transferred to the GPU side for step 8.5.

[0175] Tables 14 and 15 show the optimal conformation data obtained for the ligand file:

[0176] Table 14

[0177]

[0178] Table 15

[0179]

[0180]

[0181] REMARK VINA RESULT: -5.4 0.000 0.000. This row is the docking result. -5.4 indicates the binding affinity of the molecular conformation.

[0182] In summary, in the preparation of the small molecule library, Gypsum-DL was used to optimize the receptors, which made up for the shortcomings of Openbabel in handling small molecule charges. In the selection of receptors, the cross-docking strategy was used to select receptors with good docking accuracy and the ability to effectively restore the true conformation. In the selection of pockets, Coach-D was used to predict and optimize pockets, and the best docking pocket was further selected to improve the docking accuracy. The global grid cache was used to accelerate the intermolecular energy calculation of QuickVina2-GPU, reducing the number of communications between the device and the host, and greatly improving the speed of molecular docking.

[0183] The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A GPU-based rapid molecular docking method, characterized by: The Gypsum-DL method is used to optimize the input small molecule structure; the cross-docking method is used to select the optimal receptor structure; the Coach-D method is used to optimize the docking pocket of the receptor; and the global grid cache is used to accelerate the intermolecular energy calculation on QuickVina2-GPU. The Gypsum-DL method comprises the following steps: S1. Prepare the small molecule file to be converted; S2. Write the configuration file according to the .json configuration format requirements; S3, generate and optimize multiple conformations of small molecules; S4. Energy evaluation and ranking of conformations from the same small molecule are performed, and the one with the lowest energy is output as an .sdf file with a three-dimensional structure; S5. Use a script to convert the output .sdf file with the three-dimensional structure into the .pdbqt format required for docking; The cross-docking method specifically comprises the following steps: S6. Align and separate the prepared protein-ligand structures; S7. Each ligand is docked with different protein ligands, and the error between the top1score conformation and the ligand conformation under the protein crystal structure is calculated, expressed as the root mean square error, recorded as RMSD. If RMSD < 2, it indicates that the docking successfully restores the true conformation; S8. The docking receptor is determined by the number of RMSD < 2 and the average RMSD in each protein crystal docking result. The Coach-D method specifically comprises the following steps: S9. Upload the receptor and optional ligand to be predicted on the Coach-D server website and wait for the server to return the results; S10. After obtaining the results, find the predicted binding residues of the highest-ranked pocket; S11. Calculate the geometric center of the binding residue as the center of the docking pocket, and calculate the distance from this center to the farthest residue in the XYZ direction as the size of the pocket to contain all the binding residues; The method of using the global grid cache to accelerate the intermolecular energy calculation of QuickVina2-GPU specifically includes the following steps: S12, read receptor, ligand files and running parameter files; S13. Create and initialize the OpenCL platform; S14, generating an intramolecular energy table; S15, the GPU-side kernel 1 generates a global grid cache and stores it in the global memory area; S16, randomize the initial conformation and perform iterative search; S17. After the CPU further optimizes, scores, and sorts the conformations, the next ligand repeats the operation of S16.

2. The GPU-based rapid molecular docking method according to claim 1, characterized in that: In S2, the .json configuration format specifically includes: the path of the small molecule file to be converted; the working mode, including one of multi-process, serial, and high-performance computing; the output folder path; the maximum number of output conformations; the maximum pH, and the minimum pH.

3. The GPU-based rapid molecular docking method according to claim 1, characterized in that: The S3 specifically includes the following steps: S3.

1. Desalt small molecules to eliminate molecular fragments irrelevant to docking; S3.2, hydrogenate small molecules and adjust pH; S3.

3. Enumerate tautomers of small molecules; S3.

4. Enumerate small molecule chirality and double bonds; S3.

5. Generate and optimize 3D conformations of small molecules.

4. The GPU-based rapid molecular docking method according to claim 1, characterized in that: The S4 specifically includes the following steps: S4.

1. Use the Chem.AddHs() method in rdkit to implicitly add hydrogen to the molecule. S4.

2. Apply the MMFF force field to the small molecule using the AllChem.MMFFGetMoleculeForceField() method in rdkit.Chem. S4.

3. Use the ff.CalcEnergy() method to calculate the energy of the small molecule and determine the corresponding small molecule by finding the minimum value for output.

5. The GPU-based rapid molecular docking method according to claim 1, characterized in that: The S15 specifically includes the following steps: S15.

1. Divide the docking pocket into multiple squares. The size of each square is 0.44 ×0.44 ×0.44 , get the quantized coordinates of each grid point after division, start the grid point thread in GPU core 1, and each thread calculates the force between the atom and the receptor at its corresponding grid point; S15.

2. Convert the quantized coordinates of the grid points into their real coordinates in the receptor according to the provided docking pocket file; S15.

3. Divide a cube with the grid point as the center. The default size of the cube is 8 ×8 ×8 ; S15.4: If the cube can contain the atoms in the receptor, and the distance between the receptor atoms and the lattice atoms is less than 8 , then calculate the interaction forces between the atoms at this lattice point and these included receptor atoms, and 17 interaction forces are calculated for this lattice point for subsequent use; S15.5: After all threads complete the above operations, 17 global grid caches are obtained according to the atomic type. The global grid cache is stored in the global memory and shared by all GPU threads.

6. The GPU-based rapid molecular docking method according to claim 1, characterized in that: The S16 specifically includes the following steps: S16.1, perturbing random molecular conformations; S16.

2. For a randomly perturbed object, if the heuristic condition is not met, continue perturbing; if the condition is met, perform further local search; S16.3, locally optimize the conformation and determine whether to update the conformation based on the acceptance criterion; S16.

4. Determine whether the number of cycles of S16.1, S16.2, and S16.3 reaches the search step number; S16.

5. If the number of steps is met, the optimal configuration of each thread in GPU core 2 is transferred to the CPU-side memory.

Images