US 7797144 B2
A molecular modeling method and system to provide a dielectric function that reflects the plurality of models simultaneously used in molecular modeling. For interactions similar to those used during the development of a charge model, the dielectric is similar to the charge model dielectric. In a second embodiment, an approximation for electrostatic/electrodynamics modeling is provided that analyzes charge distribution and determines groups of charges that can be treated as a whole at moderate distances.
1. A method to molecularly model electrostatic interactions in a homogeneous portion of a chemical system with a plurality of charges from a charge model, each said charge at a defined position, said method comprising:
providing a dielectric function that is selected from a plurality of dielectric functions found in a plurality of models; and
calculating, by use of a data processing system specifically programmed with software to implement said dielectric function, an electrostatic energy between a plurality of sets of two charges in said chemical system to reduce a development time required to develop at least one new material with one or more desirable properties, wherein a plurality of distinct models apply.
2. The method of
3. The method of
4. The method of
5. The method of
6. The method of
7. The method of
8. The method of
9. The method of
10. The method of
11. The method of
12. The method of
13. The method of
14. The method of
15. The method of
16. The method of
1. Field of the Invention
This invention relates generally to providing a method and system to facilitate molecular modeling, and more specifically to providing a method and system to approximate the electrostatics and electrodynamics functions used in molecular modeling software.
2. Description of the Prior Art
Significant advances have been made in predicting, studying, and quantifying the nature of interatomic and intermolecular interactions with the use of computer models. Although such computer models can be useful in testing and validating a theory concerning the nature of these interactions, computer models find especially useful application in reducing the time required to develop new materials with desirable properties. Materials that may be developed with additional efficiency through the use of computer models include polymers, pesticides, herbicides, pharmaceuticals, semiconductor materials for integrated circuits, and petrochemicals.
Several modeling techniques are used to study chemical systems. These techniques are referred to as simulations, models, methods, approximations, and so forth. Here all techniques are referred to as models. Typically, the model selected provides a user of the model with a particular compromise between physical accuracy and the computing resources required to run the model. Quantum mechanical calculations can be performed with a high degree of accuracy, but are very expensive in terms of the computer time required to perform them. In those fields described above, where new polymers, drugs, and other products are being developed, it is more useful to use modeling techniques that requires less investment in computing resources, so that more candidate materials or molecules can be analyzed in a shorter time period for the properties desired.
A general class of models called force fields is frequently used. A force field uses analytical classical equations to describe molecules, e.g., the potential energy of a diatomic molecule can be modeled with a harmonic equation. Occasionally non-analytic terms are included, e.g., induced polarization terms. These functions are significantly simpler computationally than performing the more correct and mathematically complex quantum mechanical calculations.
Force fields are developed to reproduce experimental data or calculation results, generally quantum mechanics calculations. The selection of appropriate parameters and force field functional dependencies is a significant factor in the success of the model. Furthermore, the number of fitted parameters used in the model relative to the number of measurable or calculable values relevant to the system being modeled should be as low as possible to minimize the possibility of chance correlations.
New force field functions are continually under development to increase agreement with experiment and quantum results, to maximize the number of observables relative to the number of parameters, and to limit the necessity of performing computationally expensive calculations.
The potential energy function is the energy of a set of atoms as a function of atomic coordinates. One class of potential energy functions is the force field function. For any chosen force field function and a given set of atoms, the location of the atoms, the types of atoms, and the bonds between the atoms define a potential energy field. For a given set of atomic coordinates (configuration) the energy is well defined. On rare occasions there is bond breaking and bond making, generally explicitly defined by the user. The potential energy field is usually expressed as a sum of terms, with each term representing a particular class of interaction.
Most force fields include terms that represent bonded interactions such as bond stretching and contraction, angle stretching and contraction, rotation about a dihedral bond, out-of-plane deformation, cross terms of the above and any other term or parameter that the method developer considers to be significant. These interactions are sometimes referred to as intramolecular interactions.
Non-bonded, or intermolecular interaction potential energy generally includes a van der Waals interaction and a model using electrostatic interactions. The van der Waals (vdW) term may take a variety of forms. In general, vdW terms are mildly attractive at large separations and become very repulsive at small atomic separations. Different known parameterizations of the van der Waals term, and intramolecular force parameters, may be used.
A general expression for a potential energy field generally used for molecular modeling may thus be described as follows:
The electrostatic term is calculated using charges from a charge model. The charge model is generally developed for a set of atoms such as a molecule (composed of electrons and nuclei), an amino acid, a nucleic acid, a cell in a regular array, an aggregate, or other group of atoms (henceforth any said set of atoms is called a substructure).
The probability of electron position is provided by the quantum mechanics wavefunction, which is the solution to the Shroedinger Equation. This probability is generally considered to describe a charge distribution. Chemists frequently approximate this charge distribution with a classical description composed of point charges, dipoles, quadrupoles, and higher multipoles. The choice of charges and the methodology used to determine the charges is called a charge model.
Different techniques have been used to develop charge models. For example, the method pioneered by Kollmann et al. is disclosed in “Force Fields for Protein Simulations,” by J. W. Ponder and D. A. Case in Adv. Prot. Chem., vol. 66, pp. 27-85 (2003); and in “Molecular Dynamics Simulation of Nucleic Acids: Successes, Limitations, and Promise,” by T. E. Cheatham, III and M. A. Young in Biopolymers, vol. 56, pp. 232-256 (2001); which are hereby incorporated by reference. This method uses point charges placed at atom centers to reproduce the electrostatic field around the molecule as calculated by placing test charges in the molecule's periphery in a quantum mechanics calculation. In other methods the charge model reproduces other properties, such as hydrogen bond strength. In addition there are programs that assign charges based on the bonding structure of the molecule. For the purposes of this patent an atomistic charge model is a charge model that is any combination of the following: (1) assigns fractional electron charges to the vast majority of atoms in a chemical system, (2) assigns dipoles to the vast majority of the atoms in a chemical system, (3) assigns charges between the vast majority of atoms that are bonded, and (4) assigns dipoles between the vast majority of atoms that are bonded. Molecular modeling frequently uses this charge model to calculate the electrostatic potential between molecules and between different parts of the same molecule.
Molecular modeling simultaneously uses multiple independent models where the dielectric terms do not necessarily correspond with each other. The dielectric of a material is a macroscopic property. Poisson's equation describing electrostatics was developed to describe electrostatic interactions on a macroscopic scale, and the dielectric concept was incorporated into the equations to accommodate the dielectric properties of matter. However, there is no term for a dielectric in the Schroedinger Equation, which applies to a molecular scale. In the work of Kollman et al. noted above, the dielectric was not included in the classical equations. This has the effect of setting the classical dielectric to one, i.e., D=1. Alternately, a dielectric of 2 or 4 is frequently used for protein calculations.
There are various approaches to handle a dielectric when using classical electrostatic (macroscopic) equations to model molecules. The dielectric properties of the molecule are implicitly or explicitly incorporated into the model. When developing a charge model, the chosen dielectric is incorporated into the choice of partial charge, i.e., a charge is chosen which produces the desired properties given the dielectric used. In the charge model of Kollman et al., the charges create an electrostatic field about the molecule that approximates the quantum result given a dielectric of one. When there are two molecules, Coulomb's equation describes the interaction energy between them.
As an example of the prior art approach, consider the electrostatic potential between two charges qi and qj, at points i and j respectively, shown in
Composite Charge Model
In the prior art, algorithms for calculation of intermolecular electrostatic potential were modeled by detailed calculation of charge-charge interactions.
“The potential due to an arbitrary charge distribution at a point outside the distribution such that r>r1 is the same as: (a) the potential of a point charge equal to the net charge of the distribution, plus (b) that of a point dipole with a dipole moment equal to the dipole moment of the distribution, plus (c) that of a point quadrupole with a quadrupole moment equal to that of the distribution, and so on, all located at the origin. where r1 is the maximal distance between the origin and a charge (distribution). The potential due to a dipole moment varies as 1/r2 and the potential due to a quadrupole moment varies as 1/r3.” The impact of this result is evident in the charges found in charge models.
The point charge model has produced individual atom-atom interactions that are large. For example, consider the potential calculated using the Amber 1998 [reference] parameter set. The potential between the two aliphatic terminal carbons in two amino acid Leucine molecules that are 45 Angstroms (A) apart is 10 kcal/mole. (i.e., QiQj/rij=300(0.3)(0.3)/45=10 kcal), which is large. From an experimental perspective, the electrostatic potential between the side chains in Leucine molecules that are 45A apart is negligible. The apparent contradiction is resolved by the composite charge result noted above; in a calculation the large 10 kcal/mole term is counterbalanced by the similar terms of opposite sign in neighboring atoms and the result is small.
The large number of charges generated by the charge model introduces further problems. There may be very large number of electrostatic interactions. For example, in a protein with 3000 atoms, there are ˜n2/2 interactions or ˜4,500,000 terms. In many applications, the electrostatic potential is sampled repeatedly and the computer time required becomes quite large. Significant effort has been expended to decrease the computer time. Three techniques are noted below.
Cutoff distances have been explored where all interactions from charges greater than an “R” distance apart are set to zero. “R” is called the cutoff distance. This produces a discontinuous function at the cutoff distance, and various means have been developed to smooth this function. While smoothing makes the function differentially continuous, the cutoff approach introduces large error into the calculations.
In U.S. Pat. No. 5,612,894, which is hereby incorporated by reference, Wertz discloses an alternative procedure to reduce the number of electrostatic calculations. This procedure determines the sensitivity (importance) of a desired quantity to a given interaction and restricts the calculation of interactions with low sensitivity. This method requires techniques to recalculate whether a given interaction is important and techniques to “remember” the important ones. The interactions are separated into classes of varying importance. In some cases terms that could be in the same category to provide the counterbalancing charge terms noted above will be in different classes.
In U.S. Pat. No. 5,915,230, which is hereby incorporated by reference, Berne et al. disclose an approximation for electrostatic potential at a distance from a set of charges. A three-dimensional space is divided into cubic volumes (boxes) and the charges inside a box are handled as a set. Unfortunately the division into boxes is artificial and counterbalancing charges can be in different boxes. Charge models depend on the canceling effects at large distances as noted above. Separating the opposite, canceling charges into different boxes creates computational artifacts in the functions and makes them less stable.
As can be seen from the preceding discussion, the prior art methods for approximating potential energy functions are frequently not practical for complex molecular modeling involving many atoms. What is needed is a more accurate approximation of potential energy functions for molecular modeling, which will provide sufficient accuracy for the modeling of even complex groups of atoms.
The present invention provides a method and system for molecular modeling. Embodiments of the invention can be implemented in numerous ways, such as a method, a programmable logic device, or a program stored on electronically-readable media. Four aspects of the invention are described below.
A first aspect of the invention is directed to a method to molecularly model a chemical system. The method includes determining a plurality of charges at a plurality of charge positions in a molecular system to be modeled; providing a dielectric function that reflects a plurality of dielectric functions found in a plurality of models, wherein for at least one electrostatic/electrodynamics calculation between two or more charges (or charge distributions) having a separation distance, the separation distance is divided into regions where a plurality of different models apply; and calculating electrostatic energy between the two charges (or charge distributions).
A second aspect of the invention is again directed to a method to model the potential energy function that models a chemical system, where the invention is used for at least one calculation of the electrostatic/electrodynamics energy between two charges or charge distributions. The method includes establishing a first portion of the electrostatics/electrodynamics equation where a charge model uses a dielectric function and assigning the dielectric function to the first portion; assigning a chemically appropriate dielectric function for a second portion; and modeling an electrodynamics equation that provides acceptable accuracy for a neighborhood of the atoms associated with the first portion, wherein the modeling accommodates a switch from a first dielectric of the first portion to a second dielectric of the second portion.
A third aspect of the invention is directed to a program embodied on electronically-readable media, containing instructions to model at least one portion of a molecular system. The program includes a program code segment that establishes a first portion of the electrostatics/electrodynamics equation where a charge model assigns a dielectric function to the first portion; a program code segment to assign a chemically appropriate dielectric function for a second portion of the electrostatics/electrodynamics equation; and a program code segment to model an electrodynamics equation that provides acceptable accuracy for a neighborhood of the atoms used to develop the charge model, wherein the modeling accommodates a switch from a first dielectric of the first portion to a second dielectric of the second portion of the electrodynamics/electrostatics equation.
A fourth aspect of the invention is again directed to a method to model the potential energy function that models a chemical system, where groups of charges are chosen that have a combined charge that is substantially smaller than the individual charges and the combined charge is used, instead of the individual charges, as dictated by the details of the calculation and computer, programmer, and user resources.
A model is provided to incorporate different assumptions regarding the dielectric capacity of materials for use in molecular modeling calculations. The charge distribution in molecules is frequently modeled as a function of point charges, dipoles and/or higher multipoles. For a test charge in the periphery of the molecule, the electrostatic potential is calculated with the Poisson Equation using the charge model dielectric. When the test charge is a distance from the molecule, there may be dielectric material between the molecule and the test charge and the electrostatic potential is modified by the presence of this dielectric. The invention accommodates the need for a dielectric term that changes according to the circumstances. The dielectric function is implemented in software on a data processing system (e.g., a general-purpose computer, a personal computer, a hand-held device, or an equivalent).
The electrostatic field outside a group of charges can be modeled as a function of the net charge combined with a dipole term and higher order terms. The choice of which charges to group together effects the stability of the approximation. In embodiments reviewed below several groupings of charges are provided which take advantage of features of chemistry and charge models to increase function stability and to reduce computer time required. The resulting approximation in one embodiment is implemented in software on a data processing system.
One embodiment of the invention provides a means to model a dielectric function that takes into account the assumptions made during the development of the charge model, the ‘screening’ effect that atoms have on electrostatic potential on a molecular scale, and the macroscopic concept of a dielectric. For example, in the charge model of Kollman et al., the charges reproduce the calculated electrostatic field about the molecule for a dielectric of one. When there are two molecules, Coulomb's equation describes the interaction energy between them, again with a dielectric of one. When there is a third molecule between the two molecules, the third molecule has dielectric properties larger than one and the dielectric effect of the third molecule should be incorporated into the equation for the electrostatic potential. A model is needed that incorporates the dielectric effect of the third molecule while maintaining the equation used for two molecules with vacuum between.
The electrostatic calculation for charges within a group are determined by the details of charge model development. Electrostatic calculation involving charges in different group A of atoms 310 and group B of atoms 318 with no third group C of atoms 312, or dielectric D 314, between them is calculated using the standard equation for the chosen charge model and/or application.
For the case where there are molecules, a third group C of atoms 312, and/or dielectric D 314, between charges, there are several options. The first option is to use a standard macroscopic equation for vacuum and for the volume of the molecules/dielectric. The second option is to use a vacuum equation for the nearest neighbor atoms and dielectric for atoms further away for a chosen definition for nearest neighbor atoms. The third option is to model the dielectric with a function that is one for distances smaller than the van der Waals distance or other specific distance, where the function varies linearly from 1 at van der Waals distance to a desired bulk dielectric at a distance determined to give sufficient accuracy (for example, the van der Waals distance plus 2.0 Angstroms), and is the bulk dielectric value at distances greater than the specified distance as depicted in
Composite Charge Model
Several objects and advantages of the composite charge embodiment of invention are to (1) provide an alternate function to estimate electrostatic interaction, (2) reduce the number of electrostatic terms calculated, (3) give electrostatic approximations that are better behaved, and (4) aid the analysis of calculation results. Further objects and advantages of the invention will become apparent from a consideration of the drawings and ensuing description.
Consider the electrostatic potential between a test charge, q, and a methyl group, —CH3, where the test charge is a distance from the methyl group. The electrostatic potential between q and —CH3 using a point charge model for —CH3 is:
The electrostatic potential approximation is applied, which is described by Lorrain and Corson, as reproduced in the prior art section, where the carbon atom and three attached hydrogens are a set of charges. The charge on each hydrogen is added to the charge of the attached carbon and the point of reference is the position of the carbon nucleus. The dipole and higher multipoles are calculated for the carbon and its attached hydrogens, again the position of the carbon nucleus is the reference point. Then the potential can be modeled as
The total charge on the methyl group is near zero and qtotal˜0. The interaction between a point charge and a multipole is small since the terms drop off rapidly with distance, and the second and higher order terms are also small. The first term in equation (4) is the only significant term and the error in the composite charge model is linear in 1/r. The determination of the distance needed for the desired accuracy depends on the dielectric constant (D) of the material between the two charges and the demands of the calculation.
In this example, the invention uses one term. Compare this with the approximation formerly used in equation (3), which contains 4 terms. This substantial savings is common due to the character of computer generated charge models. The charge models generally use oppositely signed charges in close proximity to each other to model quantum mechanics effects. The charges are classical mimics of quantum results in the immediate vicinity of a molecule. The fields from these opposing charges rapidly drop off and can be approximated by the sum of the charges at moderate distances, where a moderate distance is the second or third sphere of nearest neighbors. As a result, the first term in equation (4) tends to be small and is frequently zero. In either case it can be set to zero for sufficiently large r. The remaining terms, the dipole term and higher multipole terms, tend to decrease rapidly with distance and can frequently be set to zero. However, the multipole terms occasionally play an important role in aggregate and need to be considered carefully. For example, molecules are expected to orient about an ion and the aggregate dipole becomes significant. The aggregate multipole terms could be sampled every n iterations and then handled as appropriate. However, in general, the net electrostatic potential at a distance from a molecule tends to be small as compared to the electrostatic potential due to the partial charge on a single atom.
In one embodiment of the invention, atoms are grouped by molecule, side-chain, or other division. For every group, a reference position point (x, y ,z) is chosen that is within the region of the group and the charges associated with the atoms are analyzed to produce total charge, dipole and higher multipoles centered on the reference point. Frequently dipole and/or higher multipoles can be ignored. This description of the charge distribution is then used for the electrostatic calculation for all points greater than r Angstroms away from the reference point, where r is chosen based on the dielectric and the desired accuracy. In this approximation, r must be greater than the greatest distance between the reference point and all atoms. Examples of groups are H2O, the backbone atoms in an amino acid, side-chains of amino acids, an ion and its counter ion, all neutral molecules, and all atoms in a molecule except for a subgroup (e.g., in nitropentane, all atoms except —NO2).
In one preferred embodiment, every heavy atom and all the attached hydrogen atoms are considered to be a group; the charges on hydrogen atoms are added to the charge on the attached heavy atom. A convenient reference point is the heavy atom center. In an electrostatic calculation, the complete equation (3) is used when the heavy atom is less than r Angstroms away. When the distance is greater than r then the individual hydrogen and heavy atom charges are ignored and the summed charge is used. In this embodiment, only the distances to the reference point (e.g., the heavy atoms) need be calculated when the distance is greater than r. In another embodiment, this sum is done for just the carbon heavy atoms, and hetero-atoms are handled using the classical equation (3). Other embodiments include summing all charges in a side chain of an amino acid or other polymer, summing all charges in the monomer backbone, and summing all atoms in a side chain or molecule exclusive of atoms with a formal charge.
In one embodiment, the invention is used to aid analysis of calculations. For example, the charge on hydrogen is added to the charge on the heavy atom it is attached to and the approximation is applied as described above. This step greatly reduces the number of terms calculated, easing the analysis of the remaining terms. In particular if the approximation is applied for all r, this introduces a significant error for small r, but aids analysis of the gross features of the electrostatic field.
One embodiment of the invention is used during the initial stages of a calculation, such as a dynamics calculation, and the complete equation is used at the end and/or intermittently to assure that the end coordinates are a minimum for the complete equation. One embodiment of the invention is used in a quantum mechanics calculation where a group at a sufficient distance away is treated as a point charge, dipole, and multipoles. For example, a water molecule at a distance could be treated as a collection of multipoles. In one embodiment, the invention is used in combination with the invention disclosed by U.S. Pat. No. 5,612,894 or U.S. Pat. No. 5,915,230.
Other embodiments of the invention are possible. For example, the invention may be used in quantum mechanical (QM)/molecular modeling (MM) calculations. The model for the dielectric that accommodates the switch from a charge model dielectric to bulk dielectric may be used for calculations across the QM-MM border. The composite charge methods could be used in a generalization of QM/MM methods where all interactions at a large distance are handled with MM. The invention may be used in conjunction with other models, such as generalized Born calculations, Poisson-Boltzmann calculations, and other condensed phase methods. The invention may also be used with other charge models such as a charge model based on Artificial Intelligence techniques (e.g. a knowledge based system), a charge model based on the wavefunction developed in a quantum mechanics equation using Mulliken charges, or on a model derived from Mulliken charges.
The exemplary embodiments described herein are for purposes of illustration and are not intended to be limiting. Therefore, those skilled in the art will recognize that other embodiments could be practiced without departing from the scope and spirit of the claims set forth below.