Method, device, design method and device for optimizing force field and non-bonded parameters thereof

Abstract

The present application relates to a method, device, design method and device for optimizing a force field and its non-bonded parameters. The method for optimizing the non-bonded parameters of the force field includes: obtaining M clusters; determining the quantitative intermolecular interaction energy of each polymer of at least some of the M clusters; for each cluster in the M clusters, fitting the electrostatic potential of the central molecule of the cluster based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then obtaining the charge parameters of the force field based on the atomic charge of each atom in the central molecule of the M clusters; based on the quantitative intermolecular interaction energy of each polymer of the clusters, using the atomic charge of each atom in the central molecule as the initial charge of each atom of the molecule in the polymer, fixing the bonding parameters and the initial charge unchanged, and fitting to obtain the van der Waals parameters of the molecule. The present application can improve the accuracy of the description of polarization by force field parameters.

Description

Technical Field

[0001] The present application relates to the field of computer simulation technology, and in particular to a method, device, design method and device for optimizing a force field and its non-bonded parameters. Background Art

[0002] With the rapid development of computer technology and artificial intelligence technology, computer simulation technology is being applied to more and more scenarios, such as material design, drug design, etc.

[0003] However, the applicant has found that the accuracy of the force field parameters obtained by the related technology needs to be improved. Summary of the invention

[0004] In order to solve or partially solve the problems existing in the related art, the present application provides a method, device, design method and device for optimizing a force field and its non-bonded parameters, which can effectively improve the accuracy of the obtained force field parameters.

[0005] The first aspect of the present application provides a method for optimizing the non-bonded parameters of a force field, comprising: obtaining M clusters, each cluster in the M clusters including N molecules, M and N being integers greater than or equal to 1; determining the quantified intermolecular interaction energy of each polymer of at least some of the M clusters; for each cluster in the M clusters, fitting the electrostatic potential of the central molecule of the cluster based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then obtaining the charge parameters of the force field based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is a charge determined by the atomic charge of each molecule in the cluster except the central molecule; based on the quantified intermolecular interaction energy of each polymer of the cluster, taking the atomic charge of each atom in the central molecule as the initial charge of each atom in the molecule in the polymer, fixing the bonding parameters and the initial charge unchanged, and fitting to obtain the van der Waals parameters of the molecule.

[0006] The second aspect of the present application provides a method for optimizing force field parameters, comprising: determining bonding parameters, and determining non-bonding parameters based on the above method; and obtaining force field parameters based on the bonding parameters and the non-bonding parameters.

[0007] The third aspect of the present application provides a method for calculating crystal free energy, comprising: obtaining non-bonded parameters by adopting the method for optimizing non-bonded parameters of the force field as mentioned above, or obtaining force field parameters by adopting the method for optimizing force field parameters as mentioned above; and obtaining the free energy of the crystal structure based on the non-bonded parameters and / or force field parameters.

[0008] The fourth aspect of the present application provides a method for predicting a crystal form, comprising: obtaining the free energy of a plurality of crystal structures using the above-mentioned crystal free energy calculation method; and determining a target crystal form according to the free energy of the plurality of crystal structures.

[0009] The fifth aspect of the present application provides a design method, comprising: optimizing force field parameters according to the above method; and performing drug design or material design based on the optimized force field parameters.

[0010] The sixth aspect of the present application provides a device for optimizing the non-bonded parameters of a force field, comprising: a cluster acquisition module, used to obtain M clusters, each cluster in the M clusters includes N molecules, M and N are integers greater than or equal to 1; a non-bonded parameter determination module, used to determine the quantitative intermolecular interaction energy of each polymer of at least some of the M clusters; for each cluster in the M clusters, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then the charge parameters of the force field are obtained based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is a charge determined by the atomic charge of each molecule in the cluster except the central molecule; based on the quantitative intermolecular interaction energy of each polymer of the cluster, the atomic charge of each atom in the central molecule is used as the initial charge of each atom in the molecule in the polymer, the bonding parameters and the initial charge are fixed unchanged, and the van der Waals parameters of the molecule are obtained by fitting.

[0011] The seventh aspect of the present application provides a device for optimizing force field parameters, including: a force field parameter determination module, used to determine bonding parameters, and determine non-bonding parameters based on the above device; a force field parameter optimization module, used to obtain force field parameters based on bonding parameters and non-bonding parameters.

[0012] The eighth aspect of the present application provides a crystal free energy calculation device, including: a force field parameter acquisition module, used to obtain non-bonded parameters using the above-mentioned device for optimizing the non-bonded parameters of the force field, or to obtain force field parameters using the above-mentioned device for optimizing the force field parameters; a crystal structure free energy acquisition module, used to calculate and obtain the free energy of the crystal structure based on the non-bonded parameters and / or force field parameters.

[0013] The ninth aspect of the present application provides a crystal form prediction device, comprising: a multi-free energy calculation module, used to obtain the free energy of multiple crystal structures using the above-mentioned crystal free energy calculation device; a target crystal form determination module, used to determine the target crystal form based on the free energy of multiple crystal structures.

[0014] The tenth aspect of the present application provides a design device, comprising: a parameter optimization module, used to optimize force field parameters according to the above device; and a design module, used to perform drug design or material design based on the optimized force field parameters.

[0015] The eleventh aspect of the present application provides an electronic device, comprising: a processor; and a memory on which executable codes are stored, and when the executable codes are executed by the processor, the processor executes the method.

[0016] The twelfth aspect of the present application also provides a computer-readable storage medium on which executable code is stored. When the executable code is executed by a processor of an electronic device, the processor executes the above method.

[0017] The thirteenth aspect of the present application also provides a computer program product, including executable code, which implements the above method when executed by a processor.

[0018] The method, device, design method and device for optimizing the force field and its non-bonded parameters provided in the present application obtain the atomic charge of each atom in the central molecule by fitting the electrostatic potential of the central molecule of the cluster based on the background charge of the cluster, and fully consider the polarization effect, so as to obtain a non-polarized force field capable of characterizing the polarization effect. Due to the non-bonded effect of the force field of the embodiment of the present application, the polarization effect generated by the molecule in the actual simulation is taken into account, which helps to improve the description accuracy of the interaction between molecules, thereby effectively improving the accuracy of the force field parameters.

[0019] In addition, after optimizing the non-bonding parameters in the force field parameters, considering that the change of non-bonding parameters may have a certain impact on the bonding parameters, after optimizing the non-bonding parameters, the molecular potential energy surface fitting of the non-bonding parameters and the bonding parameters is also performed to optimize the bonding parameters.

[0020] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] The above and other objects, features and advantages of the present application will become more apparent by describing in more detail exemplary embodiments of the present application in conjunction with the accompanying drawings, wherein the same reference numerals generally represent the same components in the exemplary embodiments of the present application.

[0022] Figure 1 An exemplary system architecture of the method, device, design method and device for optimizing force field and non-bonded parameters thereof according to an embodiment of the present application is schematically shown;

[0023] Figure 2 A flowchart of a method for optimizing non-bonded parameters of a force field according to an embodiment of the present application is schematically shown;

[0024] Figure 3 A schematic diagram of a rotatable key according to an embodiment of the present application is schematically shown;

[0025] Figure 4 A logic diagram of a method for determining van der Waals forces based on a global optimal molecular conformation according to an embodiment of the present application is schematically shown;

[0026] Figure 5 A logic diagram of a method for determining van der Waals forces based on experimental molecular conformation according to an embodiment of the present application is schematically shown;

[0027] Figure 6 A flow chart schematically shows a method for optimizing force field parameters according to an embodiment of the present application;

[0028] Figure 7 A logic diagram of a method for determining force field parameters according to an embodiment of the present application is schematically shown;

[0029] Figure 8 A logic diagram schematically illustrates another method for determining force field parameters according to an embodiment of the present application;

[0030] Fig. 9 A flow chart of a method for calculating crystal free energy according to an embodiment of the present application is schematically shown;

[0031] Fig.10 A flow chart of a crystal form prediction method according to an embodiment of the present application is schematically shown;

