--- trunk/tengDissertation/Introduction.tex	2006/04/18 04:11:56	2718
+++ trunk/tengDissertation/Introduction.tex	2006/04/18 22:38:19	2719
@@ -822,7 +822,7 @@ q(\Delta t)} \right]. %
 %
 q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
 q(\Delta t)} \right]. %
- \label{introEquation:positionVerlet1}
+ \label{introEquation:positionVerlet2}
 \end{align}
 
 \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
@@ -887,7 +887,21 @@ dynamical information.
 has proven to be a powerful tool for studying the functions of
 biological systems, providing structural, thermodynamic and
 dynamical information.
+
+One of the principal tools for modeling proteins, nucleic acids and
+their complexes. Stability of proteins Folding of proteins.
+Molecular recognition by:proteins, DNA, RNA, lipids, hormones STP,
+etc. Enzyme reactions Rational design of biologically active
+molecules (drug design) Small and large-scale conformational
+changes. determination and construction of 3D structures (homology,
+Xray diffraction, NMR) Dynamic processes such as ion transport in
+biological systems.
 
+Macroscopic properties are related to microscopic behavior.
+
+Time dependent (and independent) microscopic behavior of a molecule
+can be calculated by molecular dynamics simulations.
+
 \subsection{\label{introSec:mdInit}Initialization}
 
 \subsection{\label{introSec:forceEvaluation}Force Evaluation}
@@ -927,10 +941,10 @@ rotation matrix $A$ and re-formulating Hamiltonian's e
 The break through in geometric literature suggests that, in order to
 develop a long-term integration scheme, one should preserve the
 symplectic structure of the flow. Introducing conjugate momentum to
-rotation matrix $A$ and re-formulating Hamiltonian's equation, a
+rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
 symplectic integrator, RSHAKE, was proposed to evolve the
 Hamiltonian system in a constraint manifold by iteratively
-satisfying the orthogonality constraint $A_t A = 1$. An alternative
+satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
 method using quaternion representation was developed by Omelyan.
 However, both of these methods are iterative and inefficient. In
 this section, we will present a symplectic Lie-Poisson integrator
@@ -1136,8 +1150,8 @@ To reduce the cost of computing expensive functions in
    0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
 \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
 \]
-To reduce the cost of computing expensive functions in e^{\Delta
-tR_1 }, we can use Cayley transformation,
+To reduce the cost of computing expensive functions in $e^{\Delta
+tR_1 }$, we can use Cayley transformation,
 \[
 e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
 )
@@ -1213,15 +1227,15 @@ Moreover, \varphi _{\Delta t/2,V} can be divided into 
 \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
 _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
 \]
-Moreover, \varphi _{\Delta t/2,V} can be divided into two sub-flows
-which corresponding to force and torque respectively,
+Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
+sub-flows which corresponding to force and torque respectively,
 \[
 \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
 _{\Delta t/2,\tau }.
 \]
 Since the associated operators of $\varphi _{\Delta t/2,F} $ and
 $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
-order inside \varphi _{\Delta t/2,V} does not matter.
+order inside $\varphi _{\Delta t/2,V}$ does not matter.
 
 Furthermore, kinetic potential can be separated to translational
 kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
@@ -1251,129 +1265,124 @@ generalized Langevin Dynamics will be given first. Fol
 mimics a simple heat bath with stochastic and dissipative forces,
 has been applied in a variety of studies. This section will review
 the theory of Langevin dynamics simulation. A brief derivation of
-generalized Langevin Dynamics will be given first. Follow that, we
+generalized Langevin equation will be given first. Follow that, we
 will discuss the physical meaning of the terms appearing in the
 equation as well as the calculation of friction tensor from
 hydrodynamics theory.
 
-\subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
+\subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
 
+Harmonic bath model, in which an effective set of harmonic
+oscillators are used to mimic the effect of a linearly responding
+environment, has been widely used in quantum chemistry and
+statistical mechanics. One of the successful applications of
+Harmonic bath model is the derivation of Deriving Generalized
+Langevin Dynamics. Lets consider a system, in which the degree of
+freedom $x$ is assumed to couple to the bath linearly, giving a
+Hamiltonian of the form
 \begin{equation}
 H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
-\label{introEquation:bathGLE}
+\label{introEquation:bathGLE}.
 \end{equation}
-where $H_B$ is harmonic bath Hamiltonian,
+Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
+with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
 \[
-H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
-}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
+H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
+}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
+\right\}}
 \]
