Torsion angle dynamics - Revision history

CyanaAdmin: /* Forces = torques = –gradient of the target function */

2023-02-25T18:30:58Z

Forces = torques = –gradient of the target function

← Older revision		Revision as of 20:30, 25 February 2023
Line 51:		Line 51:
	Because the target function is not a physical potential, there is no unique “natural” way to define the energy and temperature scale for protein structure calculations.		Because the target function is not a physical potential, there is no unique “natural” way to define the energy and temperature scale for protein structure calculations.

	== ~~Forces~~ = torques = –gradient of the target function ==		== Generalized forces = torques = –gradient of the target function ==

	The torques about the rotatable bonds, i.e., the negative gradients of the potential energy or target function with respect to torsion angles, -∇''V''(''θ''), are calculated by a fast recursive algorithm (Abe et al., 1984). The gradient of the target function can be calculated efficiently because the target function is a sum of functions of individual interatomic distances and torsion angles. The partial derivative of the target function ''V'' with respect to a torsion angle ''θ<sub>k</sub>'' is given by		The torques about the rotatable bonds, i.e., the negative gradients of the potential energy or target function with respect to torsion angles, -∇''V''(''θ''), are calculated by a fast recursive algorithm (Abe et al., 1984). The gradient of the target function can be calculated efficiently because the target function is a sum of functions of individual interatomic distances and torsion angles. The partial derivative of the target function ''V'' with respect to a torsion angle ''θ<sub>k</sub>'' is given by

CyanaAdmin at 18:28, 25 February 2023

2023-02-25T18:28:07Z

CyanaAdmin: /* Time step */

2023-02-25T18:24:41Z

Time step

← Older revision		Revision as of 20:24, 25 February 2023
Line 179:		Line 179:
	:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.		:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.

	The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>~~\scriptstyle~~\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\scriptstyle\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.		The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\scriptstyle\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.

	=== Accuracy of acceleration calculation ===		=== Accuracy of acceleration calculation ===

Guentert: /* Potential energy = target function */

2009-02-09T13:58:35Z

Potential energy = target function

← Older revision		Revision as of 15:58, 9 February 2009
Line 47:		Line 47:
	== Potential energy = target function ==		== Potential energy = target function ==

	For torsion angle dynamics the CYANA [[target function]] (Güntert et al., 1991; Güntert et al., 1997), ''V'', has the role of the potential energy, ''E''<sub>pot</sub> = ''V'', and energies are reported in units of the target function ''V'', i.e. Å<sup>2</sup>, in the current version of CYANA. If energies and temperatures are to be reported in their standard units of kJ/mol and K, we assume that ''E''<sub>pot</sub> = ''w''<sub>0</sub>''V'', with an overall (arbitrary) scaling factor ''w''<sub>0</sub> = 10 kJ mol<sup>-1</sup>Å<sup>-2</sup>. Because the target function is not a physical potential, there is no unique “natural” way to define the energy and temperature scale for protein structure calculations.		For torsion angle dynamics the CYANA [[target function]] (Güntert et al., 1991; Güntert et al., 1997), ''V'', has the role of the potential energy, ''E''<sub>pot</sub> = ''V'', and energies are reported in units of the target function ''V'', i.e. Å<sup>2</sup>, in the current version of CYANA. If energies and temperatures are to be reported in their standard units of kJ/mol and K, we assume that ''E''<sub>pot</sub> = ''w''<sub>0</sub>''V'', with an overall (arbitrary) scaling factor ''w''<sub>0</sub> = 10 kJ mol<sup>-1</sup>Å<sup>-2</sup>.

			Because the target function is not a physical potential, there is no unique “natural” way to define the energy and temperature scale for protein structure calculations.

	== Forces = torques = –gradient of the target function ==		== Forces = torques = –gradient of the target function ==

Guentert: /* Leap-frog algorithm */

2009-02-09T13:57:14Z

Leap-frog algorithm