[0032] Fig.11 A flowchart of a design method according to an embodiment of the present application is schematically shown;

[0033] Fig.12 A block diagram of an apparatus for optimizing non-bonded parameters of a force field according to an embodiment of the present application is schematically shown;

[0034] Fig.13 A block diagram of a device for optimizing force field parameters according to an embodiment of the present application is schematically shown;

[0035] Fig.14 A block diagram of a crystal free energy calculation device according to an embodiment of the present application is schematically shown;

[0036] Fig.15 A block diagram of a crystal form prediction device according to an embodiment of the present application is schematically shown;

[0037] Fig.16 A block diagram of a design device according to an embodiment of the present application is schematically shown;

[0038] Fig.17 A block diagram of an electronic device according to an embodiment of the present application is schematically shown. DETAILED DESCRIPTION

[0039] The embodiments of the present application will be described in more detail below with reference to the accompanying drawings. Although the embodiments of the present application are shown in the accompanying drawings, it should be understood that the present application can be implemented in various forms and should not be limited by the embodiments described herein. On the contrary, these embodiments are provided to make the present application more thorough and complete, and to fully convey the scope of the present application to those skilled in the art.

[0040] The terms used in this application are only for the purpose of describing specific embodiments and are not intended to limit this application. The terms "comprise", "include", etc. used herein indicate the existence of features, steps, operations and / or components, but do not exclude the existence or addition of one or more other features, steps, operations or components.

[0041] All terms (including technical and scientific terms) used herein have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein should be interpreted as having a meaning consistent with the context of this specification, and should not be interpreted in an idealized or overly rigid manner.

[0042] It should be understood that although the terms "first", "second", "third", etc. may be used in this application to describe various information, this information should not be limited to these terms. These terms are only used to distinguish the same type of information from each other. For example, without departing from the scope of this application, the first information may also be referred to as the second information, and similarly, the second information may also be referred to as the first information. Thus, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of this application, the meaning of "multiple" is two or more, unless otherwise clearly and specifically defined.

[0043] Before describing the technical solution of the present application, some technical terms in the field involved in the present application are first explained.

[0044] The molecular force field is a set of parameterized empirical potential function parameters that describe the interactions between atoms in the system. It mainly consists of two parts: bonding parameters that describe the bonding interactions within molecules and non-bonding parameters that describe the interactions between molecules.

[0045] Bonded parameters, including bond length, bond angle, and dihedral angle parameters, respectively describe the stretching energy of chemical bonds, the bending energy of bond angles, and the torsion energy of dihedral angles.

[0046] Nonbonded parameters, including charge parameters and van der Waals parameters (in addition, some functional forms may also include hydrogen bond parameters, etc.), describe the interaction energy between non-bonded atoms or groups.

[0047] Unit Cell is the most basic geometric unit that constitutes a crystal.

[0048] A supercell is a structure that is formed by expanding a unit cell to a certain size according to a certain symmetry. The volume of a supercell is generally large and the structure inside is orderly.

[0049] Box refers to a system of a certain size established according to needs, which is filled with molecules that need to be simulated. The molecular structure inside can be ordered or disordered.

[0050] A coordinate system is a reference for describing the position and posture of an object, so it is also called a reference system or a reference system. For example, a coordinate system can be a coordinate system created during a simulation, such as a Cartesian coordinate system.

[0051] Coordinates are used to represent the absolute position of an object in a specific coordinate system. Mathematically, the essence of coordinates is ordered logarithms.

[0052] Molecular force fields are widely used in dynamics simulation, crystal prediction and other fields. The energy accuracy of molecular force fields directly determines the success or failure of computational simulation, so how to obtain high-precision force field parameters is crucial for simulation.

[0053] The electrostatic interactions of force fields in related technologies, such as AMBER force fields and CHARMM force fields, are often described by point charge models, lacking descriptions of polarization effects, so they are also called non-polarized force fields. The positive and negative charge centers of some molecules coincide, but when they are in the electric field established by adjacent polar molecules, the positive and negative charge centers of non-polar molecules will shift, thus generating induced dipoles. In these cases, the use of non-polarized force fields cannot well describe polarization effects, resulting in the inability of force field parameters to accurately describe their interaction energy.

[0054] Unlike non-polarizable force fields, polarizable force fields take into account the non-bonded effects of molecules in actual application scenarios, so they can describe the interactions between molecules with higher accuracy. For example, polarizable force fields such as Drude and AMOEBA. However, the computational cost of polarizable force fields is high, and the parameter coverage and optimization degree are not as good as those of traditional force fields in related technologies, which limits the development of polarizable force fields to a certain extent.

[0055] In the process of force field parameterization in related technologies, there is not enough experimental data, and most parameters are fitted from gas phase data obtained from quantum mechanics (QM) density functional theory calculations, including electrostatic terms that contribute greatly to intermolecular interactions. If it is applied to energy calculations or dynamic simulations of condensed phases, the lack of consideration of the polarization environment may lead to deviations in the description of electrostatic interactions, resulting in calculation errors or failures. However, the traditional force fields in related technologies have a high coverage of parameter space and strong universality through long-term optimization and iteration. They can also provide a good description of complex systems, and the computational cost is lower than that of polarization force fields, and they are more universal.

[0056] Based on this, the present application proposes a method that can better describe polarization effects based on the form of non-polarizable force fields, so as to improve the accuracy of the description of intermolecular effects.

[0057] The following will be passed Figures 1 to 17 A method, device, design method and device for optimizing a force field and its non-bonded parameters in an embodiment of the present application are described in detail.

[0058] Figure 1 The schematic diagram shows an exemplary system architecture that can be applied to the method, device, design method and device for optimizing the force field and its non-bonded parameters according to the embodiment of the present application. It should be noted that: Figure 1 What is shown is merely an example of a system architecture to which the embodiments of the present application can be applied, in order to help those skilled in the art understand the technical content of the present application, but it does not mean that the embodiments of the present application cannot be used in other devices, systems, environments or scenarios.

[0059] See also Figure 1 According to the embodiment, the system architecture 100 may include terminal devices 101, 102, 103, a network 104 and a server 105. The network 104 is used to provide a medium for communication links between the terminal devices 101, 102, 103 and the server 105. The network 104 may include various connection types, such as wired, wireless communication links or optical fiber cables, etc.

[0060] The user can use the terminal devices 101, 102, 103 to interact with other terminal devices and the server 105 through the network 104 to receive or send information, such as sending an image file, such as an algorithm program for optimizing force field parameters, such as non-bonded parameters, etc. The terminal devices 101, 102, 103 can be installed with various applications, such as drug development applications, material design applications, web browser applications, database applications, search applications, instant messaging tools, email clients, social platform software and other applications.

[0061] The terminal devices 101, 102, 103 include but are not limited to smart desktop computers, tablet computers, laptop computers, and other electronic devices that can support functions such as modeling, analytical calculations, design, and Internet access.

[0062] The server 105 may receive an image including an algorithm program, etc., and may also send force field parameters to the terminal devices 101, 102, 103. For example, the server 105 may be a background management server, a server cluster, etc.

[0063] It should be noted that the number of terminal devices, networks and servers is only illustrative, and any number of terminal devices, networks and clouds may be provided according to implementation requirements.

[0064] Figure 2 The flowchart of the method for optimizing the non-bonded parameters of the force field according to an embodiment of the present application is schematically shown.

[0065] like Figure 2 As shown, this embodiment provides a method for optimizing non-bonded parameters of a force field, and the method includes operations S210 to S230, which are specifically as follows.

[0066] In operation S210 , M clusters are obtained, each of the M clusters includes N molecules, and M and N are integers greater than or equal to 1.

[0067] A cluster is a structure composed of multiple particles. For example, a relatively stable microscopic or submicroscopic aggregate composed of several or even thousands of atoms, molecules or ions through physical or chemical binding forces. The physical and chemical properties of the cluster vary with the number of particles contained. In this embodiment, it is mainly used in a scene in which a structure composed of multiple molecules (which can be different conformations of the same molecule or different conformations of different molecules) is composed of particles.

