Correspondence between light-absorption spectrum and nonequilibrium work distribution as a mean to access free energy differences between electronic states.

The problem of recovering the free energy difference between two electronic states has been investigated by Frezzato [Chem. Phys. Lett. 533, 106 (2012)], exploring the equivalence between light-absorption spectra and work distribution, hence opening to the application of a spectroscopic version of the Jarzynski equality (JE) [Phys. Rev. Lett. 78, 2690 (1997)]. Here, assuming the validity of the time-dependent perturbation theory, we demonstrate that such equivalence does not lead to the known form of the JE. This is ascribed to the fact that light-absorption processes cannot be described as stochastic processes. To emphasize such an aspect, we devise a stochastic model for the UV-vis (ultraviolet and visible) absorption, suitable for determining the free energy difference between two generic quantum manifolds in a JE-like fashion. However, the model would require explicit knowledge of the transition dipole moments, which are in general not available. Nonetheless, we derive a spectroscopic version of the JE that allows us to recover the free energy difference between the ground and an excited electronic state when the latter state is the only one observed in the spectrum.

Binding Free Energies of Host-Guest Systems by Nonequilibrium Alchemical Simulations with Constrained Dynamics: Illustrative Calculations and Numerical Validation.

In the companion article (Giovannelli et al., 10.1021/acs.jctc.7b00594), we presented an alchemical approach, based on nonequilibrium molecular dynamics simulations, to compute absolute binding free energies of a generic host-guest system. Two alternative computational routes, called binded-domain and single-point alchemical-path schemes, have been proposed. This study is addressed to furnish numerical validation and illustrative examples of the above-mentioned alchemical schemes. Validation is provided by comparing binding free-energy data relative to two poses of a Zn(II)·anion complex with those recovered from an alternative approach, based on steered molecular dynamics simulations. We illustrate important technical and theoretical aspects for a good practice in applying both alchemical schemes, not only through the calculations on the Zn(II)·anion complex, but also estimating absolute binding free energies of 1:1 complexes of ß-cyclodextrin with aromatic compounds (benzene and naphthalene). Comparison with experimental data and previous molecular dynamics simulation studies further confirms the validity of the present nonequilibrium-alchemical methodology.

Binding Free Energies of Host-Guest Systems by Nonequilibrium Alchemical Simulations with Constrained Dynamics: Theoretical Framework.

The fast-switching decoupling method is a powerful nonequilibrium technique to compute absolute binding free energies of ligand-receptor complexes (Sandberg et al., J. Chem. Theory Comput. 2014, 11, 423-435). Inspired by the theory of noncovalent binding association of Gilson and co-workers (Biophys. J. 1997, 72, 1047-1069), we develop two approaches, termed binded-domain and single-point alchemical-path schemes (BiD-AP and SiP-AP), based on the possibility of performing alchemical trajectories during which the ligand is constrained to fixed positions relative to the receptor. The BiD-AP scheme exploits a recent generalization of nonequilibrium work theorems to estimate the free energy difference between the coupled and uncoupled states of the ligand-receptor complex. With respect to the fast-switching decoupling method without constraints, BiD-AP prevents the ligand from leaving the binding site, but still requires an estimate of the positional binding-site volume, which may not be a simple task. On the other side, the SiP-AP scheme allows avoidance of the calculation of the binding-site volume by introducing an additional equilibrium simulation of ligand and receptor in the bound state. In the companion article (DOI: 10.1021/acs.jctc.7b00595), we show that the extra computational effort required by SiP-AP leads to a significant improvement of accuracy in the free energy estimates.

Elastic Barrier Dynamical Freezing in Free Energy Calculations: A Way To Speed Up Nonequilibrium Molecular Dynamics Simulations by Orders of Magnitude.

An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called "elastic barrier dynamical freezing", according to which a stiff potential-energy barrier is enforced at the boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speedup exceeding a few orders of magnitude.

Computing Free Energy Differences of Configurational Basins.

