CN106372400A - Method and application for constructing polarized force fields and method and system for predicting drug crystal forms - Google Patents

Abstract

The invention discloses a method and application for constructing polarized force fields and a method and system for predicting drug crystal forms. The method for constructing chemical molecule polarized force fields is suitable for being executed in one or more computation apparatuses, and comprises the following steps of: carrying out optimization computation on a quantum chemical structure of a chemical molecule on the basis of an ab initio calculation method so as to obtain a locally optimized molecular structure; calculating the difference between first energy of the locally optimized molecular structure when the charge of the molecular structure is neutral and second energy of the locally optimized molecular structure when the molecular structure has a predetermined positive-valence charge, and taking the difference as a vertical ion potential corresponding to the molecular structure; calculating polarized force field parameters of the molecular structure on the basis of the vertical ion potential and the locally optimized molecular structure, wherein the polarized force field parameters comprise a multi-pole vector of atom distribution, a multi-pole polarization rate of the atom distribution and a frequency-related polarization rate; and constructing a corresponding polarized force field model on the basis of the locally optimized molecular structure and the calculated polarized force field parameters.

Description

Build the method for Polarized force field and the method and system of application, prediction drug crystal forms

Technical field

The present invention relates to quantum chemistry calculation field, more particularly, to build method and application, the prediction medicine of Polarized force field The method and system of crystal formation.

Background technology

Traditional compound design (predictions of such as medicine novel crystal forms) lacks ripe discovery approach.For example, averagely to Up to ten thousand kinds of compounds of screening just can obtain a kind of new drug.

With the development of high-performance calculation (hpc) and calculating chemical technology, compound crystal form (particularly drug crystal forms) Prediction increasingly receives publicity.However, crystal Predicting Technique is current or world-class technical barrier.Wherein, the structure in the field of force is The key technique of crystal formation prediction.

Currently, with respect to the field of force structure can by quantum chemistry Ah initiomethod or apply general force field model Lai Realize.Wherein, quantum-mechanical effects method requires a great deal of time, and receives very big restriction in practical application.Logical It is approximate according to quantum-mechanical Wo En-Oppenheimer with force field model, apply multiple approximate force field parameter (such as field of forces Parameter can be experience estimation value), its computational accuracy has much room for improvement.

For this reason, this proposes a kind of new technical scheme building Polarized force field and prediction drug crystal forms clearly.

Content of the invention

For this reason, the present invention provides a kind of new technical scheme building Polarized force field and prediction drug crystal forms, effectively solve At least one problem above having determined.

According to an aspect of the present invention, provide a kind of method building chemical molecular Polarized force field, be suitable at one or Execute in multiple computing devices.The method comprises the steps.Based on ab initio calculation (ab initio), to chemical molecular Quantum chemistry structure is optimized calculating, to obtain the molecular structure of local optimum.The molecular structure calculating this local optimum exists The first energy during neutral charge and the difference of the second energy when band makes a reservation for positive valence charge, and using this difference as this molecular structure Corresponding vertical ion gesture.Based on the molecular structure of this vertical ion gesture and this local optimum, calculate the polarization of this molecular structure Force field parameter.This Polarized force field parameter includes multistage vector, the multistage polarizability of atom distribution and the frequency dependence of atom distribution Polarizability.Molecular structure based on local optimum and calculated Polarized force field parameter, build corresponding Polarized force field model.

According to a further aspect of the invention, provide a kind of method of drug crystal forms prediction, comprise the steps.Root first Build the Polarized force field model with regard to drug molecule according to the method for the present invention.Then, according to constructed Polarized force field model, meter Calculate the crystal formation of this drug molecule.

According to a further aspect of the invention, provide a kind of application building chemical molecular Polarized force field, be suitable to be distributed in In one or more computing devices.This application include optimization calculation unit, ionic potential computing unit, force field parameter computing unit and Field of force construction unit.Based on ab initio calculation (ab initio), optimization calculation unit is suitable to the quantum chemistry knot to chemical molecular Structure is optimized calculating, to obtain the molecular structure of local optimum.Ionic potential computing unit is suitable to calculate dividing of this local optimum First energy in neutral charge for the minor structure and the difference of the second energy when band makes a reservation for positive valence charge, and using this difference as dividing Minor structure corresponding vertical ion gesture.Based on the molecular structure of this vertical ion gesture and this local optimum, force field parameter calculates single Unit is suitable to calculate the Polarized force field parameter of this molecular structure.This Polarized force field parameter includes the multistage vector of atom distribution, atom The multistage polarizability of distribution and frequency dependence polarizability.Field of force construction unit is suitable to molecular structure and calculating based on local optimum The Polarized force field parameter obtaining, builds corresponding Polarized force field model.

According to another aspect of the present invention, provide a kind of system of drug crystal forms prediction, including the structure according to the present invention The application of chemical molecular Polarized force field and crystal formation prediction application.Crystal formation prediction application is suitable to according to constructed Polarized force field mould Type, calculates the crystal formation of this drug molecule.

To sum up, different from traditional general field of force, according to the present invention structure Polarized force field technical scheme by using The partial structurtes of the chemical molecular to standard for the ab iitio mode are optimized, and build Polarized force field model on this basis, The chemical precision higher than the general field of force can be obtained.Particularly, the scheme in the structure Polarized force field of the present invention is used in During drug crystal forms prediction field, crystal formation Predicting Technique scheme can provide accurate crystalline substance based on high-precision Polarized force field Type design direction, thus be greatly improved the efficiency finding novel crystal forms and success rate.

Brief description

In order to realize above-mentioned and related purpose, herein in conjunction with explained below and accompanying drawing, some illustrative sides to be described Face, these aspects indicate can be to put into practice the various modes of principles disclosed herein, and all aspects and its equivalent aspect It is intended to fall under in the range of theme required for protection.By detailed description below be read in conjunction with the accompanying, the disclosure above-mentioned And other purpose, feature and advantage will be apparent from.Throughout the disclosure, identical reference generally refers to identical Part or element.

Fig. 1 shows the schematic diagram of some computing environment 100 implementing to arrange according to the present invention；

Fig. 2 shows the schematic diagram of the application 200 of structure chemical molecular Polarized force field according to some embodiments of the invention；

The flow chart that Fig. 3 shows the method 300 of structure chemical molecular Polarized force field according to some embodiments of the invention；

The flow chart that Fig. 4 shows the implementation 400 of step s330 in Fig. 3；And

The flow chart that Fig. 5 shows the implementation 500 of step s340 in Fig. 3.