[0068] For example, the spatial scale of clusters can range from a few angstroms to hundreds of angstroms, and many of their properties are different from those of single atoms and molecules, solids and liquids, and cannot be obtained by simple linear extension or interpolation of the properties of the two. Therefore, clusters can be regarded as a new level of material structure between atoms, molecules and macroscopic solid matter, and a transitional state of various substances from atoms or molecules to bulk matter, or in other words, they represent the initial state of condensed matter.

[0069] In operation S220 , the quantized intermolecular interaction energy of polymers of at least some of the M clusters is determined.

[0070] In this embodiment, when the cluster contains 1 or 2 molecules, the polymer is a dimer; when the cluster contains X (X ≥ 3) molecules, the polymer is a k*X polymer, where k is a positive integer.

[0071] Taking dimers as an example, there is intermolecular interaction energy between the two molecules in the dimer. The distance between the molecules can be used to determine whether the two molecules form a dimer. For example, whether the distance between the geometric centers of the two molecules is within 3.5 angstroms To determine whether a dimer is formed, Equal to 0.1 nanometer (nm).

[0072] The quantitative intermolecular interaction energy of polymers can be calculated using a variety of related techniques.

[0073] For example, the quantitative intermolecular interaction energy μ(r) is determined based on the Sutherland potential function, as shown in formula (1).

[0074]

[0075] Among them, σ represents the molecular diameter, r is the intermolecular distance, and ε is the coefficient.

[0076] For another example, the quantitative intermolecular interaction energy μ(r) is determined based on the Lennard-Jones potential function, as shown in formula (2).

[0077]

[0078] It should be noted that the above-mentioned method for calculating the intermolecular interaction energy is only an example, and a variety of other models can also be used to calculate the quantified intermolecular interaction energy, which is not limited here.

[0079] In operation S230, for each of the M clusters, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then the charge parameters of the force field are obtained based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is a charge determined by the atomic charge of each molecule in the cluster except the central molecule. The central molecule is the molecule at the center of the cluster. For example, for a single component, it is the molecule of the component. For example, for a two-component, it can be randomly selected or specified; the situation for multiple components is similar to that for two components.

[0080] Atomic charge is the most commonly used model to describe the charge distribution of a chemical system. The net charge of each atom is described by a point charge located at the nucleus. In the related art, the atomic charge can be obtained by Gaussian calculation. However, in the prior art, the atomic charge is often calculated under gas phase conditions or without the presence of background charge. This method does not or less considers the polarization effect. In this embodiment, the atomic charge of each atom of the molecule can be obtained by fitting the electrostatic potential (Electrostatic potential, referred to as ESP, also referred to as fitting electrostatic potential) of the molecule. Electrostatic potential is a very important concept in wave function analysis. With the wave function of the molecule, the electrostatic potential of each point in the three-dimensional space is determined, and it can be calculated directly based on the wave function.

[0081] The fitted electrostatic potential charge is obtained by making the atomic charge reproduce the electrostatic potential near and outside the molecular van der Waals surface calculated based on the wave function as well as possible. For example, it is obtained by minimizing the deviation between the electrostatic potential calculated based on the atomic charge and the electrostatic potential calculated based on the wave function at these fitting points by the least squares method, and the sum of the atomic charges can be constrained to be equal to the net charge of the system by the Lagrange multiplier method at the same time. The better the reproducibility of the electrostatic potential in these areas, the better the atomic charges can be used to describe the intermolecular electrostatic interaction through the classical Coulomb formula, and the fitted atomic charges obtained by fitting the electrostatic potential are in principle the atomic charges with the best reproducibility of the electrostatic potential.

[0082] Specifically, based on atomic charges, the premise for accurately describing the electrostatic interaction between molecules is that the electrostatic potential of the van der Waals surface of the molecule and some areas outside it calculated based on atomic charges is close to the electrostatic potential calculated directly based on the wave function. Fitting the electrostatic potential is to satisfy this requirement as much as possible and infer the atomic charge value that can best describe the electrostatic potential distribution on the molecular surface and nearby areas.

[0083] For example, methods for fitting electrostatic potential include but are not limited to the Merz-Kollman (MK) and CHELPG methods. The solution process can be as follows: first determine the fitting point positions near and outside the molecular van der Waals surface based on the geometric structure and the atomic van der Waals radius. Then calculate the electrostatic potential at these positions based on the wave function. Then construct the A matrix and B vector, perform the matrix operation of q=A^(-1)B, and the resulting q vector contains the atomic charges for the fitted electrostatic potential.

[0084] In a specific embodiment, the molecules other than the central molecule are first used as the "background" to fit the electrostatic potential of the central molecule. Then, the charge (calculate background charge) command of the GAUSSIAN software can be combined with the MK method to fit the electrostatic potential of the central molecule. Next, multiple rounds of iterations are performed until convergence, such as the change in the atomic charge of the central molecule <0.001). Then, these central molecules are subjected to RESP fitting to obtain the final charge parameters, such as averaging the corresponding charges of the central molecules of multiple clusters as the final charge parameters. Among them, this charge parameter is a set, which includes several atomic charges according to the number of atoms in the molecule. In this way, the final charge parameters can be used to replace the charges of each atom in the initial parameter file.

[0085] It should be noted that the initial charge of the atom can be am1bcc. In each round of iteration, the charge value corresponding to the "background" molecule can be replaced by the atomic charge of the fitted central molecule (obtained by electrostatic potential fitting).

[0086] In operation S240, based on the quantized intermolecular interaction energy of each polymer of the cluster, the atomic charge of each atom in the central molecule is used as the initial charge of each atom in the molecule in the polymer, the bonding parameters and the initial charge are fixed, and the van der Waals parameters of the molecule are obtained by fitting.

[0087] In this embodiment, the van der Waals parameters and charge parameters can be used as non-bonded parameters. The van der Waals parameters can be determined by fitting, for example, the van der Waals parameters include but are not limited to: ε, σ in formula (2).

[0088] For example, molecular dimers are extracted from the cluster (for example, the distance between the geometric centers of the two molecules is within 3.5 angstroms, and can be within 4.0 angstroms when the molecules are larger), and the quantitative intermolecular interaction energy of all dimers is calculated using a quantitative method. The force field parameters after replacing the atomic charges are used as the initial parameters, and the bonding parameters (i.e., parameters that express bonding information, such as bond parameters, angle parameters, and dihedral angle parameters) and atomic charges are fixed unchanged, and the van der Waals parameters of the molecules are fitted.

[0089] A dihedral angle is the angle between two intersecting planes. In solid geometry, it is defined as the union of a line and two half-planes, the line being the common edge of the two half-planes. In higher dimensions, a dihedral angle is the angle between two hyperplanes. In chemistry, a dihedral angle is the angle between two planes in a molecule, each consisting of three atoms, for a total of four atoms, where the common edge can be a chemical bond (two atoms) and the planes are made up of two other atoms and the chemical bond. Every polyhedron has a dihedral angle on each edge, describing the relationship between the two faces that share that edge. This dihedral angle is also called the face angle, and it is calculated as the internal angle of the polyhedron. A dihedral angle of 0° means that the two face normals are antiparallel and the two faces overlap, which means it is part of a degenerate polyhedron. A dihedral angle of 180° means that the faces are parallel, as in a tessellation.

[0090] In this embodiment, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule. The influence of polarization effect and other factors on the intermolecular interaction is taken into account, which effectively improves the accuracy of the calculation of intermolecular interactions in the condensed phase and the accuracy of the force field parameters, thereby improving the accuracy of the simulation results using the optimized force field parameters.

[0091] The following is an exemplary description of the process of obtaining M clusters. The method of obtaining M clusters in different scenarios may be different. For example, the scenario may be divided into a scenario where there is no experimental structure (non-experimental structure) and a scenario where there is an experimental structure (experimental structure). In addition, more scenarios may be divided according to actual needs.

[0092] In some embodiments, taking a non-experimental structure scenario as an example, obtaining M clusters may include the following operations.

