Forcefield Based Simulations

6       Free Energy

Who should read this chapter

Forcefield-based calculation of relative or absolute free energy should be considered an "expert" application. We recommend using the FDiscover program for relative and absolute free energy calculations. Although it is possible to run relative free energy calculations with CHARMm, that program is slower and less flexible.

This chapter explains

Relative free energy--theory and implementation

Absolute free energy

Finite difference thermodynamic integration (FDTI)

Relative free energy--methodology

To compare chemically distinct systems, the FDiscover program can be used to calculate the relative free energy difference between chemically unique species. This approach uses the finite difference thermodynamic integration (FDTI) algorithm of Mezei et al. (1987). Since FDTI does not require analytical derivatives of the Hamiltonian with respect to the coupling parameter, it is more suited for the complex coupling formalisms required in chemical perturbations. It also has been shown to possess better convergence properties than other methods.

Finite difference thermodynamic integration (FDTI)

Perturbation method (PM)

The free energy A is related to the partition function Q by the equation:

Eq. 114        

Eq. 115        

Eq. 116        

Eq. 117        

Eq. 118        

where the subscript 1 indicates that State 1 is considered the reference state (that is, the energy difference E is computed relative to an ensemble of structures for State 1). Thus, the ratio of the partition functions can be computed from an ensemble average of the energy difference between a reference state and a perturbed state.

The free energy change is then given directly as:

Eq. 119        

To calculate free energy differences larger than 2 k_BT, several sequential runs are performed, each computing the free energy change over a subinterval. Errors introduced by failure to converge at each step are propagated--each successive calculation adds its contribution to the previously accumulated total.

FDTI

Thermodynamic integration of the relative free energy assumes that the free energy change can be expressed as an integral:

Assuming that a suitable coupling parameter can be found, which adequately describes a continuous conversion between the two states, the above equation can be integrated numerically. FDTI employs the perturbation formalism to numerically compute the derivatives of the free energy function with respect to the coupling parameter. Using Eq. 119, it is possible to compute the change in free energy (A_i) for a perturbation away from the i^th point (_i ± ):

Eq. 120        

Computing A_i for many different values of spanning the interval from 0 to 1, dividing each A_i by , and then numerically integrating over the interval, the total free energy change A can be estimated. Mathematically, this is summarized as:

Eq. 121        

where k is the number quadrature points in the numerical integration. Notice that A_i \xda can be computed both for a forward (+) and a backward (-) perturbation. Computing both at the same time takes no more computer time and is a measure of convergence (both should be equal for suitable values of ). The Discover program averages the two values for the instantaneous value of A \xda at _i.

Method used in FDiscover

The value of _i in Eq. 121 for a given i depends on the numerical integration scheme used. For example, a simple trapezoidal rule approach would make each _i equal. More sophisticated integration methods may allow each _i to be different. The FDiscover program uses a Gaussian-Legendre quadrature method (Press 1986) that chooses the values needed, given the total number of intervals specified. Hence, you never explicitly specify the intervals.

Advantages

The advantages of FDTI are:

• Unlike the perturbation method (PM), large changes in free energy can be calculated in fewer steps.

• Unlike thermodynamic integration (TI), analytical derivatives of the Hamiltonian with respect to the coupling parameter are not needed.

• As has been suggested by Mezei, FDTI may converge faster than either TI or PM.

Incorporation of the coupling parameter

The first step in any free energy calculation is parameterizing the energy function to provide a continuous change between the thermodynamic states that are being compared. The relative free energy function as implemented in the FDiscover program is parameterized at the level of the forcefield parameters themselves. For example, the energy E of a harmonic bond as a function of the bond length b is:

Eq. 122        

Eq. 123        

As changes from 0 to 1, the bond description gradually changes from a State A bond to a State B bond.

Relative free energy--methodology

Background of example problem

A seminal free energy calculation reported by Jorgensen and Ravimohan (1985) showed that free energy perturbation techniques seem to accurately predict the difference in solvation energy between methanol and ethane. This calculation has since been repeated in several laboratories, using independent programs and forcefields, and has become a de facto benchmark for free energy calculations (Singh et al. 1987; Fleischman and Brooks 1987).