Specific embodiment

It is more fully described the exemplary embodiment of the disclosure below with reference to accompanying drawings.Although showing the disclosure in accompanying drawing Exemplary embodiment it being understood, however, that may be realized in various forms the disclosure and should not be by embodiments set forth here Limited.On the contrary, these embodiments are provided to be able to be best understood from the disclosure, and can be by the scope of the present disclosure Complete conveys to those skilled in the art.

Fig. 1 shows the schematic diagram of some cloud computing environments 100 implementing to arrange according to the present invention.

As shown in figure 1, computing environment 100 includes multiple computing devices.Each computing device is an equipment section in cluster Point.Here, computing environment 100 both can be enterprises lan network environment or cloud computing environment (for example surpasses in calculating The heart), the present invention does not do excessive restriction to this.Computing environment 100 for example can include calculate node 110,120,150 and 160, also Memory node 130 and 140 can be included, but not limited to this.Taking calculate node 150 and 160 as a example, the two can execute meter respectively Calculate application 151 and 161.Calculate application 151 and 161 and can be configured to the structure chemical molecular Polarized force field according to the present invention A part for application.In other words, the application of the structure chemical molecular Polarized force field of the present invention can be configured in computing environment The Distributed Application of execution in 100.With reference to Fig. 2, the application of the structure chemical molecular Polarized force field according to the present invention is carried out Describe in detail.

Fig. 2 shows the schematic diagram of the application 200 of structure chemical molecular Polarized force field according to some embodiments of the invention. Application 200 for example can execute in the computing environment 100 shown in Fig. 1, but not limited to this.

Application 200 includes optimization calculation unit 210, ionic potential computing unit 220, force field parameter computing unit 230 and power Field construction unit 240.

Based on ab initio calculation (ab initio), optimization calculation unit 210 is suitable to the quantum chemistry structure to chemical molecular It is optimized calculating, to obtain the molecular structure of local optimum.In an embodiment in accordance with the invention, ab initio calculation is by reality Apply as b3lyp/6-31g (d, p) functional mode.

Wherein, 6-31g (d, p) is the basis set of description atom,

6 is the description to inner orbit, with one slater function (sto) of 6 Gaussian functions (gto) matching,

31 represent virtual orbital are described with two groups of sto, every group of sto respectively with 3 and 1 gto matching,

D is the polarization function that non-h atom is done,

P is the polarization function that h atom is done.

In addition, multiple known (outside b3lyp/6-31g (d, p) functional) can also be applied according to embodiments of the invention Ab iitio mode, the present invention do not do excessive restriction to this.

On the basis of the above, optimization calculation unit 210 can be to the quantum chemistry under known chemical molecular gas phase configuration Structure carries out local optimum calculating.With careless acid dihydrazide (oxalyl dihydrazide, molecular formula: c2h6n4o2) it is below Example, illustrates to optimization calculation unit 210 more specifically implementation process.

The data instance of the chemical molecular structure to be optimized of careless acid dihydrazide is as follows:

Optimization calculation unit 210 can be configured to apply quantum calculation software (such as qchem, but not limited to this) and calculate The local optimum structure of chemical molecular.Here, optimization calculation unit 210 can obtain chemical molecular in local in calculating process Multiple new configuration.Optimization calculation unit 210 can select the knot of energy low (usually minimum) in obtained new configuration Structure is as local optimum structure.So, optimization calculation unit 210, by above-mentioned calculating and selection, obtains through local optimum Molecular structure.Here, the molecular structure through local optimum can be represented as cartesian coordinate data, but not limited to this.

The molecular structure of the local optimum being obtained for optimization calculation unit 210, ionic potential computing unit 220 is suitable to count Calculate the first energy in neutral charge for this molecular structure and the second energy in predetermined positive valence charge (for example positive 1 valence charge) of band The difference of amount, and using this difference as this molecular structure corresponding vertical ion gesture.For example, the local optimum of mesoxalic acid two hydrazides above The vertical ion gesture of structure is 0.370100.

The local being obtained based on the vertical ion gesture obtained by ionic potential computing unit 220 and optimization calculation unit 210 is excellent The molecular structure changed, force field parameter computing unit 230 is suitable to calculate the Polarized force field parameter of molecular structure.Here, Polarized force field Parameter includes multistage vector, the multistage polarizability of atom distribution and the frequency dependence polarizability of atom distribution, but not limited to this.

In an embodiment in accordance with the invention, force field parameter computing unit 230 will comprise the molecular structure of local optimum The data of data (the cartesian coordinate data of such as molecular structure) and vertical ion gesture is as input data.For example, field of force ginseng Number computing unit 230 can apply camcasp software (cambridge package for calculation of Anisotropic site properties, calculates Cambridge software kit of anisotropy website attribute) calculated, but do not limit In this.The input data example of the molecular structure of the local optimum of careless acid dihydrazide is as follows:

In turn below the process calculating 3 kinds of Polarized force field parameters is illustrated.

The multistage vector of atom distribution is the multiple development of electron density in molecular structure.In one embodiment, the field of force Parameter calculation unit 230 carries out spherical-harmonic expansion to electron density in molecular structure and obtains multistage vector.According to the present invention In another embodiment, the multistage vector of atom distribution of the calculated grass acid dihydrazide of force field parameter computing unit 230 Partial data form is as follows:

Multi-stage static polarizability represents the electron density of this molecular structure responding ability under External Electrical Field.Force field parameter Computing unit 230 can carry out multiple development operation to the function representing atom site in this molecular structure, many to obtain static state Level polarizability.In an embodiment in accordance with the invention, the calculated grass acid dihydrazide of force field parameter computing unit 230 is former The partial data form of the multistage polarizability of son distribution is as follows:

It should be noted that if external electrical field is not steady electric field, the response in external electric field of molecule charge is frequency Related., force field parameter computing unit 230 calculates the multistage pole of frequency dependence of the atom distribution of gained taking careless acid dihydrazide as a example Rate partial data form is as follows:

As described above, force field parameter computing unit 230 can be in optimization calculation unit 210 and ionic potential computing unit 220 On the basis of fetched data, high accuracy builds 3 kinds of force field parameters calculating Polarized force field model.On this basis, field of force structure Build unit 240 and be suitable to the molecular structure based on local optimum and calculated Polarized force field parameter, build corresponding polarization force Field model.Here, build Polarized force field and be mainly the every energy calculating with regard to molecular structure.According to one enforcement of the present invention In example, the energy that field of force construction unit 240 is calculated includes the corresponding multinomial intermolecular lattice energy of molecular crystal of molecular structure Summation, and the self-energy of molecular structure.

