Documentation of the program package MAXWELL for the
calculation of electrostatic interaction energies using multipole
expansion in the Maxwellian formalism.
Mihaly Mezei
Department of Pharmacological Sciences,
Icahn School of Medicine at Mount Sinai
New York, NY 10029
Email address:
Mihaly.Mezei@mssm.edu
and
Edwin S. Campbell
Department of Chemistry,
New York University,
New York, NY 10003
Nov. 15, 2004.
CONTENTS:
 A. OVERVIEW
 B. MOMENTS: calculating Taylor expansion coefficients
 C. MOMTRNSF: transforming Taylor expansion coefficients
 D. CHARDIR:
determinig the Mawell poles and characteristic directions
 E. MULTIPOL:
calculating the energy of a system of multipoles
 F. CRYSCON:
calculating tha lattice contribution to the energy of a crystal of multipoles
 G. CRYSTAL:
calculating the energy of a crystal of multipoles
 H. CRYSPOT:
calculating general nonbonded energies of a crystal
 I. INSTALLATION
 I. APPENDICES
 J. REFERENCES
A. OVERVIEW
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 197576 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
Taylorseries expansion coefficients
for the charge distribution.
These coefficients can be computed for
a onedeterminant wavefunction 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 abinitio package GAUSSIANxx).
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:
 (i) The
permanent multipole interaction energies are calculated with
recursion algorithms, which were developed by
Campbell and Mezei [4].
 (ii) The nonadditive induction energies are calculated in the
induced dipole approximation with equations derived by Campbell
[5].
 (iii) Additional pairwise additive interaction energies,
which have the
form of a linear combination of inverse powers of the distance
between pairs of centers on each pair of molecules can be used to
model dispersion and repulsive contributions. These algorithms were
used to develop a cooperative potential function for waterwater
interactions [6].
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 inversepower pairwise additive terms are simply
calculated on a finite lattice.
The calculations involve three different programs:
 (a) The program CRYSCON calculates
for each spherical harmonic order a set of quantities, called
crystal constants (depending only on the geometry of the
lattice sites but independent of the charge distribution)
 (b) The program CRYSTAL calculates the
electrostatic (permanent moment and induced dipole energy)
using these crystal constants and the characteristic directions
(generated by CHARDIR) in a
much shorter calculation as decribed by
Campbell and Mezei [4].
 The program CRYSPOT computes the remaining
pairwise additive contributions for each molecule in the unit cell.
It sums the
interactions over all other molecules in a crystallographic
parallelepiped about each molecule. The potential used has the form
of inverse powers of the distance between the interacting centers.
A.2. Data flow
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 (rightadjusted), 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 5digit
mantissa, e.g., 1.76921E03. The appearance of / in the format
instructs the input to start to read a new line.
B. MOMENTS
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 (onecenter
expansion) or the density can be partitioned into atomic
contributions (split expansion).
In the simplest case of a onecenter 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. 2289,
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.2346 gives the
integration formulae used.
For the sharp split the program contains a generalized
option in which the function r_{o}(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,
r_{o}(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
r_{o}(p,q).
The output can consist of one or both of the following sets of
moments:
 For single center expansion,
the moments for the entire distribution with
respect to a preselected atomic center
and, if requested, with respect to each atomic center.
 For split expansion,
the split moments of the entire charge distribution
and, if requested,
the moments for the
entire distribution , centered at a preselected atomic center.
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
 Title of the problem (20A4)
 Control constants ICON(I) (24I3)
 ICON(1): The highest order moment to be computed,
i.e., the maximum n_{1}+n_{2}+n_{3}
in Eq. 3c.
 ICON(2) 1(0): Print (do not print) the calculated
moments on unit
IUMOM
as specified by ICON(16)
in binary form (for communication between programs).
 ICON(3): When ICON(6)=0, i.e., split densities are
calculated, this determines the source of the split parameters
S(k,p,q) or w(c).
 ICON(3)=0: Set the split constants to define the Mulliken
split (see Eq.(10))
 ICON(3)=1: Read in the constants defining the split. In particular,
 If ICON(10)=0, read the S(k,p,q) (input item 4a).
S(k,p,q) gives the fraction of the overlap density
between centers p and q to be shifted to center k.
 If ICON(10)=1, read the w(c) (input item 4b).
The overlap density between centers p and q will be shared
between the two centers in the proportion w(p)/w(q).
 ICON(4) 0(1): Compute (do not compute))
the n_{j} of Eq. (A2a).
 ICON(5) 1(0): Read (do not read) the linear combination
coefficients a^{d}_{m}
which specify the type of hybrid (see item 5
below). If they are not read in, each primitive is assumed to be
pure (unhybridized), i.e., a Gaussian multiplied by a single
polynomial term in x,y and z.
 ICON(6)
 = 0: Compute a split expansion
(see ICON(3), ICON(7), and ICON(10) for additional options);
 > 0: Compute a single center expansion about the ICON(6)th
center see ICON(7) for additional option).
 ICON(7)
 When ICON(6) > 0:
 = 0: calculate the moments with respect to center ICON(6)
only;
 = 1: transform the moments with respect to center ICON(6)
to obtain moments with respect to each of the other atomic centers.
 When ICON(6)=0 (split expansion):
 = 0: Do not contract the split moments;
 > 0: contract the split moments to single center moments
about the ICON(6)th center.
 ICON(8): ID (identifier) number
of the first moment set to be
written on unit
IUMOM.
Subsequent sets' ID numbers are incremented by
one each.
 ICON(9) 0(1): compute (do not compute)
the N^{j}_{q} of Eq. A2 to
renormalize the primitives.
 ICON(10) 0(1):
split the density using the S(k,p,q) (the w(c))
 i.e., share the overlap densities in proportion defined by the
S's or w's. Used only when ICON(6)=0 and ICON(3)=1.
 ICON(11)
 = 0: real numbers < 10^{10} are to be regarded as
zero.
 > 0: real numbers <10^{ICON(11)} are to be
regarded as zero
 ICON(12)
 = 0: perform a sharp split using r_{o}(p,q).
 > 0: perform a sharp split using
r_{o}(ALPHA(p),ALPHA(q))
(cf. section 1 of this documentation.)
 ICON(14) 0(1,2): POLYATOM (GAUSSIANxx old, new) format
for molecular orbital (m.o.) coefficients.
 ICON(15) 0(1): Atomic unit (Å) units assumed for the
input coordinates.
 ICON(16): Unit number to write the calculated moments
(IUMOM) in binary form (when ICON(2)=1).
If ICON(16)=0, unit 20 will be used
 ICON(17) 0(1): Atomic units (Debye  Å^{N}) units
will be used to print the calculated multipole moments
 ICON(18) > 0: Echo each input line annotated with the
item number according to the documentation
 ICON(22)
 = 0 : Onedeterminant wave function.
 Lines specifying the center coordinates:
 3a. One line giving the number of nuclei, NOC (I3).
 3b. NOC lines (A4,8X,4F12.0) each of which contains the
symbolic name of a center, its X,Y,Z coordinates and
its nuclear charge (in this order).
 If ICON(3)=1: Split expansion information
 4a. When ICON(10)=0:
S(k,p,q) in format (14F5.0) in the following order:
all data for k=1, then all for k=2, etc.
Data for each k begins a new line.
For each k, the order is:
S(k,2,1),S(k,3,1)...S(k,n,1); S(k,3,2),
S(k,4,1),...S(k,j,1) ... S(k,j,j1), ...
 4b. When ICON(10)=1: the w(c)'s in format (14F5.0).
 If ICON(5)=1: lines defining hybrid orbital coefficients
(a^{d}_{m}).
This allows the definition of primitives where the Gaussian
is multiplied by a multiterm polynomial in x,y, and z.
 5a. One line (I3) with NOT, the number of types of functions
(i.e., number of different polynomials).
 5b. NOT sets of one, two or four lines.
 5b1. The first line gives ITYPE,NOTY,CS,CX,CY, and CZ in