Construction of the thermodynamic cycle

Consider the following thermodynamic cycle for the solvation of ethane and methanol:

This cycle requires that:

Eq. 124        

If there is no difference in intramolecular energy between ethane and methanol, then A₃ is zero and the gas phase calculation can be ignored. Although this is intuitively incorrect, a suitable forcefield could treat the intramolecular energy of ethane and methanol identically. For example, if the nonpolar hydrogens were not represented explicitly, no atom pairs would have more than two intervening bonds. Most forcefields calculate intramolecular nonbond interactions only for atom pairs separated by more than two intervening bonds. Furthermore, if the solutes are treated as rigid models, the energy would not be affected by deformations in bond lengths or angles.

In this exercise all hydrogens are explicitly included. In addition, variations in degrees of freedom are allowed for both the solvent and solute molecules.

Including all degrees of freedom requires both a vacuum and a solvent calculation. The computational penalty for an additional gas phase calculation is minor; the solvent calculation is two orders of magnitude more intensive. Moreover, because setting up the free energy calculation is virtually identical for both calculations, the gas phase calculation provides a tractable example that yields experience and confidence in performing free energy calculations without wasting an inordinate amount of computer time.

Defining the = 0 and = 1 states

To set up the calculation, you must first design the chemical perturbation for methanol going to ethane. You must first choose which of these models corresponds to the = 0 state. The = 0 state must always have enough atoms to accommodate both molecular extremes, even if invisible "dummy" atoms must be added as placeholders for atoms that appear as the perturbation proceeds.

Ethane is the logical choice for the = 0 state, because methanol can be completely subsumed within ethane. However, to illustrate that the choice is arbitrary, the calculation in the opposite direction, for methanol going to ethane, is described.

The = 0 state is provided to FDiscover as would any model for simulation, in ordinary coordinate (.car) and molecular data (.mdf) files. In the FDiscover command input file (.inp), you can specify which atoms are being perturbed and their limiting states with the Warp command. The limiting states are defined completely by the potential atom type, the partial atomic charge, and the mass of each atom.

Analyzing the output of a methanol-to-ethane free energy calculation

The free energy summary table printed out at the end of the completed free energy loop (after the last resume) contains the complete history. Throughout the run, a partial summary table is output after each resume and contains all the available data.

Assessing convergence

Graphing the Boltzmann factor [exp -(E \xda kT)] vs. time aids in assessing convergence during a given calculation. Systematic drift in exp -(E \xda kT) over time is symptomatic of nonconvergence.

To simplify analysis, various quantities are output during free energy calculations to a file with the extension .tot. Figure 32

Figure 32 . The ensemble-averaged quantity from Eq. 120 plotted vs. time for three values of

Each point plotted is actually an average of 10 points sampled over 100 fs.  

plots Column 5 of the .tot file [exp -(E \xda kT)] for the first, third, and last points for the methanol-to-ethane gas phase calculation. Note that there is a slight upward drift in the curve for the first point, although the others seem to have converged.

Absolute free energy

Theory and implementation

Example: Fentanyl

Analysis of results

This approach is a special case of the thermodynamic integration approach to free energy calculations, which is itself a general method for computing the change in free energy upon going from one thermodynamic state to another. Absolute free energy simply constrains one of these states to be a model system for which the absolute free energy is known analytically. By integrating from a known, albeit model, state to the final real state, the absolute free energy becomes the sum of the numerically computed thermodynamic integration step and the analytical absolute free energy of the model state.

Uses of absolute free energy calculations

Evaluating absolute free energies for particular conformations is an important goal for several reasons:

• Absolute free energy values of different thermodynamic states can be compared directly without having to devise a transition pathway between them, as is necessary in relative free energy calculations.

• The poor convergence properties associated with reversible structural changes in relative free energy calculations can be avoided.

Derivation for ideal systems

The absolute free energy algorithm depends on defining a model for which the partition function in Eq. 114 can be derived analytically. The model implemented in the FDiscover program is an ideal solid. That is, the atoms in the system are constrained harmonically to a lattice (analogous to a solid) and do not interact with each other (analogous to the familiar ideal gas). The Hamiltonian for such a system is:

Eq. 125        

Eq. 126        

Eq. 127        

Thermodynamic integration--derivation for real systems

As shown above, integrating the partition function analytically may be possible for simple Hamiltonians. However, for more realistic systems that include many nonbond and bond interactions between atoms, an analytic solution is impossible. As in the relative free energy calculation (see FDTI), we use thermodynamic integration to determine the change in free energy:

Eq. 128        

Eq. 129        

Without defining explicitly how the Hamiltonian depends on , we can write the difference in the energies as:

Eq. 130        

These snapshots are typically generated from either a molecular dynamics or a Monte Carlo simulation. Thus, appropriately averaging the results of a molecular dynamics trajectory enables Eq. 130 to be evaluated. This is far easier than calculating the general partition function of Eq. 129.

Eq. 130 is perfectly general for any classical system and is the fundamental equation of all thermodynamic integration methods. However, we have not yet described what the Hamiltonian is or how it can be parameterized in terms of . The Hamiltonian could be parameterized in an infinite number of ways. Because a function of state (the free energy) is being integrated, the path between the end points is arbitrary. For absolute free energy calculations, it turns out that the simplest form is perfectly adequate.

What is the Hamiltonian?

The parameterization scheme adopted by the FDiscover program is straightforward. The Hamiltonian is defined as a linear combination of the two potential energy functions that describe the extreme states:

Eq. 131        

Here, K(p) is the kinetic energy (because the two states are completely defined by their respective potential energy functions, the kinetic energy does not have to be coupled to ), V₀ is the normal potential energy function (including bonds, angles, torsions, etc.), and V_H is a harmonic oscillator site potential given by:

Eq. 132        

where K_i is the i^th atom's spring constant, r_i is the i^th atom's instantaneous coordinate, and r_i⁰ the reference lattice of the noninteracting atoms (often referred to as an Einstein solid). This choice of parameterization is motivated by the template-forcing concept in the FDiscover program and by previous free energy evaluations (Hoover and Ree 1967). In Eq. 131, is greater than 0 and less than 1, and describes an energy-space path between a system described by an unadulterated molecular forcefield ( = 0) and one represented by independent harmonic oscillators ( = 1).

Eq. 132 is thus the reference system for which an absolute free energy can be directly calculated. The potential V_H restricts the exploration of phase space to a region defined by atomic mean-squared displacements relative to the Einstein solid. Consequently, the choice of both the reference state and spring constants, in general, affects the calculated free energies. For most cases, the energy-minimized coordinates of a mechanically stable structure are used for the reference Einstein solid. Using Eq. 127 for the free energy of the Einstein solid A₁ and Eqs. 130 through 132, the absolute free energy of the real state is:

Absolute free energy of the real state

Eq. 133        

The simple form of Eq. 133 makes evaluating its ensemble average straightforward. A comprehensive dynamics trajectory for a given value of represents an ensemble of structures for that . By computing and averaging the values of (V₀ - V_H) for each structure in this ensemble, the integrand can be numerically estimated for any given value of . By performing several such calculations for many values of between 0 and 1, eventually the function can be numerically integrated.

Several practical simulation considerations are apparent from Eqs. 131 and 133.

For example, when is close to zero, the calculated ensemble of structures includes configurations far from the reference state (the ensemble is generated according to Eq. 131). Therefore, fluctuations of V_H are large when V_H is evaluated for this ensemble.

Conversely, when is close to one, the structures generated are minimally influenced by the real forcefield. These ensembles can include structures with distorted bonds, angles, and even overlapping atoms, leading to large fluctuations in V₀ in Eq. 133. Large fluctuations of (V₀ - V_H) decrease the precision with which its integral can be calculated. The overall effect is to destabilize the integral and, in extreme cases, cause it to diverge.

These undesirable effects can be minimized by reasonable choices of reference structures, spring constants, and integration algorithm. For example, the FDiscover program uses a Gaussian-Legendre quadrature algorithm to integrate Eq. 133. This has the useful property that one need not evaluate the function at the boundaries (where the divergence is worst).