[0093] First, obtain the three-dimensional (3D) structural information of the molecule. The 3D structural information includes, but is not limited to, the coordinate information of each atom in the molecule in a specified coordinate system. Specifically, the coordinate information can be the position information of the atom in a specific coordinate system represented by a string. For example, the string can include the three-dimensional conformation of the molecule and the x, y, z coordinates of each atom in the molecule.

[0094] It should be noted that the same coordinate system, in particular the spatial coordinate system, can be used in the process of optimizing the force field parameters. For example, the geometric center coordinates of the central molecule and part or all of the atomic coordinates are coordinates in the same coordinate system. It should be understood that the geometric center coordinates of the central molecule and part or all of the atomic coordinates can also be coordinates in different coordinate systems, but the coordinates in each coordinate system can be converted to each other.

[0095] Then, based on the molecular three-dimensional structure information, a conformation search is performed to obtain a specified number of low potential energy molecular conformations. The low potential energy molecular conformation refers to a molecular conformation with a relatively low potential energy or the lowest potential energy. The specified number of types can be 1, 2, 3 or more types.

[0096] In organic compound molecules, the countless specific images of atoms or groups arranged in space produced by the rotation of C-C single bonds are called molecular conformations. It refers to the spatial arrangement produced by changing the placement of atoms around single bonds in a molecule without changing the covalent bond structure. Different molecular conformations can be transformed into each other. Among various molecular conformation forms, the molecular conformation with the lowest potential energy and the most stable can be regarded as the global optimal molecular conformation. When one molecular conformation changes to another, the breaking and reformation of covalent bonds is not required. Among them, the specified number can be based on expert experience or calculation results, and is not limited here. For example, the specified number can be 1, 3, 5, 10, 15, 18, 20 or 25, etc.

[0097] Next, a molecular dynamics (MD) simulation is performed on the candidate box to obtain a stable box, which is formed by stacking multiple low-potential energy molecular conformations. For example, a candidate box with a size of 30×30×30 (in angstroms) is stacked by multiple low-potential energy molecular conformations without colliding with each other. These low-potential energy molecular conformations can be the same or different.

[0098] Then, M clusters are extracted from the stable box. In some embodiments, each of the M clusters is wrapped by at least a specific number of layers of clusters. This can avoid extracting surface clusters in the stable box, which helps to improve the accuracy of the force field parameters finally obtained. The specific number of layers can be at least two layers, such as 2 layers, 3 layers, 4 layers, 5 layers or more layers.

[0099] Through the above method, M clusters that meet the simulation requirements can be easily extracted.

[0100] In some embodiments, for non-experimental structure scenarios, the target molecular conformation may be determined in the following manner: Specifically, molecular conformation search based on molecular three-dimensional structure information to obtain a specified number of low potential energy molecular conformations may include the following operations.

[0101] First, a rotatable bond and a torsion angle corresponding to the rotatable bond are determined based on bond characteristics, wherein the bond characteristics are determined based on the chemical environment of the atoms.

[0102] Then, candidate molecular conformations are determined based on the rotatable bonds and the torsion angles corresponding to the rotatable bonds.

[0103] Next, a specified number of low-potential energy molecular conformations are determined from the candidate molecular conformations, wherein the similarity between each of the specified number of low-potential energy molecular conformations is less than a similarity threshold. The similarity threshold can be determined based on expert experience or by imitation simulation. This scheme can improve the representativeness of the determined low-potential energy molecular conformations and remove redundancy.

[0104] Figure 3 The figure schematically shows a schematic diagram of a rotatable key according to an embodiment of the present application.

[0105] like Figure 3 In a protein chain, three dihedral angles are defined as φ, ψ, and ω, as shown in the figure. The planarity of the peptide bond usually constrains ω to 180° (trans case) or 0° (cis case). The distances between the α-carbons of the trans and cis isomers are approximately 3.8 and 1.5°, respectively. Figure 3 The bonds corresponding to the three dihedral angles are rotatable bonds. Cα and Cβ are one of the six atoms that make up the peptide bond plane, and (+1) and (-1) represent different objects. A peptide bond is an amide bond formed by the dehydration condensation of the α-carboxyl group of one molecule of amino acid and the α-amino group of another molecule of amino acid, i.e. -CO-NH-. Amino acids are linked to form polypeptide chains by peptide bonds. Peptide bonds are the main covalent bonds in protein molecules and are relatively stable. Peptide bonds have the properties of partial double bonds, which are difficult to rotate freely and have a certain rigidity, thus forming a peptide bond plane. The peptide bond plane includes the spatial positions of the six atoms, C═O, N-H and two Cα, that connect the two ends of the peptide bond, which are in a relatively close plane, and the side chains R of the two adjacent amino acids form a trans configuration, thus forming a complex spatial structure of the peptide bond and the peptide chain.

[0106] Figure 4 The logical diagram of the method for determining van der Waals forces based on the global optimal molecular conformation according to an embodiment of the present application is schematically shown.

[0107] See also Figure 4 , first perform a molecular conformation search based on the three-dimensional structure information of the molecule. When the number of molecular conformations obtained is less than the specified number (such as 100+FDN×FDN×10, FDN is the number of flexible angles), you can use EDKTG of rdkit to continue to supplement the molecular conformations. Then use DFTB to perform structural optimization on the molecular conformation set. Next, remove duplicate structures through structural clustering, and perform energy ranking screening on the remaining structures to select the globally optimal molecular conformation. Then, grab the gaff2 parameter file from the globally optimal molecular conformation as a template or initial parameter for subsequent fitting. Among them, the globally optimal molecular conformation can be the molecular conformation with the lowest potential energy among multiple low-potential energy molecular conformations. The globally optimal molecular conformation can also be the weighted average of multiple low-potential energy molecular conformations.

[0108] Among them, when searching for molecular conformations, first determine the characteristics of the bonds according to the chemical environment in which the atoms are located, identify at least some of the rotatable bonds, and determine which torsion angles each bond has to consider according to the built-in library. Assuming there are four bonds, there are 6, 3, 12, and 4 torsion angles to be considered, then the total number in the molecular conformation space is 6×3×12×4=864. You can pre-set how many molecular conformations to take in the end, for example, if you set 100, then 100 random structures will be obtained by randomly selecting torsion angles. After that, you can filter out those with too high energy, and remove molecular conformations that are too similar by the RMSD of heavy atoms between molecular conformations, so as to ensure the diversity and relatively balanced weight of the generated molecular conformations. For example, the search algorithm can be the confab module of openbabel and / or ETKDG (Experimental-Torsion BasicKnowledge Distance Geometry, referred to as ETKDG).

[0109] In some embodiments, if the molecular conformation is relatively complex, the molecular conformation with the lowest potential energy that may be selected is not the global optimal molecular conformation, and a weighted average can be taken from multiple low-potential molecular conformations to reduce the error caused by a single molecular conformation. Specifically, after determining a specified number of low-potential molecular conformations from the candidate molecular conformations, the above method can also include the following operations: performing weighted averaging on the parameter values ​​of the specified number of low-potential molecular conformations to obtain the global optimal molecular conformation. Accordingly, the candidate box is formed by stacking multiple global optimal molecular conformations. At this time, there is only one molecular conformation in the candidate box.

[0110] For example, based on the multiple low potential energy molecular conformations determined in the above embodiment, a global optimal molecular conformation (such as a virtual molecular conformation) is obtained by weighted average processing. The weights of each low potential energy molecular conformation may be the same or different. For example, among multiple low potential energy molecular conformations, the greater the difference between a certain low potential energy molecular conformation and other low potential energy molecular conformations, the greater or smaller the weight of the certain low potential energy molecular conformation.

[0111] In some embodiments, for a scenario with experimental molecular conformations, M clusters may be determined in the following manner: Specifically, obtaining M clusters may include the following operations.

[0112] First, a supercell of the experimental crystal structure is obtained. For example, the size of the supercell of the experimental crystal structure is 30×30×30 (in angstroms), which includes multiple molecules with different conformations.

[0113] Then, the supercell of the experimental crystal structure is simulated by molecular dynamics to obtain a stable supercell. For example, a 5ns NVT dynamics simulation is performed on the supercell structure at 300K, and the final stable supercell structure is saved. Among them, the rmsd change of the structure in the simulation trajectory file tends to be stable, which means that the stable state is reached.