← Older revision		Revision as of 15:57, 9 February 2009
Line 125:		Line 125:

	:<math>\begin{align}		:<math>\begin{align}
	\dot q_i(t + \Delta t/2) &= \dot q_i(t - \Delta t/2) + \Delta t \ddot q_i + O(\Delta t^3)\\		\dot q_i(t + \Delta t/2) &= \dot q_i(t - \Delta t/2) + \Delta t \ddot q_i(t) + O(\Delta t^3)\\
	q_i(t + \Delta t) &= q_i(t) + \dot q_i(t + \Delta t/2) + O(\Delta t^3)		q_i(t + \Delta t) &= q_i(t) + \dot q_i(t + \Delta t/2) + O(\Delta t^3)
	\end{align}</math>		\end{align}</math>

	The degrees of freedom, ''q<sub>i</sub>'', are the Cartesian coordinates of the atoms in conventional molecular dynamics simulation, or the torsion angles in CYANA. The ''O''(Δ''t''<sup>3</sup>) terms indicate that the errors with respect to the exact solution incurred by the use of a finite time step Δ''t'' are proportional to Δ''t''<sup>3</sup>.		The degrees of freedom, ''q<sub>i</sub>'', are the Cartesian coordinates of the atoms in conventional molecular dynamics simulation, or the torsion angles in CYANA. The ''O''(Δ''t''<sup>3</sup>) terms indicate that the errors with respect to the exact solution incurred by the use of a finite time step Δ''t'' are proportional to Δ''t''<sup>3</sup>.

	=== Temperature control ===		=== Temperature control ===

Guentert: /* Equations of motion */

2009-02-09T13:56:21Z

Equations of motion

← Older revision		Revision as of 15:56, 9 February 2009
Line 74:		Line 74:
	[[Image:CPUtimeVsSize.jpg\|thumb\|300px\|Fig. 2. Dependence on the protein size of the CPU time for the calculation of one conformer using the CYANA torsion angle dynamics algorithm with 4000 time steps on two different computers (solid and dotted lines) (Güntert et al., 1997).]]		[[Image:CPUtimeVsSize.jpg\|thumb\|300px\|Fig. 2. Dependence on the protein size of the CPU time for the calculation of one conformer using the CYANA torsion angle dynamics algorithm with 4000 time steps on two different computers (solid and dotted lines) (Güntert et al., 1997).]]

	In the case of n torsion angles as degrees of freedom, the n × n mass matrix ''M''(''θ'') and the ''n''-dimensional vector <math>C(\theta,\dot\theta)</math>, can be calculated explicitly (Mazur & Abagyan, 1989; Mazur et al., 1991). To generate a trajectory this linear set of n equations would have to be solved in each time step for the torsional accelerations <math>\ddot\theta</math>, formally by		In the case of n torsion angles as degrees of freedom, the n × n mass matrix ''M''(''θ'') and the ''n''-dimensional vector <math>\scriptstyle C(\theta,\dot\theta)</math>, can be calculated explicitly (Mazur & Abagyan, 1989; Mazur et al., 1991). To generate a trajectory this linear set of n equations would have to be solved in each time step for the torsional accelerations <math>\scriptstyle\ddot\theta</math>, formally by

	:<math>\ddot\theta = M(\theta)^{-1}C(\theta,\dot\theta)</math>.		:<math>\ddot\theta = M(\theta)^{-1}C(\theta,\dot\theta)</math>.

Guentert: /* Acceleration calculation */

2009-02-09T13:30:07Z

Acceleration calculation

← Older revision		Revision as of 15:30, 9 February 2009
Line 179:		Line 179:
	The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\scriptstyle\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\scriptstyle\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.		The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\scriptstyle\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\scriptstyle\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.

	=== ~~Acceleration~~ calculation ===		=== Accuracy of acceleration calculation ===

	[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]		[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]

Guentert: /* Time step */

2009-02-09T13:29:30Z

Time step

