Outline of this lecture

1. Forcefields (15 minutes)
2. Electrostatics (10 minutes)
3. Minimization (20 minutes)
4. Quantum Mechanical Methods (5 minutes)
5. 10-minute break
6. Conformational Analyses
   1. Systematic Searches (5 minutes)
   2. Monte Carlo Methods (5 minutes)
   3. Molecular Dynamics (15 minutes)
7. Free Energy Peturbation Methods (5 minutes)

Forcefields

Description of the Classical (Newtonian) Forcefield:

Morse vs Harmonic Bond term:

Torsion Angle Term:

Gauche-Anti Isomerization of Butane

Dynamics at 300K

Dynamics at 900K

Scaling of nonbond terms:

• Urey-Bradley Terms: reduces scaling of nonbond geminal atom energies
• p1-4 Terms: reduces scaling of nonbond vicinal atom energies

Forcefield Parameters

1. A parameter should be adjusted so that the simulated system reproduces properties of the real system. It does NOT have to equal a microscopic descriptor.
2. As yet, there is no universal forcefield.
   Parameters for particular forcefields can only be used to study particular molecules.
3. Forcefield parameters must be cohesive---modelers can't just "plug & play" new parameters without testing
4. Forcefield parameters must be referenced

Types of Forcefields

• diagonal force field: no cross terms (AMBER, CHARMM)
• matrix force field: with cross terms (MM4, CFF)
• Class II force field: parameters derived from quantum mechanics (ESFF)
• "Fast" force fields: only torsion & non-bond terms (SCULPT, YETI)

1. H = harmonic term; M = Morse term
2. also known as Out-Of-Plane or Improper Dihedral Angle term
3. B = Buckingham-like term
   LJ6-9 = Lennard-Jones 6-9 term
   LJ6-12 = Lennard-Jones 6-12 term
   LJ10-12 = Lennard-Jones 10-12 term
   LJ-like 7-14 = Lennard-Jones 7-14 term with switching function at large distances
   
4. C-C: Coulombic charge-charge electrostatic term
   D-D: Dipole-Dipole electrostatic term

   Comparison of Force Field Accuracies

   

Periodic Boundary Conditions

Minimum Image Model:

Explicit Image Model:

• Much easier to set up cubic cells.
• Non-bond cutoffs must not be greater than half the cell's smallest dimension
• Must consider molecular density

Electrostatics

• Electrostatic effects dominate long-range interactions---the most important forcefield term for biomolecules.
• Electrostatic field is approximately uniform for small molecules; dielectric constant can be set to dielectric of solvent for all atoms.
• Electrostatic field is NOT uniform for atoms in large molecules; dielectric constant should be smaller for atoms in "core", and larger for atoms on surface.

  Atom location	Dielectric Constant
  in molecular "core"	3 to 5
  on molecular surface	about 10
  in bulk water	78

  Distance-dependent dielectric is used to approximate different dielectric constants for "core" and surface atoms; works particularly well for docking & ligand design.
• Can use explicit solvent for best accuracy; not all solvents are parameterized, large calculation can be computationally expensive.
• Use NO charges for conformational analyses or crude geometry optimization.
• Use the Non-bond Cutoff Method to limit number of electrostic calculations for atom pairs.
• The Cell-Multipole method and Ewald Sum Method scale as N. (N = number of atoms). The standard electrostatic Forcefield term scales as N². (N = number of atoms).

The Non-bond Cutoff Method

Non-bond interactions dominate calulcation time:

Example of a Cutoff Function:

Calculations may require a very long cutoff to retain accuracy:

Cell Multipole Method

Ewald Sum Method

Constraints

• Fixed atoms
• Tethered atoms
• Distance constraints
• Torsion constraints
• Template Force constraints

Minimization

Steepest Decents Algorithm:

Conjugate Gradients Algorithm:

• Quasi-Newton Raphson
• Truncated Newton Raphson
• VA09A

Rules for selecting Minimization Algorithms

When Derivative > 10	Use the Steepest Decents Algorithm
When Derivative < 10	For systems < 200 atoms, use the Newton Raphson Algorithm
	For systems > 200 atoms, use the Conjugate Gradients Algorithm

What is a Minimum?

A typical protocol for biomolecular Energy Minimization:

1. fix heavy atoms, minimize
2. tether heavy atoms, minimize
3. fix heavy atoms of structurally conserved regions, minimize
4. tether heavy atoms of structurally conserved regions, minimize
5. fix backbone of structurally conserved regions, minimize
6. tether backbone of structurally conserved regions, minimize
7. remove all constraints, minimize

For details about performing minimizations, see the MolViz Discover Minimization Notes.

For an example of the Discover Minimization input, see this Example Input file.

For an example of the Discover Minimization output (with NO convergence), see this Example Output file (no convergence).

For an example of the Discover Minimization output (with convergence), see this Example Output file (with convergence).

---

Quantum Mechanical Methods

Quantum Mechanical Methods must be used for applications that are beyond the capability of classical methods:

• Electronic transitions (photon absorption).
• Electron transport phenomena.
• Bond breaking/formation.
• Proton transfer (acid/base reactions).

Ab Initio Methods: Ab Initio methods are easier than you think! See the Gaussian 92 Example for more details.

No empirical (experimental) information is needed. Best for types of molecules with little experimental information. Must choose appropriate basis set (set of atomic orbitals necessary to accommodate all electrons of the atoms in the ground state). Polarization effects, electron delocalization hyperconjugate effects, D orbitals require higher-order basis sets. Results should always be evaluated. Higher-order basis sets require more computational cost.

