Implementation in ORAC 

Steered molecular dynamics in ORAC is implemented by adding an external driving potential depending on user defined internal coordinates in the form of stretching, bending, torsions. The general form of the time dependent external potential that bring the system from an initial state at $t=0$ to a different final state $t=\tau$ is given by

$$V_{ext}(t) = \frac {1}{2} \left [ \sum_{i=1}^{\rm N_r} K_i (r_i - r_{i0}(t) )^... ...(t) )^2 + \sum_{i=1}^{\rm N_\theta} K_i (\theta_i - \theta_{i0}(t) )^2 \right ]$$ (8.13)

where $r_i$, $\alpha_i$ and $\theta_i$ represents the actual $i$-th stretching , bending and torsional driven coordinate defined by arbitrarily selecting in the corresponding input definition the involved atoms. So a driven torsion or a stretching may be defined using arbitrarily chosen atoms of the solute that are not connected by any real bond. $r_{i0}(t)$, $\alpha_{i0}(t)$ and $\theta_{i0}(t)$ are time dependent parameters that defines the non equilibrium trajectory in the space of the coordinates. In ORAC , each of these parameters, given the duration $\tau$ of the non equilibrium experiment, is varied at constant speed from an initial value at time $t=0$ defining the reactants, to a final value at time $t=\tau$ defining the products :

$$r_{i}(t) = r_{i0} + \frac{r_{i\tau}-r_{i0}}{\tau} t = r_{i0} + v_{ir} t$$

$$\alpha_{i}(t) = \alpha_{i0} + \frac{\alpha_{i\tau}-\alpha_{i0}}{\tau} t = \alpha_{i0} + v_{i\alpha} t$$

$$\theta_{i}(t) = \theta_{i0} + \frac{\theta_{i\tau}-\theta_{i0}}{\tau} t = \theta_{i0} + v_{i\theta} t$$ (8.14)

As all the steering velocities are constant during the experiments, the above equations define a line

$${\bf z}(t) =\{(r_1(t),r_2(t),...\alpha_{1(t)},\alpha_2(t),..,\theta_1(t),\theta_2(t)...\}$$ (8.15)

in a reaction coordinate space at $N_r + N_\alpha + N_\theta$ dimensions

The work done by the external potential, Eq. 8.13, in the time $\tau$ of the non equilibrium driven process along the coordinate ${\bf z}$ is calculated as

$$W_0^{\tau} = \int_0^{\tau} \left [ \sum_{i=1}^{\rm N_r} K_i (r_i - r_{i0}(t) )... ...m_{i=1}^{\rm N_\theta} K_i (\theta_i - \theta_{i0}(t) ) v_{i\theta} \right ] dt$$ (8.16)

The equilibrium distribution of the starting points for independent work measurements can be determined (either by a standard equilibrium molecular dynamics simulation or by some enhanced simulation technique) by constraining the system with the harmonic constraint

$$V_{ext}(0) = \frac {1}{2} \left [ \sum_{i=1}^{\rm N_r} K_i (r_i - r_{i0})^2 + ... ...pha_{i0})^2 + \sum_{i=1}^{\rm N_\theta} K_i (\theta_i - \theta_{i0})^2 \right ]$$ (8.17)

for the reactants' state and

$$V_{ext}(\tau) = \frac {1}{2} \left [ \sum_{i=1}^{\rm N_r} K_i (r_i - r_{i\tau})^2... ...\tau})^2 + \sum_{i=1}^{\rm N_\theta} K_i (\theta_i - \theta_{i\tau})^2 \right ]$$ (8.18)

for the products' state. Having produced the work in a series of bidirectional experiments, one can then either apply the Bennett formula. Eq. 8.6, to compute the free energy differences between the reactants and the products states, or, using the intermediate work values $W_0^{t}$, apply Eq. 8.11 or Eq. 8.12 to reconstruct the entire potential of mean force along the the mono-dimensional driven trajectory in a multidimensional reaction coordinate space defined in Eq. 8.14. In order to define a non necessarily linear trajectory in a multidimensional reaction coordinate space (e.g. a putative minimum free energy path), on must be able to assign to a each steered coordinate a different steering time protocol. This can be done in ORAC by providing an auxiliary file defining the path in coordinate space. The file has the general form shown in Table 8.3.

Table: General format of the file defining of an arbitrary time protocol for a curvilinear path in a reaction coordinates space at $N_r + N_\alpha + N_\theta$ dimensions in ORAC . For a generic coordinate $\zeta =r,\alpha ,\theta$, the steering velocity between times $t_k$ and $t_{k+1}$ is constant and equal to $v_\zeta (t_k) = (\zeta (t_{k+1}) - \zeta (t_k))/(t_{k+1}-t_{k})$

$t_1$	$r_1(t_1)$	...	$r_{N_r}(t_1)$	$\alpha_{1}(t_1)$	...	$\alpha_{N_{\alpha}}(t_1)$	$\theta_{1}(t_1)$	...	$\theta_{N_{\alpha}}(t_1)$
$t_2$	$r_1(t_2)$	...	$r_{N_r}(t_2)$	$\alpha_{1}(t_2)$	...	$\alpha_{N_{\alpha}}(t_2)$	$\theta_{1}(t_2)$	...	$\theta_{N_{\alpha}}(t_2)$
	...	...	...
$t_n$	$r_1(t_n)$	...	$r_{N_r}(t_n)$	$\alpha_{1}(t_n)$	...	$\alpha_{N_{\alpha}}(t_n)$	$\theta_{1}(t_n)$	...	$\theta_{N_{\alpha}}(t_n)$

The free energy or potential of mean force obtained with the described protocols are not depurated by the jacobian terms arising form the definition of the reaction coordinates. For example, the potential of mean force, calculated with Eq. 8.11 or Eq. 8.12 along a driven distance for a freely rotating object includes the additional contribution $J(t)= 2 k_b T \ln (r_t/r_0)$ arising from the fact that the configurational probability $P(r)$, for two non interacting particles grows with the square of the distance. Moreover the PMF calculated using the driving potential given in Eq. 8.13 are in principle affected by the so-called stiff spring approximation,[148] i.e. if the constant $K_r, K\alpha, K_\theta$ in Eq. 8.13 are not large enough, then one actually computes the free energy associated to the Hamiltonian ${\cal H} = H + V_{ext}({\bf z} - {\bf z}_t)$ rather than that associated to the Hamiltonian $H({\bf z}={\bf z}_t)$. However the impact of the strength of the force constant on the computed non equilibrium average, especially if the reaction coordinate is characterized by inherently slow dynamics and/or the underlying unbiased potential of mean force is much less stiffer than the harmonic driving potential, is generally rather small even at relatively low values of force constant. With this respect, it has been shown that[148]

$$\phi({\bf z}) = F({\bf z}) + \frac{1}{2k} F' ({\bf z}) - \frac{1}{2 \beta k} F''({\bf z}) + O(1/k^2)$$ (8.19)

where $\phi({\bf z})$ is PMF of the unbiased system with the Hamiltonian $H({\bf z})$, while $F({\bf z})$ is the PMF that is actually measured in the SMD experiments, i.e. that corresponding to the biased Hamiltonian ${\cal H} = H({\bf z}) + V_{ext}({\bf z})$. From Eq. 8.19, one sees that if the derivatives of $F$ are not too high or $k$ is chosen large enough, then one can safely assume that $\phi({\bf z})\simeq F({\bf z})$.

eq:intraq:intra1