← Older revision		Revision as of 15:29, 9 February 2009
Line 157:		Line 157:
	<math>\scriptstyle\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.		<math>\scriptstyle\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.

	(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>\scriptstyle\dot\theta_i(t - \Delta t'/2)</math> by <math>\scriptstyle\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>\scriptstyle\dot\theta_e(t)</math> by <math>\scriptstyle\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)		(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>\scriptstyle\dot\theta_i(t - \Delta t'/2)</math> by <math>\scriptstyle\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>\scriptstyle\dot\theta_e(t)</math> by <math>\scriptstyle\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>\scriptstyle\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)

	:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,		:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,

Guentert: /* Time step */

2009-02-09T13:28:45Z

Time step

← Older revision		Revision as of 15:28, 9 February 2009
Line 157:		Line 157:
	<math>\scriptstyle\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.		<math>\scriptstyle\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.

	(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>\dot\theta_i(t - \Delta t'/2)</math> by <math>\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>\dot\theta_e(t)</math> by <math>\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)		(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>\scriptstyle\dot\theta_i(t - \Delta t'/2)</math> by <math>\scriptstyle\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>\scriptstyle\dot\theta_e(t)</math> by <math>\scriptstyle\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)

	:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,		:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,
Line 163:		Line 163:
	where ''T''<sup>ref</sup> is the reference value of the temperature. The instantaneous temperature, ''T''(''t''), is given by ''T''(''t'') = 2 ''E''<sub>kin</sub>(''t'')/''nk''<sub>B</sub>, where ''E''<sub>kin</sub>(''t'') is the kinetic energy at time ''t'', ''n'' denotes the number of torsion angles and ''k''<sub>B</sub> is the Boltzmann constant. Temperature control by coupling to an external heat bath (and time-step control) can be turned off by setting ''τ'' = ∞.		where ''T''<sup>ref</sup> is the reference value of the temperature. The instantaneous temperature, ''T''(''t''), is given by ''T''(''t'') = 2 ''E''<sub>kin</sub>(''t'')/''nk''<sub>B</sub>, where ''E''<sub>kin</sub>(''t'') is the kinetic energy at time ''t'', ''n'' denotes the number of torsion angles and ''k''<sub>B</sub> is the Boltzmann constant. Temperature control by coupling to an external heat bath (and time-step control) can be turned off by setting ''τ'' = ∞.

	(5) Calculate the torsional accelerations, <math>\ddot\theta(t)=\ddot\theta(\theta(t),\dot\theta_e(t))</math>, using equations [3–6].		(5) Calculate the torsional accelerations, <math>\scriptstyle\ddot\theta(t)=\ddot\theta(\theta(t),\dot\theta_e(t))</math>, using equations [3–6].

	(6) Calculate the new velocities at half time-step,		(6) Calculate the new velocities at half time-step,
Line 177:		Line 177:
	:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.		:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.

	The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.		The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\scriptstyle\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\scriptstyle\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.

	=== Acceleration calculation ===		=== Acceleration calculation ===

Guentert: /* Acceleration calculation */

2009-02-09T13:26:49Z

Acceleration calculation