Semi-empirical methods: Semi-Empirical methods are easier than you think! See the Mopac 93 geometry optimization input file and the Mopac 93 geometry optimization output file for more details.
For more information about Mopac, see also the Mopac application notes and notes about How to use Mopac with PCMODEL.

Usually considers only valence electrons. substitutes empirical parameters for some integrals. Faster than Ab Initio methods by several orders of magnitude. Results should always be carefully evaluated.

Hybrid QM/MM Methods: Treat "active site" as QM system, with surrounding "medium" as MM force at "active site"/"medium" boundary. QM system usually evaluated by semi-empirical methods. Works well only if energy of "active site" is not insubstantially dependent on "medium". Recently, a hybrid Ab Initio/semi-empirical/MM methods has been applied to systems with metals.

Conformational Analyses:
Systematic Searches

• "Bump check"
• Energy cutoff
• "Active Analog Approach": comparison to "active analogs" & inactive representations. Works epecially well if the analogs & representations ahve few degrees of freedom.

Rigid Rotor technique: no mimimization after each rotation. OK for finding geometries with energy minima. Terrible for measuring rotational energy barriers.

Dihedral Driver technique: minimization after each rotation; may have to "drive" from different directions.
Good for measuring rotational energy barriers.

Conformational Analyses:
Monte Carlo Searches

1. move atomic coordinates in a random direction.
2. minimize the structure.
3. If the conformation's absolute energy is below a pre-set energy threshold, save the results.
4. Return to the first step and repeat.

Example of a Metropolis Monte Carlo algorithm (a.k.a. Boltzmann-Weighted Monte Carlo algorithm or Weighted Monte Carlo algorithm):

1. Set the "temperature of the system" (e.g., 500K).
2. Evaluate the energy of the current conformation.
3. Move atomic coordinates in a random direction.
4. Evaluate the energy of the new conformation.
5.  
   • If the energy of the new conformation is less than the energy of the old conformation, keep the new conformation.
   • If the energy of the new conformation is more than the energy of the old conformation, calculate the Boltzman probability (between 0 and 1) that this energy increase would occur at the system's current temperature. Generate a random number (between 0 and 1). If the random number is less than the Boltzman probability, keep the new conformation. If the random number is greater than the Boltzman probability, return to the old conformation.
6. Return to step 2. Repeat steps 2-5 thousands of times.
7. Lower the temperature of the system by a small increment (e.g., 1 K).
8. Return to step 2. Repeat steps 2-6 hundreds of times, until the temperature of the system is near 1 K.

Monte Carlo methods require very long samplings.
To prove that the MC search has sampled sufficiently long enough, start from different starting points: if each run gives the same result, the caluclation has converged on the same ensemble.
Monte Carlo methods are the best techniques for ring systems.

Conformational Analyses:
Molecular Dynamics

Dynamics Timestep

Must use small timestep to avoid errors.

The RATTLE method:

Relation of standard deviation in total energy and the CPU time per step with and without RATTLE:

            No RATTLE                    RATTLE
            ---------        ------------------------------------
                              tol. = 1e-7          tol. = 1e-5
                             ----------------    ----------------
Timestep  SDV     Time/step   SDV   Time/step    SDV     Time/Step
 (fs)  (kcal mol-1) (sec)    (kcal mol-1) (sec)  (kcal mol-1)  (sec)

 0.5     0.740    0.889      0.058   0.971      0.058     0.937
 1.0     2.941    0.893      0.232   0.985      0.232     0.944
 1.5     6.352    0.884      0.525   0.988      0.526     0.953
 2.0     9.832    0.886      0.939   0.993      0.939     0.957
 2.5    10.217    0.887      1.494   0.998      1.494     0.967
 3.0            crashed      2.224   1.004      2.224     0.961
 3.5                         3.239   1.008      3.239     0.966
 4.0                         4.859   1.006      4.858     0.967
 4.5                         6.821   1.011      6.820     0.970
 5.0                        11.126   1.012     11.125     0.973
 5.5                           crashed             crashed

(``tol'' = tolerance) 

Simulated Annealing

Molecular Dynamics Analysis

Miscellaeous Notes

• Equilibrate for at least 10 psec.
• High-temerature MD may require additional restraint on peptide bond to avoid cis-trans isomerization.
• It is difficult to ensure complete sampling with Molecular Dynamics.

Molecular Dynamics with Massively Parallel Processor (MPP) Systems MPP applications must consider 3 conecpts: granularity of tasks Most MD is coarse-grained: parallelism is at level of subroutines or functions, not at fine-grained loop-level. load balancing Non-Uniform Memory Access (NUMA): atom decomposition method: atoms are assigned to specific processors Inefficient beyond 64 CPUs force decomposition method: force calculations are assigned to specific processors Inefficient beyond 64 CPUs; best method for high-performance FORTRAN and other SIMD (Single Instruction Multiple Data) formalisms (in general, MIMD formalisms are more appropriate for MD). spacial decomposition method: blocks of space are assigned to specific processors Most efficient method, because it can most naturally take advantage of cell multipole methods. Hardest method to code. Most popular modeling applications (AMBER, CHARMM, GROMOS, Discover) have bene parallelized for up to 64 CPUs. The Lennard-Jones benchmark problem is used by most researchers to test MPP MD code. However, speed is still lacking for most modeling studies.

Free Energy Peturbation method