Choice of the reference state

The choice of the reference state is the most critical step in an absolute free energy calculation. The reference state determines where the sampling of configuration space is centered. Choosing reference states that are only slightly different (for example, 0.1 Å rms in coordinates) can change the free energy significantly if the reference state has any residual strain.

Best results are obtained when the reference state represents a minimum-energy conformation on the normal energy surface. If the reference structure is a minimum for both V₀ and V_H, there is no impetus for the conformation to wander away, which, as described above, is the primary cause of divergence.

For non-minimum configurations such as the saddle point c in Figure 34, excess strain can be removed by minimizing with torsional restraints (see Example: Fentanyl).

Limitations

The absolute free energy technique is primarily used to evaluate the free energy of different conformations of the same model. This is particularly difficult to do by perturbation methods alone, since a path from one conformation to another may be difficult to construct and model (for example, a path converting an -helix into a -sheet). Although free energy methods are path independent (because what is being evaluated is the difference between thermodynamic state functions), the integration through the path must be performed reversibly, which is problematic for large structural rearrangements.

Example: Fentanyl

Figure 33 . Definitions of ₁ and ₂ for fentanyl

). Previous studies of fentanyl have examined the conformational behavior around the anilido nitrogen and identified two main classes of energy minima in the 2D energy map (Tolleaneare et al. 1986). Figure 34

Figure 34 . Two-Dimensional energy map of fentanyl as a function of the two dihedral angles defined in Figure 33

Locations a and b correspond to the two unique minima, and c is at a transition state.  

shows such a map, generated by contouring the energies obtained from a flexible-geometry minimization study performed with the Discover program. These minima are separated by approximately 10 kcal mol^-1 energy barriers and differ from each other by about 1 kcal mol^-1. The free energy of the states represented by the minima and the barriers between them will be calculated by the absolute free energy method.

The reference state has two roles in the free energy calculation. A trivial contribution is the free energy it contributes as an ideal solid. This contribution is trivial because it does not depend on the conformation of the reference structure. It depends only on the masses of its atoms and the spring constants constraining each atom to the lattice (see Eq. 127).

The nontrivial role is that it determines which part of configurational space is sampled during dynamics. In this example, the free energies of states a, b, and c (Figure 34) of fentanyl are compared. The reference states for a and b correspond to completely minimized structures at these points. Point c in Figure 34 represents the free energy at a barrier. Obviously, it would be incorrect to minimize this reference structure--it would simply roll down into one of the adjacent wells. The lowest-energy structure at the saddle point can be constructed by forcing the dihedral angles to adopt the values of the saddle point (using the force torsion command in the standalone version of the FDiscover program, or the Constraint/TorsionForce command in the Insight version) while minimizing the rest of the model (this, of course, is how the energy map in Figure 34 was generated in the first place).

It is important to choose the lowest-energy structure consistent with a given conformational state, so that the free energy calculated does not contain excess enthalpy. Constraining the sampling of configurational space to a rather small region localized around the reference structure limits the distribution of energies that is sampled. Moving the reference structure to a slightly higher-energy conformation would sample different energies and change the total free energy.

Consider Figure 35, where the history of ₁ and ₂ is superimposed on the ₁ \xda ₂ energy map for one of the ensembles. Note how only the local region around the reference structure is sampled. While this may be unsettling (what is the meaning of a free energy that does not represent all of configurational space), remember that you are interested in the difference in free energy between closely related states (a, b, and c differ mainly by rotation of two dihedral angles). To achieve this resolution, the sampling must be restricted accordingly.

Figure 35 . Time history of the dihedral angles ₁ and ₂ superimposed on the 2D energy map of Figure 34

Two trajectories are plotted, one for a minimum corresponding to state a in Figure 34, the other for the transition state (c). Note that only a small fraction of the conformational space of the model needs to be sampled to estimate the free energy of a conformer. Each trajectory corresponds to the first value, i.e., when the potential function is most like the standard potential.  

Hints for setting up the absolute free energy calculation