(A3,3X,I3,3X,4F12.0) format.
ITYPE is the name of the hybridization type. NOTY is the order
of the highest order polynomial term with nonzero coefficient (0,
1, 2 or 3 for types which have at most s, p, d, or f character,
respectively.
 5b25b4. The second, third and fourth lines (as needed,
depending on NOTY) in format (6F12.0) give the CXX, CYY, CZZ, CXY,
CXZ, CYZ, CXXX, CYYY, CZZZ, CXXZ, CXYY, CYYZ, CXZZ, CYZZ and CXYZ.
CS,CX, ...,CXYZ are the coefficients of the unnormalized pure s, x,
y,...,xyz functions in the hybrid type.
 Lines giving the basis function specifications.
Both POLYATOM and GAUSSIANxx print the coefficients and
exponents of the Gaussian primitives with the convention that first
both the primitives and the individual basis functions are to be
normalized. This will be achieved if the control constants 4 and 9
are left as 0.
 6a. N,NBFNS (2I3):
 N: the total number of primitives
 NBFNS: the total number of basis functions
 6b. One or more lines (36I2) giving the number of primitives
per basis function for each basis function.
 6c. N lines (A6,6X,A3,9X,2F12.0) each of which contains a
gaussian primitive specified by symbolic reference to its center
(CNTR) and type (TYP), and with numerical values for its exponent,
and coefficient in the basis function. The symbolic reference to
the function types (defined in 5a) must be leftadjusted in the
field and correspond to one of the types defined above. To set the
exponents and coefficients of a given basis function equal to those
of a previously readin function (complete), one line (A6,6X,A3,I3)
is needed with CNTR, TYP, and INC, the latter being the number of
the previous function.
 If ICON(14) is zero (POLYATOM m.o. format), lines giving
the symmetry specifications. (This item is not used in the moment
calculation but is included as part of the POLYATOM m.o. format.)
 7a. One line (I3) with the number, NSYMTP, of symmetry types.
(e.g., SIGMA, PI, DELTA, ...).
 7b. One line (I3,3X,A6) with the number of molecular orbitals
of a symmetry type and the name of the symmetry. Repeat this line
for each symmetry type. The sum of the numbers read in on these
NSYMTP lines must equal the number of m.o.'s read in as given
below.
 Molecular orbital coefficient matrix.
 8.a. The label of the molecular orbital, (m.o.), matrix
(20A4)
 8.b. NR, NC, WA, WB, WC, WD, IT, IFLAG, and IGNOR.
(2I3,4X,3A6,2X,A6,4X,I6,4X,I6)
WA, WB, WC, WD, IT, IFLAG, and IGNOR are only used when
ICON(14)=0.
 NR : number of mos (rows) in the m.o. matrix.
 NC : number of basis functions (columns) in the m.o.
matrix
 WA, WB, WC : 18 characters with the name of the matrix,
(e.g., closed orbitals, alpha orbitals, ... , etc.)
 WD : 6 characters with the name of the program which
produced the matrix.
 IT : the number of the iteration at which the matrix was
produced.
 IFLAG : a flag which if non zero causes the line
sequencing to be ignored and causes the individual m.o. label lines
not to be read in.
 IGNOR : a flag which if non zero causes the line
sequencing to be ignored.
If ICON(14) = 0, the m.o. coefficients are given in POLYATOM
format.
For each m.o., read in items 8.1.c8.1.e:
 8.1.c. If IFLAG = 0:
 II, K, NAM, and OE(II) (2I4,A6,6X,E15.8):
 II = the number of the m.o., the row number.
 K = the number of the m.o. within its symmetry type.
 NAM = 6 characters with the name of the symmetry type to
which the m.o. belongs.
 OE(II) = the eigenvalue (orbital energy) of the m.o.
 8.1.d. One or more lines (4E15.8) with the coefficients of
the IIth m.o.
(These are coefficients of normalized basis functions.)
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.a8.1.e are not sequenced. IT and MMM are in
cols. 6568 and 7380 respectively. Repeat 8.1.c and 8.1.d for each
m.o.
 8.1.e. One line (20A4) with the same label as in 8a.
If ICON(14) = 1 or 2, the m.o. coefficients are given in the
GAUSSIANxx format.
The input could be copied directly from the
output of GAUSSIANxx.
The m.o. coefficients in groups of 5.
For each grpup, the following data are read:
 8.2.a. Column number  just skip a line
 8.2.b. Symmetry class  just skip a line
 8.2.c. Eigenvalue  just skip a line
 8.2.d. if ICON(14) = 1 then
NR lines with the old Gaussian format (20X,5F10.0), containing the
coefficients (again, information printed by GAUSSIANxx in the
first 20 columns is discarded).
 8.2.f. if ICON(14) = 2 then
NR lines with the new Gaussian format (5d18.8), containing the
coefficients
 Lines giving the orbital occupancies.
 9a. One line (I3) with the number of occupied orbitals.
 9b. One line (I3,F12.0) for each occupied orbital with its
number and fractional occupancy (number of electrons on the m.o.,
< 2).
3. Description of the output
 The data read. The basis set data will be printed twice.
If renormalization was requested, the second printing will show the
renormalized coefficients. Otherwise the two prints will be
identical.
 The errors (unit matrix  overlap matrix) are given as a
lower triangular matrix.
 M(RHO,x_{o}, n,{e_{i}})
and n_{1},n_{2},n_{3}.
No moment smaller in absolute
value than the threshold specified by
ICON(11) is printed.
 If ICON(2) > 0 then the ID numbers
of the moments sets written on unit
IUMOM
for reference.
4. Files
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)
 NLMN(II,JJ): x,y,z exponent of the IIth moments
 PLMN(II) : value of the IIth moment.
5. Subroutines used
b. Subroutines from PA60 with major modifications.
 NORMLH: finds the harmonic order of each primitive and
normalizes primitive and basis functions;
 GGNTDE: computes the overlap integral between two primitives;
c. New subroutines for polynomials: INIT, IADDRS, MULT, SHIFT,
POLROT, described in detail
Appendix II.
d. New subroutines.
 MOMINT: computes all multipole moments;
 ONEDET: Calculates the coefficients for the contribution of
each basisfunction pair assuming that the input wave function is
a single determinant;
 SHFMOM: performs the transformation of a set of moments
expanded about a center A to a set of moments expanded about a
center B;
 PRILMN: prints the moments computed to the standard output
(unit 6) and, if ICON(2)>0, on unit
IUMOM
with unformatted write (in binary format).
 ROTMOM: transforms the moments with respect to a frame A to
moments with respect to a rotated frame;
 SHRPSP: computes the moment contributions for the overlap
densities in case of the sharp split;
 RDSPL: reads in the split coefficients, S(k,p,q) or the w(c),
the weight for center c in the sharp split;
Possible extensions and/or alterations
 a. The program could be extended for CI wave functions by
replacing the subroutine call to ONEDET (in MOMINT) to a call
preparing different coefficient sets. ICON(22) could be used to
select such new option.
 b. The storage allotted to the array ILMN is that required for
the calculation of split moments, (n_{c})*MAXLMN). In the
event that
moments of a very high order are desired for a single center
calculation only, considerable storage can be saved by
redimensioning ILMN to MAXLMN. (ICON(6) and ICON(7) must be
positive and zero, respectively).
 c. The current code computes all integrals for all pairs of
centers. As the size of the molecules increases, it may become
desirable to introduce a test to omit integrals for any pair of
centers whose distance exceeds an appropriate function of the
Gaussian exponents and contraction coefficients of the primitives.
 d. If higher order primitives are required than the current
maximum (f orbitals, MAX=3), MAX has to be increased with the
concomitant increase of the dimension of the array ETA.
C. MOMTRNSF
1. General description
This program inputs a set of cartesian multipole moments
(calculated, e.g., by the program MOMENTS):
M(RHO,x^{o},n,{e_{i}}) =
INTEGRAL_{D}
PRODUCT^{3}_{k=1}
(x_{k}x_{k}^{o})^n_{k}
RHO(x) dx
where
 {e_{i}}: the set of orthonormal basis
