拘束計算

　一般化座標を４ 章のようにユーザー定義して与えた際に、その座標を拘束した計算が一般的に可能である。

ネームリスト

　拘束計算に関するネームリストは、&constraint2, &addconstr である。その他にユーザー定義座標を与えるためのネームリスト &addCOM や &spcoord などが必要であるが、これは４-１ 節を参照のこと。 　ネームリストの仕様は、&constraint2 (３-１２ ) および &addconstr （３-１３ ）を参照すること。

　FreeFlex で拘束計算を行うにはネームリストファイルに、

=== input.nml =================================================================================
  .
  .
  .
&addCOM
  condition = "OW HW"
/

&spcoord
  condition = "Z 877 878"
/
!!! Na+ の重心が 877。878 は水全体の重心。

&addconstr
  condition = "SPCOORD 1 20.0d-10 1.0d-10"
/
================================================================================================

のように &addCOM, &spcoord, &addconstr を足してやればよい。上の２-５節のチュートリアル(sample/tutorial/05) を参照。この例では、重心877番と重心878番のz方向の距離を20 Åに固定されたトラジェクトリを得ることができる。

計算原理

　SHAKE-RATTLE法の考え方に従い、 k 番目の拘束だけを考えこれを満足させることを考える。 この時、時間発展として