In the Insight version of the FDiscover program, an absolute free energy calculation is set up using the Parameters/Absolute command. The system is first minimized, then the atoms are tethered to the minimized locations. The Parameters/Minimize and Parameters/Dynamics commands can be used to control the initial minimization, as well as the dynamics used to sample the conformational space.

In the standalone version, the reference state is defined in the FDiscover input file (.inp), with all atoms being tethered. Using the tether command allows you to use different reference structures easily. Once the tethered atom list is generated, the model could undergo dynamics to randomize the starting point just prior to starting the free energy calculation.

There is an important exception to this rule. Once a free energy calculation is initiated, the reference structure is saved into a coordinate file with the extension .cli. Thus, a permanent record is maintained of the reference state for each free energy calculation. This file is used by the FDiscover program to restart free energy calculations that may have been interrupted. The file format is identical to the .car and .cor files so that they can be used interchangeably.

One important property of the .cli file is that, if the file exists when the free energy calculation begins (for example, from an earlier or different run), the reference state is read from the file regardless of any prior tethered atom list generated. This does not mean that if a .cli file is provided, a tethered atom list command is not needed. Rather, if the .cli file is present, no matter what conformation is used for the tethered atom list command, the actual reference structure used is that in the .cli file. Be aware then, when doing new free energy calculations requiring different reference states that any existing .cli files that would override the internally generated reference state must be removed.

The magnitude of the spring constants in Eq. 132 indirectly affects the free energy calculated by controlling the volume of phase space being sampled during dynamics. Large spring constants bias the sampling to the immediate neighborhood of the reference structure. Smaller values allow the model to sample more space.

The spring constants also affect the error associated with a calculation. Ensembles that contain structures far from the reference state make greater contributions to the free energy, by virtue of the greater fluctuations in (V₀ - V_H) for this ensemble.

The FDiscover program provides two methods for assigning spring constants:

1. All atoms are assigned the same value. You may want to assign large K_i values when you are studying a mechanically unstable state, to constrain dynamics to phase space near the reference structure. In the FDiscover standalone version this choice is specified by including the keyword assigning in the free energy command along with the specification of a spring constant. In the Insight environment, the force constant is controlled by changing the Spring Constant parameter in the Parameters/Absolute command.

2. Spring constants can be estimated from the definition of mean-squared displacements under a harmonic potential. In this case, the following expression is used:

Eq. 134        

and the default is for each atom to receive a unique value of K_i.

Method 1 is useful for studying mechanically unstable states, for example, the extended state of decaglycine. Method 2 is more appropriate for studying mechanically stable states, for example, the energy-minimized -helix or hairpin structures of decaglycine. Typical values of the spring constants range from 1 to 50 kcal mol^-1 Å^-1.

If possible, you should let Method 2 choose the spring constants automatically. However, if the conformational state is not in a local minimum, you must explicitly specify a spring constant. An estimated value corresponding to a nearby minimum is a good starting point. The error associated with an explicit choice should be monitored, as well as the configurational space sampled for this calculation. Error analysis is discussed under Analysis of results.

Another choice that has to be made is how many intervals to use in the integration. This of course depends ultimately on the level of precision needed and the behavior of the integral and cannot be anticipated for all cases.

For small flexible models such as fentanyl, using a medium spring constant (50 kcal mol^-1 Å^-1), the use of 10 intervals allows a reasonably behaved integral to be estimated to with a statistical error of ± 0.5 kcal mol^-1.

For careful work, the behavior of the integral (that is, how (V₀ - V_H) varies as a function of ) must be examined. Rapid changes in (V₀ - V_H) are indicative of systematic errors. Some regions of the integral may need to be computed for longer times and with more intervals, to achieve the desired accuracy.

The number of intervals typically ranges from 4 to 12. If adequate precision cannot be achieved with 12 intervals, look for and correct systematic factors destabilizing the integral (i.e., Are the spring constants too weak? Is the reference state at too high an energy? Is the system completely equilibrated?).

Analysis of results

Although a free energy calculation is done as a single FDiscover run, the calculation actually consists of several independent cycles of dynamics equilibration (initialize dynamics) and data collection (resume dynamics)--as many cycles as intervals. The results of these dynamics calculations are stored independently for each value, in three files.

