Production of the MD trajectory with an externally driven alchemical process

In a system of $N$ particles subject to a continuous alchemical transformations, only the non-bonded potential energy function is modified because of the presence of alchemical species. The full non bonded energy of the system is given by

$$V({\bf r}_1,.. {\bf r}_N, \lambda,\eta) = \sum_{ij} [1-\lambda_{ij}(t)]\frac{Q_i Q_j}{(\lambda_{ij}\gamma_{... ...\rm erfc}(\alpha r_{ij}) - \frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i$$

$$+ \frac{1}{2 \pi V} \sum_{\bf m \ne 0} \frac{ \exp(-{\bf m}^2/\alph... ... m^2} \sum_{ij} [1-\lambda_{ij}(t)] Q_i Q_j \exp(-i 2 \pi {\bf m} \cdot r_{ij})$$

$$+ 4\epsilon_{ij} [1-\eta_{ij}(t)] \left ( \frac{1} {\left [ \gamma_... ...frac{1} {\left [ \gamma\eta_{ij}(t) + (r_{ij}/\sigma{ij})^6 \right ] } \right )$$ (9.1)

where $V$ the unit cell volume, ${\bf m}$ a reciprocal lattice vector and $\alpha$ is the Ewald convergence parameter related to the width of the Gaussian spherical charge distribution. The first term in the non-bonded energy Eq. 9.1 is limited to the zero-cell and corresponds to the electrostatic interactions in the direct lattice; $\gamma_{\rm el}$ is a a soft-core parameter for the direct-lattice electrostatic term to avoid singularities when $\lambda_{ij}$ is approaching to 1.[] the second term refers to the self interactions of the Gaussian charge distributions and the third term corresponds to the interactions between Gaussian distributions in the zero cell as well as in the infinite direct lattice, reformulated as an absolutely convergent summation in the reciprocal lattice. The last term in Eq. 9.1, finally, corresponds to the modified atom-atom Van der Waals interaction introduced in Ref. [151] incorporating a soft-core parameterization[], where the infinity in the Lennard-Jones interaction is smoothed to zero as a function of the $\eta_i$. The parameter $\gamma$ is a positive constant (usually set[152] to 0.5) that controls the smoothing to zero of the derivatives Lennard Jones function as $r$ tends to zero.[153]

Table 9.1: Combination rules for alchemical and non alchemical species. The alchemical systems may contains three species: i) alchemical growing subsystems, ii) alchemical annihilation subsystems and iii) the non alchemical solvent. The $\lambda _i(t),\eta _i{i}(t)$ atomic factors within each of this species are all identical and equal to $\lambda _{G/A/S}(t),\eta _{G/A/S} (t)$ , where the index $G,A,S$ label the growing, annihilating and solvent species.

$i$	$j$	$\lambda_{ij}(t)$	$\eta_{ij}(t)$
Alchemical	Solvent	$\lambda_i(t)$	$\eta_i(t)$
Solvent	Alchemical	$\lambda_j(t)$	$\eta_j(t)$
Solvent	Solvent	0	0
Alchemical A	Alchemical A	0	0
Alchemical A	Alchemical B	1	1
Alchemical B	Alchemical A	1	1

In the present general formulation, according to Eq. 9.1, all atoms of the systems, whether alchemical or not, are characterized by an additional, time dependent and and externally driven ``coordinate'', the $\lambda_i(t)$ parameter controlling the charging/discharging of the system and the $\eta_i(t)$ parameter for switching on or off the atom-atom Lennard-Jones potential. The time dependence of the $\eta_i(t)$, $\lambda_i(t)$ atomic factors is externally imposed using an appropriately selected time protocol. The non bonded potential energy of Eq. 9.1 coincides with the standard potential energy of a system with no alchemical species when all the alchemical atomic factors $\lambda_i(t),~\eta_i(t)$ , referring to electrostatic and Van der Waals interactions, are constant and equal to zero. At the other extreme, when $\lambda_i(t)=\eta_i(t)=1$, the alchemical species disappears according to the ``mixing'' rules for $\lambda_{ij}(t), \eta_{ij}(t)$ factors specified in Table 9.1. These rules are such that the modified alchemical potential is enforced only when one of the two interacting atoms is alchemical while atom-atom interactions within a given alchemical species are accounted for with the standard potential or simply set to zero when they do refer to atoms on different alchemical species. In general, the time protocol for the $\lambda_i$, $\eta_i$ Van der Waals and electrostatic atomic parameters may differ from each other and for different alchemical species. A simple and sufficiently flexible scheme[154] would be that, for example, of allowing only two sets of alchemical species, i.e. the species to be annihilated and the species to be created, defining hence two different time protocols for the $\lambda_i$ and two more for the $\eta_i$ atomic parameters. Such a scheme allows, for example, the determination of the energy difference when one group in a molecule is replaced by an other group in a single alchemical simulation.

As remarked by others[152], it is convenient in a, e.g., alchemical creation, to switch on first the Van der Waals parameters changing $\eta$ for the alchemical atoms from one to zero and then charge the system varying $\lambda$ from one to zero. While for soft-core Lennard Jones term and the direct lattice electrostatic term the combination rules described in Table 1 can be straightforwardly implemented at a very limited computational cost in a standardly written force routine, the same rules cannot be directly applied to the reciprocal lattice part. In common implementation of the Ewald method, for obvious reason of computational convenience, the reciprocal lattice space double sum is rewritten in terms of a squared charge weighted structure factors as

$$V_{\rm rl} = {1 \over 2 \pi V} \sum_{{\bf m} \ne 0}^{\infty} {\ex... ...\bf m} \right \vert^2 } S \left ( {\bf m} \right ) S \left ( - {\bf m} \right )$$ (9.2)

In a system subject to a continuous alchemical transformation, the charge weighted structure factor becomes a function of the atomic factors $\lambda_i(t)$:

$$S\left ( {\bf m}, \lambda \right ) = \sum_{i}^{N} (1-\lambda_i(t)) Q_i \exp \left ( - 2 \pi i {\bf m} \cdot {\bf r}_i \right )$$ (9.3)

In the PME method, the sum of Eq. 9.3 is done via FFT by smearing the atomic charges on a regular grid in the direct lattice.[35] In this approach, all charge-charge interactions between alchemical and non alchemical species are almost inextricably mixed in the PME Ewald reciprocal lattice contribution and the application of the rules reported in Table 9.1 requires an extra effort indeed, an effort that has apparently deterred many to use the full Ewald method for computing the work done during continuous alchemical transformations. To this end and with no loss of generality, it is convenient to classify the system in an alchemical ``solute'' and in a non alchemical ``solvent'', with only the former being externally driven. We then label with $q(t)$ and $Q$ the time-dependent alchemical charges and the full time-invariant atomic charges of the solute, respectively, and with $Q_S$ the charges on the solvent. The alchemical $q(t)$ and full $Q$ solute charges are related by $q(t)=(1-\lambda(t)) Q$. When evaluating the reciprocal lattice energy via Eq. 9.2, the situation for the charge-charge electrostatic interactions is in represented in Table 9.2.

Table 9.2: Charge-charge interactions in alchemical transformations using the Ewald summation. The atomic charges labeled $q(t)$, $Q$ an d $Q_s$ refer to the alchemical charge, to the full (time-invariant) solute charge and to the solvent (non alchemical) charge.

Direct Lattice (Erfc)			Reciprocal Lattice (Erf)
$\frac{EQ} {r}$	$~~~~~\frac{q(t)Q_S}{r}$	$\frac{Q_SW_S}{r}$	$\frac{q(t)q(t)}{r}$	$~~~~~\frac{q(t)Q_S}{r}$	$\frac{Q_seq_S}{r}$
Only interactions $\ge$14			All interactions

In the direct lattice, the rules reported in table 9.1 can be implemented straightforwardly by excluding in the double atomic summation of Eq. 9.1 all the so-called 12 and 13 contacts. These atom-atom contacts involve directly bonded atoms of atoms bound to a common atom for which no electrostatic charge-charge contribution should be evaluated. In the reciprocal lattice, however, because of the structure of Eq. 9.2, all intra-solute interactions are implicitly of the kind $q_i(t)q_j(t) {\rm erf} (\alpha r)/r$ and 12 and 13 pairs are automatically considered in the sum of Eq. 9.2. The latter terms may be standardly removed in the zero cell by subtracting from the energy the quantity

$$V_{\rm intra} = \sum_{ij-\textrm{excl.}} q_{i}q_{j} { {\rm erf} \left ( \alpha r_{ij} \right ) \over r_{ij}}.$$ (9.4)

Regarding the 1-4 interactions, these are fully included in the reciprocal lattice sum, while in popular force fields only a portion of them is considered via the so-called fudge factors $f$. What must be subtracted in this case is the complementary interaction $q_i(t)q_j(t) (1-f) \rm erf (\alpha r)/r$.

It should be stressed here that, when the reciprocal lattice sum is computed using Eq. 9.2, the zero cell Erf contribution of the 12, 13 and 14$(1-f)$ interactions must be removed whether the two charges are alchemical or not. So, alchemically driven simulations imply no changes on the subtraction of these peculiar self-interactions with respect to a normally implemented program with no alchemical changes. The routines that implement Eq. 9.4 must be therefore called using the atomic charges $q_i = (1-\lambda_i(t)) Q_i$ whether alchemical or not (i.e. whether $\lambda_i$ is different from zero or not). With the same spirit, the self interaction in the zero cell, i.e. the term $-\frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i$ must be computed using the same charges.

We have seen in Table 9.2 that in the direct lattice the intrasolute non bonded electrostatic interactions are computed using the full time invariant solute charges $Q$, as alchemical changes affect only solute-solvent interaction energies. To recover the bare Coulomb potential for intrasolute interaction in a system subject to an alchemical transformation one must then subtract, as done for the 12 13 and 14(1-f) pairs, the Erf $q(t)q(t)$ contribution, and add a $QQ$ Erf term to the total energy of the system, producing the alchemical correction to the electrostatic energy

$$V_{\rm alch} = \sum_{ij > 14} Q_{i}Q_{j}[1 -(1-\lambda_i(t)) (1-\lambda_j(t))] { {\rm erf} \left ( \alpha r_{ij} \right ) \over r_{ij}}.$$ (9.5)

where the summation is extended to all non bonded intrasolute interactions. It should be stressed that the energy of Eq. 9.5 is a non trivial additive term that must be included in simulations of continuous alchemical transformations. Such term stems from the time dependent alchemical charges $q(t)$ and is due to the peculiar implementation of the Ewald method. $V_{\rm alch}$ is indeed a large contribution (10-15 kJ mol${-1}$ per solute atom) and its neglect may lead to severe errors in the electrostatic energies and to incorrect MD trajectories.

We can finally re-write down the total energy of a system subject to an alchemical transformation as

$$V(x,\lambda,\eta) = \sum_{ij} [1-\lambda_{ij}(t)]\frac{Q_i Q_j}{r_{ij}}{\rm erfc}(\al... ... intra} - \frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i^2 + V_{\rm alch}$$

$$+ 4\epsilon_{ij} [1-\eta_{ij}(t)] \left ( \frac{1} {\left [ \alpha\... ...frac{1} {\left [ \alpha\eta_{ij}(t) + (r_{ij}/\sigma{ij})^6 \right ] } \right )$$ (9.6)

where $V_{\rm rl}$, $V_{\rm intra}$, $V_{\rm alch}$ are defined in Eqs. 9.2, 9.4 and 9.5, respectively. All terms in Eq. 9.6, except for the self term $\frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i^2$, contribute to the atomic forces that can be standardly computed by taking the derivatives of the energy Eq. 9.6 with respect to the atomic position ${\bf r}_i$ producing the correct trajectories for alchemically driven systems under periodic boundary conditions and treated with the Ewald sum. In the Figure 9.1 we report the time record of the intra-solute electrostatic energy during the discharging of a molecule of ethanol in water in standard conditions. In spite of the huge changes in the contributing energy energy terms, the total intrasolute energy remains approximately constant during the transformation, modulated by the intramolecular motion, exactly as it should. The changes in the solute self term $-\frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i$ compensate, at all time steps, the variation of the direct lattice and of Erf intrasolute corrections. This balance does occur provided that all terms in the energy of Eq. 9.6 are accounted for, including the intrasolute alchemical Erf correction $V_{\rm alch}$ of Eq. 9.5.

Figure: Time record for the intrasolute energy arising form electrostatic interactions during the alchemical discharging of ethanol in water at T=300 K and P=1 Atm. The simulation went on for 15 ps. The red curve is due the self term $-\frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2Q_i^2$. The green curve is due to the direct lattice contribution. The magenta curve includes the terms $-V_{\rm intra}$ (Eq. 9.4) and $V_{\rm alch}$ (Eq. 9.5).

In a multiple time scheme, the individual contributions to the non bonded forces evolve in time with disparate time scales and must be hence partitioned in appropriately defined ``integration shell'' as described in details in Chapter 3. So in condensed phases, the direct lattice term is integrated in the fast short-ranged non bonded shell, while the reciprocal lattice summations (including the Erf intramolecular correction terms in $V_{\rm intra}$) are usually assigned, with an appropriate choice of the Gaussian parameter $\alpha$, to the intermediate non bonded shell. The Lennard-Jones term, finally, is split among the short-ranged, intermediate-range and long-range integration shells. The potential subdivision for condensed phases is basically unaffected by the implementation of alchemical transformation, except for the intrasolute self term $V_{\rm alch}$ and for the now time-dependent self term $\frac{\alpha}{\pi^{1/2}}\sum_i [1-\lambda_i(t)]^2q_i$. ^9.1. The latter can be safely included in the intermediate shell, while the former (a true direct lattice term) must be integrated in the sort-range shell. The $\lambda_i(t)$ and $\eta_i(t)$) factors, finally, must be updated, according to the predefined time protocol, before the force computation of the fast short-ranged non bonded shell.