A.1 General description of the programs
The programs described here can calculate the interaction energies of either a finite set or a crystal of polarizable molecules. They were written by Mihaly Mezei in collaboration with Prof. Edwin S. Campbell at NYU in 1975-76 and are based on algorithms developed either jointly or earlier by Prof. Campbell. They have been converted to the current Fortran 77 version from 1988 to 1992 at Hunter College, CUNY and at Mount Sinai School of Medicine, NYU. It now runs on SGI, IBM, VAX/VMS, CRAY and on any generic Unix system. Inquiries can be addressed to Dr. Mihaly Mezei at Department of pharmacological Sciences at the Icahn School of Medicine at Mount Sinai, New York, NY 10029 or by electronic mail at Mihaly.Mezei@mssm.edu
The package calculates the energy of interaction of a set of charge distributions (either an arbitrary finite collection or arranged in a crystal lattice) characterized by a multipole expansion. Cooperativity is introduced through the dipole polarizability tensor of each charge distribution.
The distinguishing feature of the package is the use of the formalism of Maxwell as described by Hobson [1], which computes the interaction energy of multipoles in the form of general directional derivatives.
The calculations start with the generation of the Taylor-series expansion coefficients for the charge distribution. These coefficients can be computed for a one-determinant wave-function expanded over a Gaussian basis set with the program MOMENTS based on algorithms developed Mezei and Campbell [3] (or obtained from elsewhere, e.g., form the ab-initio package GAUSSIAN-xx). The program MOMENTS can perform the multipole expansion at the position of any of the atoms or it can partition the electron density into atomic contributions in a number of different ways and the perform the multipole expansion on each separately. If necessary, the moments can be translated and rotated by the program MOMTRNSF.
The poles (scalar multipliers for each harmonic order) and characteristic directions (used in the directional derivatives) introduced by the Maxwellian formalism are calculateded from a set of Taylor series coefficients with the program CHARDIR using the algorithm of Mezei and Campbell [2].
The program MULTIPOL computes the interaction energies between a finite set of molecules from three types of contributions:
The same three types of contributions to the lattice energy of a crystal of polarizable molecules can also be computed. The permanent multipole contribution is calculated with Campbell's generalization of the Ewald method for ionic lattices to lattices of multipoles of arbitrary order and an efficient algorithm for a sequence of calculations on the same lattice but different molecular orientations [7]. A correction to one of the equations in [7] is derived in Appendix I.) The additional inverse-power pairwise additive terms are simply calculated on a finite lattice. The calculations involve three different programs:
The program assigns each data set an identifier (ID) number that in some cases is specified by the user and in other cases by the ID number of the input data set. This allows the generation of multiple data sets in a single run and their selective utilization in subsequent runs. Current input format restricts the ID numbers to three digits.
The input format is specified using standard Fortran FORMATs. For example, (24I3,A4) describes the input of 24 integers of 3 digits in one line (right-adjusted), followed by a string of 4 characters; (5F8.0,10X,E12.5) describes the input of five real numbers, occupying 8 spaces each, with the decimal point explicitly given, followed by 10 spaces to be skipped and after that a number in scientific notation, occupying 12 spaces with a 5-digit mantissa, e.g., -1.76921E-03. The appearance of / in the format instructs the input to start to read a new line.
1. General description
This program computes the permanent multipole moments of a charge distribution using the formulae of Mezei and Campbell [3]. All notation used in this documentation is defined in [3] end the equation numbers used also refer to [3]. The moments may be calculated from the whole charge distributions (one-center expansion) or the density can be partitioned into atomic contributions (split expansion).
In the simplest case of a one-center expansion, the set of moments of the entire charge distribution is computed with respect to a single atomic center. Alternatively, a sequence of such sets of moments may be obtained, one for each of a set of atomic centers. Each set after the first is obtained by the simple transformation of the moments of the first set for the new center of expansion.
Partitioning of the density to atomic contributions only affects contributions from two different centers (overlap density). These are partitioned by either sharing it between the two centers it originates from in some proportion or separating each overlap density by a plane perpendicular to the line connecting the atoms involved (called sharp split). Section 2.1 of [3], p. 228-9, gives the multipole moments and splits of the charge density which may be used; Appendix A, p.234, gives the forms assumed for the contracted Gaussian wave functions; Appendix B, p.234-6 gives the integration formulae used.
For the sharp split the program contains a generalized option in which the function ro(p,q), which is determined by input w(p) and w(q) is replaced by a function of the Gaussian exponents of the two centers, ro(ALPHA(p),ALPHA(q)). This option places the cutoff plane at a distance XDC/ICON(12) from the center of the product Gaussian where XDC is the distance between the center of the product Gaussian and the center of the farther individual Gaussian and ICON(12) is the 12th input control constant. However, in the tests thus far made it showed no advantage over the simpler ro(p,q).
The output can consist of one or both of the following sets of moments:
At the end of the input the molecular orbitals are checked for orthonormality. This has proved extremely useful in detecting input errors.
2. Detailed description of the input
These lines are sequenced (unless IFLAG or IGNOR is non zero) with the number, II, of the m.o. and the number, MMM, of the line. Lines read according to 8.1.a-8.1.e are not sequenced. IT and MMM are in cols. 65-68 and 73-80 respectively. Repeat 8.1.c and 8.1.d for each m.o.
If ICON(14) = 1 or 2, the m.o. coefficients are given in the GAUSSIAN-xx format. The input could be copied directly from the output of GAUSSIAN-xx. The m.o. coefficients in groups of 5. For each grpup, the following data are read:
The program takes the standard input and output files as unit 5 and 6. If ICON(2)>0, the program puts out all the moments sets generated on unit IUMOM in binary form:
WRITE (IUMOM) (ILBL(J),J=1,20),(ISPEC(J),J=1,10) (description) WRITE (IUMOM) ICNTR(I),I,ISPLIT,N,NID,ICON(1) (NID: set i.d.) WRITE (IUMOM) (NLMN(II,JJ),JJ=1,3),II=1,N),(PLMN(II),II=1,N)
5. Subroutines used
The program is based on the PA60 module (properties) of POLYATOM and several subroutines are taken from there. The order of the spherical harmonics in the primitives has been changed to allow the use of the general polynomial coefficient handler routines.
1. General description
This program inputs a set of cartesian multipole moments (calculated, e.g., by the program MOMENTS):
M(RHO,xo,n,{ei}) = INTEGRALD PRODUCT3k=1 (xk-xko)^nk RHO(x) dx
where
2. Detailed description of the input
ICON(19) lines, each of which contain the name and the coordinates of one center:
NAME,X,Y,Z (A4,8X,3F12.0)
3. Output of the program
The program output consists of a list of the input data and of the moments which have been calculated. If ICON(2) = 1, then the moments will also be written on unit 30 in binary form (the format is identical to the format of unit 20 written by MOMENTS) and the ID numbers of each moments set will appear on the printed output. The output format is the same as for the program MOMENTS.
4. Subroutines used
1. General description
This program computes the multipole moments and characteristic directions of a spherical harmonic representing the cartesian moments of a charge distribution as required by the Maxwell theory of the poles of a spherical harmonics [1]. The formulae for the determination of the poles and characteristic directions are given by Mezei and Campbell [2].
For a homogeneous polynomial of degree n, Yn(x,y,z) satisfying Laplace's equation the multipole moment Pn and characteristic directions (unit vectors) s1,...,sn are determined to satisfy
Yn(x,y,z)/r2n+1 = [(-1)n * Pn /n!] * PRODni=1 (si.DEL) (1/r) |r=(x,y,z)
The program calculates the poles and characteristic directions for a set of Yn's, n=2...,N, obtained as the cartesian moments of a charge distribution. See the program MOMENTS for details.
2. Detailed description of the input
IUNFMI, IUNFMO, NREP, NOFINT, MAXORD, NID, IECHO
NOFINT lines: C, L, M, N (E15.9,5X,3I2)
where C is the coefficient of the term xLyMzN in the surface spherical harmonic of order L+M+N.
MXORD, NID, IUNFMI (3I5)
Their meaning is the same as read on line 2. Zero value for IUNFMI will retain the current value. The input continues with item 3. The calculation stops when no more data are given.
The program echoes the data read. The calculation proceeds by spherical harmonic order. For each order it echoes I(L,M,N), the coefficient of xLyMzN read, before dividing by the corresponding factorial. During the calculation some partial results are printed (in the notation of the ref. [1]). The algorithm of [1] first determines the characteristic direction set. The scalar poles are obtained as the ratio of Yn calculated from the characteristic direction set with the multipole moment set to one (called MAXWELL on the printout) and its value calculated directly the inputted coefficients (called TAYLOR on the printout) at NREP randomly chosen points. The program prints the NREP MAXWELL and TAYLOR values, the calculated value of the pole. ERROR is the calculated standard deviation of the ratios.
If IUNFMO > 0, the poles and characteristic directions are also written on unit IUNFMO with unformatted write to be read by other programs.
4. Files
The standard input and output are on units 5 and 6, respectively. If IUNFMI > 0, the moment set(s) are read from unit IUNFMI in binary form (as written by the program MOMENTS). If IUNFMO > 0, the calculated poles and characteristic directions will also be written to unit IUNFMO in binary form:
WRITE (IUNFMO) (ITITL(I),I=1,20),0,MXORD,MAXINP,NID
WRITE (IUNFMO) (POLE(I),I=1,(MXORD+1)) (the poles)
WRITE (IUNFMO) (S(I),I=1,(3*MXORD*MXORD+1)/2)) (the components
of the characteristic directions)
5. Subroutines used
1. General description
The program computes permanent multipole, and, optionally, permanent multipole-induced dipole, induced dipole-induced dipole interaction energies, supplemented with inverse power terms (called correction potential) used to represent repulsion and dispersion interactions. Electrostatic energy calculations use the Maxwellian formalism for the multipole expansion. The induced dipoles are computed according to ref. [5]. The virial sum, as defined in [11] can also be calculated.
The system is assumed to consist of a set of charge distributions (molecules). Each charge distribution is characterized by a set of sites (corresponding to the atoms of a molecule). Each site may be the center of a multipole expansion but only one of them may be polarizable. The program inputs a set of charge distribution descriptions, referred to as distribution types, given in a principle polarizability axis frame of reference local to that charge distribution (molecule). The orientation of the molecule in the global frame will determine a rotation to be applied to these charge distributions (i.e., to the characteristic directions). The same type of charge distribution can be used in different orientations, e.g., for the characterization of two water molecules in different orientations. Note, however, that the chemical equivalence of two atoms in a molecule does not necessarily imply that they can be described with the same charge distribution type.
The calculation of the permanent multipole energy will include spherical harmonics with orders less than or equal to the lower between the highest orders of the two center. This way all order interactions will be complete. Note that for this comparison the quantity MAXORD will be considered (see 5.2.1 below) that may be higher than the actually existing terms. The necessity for this may arise, for example, when in certain split density calculations it is possible that some center may have vanishing higher order moments, but the interaction of the low order (non-vanishing) moments of this center with higher orders of other centers are still sought for.
The Maxwellian poles and characteristic directions can also be supplemented with a set of Monopole fields (see CHARDIR) to describe the charge distribution belonging to that center. Note that if the partial derivatives in a Taylor series expansion are treated via the Maxwellian formalism [1] then for the pole corresponding to xLyMzN factor the characteristic direction set will consist of L copies of the unit vector in the x direction, M copies of the unit vector in the y direction and N copies of the unit vector in the z direction.
2. Detailed description of the input
The program reads information from several different files. The general information on the run is on unit 5. The description of the charge distributions is on unit 7. For a special case (water) the program can read the global cartesian coordinates of the atoms from unit 8.
The permanent multipole energy calculation can optionally be limited to orders lower than MAXINT for distant centers. When MXORDR(I) is nonzero, then for each pair of centers whose distance R falls between (I-1) Å and I Å the expansion will be limited to spherical harmonics of order MXORDR(I). For R > 10 Å, MXORDR(10) is used (or MAXINT when MXORDR(10)=0).
Interactions between molecules whose distance is more than RABMAX is set to zero. If ICINP=0, the distance between group centers will be used, if ICINP>0 the distance between the first atoms of each molecule will be used.
The principal axis dipole polarizabilities of the I-th charge distribution or zeros. For a molecule (group) (cf. geometry specification) consisting of more than one center, only one of them is allowed to have nonzero polarizabilities.
If ICDBIN > 0 , read items 5.3.4 - 5.3.6 from unit ICDBIN (written by the program CHARDIR): the binary poles and characteristic direction set with NID as the identifier number.
If ICINP=1 or 2, the input stops here (and the configurations are read from unit 8 - see input item 8 below), followed by the energy calculation.
WRITE (8) NGR,IDENT1,NCONF,NMCMIN,IPBC,EX,EY,EZ,EPWA,VPWA
WRITE (8) ((R(I,J), I=1,3), J=1, 3*NGR)
Next, the program performs the energy (and virial sum) calculation and proceeds to read the specification of an other configuration set (7 or 8 above) until MAXDAT energy calculations have been performed (or the calculation was stopped by NGR=0 in 7 above).
The output printed on unit 6 is as follows:
REPEAT from 3 for the next configuration.
4. Files
Optionally (ICINP=1), the configurations can be read from unit 8. The program puts on unit 9 a summary of the energies calculated. At continuation of the run, this file will be appended. The data is written as
WRITE (9) NEDATA,ERSLT,NMC,EPWA,VRSLT,VPWAwhere NEDATA is the number of configurations calculated and the arrays ERSLT(8,##D),VRSLT(4,##D),EPWA(##D),VPWA(##D),NMC(##D) contain the following information:
The following two subroutines are specific for water-like (c2v) molecules:
The program contains special symbols that are to be replaced by actual numbers representing the maximum dimensions of arrays, before compilation. These symbols are given below.
1. General description
This is the first in a set of two programs for the calculation of the permanent multipole and the induction contributions to the lattice energy. It prepares of the charge-distribution-independent sets of quantities, called crystal constants KNT(NU,Xc-Xj). Campbell's algorithm [7] divides the lattice sites into primitive translation lattices, which contain only corner sites. The sites in the unit cell are indexed and each primitive translation lattice is indexed by the site in the unit cell, which belongs to it. Here the superscript c is the index for the unit cell site and j is the index for the translation lattice, but both Xc and Xj are in the same unit cell. The crystal constant have to be calculated only for crystallographically nonequivalent site-lattice pairs.
As input, it requires:
The program prints the computed crystal constants as well as the two lattice sums (K1 and K2 in the notation of ref. [7]). It can also write the computed nonzero crystal constants and a sequence number onto a file in binary format.
2. Detailed description of the input
The summation of the transformed quantities is performed on a finite lattice where the lower and upper bounds of the indices are given by LiMIN, LiMAX, i = 1,2,3, respectively. Each index set {ik} represent a site <i1a1, i2a2, i3a3> where ak are the crystal basis vectors (see 5 below). The central cell is assumed to have all zero indices.
The program echoes the input information and for each translation lattice and for each nonzero crystal constant the program prints their sequence number, the components of the corresponding NU vector, the crystal constants over both the reciprocal and direct lattices, and their sum. In addition, it prints:
The (nonzero) crystal constants are written on unit 20 by the following statements:
WRITE (20) ITITLE(I),I=1,20) (alphanumeric title)
WRITE (20) NSITE,NLATT,EPS,NID,MXO,MAXLMN
WRITE (20) (RT(I),I=1,3)
WRITE (20) (KNT(J),J=1,MAXLMN) (crystal constants)
Here MAXLMN is the total number of crystal constants with orders
not exceeding MXO.
5. Subroutines used
1. General description
This program computes the permanent multipole energy of the crystal and/or the induced dipole vectors and the induction energy. A cooperative calculation may be followed by a second non-cooperative calculation. Input is assumed to be in atomic units of distance and charge. Single center multipole expansions are assumed. Its input consists of the following information.
The characteristic directions, s, are assumed to be in the principal axis system of the charge distribution (a local system) and are transformed with the rotation matrix R as defined in Appendix III.
Since the same type of charge distribution may be located at more than one site in the unit cell, an array is read in which specifies the index number for the type of charge distribution at each site. A second array specifies the index of the rotation matrix for each site. The local systems are assumed to be principal axis systems for the polarizability tensors.
K(NU,Xi-Xi)= K(NU,Xj- Xj)
K(NU,Xi-Xj)= (-1)^(NU1+NU2+NU3) K(NU,Xj-Xi)
Furthermore, identities between crystal constants are used to generate identities between the site-lattice permanent multipole energies, Up(Xc,Tj), and between the permanent multipole electric fields at sites Xc, by lattice Tj, Ep(Xc,Tj), as follows. Let C(i) be an index array such that C(i)=C(j) if and only if the type of the charge distributions and their orientations are identical.
If for two pairs <i,j> and <i',j'> ,
<C(i),C(j)> = <C(i'),C(j')> and Xi-Xj = Xi'-Xj' ,
then
Ep(Xi,Tj) = Ep(Xi',Tj') and Up(Xi,Tj) = Up(Xi',Tj') .
If for two pairs <i,j> and <i',j'>
<C(i),C(j)> = <C(j'),C(i')> and Xi-Xj = -(Xi'-Xj')
then
Up(Xi,Tj) = Up(Xi',Tj') .
After the crystal constants are read in and the redundancy search has been performed, the user has the option of providing further identities between the Ep and the Up. Since the crystal constants and the information about identities are stored on a disk file, for repeated calculation this procedure can be bypassed. Both the identities implied by the above relations and the identities read are reported in the printout.
The next part of the program permits the optional calculation of the induced dipole vectors and their interaction energies with the lattice permanent multipoles. If a cooperative calculation is performed, then a second non-cooperative calculation can also be performed.
If IUNFCD=0 then the poles and characteristic directions are read in 9c and 9d:
If IUNFCD > 0 then unit IUNFCD will be searched for the characteristic direction set with identifier NID and read from there (in the format specified by CHARDIR).
For each NU, K(NU,XNSITDUP - XNLATDUP) is computed as
K(NU,XNSITDUP-
XNLATDUP) =
K(NU, XNSITORG-
XNLATORG)
(-1)^(NUj* KTR(ITR,j))
where KTR is an internally stored matrix:
| 0 0 0 |
| 1 0 0 |
| 0 1 0 |
KTR = | 0 0 1 |
| 1 1 0 |
| 1 0 1 |
| 0 1 1 |
| 1 1 1 |
Set ITR equal to the row index for the row such that KTR(ITR,j)=1
if NUj is odd and 0 if NUj is even. [It can
be shown that a change
in the sign of the j-th component of RT introduces a factor of -1
if and only if NUj is odd.]
The program echoes the data read and prints the calculated energies, electric field components and induced dipoles, both their sum and their breakdown to sites, and, wherever applicable, site-lattice contributions.
Optionally (when ITABPR=1) the program saves the calculated energy and (when INDC=1) field tables on unit 30 either (when IUNFTB=0)in alphanumeric format or (when IUNFTB=1) in binary format with the following write statements:
WRITE (IUTAB) (ITITLE(I),I=1,20) DO ISITE=1,NUNIT WRITE (IUTAB) (SLENRG(ISITE,ILATT),ILATT=1,NUNIT) IF (INDC .EQ. 1) THEN WRITE (IUTAB) ((SLFLD(ISITE,ILATT,K),K=1,3),ILATT=1,NUNIT) END DO
5. Subroutines used
1. Introduction
This program computes lattice sums of inverse powers of distances between sites which can be used to supplement the permanent multipole and induced lattice energies with a contribution of the form:
Each molecule is described by the positions of it sites in its 'local' frame and the type index of each site. Two atoms are considered of identical types when they can be considered chemically equivalent and thus they contribute to the same inverse-distance power sums.
The orientation of the molecules is described by rotation matrices.
A set of exponents is to be specified for which the summation will be performed. For each exponent -n and each pair of site types A and B, the program computes the matrices
where type(p)=A and type(q)=B
A set of coefficients, dn<A,B>, may also be given. In this case the lattice energy contribution
will also be computed.
2. Detailed description of the input
The program echoes the data read and prints the calculated sums over the lattice for each exponent and for each pairs of atomtypes.
4. Subroutines used
The programs are written in Fortran-77. A c-shell script compile.csh is provided for the compilation of all programs on Unix systems. Before execution, compile.csh has to be edited to specify the compiler name and compilation options ( vide infra).
For non-UNIX systems the following steps are needed:
The program MULTIPOL requires a preprocessor (see Sec. E.6) to set the dimensions as needed. For Unix systems the c-shell script multipol.csh (automatically called by the compile.csh script). Typing multipol.csh multipol will read the file multipol.for containing the special chartacters and create the file multipol.f that is ready to compile. When multipol.csh is executed, the terminal output will identify the dimensions given to multipol.f. Without change, the following values will be used:
The current version uses single precision arithmetic in most places. This limited the precision of the calculated moments to 4-5 significant figures for a water molecule in a calculation with 32-bit arithmetic. The most sensitive part of the package is the calculation of the so-called poles by the program CHARDIR. If the errors printed for the calculated poles is unacceptable, or if increased precision is needed for othe reasons, the programs should be compiled with double precision. The way to do is compiler dependent - on most systems the -r8 or -d8 compiler option doubles the word sizes. Note that if one program is compiled with double precision, all the others have to be as well (since they communicate via binary files).
Appendix I. The Correction for the Self-lattice Sum of a Directional Derivative
Reference [7] derived the Eqs. (A.11,12) for the correction and the self-lattice sums of directional derivatives. However, it did not give the decomposition required by the algorithm for the study of a larger number of orientations:
RNTHETAo,2({sk},EPS) = SUM{NU} RNTHETAo,2(NU,EPS) SIGMA({sk},NU). (A1.1)
Since in Eq. (A.11) j = q = N/2, Eqs. [(38), (28d), (36), (37b,c), (23b)] imply:
RNTHETAo,2(NU,EPS) =
{[(-1)(N/2)+1 2 EPSN+1] /
[PI1/2(N+1)]}
PROD3i=1 [NUi!/(NUi/2)!]. (A1.2)
Finally, in Eq. (40) for the self-lattice sum of a directional derivative, RNTHETAo,2(NU,EPS) should replace RNTHETAo,2({sk},EPS). In addition, the equation references should read:
KNTo,1, (25b,b.1-4); KNTo,2, (39a,a.1-4).
Appendix II: Subroutines for handling polynomials of three variables
The polynomial of degree n is assumed to be of the form
SUML+M+N < n cLMNxLyMzN . (A2.1)
The coefficients cLMN are stored in a linear array, sorted by n. The number of terms of degree i is (i+1)*(i+2)/2. The length of this array is n*(n*(n+3)+2)/6, obtained from the identity
SUMni=0 i2 = (2n3+3n2+n)/6 . (A2.2)
The position of cLMN is
n'*(n'*(n'+3)+2)/6 + N*(2*n'-N+3)/2+M+1 (A2.3)
where
n'=L+M+N . (A2.4)
SUM cLMN (x-xo)L (y-yo)M (z-zo)N = SUM c'LMN xL yM zN . (A2.5)
(x-xo)L (y-yo)M
(z-zo)N =
r' = R r (A2.7)
where
Note that the LMN term of the polynomial in r' becomes
( SUM3i=1 R1ir1)L ( SUM3i=1 R2ir2)M ( SUM3i=1 R3ir3)N (A2.8)
and
( SUM3i=1 R1iri)L = SUM [L!/(l!m!n!)] R11l R12m R13n r1l r2m r3n (A2.9)
A matrix with one row will be denoted by () and a matrix with one column by []. Let ej(i) and e j(f) denote the j-th basis vector in the initial (i) and the final (f) frames, respectively. Let the matrix B, which transforms the initial set of basis vectors into the final set, be defined by the following equation:
(e1(f), e2(f), e3(f)) = (e1(i), e2(i), e3(i)) B . (A3.1)
Let xj(i) and xj(f) denote the j-th component of any vector x with respect to the initial and the final frames, respectively. Eq. (III.1) implies that the components of x with respect to the final frame transform into the components with respect to the initial frame according to the following equation:
x1(i) x1(f) [x2(i)] = B [x2(f)] (A3.2) x3(i) x3(f).The following equivalent equation for the inverse transformation of the components of x with respect to the initial frame into those with respect to the final frame is:
x1(f) x1(i) [x2(f)] = BT[x2(i)] (A3.3) x3(f) x3(i).This assumes that the frames are orthogonal so that the transpose, BT, is the inverse. Thus, the transpose, BT, is the matrix, which this program uses to transform the components of the characteristic direction vectors in the initial frame (the local frame) to the components in the final frame (the common "global" frame). The matrices are most conveniently constructed using the following equation:
(em(f).ek(i) = Bkm = (BT)mk. (A3.4)
Thus, the elements of the k-th column of BT are the components of the initial frame ek(i) in the final frame.
Appendix IV. The Definition Adopted for Eulerian Angles
Since current conventions for Eulerian angles use different orders of the successive rotations and/or notation, this appendix gives the convention adopted in these programs. In any comparison of the discussion here with other treatments, note that some authors use a row and others a column matrix for basis vectors. The notation of Appendix III for basis vectors and the components of a vector will be adopted.
Consider one of the three coordinate planes for the initial set. Then the corresponding plane for the final set is either coplanar with it or the initial and final planes intersect in a line called the line of nodes. The definition of the angles is, of course, based on the more general non-coplanar case. The two following initial choices are made here.
B1 = the rotation matrix which transforms {ej(i)} into {ej'}. (A4.1)
Since both e3'=e3(i) and e3(f) are perpendicular to e1', rotate e3(i) about the e1' axis by a positive angle THETA to e3(f). This gives the second set of intermediate basis vectors, {ej"}, such that e3"=e3(f). Let:
B2 = the rotation matrix which transforms {ej'} into {ej"}. (A4.2)
Finally, rotate e2" = e1' about the e3(f) axis by a positive angle PSI to e1(f). Let:
B3 = the rotation matrix which transforms {ej"} into {ej(f)}. (A4.3)
Each axis system is assumed to be right-handed and orthogonal and each rotation is counter-clockwise by a positive angle. (If mnp is any cyclic permutation of 123, the analytical conditions are:
| cos (ALPHA) | -s( (ALPHA) | 0 | ||
| C = | s( (ALPHA) | cos (ALPHA) | 0 | (A4.4) |
| 0 | 0 | 1 |
| 1 | 0 | 0 | ||
| D = | 0 | cos (ALPHA) | -s( (ALPHA) | (A4.5) |
| 0 | s( (ALPHA) | cos (ALPHA) |
i j 1 1 cos PSI * cos PSI - sin PSI * cos THETA * sin THETA 1 2 sin PSI * cos PSI - cos PSI * cos THETA * sin PSI 1 3 sin THETA* sin PSI 2 1 -cos PSI * sin PSI - sin PSI * cos THETA * cos PSI (A4.6) 2 2 -sin PSI * sin PSI + cos PSI * cos THETA * cos PSI 2 3 sin THETA* cos PSI 3 1 sin THETA* sin PHI 3 2 -cos PSI * sin THETA 3 3 cos THETA .