Molecular Mechanics (PDF)

Lec-5 Molecular Mechanics Dr. Uzma Habib We don’t give a damn where the electrons are. Words to the author, from the president of a well-known chemical company, emphasizing his firm’s position on basic research. “Encyclope...

Lec-5 Molecular Mechanics Dr. Uzma Habib We don’t give a damn where the electrons are. Words to the author, from the president of a well-known chemical company, emphasizing his firm’s position on basic research. “Encyclopedia of Computational Chemistry”, 5th volume, John Wiley and Sons, Inc. “Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics”, Errol Lewars, Kluwer Academic Publishers, 2004 “Introduction to Computational Chemistry”, 2nd Edition, Jensen F, John Wiley and Sons, Inc, 2007.  Molecular Modeling Molecular modeling in the broadest sense is the use of: a) Physical representations: i.e. Plastic Molecular Models b) Graphical representations: Space Filling Model Ball and Stick Model c) Mathematical representations Electrostatic Potential Map, Molecular Orbital Representation Solve Quantum Mechanical equations to Solve Quantum Mechanical determine molecular orbital equations to determine electron appearance position and atomic charge Molecular Mechanics Background The "mechanical" molecular model was developed out of a need to describe molecular structures and properties in as practical a manner as possible. ･The great computational speed of molecular mechanics allows for its use in molecules containing thousands of atoms. Molecular mechanics methods are based on the following principles: ･Nuclei and electrons are lumped into atom-like particles. ･Atom-like particles are spherical (radii obtained from measurements or theory) and have a net charge (obtained from theory). ･Interactions are based on springs and classical potentials. ･Interactions must be pre-assigned to specific sets of atoms. ･Interactions determine the spatial distribution of atom-like particles and their energies. Perspective  Molecular mechanics (MM) is based on a mathematical model of a molecule as a collection of balls (corresponding to the atoms) held together by springs (corresponding to the bonds) 4 Perspective  Within the framework of this model, the energy of the molecule changes with geometry because the springs resist being stretched or bent away from some “natural” length or angle, and the balls resist being pushed too closely together.  The mathematical model is thus conceptually very close to the intuitive feel for molecular energetics that one obtains when manipulating molecular models of plastic or metal: ◦ the model resists distortions (it may break!) from the “natural” geometry that corresponds to the bond lengths and angles imposed by the manufacturer, and in the case of space-filling models, atoms cannot be forced too closely together.  The MM model clearly ignores electrons. 5 Perspective The principle behind MM is or more precisely, to express the energy of a to use this energy equation molecule as a function of to find the bond lengths, to the various possible its resistance toward bond angles, and dihedrals potential energy surface stretching, bond bending, corresponding to the minima and atom crowding, minimum-energy geometry In other words, MM uses a conceptually mechanical model of a molecule to find its minimum-energy geometry (for flexible molecules, the geometries of the various conformers). The form of the mathematical expression for the energy, and the parameters in it, constitute a forcefield, and MM methods are sometimes called forcefield methods. 6 Force Field  A force field (FF) is a set of equations describing the potential energy surface of a chemical system. ◦ A molecular mechanics (MM) method uses a force field based on a classical mechanical representation of molecular forces to calculate static properties of a molecule (e.g., structure and energy of an energy minimum structure). ◦ Molecular dynamics (MD) also implements a force field but generates dynamic properties (e.g., evolution of an structure in time) by calculating forces and velocities of atoms.  FFs are typically constructed to yield experimentally accurate structures and relative energies. 7 Force Field  The most basic component in a FF is the atom type and one element usually contributes several atom types.  Each bond is characterized by the atom types involved and has a “natural” bond length since the variation with the chemical environment is relatively small.  Similarly, bond angles between atom types have typical values. 8 Force Field  The energy absorptions in infrared (IR) spectroscopy associated with a certain bond stretch or angle deformation also fall in narrow ranges, which demonstrates that the variation of force constants is also relatively small.  The existence of an increment system for heats of formation, for example, shows that the energy behaves additively as well.  Hence, in MM the energy is expressed classically as a function of geometric parameters. 9 Energy Terms  Advanced force fields distinguish several atom types for each element (depending on hybridization and neighboring atoms) and introduce various energy contributions to the total force field energy, 𝐸𝐹𝐹 : 𝐸𝐹𝐹 = 𝐸𝑠𝑡𝑟 + 𝐸𝑏𝑒𝑛𝑑 + 𝐸𝑡𝑜𝑟𝑠𝑖𝑜𝑛 + 𝐸𝑣𝑑𝑊 + 𝐸𝑒𝑙𝑠𝑡 + 𝐸𝑐𝑟𝑜𝑠𝑠 Where: 𝐸𝑠𝑡𝑟 and 𝐸𝑏𝑒𝑛𝑑 are energy terms due to bond stretching and angle bending, respectively; 𝐸𝑡𝑜𝑟𝑠 depends on torsional angles describing rotation about bonds; 𝐸𝑣𝑑𝑊 and 𝐸𝑒𝑙𝑠𝑡 describe (nonbonded) van der Waals and Figure 2.1 Illustration of the electrostatic interactions, fundamental force field energy respectively terms Finally 𝐸𝑐𝑟𝑜𝑠𝑠 describes coupling between the first three terms Force Fields Energy 𝐸𝐹𝐹 = 𝐸𝑠𝑡𝑟 + 𝐸𝑏𝑒𝑛𝑑 + 𝐸𝑡𝑜𝑟𝑠𝑖𝑜𝑛 + 𝐸𝑣𝑑𝑊 + 𝐸𝑒𝑙𝑠𝑡 + 𝐸𝑐𝑟𝑜𝑠𝑠  Given such an energy function of the nuclear coordinates, geometries and relative energies can be calculated by optimization.  Stable molecules correspond to minima on the potential energy surface, and they can be located by minimizing 𝑬𝑭𝑭 as a function of the nuclear coordinates.  Conformational transitions can be described by locating transition structure on the 𝑬𝑭𝑭 surface. 11 Stretch Energy  𝐸𝑠𝑡𝑟 is the energy function for stretching a bond between two atom types A and B.  The harmonic approximation gives the stretch energy of a bond between atom types A and B, 𝐸𝑠𝑡𝑟 𝐴𝐵 , as 𝐸𝑠𝑡𝑟 𝑅𝐴𝐵 − 𝑅0AB = 𝑘𝐴𝐵 𝑅𝐴𝐵 − 𝑅0𝐴𝐵 2 = 𝑘𝐴𝐵 ∆𝑅𝐴𝐵 2 Where 𝑘𝐴𝐵 is the force constant and ∆𝑅𝐴𝐵 = 𝑅𝐴𝐵 − 𝑅0AB is the bond length deviation from the natural value, 𝑅0AB, for which 𝐸𝑠𝑡𝑟 is defined to be zero 12 Stretch Energy  Further improvement can be achieved by including higher anharmonic terms to the equation.  While these expressions describe the potential well for 𝑅 close to 𝑅0, the energy goes to infinity for large distances (see Fig. 3 in the front).  In contrast, a morse potential allows the energy to approach the dissociation energy, D, as 𝑅 increases: 𝑘ൗ 𝐸𝑀𝑜𝑟𝑠𝑒 ∆𝑅 = 𝐷(1 − 𝑒 2𝐷∆𝑅 2 ) (2.6) but it is much more expensive in terms of computational cost. 13 Bending Energy  𝐸𝑏𝑒𝑛𝑑 is the energy required for bending an angle formed by three atoms A—B—C, where there is a bond between A and B, and between B and C.  Similarly to 𝐸𝑠𝑡𝑟, 𝐸𝑏𝑒𝑛𝑑 is usually expanded as a Taylor series around a “natural” bond angle and terminated at second order, giving the harmonic approximation. 𝐸𝑏𝑒𝑛𝑑 𝜃 𝐴𝐵𝐶 − 𝜃0 𝐴𝐵𝐶 = 𝑘 𝐴𝐵𝐶 𝜃 𝐴𝐵𝐶 − 𝜃0 𝐴𝐵𝐶 2 (2.7) 14 Bending Energy 𝐸𝑏𝑒𝑛𝑑 𝜃 𝐴𝐵𝐶 − 𝜃0𝐴𝐵𝐶 = 𝑘 𝐴𝐵𝐶 𝜃 𝐴𝐵𝐶 − 𝜃0𝐴𝐵𝐶 2 (2.7) While the simple harmonic expansion is adequate for most applications, there may be cases where higher accuracy is required. The next improvement is to include a third order term, analogous to 𝐸𝑠𝑡𝑟. This can give a very good description over a large range of angles, as Figure 2.4 The bending energy for CH4 illustrated in Figure 2.4 for CH4. 15 Bending Energy 𝐸𝑏𝑒𝑛𝑑 𝜃 𝐴𝐵𝐶 − 𝜃0𝐴𝐵𝐶 = 𝑘 𝐴𝐵𝐶 𝜃 𝐴𝐵𝐶 − 𝜃0𝐴𝐵𝐶 2 (2.7)  The simple harmonic approximation (P2) is seen to be accurate to about ±30ӷ from the equilibrium geometry and the cubic approximation (P3) up to ±70.ӷ Higher order terms are often included in order also to reproduce vibrational frequencies Figure 2.4 The bending energy for CH4 16 Out-of-plane Bending Energy  If the central B atom in the angle ABC is sp2- hybridized, there is a significant energy penalty associated with making the centre pyramidal, since the four atoms prefer to be located in a plane.  If the four atoms are exactly in a plane, the sum of the three angles with B as the central atom should be exactly 360°, however, a quite large pyramidalization may be achieved without seriously distorting any of these three angles. Figure 2.6 Out-of-plane variable definitions 17 Out-of-plane Bending Energy  Taking the bond distances to 1.5Å, and moving the central atom 0.2 Å out of the plane, only reduces the angle sum to 354.8° (i.e. only a 1.7° decrease per angle).  The corresponding out-of-plane angle, χ, is 7.7° for this case.  Very large force constants must be used if the ABC, ABD and CBD angle distortions are to reflect the energy cost associated with the pyramidalization. Figure 2.6 Out-of-plane variable definitions 18 Out-of-plane Bending Energy  This would have the consequence that the in-plane angle deformations for a planar structure would become unrealistically stiff.  Thus a special out-of-plane energy bend term (𝐸𝑜𝑜𝑝) is usually added, while the in-plane angles (ABC, ABD and CBD) are treated as in the general case above. Figure 2.6 Out-of-plane variable definitions 19 Out-of-plane Bending Energy  𝐸𝑜𝑜𝑝 may be written as a harmonic term in the angle χ (the equilibrium angle for a planar structure is zero)  or as a quadratic function in the distance 𝑑, as given in eq. (2.8) and shown in Figure 2.6. 𝐸𝑜𝑜𝑝(𝜒) = 𝑘𝐵 𝜒2 or (2.8) 𝐸𝑜𝑜𝑝(𝑑) = 𝑘𝐵𝑑2 Figure 2.6 Out-of-plane variable definitions 20 Out-of-plane Bending Energy 𝐸𝑜𝑜𝑝(𝜒) = 𝑘𝐵 𝜒2 or (2.8) 𝐸𝑜𝑜𝑝(𝑑) = 𝑘𝐵𝑑2  Such energy terms may also be used for increasing the inversion barrier in sp3- hybridized atoms (i.e. an extra energy penalty for being planar), and 𝐸𝑜𝑜𝑝 is also sometimes called 𝐸𝑖𝑛𝑣.  Inversion barriers are in most cases (e.g. in amines, NR3) adequately modelled without an explicit 𝐸𝑖𝑛𝑣 term, the barrier arising naturally from the increase in bond angles upon inversion. For each sp2-hybridized atom there is one additional out-of-plane force constant to be determined, 𝑘𝐵. Figure 2.6 Out-of-plane variable definitions 21 The Torsional Energy  𝐸𝑡𝑜𝑟𝑠 describes part of the energy change associated with rotation around a B—C bond in a four-atom sequence A— B—C—D, where A—B, B—C and C—D are bonded.  Looking at the B—C bond, the torsional angle is defined as the angle formed by the A—B and C—D bonds  The angle 𝜔 may be taken to be in the range [0°,360°] or [−180°,180°]. 22 The Torsional Energy  To encompass the periodicity, 𝐸𝑡𝑜𝑟𝑠 is written as a Fourier series. 𝐸𝑡𝑜𝑟𝑠 𝜔 = σ𝑛=1 𝑉𝑛𝑐𝑜𝑠(𝑛𝜔) (2.9) ◦ The n = 1 term describes a rotation that is periodic by 360°, the n = 2 term is periodic by 180°, the n = 3 term is periodic by 120°, and so on ◦ The 𝑉𝑛 constants determine the size of the barrier for rotation around the B—C bond ◦ Depending on the situation, some of these 𝑉𝑛 constants may be zero  In ethane, for example, the most stable conformation is one where the hydrogens are staggered relative to each other, while the eclipsed conformation represents an energy maximum 23 The Torsional Energy  As the three hydrogens at each end are identical, it is clear that there are three energetically equivalent staggered, and three equivalent eclipsed, conformations.  The rotational energy profile must therefore have three minima and three maxima.  In the Fourier series only those terms that have n = 3, 6, 9, etc., can therefore have 𝑉𝑛 constants different from zero. 24 The Torsional Energy  For rotation around single bonds in substituted systems, other terms may be necessary.  In the butane molecule, for example, there are still three minima, but the two gauche (torsional angle ~ ±60 ͦ)and anti (torsional angle ~180°) conformations now have different energies.  The barriers separating the two gauche and the gauche and anti conformations are also of different height.  This may be introduced by adding a term corresponding to n = 1. 25 The Torsional Energy  For the ethylene molecule, the rotation around the C=C bond must be periodic by 180°, and thus only n = 2, 4, etc., terms can enter.  The energy cost for rotation around a double bond is of course much higher than that for rotation around a single bond in ethane, which would be reflected in a larger value of the 𝑉2 constant. 26 The Torsional Energy  For rotation around the C=C bond in a molecule such as 2-butene, there would again be a large 𝑉2 constant, analogous to ethylene, but in addition there are now two different orientations of the two methyl groups relative to each other, cis and trans. 27 The Torsional Energy The Fourier series allows the representation of potentials with various minima and maxima (Fig. 4). Three terms are enough to model the most common torsional potentials. 28 The Torsional Energy  The torsional energy is fundamentally different from 𝐸𝑠𝑡𝑟 and 𝐸𝑏𝑒𝑛𝑑 in three aspects: ◦ A rotational barrier has contributions from both the non-bonded (van der Waals and electrostatic) terms, as well as the torsional energy, and the torsional parameters are therefore intimately coupled to the non-bonded parameters. ◦ The torsional energy function must be periodic in the angle 𝜔: if the bond is rotated 360° the energy should return to the same value. ◦ The cost in energy for distorting a molecule by rotation around a bond is often low, i.e. large deviations from the minimum energy structure may occur, and a Taylor expansion in 𝜔 is therefore not a good idea. 29 The van der Waals Energy  𝐸𝑣𝑑𝑊 is the van der Waals energy describing the repulsion or attraction between atoms that are not directly bonded.  Together with the electrostatic term 𝐸𝑒𝑙 , it describes the non-bonded energy.  𝐸𝑣𝑑𝑊 may be interpreted as the non-polar part of the interaction not related to electrostatic energy due to (atomic) charges.  This may for example be the interaction between two methane molecules, or two methyl groups at different ends of the same molecule. 30 The van der Waals Energy  𝐸𝑣𝑑𝑊 is zero at large interatomic distances and becomes very repulsive for short distances.  In quantum mechanical terms, the latter is due to the overlap of the electron clouds of the two atoms, as the negatively charged electrons repel each other.  At intermediate distances, however, there is a slight attraction between two such electron clouds from induced dipole–dipole interactions, physically due to electron correlation. 31 The van der Waals Energy  Even if the molecule (or part of a molecule) has no permanent dipole moment, the motion of the electrons will create a slightly uneven distribution at a given time.  This dipole moment will induce a charge polarization in the neighbor molecule (or another part of the same molecule), creating an attraction, and it can be derived theoretically that this attraction varies as the inverse sixth power of the distance between the two fragments. 32 The van der Waals Energy  The induced dipole–dipole interaction is only the leading term of such induced multipole interactions: ◦ There are also contributions from induced dipole–quadrupole, quadrupole–quadrupole, etc., interactions. ◦ These vary as 𝑅-8, 𝑅-10, etc., and the 𝑅-6 dependence is only the asymptotic behaviour at long distances.  The force associated with this potential is often referred to as a “dispersion” or “London” force.  The van der Waals term is the only interaction between rare gas atoms and it is the main interaction between non-polar molecules such as alkanes. 33 The van der Waals Energy  𝐸𝑣𝑑𝑊 is very positive at small distances, has a minimum that is slightly negative at a distance corresponding to the two atoms just “touching” each other, and approaches zero as the distance becomes large.  A general functional form that fits these conditions is given in eq. (2.11). 𝐴𝐵 𝐶 𝐸𝑣𝑑𝑊 𝑅𝐴𝐵 = 𝐸𝑟𝑒𝑝𝑢𝑙𝑠𝑖𝑜𝑛 𝑅𝐴𝐵 − 𝐴𝐵 6 (2.11) 𝑅  It is not possible to derive theoretically the functional form of the repulsive interaction, ◦ it is only required that it goes toward zero as 𝑅 goes to infinity ◦ it should approach zero faster than the 𝑅-6, term, as the energy should go towards zero from below. 34 The van der Waals Energy  A popular function that obeys these general requirements is the Lennard-Jones (LJ) potential,  where the repulsive part is given by an 𝑅 -12 dependence (𝐶1 and 𝐶2 are suitable constants). 𝐶1 𝐶2 𝐸LJ 𝑅𝐴𝐵 = 𝐴𝐵 12 − 𝐴𝐵 6 (2.12) 𝑅 𝑅  The Lennard-Jones potential can also be written as in eq. (2.13) 𝑅0 12 𝑅0 6 𝐸LJ 𝑅 = 𝜀 −2 (2.13) 𝑅 𝑅 Here 𝑹𝟎 is the minimum energy distance and 𝜺 the depth of the minimum. 35 The van der Waals Energy  The Merck Molecular Force Field (MMFF) uses a generalized Lennard-Jones potential where the exponents and two empirical constants are derived from experimental data for rare gas atoms.  The resulting buffered 14-7 potential is shown in eq. (2.14). 1.07𝑅0 7 1.12𝑅07 𝐸buf14 − 7 𝑅 = 𝜀 𝑅+0.07𝑅0 𝑅7+0.07𝑅07 −2 (2.14) 36 The van der Waals Energy  From electronic structure theory it is known that the repulsion is due to overlap of the electronic wave functions,  and furthermore that the electron density falls off approximately exponentially with the distance from the nucleus (the exact wave function for the hydrogen atom is an exponential function).  There is therefore some justification for choosing the repulsive part as an exponential function.  The general form of the “Exponential – R−6” 𝐸𝑣𝑑𝑊 function, also known as a “Buckingham” or “Hill” type potential, is given in eq. (2.15). 𝐶 𝐸H𝑖𝑙𝑙 𝑅 = 𝐴𝑒 −𝐵𝑅 − (2.15) 𝑅6 𝑨, 𝑩 𝒂𝒏𝒅 𝑪 𝒂𝒓𝒆 𝒉𝒆𝒓𝒆 𝒔𝒖𝒊𝒕𝒂𝒃𝒍𝒆 𝒄𝒐𝒏𝒔𝒕𝒂𝒏𝒕𝒔 37 The van der Waals Energy  A third functional form, which has an exponential dependence and the correct general shape, is the Morse potential, eq. (2.5). 2 𝐸𝑀𝑜𝑟𝑠𝑒 ∆𝑅 = 𝐷 1 − 𝑒 −𝛼∆𝑅 (2.5) Where, 𝐷 is the well depth and 𝛼 is the well width. Also, 𝛼 = 𝑘 2𝐷  It does not have the 𝑅-6 dependence at long range, but as mentioned before, in reality there are also 𝑅−8, 𝑅−10, etc., terms.  The 𝐷 and 𝛼 parameters of a Morse function describing 𝐸𝑣𝑑𝑊 will of course be much smaller than for 𝐸𝑠𝑡𝑟, and 𝑅0 will be longer. 38 The van der Waals Energy  For small systems, where accurate interaction energy profiles are available, it has been shown that the Morse function actually gives a slightly better description than a Buckingham potential, which again performs significantly better than a Lennard-Jones potential.  This is illustrated for the H2—He interaction in Figure 2.10, where the Buckingham and Morse parameters have been derived from the minimum energy and– distance (𝜀 and 𝑅0) and by matching the force constant at the minimum. Figure 2.10 Comparison of 𝐸𝑣𝑑𝑤 functionals for the H2—He potential39 The van der Waals Energy  The main difference between the three functions is in the repulsive part at short distances, the Lennard-Jones potential is much too hard, and the Buckingham also tends to overestimate the repulsion.  Furthermore, it has the problem of “inverting” at short distances.  For chemical purposes, however, these “problems” are irrelevant, since energies in excess of 400 kJ/mol are sufficient to break most bonds and will never be encountered in actual calculations.  The behaviour in the attractive part of the potential, which is essential for intermolecular interactions, is very similar for the three functions, as shown in Figure 2.11. 40 The van der Waals Energy Figure 2.11 Comparison of 𝐸𝑣𝑑𝑤 functionals for the attractive part of the H2—He potential 41 The van der Waals Energy  Part of the better description for the Morse and Buckingham potentials is due to the fact that they have three parameters, while the Lennard-Jones only employs two.  Since the equilibrium distance and the well depth fix two constants, there is no additional flexibility in the Lennard- Jones function to fit the form of the repulsive interaction.  The van der Waals distance, 𝑅0 𝐴𝐵, and softness parameters, 𝜀𝐴𝐵, depend on both atom types 𝐴 and 𝐵.  These parameters are in all force fields written in terms of parameters for the individual atom types.  There are several ways of combining atomic parameters to di-atomic parameters, some of them being quite complicated. 42 The van der Waals Energy  A commonly used method is to take the van der Waals minimum distance as the sum of two van der Waals radii, 𝑅0 𝐴𝐵 and the interaction parameter as the geometrical mean of the atomic “softness” constants. 𝑅0 𝐴𝐵 = 𝑅0 𝐴 + 𝑅0 𝐵 (2.18) 𝜀𝐴𝐵 = 𝜀𝐴𝜀𝐵 43 The Electronic Energy: Charges and Dipoles  The electrostatic energy, 𝐸𝑒𝑙𝑠𝑡 , is due to the electrostatic interactions arising from polarized electron distributions based on electronegativity differences.  It can be modeled by Coulomb interactions of point charges associated with individual atoms: 𝑄𝐴𝑄𝐵 𝐸𝑒𝑙𝑠𝑡 𝑅𝐴𝐵 = 𝜀𝑅𝐴𝐵  𝜀 being a dielectric constant, which can be used to model the effect of the same or other molecules present (e.g., solvent).  The atomic charges, Q, are commonly obtained by fitting to the electrostatic potential as calculated by an electronic structure method.  An 𝐸𝑒𝑙𝑠𝑡 description based on dipole–dipole interactions between polarized bonds can alternatively be employed. 44 The Electronic Energy: Charges and Dipoles  Hydrogen bonds are non-bonded interactions between a positively charged hydrogen atom and an electronegative atom with lone electron pairs (mostly oxygen or nitrogen) and can be adequately modeled by appropriately chosen atomic charges.  Although a single hydrogen bond is a very weak interaction, the large number occurring in biomolecules (e.g., proteins) makes hydrogen bonding a very important factor.  In the large size limit, the bonded interactions increase linearly with the system size, but the non-bonded interactions show a quadratic dependence and determine the computational cost. 45 The Electronic Energy  The van der Waals interactions quickly fall off with the distance (𝑅 −6 dependence) and may be neglected for large separations.  The electrostatic interaction (proportional to 𝑅 −1) is much more far reaching and needs to be considered out to very long distances.  Fast multipole methods (FMMs) can be applied to reduce the computational cost of evaluating 𝐸𝑒𝑙𝑠𝑡. 46 Other Energy Contributions  To improve performance, force fields include further parameters to take care of special cases.  Cross terms account for the interplay between different contributions (e.g., longer bonds for small angles).  Correction terms may be introduced to describe substituent effects (e.g., anomeric effect, a stereoelectronic effect that describes the tendency of heteroatomic substituents adjacent to a heteroatom within a cyclohexane ring to prefer the axial orientation instead of the less hindered equatorial orientation that would be expected from steric considerations).  Additional terms may be introduced to adequately treat special cases like pyramidalization of sp2 hybridized atoms. 47 Other Energy Contributions  Hydrogen bonding may be treated explicitly (in addition to the electrostatic interaction) with a special set of van der Waals interaction parameters.  Pseudo atoms maybe introduced to model lone pairs.  In addition, atoms in unusual bonding situations (three-membered rings, molecules with linearly conjugated 𝜋-systems, aromatic compounds, etc.), which are not described adequately by the normal parameters, can be defined as new atom types. 48 Other Energy Contributions  The force field energy, 𝐸𝐹𝐹 , corresponds to the energy relative to a molecule with non-interacting fragments.  Therefore, only energies for molecular structures built from the same fragments (conformers) can be compared directly.  So that energy between different molecules (isomers) can be compared.  The energy scale is converted to heats of formation by adding bond increments (estimated from bond dissociation energies minus the heat of formations of the atoms involved) and possibly group increments (e.g., methyl group): 𝑏𝑜𝑛𝑑𝑠 𝑔𝑟𝑜𝑢𝑝𝑠 ∆𝐻𝑓 = 𝐸𝐹𝐹 + ෍ ∆𝐻𝐴𝐵 + ෍ ∆𝐻𝐺 49 2 - Parametrization  Having settled on the functional description and a suitable number of cross terms, the problem of assigning numerical values to the parameters arises.  Determining the parameters for a force field is a substantial task.  In general, not all necessary data are available from (accurate) experiments.  Modern electronic structure computations can provide unknown data relatively easily and with sufficient accuracy.  Another problem is the large number of parameters: for a force field with N atom types, the number of different types of bonds, bond angles, and dihedral angles scales as N2, N3, and N4, respectively, each requiring several parameters. 50 2 - Parametrization  So that the number of parameters can be reduced, the atom dependency can be reduced (e.g., the torsional parameters may be treated as dependent on the B–C central bond only and not on the atom types A and D).  The parametrization effort can be reduced further by defining “generic” parameters to be used for less common bond types or when no reference data are available.  This, of course, reduces the quality of a calculation.  By deriving the di-, tri-, and tetra-atomic parameters from atomic data (atom radii, electronegativities, etc.), universal force fields (UFFs) allow one to include basically all elements.  The performance, however, is relatively poor. 51 2 - Parametrization  The kind of energy terms, their functional form, and how carefully (number, quality, and kind of reference data) the parameters were derived determine the quality of a force field.  Accurate force fields exist for organic molecules (e.g., MM2, MM3),  But more approximate force fields (e.g., with fixed bond distances) optimized for computational speed rather than accuracy are the only practical choice for the treatment of large biomolecules ◦ AMBER (assisted model building with energy refinement), ◦ CHARMM (chemistry at Harvard molecular mechanics), ◦ GROMOS (Groningen molecular simulation)  The type of molecular system to be studied determines the choice of the force field. 52 2 - Parametrization  One limitation of force field methods is that they can describe only well-known effects that have been observed for a large number of molecules (this is necessary for the parametrization).  The predictive power of these methods is limited to extrapolation or interpolation of known effects. 53 54