3.4 Free Energy Calculations
3.4.4 Replica Exchange Methods
The presented methods to obtain potentials of mean force for a one- or mulitdimensional reaction coordinate enhance only the sampling of phase space along this collective vari-ables. For the remaining 3N − f (f being the dimension of the reaction coordinate) degrees of freedom, an exhaustive sampling of phase space cannot automatically be as-sumed. Considering e.g. the adsorption of a medium-sized molecule, such as a peptide on
Figure 3.4: Section of the trajectories of the four replicas along the solute temperatures from a metadynamics+REST simulation (a). Distribution of potential energies for solute temperatures of 300K (black), 350K (red) 400K (green) and 450K (blue) (b). For com-parison the inset displays the corresponding distributions when the same temperatures are applied to all degrees of freedom as in conventional RE. The simulations refer to the RKLPDA peptide on the oxidized titanium surface, cf. chap. 6
a solid substrate, where the center-of-mass distance to the surface represents the reaction coordinate one can obtain a potential of mean force along this collective variable using one of the methods described above. However, apart from barriers in this collective variable which are supposed to be handled by the chosen method of free energy calculation, other barriers, such as rotational barriers around peptide backbone torsion angles, exist As these are not directly associated with the reaction coordinate, their sampling is generally not taken into account by the primary sampling method in a controlled manner which can, in the worst case, exclude important molecular conformations from the potential of mean force calculation [133].
In order to overcome such problems, one can associate each degree of freedom, for which one expects significant barriers, with a reaction coordinate, increasing the dimensionality of the potential of mean force. Though excess degrees of freedom can be integrated out of the final free energy profile, not only the computational effort increases dramatically, also each individual barrier has to be identified separately which renders this method feasible only for small molecules. An alternative approach which takes care of sufficient sam-pling in a natural way without a bias potential, is provided by the replica exchange (RE)
technique [160]. In this approach possible free energy barriers are overcome at elevated temperatures, increasing the probability of barrier crossing according to the Arrhenius equation of activated processes P ∼ exp(−ΔG/kBT), where ΔG denotes the height of the barrier. As, on the other hand, one is usually interested in the properties at room temperature, several replicas of the same system are simulated in parallel, each of which is assigned a certain temperature, spanning the entire range between room temperature and maximum temperature. From time to time neighboring replicas are allowed to attempt an exchange of temperatures which causes a wandering of the replicas in temperature space.
Each attempt is accepted or rejected according to a standard Monte Carlo Metropolis criterion:
P(i↔j) = exp[−(βi−βj)(Ej −Ei)] (3.43) Ideally, this leads to each replica diffusing from the bottom to the top temperature, where barrier crossing is facilitated, and back to room temperature (as displayed in Fig. 3.4 (a)).
In this process the Boltzmann based acceptance criterion ensures a canonical distribution of configurations at each temperature. The base temperature trajectory can then be used to calculate various properties, such as free energies. While the replica exchange method was originally developed to merely enhance sampling of the phase space, it has been shown that its combination with other primary sampling techniques such as meta-dynamics can significantly improve the results, e.g. for the case of protein folding [30].
Regarding surface adsorption, a combination of RE with a fixed bias potential along the reaction coordinate as primary method has been reported [133].
One of the major drawbacks of replica exchange methods is that the potential energy distribution of neighboring replicas have to overlap sufficiently to ensure a reasonable acceptance ratio of the exchange attempts. As the relative width of the distributions decrease with increasing system size, one often has to employ a small spacing of tempera-tures which makes a large numbers of replicas necessary in order to reach sufficiently high temperatures for barrier crossing. This increases the computational effort in a dramatic, sometimes even prohibitive way. Different methods have been devised in order to increase the computational efficiency of RE simulations. As the difficulties arise mainly from sol-vent contributions to the potential energy which are usually not relevant compared to conformational transitions of the solute, different methods have been devised in order to reduce the influence of the heat capacity of the solvent. The temperature intervals with
global exchange of replicas (TIGER) method [111], for example provides an empirical ap-proach to decrease the number of replicas, while maintaining a good acceptance ratio for the replica exchanges.
Another technique which appears particularly promising, as it resembles a canonical distribution while avoiding a large number of replicas, is the replica exchange with solute tempering (REST) method [114]. Here the temperature ladder is only applied to the solute molecule, whereas the solvent temperature is kept constant. By doing so, the heat capacity of the system is artificially reduced, causing the potential energy distributions of neighboring replicas to overlap even for a larger temperature spacing. This means less replicas have to be included which increases computational efficiency and at the same time the diffusion from bottom requires less intermediate steps and therefore takes place much faster. Considering the example of a hexapeptide adsorbed on a surface, the po-tential energy distributions from conventional and solute tempering (ST) simulations are compared in Fig.3.4 (b). While the ST distributions exhibit sufficient overlap, the energy distributions including all degrees of freedom into the tempering do not overlap at all.
As temperature, or more specifically thermal energy, is only defined relative to the po-tential energy, rescaling of the popo-tential energy is equivalent to rescaling of tempera-ture. Accordingly, selective tempering is achieved by only rescaling the intramolecular interactions of the solute, while keeping the total kinetic energy constant. At elevated temperatures the potential energy surface of the solute appears more shallow, conse-quently intramolecular barriers appear smaller and can be crossed easier with a proba-bility P ∼ exp(−γΔE/(kBT) = exp(−ΔE/(kBT /γ) where γ < 1.0 is the scale factor.
Due the fictitious nature of the solute tempering, the coupling of the solvent, i.e. the magnitude of its interactions with the solute, is somewhat arbitrary in this framework.
Technical reasons suggest intermediate scaling factors. A factor of √
λ has appeared par-ticularly favorable regarding electrostatic interactions [173], i.e. the potential energy of replica i is
Ei = T0
TiE0mm + T0
TiE0ms+E0ss, (3.44) where m and s denote solute molecule respectively solvent, and the subscript 0 refers to the base temperature. In contrast to common replica exchange implementations, the Hamiltonian rather than the temperature is exchanged between the replicas.
Although the Metropolis criterion in principle yields an exact canonical ensemble, one has to keep in mind that enhanced sampling applies only to the degrees of freedom included
in the tempering. Barriers involving only solvent molecules, though, cannot be overcome by this method, which might, in some unfortunate cases, lead to a distribution of states deviating from the canonical ensemble.
Though promising and computationally feasible, this technique has only been applied to a small number of examples of molecules in solution so far [33, 85]. In this work we present for the first time the application of the metadynamics method combined with the REST technique to calculate adsorption free energies for medium-sized molecules with a considerable number of internal degrees of freedom.
Chapter 4
The Oxidized Titanium Surface
I
n this chapter I describe the first principle molecular dynamics (FPMD) simulations of the oxidation of the bare Ti(0001) surface, which were performed in order to construct the reference model for the oxidized titanium surface. I analyze the resulting structure, focusing in particular on the atomic charges. In section 4.2 the development of the classical force field is addressed, which is then applied in molecular dynamics simulations presented and discussed in section 4.3.This chapter has been published in Ref. [149].