RTO care

You are here: Computation Concepts > Molecular Mechanics Theory > Methods Available in CS MOPAC > MM2

Molecular Mechanics Theory

Molecular mechanics computes the energy of a molecule in terms of a set of classical potential energy functions. The potential energy functions and the parameters used for their evaluation are known as a "force-field".

Molecular mechanical methods are based on the following principles:

•  Nuclei and electrons are lumped together and treated as unified particles (atoms).

•  Atoms are typically treated as spheres.

•  Bonds are typically treated as springs.

•  Non-bonded interactions between atoms are described using potential functions derived from classical mechanics.

•  Individual potential functions are used to describe the different interactions: bond stretching, angle bending, torsion (bond twisting), and through-space (non-bonded) interactions.

•  Potential energy functions rely on empirically derived parameters (force constants, equilibrium values) that describe the interactions between sets of atoms.

•  The sum of the interactions determines the conformation of the molecule.

•  Molecular mechanical energies have no meaning as absolute quantities. They can only be used to compare relative steric energy (strain) between two or more conformations of the same molecule.

The Force-Field

Since molecular mechanics treats bonds as springs, the mathematics of spring deformation (Hooke's Law) is used to describe the ability of bonds to stretch, bend, and twist. Non-bonded atoms (greater than two bonds apart) interact through van der Waals attraction, steric repulsion, and electrostatic attraction and repulsion. These properties are easiest to describe mathematically when atoms are considered as spheres of characteristic radii.

The total potential energy, E, of a molecule can be described by the following summation of interactions:

E= Stretching Energy + Bending Energy + Torsion Energy + Non-bonded Interaction Energy

The first three terms are the so-called bonded interactions. In general, these bonded interactions can be viewed as a strain energy imposed by a model moving from some ideal zero strain conformation. The last term, which represents the non-bonded interactions, includes the two interactions shown below.

The total potential energy can be described by the following relationships between atoms. The numbers refer to the atom positions in the figure below.

1. Bond Stretching: (1-2) bond stretching between directly bonded atoms

2. Angle Bending: (1-3) angle bending between two atoms that are adjacent to a third atom.

3. Torsion Energy: (1-4) torsional angle rotation between atoms that form a dihedral angle.

4. Repulsion for atoms that are too close and attraction at long range from dispersion forces (van der Waals interaction).

5. Interactions from charges, dipoles, quadrupoles (electrostatic interactions).

   

Different kinds of force-fields have been developed. Some include additional energy terms that describe other kinds of deformations, such as the coupling between bending and stretching in adjacent bonds, in order to improve the accuracy of the mechanical model.

The reliability of a molecular mechanical force-field depends on the parameters and the potential energy functions used to describe the total energy of a model. Parameters must be optimized for a particular set of potential energy functions, and thus are not transferable to other force fields.

MM2

Chem3D uses a modified version of Allinger's MM2 force field. For additional MM2 references see MM2

The principal additions to Allinger's MM2 force field are:

•  A charge-dipole interaction term

•  A quartic stretching term

•  Cutoffs for electrostatic and van der Waals terms with 5th order polynomial switching function

•  Automatic pi system calculations when necessary

•  Torsional and non-bonded constraints.

Chem3D stores the parameters for the potential energy function in tables. To view these tables, go to View>Parameter Tables.

Each parameter is classified by a Quality number. This number indicates the reliability of the data. The quality ranges from 4, where the data is derived completely from experimental data (or ab initio data), to 1, where the data is guessed by Chem3D.

The parameter table, MM2 Constants, contains adjustable parameters that correct for failings of the potential functions in outlying situations.

Note: Editing of MM2 parameters should only be done with the greatest of caution by expert users. Within a force-field equation, parameters operate interdependently; changing one normally requires that others be changed to compensate for its effects.

Bond Stretching Energy

The bond stretching energy equation is based on Hooke's law:

   

The Ks parameter controls the stiffness of the spring's stretching (bond stretching force constant), while ro defines its equilibrium length (the standard measurement used in building models). Unique Ks and ro parameters are assigned to each pair of bonded atoms based on their atom types (C-C, C-H, O-C). The parameters are stored in the Bond Stretching parameter table. The constant, 71.94, is a conversion factor to obtain the final units as kcal/mole.

The result of this equation is the energy contribution associated with the deformation of the bond from its equilibrium bond length.

This simple parabolic model fails when bonds are stretched toward the point of dissociation. The Morse function would be the best correction for this problem. However, the Morse Function leads to a large increase in computation time. As an alternative, cubic stretch and quartic stretch constants are added to provide a result approaching a Morse-function correction.

The cubic stretch term allows for an asymmetric shape of the potential well, allowing these long bonds to be handled. However, the cubic stretch term is not sufficient to handle abnormally long bonds. A quartic stretch term is used to correct problems caused by these very long bonds. With the addition of the cubic and quartic stretch term,

the equation for bond stretching becomes:

   

Both the cubic and quartic stretch constants are defined in the MM2 Constants table.

To precisely reproduce the energies obtained with Allinger's force field: set the cubic and quartic stretching constant to "0" in the MM2 Constants tables.

Angle Bending Energy

   

.i.Angle bending energy;

The bending energy equation is also based on Hooke's law. The Kb parameter controls the stiffness of the spring's bending (angular force constant), while defines the equilibrium angle. This equation estimates the energy associated with deformation about the equilibrium bond angle. The constant, 0.02191418, is a conversion factor to obtain the final units as kcal/mole.

Unique parameters for angle bending are assigned to each bonded triplet of atoms based on their atom types (C-C-C, C-O-C, C-C-H). For each triplet of atoms, the equilibrium angle differs depending on what other atoms the central atom is bonded to. For each angle there are three possibilities: XR2, XRH or XH2. For example, the XH2 parameter would be used for a C-C-C angle in propane, because the other atoms the central atom is bonded to are both hydrogens. For isobutane, the XRH parameter would be used, and for 2,2-dimethylpropane, the XR2 parameter would be used.

The effect of the Kb and parameters is to broaden or steepen the slope of the parabola. The larger the value of Kb, the more energy is required to deform an angle from its equilibrium value. Shallow potentials are achieved with Kb values less than 1.0.

A sextic term is added to increase the energy of angles with large deformations from their ideal value. The sextic bending constant, SF, is defined in the MM2 Constants table. With the addition of the sextic term, the equation for angle bending becomes:

Note: The default value of the sextic force constant is 0.00000007. To precisely reproduce the energies obtained with Allinger's force field, set the sextic bending constant to "0" in the MM2 Constants tables.

There are three parameter tables for the angle bending parameters:

•  Angle Bending parameters

•  3-Membered Ring Angle Bending parameters

•  4-Membered Ring Angle Bending parameters

There are three additional angle bending force constants available in the MM2 Constants table. These are the "-CHR-Bending" constants, specifically for carbons with one or two attached hydrogens.

The -CHR- Bending Kb for 1-1-1 angles1allows more accurate force constants to be specified for Type 1 (CHR) and Type 2 (CHR) interactions.

The -CHR-Bending Kb for 1-1-1 angles in 4membered rings and the -CHR- Bending Kb for 22-22-22 angles in 3-membered rings require separate constants for accurate specification.

Torsion Energy

   

This term accounts for the tendency for dihedral angles (torsionals) to have an energy minimum occurring at specific intervals of 360/n. In Chem & Bio 3D 11.0, n can equal 1, 2, or 3.

   

   

The Vn/2 parameter is the torsional force constant. It determines the amplitude of the curve. The n signifies its periodicity. n shifts the entire curve about the rotation angle axis. The parameters are determined through curve-fitting techniques. Unique parameters for torsional rotation are assigned to each bonded quartet of atoms based on their atom types (C-C-C-C, C-O-C-N, H-C-C-H).

Chem & Bio 3D 11.0 provides three torsional parameters tables:

•  Torsional parameters

•  4-Membered ring torsions

•  3-Membered ring torsions.

Non-Bonded Energy

The non-bonded energy represents the pairwise sum of the energies of all possible interacting nonbonded atoms within a predetermined "cutoff" distance.

The non-bonded energy accounts for repulsive forces experienced between atoms at close distances, and for the attractive forces felt at longer distances. It also accounts for their rapid falloff as the interacting atoms move farther apart by a few Angstroms.

van der Waals Energy

Repulsive forces dominate when the distance between interacting atoms becomes less than the sum of their contact radii. In Chem3D repulsion is modeled by an equation which combines an exponential repulsion with an attractive dispersion interaction (1/R6):

   

The parameters include:

•  Ri* and Rj*—the van der Waals radii for the atoms

•  Epsilon (—determines the depth of the attractive potential energy well and how easy it is to push atoms together

•  rij—which is the actual distance between the atoms

At short distances the above equation favors repulsive over dispersive interactions. To compensate for this at short distances (R=3.311) this term is replaced with:

   

The R* and Epsilon parameters are stored in the MM2 Atom Types table.

For certain interactions, values in the van der Waals interactions parameter table are used instead of those in the MM2 atom types table. These situations include interactions where one of the atoms is very electronegative relative to the other, such as in the case of a water molecule.

Cutoff Parameters for van der Waals Interactions

The use of cutoff distances for van der Waals terms greatly improves the computational speed for large molecules by eliminating long range, relatively insignificant, interactions from the computation.

Chem3D uses a fifth-order polynomial switching function so that the resulting force field maintains second-order continuity. The cutoff is implemented gradually, beginning at 90% of the specified cutoff distance. This distance is set in the MM2 Constants table.

The van der Waals interactions fall off as 1/r6, and can be cut off at much shorter distances, for example 10Å. This cut off speeds the computations significantly, even for relatively small molecules.

Note: To precisely reproduce the energies obtained with Allinger's force field: set the van der Waals cutoff constants to large values in the MM2 Constants table.

Electrostatic Energy

The electrostatic energy is a function of the charge on the non-bonded atoms, q, their interatomic distance, rij, and a molecular dielectric expression, D, that accounts for the attenuation of electrostatic interaction by the environment (solvent or the molecule itself).

In Chem3D, the electrostatic energy is modeled using atomic charges for charged molecules and bond dipoles for neutral molecules.

There are three possible interactions accounted for by Chem3D:

•  charge/charge

•  dipole/dipole

•  dipole/charge

Each type of interaction uses a different form of the electrostatic equation as shown below:

where the value 332.05382 converts the result to units of kcal/mole.

dipole/dipole contribution

   

where the value 14.388 converts the result from ergs/mole to kcal/mole, is the angle between the two dipoles i and j, i and j are the angles the dipoles form with the vector, rij, connecting the two at their midpoints, and D is the (effective) dielectric constant.

dipole/charge contribution

   

where the value 69.120 converts the result to units of kcal/mole.

Bond dipole parametersfor each atom pair are stored in the bond stretching parameter table. The charge, q, is stored in the atom types table. The molecular dielectric expression is set to a constant value between 1.0 and 5.0 in the MM2 Atom types table.

Note: Chem & Bio 3D 11.0 does not use a distance-dependent dielectric.

Cutoff Parameters for Electrostatic Interactions

The use of cutoff distances for electrostatic terms, as for van der Waals terms, greatly improves the computational speed for large molecules by eliminating long-range interactions from the computation.

As in the van der Waals calculations, Chem3D uses a fifth-order polynomial switching function to maintain second-order continuity in the force-field. The switching function is invoked as minimum values for charge/charge, charge/dipole, or dipole/dipole interactions are reached. These cutoff values are located in the MM2 Constants parameter table.

Since the charge-charge interaction energy between two point charges separated by a distance r is proportional to 1/r, the charge-charge cutoff must be rather large, typically 30 to 40Å, depending on the size of the molecule. The charge-dipole, dipole-dipole interactions fall off as 1/r2, 1/r3 and can be cutoff at much shorter distances, for example 25 and 18Å respectively. To precisely reproduce the energies obtained with Allinger's force field: set the cutoff constants to large values (99, for example) in the MM2 Constants table.

OOP Bending

Atoms that are arranged in a trigonal planar fashion, as in sp2 hybridization, require an additional term to account for out-of-plane (OOP) bending. MM2 uses the following equation to describe OOP bending:

The form of the equation is the same as for angle bending, however, the value used is angle of deviation from coplanarity for an atom pair and  is set to zero.

   

The special force constants for each atom pair are located in the Out of Plane bending parameters table. The sextic correction is used as previously described for Angle Bending. The sextic constant, SF, is located in the MM2 Constants table.

Pi Bonds and Atoms with Pi Bonds

For models containing pi systems, MM2 performs a Pariser-Parr-Pople pi orbital SCF computation for each system. A pi system is defined as a sequence of three or more atoms of types which appear in the Conjugate Pi system Atoms table. Because of this computation, MM2 may calculate bond orders other than 1, 1.5, 2, and so on.

Note: The method used is that of D.H. Lo and M.A. Whitehead, Can. J. Chem., 46, 2027(1968), with heterocycle parameter according to G.D. Zeiss and M.A. Whitehead, J. Chem. Soc. (A), 1727 (1971). The SCF computation yields bond orders which are used to scale the bond stretching force constants, standard bond lengths and twofold torsional barriers.

The following is a step-wise overview of the process:

1. A Fock matrix is generated based on the favorability of electron sharing between pairs of atoms in a pi system.

2. The pi molecular orbitals are computed from the Fock matrix.

3. The pi molecular orbitals are used to compute a new Fock matrix, then this new Fock matrix is used to compute better pi molecular orbitals.

Step 2 and 3 are repeated until the computation of the Fock matrix and the pi molecular orbitals converge. This method is called the self-consistent field technique or a pi-SCF calculation.

4. A pi bond order is computed from the pi molecular orbitals.

5. The pi bond order is used to modify the bond length(BLres) and force constant (Ksres) for each bond in the pi system.

6. The modified values of Ksres and BLres are used in the molecular mechanics portion of the MM2 computation to further refine the molecule.

Stretch-Bend Cross Terms

Stretch-bend cross terms are used when a coupling occurs between bond stretching and angle bending. For example, when an angle is compressed, the MM2 force field uses the stretch-bend force constants to lengthen the bonds from the central atom in the angle to the other two atoms in the angle.

   

The force constant (Ksb) differs for different atom combinations.

The seven different atom combinations where force constants are available for describing the situation follow:

•  X-B, C, N, O-Y

•  B-B, C, N, O-H

•  X-Al, S-Y

•  X-Al, S-H

•  X-Si, P-Y

•  X-Si, P-H

•  X-Ga, Ge, As, Se-Y, P-Y

where X and Y are any non-hydrogen atom.

User-Imposed Constraints

Additional terms are included in the force field when constraints are applied to torsional angles and non-bonded distances by the Optimal field in the Measurement table. These terms use a harmonic potential function, where the force constant has been set to a large value (4 for torsional constraints and 106 for non-bonded distances) in order to enforce the constraint.

For torsional constraints the additional term and force constant is described by:

   

For non-bonded distance constraints the additional term and force constant is:

.

Molecular Dynamics Simulation

Molecular dynamics simulates molecular motion. This simulation is useful because motion is inherent to all chemical processes: vibrations, like bond stretching and angle bending, give rise to IR spectra; chemical reactions, hormone-receptor binding, and other complex processes are associated with many kinds of intramolecular and intermolecular motions.

The MM2 method of molecular dynamics simulation uses Newton's equations of motion to simulate the movement of atoms.

Conformational transitions and local vibrations are the usual subjects of molecular dynamics studies. Molecular dynamics alters the values of the intramolecular degrees of freedom in a stepwise fashion. The steps in a molecular dynamics simulation represent the changes in atom position over time, for a given amount of kinetic energy.

The Molecular Dynamics (MM2) command in the Calculations menu can be used to compute a molecular dynamics trajectory for a molecule or fragment in Chem3D. A common use of molecular dynamics is to explore the confor-ma-tion-al space accessible to a molecule, and to prepare sequences of frames representing a molecule in motion. For more information on Molecular Dynamics see MM2.

Methods Available in CS MOPAC

CS MOPAC methods include both open shell and closed shell Hartree-Fock approximations. Configuration interaction may be estimated by either iterative self-consistent field (SCF) or multi-electron (MECI) methods. For semi-empirical methods, you may choose between five Hamiltonian approximations.

RHF

The default Hartree-Fock method assumes that the molecule is a closed shell and imposes spin restrictions. The spin restrictions allow the Fock matrix to be simplified. Since alpha (spin up) and beta (spin down) electrons are always paired, the basic RHF method is restricted to even electron closed shell systems.

Further approximations are made to the RHF method when an open shell system is presented. This approximation has been termed the 1/2 electron approximation by Dewar. In this method, unpaired electrons are treated as two 1/2 electrons of equal charge and opposite spin. This allows the computation to be performed as a closed shell. A CI calculation is automatically invoked to correct errors in energy values inherent to the 1/2 electron approximation. For more information see Configuration Interaction.

With the addition of the 1/2 electron approximation, RHF methods can be run on any starting configuration.

UHF

The UHF method treats alpha (spin up) and beta (spin down) electrons separately, allowing them to occupy different molecular orbitals and thus have different orbital energies. For many open and closed shell systems, this treatment of electrons results in better estimates of the energy in systems where energy levels are closely spaced, and where bond breaking is occurring.

UHF can be run on both open and closed shell systems. The major caveat to this method is the time involved. Since alpha and beta electrons are treated separately, twice as many integrals need to be solved. As your models get large, the time for the computation may make it a less satisfactory method.

Configuration Interaction

The effects of electron-electron repulsion are underestimated by SCF-RHF methods, which results in the overestimation of energies.

SCF-RHF calculations use a single determinant that includes only the electron configuration that describes the occupied orbitals for most molecules in their ground state. Further, each electron is assumed to exist in the average field created by all other electrons in the system, which tends to overestimate the repulsion between electrons. Repulsive interactions can be minimized by allowing the electrons to exist in more places (i.e. more orbitals, specifically termed virtual orbitals). The multi-electron configuration interaction (MECI) method in MOPAC addresses this problem by allowing multiple sets of electron assignments (i.e., configurations) to be used in constructing the molecular wave functions. Molecular wave functions representing different configurations are combined in a manner analogous to the LCAO approach.

For a particular molecule, configuration interaction uses these occupied orbitals as a reference electron configuration, then promotes the electrons to unoccupied (virtual) orbitals. These new states, Slater determinants or microstates in MOPAC, are then linearly combined with the ground state configuration. The linear combination of microstates yields an improved electronic configuration and hence a better representation of the molecule.

   

'상태와 변화' 카테고리의 다른 글

Combustion  (0) 2016.06.27
Oxygen  (0) 2016.06.27
Computational Chemistry  (1) 2016.06.27
오존공정의 기본 메카니즘  (0) 2016.06.27
POCP  (0) 2016.06.27