[0114] In addition, M clusters can be extracted from the stable supercell. For example, each of the M clusters is surrounded by at least two layers of clusters. For example, several clusters with a radius of 10 angstroms are extracted from the supercell.

[0115] Figure 5 The logical diagram of the method for determining van der Waals forces based on experimental molecular conformation according to an embodiment of the present application is schematically shown.

[0116] See also Figure 5 ,and Figure 4 The difference is that a supercell of approximately 30×30×30 (in angstroms) is directly obtained based on the experimental crystal structure, and a 5ns MD simulation is performed without determining the global optimal molecular conformation. Figure 4 The same contents will not be described in detail here.

[0117] The above method can be used to optimize the force field parameters of experimental molecules with specific molecular conformations.

[0118] In some embodiments, for non-experimental structural scenarios, the above method may further include the following operations: for a multi-component system, the molecules of each component system are stacked according to the component ratio of the multi-component system to obtain a box.

[0119] Accordingly, obtaining the M clusters includes: sampling the box according to the components of the multi-component system to obtain the M clusters.

[0120] Specifically, if it is a multi-component system, packing (i.e., the process of forming a molecular structure stack) is performed in the process of determining the non-bonded parameters, and each component is stacked into a box at the same time according to the required ratio (similar to a supercell, the structure in the box in this application can be disordered), and MD sampling is performed. The extraction of clusters must maintain the ratio of the components, and the extraction of polymers must be based on the ratio of the components. Among them, the polymer can be a dimer or a trimer, etc.

[0121] In some embodiments, for the experimental structure scenario, obtaining M clusters in the above method includes: performing dynamic sampling on the supercell to obtain M clusters.

[0122] In the embodiments of the present application, the background charge of the cluster is used to fit the electrostatic potential of the central molecule of the cluster to obtain the atomic charge of each atom in the central molecule, which effectively improves the accuracy of the calculation of intermolecular interactions in the condensed phase and improves the success rate of the simulation.

[0123] It should be noted that when the cluster contains only one type of molecule, the central molecule is a molecule of a certain conformation in the cluster; when the cluster contains two or more types of molecules, the central molecule can randomly select a molecule of a certain conformation of one of the two or more molecules as needed, and based on the different molecular types selected by the central molecule, the non-bonded parameters obtained by subsequent fitting are the non-bonded parameters corresponding to that molecular type.

[0124] Another aspect of the present application also provides a method for optimizing force field parameters.

[0125] Figure 6 The flowchart of the method for optimizing force field parameters according to an embodiment of the present application is schematically shown.

[0126] like Figure 6 As shown, the method for optimizing force field parameters may include operations S610 to S620.

[0127] In operation S610, bonding parameters are determined, and non-bonding parameters are determined based on the above method. Specifically, the non-bonding parameters may include charge parameters and van der Waals parameters. The specific calculation process is as shown above and will not be described in detail here.

[0128] In operation S620, force field parameters are obtained based on the bonding parameters and the non-bonding parameters.

[0129] In this embodiment, the bonding parameters include at least one of a bond parameter, an angle parameter, and a dihedral angle parameter. The bonding parameters are physical quantities that can characterize the properties of chemical bonds and are called bond parameters. For example, the bond parameters mainly include force constants, bond lengths, etc. The angle parameters may include: force constants, angles between bonds, etc. The dihedral angle and its parameters may be as shown in the above related content and will not be described in detail here.

[0130] It should be noted that a variety of related technologies can be used to determine the bonding parameters, which are not limited here.

[0131] The field parameters determined in this embodiment take into account the influence of the polarization effect on the force field parameters, so that the force field parameters can better describe the force field, which helps to improve the accuracy of the results obtained by simulation based on the force field parameters.

[0132] In some embodiments, in order to further improve the accuracy of the force field parameters, after obtaining the optimized non-bonded parameters, the bonding parameters may be optimized based on the optimized non-bonded parameters. Specifically, after obtaining the non-bonded parameters, the above method may further include the following operations.

[0133] First, the quantized energy of each scanning point is calculated to obtain potential energy surface data. The scanning point is a point scanned during the potential energy surface scanning of at least a part of the flexible angle of the target molecular conformation.

[0134] Then, the potential energy surface data are used as training data to fit the non-bonding parameters and bonding parameters to obtain the fitted dihedral angle parameters.

[0135] Then, the bond parameters, angle parameters and fitted dihedral parameters are taken as the optimized bonding parameters.

[0136] Among them, the dihedral angle parameter has a great influence on the molecular conformation of a single molecule, and the dihedral angle assigned by the force field through the atomic type may not be accurate. Therefore, when simulating flexible molecules, it is usually necessary to optimize the dihedral angle to obtain a more accurate molecular conformation.

[0137] For example, a one-dimensional potential energy surface scan is performed on each flexible angle of the global optimal molecular conformation or experimental molecular conformation obtained in the above embodiment, and the calculation results of each scanning point under a given quantization method are calculated. These structures and their energies are used as a fitting set of force field dihedral angle parameters. Wherein, for example, the quantization energy can be determined by quantization software, such as using GAUSSIAN09. The quantization method can use b3lyp / 6-31g (d).

[0138] Then, the scanned structure and QM energy can be used as a training set to fit the molecular potential energy surface to the force field parameters determined in the above embodiment to obtain fitting dihedral angle parameters. The fitting dihedral angle parameters can more realistically describe the dihedral angle potential energy surface in the condensed phase system.

[0139] In certain embodiments, each of the M clusters is formed by stacking target molecular conformations, and the target molecular conformations include: experimental molecular conformations and / or global optimal molecular conformations. Among them, the experimental molecular conformations can be directly input by the user or called from a database. The global optimal molecular conformation can be the molecular conformation with the lowest potential energy among multiple possible molecular conformations. The determination method can refer to the relevant content in the above embodiment, which will not be described in detail here.

[0140] The following is an exemplary description of a method for determining bonding parameters.

[0141] In some embodiments, the bonding parameters and angle parameters in the bonding parameters may be determined as follows.

[0142] First, the Hessian matrix of the experimental molecular conformation or the global optimal molecular conformation is determined. For example, the Hessian matrix of the global optimal molecular conformation or the experimental molecular conformation is calculated using a quantitative method.

[0143] Among them, the Hessian matrix, also known as the Hessian matrix, the Hessian matrix, the Hessian matrix, etc., is a square matrix composed of the second-order partial derivatives of a multivariate function, which describes the local curvature of the function. The Hessian matrix is ​​often used to solve optimization problems using the Newton method. The Hessian matrix can be used to determine the extreme value problem of a multivariate function. In the optimization scenario of force field parameters, the listed objective functions are often very complex. In order to simplify the problem, the objective function is often expanded into a Taylor polynomial in the neighborhood of a certain point to approximate the original function. At this time, the Hessian matrix will be involved in the matrix form of the Taylor expansion of the function at a certain point.

[0144] Then, the Hessian matrix is ​​converted into the force constant of each bond and the force constant of each angle of the target molecular conformation; and the molecular bond length and molecular bond angle of the target molecular conformation are determined.

[0145] For example, the Seminario method can be used to convert the Hessian matrix into force constants for each bond and each angle of the molecule, and the molecular bond length and bond angle optimized by the quantitative method can be used as the equilibrium value of the bond parameter and the angle parameter. It should be noted that when there are multiple molecular conformations, the average value of the bonds and angles of multiple molecular conformations can be selected as the equilibrium value of the force field parameter. This method can reduce errors when the molecular conformation is more complex.

[0146] Next, the force constant of each bond of the target molecular conformation and the corresponding molecular bond length are taken as bond parameters, and the force constant of each angle of the target molecular conformation and the corresponding molecular bond angle are taken as angle parameters.

[0147] For example, the initial parameter value is replaced with the optimized value, wherein the bond parameter includes the bond force constant and the bond equilibrium value, and the angle parameter includes the angle force constant and the angle equilibrium value.

