Quadrupole moment

Just as the dipole moment yields a reasonable origin of coordinates, the quadrupole moment establishes a plausible scale of coordinates. This use is based on Tchebychev's theorem, which states that the fraction of the data lying more than $n$ standard deviations from the mean is $1/n^2$. This is an absolute guarantee, not depending on a Gaussian or any other distribution; consequently it tends to be fairly conservative. In rough terms, less than 11% of the data lies more than 3 standard deviations away from the average.

The theorem is so easy to prove, that it might as well be done right here. By definition, the square of the standard deviation is the second moment of the data with respect to the first momment (suppose it is zero to simplify the calculation):

$$\begin{eqnarray*} \sigma^2 & = & \frac{1}{N}\sum_{i=1}^N x_i^2. \end{eqnarray*}$$

Divide the data into two subsets, $A$ for which $\vert x_i\vert\leq \rho \sigma$ and $B$ for the remainder:

$$\begin{eqnarray*} \sigma^2 & = & \frac{1}{N}\sum_A x_i^2 + \frac{1}{N}\sum_B x_i^2. \end{eqnarray*}$$

Create an inequality by discarding data from set $A$ and underestimating each term in the sum for $B$:

$$\begin{eqnarray*} \sigma^2 & \geq & \frac{1}{N}\sum_B (\rho\sigma)^2. \end{eqnarray*}$$

Suppose there were $M$ items in $B$, the set with the larger data. Then

$$\begin{eqnarray*} \sigma^2 & \geq & \frac{M}{N}(\rho\sigma)^2, \end{eqnarray*}$$

which readily translates into

$$\begin{eqnarray*} {\rm Probability}(x\geq \rho\sigma) & \leq & \frac{1}{\rho^2}. \end{eqnarray*}$$

Applied to multidimensional data, a similar case could be made for quadratic forms constructed from the components of the data, but it is preferable to go ahead and use linear algebra to realize that the eigenvalues of the second moment, or quadrupole, matrix have a similar significance, providing three separate scales for three orthogonal directions. An average, provided by the trace of the quadrupole matrix, is adequate when isotropic coordinates are preferred. The plane of the largest two eigenvalues would be the most revealing, in the case that the third eigenvalue were significantly smaller than the other two.

/* calculate quadrupole moment */

void geoqm(q,s,n) double q[3][3], *s[7]; int n; {
double en; int i, j, k;
en=(double)n;
for (i=0; i<3; i++) for (j=0; j<3; j++) q[i][j]=0.0;
for (k=0; k<n; k++) 
  for (i=0; i<3; i++) for (j=0; j<3; j++) q[i][j]+=s[k][i]*s[k][j];
for (i=0; i<3; i++) for (j=0; j<3; j++) q[i][j]/=en;
}

Once the quadrupole matrix exists, we require a means of diagonalizing it. Many diagonalization programs exist, especially for symmetric matrices. For pedagogical reasons, GEOM uses the Jacobi method [1], which is not necessarily the most effective procedure for very large matrices. Nevertheless, it is simple and reliable, so there is no reason not to use it for the matrices which we are likely to encounter.

The Jacobi method consists in using planar rotations to reduce the matrix element belonging to that plane to zero. However, trying to eliminate all the off-diagonal elements is like swatting flies. As a new element is eliminated, many of those previously destroyed partially revive, so they must be confronted over and over again.

In the history of computing methods, different strategies for deciding the order in which to attack the off-diagonal elements, as well as the validity of different approximations to the angle of rotation and the rotation matrices have been studied in considerable detail; all of this goes far beyond the requirements of GEOM.

Amongst the idiosyncracies of the Jacobi method, other than the fact that it was far too computation intensive to be practical before the advent of electronic computers, is its behavior with respect to degeneracies. Two diagonal elements which are more nearly equal than even the tiniest off-diagonal element connecting them will provoke a $45^o$ rotation trying to break the degeneracy. This can be an annoyance in symmetric molecules, where symmetry implies degeneracy and near symmetry implies near degeneracy.

The formulas required by the Jacobi method can be obtained from $2 \times 2$ matrices. Let us call $c$ and $s$, respectively, the sine and cosine of an angle of rotation, $\theta$. Then 0.4em

$$\begin{eqnarray*} \left[ \begin{array}{cc} c & s \\ -s & c \end{array} \ri... ... \\ sc(d-a)+(c^2-s^2)b & s^2a-2scb+c^2d \end{array} \right] \end{eqnarray*}$$

Some algebra and use of standard trigonometric identities shows that the condition for vanishing of the new off-diagonal element is

$$\begin{eqnarray*} \tan 2 \theta & = & \frac{2b}{a-d}, \end{eqnarray*}$$

which tells us the angle of rotation to use. The trace of the new matrix is the same as that of the old, as it should be.

However, the trace of the square of the matrix is also preserved, a bit of algebra that we would like to avoid performing in detail. In the process, the square-sum of the diagonal elements grows, at the expense of the off-diagonal elements. By the nature of the transformation, the same is true for any larger matrix in which this $2 \times 2$ submatrix is embedded, providing something of a convergence criterion.

If all the square-sum of a matrix is on the diagonal, it itself is diagonal. So at any moment we can sweep the largest off-diagonal element onto the diagonal, but without knowing the exact way the remaining elements will rearrange themselves, or which of them will be the largest remaining. So there are strategies for selecting the order in which to crush off-diagonal elements.

/* make a jacobi rotation of the rows of the    */ 
/* matrix z through angle t in plane i, j       */
/* two columns will be affected                 */

void sphjr(z,i,j,th) double z[3][3]; int i, j; double th; {
int    k;
double c, s, t;
s=sin(th); c=cos(th);
for (k=0; k<3; k++) {
  t=z[k][i];
  z[k][i]= c*t+s*z[k][j];
  z[k][j]=-s*t+c*z[k][j];
  }
}

/* make a jacobi rotation of the columns of the */ 
/* matrix z through angle t in plane i, j       */
/* two rows will be affected                    */

void sphjc(z,i,j,th) double z[3][3]; int i, j; double th; {
int    k;
double c, s, t;
s=sin(th); c=cos(th);
for (k=0; k<3; k++) {
  t=z[i][k];
  z[i][k]= c*t+s*z[j][k];
  z[j][k]=-s*t+c*z[j][k];
  }
}

/* diagonalize a symmetric matrix and record its */
/* eigenvectors (using Jacobi's method).         */
/* u[3][3] = matrix of eigenvectors              */
/* v[3][3] = matrix being diagonalized           */

void sphji(u,v) double u[3][3], v[3][3]; {
int    i, j, k, l, ll;
double e, th;
for (i=0; i<3; i++) for (j=0; j<3; j++) u[i][j]=0.0;
for (i=0; i<3; i++) u[i][i]=1.0;
  ll=0;
for (k=0, e=1.0; k<10; k++, e/=3.0) {
  l=0; do {
  for (i=0; i<2; i++) for (j=i+1; j<3; j++) {
 ll++; if (ll>5000) return;
    if(v[i][j]>-e&&v[i][j]<e) continue; 
    th=0.5*atan2(v[i][i]-v[j][j],v[i][j]+v[j][i]);
    sphjr(u,i,j,th); sphjr(v,i,j,th); sphjc(v,i,j,th);
    l=1;
    }}  while (l==1);
  }
 }