Pucker座標

# Cremer Pople のPucker 座標変換
参考資料
- [D.Cremer and J.A.Pople, "A General Definition of Ring Puckering Coordinates",  
(1975), *JACS*, **97**, 1354-1358](https://pubs.acs.org/doi/10.1021/ja00839a011)  
- $AMBERHOME/AmberTools/src/cpptraj/src/TorsionRoutines.cpp

## 反応座標における Value
六員環を構成する原子の座標を位置ベクトル $ \vec{R}_j (j=1,2,\ldots,5,6) $ とする。
座標系の原点は以下のように取る。（原点は六員環の中心）

```math
\sum_{j=1}^{6}\vec{R}_j = \vec{0}
```
CPPTRAJの計算では6原子の座標から中心を計算し、  
各原子の座標から中心座標を引いた値を実際に計算に用いる座標としている。  
ただしこの処理を行わず、元の座標のまま $\vec{R'} 、 \vec{R''}$ を計算しても  
同じPucker座標に変換可能。  

次に以下のようなベクトルを定義する。

```math
\begin{align}
\vec{R'} &= \sum_{j=1}^{6}\vec{R}_j\sin\bigg[\frac{\pi(j-1)}{3}\bigg]\\
\vec{R''} &= \sum_{j=1}^{6}\vec{R}_j\cos\bigg[\frac{\pi(j-1)}{3}\bigg]\\
\vec{n} &= \frac{\vec{R'} \times \vec{R''}}{|\vec{R'} \times \vec{R''}|}
\end{align}
```

ここで $\times$ はクロス積を表し、
$\vec{A} = (A_x,A_y,A_z) , \vec{B} = (B_x,B_y,B_z)$ とすると
$\vec{A}\times\vec{B} = ((A_yB_z-A_zB_y),(A_zB_x-A_xB_z),(A_xB_y-A_yB_x))$
となる  
この操作により、2つのベクトル $\vec{R'},\vec{R''}$ が成す平面が
六員環の平均面(Average plane, Mean plane)となり、  
ベクトル $\vec{n}$ は平均面の法線ベクトルとなる。
さらに
$$z_j=\vec{R_j}\cdot\vec{n}$$
を計算することで平均的な面からの距離$z_j$が求まる。  

```math
\begin{align}
q_2\sin(\phi) &= \sqrt\frac{1}{3}\sum_{j=1}^{6}z_j\sin\bigg[\frac{2\pi}{3}(j-1)\bigg] &&\quad (0 \leq \phi \leq 2\pi)\\
q_2\cos(\phi) &= - \sqrt\frac{1}{3}\sum_{j=1}^{6}z_j\cos\bigg[\frac{2\pi}{3}(j-1)\bigg]&&\quad (0 \leq \phi \leq 2\pi) \\
q_3 &= \sqrt\frac{1}{6}\sum_{j=1}^{6}(-1)^{j-1}z_j \\
{q_2}^2 &= (q_2\sin(\phi))^2 + (q_2\cos(\phi))^2 \\
Q^2 &= {q_2}^2 + {q_3}^2\\
q_2 &= Q\cos(\theta)&&\quad (0 \leq \theta \leq \pi)\\
q_3 &= Q\sin(\theta)&&\quad (0 \leq \theta \leq \pi)\\
\end{align}
```

上記の数式によりPucker座標（$Q,\theta,\phi$）が定義されている  
またCPPTRAJでは実際には以下の計算により $\phi,\theta$ を計算している。

```math
\theta = \left\{
\begin{array}{ll}
\arctan2(q_2,q_3)  & (\theta \geq 0)\\
\arctan2(q_2,q_3) + \pi  & (\theta \lt 0)
\end{array}
\right.\\
\phi = \left\{
\begin{array}{ll}
\pi - \arcsin\bigg(\frac{q_2\sin(\phi)}{q_2}\bigg)  & (q_2\cos(\phi) \lt 0)\\
\arcsin\bigg(\frac{q_2\sin(\phi)}{q_2}\bigg) + 2\pi  & (q_2\cos(\phi) \geq 0 , \phi \lt 0)\\
\\arcsin\bigg(\frac{q_2\sin(\phi)}{q_2}\bigg) & (q_2\cos(\phi) \geq 0 , \phi \geq 0)
\end{array}
\right.\\
```

C++ の関数arctan2(y,x) の返り値の範囲は $[-\pi,\pi]$ であるが、  
$q_2 \geq 0$ のため、実際の $\theta$ の範囲は $[-\frac{\pi}{2},\frac{\pi}{2}]$ である。

## 反応座標のValueに対する各原子の座標微分からForceを計算する。
まず $\theta$ についての座標微分を計算する。  
ここでは $R_1$ について計算していく。

```math
\theta = \arctan\bigg(\frac{q_2}{q_3}\bigg)\\
```

```math
\begin{align}
\theta を \vec{R_1}_x について偏微分すると\\
\frac{ \partial \theta }{\partial \vec{R_1}} &= \frac{ \partial \theta }{ \partial (\frac{q_2}{q_3})} \frac{ \partial (\frac{q_2}{q_3}) }{ \partial \vec{R_1}} \\
\frac{ \partial \theta }{ \partial (\frac{q_2}{q_3})} &= \frac{1}{1+(\frac{q_2}{q_3})^2}\\
&= \frac{{q_3}^2}{{q_3}^2+{q_2}^2}\\
\frac{ \partial (\frac{q_2}{q_3}) }{ \partial \vec{R_1}} &= \frac{\frac{\partial q_2}{\partial \vec{R_1}}q_3 - \frac{\partial q_3}{\partial \vec{R_1}}q_2} {{q_3}^2}\\
\end{align}
```

さらに $q_2,q_3$ についてそれぞれ $\vec{R_1}_x$ の偏微分を求める。

```math
\begin{align}
q_2 &= \sqrt{{q_2\sin(\phi)}^2 + {q_2\cos(\phi)}^2}\quad より\\
Q_2 &= {q_2\sin(\phi)}^2 + {q_2\cos(\phi)}^2 \quadとおくと\\
\frac{\partial q_2}{\partial \vec{R_1}} &=  \frac{\partial q_2}{\partial Q_2}  \frac{\partial Q_2}{\partial \vec{R_1}}\\
q_2 &= \sqrt{Q_2}  \quadより\\
\frac{\partial q_2}{\partial Q_2} &= \frac{1}{2\sqrt{Q_2}}\\
\frac{\partial Q_2}{\partial \vec{R_1}} &= 2q_2\sin{\phi}\frac{\partial q_2\sin(\phi)}{\partial \vec{R_1}}+2q_2\cos{\phi}\frac{\partial q_2\cos(\phi)}{\partial \vec{R_1}}\\
\end{align}
```


さらに $q_2\sin(\phi) , q_2\cos(\phi)$ について $\vec{R_1}x$ の偏微分を求める。  
ここで $q_2\sin(\phi),q_2\cos(\phi)$ ともに定義式を代入し偏微分を行うことに注意。  
また、 $\vec{R'},\vec{R''}$ の $\vec{R_j}$ の係数にあたる、  
$\sin\bigg[\frac{\pi(j-1)}{3}\bigg],\cos\bigg[\frac{\pi(j-1)}{3}\bigg]$ をそれぞれ、 $k_{s1\_j},k_{c1\_j}$ と、  
また、 $q_2\sin(\phi),q_2\cos(\phi)$ の $z_j$ の係数にあたる、  
$\sin\bigg[\frac{2\pi(j-1)}{3}\bigg],\cos\bigg[\frac{2\pi(j-1)}{3}\bigg]$ をそれぞれ、 $k_{s2\_j},k_{c2\_j}$ とおく。    

```math
\begin{align}
q_2\sin(\phi) &= \sqrt\frac{1}{3}\sum_{j=1}^{6}z_j\sin\bigg[\frac{2\pi}{3}(j-1)\bigg]\\
&= \sqrt\frac{1}{3}\bigg[\sin(0)z_1 + \sin(\frac{2\pi}{3})z_2 + \ldots \sin(\frac{4\pi}{3})z_6\bigg]\\
&定数項を k_{s2\_j} で置換すると\\
q_2\sin(\phi) &= k_{s2\_1} * z_1 + k_{s2\_2} * z_2 + \ldots + k_{s2\_6} * z_6 \\
\frac{\partial q_2\sin(\phi)}{\partial \vec{R_1}} &= k_{s2\_1} * \frac{\partial z_1}{\partial \vec{R_1}} + k_{s2\_2} * \frac{\partial z_2}{\partial \vec{R_1}} + \ldots + k_{s2\_6} * \frac{\partial z_6}{\partial \vec{R_1}} \\
\end{align}
```

$z_1,z_2,\ldots$ それぞれについて $\vec{R_1}_x$ での偏微分を求める。  
この後の方針としては、$z_j$ についての $\vec{R_1}_x$ での偏微分を、  
$\vec{R_1}_x=(R_1x,R_1y,R_1z)$ での各成分での偏微分を計算する。  
その過程で、単位法線ベクトル $\vec{n}$ の各成分に対する偏微分が必要となり、  
これを求めるのに、 $\vec{R'}\times\vec{R''}$ の各成分、および絶対値の偏微分を計算する。  

```math
\begin{align}
z_1 &= \vec{R_1} \cdot \vec{n} \\
&= \vec{R_1}_x * \vec n_x + \vec{R_1}_y * \vec n_y + \vec{R_1}_z * \vec n_z\\
\frac{\partial z_1}{\partial \vec{R_1}} &= \bigg( \frac{\partial z_1}{\partial \vec{R_1}_x} , \frac{\partial z_1}{\partial \vec{R_1}_y} ,\frac{\partial z_1}{\partial \vec{R_1}_z} \bigg)\\
&ここでまず\vec{R_1} のx成分の微分を求める。\\
\frac{\partial z_1}{\partial \vec{R_1}_x} &= \frac{\partial(\vec{R_1}_x * \vec n_x + \vec{R_1}_y * \vec n_y + \vec{R_1}_z * \vec n_z)}{\partial \vec{R_1}_x}\\
&= \vec n_x + \vec{R_1}_x \frac{\partial \vec{n}_x}{\partial \vec{R_1}_x } +  \vec{R_1}_x \frac{\partial \vec{n}_y}{\partial \vec{R_1}_x } + \vec{R_1}_x \frac{\partial \vec{n}_z}{\partial \vec{R_1}_x }\\
\vec{R'R''} &= \vec{R'} \times \vec{R''} = (RR_x , RR_y , RR_z) とおくと\\
|\vec{R'R''}| &= \sqrt{(RR_x)^2 + (RR_y)^2 + (RR_z)^2}\\
\vec{n} &= \frac{\vec{R'R''}}{|\vec{R'R''}|} より、\\

\frac{\partial \vec{n}_x}{\partial \vec{R_1}_x } &= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\frac {RR_x}{|\vec{R'R''}|}\bigg)\\
&= \frac{1}{|\vec{R'R''}|^2}\bigg( |\vec{R'R''}| \frac{\partial RR_x}{\partial \vec{R_1}_x} - RR_x \frac{\partial |\vec{R'R''}|}{\partial \vec{R_1}_x}  \bigg)\\
\\
\frac{\partial RR_x}{\partial \vec{R_1}_x} &= \frac{\partial (\vec{R'}_y\vec{R''}_z - \vec{R'}_z\vec{R''}_y)}{\partial \vec{R_1}_x} \\
&= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\sum_{j=1}^6\vec{R_j}_y\sin\bigg[\frac{\pi(j-1)}{3}\bigg] * \sum_{j=1}^6\vec{R_j}_z\cos\bigg[\frac{\pi(j-1)}{3}\bigg] - \sum_{j=1}^6\vec{R_j}_z\sin\bigg[\frac{\pi(j-1)}{3}\bigg] * \sum_{j=1}^6\vec{R_j}_y\cos\bigg[\frac{\pi(j-1)}{3}\bigg]\bigg)\\
&= 0\\
\end{align}
```

※  $\vec{R_1}_x$ 成分を含まないため。  

```math
\begin{align}
\\
\frac{\partial RR_y}{\partial \vec{R_1}_x} &= \frac{\partial (\vec{R'}_z\vec{R''}_x - \vec{R'}_x\vec{R''}_z)}{\partial \vec{R_1}_x}\\
&= \frac{\partial \vec{R'}_z\vec{R''}_x}{\partial \vec{R_1}_x} - \frac{\partial \vec{R'}_x\vec{R''}_z}{\partial \vec{R_1}_x}\\
\frac{\partial \vec{R'}_z\vec{R''}_x}{\partial \vec{R_1}_x} &= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\vec{R'}_z * \sum_{j=1}^6\vec{R_j}_x\cos\bigg[\frac{\pi(j-1)}{3}\bigg]\bigg)\\
&=k_{c1\_1} * \vec{R'}_z\\
\frac{\partial \vec{R'}_x\vec{R''}_z}{\partial \vec{R_1}_x} &= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\sum_{j=1}^6\vec{R_j}_x\sin\bigg[\frac{\pi(j-1)}{3}\bigg] * \vec{R''}_z \bigg)\\
&= k_{s1\_1} * \vec{R''}_z \\
\frac{\partial RR_y}{\partial \vec{R_1}_x} &= k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z \\
同様に\\
\frac{\partial RR_z}{\partial \vec{R_1}_x} &= \frac{\partial (\vec{R'}_x \vec{R''}_y - \vec{R'}_y\vec{R''}_x)}{\partial \vec{R_1}_x}\\
&= k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y\\
\\
\frac{\partial |\vec{R'R''}|}{\partial \vec{R_1}_x} &= \frac{\partial \sqrt{(RR_x)^2 + (RR_y)^2 + (RR_z)^2}}{\partial \vec{R_1}_x}\\
&= \frac {\partial\sqrt{(RR_x)^2 + (RR_y)^2 + (RR_z)^2}}{\partial \bigg( (RR_x)^2 + (RR_y)^2 + (RR_z)^2\bigg)} \frac{\partial \bigg( (RR_x)^2 + (RR_y)^2 + (RR_z)^2\bigg)}{\partial \vec{R_1}_x}\\
&= \frac{1}{2\sqrt{(RR_x)^2 + (RR_y)^2 + (RR_z)^2}}\bigg(2RR_x\frac{\partial RR_x}{\partial \vec{R_1}_x} + 2RR_y\frac{\partial RR_y}{\partial\vec{R_1}_x} + 2RR_z\frac{\partial RR_z}{\partial \vec{R_1}_x} \bigg)\\
&=\frac{1}{2|\vec{R'R''}|}\bigg( 2RR_x * 0 + 2RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + 2RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )   \bigg) \\
&= \frac{1}{|\vec{R'R''}|}\bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  )  \bigg)\\
\\
\frac{\partial \vec{n}_x}{\partial \vec{R_1}_x } &= \frac{1}{|\vec{R'R''}|^2}\bigg( |\vec{R'R''}| \frac{\partial RR_x}{\partial \vec{R_1}_x} - RR_x \frac{\partial |\vec{R'R''}|}{\partial \vec{R_1}_x}  \bigg)\\
&= \frac{-RR_x}{|\vec{R'R''}|^2} \frac{1}{|\vec{R'R''}|}\bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  )  \bigg)\\
&= \frac{- RR_x}{|\vec{R'R''}|^3} \bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  )  \bigg)\\

\\
\\
\frac{\partial \vec{n}_y}{\partial \vec{R_1}_x } &= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\frac{RR_y}{|\vec{R'R''}|}  \bigg)\\
&= \frac{1}{|\vec{R'R''}|^2} \bigg( |\vec{R'R''}| \frac{ \partial RR_y}{\partial \vec{R_1}_x} - RR_y \frac{\partial |\vec{R'R''}|}{\partial \vec{R_1}_x} \bigg)\\
&= \frac{1}{|\vec{R'R''}|^2} \bigg[ |\vec{R'R''}| * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) \\
&- RR_y * \frac{1}{|\vec{R'R''}|}\bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )    \bigg) \bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ |\vec{R'R''}|^2 * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) \\
&- RR_y * \bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  )  \bigg) \bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ (|\vec{R'R''}|^2 - {RR_y}^2 ) * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) - RR_y * RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  \bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ ({RR_x}^2 + {RR_z}^2 ) * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) - RR_y *  RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y  ) \bigg]\\
\\
\\
\frac{\partial \vec{n}_z}{\partial \vec{R_1}_x } &= \frac{\partial}{\partial \vec{R_1}_x}\bigg(\frac{RR_z}{|\vec{R'R''}|}  \bigg)\\
&= \frac{1}{|\vec{R'R''}|^2} \bigg( |\vec{R'R''}| \frac{ \partial RR_z}{\partial \vec{R_1}_x} - RR_z \frac{\partial |\vec{R'R''}|}{\partial \vec{R_1}_x} \bigg)\\
&= \frac{1}{|\vec{R'R''}|^2} \bigg[ |\vec{R'R''}|* (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  \\
&- RR_z * \frac{1}{|\vec{R'R''}|}\bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  \bigg)\bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ |\vec{R'R''}|^2 * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  \\
&- RR_z * \bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  \bigg)\bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ (|\vec{R'R''}|^2 - {RR_z}^2 )* (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  - RR_z * RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)    \bigg]\\
&= \frac{1}{|\vec{R'R''}|^3} \bigg[ ( {RR_x}^2+{RR_y}^2 )* (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  - RR_z * RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)    \bigg]\\
\\
\\
\end{align}
```


**ここまでチェック終了**  
これまでの結果をここで一度まとめると以下のようになる。  

```math
\begin{align}

\frac{\partial \vec{n}_x}{\partial \vec{R_1}_x } &=\frac{- RR_x}{|\vec{R'R''}|^3} \bigg( RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)  + RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  )  \bigg)\\
\frac{\partial \vec{n}_y}{\partial \vec{R_1}_x } &=\frac{1}{|\vec{R'R''}|^3} \bigg[ ({RR_x}^2 + {RR_z}^2 ) * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) - RR_y *  RR_z * (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y  ) \bigg]\\
\frac{\partial \vec{n}_z}{\partial \vec{R_1}_x } &=\frac{1}{|\vec{R'R''}|^3} \bigg[ ( {RR_x}^2+{RR_y}^2 )* (k_{s1\_1} * \vec{R''}_y - k_{c1\_1} * \vec{R'}_y )  - RR_z * RR_y * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z)    \bigg]\\
\frac{\partial z_1}{\partial \vec{R_1}_x} &= \vec n_x + \vec{R_1}_x \frac{\partial \vec{n}_x}{\partial \vec{R_1}_x } +  \vec{R_1}_x \frac{\partial \vec{n}_y}{\partial \vec{R_1}_x } + \vec{R_1}_x \frac{\partial \vec{n}_z}{\partial \vec{R_1}_x }\\
&= \frac{RR_x}{|\vec{R'R''}|} + \frac{\vec{R_1}_x}{|\vec{R'R''}|^3}\bigg[ \bigg( (-RR_xRR_y+{RR_x}^2+{RR_z}^2-RR_zRR_y) * (k_{c1\_1} * \vec{R'}_z - k_{s1\_1} * \vec{R''}_z) \bigg) + \bigg( \bigg)  \bigg]\\
\end{align}
```


## Pucker 座標計算の実例
以下の例について求める。

| Atom |  j  |      x |      y |     z |
|:----:|:---:| ------:| ------:| -----:|
|  O5  |  1  |  0.170 |  1.518 | 1.980 |
|  C1  |  2  | -1.345 |  1.518 | 1.980 |
|  C2  |  3  | -1.898 |  0.108 | 1.979 |
|  C3  |  4  | -1.348 | -0.698 | 0.820 |
|  C4  |  5  | -0.367 |  0.116 |     0 |
|  C5  |  6  |  0.719 |  0.713 | 0.820 |

原点は (-0.678,0.546,1.263)
原点移動後の座標は

| Atom |  j  |      x |      y |      z |
|:----:|:---:| ------:| ------:| ------:|
|  O5  |  1  |  0.848 |  0.972 |  0.717 |
|  C1  |  2  | -0.669 |  0.972 |  0.717 |
|  C2  |  3  | -1.220 | -0.438 |  0.716 |
|  C3  |  4  | -0.669 | -1.244 | -0.443 |
|  C4  |  5  |  0.311 | -0.430 | -1.263 |
|  C5  |  6  |  1.397 |  0.167 | -0.443 |

$\sin\bigg[\frac{\pi}{3}(j-1)\bigg] , \cos\bigg[\frac{\pi}{3}(j-1)\bigg],\sin\bigg[\frac{2\pi}{3}(j-1)\bigg] , \cos\bigg[\frac{2\pi}{3}(j-1)\bigg]$ の  
値をそれぞれ計算して表に書き加えると、  

| Atom |  j  |      x |      y |      z | sin1_j | cos1_j | sin2_j | cos2_j |
|:----:|:---:| ------:| ------:| ------:| ------:| ------:| ------:| ------:|
|  O5  |  1  |  0.848 |  0.972 |  0.717 |    0.0 |    1.0 |    0.0 |    1.0 |
|  C1  |  2  | -0.669 |  0.972 |  0.717 |  0.866 |    0.5 |  0.866 |   -0.5 |
|  C2  |  3  | -1.220 | -0.438 |  0.716 |  0.866 |   -0.5 | -0.866 |   -0.5 |
|  C3  |  4  | -0.670 | -1.244 | -0.443 |    0.0 |   -1.0 |    0.0 |    1.0 |
|  C4  |  5  |  0.311 | -0.430 | -1.263 | -0.866 |   -0.5 |  0.866 |   -0.5 |
|  C5  |  6  |  1.397 |  0.167 | -0.443 | -0.866 |    0.5 | -0.866 |   -0.5 |

上の表の値を $\vec{R'},\vec{R''}$ を求める式に代入してそれぞれ計算すると、
  
```math
\begin{align}
\vec{R'} &= (-3.113,0.690,2.718)\\
\vec{R''} &= (2.338,3.220,1.571)
\end{align}
```

さらに法線ベクトル $\vec n$ を求めると

```math
\vec n = ( -0.428 , 0.628 ,-0.650 )
```

法線ベクトル $\vec n$ と $\vec{R_j}$ の内積により、平均面からの距離 $z_j$ を計算し表に挿入すると、
  
| Atom | j   | X      | Y      | Z      | sin1_j | cos1_j | sin2_j | cos2_j | $z_j$  |
| ---- | --- | ------ | ------ | ------ | ------ | ------ | ------ | ------ | ------ |
| O5   | 1   | 0.848  | 0.972  | 0.717  | 0.000  | 1.000  | 0.000  | 1.000  | -0.219 |
| C2   | 2   | -1.220 | -0.438 | 0.716  | 0.866  | -0.500 | -0.866 | -0.500 | -0.218 |
| C1   | 3   | -0.667 | 0.972  | 0.717  | 0.866  | 0.500  | 0.866  | -0.500 | 0.430  |
| C3   | 4   | -0.670 | -1.244 | -0.443 | 0.000  | -1.000 | -0.000 | 1.000  | -0.206 |
| C4   | 5   | 0.311  | -0.430 | -1.263 | -0.866 | -0.500 | 0.866  | -0.500 | 0.418  |
| C5   | 6   | 1.397  | 0.167  | -0.443 | -0.866 | 0.500  | -0.866 | -0.500 | -0.205 |

さらに $q_2\cos(\phi) , q_2\sin(\phi) , q_3,$ および Pucker 座標を計算すると  

```math
\begin{align}
q_2\cos(\phi) &= -0.368\\
q_2\sin(\phi) &= -0.635\\
q_3 &= -0.015\\
Q &= 0.734 \\
\theta &= 91.19\\
\phi &=  239.93\\
\end{align}
```

注意点として、Pucker を計算する際に六員環を構成する原子を指定しており、  
指定の際は {C1,C2,C3,C4,C5,O5} の順で入力しているが、内部では上記の表のとおり  
ナンバリングが修正されている。  
新たに反応座標として Pucker を定義する際に CPPTRAJ での定義と齟齬がないようにする必要がある。  


```TorsionRoutines.cpp
# CPPTRAJ による計算
// Pucker_CP()
/** Return the pucker (in radians) of coords in a1-a5 based on method of
  * Cremer & Pople.
  */
double Pucker_CP(const double* a1, const double* a2, const double* a3,
                 const double* a4, const double* a5, const double* a6,
                 int N, double& amplitude, double& theta)
{
  double dN = (double)N;
  double twopi_over_N = Constants::TWOPI / dN;
  Vec3 XYZ[6];
  XYZ[1].Assign( a1 );
  XYZ[2].Assign( a2 );
  XYZ[3].Assign( a3 );
  XYZ[4].Assign( a4 );
  if (N == 5) {
    XYZ[0].Assign( a5 ); // Ring apex
  } else if (N == 6) {
    XYZ[0].Assign( a6 ); // Ring apex
    XYZ[5].Assign( a5 );
  } else
    return -1.0; // Internal error
  // Calculate geometric center
  Vec3 Center(0.0, 0.0, 0.0);
  for (int i = 0; i < N; i++)
    Center += XYZ[i];
  Center /= dN;
  // Translate to center
  for (int i = 0; i < N; i++)
    XYZ[i] -= Center;
  // Calculate position vectors
  Vec3 R1(0.0, 0.0, 0.0);
  Vec3 R2(0.0, 0.0, 0.0);
  for (int i = 0; i < N; i++) {
    double factor = twopi_over_N * (double)i;
    double sin_val = sin( factor );
    double cos_val = cos( factor );
    R1[0] += XYZ[i][0] * sin_val;
    R2[0] += XYZ[i][0] * cos_val;
    R1[1] += XYZ[i][1] * sin_val;
    R2[1] += XYZ[i][1] * cos_val;
    R1[2] += XYZ[i][2] * sin_val;
    R2[2] += XYZ[i][2] * cos_val;
  }
  // Calculate vector normal to plane
  Vec3 NXYZ = R1.Cross( R2 );
  // Normalize
  NXYZ.Normalize();
  // Calculate sums for determining q2/phi2
  // NOTE: Total amplitude Q should = sqrt(SUM[ Zn*Zn ])
  double sum1 = 0.0;
  double sum2 = 0.0;
  double Zn[6];
  for (int i = 0; i < N; i++) {
    double factor = 2.0 * twopi_over_N * (double)i;
    double cos_val = cos( factor );
    double sin_val = sin( factor );
    Zn[i] = XYZ[i] * NXYZ;
    sum1 += Zn[i] * cos_val;
    sum2 -= Zn[i] * sin_val;
  }
  // Calculate amplitude (q2, ==Q for N==5)
  double norm = sqrt(sum1*sum1 + sum2*sum2);
  amplitude = norm * sqrt( 2.0 / dN );
  // For even # coords (only 6 currently) calc extra pucker coord (q3)
  if (N == 6) {
    double q3 = 0.0;
    double mult = 1.0;
    for (int i = 0; i < N; i++) {
      q3 += mult * Zn[i]; // mult ~ pow( -1.0, i )
      mult = -mult;
    }
    q3 /= sqrt( dN );
    // Calculate theta
    theta = atan2( amplitude, q3 );
    // Update amplitude (Q for N==6)
    amplitude = sqrt( amplitude*amplitude + q3*q3 );
  }
  // Calculate pucker (phi2)
  double pucker = asin( sum2 / norm );
  if (sum1 < 0.0)
    pucker = Constants::PI - pucker;
  else if (pucker < 0.0)
    pucker += Constants::TWOPI;

  return pucker;
}
```