vectors;
 n =
<n_{1},n_{2},n_{3}>, n_{i} a
nonnegative integer;
 RHO(x) : a density defined over the region D;
 x^{o} : the reference center of the
multipole expansion.
and can perform any of the following transformations
 (a) All multicenter expansion can be converted to a single
center expansion as follows. Given the multipole moments of a
distribution RHO^{j} with respect to a reference center
x^{oj}, i<j<N ,
and any preselected one of the reference centers,
x^{ok}, for each j<>k the moment
M(RHO^{j},x^{oj},n,{e_{i}})
can be transformed to M(RHO^{j},x^{oj},
n,{e_{i}}).
Then the results are summed to give
SUM M(RHO^{j},x^{ok},n,
{e_{i}}),
the moment n of the density
SUM_{j=1} RHO^{j}(x)
with respect to x^{ok}.
 (b) The set of moments
M(RHO^{j},x^{oj}, n,{e_{i}})
of RHO(x) with respect to
one given center can be transformed to the set
M(RHO^{j},x^{oj}, n,{e_{i}})
with respect to a second center,
x^{o}.
This can be repeated for a sequence of
x^{o}'.
 (c) The set of moments
M(RHO^{j},x^{oj}, n,{e_{i}})
of RHO( x) in a coordinate
system whose basis vectors are { e_{i}}
can be transformed to moments
in a rotated frame whose basis vectors are
{ e_{i}'}.
In this case
(e_{1},e_{2},e_{3}) =
(e_{1}',e_{2}',e_{3}') R
where R is a
rotation matrix.
2. Detailed description of the input
 Line 1: Identification of the run. (20A4)
 Line 2: ICON(I),I=1,24 (24I3)
Control constants of the run:
 ICON(1)
Highest value of n=n_{1}+n_{2}+n_{3}
to be used in the run; n<l4 (this limit
can be changed by redimensioning the program(s)).
 ICON(2)
 ICON(2) = 0 : Computed moments are written to unit 6 only.
 ICON(2) = 1 : Computed moments are also written on unit 30
in binary form (see MOMENTS for the
format).
 ICON(8)
Identifier number (ID) of the first output
moment set written on unit 30. ID numbers for
subsequent sets are incremented by one each.
 ICON(11) > 0: Real numbers less than 10^(ICON(11)) will be
treated as zero.
 ICON(13)
 ICON(13) = 0 : Read the moments from unit 5 as described in
5 below.
 ICON(13) = 1 : Read the moments from unit 20 in
binary form
as written by the program MOMENTS.
 ICON(14) 1(0): Shift (do not shift) the moments to a
different expansion center.
 ICON(15)
0(1): Atomic unit (Å) units assumed for the input
coordinates.
 ICON(16) : If ICON(14)=1, shift the moments to ICON(16)
different centers, x^{o}', to be read in (see
5 below).

ICON(17)
 = 0: Do not calculate moments in a rotated
coordinate system.
 > 0: Calculate the moments in a rotated coordinate
system.
 = 1: The elements of the rotation matrix R are
given in the input.
 = 2: The elements of the
matrix R are obtained from the
Eulerian angles
read (see Appendices
III and IV.
 ICON(19)
 = 0: Do not contract moments for a sequence of
centers to a single center.
 > 0: The number of sets of moments which are to be
contracted.
 ICON(21)
The index number of x^{ok} in
Section (1a) above (when ICON(19)>0).
 If ICON(19) > 0
ICON(19) lines, each of which contain the name and the
coordinates of one center:
NAME,X,Y,Z (A4,8X,3F12.0)
 ICON(19) (or just 1 if ICON(19)=0) sets consisting of the
following
 (i) NOFMOM, NID (2I5)
 (iia) If ICON(13)=0 then NOFMOM lines containing
M(RHO,x^{oj},n,{e_{i}}),
n_{1}, n_{2}, n_{3} (E15.9,5X,3I2)
 (iib) If ICON(13)>0, then unit 20 is scanned for the moment
set with ID number NID
and the moments are read from there in the
format established by the program MOMENTS.
 If ICON(14) = 1 then
 (i) Coordinates for the original center of expansion:
NAME,X,Y,Z (A4,8X,3F12.0)
 (ii) ICON(16) lines each containing the coordinates for
one new expansion center:
NAME,X,Y,Z (A4,8X,3F12.0)

 (i) If ICON(l7) = 1: the elements of
R rowwise. (3E15.9)
 (ii) If ICON(17) = 2: the Eulerian angles PHI, THETA, PSI in
degrees. (3F14.0)
Note, that if both ICON(14) and ICON(17) are nonzero,
they both will operate on the input moment set.
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
a) From the program MOMENTS:
LAB, SHFMOM, PRILMN, INIT, IADDRS, SHIFT, POLROT
b) EULER: a subroutine to calculate the
R given as input the
Eulerian angles as defined by Goldstein
[9] and described in
III and
IV.
D. CHARDIR
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, Y_{n}(x,y,z)
satisfying
Laplace's equation the multipole moment P^{n} and
characteristic directions (unit vectors)
s_{1},...,s_{n} are
determined to satisfy
Y_{n}(x,y,z)/r^{2n+1} =
[(1)^{n} * P^{n} /n!] *
PROD^{n}_{i=1}
(s_{i}.DEL) (1/r)
_{r=(x,y,z)}
The program calculates the poles and characteristic directions for
a set of Y_{n}'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
 Title of the problem (20A4)
 Control constants (options) (7I5) :
IUNFMI, IUNFMO, NREP, NOFINT, MAXORD, NID, IECHO
 IUNFMI : Unit number of the binary input of the cartesian
moments. If given as zero, formatted input is assumed from unit 5.
 IUNFMO : Unit number for the binary output of the calculated
characteristic directions. If given as zero, only the formatted
output on unit 6 will be produced.
 NREP : The number of random points to generate for the
determination of the poles. Recommended value: 25.
 NOFINT : The number of nonzero polynomial coefficients for the
spherical harmonics of all orders to be read from unit 5
(IUNFMI=0).
 MAXORD : The highest spherical harmonic order for which poles
and characteristic directions are to be calculated.
 NID: Identifier number of the dataset. For unformatted
(binary) input, it has to match the
identifier (ID) number of the
moment set for which the calculation is to be run. The calculated
poles and characteristic directions will inherit this number.
 IECHO: If > 0, items read from unit 5 will be echoed
annotated with the item number (except the moment values that are
echoed automatically).
 Spherical harmonic coefficients
 3a. If IUNFMI=0 then the spherical harmonic coefficients
are read from unit 5:
NOFINT lines: C, L, M, N (E15.9,5X,3I2)
where C is the coefficient of the term
x^{L}y^{M}z^{N} in the surface
spherical harmonic of order L+M+N.
 3b. If IUNFMI>0, the surface spherical harmonic coefficients
on unit IUNFMT as written by the program MOMENTS.
The program will use the moments set whose
ID number coincides with NID read in the
second line.
 Data for a new moment set:
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.
3. The output of the program
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 x^{L}y^{M}z^{N} 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 Y_{n} 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
 ITITL: description
 MXORD: highest order of pole calculated
 MAXINP: highest order of monopole field read from
IUNFMI.
 NID: dataset
identifier number (same as the
identifier of the moment set used)
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
 INIT: See
Appendix II.
Also generates PROD [(2*j1)].
 RMPS: Calculates the directional derivatives in an arbitrary
direction. Adapted from RMP of the program MULTIPOL.
 MAXWLL: Determines the characteristic directions.
 TAYLOR: Evaluates Y_{n} based on the coefficients read
and cartesian characteristic directions.
 RMPOLE: Calculates the value of a multipole moment from the
ratio of the calculated values of TAYLOR and RMPS  the latter
using P^{n}=1.
 PRNTCD: Prints the calculated multipole moments and
characteristic directions.
 LAGUER, ZROOT: Polynomial root finder routine pair from