[0148] Figure 7 A logic diagram of a method for determining force field parameters according to an embodiment of the present application is schematically shown.

[0149] Please also read Figure 5 and Figure 7 As shown, Figure 5 The difference is, Figure 7 After obtaining the experimental molecular 3D conformation, the Hessian matrix corresponding to the experimental molecular 3D conformation is also calculated, and then the Hessian matrix is ​​converted into bond parameters and angle parameters by the Seminario method. In addition, the QM dihedral angle scanning point training set is obtained by flexible angle scanning to fit the dihedral angle parameters. In this way, a complete force field parameter file can be obtained based on non-bonded parameters and bonding parameters.

[0150] Figure 8 A logic diagram of another method for determining force field parameters according to an embodiment of the present application is schematically shown.

[0151] Please also read Figure 4 and Figure 8 As shown, Figure 4 The difference is, Figure 8 After obtaining the global optimal molecular conformation (such as the global minimum potential energy), the Hessian matrix corresponding to the global optimal molecular conformation is also calculated, and then the Hessian matrix is ​​converted into bond parameters and angle parameters by the Seminario method. In addition, a QM dihedral angle scanning point training set is obtained by flexible angle scanning to fit the dihedral angle. In this way, a complete force field parameter file can be obtained based on non-bonded parameters and bonding parameters.

[0152] In this embodiment, a process for quickly determining the parameterization of a "proprietary force field" is provided. After determining the global optimal molecular conformation, the global optimal molecular conformation is selected, and the gas-phase RESP fitting of the atomic charges is performed while keeping the initial van der Waals parameters unchanged. Finally, the dihedral angle parameters are fitted to obtain a set of optimized force field parameters.

[0153] Molecules in the condensed phase are affected by the surrounding chemical environment and will produce a polarization effect. Therefore, in order to more accurately describe the charge distribution of the molecules, this embodiment takes into account the chemical environment around the molecules in the condensed phase when fitting the charge, and the description of the electrostatic interaction between molecules will be more accurate. Specifically, the electrostatic potential of organic molecules is fitted in the presence of background charge to obtain the charge parameters of the force field, which effectively improves the accuracy of the force field parameters.

[0154] In addition, the parameters of the force field items of the general force field are interdependent. If the charge parameters are directly replaced, the overall energy or results calculated or simulated will be biased, especially the electrostatic parameters and van der Waals parameters. In this embodiment, after determining the atomic charge, other parameters are also regenerated or fitted. This embodiment effectively improves the accuracy of the force field parameters.

[0155] Another aspect of the present application provides a method for calculating crystal free energy.

[0156] Fig. 9 The figure schematically shows a flow chart of a method for calculating crystal free energy according to an embodiment of the present application.

[0157] See also Fig. 9 , the crystal free energy calculation method may include the following operations.

[0158] In operation S910, the non-bonded parameters are obtained by using the above method for optimizing the non-bonded parameters of the force field, or the force field parameters are obtained by using the above method for optimizing the force field parameters.

[0159] In operation S920 , free energy of the crystal structure is calculated based on non-bonded parameters and / or force field parameters.

[0160] Among them, the free energy prediction of crystal structure aims to theoretically evaluate the relative stability of different crystal forms of materials (such as drugs) with temperature changes, evaluate experimental crystal forms, predict crystal forms, and the relative stability changes of hydrates, co-crystals, free forms, etc., so as to quickly select the dominant crystal form of the target drug.

[0161] Another aspect of the present application provides a method for calculating crystal free energy.

[0162] Fig.10 A flow chart of a crystal form prediction method according to an embodiment of the present application is schematically shown.

[0163] See also Fig.10 The above-mentioned crystal form prediction method may include operations S1010 to S1020.

[0164] In operation S1010, the free energies of various crystal structures are obtained using the above crystal free energy calculation method.

[0165] In operation S1020, a target crystal form is determined according to free energy of the plurality of crystal structures.

[0166] Crystal structure research is one of the key points of compound research. Crystal structure can strongly affect the bioavailability of drugs and ultimately affect the therapeutic performance of drugs. The crystal form prediction method can realize the crystal structure analysis of drug molecules through computer simulation, which can not only help researchers simplify experimental steps, but also screen out a few potential crystal structures with clear targets for experimental crystallization.

[0167] Another aspect of the present application provides a design method.

[0168] Fig.11 The flowchart of the design method according to the embodiment of the present application is schematically shown.

[0169] See also Fig.11 , the design method may include operation S1110 and operation S1120.

[0170] In operation S1110 , the force field parameters are optimized according to the method described above.

[0171] The method for optimizing the force field parameters can be as shown above and will not be described in detail here.

[0172] In operation S1120 , drug design or material design is performed based on the optimized force field parameters.

[0173] The force field parameters using background charge can be applied in the calculation of condensed phase thermodynamic properties, such as energy calculation. Drug design or material design can then be performed based on condensed phase thermodynamic properties, thus improving the design success rate and development efficiency.

[0174] Another aspect of the present application also provides a device for optimizing non-bonded parameters of a force field.

[0175] Fig.12 A block diagram of an apparatus for optimizing non-bonded parameters of a force field according to an embodiment of the present application is schematically shown.

[0176] See also Fig.12 The non-bonded parameter optimization device 1200 for the force field may include: a cluster obtaining module 1210 and a non-bonded parameter determining module 1220 .

[0177] The cluster obtaining module 1210 is used to obtain M clusters, each of the M clusters includes N molecules, and M and N are integers greater than or equal to 1.

[0178] The non-bonded parameter determination module 1220 is used to determine the quantitative intermolecular interaction energy of each polymer of at least some of the M clusters; for each cluster in the M clusters, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then the charge parameters of the force field are obtained based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is the charge determined by the atomic charge of all molecules in the cluster except the central molecule; based on the quantitative intermolecular interaction energy of each polymer of the cluster, the atomic charge of each atom in the central molecule is used as the initial charge of each atom in the molecule in the polymer, the bonding parameters and the initial charge are fixed unchanged, and the van der Waals parameters of the molecule are obtained by fitting.

[0179] In some embodiments, the cluster acquisition module 1210 includes: a three-dimensional structure information acquisition unit, a conformation search unit, a first molecular dynamics simulation unit, and a first cluster extraction unit.

[0180] Wherein, the three-dimensional structure information obtaining unit is used to obtain the three-dimensional structure information of the molecule.

[0181] The conformation search unit is used to search for molecular conformations based on the three-dimensional structure information of the molecule, and obtain a specified number of types of low-potential energy molecular conformations.

[0182] The first molecular dynamics simulation unit is used to perform molecular dynamics simulation on a candidate unit cell to obtain a stable unit cell. The candidate unit cell is formed by stacking a plurality of low potential energy molecular conformations, and each low potential energy molecular conformation may be the same or different.

[0183] The first cluster extraction unit is used to extract M clusters from the stable unit cell.

[0184] In certain embodiments, the conformation search unit includes: a rotatable bond and angle determination subunit, a candidate molecular conformation determination subunit, and a low potential energy molecular conformation determination subunit.

[0185] The rotatable bond and angle determination subunit is used to determine the rotatable bond and the torsion angle corresponding to the rotatable bond based on the characteristics of the bond, wherein the characteristics of the bond are determined based on the chemical environment of the atom.

[0186] The candidate molecular conformation determining subunit is used to determine the candidate molecular conformation based on the rotatable bonds and the torsion angles corresponding to the rotatable bonds.

[0187] The low potential energy molecular conformation determination subunit is used to determine a specified number of low potential energy molecular conformations from candidate molecular conformations, wherein the similarity between each of the specified number of low potential energy molecular conformations is less than a similarity threshold.

[0188] In some embodiments, the above-mentioned device 1200 may also include: a global optimal molecular conformation determination module, including a first optimal conformation obtaining unit and a second optimal conformation obtaining unit.

[0189] The first optimal conformation obtaining unit is used to select the molecular conformation with the lowest potential energy from the candidate molecular conformations as the global optimal molecular conformation when the designated number of types is 1.