The standard output file (extension .out) is the most important. Among other things, it includes a history of average energy values (potential and kinetic) during each dynamics cycle. Data needed for the actual free energy calculation are collected only during the resume dynamics phase. The actual free energy results are in the SUMMARY OF FREE-ENERGY CALCULATION table included at the end of each resume dynamics output. The total free energy can be calculated only at the conclusion of all the cycles. However, the incomplete table is output at the end of each interval as an intermediate result.

A second output file with the extension .tot contains only the instantaneous values of the kinetic energy plus the parameterized potential energy for the entire run (all values of ). The frequency of output into this file is controlled by the sampling option in the absolute free energy command and is, by default, output every step. The .tot file is particularly convenient as input to plotting and statistical analysis programs (such as SAS or RS1).

The third source of analysis information is the dynamics history file (extension .his). This file contains all the coordinates and velocities for each run. Because a free energy calculation consists of several cycles, each with its own initialize command, there can be several history files. To allow for several history files with the same root name, an option has been added to the initialize command to change the extension from .his to a number (e.g., FENTANYL.1, FENTANYL.2, etc.). This is specified by including the keyword save in the initialize command. If save is not included, the Discover program keeps only the history file corresponding to the last run.

Assessing and minimizing statistical and systematic errors

A major source of systematic error in these calculations is lack of convergence (that is, failure to equilibrate long enough to achieve thermodynamic equilibrium at each value) and insufficient sampling of configurational space. Other sources of systematic error include inaccuracies in the force field (both in the functional form and the parameters) and quantum mechanical effects.

Random errors are a natural consequence of free energy calculations. The statistical distribution of states available to a molecule at a given temperature is precisely what defines its entropy. Measuring entropy is an inherently statistical process that can be quantified with standard random error analysis procedures. Statistical errors are calculated for the average of (V₀ - V_H) at each value.

Systematic error/convergence

Careful analysis of the random errors reported by the FDiscover program helps point to likely systematic errors. For example, rapid changes in error as a function of could be caused by a failure to completely equilibrate.

Although repeating the entire run with longer equilibration and sampling times is an obvious solution, it is not a particularly efficient one. The integral computed under both short and long conditions may be indistinguishable up to = 0.8. It is wasteful to compute all intervals at the worst-case level.

For such situations, the absolute free energy command in the standalone version of the FDiscover program allows you to perform a free energy calculation over just a part of the interval. If you are using the Insight interface to the FDiscover program, the parameters Lower_lambda, Upper_lambda, and Quadrature Points in the Parameters/Absolute command control the points sampled, allowing you to break the calculations into ranges.

You could perform two independent free energy calculations. The first calculation would use short equilibration and sampling times to compute the free energy from = 0 to 0.8. The second calculation would compute the free energy from = 0.8 to 1.0 with longer equilibration and sampling times.

The total free energy is then just the sum of the two. Note that, when summing the free energies from these two runs, the ideal solid contribution must be included only once. Mathematically, the integral is simply being broken into pieces:

Eq. 135        

The .cli file contains the free energy comparison list, i.e., the template coordinates prevalent during tethering. This file can be used as input for further simulations with the same template coordinates, by specifying file filename with the calculate command in the standalone version of the FDiscover program.

The .tot file contains a history of instantaneous parameterized energies from the dynamics simulations.

Table of spring constants

A table of spring constants is printed before commencing the first dynamics simulation for the integration over . The data tabulated are absolute atom #, harmonic force constant, and (optionally) mean-square displacement and mean harmonic energies as deduced from spring-constant-generation procedures utilizing molecular dynamics.

Dynamics running averages

Dynamics running averages contain additional categories for (V₀ - V_H) and <V₀>.

Summary of averages per

Each table of dynamics averages generated by resume is followed by a table summarizing contributions to the free energy integral, including statistical errors and variances on the value dF for each lambda value.

The final page of the output contains a SUMMARY OF FREE ENERGY CALCULATIONS. It also includes a total free energy summary, giving the absolute free energy of the ideal solid (the = 1 state of the template). Each value has an associated output table in the .out file.