-and $\Delta U$ is bilinear system-bath coupling,
+where the index $\alpha$ runs over all the bath degrees of freedom,
+$\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
+the harmonic bath masses, and $\Delta U$ is bilinear system-bath
+coupling,
 \[
 \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
 \]
-Completing the square,
+where $g_\alpha$ are the coupling constants between the bath and the
+coordinate $x$. Introducing
 \[
-H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
-{\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
-w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
-w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
-1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
-\]
-and putting it back into Eq.~\ref{introEquation:bathGLE},
+W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
+}}{{2m_\alpha  w_\alpha ^2 }}} x^2
+\] and combining the last two terms in Equation
+\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
+Hamiltonian as
 \[
 H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
 {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
 w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
 w_\alpha ^2 }}x} \right)^2 } \right\}}
 \]
-where
-\[
-W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
-}}{{2m_\alpha  w_\alpha ^2 }}} x^2
-\]
 Since the first two terms of the new Hamiltonian depend only on the
 system coordinates, we can get the equations of motion for
 Generalized Langevin Dynamics by Hamilton's equations
 \ref{introEquation:motionHamiltonianCoordinate,
 introEquation:motionHamiltonianMomentum},
-\begin{align}
-\dot p &=  - \frac{{\partial H}}{{\partial x}}
-       &= m\ddot x
-       &= - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)}
-\label{introEquation:Lp5}
-\end{align}
-, and
-\begin{align}
-\dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
-                &= m\ddot x_\alpha
-                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
-\end{align}
+\begin{equation}
+m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
+\sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
+\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
+\label{introEquation:coorMotionGLE}
+\end{equation}
+and
+\begin{equation}
+m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
+\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
+\label{introEquation:bathMotionGLE}
+\end{equation}
 
-\subsection{\label{introSection:laplaceTransform}The Laplace Transform}
+In order to derive an equation for $x$, the dynamics of the bath
+variables $x_\alpha$ must be solved exactly first. As an integral
+transform which is particularly useful in solving linear ordinary
+differential equations, Laplace transform is the appropriate tool to
+solve this problem. The basic idea is to transform the difficult
+differential equations into simple algebra problems which can be
+solved easily. Then applying inverse Laplace transform, also known
+as the Bromwich integral, we can retrieve the solutions of the
+original problems.
 
+Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
+transform of f(t) is a new function defined as
 \[
-L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
+L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
 \]
+where  $p$ is real and  $L$ is called the Laplace Transform
+Operator. Below are some important properties of Laplace transform
+\begin{equation}
+\begin{array}{c}
+ L(x + y) = L(x) + L(y) \\
+ L(ax) = aL(x) \\
+ L(\dot x) = pL(x) - px(0) \\
+ L(\ddot x) = p^2 L(x) - px(0) - \dot x(0) \\
+ L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p) \\
+ \end{array}
+\end{equation}
 
+Applying Laplace transform to the bath coordinates, we obtain
 \[
-L(x + y) = L(x) + L(y)
+\begin{array}{c}
+ p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) \\
+ L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
+ \end{array}
 \]
-
+By the same way, the system coordinates become
 \[
-L(ax) = aL(x)
+\begin{array}{c}
+ mL(\ddot x) =  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
+  - \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x) - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0) - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}  \\
+ \end{array}
 \]
 
+With the help of some relatively important inverse Laplace
+transformations:
 \[
-L(\dot x) = pL(x) - px(0)
-\]
-
-\[
-L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
-\]
-
-\[
-L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
-\]
-
-Some relatively important transformation,
-\[
-L(\cos at) = \frac{p}{{p^2  + a^2 }}
-\]
-
-\[
-L(\sin at) = \frac{a}{{p^2  + a^2 }}
-\]
-
-\[
-L(1) = \frac{1}{p}
+\begin{array}{c}
+ L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
+ L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
+ L(1) = \frac{1}{p} \\
+ \end{array}
 \]
-
-First, the bath coordinates,
-\[
-p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
-_\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
-}}L(x)
-\]
-\[
-L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
-px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
-\]
-Then, the system coordinates,
+, we obtain
 \begin{align}
-mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
-\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
-}}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
-(0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
-}}\omega _\alpha ^2 L(x)} \right\}}
-%
- &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
- \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
- - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
- - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
-\end{align}
-Then, the inverse transform,
-
-\begin{align}
 m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
 }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
@@ -1392,61 +1401,116 @@ t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  =
 (\omega _\alpha  t)} \right\}}
 \end{align}
 
+Introducing a \emph{dynamic friction kernel}
 \begin{equation}
