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

E-mail 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

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:

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 (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.

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 (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

  1. Title of the problem    (20A4)
  2. Control constants ICON(I)    (24I3)
  3. Lines specifying the center coordinates:
  4. If ICON(3)=1: Split expansion information
  5. If ICON(5)=1: lines defining hybrid orbital coefficients (adm). This allows the definition of primitives where the Gaussian is multiplied by a multiterm polynomial in x,y, and z.
  6. Lines giving the basis function specifications. Both POLYATOM and GAUSSIAN-xx 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.
  7. 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.)
  8. Molecular orbital coefficient matrix.
  9. Lines giving the orbital occupancies.
3. Description of the output
  1. 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.
  2. The errors (unit matrix - overlap matrix) are given as a lower triangular matrix.
  3. M(RHO,xo, n,{ei}) and n1,n2,n3. No moment smaller in absolute value than the threshold specified by ICON(11) is printed.
  4. 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)

5. Subroutines used

  • b. Subroutines from PA60 with major modifications.
  • c. New subroutines for polynomials: INIT, IADDRS, MULT, SHIFT, POLROT, described in detail Appendix II.

  • d. New subroutines.
  • Possible extensions and/or alterations

    C. MOMTRNSF

    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

    and can perform any of the following transformations

    2. Detailed description of the input

    1. Line 1: Identification of the run.    (20A4)
    2. Line 2: ICON(I),I=1,24    (24I3) Control constants of the run:
      • ICON(1) Highest value of n=n1+n2+n3 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, xo', 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 xok in Section (1a) above (when ICON(19)>0).
    3. 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)

    4. 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,xoj,n,{ei}), n1, n2, n3    (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.
    5. 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 row-wise.    (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, 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

    1. Title of the problem    (20A4)
    2. 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).
    3. 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 xLyMzN 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.
    4. 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 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

    E. MULTIPOL

    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.

    1. Title of the problem    (20A4)

    2. 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 C2v 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 solute-solvent 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)

    3. (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 (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).

    4. 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.

    5. 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.)
      • 5.1. NT    (4I3)
          NT: The number of different types of charge distributions.
      • 5.2. For each charge distribution type I (i.e., NT times):
        • 5.2.1. MAXS(I),MAXORD(I),NID,ICHRRT(I),ICDBIN,ICRDRT(i)    (6I3)
          • MAXS(I): the highest order of the multipoles and characteristic direction for the I-th charge distribution type to be read in. Ignored when IDBIN > 0 .
          • MAXORD(I): the highest harmonic order to use for the permanent multipole energy calculations involving this charge distribution type. This can be more or less than MAXS(I) since the order of an interaction is the sum of the orders of the twointeracting multipoles.
          • NID: Identifier number of the characteristic direction set on unit ICDBIN.
          • ICHRRT(I) > 0 : The characteristic directions will be transformed into a rotated frame before use (referred to as 'prerotation'). See Appendices III and IV for the conventions.
          • ICDBIN > 0 : The multipole moments and characteristic directions are read from unit ICDBIN in binary form (as written by CHARDIR).
          • ICRDRT(I) > 0 : The 'prerotation' will also applied to the coordinates. Ignored when ICHRRT(I)=0.
        • 5.2.2. If ICHRRT(I)>0 : prerotation information for the I-th characteristic direction set.
          • If ICHRRT(I)=1 then the rotation is specified by a rotation matrix, read in row-wise.    (3(3E15.9),/)
          • If ICHRRT(I)=2 then the rotation is specified by the Eulerian angles PHI, THETA, PSI.    (3F14.0)
        • 5.2.3. ALPHA(I,1), ALPHA(I,2), ALPHA(I,3)    (3E15.9)

          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 ICDBIN = 0:
          • 5.3.4. CHRGE(I) The total charge of the I-th distribution.
          • 5.3.5. MAXS(I) records: POLE (I,J)
              POLE(I,J) is the multipole moment of the spherical harmonic of order J for charge distribution type I (in increasing order).    (E15.9)
          • 5.3.6. (MAXS (I))*(MAXS(I)+1)/2) times:S1, S2, S3,    (3E15.9)
              the three components of each characteristic direction. The characteristic directions belonging to the same harmonic order are to be grouped together and read in increasing order.
        • 5.3.7. NH: the number of monopole fileds    (I3)
        • 5.3.8. NH times: HPOLE,L,M,N    (E15.9,3I2)
            HPOLE is the value of a Cartesian multipole moment of order . Each monopole field is calculated as HPOLE * [a directional derivative whose characteristic direction set consists of N, L, M copies of e1, e2, e3, respectively]
    6. 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.

    7. 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: 0-1-2: 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: 0-1: Rotation is specified by Eulerian angles (see 7.1.2.1. below) or a rotation matrix.
        • ITWROT: 0-1: 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 charge-distribution 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 quantum-mechanical 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: ALPHA11, ALPHA22, ALPHA33, TYPE,   (3E15.9,I3)
            The polarizabilities belonging to the distribution type TYPE in the subsequent calculations will be ALPHA11, ALPHA22, ALPHA33.
    8. Configuration specification: If ICINP > 0 then configuration specification by atomic coordinates is read from unit 8.
      • 8.1. If ICINP=1 or 3 then the data is in alphanumeric format:
        • 8.1.1. NGR, IAU, ILAB    (2I3,18A4)
          • NGR: the number of molecules,
          • IAU=0: input coordinates are already in atomic units,
          • IAU=1: coordinates are in (and will be converted to a.u.).
          • ILAB: Dataset descriptor.
        • 8.1.2. IPBC,EX,EY,EZ,EPWA,VPWA    (I3,3F15.0,2E15.8)
          • IPBC: Periodic cell type. Only IPBC=0 (rectangular cell) is used.
          • EX, EY, EZ: The edges of the periodic cell in .
          • EPWA, VPWA: The pairwise additive energy and virial sum (in a.u.) to be used for comparison and correlation.
        • 8.1.3. 3*NGR times : X(I),Y(I),Z(I)    (3F20.9)
            X(I),Y(I),Z(I): Coordinates of the I-th atom. It is assumed that the atoms are in the same order as in the charge-distribution file, (i.e., O,H,H corresponding to the file attached for water).
      • 8.2. If ICINP=2 or 4 then the data is in binary format and was written (by a Monte Carlo program, see [11]) as follows:
        • 8.2.1.
          WRITE (8) NGR,IDENT1,NCONF,NMCMIN,IPBC,EX,EY,EZ,EPWA,VPWA
          
          • For IPBC,EX,EY,EZ,EPWA,VPWA, see 8.1.2 above.
          • IDENT1: 4-character alphanumeric identifier.
          • NCONF: Configuration number.
          • NMCMIN: Monte Carlo stepnumber (discarded).
        • 8.2.2.
          WRITE (8) ((R(I,J), I=1,3), J=1, 3*NGR)
          
            R(I,J): The I-th coordinate of atom J in a.u.
      • If IREADFLD > 0: NGR records specifying the components of the electric field at the group centers    (3F20.9)

        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).

    3. The output of the program

    The output printed on unit 6 is as follows:

    1. Program name, title, date.
    2. Charge distribution data read in. The characteristic directions are printed in the local (principal axis) coordinate system, that is, after the pre-rotation (if ICHRRT(I)>0)
    3. Geometry data read in (coordinates, rotation matrices).
    4. 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.
    5. 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;
    6. If IEOUTL=0 then the contribution of each spherical harmonic order to the permanent moment energy.
    7. 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.

    8. 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: 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.

    A c-shell 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 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

    1. Title of the first translation lattice. (20A4)
    2. 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.  
      • 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 {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.

    4. ((ACOMP(I,J), J=1,3), 1=1,3)    (3F20.0)
        ACOMP(I,J): the J-th component of the I-th unit cell basis vector. These can either be scaled or unscaled [7]
    5. 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 = XNSITE-XNLATT .
      • 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 site-lattice 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:

    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

    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 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.

    2. Detailed description of the input
    1. Title    (20A4)
    2. 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 non-cooperatively as well as cooperatively.
    3. R0    (F20.0) : Scaling factor for the distance.
        Each crystal constant K(NU,Xc-Xj) is divided by R0^(NU1+NU2+NU3+1). This adjusts for any scaling of the input crystal constants.
    4. IRTYP(I), I=1,NUNIT    (24I3) : Index numbers for the different types of rotation matrices by site.
    5. ICTYP(I), I=1,NUNIT    (24I3) : Index which defines the type of charge distribution for site I.
    6. NRTYP    (I3) : Number of distinct rotation matrices.
    7. (((ROT (ITYPE,I,J), J=1,3), I=1,3), ITYPE=1,NRTYP)    (3F20.0) :
        The rotation matrices by type, listed row-wise.
    8. NTYP    (I3) : Total number of types of (i.e., different) charge distributions.
    9. 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 J-1;
      • (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).

    10. NDUPEN    (I3) : Number of elements in the list (item 11) which specify duplicate site-lattice interaction energies.
    11. NDUPEN lines.
        ISD, ILD, ISC, ILC    (4I3) : indices such that
        • Up(ISB, ILD) = Up(ISC, ILC) where the first variables, ISC and ISD, are indices for sites and the second, ILC, ILD are indices for lattices.
    12. NDUPFLD    (I3) : Number of elements in the list (item 13) which specifies duplicate electric field components.
    13. NDUPFLD lines.
        ISD, ILD, K, ISC, ILC, ISG    (6I3)
        • ISG is a constant for which the Kth components of the permanent multipole fields,
        • Ep defined by the translation lattice ILD(ILC) at the site ISD(ISC) satisfy the equation Ep(ISD,ILD,K) = ISG* Ep(ISC,ILC,K).
    14. 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:
        • (i) Title of the problem    (20A4)
        • (ii) ISITE, ILATT, EPS    (5X,I4,12X,I4,6X,F10.5) :
          • ISITE: Lattice site index;
          • ILATT: Translation lattice index;
          • EPS: Ewald parameter;
        • (iii) RT(K), K=1,3: RT = XISITE - XILATT ; (3X,3E20.12)
        • (iv) NZ: Number of non-zero crystal constants for < ISITE, ILATT > (I3)
        • (v) RKN(J) , J=1, NZ:    (3(E20.12,I3,1X))
          • RKN(J): the crystal constant corresponding to NU which has been calculated as in CRYSCON and assigned an index INX(J) as a function of <NU1,NU2,NU3> ;
        • (vi) NCRTTR    (I3) : Number of crystal constant sets that are to be obtained by transformations from the constants for one <ISITE,ILATT> pair to another <ISITE,ILATT> pair;
        • (vii) NCRTTR lines:
            NSITDUP, NLATDUP, NSITORG, NLATORG, ITR    (5I3)

            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.]
    15. 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.
    16. If ITABRD > 0 AND IUNFTB = 0 :
      • a) Title of the Up and Ep tables;    (20A4)
      • b) The Up(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 Ep(ISITE,ILATT):
          ((SLFIELD (ISITE,ILATT,K),K=1,3),ILATT=1,NUNIT)    (3E24.12) : the K-th component of the electric field defined at the unit cell site ISITE by the permanent multipoles at the sites of the translation lattice ILATT.
    17. 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, site-lattice 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

    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
    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 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

    Sn<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, dn<A,B>, may also be given. In this case the lattice energy contribution

    SUM{n} SUM{<A,B>} dn<A,B> Sn<A,B>

    will also be computed.

    2. Detailed description of the input

    1. Program Title.    (20A4)
    2. NUNIT    (24I3) : Number of sites in the unit cell.
    3. 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 {ik} represent a site <i1a1, i2a2, i3a3> where ak are the crystal basis vectors (see 5 below).
      • a) L1MIN, L2MIN, L3MIN    (24I3);
      • b) L1MAX, L2MAX, L3MAX    (24I3).
    4. R0: A scaling factor for distance    (3F20.0):
        The basis vector components in item 5 will be multiplied by R0.
    5. ((ACOMP(I,J), J=1,3), I=1,3)    (3F20.0) :
        ACOMP(I,J); The J-th component of the I-th unit cell basis vector.
    6. IMOLTYP(I),I=1,NUNIT    (24I3) : The type of the molecule at site I.
    7. (RT(I,J),J=1,3), I=1, NUNIT):
        RT(I,J): The J-th component of the position of the I-th unit cell site.
    8. IRTYP(I), I=1, NUNIT    (24I3) : The index which identifies the rotation matrix for unit cell set I.
    9. NRTYP    (I3) : The number of distinct rotation matrices.
    10. (((ROT(ITYPE,I,J), J=1,3) I=1,3), ITYPE=1,NRTYP)    (3E15.9) :
        ROT(ITYPE,I,J): The element in the I-th row and J-th column of the rotation matrix whose index is ITYPE.
    11. NGR: the number of different types of groups (molecules)
    12. 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 J-th local coordinate of the I-th 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 (dn).
      • 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

    I. INSTALLATION

    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:

    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 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).

    J. APPENDICES

    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)

    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 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.

    Let:

    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:

    The forms of the three rotation matrices follow immediately. First, the elements of a matrix C for a rotation about e3(i) by a positive angle ALPHA are:
    cos  (ALPHA) -s( (ALPHA)0
    C = s( (ALPHA) cos  (ALPHA)0(A4.4)
    001
    Both B1 and B3 have the form C. Second, the elements of a matrix D for a rotation about e1(i) by a positive angle THETA are:
    100
    D = 0cos  (ALPHA) -s( (ALPHA)(A4.5)
    0s( (ALPHA) cos  (ALPHA)
    Thus B2 = D. Finally, since the matrix B of equation (III.1) is a right matrix multiplier, B = B1B2B3. According to equation (III.3), BT = B1T B2T B3T 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. BT, 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

    1. Description of the Maxwell formalism:
      E.W. Hobson, "The Theory of Spherical and Ellipsoidal Harmonics", Cambridge Univ. Press, (1931).

    2. Algoritm for calculating the characteristic directions:
      M. Mezei and E.S. Campbell, J. Comp. Phys, 20, 110 (1976).       Manuscript (PDF)       Manuscript (PS)

    3. Methods of partitioning the electon density:
      M. Mezei and E.S. Campbell, Theoret. Chim. Acta, 43, 227 (1977).

    4. Recursion formulae for calculating the directional derivatives:
      E.S. Campbell and M. Mezei, J. Comp. Phys., 21, 114 (1976).       Manuscript (PDF)       Manuscript (PS)

    5. System of equations to solve to obtain induced dipoles:
      E.S. Campbell, Helv. Phys. Acta, 40, 367 (1967)).       Manuscript (PDF)       Manuscript (PS)

    6. 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)

    7. 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)

    8. Brief description of the Maxwell package:
      M. Mezei and E.S. Campbell, J. Comp. Phys., 29, 297 (1978).

    9. Herbert Goldstein, "Classical Mechanics", Addison-Wesley Publishing Company, Reading, Mass., First edition, 1950, PP. 107-9; second edition, 1980, pp 145-7..

    10. W.H. Press, B.P. Flannery, S.A. Teukolsky & W.T. Vettering, 'Numerical Recipes', Cambridge University Press (1986).

    11. Effective cooperative water model:
      M. Mezei, J. Phys. Chem., 95, 7042 (1991).

    12. H. Popkie, H. Kistenmacher, and E. Clementi, J. Chem. Phys., 59, 1325 (1973).

    13. S.P. Liebmann, and J. W. Moskowitz, J. Chem. Phys., 54, 3622 (1971).

    Return to contents page