Geometry Optimization

1
Geometry Optimization
Computational Chemistry 5510 Spring 2006 Hai Lin
2
Potential Energy Surface (PES)
A hypersurface in 3N-dimensional hyperspace,
where N is the number of atoms.
3
Typical Tasks for Geometry Optimization

• Find a Local Minimum Structure
• Find the Global Minimum Structure
• Find the Transition State Structure
• Find the Reaction Pathway Connecting Two Minima
  and Passing through the Transition State
  Structure (will be talked about later this
  semester)

4
Gradients

• Gradient gx ?E/?x
• First-order derivatives of energy w/r variables
  (e.g., Cartesian coordinates X, Y, Z, or
  internal coordinates such as bondlength and angle
  displacements)
• The negative of the gradient is called force.
• Stationary Point
• A point on the PES where gradient is zero
• Including minimum, maximum, and saddle points
• (A transition state is a one-order saddle point.)

E
x
5
Hessians

• Hessian Hxy ?2E / ?x?y
• A matrix of second-order derivatives of the
  energy w/r variables (e.g., Cartesian or internal
  coordinates)
• Sometimes called force matrix
• Normal-mode Analysis
• Diagonalization of the Hessian matrix to obtain
  vibrational frequencies (related to the
  eigenvalues) and normal modes of vibrations
  (eigenvectors)
• Note Strictly speaking, normal-mode analysis is
  only meaningful at stationary points.

6
Minima, Maxima, Saddle Points

• A normal-mode analysis (also known as frequency
  calculation) can be used to identify a stationary
  point.
• Minimum All frequencies are real (eigenvalues of
  the Hessian are positive).
• Maximum All frequencies are imaginary
  (eigenvalues of the Hessian are negative).
• Saddle point Some frequencies are imaginary.
• In particular, a saddle point with only one
  imaginary frequency is called a transition state.

7
Locate a Local Minimum

• Zero-order methods (energy only seldom used)
• First-order methods (energy and gradient)
• Second-order methods (energy, gradient, and
  Hessian)
• Use of more energetic information usually makes
  it easier to find a minimum at a cost of higher
  computational effort.
• In numerical calculations, due to the finite
  precision, a zero gradient is not exactly zero
  and is smaller than a pre-defined cut-off value.
  (Loose, tight, very tight ...)

8
Steepest Descent Method

• A first-order method
• Take a step in the direction of negative gradient
  
• d - g
• The step size can be fixed or depend on the norm
  of gradient.

• Once the energy increases, perform a line search
  using the points in the last three steps to
  determine where to go for the next step.

4
5
9
Steepest Descent Method (2)

• Guarantee lowing the energy (when step size
  dependeds on gradient)
• Oscillating around the minimum energy path
• Never reach the minimum
• An ever descreasing speed approaching the minimum

• Relax a poor initial geometry quickly, and is
  very useful for preliminary optimization

10
Conjugate Gradient Methods

• Instead of follow the direction current gradient,
  it follows a direction conjugate to previous
  searching directions
• di - gi bi di 1
• where bi is determined by gradients at the
  current and the last or last few points.
• Much better convergence characteristics than the
  steepest descent method
• Need to store the history Implicitly use
  information about second-order derivatives
• Used quite often, especially for biomolecules

11
Newton-Raphson Methods

• Second-order methods
• Taylor expansion around a given point
• f(x) ? f(x0) gT(x x0) ½ (x x0)T H (x
  x0)
• ?
• (x x0) g / H - H1 g
• The displacements of the coordinates for the next
  step follow the eigenvectors of the current
  Hessian, and the step size following the i-th
  eigenvector is determined as
• Dxi fi / ei
• where fi is the projection of gradient along the
  i-th eigenvector, whose eigenvalue is ei.

12
Newton-Raphson Methods (2)

• Converge very efficiently, especially when close
  to stationary points
• Converge to the nearest stationary point,
  regardless it is a minimum, maximum, or saddle
  point.
• A Trust Radius is used to prevent too large step
  sizes that move the system out of the PES region
  where the Tylaor expansion is valid.

• A shift parameter l can also be introduced to
  modify the Hessian in order to control step
  sizes
• Dxi fi / (ei - l)

13
Newton-Raphson Methods (3)

• Hessian can be calculated exactly, but usually
  very expensively. Alternative solution is to use
  an approximated Hessian that is updated (refined)
  with energy and gradient informations on the fly
  during the optimization.
• The NR methods are not practical for optimizing
  large-size systems, because calculating and
  handling very large Hessians is difficult. (How
  many Hessian elements do you have for a protein
  of 10,000 atoms?)

14
Restrained Optimization

• Keep the variable close to a desired value.
• Impose a penalty function, often a harmonic
  potential function, to the energy expression
• E E kp(x x0)2
• The force constant kp can be initially set very
  large and descreases gradually.
• Example
• Relaxation of a X-ray
• structure of a protein.

E
x0
x
15
Constrained Optimization