A simulation-based approach is proposed to estimate free energy differences between configurational states A and B, defined in terms of collective coordinates of the molecular system. The computational protocol is organized into three stages that can be carried on simultaneously. Two of them consist of independent simulations aimed at sampling, in turn, A and B states. In order to limit the evolution of the system around A and B, biased sampling simulations such as umbrella sampling can be employed. These simulations allow us to estimate local configuration integrals associated with A and B, which can be viewed as vibrational contributions to the free energy. Free energy evaluation is completed by the linking-path stage, in which the potential of mean force difference is estimated between two arbitrary points of the configurational surface, located the first around A and the second around B. The linking path in the space of the collective coordinates is arbitrary and can be computed with any method, starting from adaptive biasing potential/force approaches to nonequilibrium techniques. As an illustrative example, we present the calculation of free energy differences between conformational states of the alanine dipeptide in the space of backbone dihedral angles. The basic advantage of this method, that we term "path-linked domains" scheme, is to prevent accurate calculation of the whole free energy hypersurface in the space of the collective coordinates, thus limiting the statistical sampling to a minimum. Path-linked domains schemes can be applied to a variety of biochemical processes, such as protein-ligand complexation or folding-unfolding interconversion.

Simulations in generalized ensembles through noninstantaneous switches.

Generalized-ensemble simulations, such as replica exchange and serial generalized-ensemble methods, are powerful simulation tools to enhance sampling of free energy landscapes in systems with high energy barriers. In these methods, sampling is enhanced through instantaneous transitions of replicas, i.e., copies of the system, between different ensembles characterized by some control parameter associated with thermodynamical variables (e.g., temperature or pressure) or collective mechanical variables (e.g., interatomic distances or torsional angles). An interesting evolution of these methodologies has been proposed by replacing the conventional instantaneous (trial) switches of replicas with noninstantaneous switches, realized by varying the control parameter in a finite time and accepting the final replica configuration with a Metropolis-like criterion based on the Crooks nonequilibrium work (CNW) theorem. Here we revise these techniques focusing on their correlation with the CNW theorem in the framework of Markovian processes. An outcome of this report is the derivation of the acceptance probability for noninstantaneous switches in serial generalized-ensemble simulations, where we show that explicit knowledge of the time dependence of the weight factors entering such simulations is not necessary. A generalized relationship of the CNW theorem is also provided in terms of the underlying equilibrium probability distribution at a fixed control parameter. Illustrative calculations on a toy model are performed with serial generalized-ensemble simulations, especially focusing on the different behavior of instantaneous and noninstantaneous replica transition schemes.

Annealed importance sampling with constant cooling rate.

Annealed importance sampling is a simulation method devised by Neal [Stat. Comput. 11, 125 (2001)] to assign weights to configurations generated by simulated annealing trajectories. In particular, the equilibrium average of a generic physical quantity can be computed by a weighted average exploiting weights and estimates of this quantity associated to the final configurations of the annealed trajectories. Here, we review annealed importance sampling from the perspective of nonequilibrium path-ensemble averages [G. E. Crooks, Phys. Rev. E 61, 2361 (2000)]. The equivalence of Neal's and Crooks' treatments highlights the generality of the method, which goes beyond the mere thermal-based protocols. Furthermore, we show that a temperature schedule based on a constant cooling rate outperforms stepwise cooling schedules and that, for a given elapsed computer time, performances of annealed importance sampling are, in general, improved by increasing the number of intermediate temperatures.

Combining path-breaking with bidirectional nonequilibrium simulations to improve efficiency in free energy calculations.

An important limitation of unidirectional nonequilibrium simulations is the amount of realizations of the process necessary to reach suitable convergence of free energy estimates via Jarzynski's relationship [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. To this regard, an improvement of the method has been achieved by means of path-breaking schemes [R. Chelli et al., J. Chem. Phys. 138, 214109 (2013)] based on stopping highly dissipative trajectories before their normal end, under the founded assumption that such trajectories contribute marginally to the work exponential averages. Here, we combine the path-breaking scheme, called probability threshold scheme, to bidirectional nonequilibrium methods for free energy calculations [G. E. Crooks, Phys. Rev. E 61, 2361 (2000); R. Chelli and P. Procacci, Phys. Chem. Chem. Phys. 11, 1152 (2009)]. The method is illustrated and tested on a benchmark system, i.e., the helix-coil transition of deca-alanine. By using path-breaking in our test system, the computer time needed to carry out a series of nonequilibrium trajectories can be reduced up to a factor 4, with marginal loss of accuracy in free energy estimates.