In an embodiment in accordance with the invention, field of force construction unit 240 can adopt Density functional tight-binding (dftb) Mode calculates the self-energy of molecular structure.In addition, can also be calculated using mode known to other according to embodiments of the invention The self-energy of molecular structure, repeats no more here.

In addition, field of force construction unit 240 can calculate the molecular structure of local optimum according to above-mentioned Polarized force field parameter Corresponding molecular crystal, the summation of the multinomial intermolecular lattice energy under periodic boundary condition.According to the present invention one In embodiment, multinomial intermolecular lattice energy includes representing the energy e of intermolecular static state electrostatic interaction_lattice-es, represent point Polarize between son the energy e of interaction_lattice-indWith the energy e representing that intermolecular dispersion interacts_lattice-disp.

Wherein, field of force construction unit 240 can calculate e according to following formula (1), (2) and (3)_lattice-es.

$\begin{array}{c}{u}_{lattice-es}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2\end{array};---\left(1\right)$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);---\left(2\right)$

Wherein,It is based on lower Ewald summation (ewald summation) of the multistage vector of atom distribution Turned letter interphase interaction matrix.Generally, the dimension of interaction matrix can be for example 25*25, but not limited to this.A, b representThe index of two dimensions of matrix.

q_tAnd q_uThe atom exactly representing the molecule of two interactions in lattice is distributed multistage vector.T and u represent be The index of multistage vector.N is the index of the cyclic extension of ewald summation.K is the vector basis based on reciprocal lattice vector.S square What the direction type that the atom that battle array represents in two molecules in lattice is formed when interacting was constituted is spatially steered matrix.F is down Convergence parameter equation under space.L is the exponent number of multiple development.It is based on the lower Ai Waer of the multistage vector of atom distribution The real space interaction matrix of moral summation.R is the vector basis based on real space lattice vector.Expression be Interaction matrix under construction is between turned letter during reciprocal lattice vector kn=0.It is that in the multistage vector of atom distribution, dimension indexes as 0 When value,Represent self-interaction.Represent the atom being derived from two molecules in lattice Form the arrow composition of direction when interacting is spatially steered matrix.S is the side based on two atom position constructions in molecule To vector.ω_a, ω_b, it is the vector matrix based on selected molecule local coordinate system construction for two atom positions respectively.It is based on Reciprocal lattice vector under between turned letter.It is the vector arrow that two atom positions in two molecules are constituted.It is in real space next week Real space lattice vector after the extension of phase property boundary condition.Convergent function under between for turned letter. For the convergent function under the real space.

In addition, representing the energy e that intermolecular polarization interacts_lattice-indIncluding 2 parts, that is, the polarization restraining (induces) Many body polarizations energy in multistage vector sum lattice.Below the calculating process of 2 portion of energy is illustrated.

On the one hand, field of force construction unit 240 can iterate to calculate the polarization restraining multistage according to following formula (4) and (5) Vector.

${\mathrm{\δq}}_t^a=-\underset{b\&notequal;a}{\σ}\underset{b}{\σ}\underset{a^{\′}}{\σ}\underset{t^{\′}u}{\σ}{\mathrm{\α}}_{{\mathrm{tt}}^{\′}}^{{\mathrm{aa}}^{\′}}{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}(q_u^b+{\mathrm{\δq}}_u^b)---\left(4\right)$

$e_{ind}=\frac{1}{2}\underset{a}{\σ}\underset{b\&notequal;a}{\σ}{\mathrm{\δq}}_t^a{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}{q}_u^b---\left(5\right)$

Wherein,Respectively induce multistage vector,Then for the multistage polarizability of atom distribution.b≠a Represent that processing intramolecular polarization without formula (4) and (5) interacts.Because intramolecular polarization was processed with dftb, Therefore formula (4) and (5) are used for calculating intermolecular polarization interaction.f_n(r, β) is tang-toennies convergence equation, Specific formula for calculation is:

$f_n(r,\mathrm{\β})=1-\underset{k=0}{\overset{n}{\σ}}\left(\frac{{\left(\mathrm{\β}r\right)}^k}{k!}\right)e^{-\mathrm{\β}r}$

R is the vector field homoemorphism between two atom positions.N is convergence order of equation number, and the index of usual exponent number is the model of 1-5 Enclose and be obtained with restraining, but not limited to this.β is the related parameter of the property to molecule, can be by building three body groups The mode that cluster is fitted obtains.For example, β may be at the interval of [1.2-2.0].

On the other hand, field of force construction unit 240 can be by following calculating operations come many body polarizations energy in counting lattice.Tool For body, it is predetermined angstrom of degree (such as 25 that field of force construction unit 240 is based on molecular structure corresponding molecular crystal construction truncation radius Angstrom) cluster.Then, field of force construction unit 240 calculates the induction of the decay of each atom in this cluster by way of being certainly in harmony Multistage vectorWith the multistage vector of unattenuated inductionOn this basis, power Field construction unit 240 calculates according to formula (6) decays and the multistage vectorial difference of unattenuated induction:

${\mathrm{\δe}}_{damp}^{cluster}=e_{ind}^{cluster}\left(damped\right)-e_{ind}^{cluster}\left(undamped\right)---\left(6\right)$

In addition, field of force construction unit 240 can calculate according to formula (7), (8) and (9) induces multistage vectorLattice energy in Ewald sum formula

$\begin{array}{c}{e}_{ind}^{ewald}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2\end{array};---\left(7\right)$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);---\left(8\right)$

The parameter declaration of formula (7)-(9) refers to the parameter declaration of above formula (1), (2) and (3).

Furthermore, field of force construction unit 240 can be according to many body polarizations energy in formula (10) counting lattice:

$e_{ind}^{lattice}\left(damped\right)=e_{ind}^{ewald}\left(undamped\right)+{\mathrm{\δe}}_{damp}^{cluster}---\left(10\right)$

Above illustrate field of force construction unit 240 to the energy representing intermolecular static state electrostatic interaction e_lattice-esCalculating process (using formula (1), (2) and (3)) and represent intermolecular polarization interact energy e_lattice-indCalculating process (using formula (4) to (10)).

Continue with and the energy e representing that intermolecular dispersion interacts is calculated to field of force construction unit 240_lattice-disp's Process illustrates.Field of force construction unit 240 can calculate two bulk dispersion respectively can be with three bulk dispersion energy.Specifically, field of force structure Build unit 240 and build the truncation radius that dispersion can calculate first.For example, truncation radius are 15-25 angstrom, but not limited to this.It is based on Frequency dependence polarizability, field of force construction unit 250 can calculate abbe number c6 and c8 according to formula (11) and (12).