Numerical Recipes
[10]. The input array element A(I) is the
coefficient of x^{i1}.
 RANDOM: Congruential random number generator with 32bit
integer arithmetic, using Forsythe's constants.
E. MULTIPOL
1. General description
The program computes permanent multipole, and, optionally,
permanent multipoleinduced dipole, induced dipoleinduced 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 (nonvanishing)
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
x^{L}y^{M}z^{N} 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.
 Title of the problem (20A4)
 INDC, MAXINT, ICINP, ICORR, IFLD, IEOUTL, IPBC,
IREST, MAXDAT, LISTEN, NDSKIP, NDATX, IVIR, NCTRR, ISLTI, IREADFLD,
IECHO
(20I3)
 INDC: = 0 : Do not perform induced calculations.
 MAXINT: Interactions of at most MAXINT order will be
computed. If MAXINT = 0, the order highest allowed by the charge
distribution data (subject to MXORDR and RABMAX, see
4 and 5 below).
 ICINP:
 0: General coordinate input (see
7.1
below) only.
 1,2: Atomic coordinates will be read from unit 8.
The current version also assumes that all molecules are the same,
have three atoms and C_{2v} symmetry (e.g., waters, see
8
below).
 3,4: General coordinate input is followed by input
from unit 8.
 ICORR: =0: Do not use correction potential (see 6 below).
 IFLD:
 = 0: Do not compute the electric field components
unless INDC > 0.
 = 1: Compute electric fields.
 = 2: Compute electric fields and field gradients.
 IEOUTL:
 = 0: Print full output (energy contributions by
interaction order, electric field components, dipoles, etc.).
 = 1: Print only the permanent, induced and correction
contributions to the total energy.
 = 2: Don't print the results (only save them on unit 9)
 IPBC:
 = 0: Do not use periodic boundary conditions.
 = 1: Rectangular periodic boundary conditions will
be applied for all interactions based on the position of the first
atom of the molecule. See also
3 and 4 below.
 IREST:
 = 0: New calculation, save results on unit 9.
 = 1: Run is a continuation calculation (i.e. results of previous
calculations were saved on unit 9).
 = 2: New calculation, but don't save the results on unit 9.
 MAXDAT: Number of energy calculations to be performed.
 LISTEN: >0 : Do not print the previously calculated energies
at the end of the run.
 NDSKIP: Skip NDSKIP configurations both on unit 9 and on the
configuration input (to allow reinsertion of data).
 IVIR: >0 : Do not compute virial sums
[11].
 NDATX: When > 0, for the calculation of average energies
and the correlation coefficients only NDATX configurations will be
used.
 NCTRR: >0 : When a configuration is read from unit 8
(ICINP > 0) read NCTRR atoms per molecule, even though less is
required for the calculation (used to skip dummy atoms).
 ISLTI:
 = 1 : After the calculation of the energy of a
complex, the calculation is repeated for all pairs involving the
first molecule (assumed to be the 'solute')
and the contributions are summed. to give the solutesolvent pairwise
additive energy.
 = 2 : After the calculation of the energy of a
complex, the calculation is repeated for all pairs and the energy
contributions are summed. This will give the purely pairwise additive energy
and thus allows the calculation of the cooperative component.
 IREADFLD:: when > 0, the program reads in the
electric fields at the group centers and ingnores the calculated
values
 IECHO:
 > 0: Echo each inputted number annotated with the item's number
