Data used for parameterization:
1. Variety of Sources
2. Measured via Different Kinds of Experiments
3. Different Units
4. Different Error Bars; Accuracy

Based On:
1. Structural Data
2. Structural & Energetic
3. Structural,E's & 8's
4. Strucure,E's,V's,Thermodynamic

Procedures:

1. Opt by Inspection
2. Least Squares
• y (improved) = y (initial) + dR(bo , kb , ao, ko ,...Alj, Blj, ql , qj...)
• Weight Function:Different Units
• Linearize via Taylor Series Expansion.
• Intuition

General Application of Force Field

• Structure, Energy, and Properties of Isolated Molecules
• Isomer Comparison
  
• Substituent Effects
• Not good for weakly bound systems.
• Not good for bond making/bond breaking processes
• Care needs to be taken in aromatic systems. These systems usually involve >1 Kekule description.

Conjugated Systems

Bonding Less Easy to Define
A. Simple Linear Polyenes:
• Methods as described work reasonably well.
• 'Essentially' double bonds (e.g. 1,2 bond in butadiene) are treated as ordinary double bonds, while 'essentially' single bonds (e.g., 2,3 bond in butadiene) are treated as single bonds
• It is known that the resonance interaction in a linear conjugated pyrene is not great.
• The approximation of taking the bond E's of such bonds as additiv quantities appears adequate (Denar, 1969).

B. Cyclic Polyenes Late 60's: Boyd, Chang, Allinger

Methods as described do not work.
• Benzene
  Can we find a 'natural' bond length and force constant for bond stretching which would be somewhere between those for C-C and C=C? Will they reproduce benzene structure?

  Yes: Allows simple benzene derivatives to be dealt with using these parameters by methods previously discussed.

  New parameters to add: Car - Cca ; kar

  Good results for benzenoid hydrocarbons using these parameters assigned to benzene ring.

• Napthalene Allinger, Sprague 172-3
  • Using simple benzenoid parameters -> all bonds=
  • Using 'single' and 'double' bonds ->alternating bonds =>exptal
  • Why? In napthalene, the bond lengths are determined largely by bond orders.
    Benzene bond orders are all=and included as such in the parameters
  • => Can not apply FF calculation as discussed previously except on a case by case basis
  • What about more generally?

  Bond Order: Measure of E density between pi centers
  

  • Would need potentials depending on pi-electron distribution to calibrate this effect.
    
  • In particular, expect bond potential (force) to vary as the pi electron density between atoms changes.
    
  • Requires some sort of quantum mechanical techniques
  • Full quantum mechanical techniques defeat the purpose.
  • Instead - pi-quantum mechanical treatment only
    • Partition sigma and pi frame
    • Essentially orthogonal & therefore independent
    • Good assumption for planar systems



Two approaches

• I. Warshel and Karplus
  • Initial structure
  • Do SCF on PI system E=F(coords) Analytic
  • Do MM on sigma system (Morse potential); conjugated bond distortions are outside Hooke's law.

    Total Energy = MM Energy (sigma) + QM Energy (pi)

  • Do structural optimization to minimize this sum
  • Good results for planar systems
  • Nonplanar: SCF calculations modified to include nearest neighbor overlap

  II. Allinger and Sprague
  • Initial geometry
  • SCF done on pi system. If nonplanar, flatten to do SCF
  • Calculate bond orders from coefficients of wave function.
  • Assume linear relationship between bond orders of a bond and
    Rb
    bo
    torsional barrier
    All necessary Ri and bi calculations using these structural MM proceeds with these new parameters.
  • Works well with planar systems:C-pi approximation is good.
    Breaks down with nonplanar systems.
  • Bond Order
    Pi's=(Sum)N(cicj)
    N: Number of electrons in each occupied orbital.
    Ci: coefficient of atomic orbital i.
    1. bo=l=l1-rPis l1 r empirically determined
       (Hydrocarbons: l0=1.51-0.17JPlj)
    2. kb=5.0+4.6Plj, for hydrocarbons plot exptal k vs. Pij
    3. Worst assumption is the assumption that there is a linear relationship between bond order and torsional potential.
       V2y=15.0f(PljB;j)
       Blj:resonance inegral; calculate heat of formation
       2-fold torsional constant across bond lij
       ________
       Plj bond order between cij |
       Bij: resonance integral between cij |
       .....................................Z: doubly occ. orbs | |
       As bond is twisted, Pij ->0
       Bij ->0
       Vz -Pij2
    • Rotational barriers for ethylene, stilbene and butadiene fit against this function and the proportionality constant obtained.
    • Works well for many nonplanar aromatics
      More Recent modifications (MMP2)
      • MMP1 gives error Ethylene: pi pond order always = 1.0 regardless of distance between atoms.
        +> Conjugaton E would be independent of dist between atoms!
        
      • Want El->O as atoms separate
        => Ec = RPljBlj as perturbution theory requires k. parameter
More recent controversy involving conjugated systems:
1. Through bond coupling: dominance in filled-filled orbitals through-bond without bond lengthening
2. [18] Annulenes {Many errors due to faulty experimental results=>bad parameterizations}

Minimization Techniques