+\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
+}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
+\label{introEquation:dynamicFrictionKernelDefinition}
+\end{equation}
+and \emph{a random force}
+\begin{equation}
+R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
+- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
+\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
+\label{introEquation:randomForceDefinition}
+\end{equation}
+the equation of motion can be rewritten as
+\begin{equation}
 m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
 (t)\dot x(t - \tau )d\tau }  + R(t)
 \label{introEuqation:GeneralizedLangevinDynamics}
 \end{equation}
-%where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
-%$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
+which is known as the \emph{generalized Langevin equation}.
+
+\subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
+
+One may notice that $R(t)$ depends only on initial conditions, which
+implies it is completely deterministic within the context of a
+harmonic bath. However, it is easy to verify that $R(t)$ is totally
+uncorrelated to $x$ and $\dot x$,
 \[
-\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
-}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
+\begin{array}{l}
+ \left\langle {x(t)R(t)} \right\rangle  = 0, \\
+ \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
+ \end{array}
 \]
-For an infinite harmonic bath, we can use the spectral density and
-an integral over frequencies.
+This property is what we expect from a truly random process. As long
+as the model, which is gaussian distribution in general, chosen for
+$R(t)$ is a truly random process, the stochastic nature of the GLE
+still remains.
 
+%dynamic friction kernel
+The convolution integral
 \[
-R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
-- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
-\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
-(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }
 \]
-The random forces depend only on initial conditions.
+depends on the entire history of the evolution of $x$, which implies
+that the bath retains memory of previous motions. In other words,
+the bath requires a finite time to respond to change in the motion
+of the system. For a sluggish bath which responds slowly to changes
+in the system coordinate, we may regard $\xi(t)$ as a constant
+$\xi(t) = \Xi_0$. Hence, the convolution integral becomes
+\[
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
+\]
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
+\[
+m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
+\frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
+\]
+which can be used to describe dynamic caging effect. The other
+extreme is the bath that responds infinitely quickly to motions in
+the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
+time:
+\[
+\xi (t) = 2\xi _0 \delta (t)
+\]
+Hence, the convolution integral becomes
+\[
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
+{\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
+\]
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
+\begin{equation}
+m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
+x(t) + R(t) \label{introEquation:LangevinEquation}
+\end{equation}
+which is known as the Langevin equation. The static friction
+coefficient $\xi _0$ can either be calculated from spectral density
+or be determined by Stokes' law for regular shaped particles.A
+briefly review on calculating friction tensor for arbitrary shaped
+particles is given in section \ref{introSection:frictionTensor}.
 
 \subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
-So we can define a new set of coordinates,
+
+Defining a new set of coordinates,
 \[
 q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
 ^2 }}x(0)
-\]
-This makes
+\],
+we can rewrite $R(T)$ as
 \[
-R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
+R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
 \]
 And since the $q$ coordinates are harmonic oscillators,
 \[
-\begin{array}{l}
+\begin{array}{c}
+ \left\langle {q_\alpha ^2 } \right\rangle  = \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
  \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
+ \left\langle {R(t)R(0)} \right\rangle  = \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
+  = \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
+  = kT\xi (t) \\
  \end{array}
 \]
-
-\begin{align}
-\left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
-{\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
-(t)q_\beta  (0)} \right\rangle } }
-%
-&= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
-\right\rangle \cos (\omega _\alpha  t)}
-%
-&= kT\xi (t)
-\end{align}
-
+Thus, we recover the \emph{second fluctuation dissipation theorem}
 \begin{equation}
 \xi (t) = \left\langle {R(t)R(0)} \right\rangle
-\label{introEquation:secondFluctuationDissipation}
+\label{introEquation:secondFluctuationDissipation}.
 \end{equation}
+In effect, it acts as a constraint on the possible ways in which one
+can model the random force and friction kernel.
 
 \subsection{\label{introSection:frictionTensor} Friction Tensor}
 Theoretically, the friction kernel can be determined using velocity
@@ -1622,7 +1686,7 @@ where \delta _{ij} is Kronecker delta function. Invert
 B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
 )T_{ij}
 \]
-where \delta _{ij} is Kronecker delta function. Inverting matrix
+where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
 $B$, we obtain
 
 \[
@@ -1666,7 +1730,7 @@ translation-rotation coupling resistance tensor do dep
 Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
 we can easily find out that the translational resistance tensor is
 origin independent, while the rotational resistance tensor and
-translation-rotation coupling resistance tensor do depend on the
+translation-rotation coupling resistance tensor depend on the
 origin. Given resistance tensor at an arbitrary origin $O$, and a
 vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
 obtain the resistance tensor at $P$ by