[0190] The second optimal conformation obtaining unit is used to perform weighted average processing on the parameter values ​​of the low potential energy molecular conformations of the specified number of types when the specified number of types is multiple, so as to obtain the global optimal molecular conformation.

[0191] In certain embodiments, the candidate box is stacked from multiple global optimal molecular conformations.

[0192] In some embodiments, the cluster acquisition module 1210 includes: an experimental supercell acquisition unit and a second molecular dynamics simulation unit.

[0193] The experimental supercell acquisition unit is used to obtain supercells of experimental crystal structures.

[0194] The second molecular dynamics simulation unit is used to perform molecular dynamics simulation on the supercell of the experimental crystal structure to obtain a stable supercell structure, so as to obtain M clusters based on the stable supercell structure.

[0195] In some embodiments, the cluster acquisition module 1210 is specifically configured to extract M clusters from a stable unit cell, and each of the M clusters is respectively wrapped by at least two layers of clusters.

[0196] In some embodiments, the above-mentioned device further includes: a first molecular stacking module, which is used to stack the molecules of each component system according to the component ratio of the multi-component system to obtain a box.

[0197] Accordingly, the cluster obtaining module 1210 is specifically used to sample the box according to the components of the multi-component system to obtain M clusters.

[0198] In some embodiments, the cluster acquisition module 1210 is specifically used to perform dynamic sampling on the supercell to obtain M clusters.

[0199] In certain embodiments, the M clusters form at least two layers to wrap the central molecule.

[0200] Another aspect of the present application also provides a method for optimizing force field parameters.

[0201] Fig.13 A block diagram of an apparatus for optimizing force field parameters according to an embodiment of the present application is schematically shown.

[0202] See also Fig.13 The above-mentioned device 1300 for optimizing force field parameters includes: a force field parameter determination module 1310 and a force field parameter optimization module 1320.

[0203] The force field parameter determination module 1310 is used to determine bonding parameters, and to determine non-bonding parameters based on the above method.

[0204] The force field parameter optimization module 1320 is used to obtain force field parameters based on bonding parameters and non-bonding parameters.

[0205] In certain embodiments, the bonding parameter includes at least one of a bonding parameter, an angle parameter, and a dihedral parameter.

[0206] The above-mentioned device 1300 may further include: a potential energy surface data acquisition module, a molecular potential energy surface fitting module and a bonding parameter acquisition module.

[0207] The potential energy surface data acquisition module is used to calculate the quantized energy of each scanning point after obtaining the non-bonded parameters to obtain potential energy surface data. The scanning point is the point scanned during the potential energy surface scanning of at least part of the flexible angle of the target molecular conformation.

[0208] The molecular potential energy surface fitting module is used to use the potential energy surface data as training data, perform molecular potential energy surface fitting on non-bonding parameters and bonding parameters, and obtain fitting dihedral angle parameters.

[0209] The bonding parameter acquisition module is used to use the bonding parameters, angle parameters and fitted dihedral angle parameters as optimized bonding parameters.

[0210] In certain embodiments, each of the M clusters is formed by stacking target molecular conformations, and the target molecular conformations include: experimental molecular conformations and / or global optimal molecular conformations.

[0211] In some embodiments, the force field parameter determination module 1310 includes: a Hessian matrix determination unit, a force constant determination unit, and a bond angle parameter determination unit.

[0212] The Hessian matrix determination unit is used to determine the Hessian matrix of the experimental molecular conformation or the global optimal molecular conformation.

[0213] The force constant determination unit is used to convert the Hessian matrix into the force constants of each bond and each angle of the target molecular conformation; and determine the molecular bond length and molecular bond angle of the target molecular conformation.

[0214] The bond angle parameter determination unit is used to use the force constant of each bond of the target molecular conformation and the corresponding molecular bond length as bond parameters, and use the force constant of each angle of the target molecular conformation and the corresponding molecular bond angle as angle parameters.

[0215] Another aspect of the present application provides a crystal free energy calculation device.

[0216] Fig.14 The block diagram of the crystal free energy calculation device according to an embodiment of the present application is schematically shown.

[0217] See also Fig.14 The crystal free energy calculation device 1400 may include a force field parameter acquisition module 1410 and a crystal structure free energy acquisition module 1420 .

[0218] The force field parameter obtaining module 1410 is used to obtain non-bonded parameters by using the above-mentioned device 1200 for optimizing non-bonded parameters of the force field, or to obtain force field parameters by using the above-mentioned device 1300 for optimizing force field parameters.

[0219] The crystal structure free energy obtaining module 1420 is used to calculate and obtain the free energy of the crystal structure based on non-bonded parameters and / or force field parameters.

[0220] Another aspect of the present application also provides a crystal form prediction device.

[0221] Fig.15 A block diagram of a crystal form prediction device according to an embodiment of the present application is schematically shown.

[0222] See also Fig.15 The above-mentioned device 1500 includes: a multi-free energy calculation module 1510 and a target crystal form determination module 1520.

[0223] The multi-free energy calculation module 1510 is used to obtain the free energy of multiple crystal structures using the crystal free energy calculation device 1400 .

[0224] The target crystal form determination module 1520 is used to determine the target crystal form according to the free energy of multiple crystal structures.

[0225] Another aspect of the present application also provides a design device.

[0226] Fig.16 A block diagram of a design device according to an embodiment of the present application is schematically shown.

[0227] See also Fig.16 , the device 1600 may include: a parameter optimization module 1610 and a design module 1600 .

[0228] The parameter optimization module 1610 is used to optimize the force field parameters according to the above method.

[0229] The design module 1620 is used to perform drug design or material design based on the optimized force field parameters.

[0230] Regarding the devices 1200, 1300, 1400, 1500, and 1600 in the above embodiments, the specific manner in which each module and unit performs operations has been described in detail in the embodiments of the method, and will not be elaborated again here.

[0231] Another aspect of the present application provides an electronic device.

[0232] Fig.17 A block diagram of an electronic device according to an embodiment of the present application is schematically shown.

[0233] See also Fig.17 , the electronic device 1700 includes a memory 1710 and a processor 1720.

[0234] The processor 1720 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or the processor may be any conventional processor, etc.

[0235] The memory 1710 may include various types of storage units, such as system memory, read-only memory (ROM), and permanent storage. Among them, ROM may store static data or instructions required by the processor 1720 or other modules of the computer. The permanent storage may be a readable and writable storage device. The permanent storage may be a non-volatile storage device that does not lose the stored instructions and data even after the computer is powered off. In some embodiments, the permanent storage device uses a large-capacity storage device (such as a magnetic or optical disk, flash memory) as a permanent storage device. In some other embodiments, the permanent storage device may be a removable storage device (such as a floppy disk, optical drive). The system memory may be a readable and writable storage device or a volatile readable and writable storage device, such as a dynamic random access memory. The system memory may store some or all instructions and data required by the processor at run time. In addition, the memory 1710 may include any combination of computer-readable storage media, including various types of semiconductor memory chips (such as DRAM, SRAM, SDRAM, flash memory, programmable read-only memory), and disks and / or optical disks may also be used. In some embodiments, the memory 1710 may include a readable and / or writable removable storage device, such as a laser disc (CD), a read-only digital versatile disc (such as a DVD-ROM, a double-layer DVD-ROM), a read-only Blu-ray disc, an ultra-density optical disc, a flash memory card (such as an SD card, a mini SD card, a Micro-SD card, etc.), a magnetic floppy disk, etc. The computer-readable storage medium does not include carrier waves and transient electronic signals transmitted wirelessly or wired.

[0236] The memory 1710 stores executable codes, and when the executable codes are processed by the processor 1720 , the processor 1720 may execute part or all of the methods described above.

[0237] In addition, the method according to the present application may also be implemented as a computer program or a computer program product, which includes computer program code instructions for executing some or all of the steps in the above method of the present application.

[0238] Alternatively, the present application can also be implemented as a computer-readable storage medium (or non-transitory machine-readable storage medium or machine-readable storage medium) on which executable code (or computer program or computer instruction code) is stored. When the executable code (or computer program or computer instruction code) is executed by a processor of an electronic device (or server, etc.), the processor executes part or all of the steps of the above-mentioned method according to the present application.