• V=Vb + Vo + Vt + Vnb +...
• Once suitable potential functions and parameters are chosen, we can calculate the potential energy for any arrangement of the component atoms that do not vary much from the reference geometry.
  Intermolecular interaction potentials
  Parameters
  Pair interactions
  Site-site, group-group or atom-atom cutoffs
  Nonbonded' interactions
  Switching functions
  Periodic B.C.

  Objective: Find a set of coordinates corresponding to a minimum on the PE's
  • Minimum for macromolecules implies eliminating the worst steric conflicts
    
  • Provide relative energy information.
    
  • Requires taking the derivative of the potential with respect to coordinates and setting this functional equal to zero. i.e., derivative of potential function with respect to all degrees of freedom in the molecule.
    
  • Characterization of any minimum structure found is needed. This is done by calculating the second derivative of the potential energy with respect to coordinates (i.e., calculating the hessian).
    
  Use:
  1. Comparison of Structures/Properties
  2. Binding Energies
  3. Harmonic Analysis
  4. Comparing/Fitting FF's
  5. Template forcing
  6. Systematic mapping of E space
  7. Limiting scope of MC & MD
Algorithms
• Many methods are available for minimizing 5 - 10 variables.
• These methods are not necessarily suitable for MM calculations involving greater than 100 atoms.
• Methods are either:
  1. Convergent: Stop when they have achieved a minimum of confirmation.
  2. Nonconvergent: Never stop wondering around on the PES.
• Methods usually used:
  1. Steepest descents
  2. Conjugate gradients
  3. Newton Raphson (NR)
     Quasi NR
     Full NR
  4. Quenched dynamics
• Steps
  • Energy equation defined and evaluated for a given conformation
    Possible restraining functions. (bias optimization towards some particular minimum.)
  • Conformation adjusted to lower the energy of the potential energy
    Number of iterations depends on
    • Nature of algorithm
    • Form of potential energy function
    • Size of molecule
Global versus local minima:
• Small molecules - global
• Large molecules
  • Local minimum is usually determined.
  • Difficult to find global minimum of a general nonlinear function of greater than 10 independent variables

  Solving Nonlinear Optimization Problem:
  f(x)=f(xj)+(x-xo)f1(xo)+(x-xo)2f''(xo)/2+....
  

  • x: vector of coordinates
    
  • derivatives are matrices
    
  • order:highest order derivative used on method.
    1. zeroth order: grid search
    2. first order: Steepest descents
    3. conjugate gradients
    4. second order: Newton Ruphson

Steepest Descents:

• 'Walk along gradient'
  
• f(x0=F(XJ)+(X-XO)F'(XO) Force = negative gradient of potential.
• Useful for macromolecular structure
• Steps:
  1. Choose descent direction: 3N-D vector of unit lenghth
  2. Determine descent step size
  3. Descent step taken:
     Direction parallel to net force; points downhill steepest descent
• Step size Initial Step Chosen
  Step taken
  E lowered?
  yes....................................................................no
  Increase step see(1.2).....................................Decrease step size (0.5)
• With and without line minimization
• Converges quickly when derivatives large-Far min.
• Converges slowly when derivatives small-near min.
Assessment

• Sensitive to step size
  • Too large of a step: unreliable structures
  • Too small a step: large amounts of time

Conjugate Gradients:

• non-interfering directions
• 
  Xr=Xr-1 T ZrSr
• Descent method based on gradient at current step, and on gradients at previous steps.
• S1=-Gi
  Sr=-Gr+brSr1
  br=Gr2/(Gr-1)2
• first line search = Steepest Descents with line minimization

Newton Raphson:

• Assumption: Energy is a quadratic function in the vicinity of the minimum.
  Thus,
  
  • v9x0=a+bx+cx2
  • v1(x)=b+2cx
  • v11(x)=ac
  At min,
  1. v1(x*)=b+2cx*=0
     x*=-b/ac
  2. v1(x)=b+2cx
  3. v1(x)/ac=b/ac+x
  4. v1(x)/v11(x)=-x*+x x*=x-V1(x)/V11(x) x*=x-V1(x)[V11(x)]
     or
     V1(x)+V11(x) X=0 i=l1,...,3n-6
     V9x)+F(x) X=0
     x1=x0-F-1 V(X)
     F-1:difficult!
     can increase speed if min by factor of 10-100

3n x 3n matrix with all coordinates is 6-fold singular due to 3 translational and 3 rotational degrees of freedom (only 3N-6 nonzero vibrational frequencies) 3 ways around this:

1. Reduced matrix method.
   translation & rotation prevention by fixing
   1 atom@x=y=z=0
   2nd atom@ x=y=0........along 1 axis 3rd in a plane x=0
   =>can delete 6 egns,=>removes 6 rows and 6 ncolumns if F matrix
   =>no more singularities
   =>Problem w. oscillating when trial geom pcor
   Care must be taken not to restrict coordiantes
   F a bond length or< is restricted
2. Generalized inverse
   • Diagonalize 3n x 3n F matrix
   • Remove 6 0 eigenvalues and corresponding eigenvectors.
   • Obtain inverse
     .......................eigenvector
     >F3n,3n-1=A3n,#n-6,3n-6 (-1) A3n-6,3n
3. Apply Eckart constraints
   Eckart conditions fix the center if mass in space and constrain infinitesimal rotations of the molecule.
   In this method, a (3n+6) x (3n+6) matrix must be inverted, but E min is faster than other methods.
   This is just the opposite of steepest descents which
   converges slow-or diverges- when one is far from the minimum.
   ...........fast near minimum
   Characterictics of a block-diag min is like steepest descents works well for that surface <