${\begin{aligned}\mathbf {v} _{i}'(t+{\frac {1}{2}}\Delta t)&=\mathbf {v} _{i}(t)+{\frac {\Delta t}{2m_{i}}}\mathbf {F} _{i}(t)\\\mathbf {x} _{i}'(t+\Delta t)&=\mathbf {x} _{i}(t)+\Delta t\mathbf {v} _{i}'(t+{\frac {1}{2}}\Delta t)\\(*)\quad \mathbf {x} _{i}(t+\Delta t)&=\mathbf {x} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{x}(t)}{m_{i}}}{\frac {\partial g_{k}(\{\mathbf {x} _{i}(t)\})}{\partial \mathbf {x} _{i}}}\\(*)\quad \mathbf {v} _{i}(t+{\frac {1}{2}}\Delta t)&=\mathbf {v} _{i}'(t+{\frac {1}{2}}\Delta t)+{\frac {\gamma _{k}^{x}(t)}{m_{i}\Delta t}}{\frac {\partial g_{k}(\{\mathbf {x} _{i}(t)\})}{\partial \mathbf {x} _{i}}}\\\mathbf {v} _{i}'(t+\Delta t)&=\mathbf {v} _{i}(t+{\frac {1}{2}}\Delta t)+{\frac {\Delta t}{2m_{i}}}\mathbf {F} _{i}(t+\Delta t)\\(*)\quad \mathbf {v} _{i}(t+\Delta t)&=\mathbf {v} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{v}(t+\Delta t)}{m_{i}}}{\frac {\partial g_{k}(\{\mathbf {x} _{i}(t+\Delta t)\})}{\partial \mathbf {x} _{i}}}\end{aligned}}$

(1)

を使用する。$\gamma _{k}^{x},\gamma _{k}^{v}$ は、時間発展後の座標が拘束条件を満たすという条件

$(*)\quad g_{k}\left(\{\mathbf {x} _{i}(t)\},\{\mathbf {x} _{i}'(t+\Delta t)\},\gamma _{k}^{x}\right)=0$

(2)

$(*)\quad {\frac {dg_{k}}{dt}}\left(\{\mathbf {x} _{i}(t+\Delta t)\},\{\mathbf {v} _{i}'(t+\Delta t)\},\gamma _{k}^{v}\right)=0$

(3)

を解くことにより求める。 拘束が複数存在する場合には上記の手順を繰り返すことにより全ての拘束条件を満たす構造を求める。 収束判定は $(g_{k})^{2}<(tolerance)^{2},\left({\frac {dg_{k}}{dt}}\right)^{2}<\left({\frac {tolerance}{\Delta t}}\right)^{2}$

により行っている。拘束の種類によらず同じ tolerance を用いる為 $g_{k}$ は無次元に取る必要がある。

　実際の計算では、拘束条件によって値が変化する (*) の部分は constraint2_core_struct の拡張構造体が担い、その他の拘束条件の種類によらず形が変わらない部分は constraint2_struct が担う。

SHAKE法による補正前 $|g_{k}|$ が大きい場合

　初期構造が固定したい状態と離れている場合、拘束計算が収束しない可能性がある。 これに対処するため、SHAKE計算のループの前に$(g_{k})^{2}>(threshold)^{2}$ かどうかの判定を行い、$|g_{k}|$が threshold よりも大きい場合、拡張構造体内のchange_condition(g) サブルーチンによって拘束条件を一時的に変更する。 具体的な例として、$g=q(\{x_{i}\})-q_{0}$ の場合を考えると、$q_{0}$ を $q_{0}'=q-sign(g)*d$ (d は適当な定数） に変更することになる。 なお、変更時には警告メッセージを出力する。 SHAKE 計算終了後、拘束条件を元に戻すサブルーチン restore_condition() を実行し変更を元に戻す。

プログラムの構造と拘束条件の追加

　MD計算の際にどのような設定を行うかは５-４ 章及びチュートリアル２-５ を参照すること。 ここではソースコードの構造について述べる。

　拘束関連のプログラムの構造は以下のようになっている。

　拘束計算のコアとなる処理はconstraint2_core/ 以下のファイルで行っている。 このディレクトリ下にはconstraint_core.F に定義されているconstraint_core_struct の拡張構造体が置かれている。 この拡張構造体は、少なくとも <source lang="fortran">

   subroutine init(condition)
   !!!
   !!! 拘束条件を指定した文字列 condition から、
   !!! 必要なパラメータを取り出すサブルーチン
   !!!
     character(len=*), intent(in) :: condition
     !!!
     !!! 拘束条件を指定した文字列
     !!!
   endsubroutine

</source>

<source lang="fortran">

   subroutine SHAKE (coord,dt)
   !!!
   !!! 距離に関する拘束の式（テキスト式(2)）を解き、
   !!! coord に新しい座標を返すサブルーチン
   !!!
     type(coord_struct), intent(inout) :: coord
     !!!
     !!!  input: 拘束考慮前の座標
     !!!  output: 拘束考慮後の座標
     !!!
     !!!   REAL8, intent(in) :: dt
     !!!     時間刻み
     !!!
   endsubroutine

</source>

<source lang="fortran">

   subroutine RATTLE(coord)
   !!!
   !!! 速度に関する拘束の式(3)を解き、
   !!! coord に新しい速度を返すサブルーチン
   !!!
     type(coord_struct), intent(inout) :: coord
     !!!
     !!!  input: 拘束考慮前の座標
     !!!  output: 拘束考慮後の座標
     !!!
   endsubroutine

</source>

<source lang="fortran">

   REAL8 function g(coord)
   !!!
   !!! 拘束条件を規定する式 g を返す関数
   !!! ( g = 0 のとき拘束条件が満たされる）
   !!!
   !!!   type(coord_struct), intent(inout) :: coord
   !!!

</source>

<source lang="fortran">

   REAL8 function dgdt (coord)
   !!!
   !!! dg/dt を返す関数
   !!! ( dg/dt = 0 のとき拘束条件が満たされる）
   !!!
   !!!   type(coord_struct), intent(inout) :: coord
   !!!

</source>

<source lang="fortran">

   subroutine change_condition(g,coord,oldcoord)
   !!!
   !!! 系の構造が拘束条件と大きく離れている場合の処理
   !!!
   !!! REAL8,intent(in) :: g    !!! 拘束条件の式の値
   !!! type(coord_struct),intent(inout) :: coord,oldcoord
   !!!   時間発展前後の座標
   !!!

</source>

<source lang="fortran">

   subroutine restore_condition()
   !!!
   !!! change_condition で変更した値を元に戻すサブルーチン

</source> の7つのプロシージャを定義する必要がある。 その他追加のパラメータや変数が必要な場合は読み込み用のサブルーチンを適宜作成する。 なお、これらのプロシージャの実行によって座標や速度が更新される場合、その実行後に重心やユーザー定義座標を用いる場合には、それらの更新も別途行って同期する必要があることに注意すること。

　実際の追加には拡張構造体の作成のほか、constraint2_condition.F 中で分岐処理を追加する必要がある。select case 文に新しい case として、 <source lang="fortran">

     case("SAMPLE")
       allocate(sample)
       self%core => sample

</source> のようなものを追加するだけなのでソースを参考にしてほしい。

拘束条件の指定と実際の計算式

　MD計算で拘束条件を追加するには、FreeFlexに読み込ませるネームリストファイルに

&addconstr
  condition = “keyword options”
/

の形で新しい条件を指定すればよい。 また、&constraint2 ネームリストにより拘束計算全体に関わる設定も変更できる。 詳しくは tutorial.docx を見ること。

DISTANCE

　拘束条件を指定する文字列として

condition = “DISTANCE  i  j  r0  runit  maxmove”

とした場合、サイト-サイト間距離の拘束を行う。i, j はサイトの番号、r0 は拘束距離（m）を意味する。 また、runit, maxmoveは計算に用いるパラメータで下記の説明の通りである。 これらのパラメータは指定しなくても良く、その場合デフォルトの値として r^unit=1 Å,maxmove=0.01 Å が使用される。 拘束条件にDISTANCE を指定した際の具体的な拘束条件の式は

$g={\frac {(\mathbf {x} _{i}-\mathbf {x} _{j})^{2}-r_{0}^{2}}{(r^{unit})^{2}}}=0$

である。

　上記拘束条件のもとラグランジュ未定乗数を求めると、

${\begin{aligned}\gamma _{k}^{x}&={\frac {(r^{unit})^{2}\left[\left\{\mathbf {x} _{i}'(t+\Delta t)-\mathbf {x} _{j}'(t+\Delta t)\right\}^{2}-r_{0}^{2}\right]}{2\left({\frac {1}{m_{i}}}+{\frac {1}{m_{j}}}\right)\left\{\mathbf {x} _{i}(t)-\mathbf {x} _{j}(t)\right\}\cdot \left\{\mathbf {x} _{i}'(t+\Delta t)-\mathbf {x} _{j}'(t+\Delta t)\right\}}}\\\gamma _{k}^{v}&={\frac {(r^{unit})^{2}\left\{\mathbf {v} _{i}'(t+\Delta t)-\mathbf {v} _{j}'(t+\Delta t)\right\}\cdot \left\{\mathbf {x} _{i}(t+\Delta t)-\mathbf {x} _{j}(t+\Delta t)\right\}}{\left({\frac {1}{m_{i}}}+{\frac {1}{m_{j}}}\right)\left\{\mathbf {x} _{i}(t+\Delta t)-\mathbf {x} _{j}(t+\Delta t)\right\}^{2}}}\\\end{aligned}}$

となる。

　系の構造が拘束条件から大きく外れている場合、$r_{0}$ の代わりに

$r_{0}'={\sqrt {\left[\mathbf {x} _{i}(t)-\mathbf {x} _{j}(t)\right]^{2}}}-sign(g_{k})*maxmove$

を用いる。

SPCOORD

拘束条件を指定する文字列として

condition = “SPCOORD  n  nq  q0(1:nq)  qunit  maxmove”

とした場合、ユーザー定義座標の拘束を行う。 nはユーザー定義座標（ネームリストの上から1,2,…となる）の番号、nq は指定したユーザー定義座標の次元数 (※ nq は spcoord_core_struct の関数(get_nq)で設定されるので、ネームリストの指定は不要。)、 q0(1:nq) は拘束値、qunit, maxmove は計算に用いるパラメータである。具体的な拘束条件の式としては、

$g={\frac {q\left(\{\mathbf {x} _{i}\}\right)-q_{0}}{q^{unit}}}=0$

を用いる。 ここでユーザー定義座標を q で表記している。この時、

$\mathbf {x} _{i}(t+\Delta t)=\mathbf {x} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{q}}{m_{i}}}{\frac {1}{q^{unit}}}{\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}$

であるから、

${\begin{aligned}&{\frac {1}{q^{unit}}}\left[q\left(\left\{\mathbf {x} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{q}}{m_{i}}}{\frac {1}{q^{unit}}}{\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}\right\}\right)-q_{0}\right]\\\approx &{\frac {1}{q^{unit}}}\left[q'(t+\Delta t)+\sum _{i}{\frac {\gamma _{k}^{q}}{m_{i}}}{\frac {1}{q^{unit}}}{\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}\cdot {\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}-q_{0}\right]=0\end{aligned}}$

$\gamma _{k}^{q}=-q^{unit}\left[q(t+\Delta t)-q_{0}\right]/\sum _{i}{\frac {1}{m_{i}}}{\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}\cdot {\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}$

となる。 同様にして速度については、

${\frac {dg_{k}}{dt}}={\frac {1}{q^{unit}}}\sum _{i}{\frac {\partial q}{\partial \mathbf {x} _{i}}}\cdot \mathbf {v} _{i}=0$

が拘束条件となる。 従って速度の時間発展が、

$\mathbf {v} _{i}(t+\Delta t)=\mathbf {v} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{v}}{m_{i}}}{\frac {1}{q^{unit}}}{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}$

と表されることから、

${\frac {1}{q^{unit}}}\sum _{i}{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}\left[\mathbf {v} _{i}'(t+\Delta t)+{\frac {\gamma _{k}^{v}}{m_{i}}}{\frac {1}{q^{unit}}}{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}\right]=0$

$\gamma _{k}^{v}=-\xi ^{unit}\sum _{i}{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}\cdot \mathbf {v} _{i}'(t+\Delta t)/\sum _{i}{\frac {1}{m_{i}}}\left[{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}\right]^{2}$

が得られる。

　系の構造が拘束条件から大きく外れている場合、q_0 の代わりに

$q_{0}'=q(t)-sign(g_{k})*maxmove$

を用いる。

${\frac {\partial q}{\partial \mathbf {x} _{i}}}$ が不連続になる場合の処理

　water finger 座標等、時間発展やSHAKE法のループ中に注目する粒子の入れ替わりが起こる場合、${\frac {\partial q}{\partial \mathbf {x} _{i}}}$ が不連続になる。 ${\frac {\partial q(t)}{\partial \mathbf {x} _{i}}},{\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}$ の値として入れ替わり前後の値を用いてしまうと通常のSHAKE法のループが収束する保証はない。 そこでSHAKE計算においては粒子の入れ替わりを無視し、時刻 t+Δt における注目粒子を用いた場合の ${\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}$ に対応する値を使用することにする。 この処理は、spcoord_core_struct 内に注目粒子を交換せずに ${\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}$ を更新するサブルーチン update_q_constr(coord) を用意することで、

• constraint2_struct 内のサブルーチン save_oldcoordにおいて時刻 t における座標を oldcoord として保存。
• constraint2_core_struct 内のサブルーチン SHAKE(coord,oldcoord,dt) において、
  • coord%spcoord(n)%update_q_constr(coord) を実行。
  • coord%spcoord(n)%get_dqdx(dqdx) により${\frac {\partial q(t+\Delta t)}{\partial \mathbf {x} _{i}}}$ を取得。
  • coord%spcoord(n)%update_q_constr(oldcoord) を実行。
  • coord%spcoord(n)%get_dqdx(old_dqdx) により ${\frac {\partial q(t)}{\partial \mathbf {x} _{i}}}$ を取得。

のような処理により実装が可能となる。

FREEZE

• 古い仕様

　条件を

condition = “FREEZE  i”

とした場合、i 番目の重心グループの位置を固定する。 プログラム内では、位置・速度の補正時に

${\begin{aligned}\mathbf {x} _{i}(t+\Delta t)&=\mathbf {x} _{i}(t)\\\mathbf {v} _{i}(t+\Delta t)&=0\end{aligned}}$

としているだけである。

• 新しい仕様

$A$ 番目の重心グループの位置を空間固定系に対して固定する。 ただし単純に重心を空間固定系に対して拘束しようとすると1体の拘束力をかけることとなり、 内力場でなくなるため重心が移動する問題がある。 ここではこれを解決する方法として、$A$ の重心と系全体の重心間に拘束をかけることを考える。

$A$ 番目のグループの質量の和を $M_{A}$、系に含まれる全粒子の質量の和を $M$ とする。

${\begin{aligned}M_{A}=\sum _{i\in A}m_{i}\\M=\sum _{i}m_{i}\end{aligned}}$

系の相互作用が内力のみから成る場合、系の重心は空間固定系に対して移動しない。 従って $A$ の重心を空間固定系に対して固定したい場合、 系の重心に対する相対位置 $\Delta R$ を固定すればこれは達成される。 このとき拘束条件 $g$ およびその微分 $\partial g/\partial r_{i}$ は、次の式で書ける。

${\begin{aligned}g&={\frac {\sum _{i\in A}m_{i}r_{i}}{M_{A}}}-{\frac {\sum _{j}m_{j}r_{j}}{M}}-\Delta R\\&=\sum _{i\in A}m_{i}r_{i}{\Bigl (}{\frac {1}{M_{A}}}-{\frac {1}{M}}{\Bigr )}-\sum _{j\in {\bar {A}}}m_{j}r_{j}{\frac {1}{M}}-\Delta R=0\end{aligned}}$

${\begin{aligned}{\frac {\partial g}{\partial r_{i}}}={\begin{cases}m_{i}{\Bigl (}{\frac {1}{M_{A}}}-{\frac {1}{M}}{\Bigr )}&i\in A\\-{\frac {m_{i}}{M}}&i\in {\bar {A}}\\\end{cases}}\end{aligned}}$

拘束をかける前の座標および速度を $r_{i}^{0}$、$v_{i}^{0}$ とし、 拘束力による変位を座標および速度に対してそれぞれ $\delta r_{i}$、$\delta v_{i}$ とする。 ここで $r_{i}^{0}$ および $v_{i}^{0}$ は通常の力 $-\partial U/\partial r_{i}$ による変位を含むため、 必ずしも拘束条件は満たされていない。

${\begin{aligned}r_{i}&=r_{i}^{0}+\delta r_{i}\\v_{i}&=v_{i}^{0}+\delta v_{i}\end{aligned}}$

しかし通常の力が内力のみから構成されるとき、系の重心は初期配置での重心 $R_{0}$ に常に固定される。

${\begin{aligned}{\frac {\sum _{i}m_{i}r_{i}^{0}}{M}}=R_{0}\end{aligned}}$

拘束力による変位は未定乗数 $\lambda$ を用いて次式で書ける。

${\begin{aligned}\delta r_{i}=\lambda {\frac {\partial g}{\partial r_{i}}}{\frac {(\Delta t)^{2}}{2m_{i}}}={\begin{cases}\lambda {\Bigl (}{\frac {1}{M_{A}}}-{\frac {1}{M}}{\Bigr )}{\frac {(\Delta t)^{2}}{2}}&\equiv \delta r_{A}&i\in A\\-{\frac {\lambda }{M}}{\frac {(\Delta t)^{2}}{2}}&=\delta r_{A}{\frac {M_{A}}{M_{A}-M}}&i\in {\bar {A}}\\\end{cases}}\end{aligned}}$

このとき $\lambda$ は拘束条件 $g$ を満たすように決まるため、 $\lambda$ を求めるために解くべき方程式は次式となる。

${\begin{aligned}g&={\frac {\sum _{i\in A}m_{i}(r_{i}^{0}+\delta r_{i})}{M_{A}}}-{\frac {\sum _{j}m_{j}(r_{j}^{0}+\delta r_{j})}{M}}-\Delta R\\&={\frac {\sum _{i\in A}m_{i}r_{i}^{0}}{M_{A}}}+\delta r_{A}^{0}-R_{0}-\Delta R=0\end{aligned}}$

従って、グループ $A$ に属する粒子の拘束力による変位は

${\begin{aligned}\delta r_{A}&=\Delta R-{\frac {\sum _{i\in A}m_{i}r_{i}^{0}}{M_{A}}}+R_{0}\\&=R_{A}-{\frac {\sum _{i\in A}m_{i}r_{i}^{0}}{M_{A}}}\end{aligned}}$

となる。また、グループ $A$ に属さない粒子の変位はこれに $M_{A}/(M_{A}-M)$ をかけたものとなる。 ここまで議論した変位は MoveA における $t\rightarrow t+\Delta t$の時間発展である。 これに対応する速度の変化は $t\rightarrow t+\Delta t/2$ の時間発展であるから、

${\begin{aligned}\delta v_{i}={\frac {\delta r_{i}}{\Delta t}}\end{aligned}}$

となる。 また、$g$ の時間微分も拘束条件として課される。

${\begin{aligned}{\frac {dg}{dt}}={\frac {\sum _{i\in A}m_{i}v_{i}}{M_{A}}}-{\frac {\sum _{j}m_{j}v_{j}}{M}}=0\end{aligned}}$

拘束力による速度変化は、$t+\Delta t/2\rightarrow t+\Delta t$ における時間発展であるから、

${\begin{aligned}\delta v_{i}=\lambda {\frac {\partial g}{\partial r_{i}}}{\frac {\Delta t}{2m_{i}}}={\begin{cases}\lambda {\Bigl (}{\frac {1}{M_{A}}}-{\frac {1}{M}}{\Bigr )}{\frac {\Delta t}{2}}&\equiv \delta v_{A}&i\in A\\-{\frac {\lambda }{M}}{\frac {\Delta t}{2}}&=\delta v_{A}{\Bigl (}{\frac {M_{A}}{M_{A}-M}}{\Bigr )}&i\in {\bar {A}}\end{cases}}\end{aligned}}$

となる。ただし重心の速度は0である。

${\begin{aligned}{\frac {\sum _{i}m_{i}v_{i}^{0}}{M}}=0\end{aligned}}$

従って拘束力を求めるために解くべき方程式は

${\begin{aligned}{\frac {dg}{dt}}&={\frac {1}{M_{A}}}\sum _{i\in A}m_{i}(v_{i}^{0}+\delta v_{i})-{\frac {1}{M}}\sum _{j}m_{j}(v_{j}^{0}+\delta v_{j})\\&={\frac {\sum _{i\in A}m_{i}v_{i}^{0}}{M_{A}}}+\delta v_{A}=0\end{aligned}}$

であり、速度変化はこれを解いて次式となる。

${\begin{aligned}\delta v_{A}=-{\frac {\sum _{i\in A}m_{i}v_{i}^{0}}{M_{A}}}\end{aligned}}$

平均力と拘束力

結構ムズイのであきらめた。

Giovanni Ciccotti, Raymond Kapral, and Eric Vanden-Eijnden, "Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics," ChemPhysChem 6 (2005) 1809-1814.

Michiel Sprik and Giovanni Ciccotti, "Free energy from constrained molecular dynamics," J. Chem. Phys. 109 (1998) 7737-7744.

　FreeFlex version 1.0以降では、通常のRATTLE法に加え、質量のないサイトに対する拘束や特異点のある拘束条件に対する拘束アルゴリズム（SF-RATTLE）をサポートするよう改良がなされた。 　また、これに伴い通常のRATTLE法に対しても若干の修正が加えられており、これらについて説明する。

SF-RATTLE

FreeFlexでは一般化座標を用いて、自由な拘束条件を課すことができる。 ただし拘束条件の種類によってはRATTLE拘束に対して特異点となる場合があり、これを回避するために特別なアルゴリズム（SF-RATTLE）を 用いる必要がある。 例えば粒子1, 2, 3のなす角 $\theta$ を $\theta _{0}$ に固定するような拘束条件

${\begin{aligned}g=\cos {\theta }-\cos {\theta _{0}}=0\end{aligned}}$

に対して、$\cos {\theta }=0,\pi$ は g の微分が 0 になるために拘束力の向きが決まらず、特異点となっている。

質量のないサイトを含む拘束条件

TIP4P、TIP5PやCRK水モデルをはじめとして、ダミーサイトを含む分子モデルがしばしば用いられる。 通常のRATTLE法をこれらの拘束条件に直接適用することはできないが（拘束力 $F_{i}^{\mathrm {con} }$ による 変位 $\delta r_{i}={\frac {(\Delta t)^{2}}{2m_{i}}}F_{i}^{\mathrm {con} }$ が $m_{i}\rightarrow 0$で発散する） 、ダミーサイトに働く力を質量のあるサイトに分配することで通常のRATTLE法と同様に取り扱うことが可能となる。

ダミーサイトの座標は、質量を持つ粒子との拘束条件によって決められる。 従って、ダミーサイトには通常3つの（独立な）拘束条件がかかるべきであり、これ以上の拘束は冗長である。 現行のFreeFlexでは4つ以上の拘束条件を許容しないので拘束条件は冗長に置かないように注意する。