← Older revision		Revision as of 15:26, 9 February 2009
Line 182:		Line 182:
	[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]		[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]

	Unlike the situation in Cartesian dynamics, where the accelerations are a function of the positions only, the torsional accelerations also depend on the velocities. These, however, are only known at half time-steps, whereas the positions and accelerations are required at full time-steps. In step 7 the presently used algorithm employs linear extrapolation to obtain an estimate of the velocity after the full time-step, , which is used in the next integration step to calculate the torsional accelerations (step 5). They are therefore beset with an additional error that could be eliminated by iterating the steps 5–7. In general, no significant improvement can be observed after one iteration (Mathiowetz et al., 1994). Using instead of introduces an additional error of order O(~~t3~~) into the velocity calculation of the leap-frog algorithm in step 6. However, the intrinsic accuracy of the velocity step is also O(~~t3~~). Thus, using the estimated velocities, , does not change the order of accuracy of the integration algorithm. Each iteration of steps 5–7 requires the calculation of the torsional accelerations which, in general, takes as long as the calculation of the target function and its gradient. It is therefore more efficient to slightly decrease the time step, t, than to perform steps 5–7 twice in each integration step. The situation could be different if a full physical force field were used, since then the calculation of torsional accelerations would require only a negligible fraction of the total computation time.		Unlike the situation in Cartesian dynamics, where the accelerations are a function of the positions only, the torsional accelerations also depend on the velocities. These, however, are only known at half time-steps, whereas the positions and accelerations are required at full time-steps. In step 7 the presently used algorithm employs linear extrapolation to obtain an estimate of the velocity after the full time-step, <math>\scriptstyle\dot\delta_e(t)</math>, which is used in the next integration step to calculate the torsional accelerations (step 5). They are therefore beset with an additional error that could be eliminated by iterating the steps 5–7. In general, no significant improvement can be observed after one iteration (Mathiowetz et al., 1994). Using <math>\scriptstyle\dot\theta_e(t)</math> instead of <math>\scriptstyle\dot\theta(t)</math> introduces an additional error of order ''O''(Δ''t''<sup>3</sup>) into the velocity calculation of the leap-frog algorithm in step 6. However, the intrinsic accuracy of the velocity step is also ''O''(Δ''t''<sup>3</sup>). Thus, using the estimated velocities, <math>\scriptstyle\dot\delta_e(t)</math>, does not change the order of accuracy of the integration algorithm. Each iteration of steps 5–7 requires the calculation of the torsional accelerations which, in general, takes as long as the calculation of the target function and its gradient. It is therefore more efficient to slightly decrease the time step, Δ''t'', than to perform steps 5–7 twice in each integration step. The situation could be different if a full physical force field were used, since then the calculation of torsional accelerations would require only a negligible fraction of the total computation time.

	Since in structure calculations with torsion angle dynamics the time-steps are made as long as possible a safeguard against occasional strong violations of energy conservation was introduced: If the total energy, E, has changed by more than 10% in a single time-step, this time-step is rejected and replaced by two time-steps of half length.		Since in structure calculations with torsion angle dynamics the time-steps are made as long as possible a safeguard against occasional strong violations of energy conservation was introduced: If the total energy, ''E'', has changed by more than 10% in a single time-step, this time-step is rejected and replaced by two time-steps of half length.

	== References ==		== References ==

← Older revision		Revision as of 20:28, 25 February 2023
Line 29:		Line 29:
	[[Image:TreeStructureRigidBodies.jpg\|thumb\|500px\|Fig. 1. (a) Tree structure of torsion angles for the tripeptide Val-Ser-Ile. Circles represent rigid units. Rotatable bonds are indicated by arrows that point towards the part of the structure that is rotated if the corresponding torsion angle is changed. (b) Excerpt from the tree structure formed by the torsion angles of a molecule, and definition of quantities required by the CYANA torsion angle dynamics algorithm.]]		[[Image:TreeStructureRigidBodies.jpg\|thumb\|500px\|Fig. 1. (a) Tree structure of torsion angles for the tripeptide Val-Ser-Ile. Circles represent rigid units. Rotatable bonds are indicated by arrows that point towards the part of the structure that is rotated if the corresponding torsion angle is changed. (b) Excerpt from the tree structure formed by the torsion angles of a molecule, and definition of quantities required by the CYANA torsion angle dynamics algorithm.]]

	The following quantities are defined for each cluster ''k'' (see Fig. 1): the “reference point”, ''r<sub>k</sub>'', which is the position vector of the end point of the bond between the clusters ''p''(''k'') and ''k''; <math>~~\scriptstyle~~ v_k = \dot{r}_k</math>, the velocity of the reference point; ''ω<sub>k</sub>'', the angular velocity vector of the cluster; ''Y<sub>k</sub>'', the vector from the reference point to the center of mass of the cluster; ''m<sub>k</sub>'', the mass of the cluster ''k''; ''I<sub>k</sub>'', the inertia tensor of the cluster ''k'' with respect to the reference point, given by		The following quantities are defined for each cluster ''k'' (see Fig. 1): the “reference point”, ''r<sub>k</sub>'', which is the position vector of the end point of the bond between the clusters ''p''(''k'') and ''k''; <math>v_k = \dot{r}_k</math>, the velocity of the reference point; ''ω<sub>k</sub>'', the angular velocity vector of the cluster; ''Y<sub>k</sub>'', the vector from the reference point to the center of mass of the cluster; ''m<sub>k</sub>'', the mass of the cluster ''k''; ''I<sub>k</sub>'', the inertia tensor of the cluster ''k'' with respect to the reference point, given by

	:<math>\textstyle I_k = \sum_\alpha m_\alpha I(y_\alpha)</math>		:<math>\textstyle I_k = \sum_\alpha m_\alpha I(y_\alpha)</math>