(to help debugging the input)
 (MXORDR(I),I=1,10) (15I3)
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 (I1) Å 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).
 RABMAX (F15.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 charge distribution description. (For water, described
by a wave function calculated by Popkie, Kistenmacher and Clementi
[12] and by the polarizability tensor calculated by Liebmann and
Moskowitz
[13], used in
[6] and
[11], a file is available with the
complete input.)
 If ICORR > 0, then the correction potential description is
read:
 6.1. NEXP (I3):r
the number of different exponents in the potential
 6.2. IEXP(I), I=1, NEXP (4I3) :
the exponents in the potential
 6.3. NCOEF (I3) :
the number of coefficients to read in
 6.4. NCOEF times:
COEF, TYPE1, TYPE2, IE (E20.12,3I3)
It will result in a contribution of COEF*R(TYPE1,TYPE2)^(IEXP(IE))
to the energy where TYPE1 and TYPE2 refer to atoms with
distribution type TYPE1 and TYPE2, respectively and R(TYPE1,TYPE2)
is the distance between the interacting centers. Notice that
generally IEXP has to be negative.
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.
 If ICINP=0, 3 or 4: Configuration data or new correction
potential or new polarizabilities. Each set of energy calculations
has to read in the information about the configuration. Optionally,
the correction potential or the polarizabilities can also be
modified. The choice is controlled with the variables NGR and
INPTYP.
NGR, INPTYP, IMATRD, ITWROT (4I3)
 If NGR = 0, STOP the calculation.
 If INPTYP < 3: Geometry input  items 7.1.
 If INPTYP = 3: New correction potential input, as
described under items 6, followed by another input of items 7.
 If INPTYP = 4: New polarizability tensor input, described
under items 7.2. below, followed by an other input of items 7.
 IMATRD, ITWROT: See below for case when NGR>0 and
INPTYP<2.
 7.1. (INPTYP < 3) Geometry specification input. The meaning
of the variables read in item 7 in this special case is as follows:
 If NGR < 0, proceed to the energy calculation (i.e., no more
input is needed).
 If NGR > 0, it is the number of groups (molecules to be
read in.
 INPTYP: 012: controls the input items. For INPTYP<2, part
ofthe input can be skipped (useful when a second configuration is
just a modification of the first).
 IMATRD: 01: Rotation is specified by
Eulerian angles (see
7.1.2.1. below) or a rotation matrix.
 ITWROT: 01: If ITWROT=0 the first group (molecule) will not
be rotated, but will be assumed to have its principal axis system
parallel to the global system.
All items under 7.1. are repeated NGR times:
 7.1.1. Molecular geometry description. For the first
configuration, INPTYP has to be 2.
Each molecule is described by specifying the atoms in a local
frame (frequently the principal axis frame). This should be
the same frame the characteristic directions are calculated and
described.
The full configuration is specified by giving the position
and orientation of the local frame of each molecule in a common
frame, usually referred to as the laboratory or global frame.
 7.1.1.1. If INPTYP > 0: R1,R2,R3 (3E15.9) :
the coordinates of the center of the local coordinate system for
this molecule (in a.u.) in the laboratory frame.
 7.1.1.2. If INPTYP=2 :.NCTR (I3) : the
number of centers in the molecule.
 7.1.1.3. If INPTYP=2 then NCTR numbers: Type1, ..., TypeNCTR
(24I3)
The type codes of the different centers in the
molecule. The type codes are the index of the charge distribution
in the chargedistribution list.
 7.1.1.4. NCTR times: RC1, RC2, RC3 (3E15.9) :
the coordinates of each center in the local system (in a.u.)
 7.1.2. Rotation description (see Appendices
III and IV
for the conventions) and quantummechanical energy of the molecule.
If ITWOROT=0 then omit it for the first molecule.
This item describes the orientation of each molecule's
local frame in the laboratory system.
 7.1.2.1. If IMATRD = 0: PHI, THETA, PSI, Energy
(3F14.0,E15.9)
PHI, THETA, PSI: the three
Eulerian angles
describing the orientation of the molecule.
 7.1.2.2. If IMATRD = 1:
Rotation matrix by rows, Energy
(3E15.9,/,3E15.9,/,3E15.9,E15.9)
 7.2 If INPTYP=4 then new polarizability tensor
ALPHA input (to be followed by an other input of item 7):
 7.2.1. NALP (I3) : Number of polarizability sets to be read
 7.2.2. NALP times: ALPHA_{11}, ALPHA_{22},
ALPHA_{33}, TYPE, (3E15.9,I3)
The polarizabilities belonging to the distribution type
TYPE in the subsequent calculations will be ALPHA_{11},
ALPHA_{22}, ALPHA_{33}.
 Configuration specification:
If ICINP > 0 then configuration specification by atomic
coordinates is read from unit 8.
3. The output of the program
The output printed on unit 6 is as follows:
 Program name, title, date.
 Charge distribution data read in. The
characteristic directions are printed in the local (principal axis)
coordinate
system, that is, after the prerotation (if ICHRRT(I)>0)
 Geometry data read in (coordinates, rotation matrices).
 For IEOUTL=0: If INDC=0 or IFLD>0 then the calculated
components of electric fields and, if IFLD>1, field gradients.
The components will be given with respect to the local frame at the
polarizable center (i.e., the center where the charge distribution
type involved nonzero polarizabilities). The field gradients are
grouped into triples, each representing the components of the
gradient of one field component.
 If IEOUTL=0 and INDC=1 (induced calculation) then for each
molecule:
 5.a The components and the magnitude of the electric fields
in the global frame at the polarizable centers;
 5.b The components of the induced dipoles in both the local
and the global fields and their magnitudes;
 5.c The components of the total molecular dipole vectors
(permanent plus induced) in both the local and the global fields
and their magnitudes;
 5.d The angles between the permanent and induced dipole
vectors;
 5.e The angles between the permanent dipole vectors and the
electric field;
 If IEOUTL=0 then the contribution of each spherical
harmonic order to the permanent moment energy.
 If IEOUTL > 0 then various energies and, if IVIR > 0,
virial sums are printed:
 7.a the cumulative contributions of each harmonic order to
the total permanent multipole energy and to the virial sum.
 7.b if distance cutoff was used (see
item 4 of the input), the number
of molecule pairs and their average distance whose interactions
were computed;
 7.c the average induced and total molecular dipoles, the
range and the average electric fields at the polarizable centers;
 7.d the permanent multipole energy and virial sum;
 7.e if INDC=1 then the induced dipole contribution to the
energy and the virial sum;
 7.f if ICORR>0 then the correction potential energy;
 7.g the total energy.
REPEAT from 3 for the next configuration.
 At the end of the run, the energies and virial sums
calculated for each configuration will be listed in a tabular form,
their average and standard deviation will be calculated as well as
the correlation coefficient between the calculated total energy and
the original pairwise energy (if available).
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,VPWA
where 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:
 ERSLT(1,I), VRSLT(1,I) : calculated total energy and virial
sums of configuration I
 ERSLT(2,I), VRSLT(2,I): Correction potential contributions of
configuration I
 ERSLT(3,I), VRSLT(3,I): Induced contributions of configuration
I.
 ERSLT(4,I), VRSLT(4,I): Permanent moment contributions of
configuration I.
 ERSLT(4+J,I) VRSLT(4+J,I): When ISLTI > 0, the sum of the
pairwise calculated total, correction, induced and permanent moment
energies for J=1,2,3,4, respectively.
 EPWA(I),VPWA(I),NMC(I): If the input configurations were
obtained from a Monte Carlo simulation (read from unit 8), the
pairwise energy, pairwise virial sum and the Monte Carlo
configuration number of configuration I, respectively.
5. Subroutines used
6. The dimensioning of the program
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.
 @@@ : Highest order interactions allowed (MXORD) + 2
 @@3 : 3*MXORD = 3*(@@@2)
 @3S : Maximum number of characteristic direction components,
[MXORD*(MXORD+1)/2]*3
 @@M : Maximum number of moments,
(MXORD^{3}+6*MXORD^{2}+11*MXORD+6)/6
 @@T : Maximum number of monopole fields.
 ### : Maximum number of molecules (groups).
 $$$ : Maximum number of atoms per molecule (centers per group).
 &&&: Maximum number of atoms (centers) and
molecules (groups).
 @@S : Maximum number of different charge distributions.
 @SS : Total number of characteristic direction components,
@@T * @3S .
 @@E : Maximum number of correction potential exponents.
 ##D : Maximum number of energies and virial sums saved.
A cshell
script and a VAX/VMS command file are given
provided with the distribution to perform the dimensioning,
as described in Sec. I.
F. CRYSCON
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 chargedistributionindependent
sets of quantities, called crystal constants
K^{N}_{T}(NU,X^{c}X^{j}).
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
X^{c} and
X^{j} are in the
same unit cell. The crystal constant have to be calculated only for
crystallographically nonequivalent sitelattice pairs.
As input, it requires:
 (a) the components of basis vectors for
the crystal with respect to an orthogonal Cartesian frame;
 (b) the
size of the finite lattice over which the required summations are
to be performed;
 (c) for each translational lattice the
X^{c}X^{j}
vector;
 (d) the arbitrary Ewald parameter EPS. Calculations can be
performed with more than one EPS value for checking purposes. The
calculation stops when the list of translational lattices is
exhausted.
The program prints the computed crystal constants as well as
the two lattice sums (K_{1} and K_{2} 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
 Title of the first translation lattice. (20A4)
 MX0, MAXRRT, IBINCC, NID (12I3)
 MXO: highest spherical harmonic order crystal constant to
be computed;
 MAXPRT: Highest spherical harmonic order crystal constant
to be printed;
 IBINCC 1 (0) : Write (do not write) the crystal constants on
unit 20 in binary form;
 NID:
Identifier number
of the first crystal constant set generated
(subsequent sets will have their NID incremented by one each).

 3.a. L1MIN, L2MIN, L3MIN; (12I3)
 3.b. L1MAX, L2MAX, L3MAX. (12I3)
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
{i_{k}} represent a site
<i_{1}a_{1},
i_{2}a_{2},
i_{3}a_{3}> where
a_{k} are the crystal
basis vectors (see 5 below).
The central cell is assumed to have all zero indices.
 ((ACOMP(I,J), J=1,3), 1=1,3) (3F20.0)
ACOMP(I,J): the Jth component of the Ith unit cell basis
vector. These can either be scaled or unscaled
[7]
 The following sequence of lines for each unit cell site
NSITE and translation lattice NLATT.
 a). The first line is
RT(1), RT(2), RT(3), NSITE, NLATT. (3F20.0,2I3)
If NSITE=0, the calculation will end.
Here RT =
X^{NSITE}X^{NLATT} .
 b). EPS: EPS (3F20.0,I3) : the Ewald parameter.
The first of a possible sequence of EPS values must be
positive. Then the crystal constant will be computed
for the sitelattice pair described in line (5a).
If this EPS is positive the crystal constant calculation for (5a)
will be repeated. If this EPS is negative, read in line (5c).
 c). Next translation lattice calculation's title. (20A4)
Then proceed to (5a) for the next translation lattice.
3. Output of the program
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:
 (i) the reciprocal of the unit cell volume
AINV=a_{1}.(a_{2}x
a_{3});
 (ii) the relative errors for the crystal constants over the
direct and the reciprocal lattices. Sums are computed for distance
classes until the contribution lies below a threshold. The relative
error for the each lattice is the largest of the ratios
(contribution from the last class)/(value of the constant) for all
NU.
4. Files
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
 INITCRC: Initializes constants.
 GENRQVC: Generates the auxiliary QLT vectors and their
magnitudes.
 SORTIT, SORTRC: Together they sort the lattice distances into
classes of identical magnitude.
 CRSTC1: Calculates the crystal constant contributions KNT1.
 CRSTC2: Calculates the crystal constant contributions KNT2.
 SELFCOR: Calculates and applies the correction for the
selftranslation lattice,
i.e., for the case when the site of index NSITE
belongs to the translation lattice of index NLATT. This uses a
correction to an equation of
[7], which is derived in
Appendix II.
 SAVECT: Writes the calculated crystal constants on unit 6 and,
optionally, to unit 20.
G. PROGRAM CRYSTAL
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
noncooperative 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.
 a. Options.
 b. Charge distributions. Each type of charge distribution
(molecule) is specified by the list of characteristic directions
and the corresponding multipole moments (as computed, e.g., by the
program CHARDIR).
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.
 c. Crystal constants. The crystal constants for each
