Classical Forcefield derivation - ProkopHapala/FireCore GitHub Wiki

Torsion

Let's consider torsion between atoms A,B,C,D where A,B are end-points and bond CD is the axis. Atom C is bonded to A, and D is bonded to B. Let's define vectors along the bonds ${\vec a}=A-C$, ${\vec b}=B-D$ and ${\vec c}=D-C$. It is also useful to define unitary vector along axis ${\hat c} = {\vec c}/|\vec c|$.

Energy

The energy of torsion $E(\phi)$ should depend on only on angle $\phi$ between vectors $\vec a$ and $\vec b$ projected to plane perpendicular to the axis ${\hat c}$. But the angle ACD and BDC can be arbitrary, therefor it is useful to define also these projected vectors which are perpendicular to $\hat c$. That is ${\vec u}={\vec a}-{\hat c}({\vec a} \cdot {\hat c})$ and ${\vec v}={\vec b}-{\hat c}({\vec b} \cdot {\hat c})$. The length of these vectors $l_u=|\vec u|$ and $l_v=|\vec v|$ and their normalized counterparts $\hat u = {\vec u}/l_u$ and $\hat v={\vec v}/l_v$.

Forces

Now we know that forces on the endpoints can be expressed using torque $\tau_c$ along axis ${\hat c}$ which is a derivative of of energy which respect to angle $\phi$. Formally it useful to express it as a vector along axis $\tau_c$, therefore ${\vec \tau} = {\hat c} (\partial E(\phi) / \partial \phi)$. Now the forces on the endpoints (atoms A,B) can be expressed as ${f}_A = {\vec \tau}_c \times {\hat u}/(l_u)$ and ${f}_B = {\vec \tau}_c \times {\hat v}/l_v$. Or also ${\vec f}_A = {\vec \tau}_c \times {\vec u}/|\vec u|^2$ and ${\vec f}_B = {\vec \tau}_c \times {\vec v}/|\vec v|^2$.

Recoil forces for torsion

Any isolated physical systems must conserve linear and angular momentum. Therefore also the forces acting between particles (e.g. atoms) cannot introduce any impulse on the center of mass, nor torque. This means that the forces and torques introduced by ${\vec f}_A$ and ${\vec f}_B$ acting in end-points A,B must be compensated by ${\vec f}_C$ and ${\vec f}_D$ acting at axial points C,D. Notice that by definition the forces ${\vec f}_A$ and ${\vec f}_B$ has zero component along axis ${\hat c}$. Therefore also ${\vec f}_C$ and ${\vec f}_D$ must be zero along this axis (note: just to be precise they can still have non-zero but exactly opposite component along this axis, but there is no physical reason why there should be such force, also there is bond between CD which fix their distance). Without loss of generality we may consider just 2D system where all the forces act just in plane perpendicular to ${\hat c}$ (i.e. we have only components x,y assuming ${\hat c}$ is along z-axis). Now we have just 4 free degrees of freedom for ${\vec f}_C = (f^x_C,f^y_C)$ and ${\vec f}_D=(f^x_D,f^y_D)$ and already two conditions following from conservation of linear momentum ${\vec f}_C + {\vec f}_D = {\vec f}_G = -{\vec f}_A - {\vec f}_B$ where ${\vec f}_G = (f^x_G,f^y_G)$ is linear recoil on center of mass due to ${\vec f}_C$ and ${\vec f}_D$.

So we have conditions $f^x_C + f^x_C = f^x_G$ and $f^y_C + f^y_C = f^y_G$. We have to envroce additional two conditions to determine fullty the recoil forces ${\vec f}_C$ and ${\vec f}_D$. This can be done easily considering that total torque must be zero along arbitery axis.

Notice, that it may seem that torque changes drastically depending on placement of the axis. But actually shift o the axis of rotation does not change the torque as long as force on center of mass is zero. This is due to linearity of torque formulas which can be shown component-wise.