Line 76:		Line 76:
	[[Image:CPUtimeVsSize.jpg\|thumb\|300px\|Fig. 2. Dependence on the protein size of the CPU time for the calculation of one conformer using the CYANA torsion angle dynamics algorithm with 4000 time steps on two different computers (solid and dotted lines) (Güntert et al., 1997).]]		[[Image:CPUtimeVsSize.jpg\|thumb\|300px\|Fig. 2. Dependence on the protein size of the CPU time for the calculation of one conformer using the CYANA torsion angle dynamics algorithm with 4000 time steps on two different computers (solid and dotted lines) (Güntert et al., 1997).]]

	In the case of n torsion angles as degrees of freedom, the n × n mass matrix ''M''(''θ'') and the ''n''-dimensional vector <math>~~\scriptstyle~~ C(\theta,\dot\theta)</math>, can be calculated explicitly (Mazur & Abagyan, 1989; Mazur et al., 1991). To generate a trajectory this linear set of n equations would have to be solved in each time step for the torsional accelerations <math>~~\scriptstyle~~\ddot\theta</math>, formally by		In the case of n torsion angles as degrees of freedom, the n × n mass matrix ''M''(''θ'') and the ''n''-dimensional vector <math>C(\theta,\dot\theta)</math>, can be calculated explicitly (Mazur & Abagyan, 1989; Mazur et al., 1991). To generate a trajectory this linear set of n equations would have to be solved in each time step for the torsional accelerations <math>\ddot\theta</math>, formally by

	:<math>\ddot\theta = M(\theta)^{-1}C(\theta,\dot\theta)</math>.		:<math>\ddot\theta = M(\theta)^{-1}C(\theta,\dot\theta)</math>.
Line 124:		Line 124:
	=== Leap-frog algorithm ===		=== Leap-frog algorithm ===

	The integration scheme for the equations of motion in torsion angle dynamics (Mathiowetz et al., 1994) is a variant of the “leap-frog” algorithm used in Cartesian space molecular dynamics (Allen & Tildesley, 1987). To obtain a trajectory, the equations of motion are numerically integrated by advancing the ''i'' = 1,…,''n'' (generalized) coordinates ''q<sub>i</sub>'' and velocities <math>~~\scriptstyle~~\dot q_i</math> that describe the system by a small but finite time step Δ''t'':		The integration scheme for the equations of motion in torsion angle dynamics (Mathiowetz et al., 1994) is a variant of the “leap-frog” algorithm used in Cartesian space molecular dynamics (Allen & Tildesley, 1987). To obtain a trajectory, the equations of motion are numerically integrated by advancing the ''i'' = 1,…,''n'' (generalized) coordinates ''q<sub>i</sub>'' and velocities <math>\dot q_i</math> that describe the system by a small but finite time step Δ''t'':

	:<math>\begin{align}		:<math>\begin{align}