$c_6^{ab}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right)d\mathrm{\υ}---\left(11\right)$

$c_8^{ab}=\frac{15}{2\mathrm{\π}}\left({\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\right(i\mathrm{\υ}\left){\mathrm{\α}}_{22}^b\right(i\mathrm{\υ})d\mathrm{\υ}+{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{22}^a(i\mathrm{\υ}\left){\mathrm{\α}}_{11}^b\right(i\mathrm{\υ}\left)d\mathrm{\υ}\right)---\left(12\right)$

WhereinWithIt is the frequency dependence polarizability of the even level of induction of atom from two molecules respectively,WithIt is the frequency dependence polarizability of the three-level induction three-level of atom from two molecules respectively.

Based on calculated abbe number c6 and c8, field of force construction unit 240 can calculate disome according to formula (13) Dispersion energy:

$e_{2\·bodydisp}=-\underset{a}{\σ}\underset{b\&notequal;a}{\σ}\underset{a}{\σ}\underset{b}{\σ}(f_6\frac{c_6^{ab}}{r_{ab}^6}+f_8\frac{c_8^{ab}}{r_{ab}^8}+\mathrm{...})---\left(13\right)$

Wherein, f₆And f₈Represent above it is stated that tang-toennies convergence equation.6 and 8 is order of convergence number. r_abIt is the distance of the atom interdigit in two molecules.

Field of force construction unit 240 can calculate three bulk dispersion coefficient c9 according to formula (14).

$c_9^{abc}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^c\left(i\mathrm{\υ}\right)d\mathrm{\υ}---\left(14\right)$

Represent the frequency of three not homoatomic even level induction idol levels in two molecules Rate correlation polarizability.

Based on calculated three bulk dispersion coefficient c9, field of force construction unit 240 can calculate three bodies according to formula (15) Dispersion energy.

$e_{3\·bodydisp}^{abc}=\underset{a\&element;a}{\σ}\underset{b\&element;b}{\σ}\underset{c\&element;c}{\σ}{f}_9{c}_9^{abc}\frac{(1+3\cos \hat{a}\cos \hat{b}\cos \hat{c})}{r_{ab}^3{r}_{bc}^3{r}_{ac}^3}---\left(15\right)$

Wherein a, b and c refer to three molecules, and a, b and c refer to the atom in molecule,WithRefer to three atoms and form triangle The angle at three angles after shape.

, the Polarized force field result constructed by application 200 of the present invention is as follows taking careless acid dihydrazide as a example:

Qm 1-body=-896.247339368hartrees

Qm 2-body=-0.142601938hartrees

Mm full=-763.941kj/mol

Mm 1-body=0.000kj/mol

Mm 2-body=-775.516kj/mol

Mm many-body=11.576kj/mol

Qm 1+2-body energy=-896.389941306hartrees

Mm many-body energy=11.576kj/mol (5.788kj/mol per monomer)

Hmbi energy=-896.385532326hartrees

The above results show, the reference energy of experiment is -896.3934378hartrees, and error, in several kj/mol, reaches Chemical precision.

To sum up, different from traditional general field of force, according to the present invention structure Polarized force field application 200 by using from The partial structurtes of the head chemical molecular to standard for the calculation are optimized, and build Polarized force field model on this basis, can To obtain the chemical precision higher than the general field of force.Particularly, when applying 200 to be used in drug crystal forms prediction field, brilliant Type Predicting Technique scheme can provide accurate crystal formation design direction, thus greatly carrying based on high-precision Polarized force field The efficiency of high discovery novel crystal forms and success rate.

The flow chart that Fig. 3 shows the method 300 of structure chemical molecular Polarized force field according to some embodiments of the invention. Method 300 is suitable to execution in one or more computing devices.In other words, method 300 can be in cloud computing environment or internal collection Execute in group.

As shown in figure 3, method 300 starts from step s310.In step s310, based on ab initio calculation (ab initio), Calculating is optimized to the quantum chemistry structure of chemical molecular, to obtain the molecular structure of local optimum.According to the present invention one Embodiment, step s310 can utilize b3lyp/6-31g (d, p) functional mode, and the quantum chemistry structure of this chemical molecular is entered Row optimizes and calculates, and using the low molecular structure of local energy calculating gained as described local optimum molecular structure.Step S310 more specifically embodiment is consistent with above optimization calculation unit 210, repeats no more here.

Obtain the molecular structure of local optimum in step s310 after, method 300 execution step s320.In step s320 In, calculate first energy in neutral charge for the molecular structure of local optimum and make a reservation for positive valence charge (for example positive 1 valency electricity in band Lotus) when the second energy difference, and using this difference as this molecular structure corresponding vertical ion gesture.

Molecular structure based on the local optimum obtaining in step s310 and step s320 calculate gained vertical ion gesture, side Method 300 execution step s330.In step s330, calculate the Polarized force field parameter of molecular structure.Here, Polarized force field parameter example As multistage vector, the multistage polarizability of atom distribution and the frequency dependence polarizability of atom distribution can be included.According to the present invention One embodiment, step s330 may be implemented as method 400 as shown in Figure 4.

As shown in figure 4, method 400 starts from step s410, spherical-harmonic expansion is carried out to electron density in this molecular structure, To obtain the multistage vector of this atom distribution.Method 400 also includes step s420, to the letter representing atom site in molecular structure Number carries out multiple development operation, with obtain represent the electron density of this molecular structure under External Electrical Field responding ability multistage Polarizability.In addition, method 400 also execution step s430, obtain the frequency dependence polarizability of atom distribution in molecular structure.Step S330 more specifically embodiment is consistent with above force field parameter computing unit 230, repeats no more here.

Molecular structure based on step s310 gained local optimum and step s330 gained Polarized force field parameter, method 300 Execution step s340.In step s340, build the analytical structure corresponding Polarized force field model of local optimum.Specifically, exist Calculate the self-energy of the molecular structure of local optimum using Density functional tight-binding (dftb) mode in step s340.In addition, In step s340 according to Polarized force field parameter calculate local optimum molecular structure corresponding to molecular crystal, in periodic boundary Under the conditions of multinomial intermolecular lattice energy summation.

In an embodiment in accordance with the invention, multinomial intermolecular lattice energy includes representing intermolecular static state electrostatic phase interaction Energy e_lattice-es, represent intermolecular polarization interact energy e_lattice-indIntermolecular dispersion interacts with expression Energy e_lattice-disp.Step s340 may be implemented as the method 500 shown in Fig. 5.

As shown in figure 5, method 500 starts from step s510.In step s510, according to atom distribution multistage vector, with And in lattice the mutual molecule of any two position coordinateses, calculate represent intermolecular static state electrostatic interaction energy e_lattice-es.Here, step s510 can calculate e according to above formula (1), (2) and (3)_lattice-es, repeat no more here.

Method 500 also includes step s520, calculates the energy representing that intermolecular polarization interacts.Here, represent molecule Between polarize the energy e of interaction_lattice-indEnergy include 2 parts, that is, the multistage vector sum lattice of polarization (induction) restrained Interior many body polarizations energy.Step s520 can calculate the multistage vector of polarization of convergence according to above formula (4) and (5), here not Repeat again.Step s520 can repeat no more here according to many body polarizations energy in above formula (6) to (10) counting lattice.

Method 500 also includes step s530, calculates the energy e representing that intermolecular dispersion interacts_lattice-disp.According to One embodiment of the invention, step s530 calculates two bulk dispersion respectively can be with three bulk dispersion energy.In order to calculate two bulk dispersion energy, The truncation radius that dispersion can calculate are built first in step s530, such as 15-25 angstrom.Then, step s530 can be according to above Middle formula (11), (12) and (13) calculates two bulk dispersion energy, repeats no more here.In addition, step s530 can be according to above Formula (14) and (15) calculate three bulk dispersion energy, repeat no more here.It should be noted that step s340 more specifically embodiment is permissible See above middle field of force construction unit 240, repeat no more here.

A8, the method as any one of a5-a7, wherein, the described energy calculating expression intermolecular dispersion interaction Amount e_lattice-dispStep include:

Calculate two bulk dispersion energy；

Calculate three bulk dispersion energy.

A9, the method as described in a8, wherein, the described step calculating two bulk dispersion energy includes:

Build the truncation radius that dispersion can calculate；

Based on described frequency dependence polarizability, and according to following manner calculating abbe number c6 and c8:

$c_6^{ab}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right)d\mathrm{\υ}$

$c_8^{ab}=\frac{15}{2\mathrm{\π}}\left({\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\right(i\mathrm{\υ}\left){\mathrm{\α}}_{22}^b\right(i\mathrm{\υ})d\mathrm{\υ}+{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{22}^a(i\mathrm{\υ}\left){\mathrm{\α}}_{11}^b\right(i\mathrm{\υ}\left)d\mathrm{\υ}\right)$

Wherein.WithRepresent the frequency dependence polarization of the even level of induction of atom from two molecules respectively Rate.WithRepresent that the three-level of the atom from two molecules induces the frequency dependence polarizability of three-level respectively,

Based on abbe number c6 and c8, and according to following formula calculating two bulk dispersion energy:

$e_{2\·bodydisp}=-\underset{a}{\σ}\underset{b\&notequal;a}{\σ}\underset{a}{\σ}\underset{b}{\σ}(f_6\frac{c_6^{ab}}{r_{ab}^6}+f_8\frac{c_8^{ab}}{r_{ab}^8}+\mathrm{...})$

Wherein, f₆And f₈Represent tang-toennies convergence equation.6 and 8 is order of convergence number, r_abRepresent and be derived from two The distance of the atom interdigit in molecule.

A10, the method as described in a8 or a9, wherein, the step calculating three bulk dispersion energy includes:

According to following formula calculating three bulk dispersion coefficient c9:

$c_9^{abc}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^c\left(i\mathrm{\υ}\right)d\mathrm{\υ}$

Wherein,Represent that three not homoatomic even level inductions are even in two molecules The frequency dependence polarizability of level；

Based on calculated three bulk dispersion coefficient c9, according to following manner calculating three bulk dispersion energy:

$e_{3\·bodydisp}^{abc}=\underset{a\&element;a}{\σ}\underset{b\&element;b}{\σ}\underset{c\&element;c}{\σ}{f}_9{c}_9^{abc}\frac{(1+3\cos \hat{a}\cos \hat{b}\cos \hat{c})}{r_{ab}^3{r}_{bc}^3{r}_{ac}^3}$

Wherein, a, b and c refer to three molecules, and a, b and c refer to the atom in molecule,WithRepresent that three atoms are formed The angle at three angles after triangle.

B13, the application as described in b12, wherein, described optimization calculation unit is suitable to be based on ab iitio according to following manner Method (ab initio), is optimized calculating to the quantum chemistry structure of chemical molecular, to obtain the molecular structure of local optimum: Using b3lyp/6-31g (d, p) functional mode, calculating is optimized to the quantum chemistry structure of this chemical molecular, and will calculate The low molecular structure of the local energy of gained is as the molecular structure of described local optimum.B14, answering as described in b12 or b13 With wherein, described force field parameter computing unit is suitable to calculate the Polarized force field parameter of described molecular structure according to following manner: right In this molecular structure, electron density carries out spherical-harmonic expansion, to obtain the multistage vector of this atom distribution；To this molecule of expression In structure, the function of atom site carries out multiple development operation, is made in external electric field with obtaining the electron density representing this molecular structure Multistage polarizability with the atom distribution of lower responding ability；Obtain the frequency dependence polarizability of atom distribution in this molecular structure. B15, the application as any one of b12-b14, wherein, the described molecular structure based on local optimum and calculated pole Change force field parameter, the step building corresponding Polarized force field model includes: calculates using Density functional tight-binding (dftb) mode The self-energy of the molecular structure of described local optimum；Calculate the molecular structure of described local optimum according to described Polarized force field parameter Corresponding molecular crystal, the summation of the multinomial intermolecular lattice energy under periodic boundary condition.B16, answering as described in b15 With wherein, described multinomial intermolecular lattice energy includes representing the energy e of intermolecular static state electrostatic interaction_lattice-es, represent The energy e that intermolecular polarization interacts_lattice-indWith the energy e representing that intermolecular dispersion interacts_lattice-disp, described According to the molecular crystal corresponding to the molecular structure that described Polarized force field parameter calculates described local optimum in periodic boundary bar The step of the summation of multinomial intermolecular lattice energy under part includes: according to the multistage vector of described atom distribution, and in lattice The position coordinateses of the mutual molecule of any two, calculate the energy e representing intermolecular static state electrostatic interaction_lattice-es；Meter Calculate the energy e representing that intermolecular polarization interacts_lattice-ind；Calculate the energy representing that intermolecular dispersion interacts e_lattice-disp.B17, the application as described in b16, wherein, the described energy calculating expression intermolecular static state electrostatic interaction e_lattice-esStep include following manner:

$\begin{array}{c}{u}_{lattice-es}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2\end{array};$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);$

Wherein,It is based on lower Ewald summation (ewald summation) of the multistage vector of atom distribution Turned letter interphase interaction matrix, a, b representThe index of two dimensions of matrix, q_tAnd q_uExactly represent two in lattice The atom of the molecule of individual interaction is distributed multistage vector, and what t and u represented is the index of multistage vector, and n is ewald The index of the cyclic extension of summation,

K is the vector basis based on reciprocal lattice vector, and the atom in two molecules in s matrix representative lattice is formed when interacting Direction type constitute be spatially steered matrix, f be turned letter between under convergence parameter equation, l is the exponent number of multiple development,Convergent function under between for turned letter,It is based on the lower Ewald summation of the multistage vector of atom distribution Real space interaction matrix, r is the vector basis based on real space lattice vector,Expression is in reciprocal lattice vector kn Interaction matrix under construction is between turned letter when=0,It is that in the multistage vector of atom distribution, dimension indexes as value when 0,Represent self-interaction,Represent the atom being derived from two molecules in lattice to be formed mutually What direction arrow during effect was constituted is spatially steered matrix, and s is the direction vector based on two atom position constructions in molecule, ω_aAnd ω_bIt is the vector matrix based on selected molecule local coordinate system construction for two atom positions respectively,Under between based on turned letter Reciprocal lattice vector,It is the vector arrow that two atom positions in two molecules are constituted,It is then on real space periodical side Real space lattice vector after boundary's condition extension,For the convergent function under the real space.

B18, the application as described in b16 or b17, wherein, described calculating represents the intermolecular energy polarizing and interacting e_lattice-indStep include:

The multistage vector of polarization according to the iterative calculation convergence of following formula:

${\mathrm{\δq}}_t^a=-\underset{b\&notequal;a}{\σ}\underset{b}{\σ}\underset{a^{\′}}{\σ}\underset{t^{\′}u}{\σ}{\mathrm{\α}}_{{\mathrm{tt}}^{\′}}^{{\mathrm{aa}}^{\′}}{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}(q_u^b+{\mathrm{\δq}}_u^b)$

$e_{ind}=\frac{1}{2}\underset{a}{\σ}\underset{b\&notequal;a}{\σ}{\mathrm{\δq}}_t^a{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}{q}_u^b$

Respectively induce multistage vector,Then for the multistage polarizability of atom distribution.f_n(r, β) table Show tang-toennies convergence equation；With

According to many body polarizations energy in following manner counting lattice:

The cluster that truncation radius are predetermined angstrom of degree is constructed based on corresponding molecular crystal, and is calculated by way of being certainly in harmony The multistage vector of the induction of the decay of each atom in this clusterWith the multistage vector of unattenuated induction

Calculated according to following manner and decay and the multistage vectorial difference of unattenuated induction:

${\mathrm{\δe}}_{damp}^{cluster}=e_{ind}^{cluster}\left(damped\right)-e_{ind}^{cluster}\left(undamped\right);$

Calculated according to following manner and induce multistage vectorLattice energy in Ewald sum formula

$\begin{array}{c}{e}_{ind}^{ewald}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2\end{array};$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);$

According to many body polarizations energy in following formula counting lattices:

$e_{ind}^{lattice}\left(damped\right)=e_{ind}^{ewald}\left(undamped\right)+{\mathrm{\δe}}_{damp}^{cluster}$

B19, the application as any one of b16-b18, wherein, described calculating represents what intermolecular dispersion interacted Energy e_lattice-dispStep include:

Calculate two bulk dispersion energy；

Calculate three bulk dispersion energy.

B20, the application as described in b19, wherein, the described step calculating two bulk dispersion energy includes:

Build the truncation radius that dispersion can calculate；

Based on described frequency dependence polarizability, and according to following manner calculating abbe number c6 and c8:

$c_6^{ab}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right)d\mathrm{\υ}$