translation lattice can be specified by their numerical values or
by information which tells the program how to obtain them from the
crystal constants of another translational lattice that have been
read in already, as described under item 14 in the input list. The
program automatically recognizes the trivial identities for all
<i,j>
K(NU,X_{i}X_{i})=
K(NU,X_{j} X_{j})
K(NU,X_{i}X_{j})=
(1)^(NU_{1}+NU_{2}+NU_{3})
K(NU,X_{j}X_{i})
Furthermore, identities between crystal constants are used to
generate identities between the sitelattice permanent multipole
energies, U_{p}(X_{c},T_{j}),
and between the permanent multipole
electric fields at sites
X_{c}, by lattice T_{j},
E_{p}(X_{c},T_{j}), 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
X^{i}X^{j} =
X^{i'}X^{j'} ,
then
E_{p}(X_{i},T_{j}) =
E_{p}(X_{i'},T_{j'}) and
U_{p}(X_{i},T_{j}) =
U_{p}(X_{i'},T_{j'}) .
If for two pairs <i,j> and <i',j'>
<C(i),C(j)> = <C(j'),C(i')> and
X^{i}X^{j} =
(X^{i'}X^{j'})
then
U_{p}(X_{i},T_{j}) =
U_{p}(X_{i'},T_{j'}) .
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
E_{p} and the U_{p}. 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.
 d. There is an option, which permits the user to bypass the
input steps bc and read in the
E_{p} and the U_{p} values computed
previously. Other options permit the printing of the
E_{p} and the U_{p}
values and/or their saving on a disk file.
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 noncooperative calculation can also be
performed.
2. Detailed description of the input
 Title (20A4)
 NUNIT, MXOEN, INDC, IUNFCD, IUNFCC, IUNFTB,IREPRT, MXPRT,
ITABPR, ITABUF, ITABRD, ICCFIL,ISTNPR, IFLDPR, IDIPPR, ITDIPP,
ISMPIN (24I3)
 NUNIT: Number of sites in the unit cells
 MXOEN: Highest order interaction to be computed;
 INDC =1(0): Compute (do not compute) induced moment energies;
 IUNFCD=0 : Characteristic directions and poles are read in
decimal format [cf.
(9c and 9d)] from unit 5;
 IUNFCD=1 : Characteristic directions and poles are read in
binary format from unit 10;
 IUNFCC=0 : Crystal constants are read in decimal format
[cf 14. below] from unit 5;
 IUNFCC=1 : Crystal constants are read in
binary format from unit 20;
 IUNFTB=1(0): use binary
(alphanumeric) [cf. (16) below] format
when permanent multipole energy and field tables are
read (ITABRD=1) or written(ITABPR=1),
 IREPRT=1(0): Print (do not print) the characteristic
directions, rotation matrices and, if they are included as input,
tables of permanent multipoleelectric fields and energies;
 MXPRT : Highest order crystal constant that is to be
printed;
 ITABPR=1(0): Print (do not print) permanent multipole electric
fields and energies;
 ITABUF=1(0): Write (do not write) the permanent multipole
electric fields and energies un unit 30;
 ITABRD=1(0): Read (compute) the permanent multipole electric
fields and energies;
 ICCFIL=0(1): Save the crystal constants and duplicate
information on unit 4, (Read the crystal constants and duplicate
information from unit 4)
 ISTNPR=1(0): Print (do not print) the site energies;
 IFLDPR=1(0): Print (do not print) the components of the
electric field at each site;
 IDIPPR=1(0): Print (do not print) the components of theinduced
dipole vector (global frame);
 ITDIPP=1(0): Print (do not print) the total dipole vector of
the unit cell;
 ISMPIN=1(0): Compute (do not compute) the induced dipole
vectors noncooperatively as well as cooperatively.
 R0 (F20.0) : Scaling factor for the distance.
Each crystal constant
K(NU,X_{c}X_{j})
is divided by
R0^(NU_{1}+NU_{2}+NU_{3}+1). This
adjusts for any scaling of the input crystal constants.
 IRTYP(I), I=1,NUNIT (24I3) :
Index numbers for the different types of rotation matrices by site.
 ICTYP(I), I=1,NUNIT (24I3) :
Index which defines the type of charge distribution for site I.
 NRTYP (I3) : Number of distinct rotation matrices.
 (((ROT (ITYPE,I,J), J=1,3), I=1,3), ITYPE=1,NRTYP) (3F20.0) :
The rotation matrices by type, listed rowwise.
 NTYP (I3) :
Total number of types of (i.e., different) charge distributions.
 For each ITYP 1 < ITYP < NTYP:
 (a). MXO, NID, IUCD (I2,I8,I10)
 MXO: Highest spherical harmonic order included in theinput,
 NID: Characteristic direction set
identifier (used when IUNFCD>0),
 IUCD: When neither IUCD nor IUNFCD are zero, IUNFCD is replaced
by IUCD.
 (b). ALPHA(I), 1=1,3 (3E15.9) :
The diagonal component of the dipole
polarizability tensor in the principal polarizability axis frame,
in (a.u.)^{3};
If IUNFCD=0 then the poles and characteristic directions
are read in 9c and 9d:
 (c); POLE(J), J = 1, [MXO+1] (3E15.9) :
multipole moment of order
J1;
 (d). <s(1), s(2), s(3) > (3E15.9) :
Components of the characteristic directions. The MX0(MX0+1)/2
s must be read in sets of increasing multipole order
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).
 NDUPEN (I3) : Number of elements in the list (item 11)
which specify duplicate sitelattice interaction energies.
 NDUPEN lines.
ISD, ILD, ISC, ILC (4I3) : indices such that
 U_{p}(ISB, ILD) = U_{p}(ISC, ILC) where
the first variables, ISC and ISD, are indices for sites and the
second, ILC, ILD are indices for lattices.
 NDUPFLD (I3) :
Number of elements in the list (item 13) which
specifies duplicate electric field components.
 NDUPFLD lines.
ISD, ILD, K, ISC, ILC, ISG (6I3)
 ISG is a constant for which the Kth components of the permanent
multipole fields,
 E_{p} defined by the translation lattice
ILD(ILC) at the site ISD(ISC) satisfy the equation
E_{p}(ISD,ILD,K) = ISG*
E_{p}(ISC,ILC,K).
 If (ICCFILE .EQ. 0 .AND. IUNFCC .EQ. 0), then read the
following data (in the subroutine READKNT).
 (a) NCRCT (I3) :
Number of <ISITE,ILATT> pairs for which crystal
constants are to be read;
 (b) For each of the NCRCT pairs for which crystal constants are
to be read in, read the following line sequence:
If (ICCFILE .EQ. 0 .AND. IUNFCC .EQ. 1), then read first from
Unit 5:
 a) NCRCT (I3) :
Number of <ISITE,ILATT> pairs for which crystal
constants are to be read [cf.(14b)];
 (b) For each of the NCRCT pairs
read first from Unit 5 NID, (I3) the
identifier number of the crystal constant set.
After that, the program reads from unit 20 the crystal
constants in
binary format as written by the program
CRSYCON.
If ITABRD > 0 AND IUNFTB = 0 :
 a) Title of the U_{p} and
E_{p} tables; (20A4)
 b) The U_{p}(ISITE,ILATT)
values:
(SLENRGY(ISITE,ILATT),ILATT=1 NUNIT) (3E24.12) :
The permanent
multipole interaction energy between the unit cell site ISITE and
the translation lattice ILATT.
 c) If INDC = 1 then the
components of the E_{p}(ISITE,ILATT):
((SLFIELD (ISITE,ILATT,K),K=1,3),ILATT=1,NUNIT) (3E24.12) :
the Kth component of the electric field defined at the unit cell site
ISITE by the permanent multipoles at the sites of the translation
lattice ILATT.
If ITABRD > 0 AND IUNFTB = 1 :
The energy and (if INDC=1) the field table in
binary form
read from unit 30 as written by an earlier run.
3. The output of the program
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,
sitelattice contributions.
4. File
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
 INITCRC: Initializes constants.
 GENSIG: Generates the SIGMA vector