Convective Replica-Exchange in Ergodic Regimes.

In a recent article (J. Comput. Chem. 2013, 34, 132-140), convective replica-exchange (convective-RE) has been presented as an alternative to the standard even-odd transition scheme. Computations on systems of various complexity have shown that convective-RE may increase the number of replica round-trips in temperature space with respect to the standard exchange scheme, leading to a more effective sampling of energy basins. Moreover, it has been shown that the method may prevent the formation of bottlenecks in the diffusive walk of replicas through the space of temperature states. By using an ideal temperature-RE model and a classical harmonic-oscillator RE scheme, we study the performances of convective-RE when ergodicity is not broken and convergence of acceptance probabilities is attained. In this dynamic regime, the round-trip ratio between convective and standard-RE is at maximum ∼ 1.5, a value much smaller than that observed in nonergodic simulations. For large acceptance probabilities, the standard-RE outperforms convective-RE. Our observations suggest that convective-RE can safely be used in either ergodic or non-ergodic regimes; however, convective-RE is advantageous only when bottlenecks occur in the state-space diffusion of replicas, or when acceptance probabilities are globally low. We also show that decoupling of the state-space dynamics of the stick replica from the dynamics of the remaining replicas improves the efficiency of convective-RE at low acceptance probability regimes.

Nonequilibrium Candidate Monte Carlo Simulations with Configurational Freezing Schemes.

Nonequilibrium Candidate Monte Carlo simulation [Nilmeier et al., Proc. Natl. Acad. Sci. U.S.A. 2011, 108, E1009-E1018] is a tool devised to design Monte Carlo moves with high acceptance probabilities that connect uncorrelated configurations. Such moves are generated through nonequilibrium driven dynamics, producing candidate configurations accepted with a Monte Carlo-like criterion that preserves the equilibrium distribution. The probability of accepting a candidate configuration as the next sample in the Markov chain basically depends on the work performed on the system during the nonequilibrium trajectory and increases with decreasing such a work. It is thus strategically relevant to find ways of producing nonequilibrium moves with low work, namely moves where dissipation is as low as possible. This is the goal of our methodology, in which we combine Nonequilibrium Candidate Monte Carlo with Configurational Freezing schemes developed by Nicolini et al. (J. Chem. Theory Comput. 2011, 7, 582-593). The idea is to limit the configurational sampling to particles of a well-established region of the simulation sample, namely the region where dissipation occurs, while leaving fixed the other particles. This allows to make the system relaxation faster around the region perturbed by the finite-time switching move and hence to reduce the dissipated work, eventually enhancing the probability of accepting the generated move. Our combined approach enhances significantly configurational sampling, as shown by the case of a bistable dimer immersed in a dense fluid.

Path-breaking schemes for nonequilibrium free energy calculations.

We propose a path-breaking route to the enhancement of unidirectional nonequilibrium simulations for the calculation of free energy differences via Jarzynski's equality [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. One of the most important limitations of unidirectional nonequilibrium simulations is the amount of realizations necessary to reach suitable convergence of the work exponential average featuring the Jarzynski's relationship. In this respect, a significant improvement of the performances could be obtained by finding a way of stopping trajectories with negligible contribution to the work exponential average, before their normal end. This is achieved using path-breaking schemes which are essentially based on periodic checks of the work dissipated during the pulling trajectories. Such schemes can be based either on breaking trajectories whose dissipated work exceeds a given threshold or on breaking trajectories with a probability increasing with the dissipated work. In both cases, the computer time needed to carry out a series of nonequilibrium trajectories is reduced up to a factor ranging from 2 to more than 10, at least for the processes under consideration in the present study. The efficiency depends on several aspects, such as the type of process, the number of check-points along the pathway and the pulling rate as well. The method is illustrated through radically different processes, i.e., the helix-coil transition of deca-alanine and the pulling of the distance between two methane molecules in water solution.