$c_8^{ab}=\frac{15}{2\mathrm{\π}}\left({\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\right(i\mathrm{\υ}\left){\mathrm{\α}}_{22}^b\right(i\mathrm{\υ})d\mathrm{\υ}+{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{22}^a(i\mathrm{\υ}\left){\mathrm{\α}}_{11}^b\right(i\mathrm{\υ}\left)d\mathrm{\υ}\right)$

WhereinWithRepresent the frequency dependence polarization of the even level of induction of atom from two molecules respectively Rate,WithRepresent that the three-level of the atom from two molecules induces the frequency dependence polarizability of three-level respectively,

Based on abbe number c6 and c8, and according to following formula calculating two bulk dispersion energy:

$e_{2\·bodydisp}=-\underset{a}{\σ}\underset{b\&notequal;a}{\σ}\underset{a}{\σ}\underset{b}{\σ}(f_6\frac{c_6^{ab}}{r_{ab}^6}+f_8\frac{c_8^{ab}}{r_{ab}^8}+\mathrm{...})$

Wherein, f₆And f₈Represent tang-toennies convergence equation.6 and 8 is order of convergence number, r_abRepresent and be derived from two The distance of the atom interdigit in molecule.

B21, the application as described in b19 or b20, wherein, the step calculating three bulk dispersion energy includes:

According to following formula calculating three bulk dispersion coefficient c9:

$c_9^{abc}=\frac{3}{\mathrm{\π}}{\&integral;}_0^{\mathrm{\∞}}{\mathrm{\α}}_{11}^a\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^b\left(i\mathrm{\υ}\right){\mathrm{\α}}_{11}^c\left(i\mathrm{\υ}\right)d\mathrm{\υ}$

Wherein,WithRepresent that three not homoatomic even level inductions are even in two molecules The frequency dependence polarizability of level；

Based on calculated three bulk dispersion coefficient c9, according to following manner calculating three bulk dispersion energy:

$e_{3\·bodydisp}^{abc}=\underset{a\&element;a}{\σ}\underset{b\&element;b}{\σ}\underset{c\&element;c}{\σ}{f}_9{c}_9^{abc}\frac{(1+3\cos \hat{a}\cos \hat{b}\cos \hat{c})}{r_{ab}^3{r}_{bc}^3{r}_{ac}^3}$

Wherein a, b and c refer to three molecules, and a, b and c refer to the atom in molecule,WithRefer to three atom shapes The angle at three angles after triangularity.

In description mentioned herein, illustrate a large amount of details.It is to be appreciated, however, that the enforcement of the present invention Example can be put into practice in the case of not having these details.In some instances, known method, knot are not been shown in detail Structure and technology, so as not to obscure the understanding of this description.

Similarly it will be appreciated that in order to simplify the disclosure and help understand one or more of each inventive aspect, Above in the description to the exemplary embodiment of the present invention, each feature of the present invention is grouped together into single enforcement sometimes In example, figure or descriptions thereof.However, the method for the disclosure should be construed to reflect an intention that i.e. required guarantor The application claims of shield are than the feature more features being expressly recited in each claim.More precisely, as following As claims are reflected, inventive aspect is all features less than single embodiment disclosed above.Therefore, abide by The claims following specific embodiment are thus expressly incorporated in this specific embodiment, wherein each claim itself Separate embodiments as the present invention.

Those skilled in the art should be understood module or unit or the group of the equipment in example disclosed herein Part can be arranged in equipment as depicted in this embodiment, or alternatively can be positioned at and the equipment in this example In different one or more equipment.Module in aforementioned exemplary can be combined as a module or be segmented into multiple in addition Submodule.

Those skilled in the art are appreciated that and the module in the equipment in embodiment can be carried out adaptively Change and they are arranged in one or more equipment different from this embodiment.Can be the module in embodiment or list Unit or assembly be combined into a module or unit or assembly, and can be divided in addition multiple submodule or subelement or Sub-component.In addition to such feature and/or at least some of process or unit exclude each other, can adopt any Combination is to all features disclosed in this specification (including adjoint claim, summary and accompanying drawing) and so disclosed Where method or all processes of equipment or unit are combined.Unless expressly stated otherwise, this specification (includes adjoint power Profit requires, summary and accompanying drawing) disclosed in each feature can carry out generation by the alternative features providing identical, equivalent or similar purpose Replace.

Although additionally, it will be appreciated by those of skill in the art that some embodiments described herein include other embodiments In included some features rather than further feature, but the combination of the feature of different embodiment means to be in the present invention's Within the scope of and form different embodiments.For example, in the following claims, embodiment required for protection appoint One of meaning can in any combination mode using.

Additionally, some heres in described embodiment be described as can be by the processor of computer system or by executing Method or the combination of method element that other devices of described function are implemented.Therefore, have for implementing methods described or method The processor of the necessary instruction of element forms the device for implementing the method or method element.Additionally, device embodiment This described element is the example of following device: this device is used for implementing performed by the element of the purpose in order to implement this invention Function.

As used in this, unless specifically stated so, come using ordinal number " first ", " second ", " the 3rd " etc. Description plain objects are merely representative of the different instances being related to similar object, and are not intended to imply that the object being so described must Must have the time upper, spatially, sequence aspect or given order in any other manner.

Although the present invention is described according to the embodiment of limited quantity, benefit from above description, the art Interior it is clear for the skilled person that it can be envisaged that other embodiments in the scope of the present invention thus describing.Additionally, it should be noted that Language used in this specification primarily to the purpose of readable and teaching and select, rather than in order to explain or limit Determine subject of the present invention and select.Therefore, in the case of without departing from the scope of the appended claims and spirit, for this For the those of ordinary skill of technical field, many modifications and changes will be apparent from.For the scope of the present invention, to this It is illustrative and not restrictive for inventing done disclosure, and it is intended that the scope of the present invention be defined by the claims appended hereto.

Claims (10)

1. a kind of method building chemical molecular Polarized force field, is suitable to execution, the method bag in one or more computing devices Include:

Based on ab initio calculation (ab initio), calculating is optimized to the quantum chemistry structure of chemical molecular, to obtain local The molecular structure optimizing；

Calculate the molecular structure of this local optimum the first energy in neutral charge and when band makes a reservation for positive valence charge second The difference of energy, and using this difference as this molecular structure corresponding vertical ion gesture；

Based on the molecular structure of this vertical ion gesture and this local optimum, calculate the Polarized force field parameter of this molecular structure, this pole Change multistage vector, the multistage polarizability of atom distribution and the frequency dependence polarizability that force field parameter includes atom distribution；And

Molecular structure based on local optimum and calculated Polarized force field parameter, build corresponding Polarized force field model.

2. the method for claim 1, wherein described based on ab initio calculation (ab initio), the amount to chemical molecular Sub- chemical constitution is optimized calculating, is included with the step obtaining the molecular structure of local optimum:

Using b3lyp/6-31g (d, p) functional mode, calculating is optimized to the quantum chemistry structure of this chemical molecular, and will Calculate the molecular structure as described local optimum for the low molecular structure of local energy of gained.

3. method as claimed in claim 1 or 2, wherein, described based on this vertical ion gesture, calculate the polarization of this molecular structure Force field parameter, this Polarized force field parameter includes multistage vector, the multistage polarizability of atom distribution and the frequency dependence of atom distribution The step of polarizability includes:

Spherical-harmonic expansion is carried out to electron density in this molecular structure, to obtain the multistage vector of this atom distribution；

Multiple development operation is carried out to the function representing atom site in this molecular structure, to obtain the electricity representing this molecular structure The multistage polarizability of sub- density responding ability under External Electrical Field；And

Obtain the frequency dependence polarizability of atom distribution in this molecular structure.

4. the method as any one of claim 1-3, wherein, the described molecular structure based on local optimum and calculating The Polarized force field parameter arriving, the step building corresponding Polarized force field model includes:

Calculate the self-energy of the molecular structure of described local optimum using Density functional tight-binding (dftb) mode；

According to the molecular crystal corresponding to the molecular structure that described Polarized force field parameter calculates described local optimum, on periodicity side The summation of the multinomial intermolecular lattice energy under the conditions of boundary.

5. method as claimed in claim 4, wherein, described multinomial intermolecular lattice energy includes representing intermolecular static state electrostatic phase The energy e of interaction_lattice-es, represent intermolecular polarization interact energy e_lattice-indMutual with the intermolecular dispersion of expression The energy e of effect_lattice-disp, described according to described Polarized force field parameter calculate described local optimum molecular structure corresponding to The step of the summation of multinomial intermolecular lattice energy under periodic boundary condition for the molecular crystal include:

According to described atom distribution multistage vector, and in lattice the mutual molecule of any two position coordinateses, computational chart Show the energy e of intermolecular static state electrostatic interaction_lattice-es；

Calculate the energy e representing that intermolecular polarization interacts_lattice-ind；And

Calculate the energy e representing that intermolecular dispersion interacts_lattice-disp.

6. method as claimed in claim 5, wherein, the described energy calculating expression intermolecular static state electrostatic interaction e_lattice-esStep include following manner:

$\begin{array}{c}{u}_{lattice-es}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2;\end{array}$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);$

Wherein,It is the turned letter based on lower Ewald summation (ewald summation) of the multistage vector of atom distribution Interphase interaction matrix, a and b representsThe index of two dimensions of matrix, q_tAnd q_uExactly represent two in lattice The atom of the molecule interacting is distributed multistage vector, and what t and u represented is the index of multistage vector, and n is ewald The index of the cyclic extension of summation, k is the vector basis based on reciprocal lattice vector, in two molecules in s matrix representative lattice What the direction type that atom is formed when interacting was constituted is spatially steered matrix, f be between turned letter under convergence parameter equation, l is The exponent number of multiple development,Convergent function under between for turned letter,It is that multistage vector is distributed based on atom The real space interaction matrix of lower Ewald summation, r is the vector basis based on real space lattice vector, Expression is the interaction matrix under constructing in reciprocal lattice vector kn=0 between turned letter,It is to tie up in the multistage vector of atom distribution Degree index is value when 0,Represent self-interaction,Represent and in lattice, be derived from two molecules In atom formed interact when direction arrow constitute be spatially steered matrix, s is based on two atom positions in molecule The direction vector of construction, ω_aAnd ω_bIt is the vector matrix based on selected molecule local coordinate system construction for two atom positions respectively,Reciprocal lattice vector under between based on turned letter,It is the vector arrow that two atom positions in two molecules are constituted,It is then in reality Real space lattice vector after the periodical boundary condition extension of space,For the convergent function under the real space.

7. the method as described in claim 5 or 6, wherein, described calculating represents the intermolecular energy polarizing and interacting e_lattice-indStep include:

The multistage vector of polarization according to the iterative calculation convergence of following formula:

${\mathrm{\δq}}_t^a=-\underset{b\&notequal;a}{\σ}\underset{b}{\σ}\underset{a^{\′}}{\σ}\underset{t^{\′}u}{\σ}{\mathrm{\α}}_{{\mathrm{tt}}^{\′}}^{{\mathrm{aa}}^{\′}}{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}(q_u^b+{\mathrm{\δq}}_u^b)$

$e_{ind}=\frac{1}{2}\underset{a}{\σ}\underset{b\&notequal;a}{\σ}{\mathrm{\δq}}_t^a{f}_n(r,\mathrm{\β})t_{t^{\′}u}^{a^{\′}b}{q}_u^b$

Respectively induce multistage vector,Then for the multistage polarizability of atom distribution, f_n(r, β) represents Tang-toennies restrains equation；With

According to many body polarizations energy in following manner counting lattice:

The cluster that truncation radius are predetermined angstrom of degree is constructed based on corresponding molecular crystal, and is calculated by way of being certainly in harmony in this cluster The multistage vector of induction of the decay of each atomWith the multistage vector of unattenuated induction

Calculated according to following manner and decay and the multistage vectorial difference of unattenuated induction:

${\mathrm{\δe}}_{damp}^{cluster}=e_{ind}^{cluster}\left(damped\right)-e_{ind}^{cluster}\left(undamped\right);$

Calculated according to following manner and induce multistage vectorLattice energy in Ewald sum formula

$\begin{array}{c}{r}_{ind}^{ewald}=\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a({\mathrm{recipt}}_{tu}^{ab}+{\mathrm{direct}}_{tu}^{ab}){\tilde{q}}_u^b\\ +\frac{1}{2}\underset{a}{\σ}\underset{b}{\σ}\underset{n=0}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^{a}recip\_kn0t_{tu}^{ab}{\tilde{q}}_u^b\\ -\frac{1}{2}\underset{m}{\σ}\underset{a\&element;m}{\σ}\underset{b\&element;m}{\σ}\underset{tu}{\σ}\mathrm{\δ}{\tilde{q}}_t^a{t}_{tu}^{ab}{\tilde{q}}_u^b-\frac{\mathrm{\γ}}{\sqrt{\mathrm{\π}}}\underset{n}{\σ}{\tilde{q}}_{00}^2;\end{array}$

${\mathrm{recipt}}_{tu}^{ab}=\left(\begin{array}{c}{l}_1+l_2\\ l_1\end{array}\right){\stackrel{\&overbar;}{s}}_{l_1,l_2,l}^{{\mathrm{\κ}}_1,{\mathrm{\κ}}_2}({\mathrm{\ω}}_a,{\mathrm{\ω}}_b,{\hat{k}}_n)k_n^{l}f({\underset{~}{k}}_n,{\underset{~}{r}}_{ab},l);$

According to many body polarizations energy in following formula counting lattices:

$e_{ind}^{lattice}\left(damped\right)=e_{ind}^{ewald}\left(undamped\right)+{\mathrm{\δe}}_{damp}^{cluster}.$

8. a kind of method of drug crystal forms prediction, comprising:

According to the method as any one of claim 1-7, build the Polarized force field model with regard to drug molecule；

According to constructed Polarized force field model, calculate the crystal formation of this drug molecule.

9. a kind of application building chemical molecular Polarized force field, is suitable to be distributed in one or more computing devices, this application bag Include:

Optimization calculation unit, based on ab initio calculation (ab initio), is suitable to the quantum chemistry structure of chemical molecular is carried out excellent Change and calculate, to obtain the molecular structure of local optimum；

Ionic potential computing unit, is suitable to calculate the first energy in neutral charge for the molecular structure of this local optimum and is carrying in advance The difference of the second energy during fixed positive valence charge, and using this difference as this molecular structure corresponding vertical ion gesture；

Force field parameter computing unit, based on the molecular structure of this vertical ion gesture and this local optimum, is suitable to calculate this molecule knot The Polarized force field parameter of structure, this Polarized force field parameter include atom distribution multistage vector, atom distribution multistage polarizability and Frequency dependence polarizability；And

Field of force construction unit, is suitable to the molecular structure based on local optimum and calculated Polarized force field parameter, builds corresponding Polarized force field model.

10. a kind of system of drug crystal forms prediction, comprising:

According to the application building chemical molecular Polarized force field as claimed in claim 9；And

Crystal formation prediction application, is suitable to, according to constructed Polarized force field model, calculate the crystal formation of this drug molecule.