Line 157:		Line 157:
	:<math>\epsilon(t) = \left\|\frac{E(t)-E(t-\Delta t')}{E(t)}\right\|</math>.		:<math>\epsilon(t) = \left\|\frac{E(t)-E(t-\Delta t')}{E(t)}\right\|</math>.

	<math>~~\scriptstyle~~\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.		<math>\lambda_\epsilon^\max</math> is the maximally allowed value of the scaling factor, which is set to 1.025 in CYANA. The time constant, ''τ'' >> 1, is a user-defined parameter that is measured in units of the time-step. Typically, a value of ''τ'' = 20 is used for structure calculations with CYANA. To calculate ''ε''(''t''), the total energy ''E''(''t'') is evaluated before velocity scaling is applied (step 4 below), whereas for ''E''(''t'' - Δ''t''') the value after velocity scaling in the preceding time-step is used. Thus, the measurement of the accuracy of energy conservation is not affected by the scaling of velocities. An exact algorithm would yield ''E''(''t'') = ''E''(''t'' - Δ''t''') and consequently ''ε''(''t'') = 0.

	(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>~~\scriptstyle~~\dot\theta_i(t - \Delta t'/2)</math> by <math>~~\scriptstyle~~\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>~~\scriptstyle~~\dot\theta_e(t)</math> by <math>~~\scriptstyle~~\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>~~\scriptstyle~~\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)		(4) Adapt the temperature by scaling of the torsional velocities, i.e., replace <math>\dot\theta_i(t - \Delta t'/2)</math> by <math>\lambda_T\dot\theta_i(t - \Delta t'/2)</math> and <math>\dot\theta_e(t)</math> by <math>\lambda_T\dot\theta_e(t)</math> (see step 7 below for the definition of <math>\dot\theta_e(t)</math>). The velocity scaling factor, ''λ<sub>T</sub>'', is given by (Berendsen et al., 1984)

	:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,		:<math>\lambda_T = \sqrt{1+\frac{T^\text{ref}-T(t)}{\tau T(t)}}</math>,
Line 165:		Line 165:
	where ''T''<sup>ref</sup> is the reference value of the temperature. The instantaneous temperature, ''T''(''t''), is given by ''T''(''t'') = 2 ''E''<sub>kin</sub>(''t'')/''nk''<sub>B</sub>, where ''E''<sub>kin</sub>(''t'') is the kinetic energy at time ''t'', ''n'' denotes the number of torsion angles and ''k''<sub>B</sub> is the Boltzmann constant. Temperature control by coupling to an external heat bath (and time-step control) can be turned off by setting ''τ'' = ∞.		where ''T''<sup>ref</sup> is the reference value of the temperature. The instantaneous temperature, ''T''(''t''), is given by ''T''(''t'') = 2 ''E''<sub>kin</sub>(''t'')/''nk''<sub>B</sub>, where ''E''<sub>kin</sub>(''t'') is the kinetic energy at time ''t'', ''n'' denotes the number of torsion angles and ''k''<sub>B</sub> is the Boltzmann constant. Temperature control by coupling to an external heat bath (and time-step control) can be turned off by setting ''τ'' = ∞.

	(5) Calculate the torsional accelerations, <math>~~\scriptstyle~~\ddot\theta(t)=\ddot\theta(\theta(t),\dot\theta_e(t))</math>, using equations [3–6].		(5) Calculate the torsional accelerations, <math>\ddot\theta(t)=\ddot\theta(\theta(t),\dot\theta_e(t))</math>, using equations [3–6].

	(6) Calculate the new velocities at half time-step,		(6) Calculate the new velocities at half time-step,
Line 179:		Line 179:
	:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.		:<math>\theta(t + \Delta t) = \theta(t) + \Delta t \dot\theta(t+\Delta t/2)</math>.

	The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>~~\scriptstyle~~\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.		The algorithm is initialized by setting ''t'' = 0, Δ''t''' = Δ''t'' and <math>\dot\theta_e(0)=\dot\theta(-\Delta t/2)</math>. The initial torsional velocities, <math>\dot\theta(-\Delta t/2)</math>, are chosen randomly from a normal distribution with zero mean value and a standard deviation which ensures that the initial temperature has a predefined value, ''T''(0). Once the time step ''t'' → ''t'' + Δ''t'' is completed by going through the operations 1 to 8, ''t'' is replaced by ''t'' + Δ''t'', and Δ''t''' by Δ''t''.

	=== Accuracy of acceleration calculation ===		=== Accuracy of acceleration calculation ===
Line 185:		Line 185:
	[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]		[[Image:Schedule.jpg\|thumb\|300px\|Fig. 3. Plots versus the number of torsion angle dynamics steps for a structure calculation with simulated annealing for the protein cyclophilin A. (a) Temperature of the heat bath to which the system is weakly coupled. (b) Integration time step length, Δ''t''. (c) RMS torsion angle change, ''σ<sub>θ</sub>'', between successive time steps. Data are averaged over 50 time steps and 80 independently calculated conformers (Güntert et al., 1997).]]

	Unlike the situation in Cartesian dynamics, where the accelerations are a function of the positions only, the torsional accelerations also depend on the velocities. These, however, are only known at half time-steps, whereas the positions and accelerations are required at full time-steps. In step 7 the presently used algorithm employs linear extrapolation to obtain an estimate of the velocity after the full time-step, <math>~~\scriptstyle~~\dot\delta_e(t)</math>, which is used in the next integration step to calculate the torsional accelerations (step 5). They are therefore beset with an additional error that could be eliminated by iterating the steps 5–7. In general, no significant improvement can be observed after one iteration (Mathiowetz et al., 1994). Using <math>~~\scriptstyle~~\dot\theta_e(t)</math> instead of <math>~~\scriptstyle~~\dot\theta(t)</math> introduces an additional error of order ''O''(Δ''t''<sup>3</sup>) into the velocity calculation of the leap-frog algorithm in step 6. However, the intrinsic accuracy of the velocity step is also ''O''(Δ''t''<sup>3</sup>). Thus, using the estimated velocities, <math>~~\scriptstyle~~\dot\delta_e(t)</math>, does not change the order of accuracy of the integration algorithm. Each iteration of steps 5–7 requires the calculation of the torsional accelerations which, in general, takes as long as the calculation of the target function and its gradient. It is therefore more efficient to slightly decrease the time step, Δ''t'', than to perform steps 5–7 twice in each integration step. The situation could be different if a full physical force field were used, since then the calculation of torsional accelerations would require only a negligible fraction of the total computation time.		Unlike the situation in Cartesian dynamics, where the accelerations are a function of the positions only, the torsional accelerations also depend on the velocities. These, however, are only known at half time-steps, whereas the positions and accelerations are required at full time-steps. In step 7 the presently used algorithm employs linear extrapolation to obtain an estimate of the velocity after the full time-step, <math>\dot\delta_e(t)</math>, which is used in the next integration step to calculate the torsional accelerations (step 5). They are therefore beset with an additional error that could be eliminated by iterating the steps 5–7. In general, no significant improvement can be observed after one iteration (Mathiowetz et al., 1994). Using <math>\dot\theta_e(t)</math> instead of <math>\dot\theta(t)</math> introduces an additional error of order ''O''(Δ''t''<sup>3</sup>) into the velocity calculation of the leap-frog algorithm in step 6. However, the intrinsic accuracy of the velocity step is also ''O''(Δ''t''<sup>3</sup>). Thus, using the estimated velocities, <math>\dot\delta_e(t)</math>, does not change the order of accuracy of the integration algorithm. Each iteration of steps 5–7 requires the calculation of the torsional accelerations which, in general, takes as long as the calculation of the target function and its gradient. It is therefore more efficient to slightly decrease the time step, Δ''t'', than to perform steps 5–7 twice in each integration step. The situation could be different if a full physical force field were used, since then the calculation of torsional accelerations would require only a negligible fraction of the total computation time.

	Since in structure calculations with torsion angle dynamics the time-steps are made as long as possible a safeguard against occasional strong violations of energy conservation was introduced: If the total energy, ''E'', has changed by more than 10% in a single time-step, this time-step is rejected and replaced by two time-steps of half length.		Since in structure calculations with torsion angle dynamics the time-steps are made as long as possible a safeguard against occasional strong violations of energy conservation was introduced: If the total energy, ''E'', has changed by more than 10% in a single time-step, this time-step is rejected and replaced by two time-steps of half length.