[7] from a set of
characteristic directions.
 INDUCED: Calculates the induced dipoles and induced energies.
 READSP: Reads in the characteristic directions.
 DISKTR: reads or writes a crystal constant set to or from unit
4.
 RDKNT: Reads in the crystal constants.
 LUDCMP, LUBKSB: LU decomposition solver of a linear equation
system
[10]
H. CRYSPOT
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:
SUM_{{n}} SUM_{{<i,j>}}
SUM_{{<p,q>}}
X_{<i,p>} 
X_{<j,q>}^{n}
 i,j: indices for sites in the crystal,
 p,q: sites within the charge distributions (molecules) at
lattice sites i and j, respectively.
The crystal geometry is specified by the three orthogonal
(Cartesian) components of the crystal axes, the coordinates of the
sites in the unit cell and a scaling factor.
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
inversedistance 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
S^{n}_{<A,B>} =
SUM_{{n}} SUM_{{<i,j>}}
SUM_{{<p,q>}}
X_{<i,p>} 
X_{<j,q>}^{n}
where type(p)=A and type(q)=B
A set of coefficients, d^{n}_{<A,B>}, may
also be given. In this
case the lattice energy contribution
SUM_{{n}} SUM_{{<A,B>}}
d^{n}_{<A,B>}
S^{n}_{<A,B>}
will also be computed.
2. Detailed description of the input
 Program Title. (20A4)
 NUNIT (24I3) :
Number of sites in the unit cell.
 The summation of the powers of distance is performed over
a finite lattice where the upper and lower bounds of the indices
are given by LiNIN, LiMAX, i=1,2,3, respectively. Each index set
{i_{k}} represent a site
<i_{1}a_{1},
i_{2}a_{2},
i_{3}a_{3}>
where a_{k} are the crystal
basis vectors (see 5 below).
 a) L1MIN, L2MIN, L3MIN (24I3);
 b) L1MAX, L2MAX, L3MAX (24I3).
 R0: A scaling factor for distance (3F20.0):
The basis vector components in item 5 will be multiplied by R0.
 ((ACOMP(I,J), J=1,3), I=1,3) (3F20.0) :
ACOMP(I,J); The Jth component of the Ith unit cell basis
vector.
 IMOLTYP(I),I=1,NUNIT (24I3) :
The type of the molecule at site I.
 (RT(I,J),J=1,3), I=1, NUNIT):
RT(I,J): The Jth component of the position of the Ith
unit cell site.
 IRTYP(I), I=1, NUNIT (24I3) :
The index which identifies the rotation
matrix for unit cell set I.
 NRTYP (I3) :
The number of distinct rotation matrices.
 (((ROT(ITYPE,I,J), J=1,3) I=1,3), ITYPE=1,NRTYP) (3E15.9) :
ROT(ITYPE,I,J): The element in the Ith row and Jth column
of the rotation matrix whose index is ITYPE.
 NGR: the number of different types of groups (molecules)
 IG=1,NGR
 (a) NCTR (I3) : The number of sites in the molecule of type IG
 (b) CTYP(I) I=1,NCTR (24I3) :
The index specifying the type of site I
in group IG ;
 (c) ((RC(IG,I,J), J=1,3), I=1,NCTR) (3E15.9) :
RC(IG,I,J): The Jth local coordinate of the Ith site of molecule
IG.
 NEXP (I3) :
The number of exponents for the potential.
(IEXP(IE), IE=1, NEXP) (12I3) :
The negative of the exponent to be used.
 NCOEF (I3) :
If NCOEF = 0, only the sums of the inverse powers will
be computed. Otherwise, NCOEF is the number of non zero
coefficients in the potential (d^{n}).
 C, IT1, IT2, IE (E20.12,3I3) :
 C=COEF(IT2,IT1,IE) where
 IE is the index of the exponent.
1 < IT1 < NTYP and
1 < IT2 < NTYP are the indices which
specify the types of sites as identified for a particular group in
(12b).
3. The output of the program
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
 INITCRC: Initializes constants.
 CORRPOT: Reads in the parameters of the correction potential.
I. INSTALLATION
The programs are written in Fortran77.
A cshell 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 nonUNIX systems the following steps are needed:
 Eliminate all CALL SYSTEM statements from the source code
(they are nonessential).
 Lines required by IBM (MVS or VM) are given prefixed with
CIBM and by VAX/VMS with CVAX
 these comments have to be blanked out, as
appropriate.
 The programs have to be compiled one by one with the
requisite Fortran compilation commands.
The program MULTIPOL requires a
preprocessor (see Sec. E.6)
to set the dimensions as needed. For Unix systems
the cshell 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:
 Maximum harmonic order: 10,
 Maximum number of molecules: 125,
 Maximum number of atomtypes: 3,
 Maximum number of configurations: 500,
 Maximum number of correction potential exponents: 4.
For VAX/VMS systems
the command file compile.com is provided.
The current version uses single precision arithmetic in most
places. This limited the precision of the calculated moments to 45
significant figures for a water molecule in a calculation with
32bit arithmetic.
The most sensitive part of the package is the calculation of the
socalled 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).
J. APPENDICES
Appendix I. The Correction for the Selflattice Sum of a
Directional Derivative
Reference
[7] derived the Eqs. (A.11,12) for the correction
and the selflattice sums of directional derivatives. However, it
did not give the decomposition required by the algorithm for the
study of a larger number of orientations:
R^{N}_{THETAo,2}({s_{k}},EPS)
= SUM_{{NU}}
R^{N}_{THETAo,2}(NU,EPS)
SIGMA({s_{k}},NU). (A1.1)
Since in Eq. (A.11) j = q = N/2, Eqs. [(38), (28d), (36), (37b,c),
(23b)] imply:
R^{N}_{THETAo,2}(NU,EPS) =
{[(1)^{(N/2)+1} 2 EPS^{N+1}] /
[PI^{1/2}(N+1)]}
PROD^{3}_{i=1} [NU_{i}!/(NU_{i}/2)!]. (A1.2)
Finally, in Eq. (40) for the selflattice sum of a directional
derivative,
R^{N}_{THETAo,2}(NU,EPS) should replace
R^{N}_{THETAo,2}({s_{k}},EPS).
In addition, the equation references should read:
K^{N}_{To,1}, (25b,b.14);
K^{N}_{To,2}, (39a,a.14).
Appendix II: Subroutines for handling polynomials of
three variables
The polynomial of degree n is assumed to be of the form
SUM_{L+M+N < n}
c_{LMN}x^{L}y^{M}z^{N} . (A2.1)
The coefficients c_{LMN} 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
SUM^{n}_{i=0} i^{2} =
(2n^{3}+3n^{2}+n)/6 . (A2.2)
The position of c_{LMN} is
n'*(n'*(n'+3)+2)/6 + N*(2*n'N+3)/2+M+1 (A2.3)
where
n'=L+M+N . (A2.4)
 INIT sets up the following constant arrays:
 FACT(I): (I1) factorial;
 NOFMEM(I): SUM^{I1}_{j=0}
(j+1)*(j+2)/2 , the number of terms in a polynomial of degree I1.
 EXPOS(1,I), EXPOS(2,I), EXPOS(3,I): L, M, N, for the
coefficient c_{LMN} at the Ith place in the linear array;
 BINOM(I,J): The binomial coefficient
B_{I,J}=I!/[J!*(IJ)!].
 IADDRS(L,M,N) returns the address of the polynomial
coefficient c_{LMN}.
 MULT(FA,NA,FB,NB,PROD): calculates the coefficients of a
product of two polynomials. The arrays FA and FB hold the
coefficients of the two input polynomials of degree NA and NB. The
product's coefficients will be in PROD.
 SHIFT(SOURCE,DEST,NDEG,X0,Y0,Z0): computes the coefficients
