The Calculation of Enthalpy and Entropy Differences

1. The Calculation of Enthalpy and Entropy Differences???(Housekeeping Details for the Calculation of Free Energy Differences)first edition: p. 493-502second edition: p. 574-585

2. The Calculation of Enthalpy and Entropy Differences • Free energies can now be calculated with errors of less than 1 kcal/mol in favorable cases. • Enthalpy and entropy differences for solvation could be calculated by simulating the two systems separately and taking the differences in the total. • Leads to much larger errors than for the free energy since the free energy reduces to interaction terms only involving the solute. • For example: The solvent-solvent interaction term which contributes the so-called cavity (solvent reorganization) term to the energy is said to be canceled exactly by a corresponding term in the entropy. • Yu, H.-A.; Karplus, M. J. Chem. Phys. 1988, 89, 2366-2379.

3. Partitioning the Free Energy • If the thermodynamic integration method is used the overall free energy can be partitioned into individual contributions. • However, while the total free energy is a state function the individual contributions are not.

4. Partitioning the Free Energy • Calculation of the free energy differences by thermodynamic integration: • When performing this procedure on individual contributions, energy is transferred between the contributors. For example relieving strain in a bond angle may increase the potential energy in certain bond distances.

5. Partitioning the Free Energy • A common practice is to partition the free energy into contributions from the van der Waals and electrostatic interactions. • The biotin/streptavidin complex has an extremely strong association constant (-18.3 kcal/mol). • The favorable electrostatic interaction, from H-bonding, was canceled by the free energy of interaction of biotin with water. • However there was a very large van der Waals interaction. biotin

6. Partitioning the Free Energy • Strong van der Waals interaction.

7. Potential Pitfalls with Free Energy Calculations • Two major sources of error: 1) Hamiltonian. 2) insufficient sampling of phase space. • Inadequate sampling: • errors may be identified by running the simulation for longer periods of time (molecular dynamics (MD)) or for more iterations (Monte Carlo (MC)). • The perturbation may be run in the forward and reverse directions. The difference in the calculated energy values, hysteresis, gives a lower bounds estimate of the error. • Note, very short simulations may give almost 0 hysteresis while the errors may still be large.

8. Implementation Aspects • Simulation Method: • Molecular dynamics is almost always used for systems with significant degree of conformational flexibility. • Monte Carlo gives good results for small rigid molecules. • Thermodynamic perturbation or integration preferred over the slow growth methods. • Slow growth suffers from “Hamiltonian lag” and adding additional values of  requires rerunning the simulation from scratch.

9. Implementation Aspects • Coupling Parameter () •  does not have to be a constant value. Could use small values when the free energy is changing quickly and large values when the free energy is changing slowly. (Dynamically modified windows) • Choice of Pathway • A change that involves high energy barriers will require much smaller increments in  to insure reversibility than a pathway that proceeds via a lower barrier.

10. Implementation Aspects • The molecular topology at all stages is a union of the initial and final states, using dummy atoms where necessary. • Single-topology: • Dual Topology: • Both molecular topologies are maintained, but do not interact with each other. • The simplest Hamiltonian that describes the interaction between these groups and the environment is a linear relationship:

11. Implementation Aspects • Dual Topology: • Can result in a singularity in the function for which an ensemble average is to be formed. • Problem with thermodynamic integration where the derivative of the parameterized Hamiltonian with respect to  is the observable. • One solution when performing MC is to change the scaling factors: • When n  4 the free energy difference is always finite and can be integrated numerically. • However this results in difficulties in calculating the first and second derivatives of the potential energy function required for MD. • Solution: Soft Core Potentials.

12.  = 1 Potential energy  = 0 Atom-atom separation Implementation Aspects • Where  determines the softness of the interaction, removing the singularity. • Soft Core Potential: • The traditional Lennard-Jones interaction can be replaced: • Similar expressions can be developed for for electrostatic interactions.

13. Potentials of Mean Force • May wish to examine the Free Energy as a function of some inter- or intramolecular coordinate. (ie. Distance, torsion angle etc.) • The free energy along the chosen coordinate is known as the Potential of Mean Force (PMF). • Calculated for physically achievable processes so the point of highest energy corresponds to a TS. • Simplest type of PMF is the free energy change as the separation (r) between two particles is varied. • PME can be calculated from the radial distribution function (g(r)) using: • Recall: g(r) is the probability of finding an atom at a distance r from another atom.

14. Potentials of Mean Force • Problem: The logarithmic relationship between the PMF and g(r) means a relatively small change in the free energy (small multiple of kBT may correspond to g(r) changing by an order of magnitude. • MC and MD methods do not adequately sample regions where the radical distribution function differs drastically from the most likely value. • Solution: Umbrella Sampling. • The coordinates of interest are allowed to vary over their range of values throughout the simulation. (Subject to a potential modified using a forcing function.)

15. Umbrella Sampling • The Potential Function can be written as a perturbation: • Where W(rN) is a weighting function which often takes a quadratic form: • Result: For configurations far from the equilibrium state, r0N, the weighting function will be large so the simulation will be biased along some relevant reaction coordinate. • The Boltzmann averages can be extracted from the non-Boltzmann distribution using: • Subscript W indicates that the average is based on the probability PW(rN), determined from the modified energy function V ‘(rN).

16. Umbrella Sampling • This free energy perturbation method can be used with MD and MC. • Calculation can be broken into a series of steps characterized by a coupling parameter  with holonomic constraint methods used to fix the desired coordinates. (ie. SHAKE)