• Satisfy certain constraints.
• G(x1, x2, ..., xN) 0
• Example keep the reaction coordinate defined by
  a bondlength to a desired value (r r0 0).
• Lagrange undetermined multipliers
• L(x1, x2, ..., xN, l) E(x1, x2, ..., xN) l
  G(x1, x2, ..., xN)
• where E(x1, x2, ..., xN) is the energy function
  of the system.
• Now we require
• ?L / ?xi 0, ?L / ?l 0

r
16
Choice of Coordinates

• Cartesian Coordinates
• Straightforward, non-redundant, may not be
  efficient due to high correlations, but sometimes
  the only choice for large-size molecules
• Internal Coordinates
• Close to chemical intuitions, usually redundant,
  cause difficulty for certain geometries (e.g., an
  angle accidentally close to 180 deg), need a
  scheme to automatically generate for large
  molecules
• Mixed Cartesian and Internal Coordinates
• Convenient to use, not available in most programs

17
Combinatorial Explosion Problem
The number of local minima typically goes
exponentially with the number of variables
(degrees of freedom).
Possible Conformations (3n) for linear alkanes
CH3(CH2)n1CH3
f
n 1 3 n 2 9 n 5 243 n 10 59,049 n
15 14,348,907 n 100 ?
f1
f2
18
Locate the Global Minimum

• Systematic (Grid) Search (Manual or automatic)
• Monte Carlo (Stochastic) Methods
• Randam jumps on the PES
• Molecular dynamics
• Moving around based on Newtons laws of motion
• Simulated Annealing
• Cooling a melt slowly to form a crystal
• Genetic Algorithms
• Darwins principle survival of the fittest and
  carry on to the next generation

19
Monte Carlo Methods
20
Monte Carlo Methods (2)

• Sample the PES only.
• Provide no time-dependent information.
• Special techniques are often used to promote a
  reasonable acceptance ratio.
• Energy minimization following distorted geometry
  can be used to increase efficiency.
• Geometry distortion is based on the last geometry
  or on a set of selected previous geometries.

21
Molecular Dynamics

• Sample the 6N-dimensional (kinetic potential)
  hyperspace.
• Provide time-dependent information.
• Special techniques can be used to promote a large
  step size of time.
• The total enegy set the limit of barrier heights
  that can be overcomed.
• Energy minimization is done for geometry
  snapshots extracted from the MD trajectory.

22
Simulated Annealing

• High temperature
• Overcome high barriers
• Sample large area

Gradually decease the temperature in the MC or MD
run

• Low temperature
• Overcome low barriers
• Sample small area

• If the annealing process is infinitely slow (not
  possible in practice), the system reaches the
  global minimum.

23
Genetic Algorithms
Select a number of optimized conformations
(parents), calculate their fitness value
Define a fitness function f
Use the survived conformations as the next
generation (new parents)
Generate new conformations (children) by
crossover, mutation, ..., followed by
optimizations
Select conformations of highest fitness values
from parents and children
Calculate fitness values for children
Seems robust in finding minima pretty close to
the global minimum becomes popular in recent
years.
24
Summary

• Typical Tasks for Geometry Optimizations
• Describe the PES
• Gradient, Hessian, Minima, Maxima, Saddle Points
• Find a local minimum
• Steepest Descent, Conjugated gradient,
  Newton-Rasphon
• Restrained and Constrained Optimizations
• Choice of Coordinates
• Find the global minimum
• Monte Carlo, Molecular Dynamics, Simulated
  Annealing, Genetic Algorithms

25
Your Homework

• Read Textbook (Remember to take notes!)
• 14.1, 14.2. 14.3
• 14.4
• 14.6
• 14.7.1, 14.7.2, 14.7.3, 14.7.4

26
Practise with Tinker

• Optimization Programs Provided by Tinker
• For local minimum
• minimize (use energy and gradient)
• Newton (use energy, gradient, and Hessian)
• For global minimum
• Anneal (simulated annealing)
• Monte (Monte Carlo)

27
Practise with Tinker (2)

• A small piece of peptide with known c5, c7a, and
  c7e conformations.
• (Coordinate files under the test subdirectory)
• Calculate the single-point energy for all three
  geometries. Which one is of the lowest energy?
• Optimize all three geometries using programs
  minimize and Newton. Which one is of the lowest
  energy?
• Distort the geometries by editing manually the
  coordinate files (.xyz), repeat step 2.

28
Practise with Tinker (3)

• A small piece of peptide with known c5, c7a, and
  c7e conformations.
• Optimize all three geometries (c5, c7a, c7e)
  using programs monte and anneal. Can you get the
  global minimum?
• Distort the geometries by editing manually the
  coordinate files (.xyz), repeat step 4.

29
Your Home work

• Try to find the global minimum or a minimum that
  is close to the global minimum for
• methane
• ethane
• butane
• octane
• long-chain n-alkane
• The files are named methane, ethane, butane,
  octane, and alkane under the test subdirectory.
  Use the MM3 force field.