for a polynomial when the center of the coordinate system is
shifted, i.e.
SUM c_{LMN} (xx_{o})^{L}
(yy_{o})^{M} (zz_{o})^{N} =
SUM c'_{LMN} x^{L} y^{M} z^{N} .
(A2.5)
 SOURCE contains the coefficient c_{LMN} and DEST
contains the
coefficients c'_{LMN}. The degree of the polynomial is
NDEG.
Note that
(xx_{o})^{L} (yy_{o})^{M}
(zz_{o})^{N} =
SUM_{l<L,m<M,n<N}
B_{L,l}B_{M,m}B_{N,n}
(x_{o})^{Ll} (y_{o})^{Mm}
(z_{o})^{Nn} x^{l} y^{m}
z^{n} (A2.6)
 POLROT(R,SOURCE,DEST,NDEG): computes the coefficients for a
polynomial after rotation of the coordinate system, i.e., after the
transformation
r' = R r (A2.7)
where
 r and r' are column vectors and
 R is the rotation matrix (see
Appendix III.
for the conventions);
 SOURCE holds the coefficients
with respect to r;
 DEST with respect to r',
 NDEG is the degree of the polynomial and
 R is a 3x3 matrix representing the rotation (see
Appendix III.
for the conventions).
Note that the LMN term of the polynomial in
r' becomes
( SUM^{3}_{i=1}
R_{1i}r_{1})^{L}
( SUM^{3}_{i=1}
R_{2i}r_{2})^{M}
( SUM^{3}_{i=1}
R_{3i}r_{3})^{N} (A2.8)
and
( SUM^{3}_{i=1}
R_{1i}r_{i})^{L} =
SUM_{}
[L!/(l!m!n!)] R_{11}^{l}
R_{12}^{m}
R_{13}^{n}
r_{1}^{l} r_{2}^{m}
r_{3}^{n} (A2.9)
Appendix III. The Convention for the Matrices for the
Transformation of Basis Vectors and the Components of a Vector
A matrix with one row will be denoted by () and a matrix with one
column by []. Let e_{j}(i) and e
_{j}(f) denote the jth 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:
(e_{1}(f), e_{2}(f),
e_{3}(f)) = (e_{1}(i),
e_{2}(i), e_{3}(i)) B
. (A3.1)
Let x_{j}(i) and x_{j}(f) denote the jth
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:
x_{1}(i) x_{1}(f)
[x_{2}(i)] = B [x_{2}(f)] (A3.2)
x_{3}(i) x_{3}(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:
x_{1}(f) x_{1}(i)
[x_{2}(f)] = B^{T}[x_{2}(i)] (A3.3)
x_{3}(f) x_{3}(i).
This assumes that the frames are orthogonal so that the transpose,
B^{T}, is the inverse. Thus, the transpose,
B^{T}, 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:
(e_{m}(f).e_{k}(i) =
B^{km} = (B^{T})^{mk}. (A3.4)
Thus, the elements of the kth column of
B^{T} are the components of
the initial frame e_{k}(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 noncoplanar case. The two
following initial choices are made here.
 (i) The coordinate plane is the
(e_{1}(i),e_{2}(2)) plane.
 (ii) First rotate e_{1}(i) about the
e_{3}(i) axis by a positive angle PHI to the
line of intersection of the initial and final
(e_{1}(i),e_{2}(2)) planes
(their line of nodes). This gives the first set of intermediate
basis vectors, {e_{j}'}.
Let:
B_{1} =
the rotation matrix which transforms {e_{j}(i)}
into {e_{j}'}. (A4.1)
Since both e_{3}'=e_{3}(i)
and e_{3}(f) are perpendicular to
e_{1}', rotate
e_{3}(i) about the
e_{1}' axis by a positive angle THETA to
e_{3}(f). This gives
the second set of intermediate basis vectors,
{e_{j}"}, such that
e_{3}"=e_{3}(f). Let:
B_{2} =
the rotation matrix which transforms
{e_{j}'} into
{e_{j}"}. (A4.2)
Finally, rotate
e_{2}" = e_{1}' about the
e_{3}(f) axis by a positive angle PSI to
e_{1}(f). Let:
B_{3} =
the rotation matrix which transforms
{e_{j}"} into
{e_{j}(f)}. (A4.3)
Each axis system is assumed to be righthanded and orthogonal
and each rotation is counterclockwise by a positive angle. (If mnp
is any cyclic permutation of 123, the analytical conditions are:
 (i) for a righthanded system,
e_{p} =
e_{m} x e_{n} ;
 (ii) for a counterclockwise rotation about a basis
vector e_{m}(i), e_{n}(f) has
positive components with respect to both
e_{n}(i) and
e_{p}(i), if the angle is < 90^{o}.
The forms of the three rotation matrices follow
immediately.
First, the elements of a matrix C for a rotation about
e_{3}(i) by a positive angle ALPHA are:
 cos (ALPHA) 
s( (ALPHA)  0  
C =  s( (ALPHA) 
cos (ALPHA)  0  (A4.4) 
 0  0  1  
Both B_{1} and B_{3} have the form
C.
Second, the elements of a matrix D for a rotation about
e_{1}(i) by a positive angle THETA are:
 1  0  0  
D =  0  cos (ALPHA) 
s( (ALPHA)  (A4.5) 
 0  s( (ALPHA) 
cos (ALPHA)  
Thus B_{2} = D. Finally, since the matrix
B of equation (III.1) is a right matrix multiplier, B
= B_{1}B_{2}B_{3}.
According to equation (III.3),
B^{T} =
B_{1}^{T} B_{2}^{T}
B_{3}^{T} is the matrix for the
transformation of the components
of any vector x with respect to the initial frame into
those with respect to the final frame.
B^{T}, which the program uses to transform
the components of the characteristic directions from the initial
(local) frame into those for the final (global) frame has the
elements:
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 .
K. REFERENCES
 Description of the Maxwell formalism:
E.W. Hobson, "The Theory of Spherical and Ellipsoidal
Harmonics", Cambridge Univ. Press, (1931).
 Algoritm for calculating the characteristic directions:
M. Mezei and E.S. Campbell, J. Comp. Phys,
20, 110 (1976).
Manuscript (PDF)
Manuscript (PS)
 Methods of partitioning the electon density:
M. Mezei and E.S. Campbell, Theoret. Chim. Acta,
43, 227 (1977).
 Recursion formulae for calculating the directional
derivatives:
E.S. Campbell and M. Mezei, J. Comp. Phys.,
21, 114 (1976).
Manuscript (PDF)
Manuscript (PS)
 System of equations to solve to obtain induced dipoles:
E.S. Campbell, Helv. Phys. Acta,
40, 367 (1967)).
Manuscript (PDF)
Manuscript (PS)
 Cooperative water model based on the multipole expansion
using dipole polarizibility:
E.S. Campbell and M. Mezei, J. Chem. Phys.,
67, 2338 (1977).
Manuscript (PDF)
 Extension of the Ewald summation to multipole lattices:
E.S. Campbell, J. Phys.Chem. Solids,
24, 197 (1963).
Manuscript (PDF)
Manuscript (PS)
E.S. Campbell, J. Phys.Chem. Solids,
26, 1395 (1965).
Manuscript (PDF)
Manuscript (PS)
 Brief description of the Maxwell package:
M. Mezei and E.S. Campbell, J. Comp. Phys.,
29, 297 (1978).
 Herbert Goldstein, "Classical Mechanics", AddisonWesley
Publishing Company, Reading, Mass., First edition, 1950, PP. 1079;
second edition, 1980, pp 1457..
 W.H. Press, B.P. Flannery, S.A. Teukolsky & W.T.
Vettering, 'Numerical Recipes', Cambridge University Press (1986).
 Effective cooperative water model:
M. Mezei, J. Phys. Chem.,
95, 7042 (1991).
 H. Popkie, H. Kistenmacher, and E. Clementi, J. Chem.
Phys., 59, 1325 (1973).
 S.P. Liebmann, and J. W. Moskowitz, J. Chem. Phys.,
54, 3622 (1971).
Return to contents page