[0239] The embodiments of the present application have been described above, and the above description is exemplary, not exhaustive, and is not limited to the disclosed embodiments. Many modifications and changes will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The selection of terms used herein is intended to best explain the principles of the embodiments, practical applications, or improvements to the technology in the market, or to enable other persons of ordinary skill in the art to understand the embodiments disclosed herein.

Claims

1. A method for optimizing non-bonded parameters of a force field, characterized in that: include: Obtaining M clusters, each of the M clusters includes N molecules, where M and N are integers greater than or equal to 1; Determining the quantified intermolecular interaction energy of each polymer of at least some of the M clusters; For each cluster in the M clusters, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then the charge parameters of the force field are obtained based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is a charge determined by the atomic charge of each molecule in the cluster except the central molecule, and the central molecule is a molecule at the center of the cluster; Based on the quantified intermolecular interaction energy of each polymer of the cluster, the atomic charge of each atom in the central molecule is used as the initial charge of each atom in the molecule in the polymer, the bonding parameters and the initial charge are fixed, and the van der Waals parameters of the molecule are obtained by fitting.

2. The method according to claim 1, characterized in that: The obtaining of M clusters comprises: Obtain molecular three-dimensional structure information; Performing conformation search based on the three-dimensional structure information of the molecule to obtain a specified number of low potential energy molecular conformations; Performing molecular dynamics simulation on the candidate box to obtain a stable box, wherein the candidate box is formed by stacking a plurality of low potential energy molecular conformations of the specified number of types; The M clusters are extracted from the stable box.

3. The method according to claim 2, characterized in that The conformation search based on the three-dimensional structure information of the molecule is performed to obtain a specified number of types of low potential energy molecular conformations, including: Determining a rotatable bond and a torsion angle corresponding to the rotatable bond based on a characteristic of the bond, wherein the characteristic of the bond is determined based on a chemical environment in which an atom is located; Determining a candidate molecular conformation based on the rotatable bond and the torsion angle corresponding to the rotatable bond; The designated number of low potential energy molecular conformations are determined from the candidate molecular conformations, wherein the similarities between any two of the designated number of low potential energy molecular conformations are less than a similarity threshold.

4. The method according to claim 1, characterized in that: The obtaining of M clusters comprises: Obtain supercells of experimental crystal structures; A molecular dynamics simulation is performed on the supercell of the experimental crystal structure to obtain a stable supercell structure, so as to obtain M clusters based on the stable supercell structure.

5. The method according to claim 3, characterized in that: After determining the designated number of low potential energy molecular conformations from the candidate molecular conformations, the method further comprises: When the specified number of types is 1, the molecular conformation with the lowest potential energy is selected from the candidate molecular conformations as the global optimal molecular conformation; When the specified number of types is multiple, the parameter values ​​of the low potential energy molecular conformations of the specified number of types are weighted averaged to obtain the global optimal molecular conformation.

6. The method according to claim 5, characterized in that ； The candidate box is formed by stacking a plurality of the global optimal molecular conformations.

7. The method according to claim 2 or 4, characterized in that: The M clusters form at least two layers to wrap the central molecule.

8. The method according to claim 5, characterized in that Also includes: For a multi-component system, the molecules of each component system are stacked according to the proportion of each component of the multi-component system to obtain a box; The obtaining of M clusters includes: performing dynamic sampling on the box to obtain M clusters.

9. The method according to claim 4, characterized in that The obtaining of M clusters includes: performing dynamic sampling on the supercell to obtain M clusters.

10. A method for optimizing force field parameters, characterized in that: include: Determining bonding parameters, and determining the non-bonding parameters according to the method according to any one of claims 1 to 9; The force field parameters are obtained based on the bonding parameters and the non-bonding parameters.

11. The method according to claim 10, characterized in that: The bonding parameters include at least one of a bonding parameter, an angle parameter and a dihedral angle parameter; After obtaining the non-key parameter, the method further includes: Calculating the quantized energy of each scanning point to obtain potential energy surface data, wherein the scanning point is a point scanned during a potential energy surface scanning of at least a portion of the flexible angle of the target molecular conformation; Using the potential energy surface data as training data, performing molecular one-dimensional potential energy surface fitting based on the non-bonding parameters and the bonding parameters to obtain fitting dihedral angle parameters; The bond parameter, the angle parameter and the fitted dihedral angle parameter are used as optimized bonding parameters.

12. The method according to claim 10, characterized in that Each cluster in the M clusters is formed by stacking target molecular conformations, and the target molecular conformations include: experimental molecular conformations and / or global optimal molecular conformations.

13. The method according to claim 12, characterized in that The bond parameters and angle parameters in the bonding parameters are determined as follows: Determining a Hessian matrix of the experimental molecular conformation or the global optimal molecular conformation; Converting the Hessian matrix into force constants of each bond and each angle of the target molecular conformation; and determining the molecular bond length and molecular bond angle of the target molecular conformation; The force constant of each bond of the target molecular conformation and the corresponding molecular bond length are used as the bond parameters, and the force constant of each angle of the target molecular conformation and the corresponding molecular bond angle are used as the angle parameters.

14. A method for calculating crystal free energy, characterized in that: The method comprises: The non-bonded parameters are obtained by using the method for optimizing the non-bonded parameters of the force field according to any one of claims 1 to 9, or the force field parameters are obtained by using the method for optimizing the force field parameters according to any one of claims 10 to 13; The free energy of the crystal structure is calculated based on the non-bonded parameters and / or the force field parameters.

15. A method for predicting a crystal form, characterized in that: The method comprises: The free energy of various crystal structures is obtained by using the crystal free energy calculation method described in claim 14; The target crystal form is determined based on the free energy of the multiple crystal structures.

16. A design method, characterized in that: The method comprises: The method according to any one of claims 10 to 13, optimizing force field parameters; Drug design or material design is performed based on the optimized force field parameters.

17. A device for optimizing non-bonded parameters of a force field, characterized in that: include: A cluster acquisition module, used to obtain M clusters, each of the M clusters includes N molecules, and M and N are integers greater than or equal to 1; A non-bonded parameter determination module is used to determine the quantitative intermolecular interaction energy of each polymer of at least some of the M clusters; for each cluster in the M clusters, the electrostatic potential of the central molecule of the cluster is fitted based on the background charge of the cluster to obtain the atomic charge of each atom in the central molecule, and then the charge parameters of the force field are obtained based on the atomic charge of each atom in the central molecule of the M clusters; wherein the background charge is a charge determined by the atomic charge of all molecules in the cluster except the central molecule, and the central molecule is a molecule at the center of the cluster; based on the quantitative intermolecular interaction energy of each polymer of the cluster, the atomic charge of each atom in the central molecule is used as the initial charge of each atom in the molecule in the polymer, the bonding parameters and the initial charge are fixed unchanged, and the van der Waals parameters of the molecule are obtained by fitting.

18. A device for optimizing force field parameters, characterized in that: include: a force field parameter determination module, configured to determine bonding parameters, and to determine the non-bonding parameters according to the apparatus of claim 17; A force field parameter optimization module is used to obtain the force field parameters based on the bonding parameters and the non-bonding parameters.

19. A crystal free energy calculation device, characterized in that: The device comprises: A force field parameter acquisition module, used to obtain the non-bonded parameters by using the device for optimizing the non-bonded parameters of the force field according to claim 17, or to obtain the force field parameters by using the device for optimizing the force field parameters according to claim 18; A crystal structure free energy obtaining module is used to calculate and obtain the free energy of the crystal structure based on the non-bonded parameters and / or the force field parameters.

20. A crystal form prediction device, characterized in that: The device comprises: A multi-free energy calculation module, used to obtain the free energy of multiple crystal structures using the crystal free energy calculation device according to claim 19; The target crystal form determination module is used to determine the target crystal form according to the free energy of the multiple crystal structures.

21. A design device, characterized in that: include: A parameter optimization module, configured to optimize force field parameters of the apparatus according to claim 18; The design module is used for drug design or material design based on optimized force field parameters.

Images