II. A program for computing normal modes of molecules, crystal phonon dispersion relations and structure factors for neutron inelastic scattering

Computer Physics Communications 14 (1978) 219—243 © North-Holland Publishing Company

II. A PROGRAM FOR COMPUTING NORMAL MODES OF MOLECULES, CRYSTAL PHONON DISPERSION RELATIONS AND STRUCTURE FACTORS FOR NEUTRON INELASTIC SCATTERING Flemming Yssing HANSEN Fysisk-Kemisk Institut, The technical University of Denmark, DK 2800 Lyngby, Denmark Received 21 July 1977 Revised manuscript received 12 December 1977

PROGRAM SUMMARY Title ofprogram: FYFRE

theory, symmetry, dynamical matrix, phonon, annihilation, creation

Catalogue number: ACXR Nature of physical problem The calculation of normal modes in molecules and crystals based on a harmonic potential field defined in terms of any set of coordinates. For crystals the structure factors and

Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland. (see application form in this issue)

partial differential cross sections for inelastic one phonon

Computer: IBM 370/165 Installation: NEUCC, The Technical University of Denmark

annihilation processes and one phonon creation processes are calculated taking temperature effects into account. Also,

Operating system: OS/MVT, Release 21.7

the polarization factors for each atom in the unit cell is calculated for each normal mode.

Programming language used: FORTRAN IV

Method of solution

High speed storage required: 57500 words

The calculations basedgenerates on the results of program described in part are I, which a complete set ofFYCOOR coor-

No. of bits in a word: 32

dinates and force from a very inputunits of basis coordinates. Thesematrices results are stored onlimited permanent and

Overlay structure: none

used in the present program, which calculates the dynamical matrix by simple matrix operations. Symmetry coordinate

No. of magnetic tapes required: none

matrices used frequencies to block diagonalize the dynamical matrix. Once theare normal and eigenvectors have been cal-

Other peripherals used: card reader, line printer, magnetic

culated the inelastic one phonon structure factors and related quantities are calculated.

•disc (6 files) No. of cards in combined program and test deck: 1284

Typical running time

Card punchingcode: EBCDIC

The running time depends on the problem. For the two test runs the running times including the compilation time were 16 sand 18 s. The compilation time is 12 s.

CPC Library subprograms used Cat. No Title

Keywords: solid state physics, crystal, molecule, normal

Unusual features of the program The storage requirement depends on the problem that is the number of atoms in the unit cell or in the molecule, the number of basis coordinates and the number of updated force constants. In the present set up the storage requirement

mode, neutron, structure factor, dispersion relation, group

is set to 57500 words, which is estimated to be sufficient for

ACXQ ACXS

FYCOOR FYADJ

Ref in CPC This issue This issue

219

F. Y. Hansen / FYFRE

220

a number of cases. Up to 20 basis coordinates, 10 atoms per unit cell or per molecule and 30 updates of force constants may be considered. For a more precise estimate of the storage requirement one should consult the long write up. Detailed information on how to change storage requirements

may also be found there. Restrictions on the complexity of problems The only restrictions of the program is the high speed storage available.

LONG WRITE-UP 1. Introduction The FORTRAN IV program, FYFRE, computes certain quantities which are of interest in the theory of harmonic oscillations of a molecule, a crystal and in the theory of inelastic neutron scattering from a crystal. These quantities are the dynamical matrix, eigenvalues and eigenvectors of this matrix, and corresponding structure factors for coherent single phonon scattering of neutrons. The eigenvalues lead directly to the normal modes of a molecule and to the dispersion relations for a crystal. The structure factors are useful in planning and interpreting certain neutron scattering experiments. It is also an important part of the calculation of cross sections for inelastic scattering of neutrons from crystalline or polycrystalline materials. The structure factor depends on the ternperature through the Debye Wailer factor and the partial differential cross section depends on temperature partly through the Debye Wailer factor and partly through the population factor. For the calculation of the Debye Waller factor, the complete dynam. ical spectrum of the crystal has to be known, but in practice, for high temperatures, that is compared to the Debye frequency of the crystal, an approximate expression for the Debye Waller factor may be given in terms of the Debye frequency and the temperature. The calculations are based on the results of the program FYCOOR [2]. Some of the main features of the program will be given here, while a more detailed description may be found in the subsequent sections. _l. The program permits entry of complex wave vector dependent symmetry coordinates used to factor the secular equation and to determine the symmetry of the various modes. The only requirement on the symmetry coordinates is linear independence. The program orthonormalizes the symmetry coordinates by a Gram—Schmidt procedure beginning with the last coordinate. Complex wave vector dependent symme-

try coordinates may be needed in certain directions of non-symmorphic crystals. These coordinates may be found by the projection technique using multiplier representations of the space groups [3]. It is evident that the program also permits entry of real wave vector independent symmetry coordinates. 2. An initial set of force constants must be supplied for the calculation. The force constants are defined in terms of the virtual zero dummy values assigned to each force constant in program FYCOOR and the input force matrices are updated accordingly. The program may also be used to generate the “analytic form” of the dynamical matrix in terms of the contributions from each force constant to the matrix elements. This is done by updating one force constant at a time to the value 1. In that case the program does not calculate the eigenvalues but only prints the dynamical matrix and write the results simultaneously on unit 1 to be defined in the JCL cards, for later use in the force constant adjuster program (part III of this series) [4]. 3. The dynamical matrix is Hermitian, but numerical approximations and round-offerrors may obscure the Hermiticity. Measures are therefore taken to make the matrix Hermitian before the calculation of the eigenvalues and the eigenvectors. The dynamical matrix is in general a complex matrix. Since the subroutine for calculating eigenvalues and eigenvectors is adapted to real symmetric matrices only, one has to use the following expedient in order to obtain the eigenvalues and the eigenvectors of the dynamical matrix. The complex (n X n) matrix is expanded to a (2n X 2n) real symmetric matrix, which can be shown to have the same eigenvalues and eigenvectors as the original matrix, only the multiplicity is twice the one found for the original matrix. 4. Each normalized eigenvector of a Hermitian matrix contains an arbitrary factor e~,where a is real. The program selects this factor to make the first nonzero entry in each eigenvector to be real and positive.

F.Y. Hansen /FYFRE

The complex entries in the eigenvectors are printed as r and 0 corresponding to the complex number r~eW with 0 ~ 0 < 180. 5. The structure factor for each normal mode may be calculated with and without the influence of temperature through the Debye Wailer factor for both single phonon annthilation and creation processes. The polarization factors for each atom in the origin unit cell are also calculated. If temperature effects are considered, a calculation of the partial differential cross section for neutron scattering may be performed both for the annihilation process and the creation process. This program may be considered as an extension and generalization of earlier programs [1].

2. Theory

221

nates generated from the basis coordinates belonging to cell n~,D~P,for a displacement of the relevant atoms. For each set of values off, I, k, n~,a sequence of matrices

[S~It[S,It,ib° 1]

(11.3)

,

[S1It~] [SkItk]D”~

and

(11.4) (11.5)

F’0~~ ,

have been generated and written on units 2, 3 and 4 respectively by the program FYCOOR ready for use in the present program. According to the Born v. Kármán conditions used in ordinary normal mode calculations on crystals, the displacements of the atoms in cells n~,and flq are related by

2.1. Dynamical matrtc In part I [2] it was shown that the potential energy per unit cell may be written (eq. 1.7)

~‘flp = exp(iq. Rnnnq) Xnq~ (11.6) q is a wavevector in the first Brillouin zone and Rnpnq = Rn,, Rnq where Rn,, and R~ are the posi—

2E~0~ =

~ ([S~It~J [s11t,]Do

-

j,k

flp

I

Xo)cTFik Oflp

7P~ X). (11.1) ([SIIt,][SkItk]b’ In this expression [S 11t,]denotes the ith symmetry 1P symbolizes a set of operation of the group and D’ D-matrices, whose rows contain the cartesian displacement coordinates of the atoms for the various basis

x

1’

tion vectorsfor of the cellsdot n~,product and flq. In [2] an explicit expression in part (11.2) was given (1.3) that and from that expression and eq. (11.6) it is seen

[SjIti][SkItk]~~X~ = I =

~

D~exp(iq.Rnq) exp(iq R,~,)X .

0 coordinates in terms of which the force field is defined. The suffix n~indicates that the basis coordinates belong to cell n~,and the bar reflects that the symbol represents a set of displacement patterns, Xis a set of matrices whose elements are arbitrary cartesian displacement coordinates of the atoms in jk the crystal. The elements of the force matrix F0~~ are the force constants coupling coordinates, which are generated from the basis coordinates by respectively the /th symmetry operation and belong to the origin unit cell, and by the kth symmetry operation and belong to cell ne,. The dot product (see eq. 1.2)

[S11r1] [SkItk]b~P’

X,

(11.2)

simply reflects the excitation of the various coordi-

~l1

.

(11.7)

-

D~kdenotes the D-matrix associated with the displacements of the atoms in cell flq and generated from the set in eq. (1.3) by the combined action of the kth and the ith symmetry operation. I is the number of cells involved in the definition of the set of internal or basis coordinates. From (11.7) it is now seen that the expression for the potential energy (11.1) briefly may be written 2E~t= X~TM(q)X,,,

(11.8)

where M(q) is a (3r X 3r) Hermitian matrix, whose elements depend upon the direction and magnitude of the wavevector q, and r is the number of atoms in the origin unit cell of the crystal or in the molecule. In the case of a molecule, the wave vector is of course

F. Y. Hansen / FYFRE

222

always zero. The kinetic energy per unit cell, expressed in terms of the cartesian displacement coordinates of the atoms in the origin unit cell, is

method is IGF Ew2 I = 0

(II 11 ‘

where F is the force matrix and G the inverse kinetic 2Eirjn

=

XOCT mX

0

(11.9)

.

-

.

mis a (3r X 3r) diagonal matrix with m1 m~, m,. occurring successively three times each down the diagonal. m1 is the atomic mass of the ith atom in the cell. Then, by Lagrange s equation and the expressions for the kinetic and potential energies, the secular equation of the problem becomes 2I = 0 or IM(q) mw —1 2 Im M(q) Ew I = 0 (11.10) which gives the dispersion relations. E is the (3r X 3r) unit matrix and it is seen that there will be 3r eigenvalues and 3r eigenvectors. It is important to notice that both the potential energy and the kinetic energy have been expressed in terms of the cartesian displacement coordinates of the atoms in the unit cell, and not in terms of the internal coordinates (e.g. valence coordinates etc.). We therefore always find 3r frequencies no matter how many coordinates havebeen considered in the definition of the force field. The force constants are all defined in terms of such basis coordinates, and in part I [2] it was described how they have to be chosen in order to reflect the symmetry of the problem correctly. In the case of a molecule with r atoms, we therefore always have to find 6 zero frequencies, three for the simple translation and three for the simple rotational motion of the molecule. This may be accomplished by always chosing the basis coordinates to be orthogonal against simple translation and rotation. In the case of a crystal we always have to find three zero frequencies corresponding to simple translation of the crystal, while the simple rotation of the atoms in the origin unit cell has been transformed to internal coordinates with non-zero frequency. Our method of -expressing the potential and kinetic energy [5] differs from the one usually applied and described by ...,

-

.

-

.

—

.

.

.

...,

.

.

-

—

energy matrix. The method applies mass adjusted displacements. For a molecule with r atoms both G and F will be (3r 6) X (3r 6) matrices giving (3r 6) non-zero frequencies, and in the case of a crystal the matrices will be (3r 3) X (3r 3) matrices giving (3r 3) non-zero frequencies. However, often more than (3r 6) respectively and are coordinates considered, in so the linear (3r a crystal, 3) internal casedependencies of a molecule among the internal coordinates exist. In the Wilson method this gives rise to relations between force constants associated with linear dependent coordinates, and this redundancy problem is considered in detail by Wilson et al. [6]. In summary, we avoid considering redundancies by always expressing the potential and kinetic energies in terms of the 3r linear independent displacement coordinates at the expense of having a secular equation of the order 3r instead of one of the order (3r 6) in the case of a molecule and (3r 3) in the case of a crystal. In fact we may also reduce the dimensionality of the secular determinant (11.10), if the potential energy and the kinetic energy are expressed in terms of the 3r linear independent symmetry coordinates for the wave vector considered. A procedure for construction of symmetry coordinates, which may be used for all groups and all wave vectors q is described elsewhere [3]. The symmetry coordinates are displacement patterns of the atoms in the origin unit cell or in the molecule, which span the group and transform according to the various irreducible representations. Let us denote by S, the ith symmetry coordinate with the bar reflecting that the symbol represents a displacement pattern of the atoms. The representation of the symmetry coordinates in the local cartesian coordinate systems (see for details ref. 3) is written as A 1,1 r~ ~ ~i 1 r•11’ .~,arJ , .

—

.

—

—

.

—

—

—

—

-

—

—

—

.

Wilsonthe et kinetic al. [6]. and In the FG-matrix both the Wilson potential energies method are ex- [6] pressed in terms of internal coordinates such as valence coordinates, and not in terms of the cartesian displacement coordinates of the atoms in the origin unit cell. The secular equation in the Wilson FG-

. .

.

‘-‘I

“.,

‘-~3rJ

L

A3r,1.

[I

...

IA

k

.

.

A3r,3r

(11.12)

r

where the elements of the ith column of A are the

223

F.)’. Hansen /FYFRE

components of S~in the local cartesian coordinate systems. The arbitrary displacements of the atoms in the unit cell, denoted briefly X0, may be written as

3r to 3r 3. Another important consequence of expressing the potential and kinetic energies in terms —

r 1

U1,

=

kr]

...,

=

[l~

...

kr] X0

.

(11.13)

IX3r —

X0 may now be expressed in terms of the 3r linear independent displacement patterns S. and from (11.12) and (11.13) it is seen that —

—

Xo

—

— ...~

=

...,

(11.14)

problems are set up more or less manually, one may easily make mistakes, which are very difficult to trace, and this has obstructed the application of other programs. The program FYCOOR has eliminated these problems, and we have always obtained

q~

a perfect no matter which basis coordinates have beenblocking considered.

—i

—

ES~, ~

A Xo S3~] s, 53r]

-

50 =

As

/

The elements of s are the components of ~, in terms of the symmetry coordinates. Introduction of (11.15) into the expressions for the potential energy (11.8) and for the kinetic energy (11.9), gives 2E

.CTACT~. S S,

~II 16’

—

=

sC’TAC’TM(q) As.

(11.17)

-

kin

‘

/

and 2E

ot I’

-

The secular equation now becomes 2AC’TmAI = 0 or AC’TM(q) A W —

IACTm_1M(q) A

—

Ew21 = 0,

of symmetry coordinates is that the dynamical matrix m~M(q) is blocked into smaller matrices according to the symmetry of the wave vector considered, so coordinates of different symmetry do not couple. This simplifies the computation of the eigenvalues and the eigenvectors, and serves as a very critical check of the problem set up. That is, if the symmetry of the potential energy is wrong, it will show up as a wrong blocking of the dynamical matrix. When the

Long range and intermediate range Coulomb interactions may also be considered and in the “dipole coordinate” approximation of A.W. Solbrig [7], for example, the contribution may be expressed in terms of a matrix, which merely has to be added to the dy..

.

namical matrix, so it may easily be added to this program. One may also be interested in frequencies of certam symmetries rather than in all frequencies. In such cases the A matrix has to include only symmetry coordinates of the proper symmetry, so the number of columns in Amolecular becomes crystals less thanare 3r.studied This is for often the case, when example and only intermolecular modes are of interest. -

(11.18) 2.2. The stnicture factor

where A is unitary. In the case of a molecule one may now project 3 symmetry coordinates from simple translational displacement of the atoms and 3 symmetry coordinates from simple rotational displacement of the atoms. Since the coordinates, as defined in terms of the D-matrices, are orthogonal against simple translation and rotation, the secular determinant will block out into a (6 X 6) block giving 6 zero eigenvalues and a (3r 6)X(3r 6) block with non-zero eigenvalues, this in effect having reduced the dimensionality of the secular determinant from 3r to 3r 6. For a crystal the dimensionality of the secular determinant may in the same way be reduced from —

—

In the inelastic phonon scattering of neutrons either a phonon may be annihilated (neutron energy gain) or created (neutron energy loss). The inelastic one phonon scattering cross section for a general crystal lattice for an annihilation process may be written [8] I_d2o_\ ~d~2dE’)

Ik’I ~

~m.

—

x

r ~ 11

/

1 w.(q) I

b

exp [—W,(K)] m~

224

F.)’. Hansen /FYFRE

X exp(is~‘r.) Ke~(q) 2 n1(q) ~ [w + w1(q)] Ik’I =N~—1./ F1(q, ~)2 n,(q) ~ [w + w1(q)] with

,

(11.19)

cross sections for phonon annihilation and phonon creation processes satisfy the detailed balance condition. N 2is isthe number of atoms in the sample and called the structure factor. hw is the F1(q,ic)change of the neutrons. energy The i-summation in the structure factor extends

—11w = 11w 1(q), and —icq—

r,

or

E’ = E +Ilw.(q) k’ k 0

+

q

—

(11.20)

and

over all atoms (r) in the origin unit cell, b~is the coherent scattering length of a bound atom i, m1 the atomic mass of atom i and r• the position vector of .

atom i. W1Ø)that is the Debye factor atom i. It is noticed Fj(q ~ )2 Waller depends uponofthe temper-

r

(11 21)

ature through factor onlyasand furthermore, thatthe thisDebye factorWailer does not appear a simple

The corresponding expressions for the phonon creation process are I d~cr \

scaling factor. The relative contributions from the various atoms may change with temperature. The Debye Waller factor may, in principle, only de deter-

L~2dE’)

(11.22)

mined when the complete frequency spectrum of the crystal is known. However, it can be shown that for all practical purposes it only has a single parameter dependence on theAtspectrum [8], that isW~(K) on the Debye frequency. low temperatures is small

(11.23)

and unimportant and in the high temperature limit an expansion in (kB7)~gives [8] 21K12 3T W(ic)= 11 2 2mav kBOD

X [w ~th ‘Wi ~

11w

=

—

N~—1~1~j(q,x)2(nj(q) + 1) w -(q)]

hw,{q),

,

and

lc=q+t, or E’=E—hw-(q) =

—

q

—

.

and (11.24)

.

In these expressions k’ is the wave vector of the scattered neutrons, k 0 the wave vector of the incident monocromatic neutrons, q a wavelattice vectorvector. in the In first Brillouin zone and t a reciprocal an experiment k’, k 0 and w are known and because q is a vector in the first Briilouin zone both q and t may be determined from (11.24). is theis frequency of the jth (11.21) normal or mode at q w1.(q) and e~(q) the displacement vector of atom i associated with the jth normal mode with frequency w1(q). It is understood that for each /, q and rare determined either from (11.21) or (11.24). n1(q) is the statistical population factor [8,9] given by -

=

exp[hw/(q)/(k~T)]

—

i

(11.25)

where kB is the Boltzmann constant and T the temperature. It is noticed that the partial differential

[

X [I

+

1 /O~\~ 1 ~ ~ —

y~-j

(11.26)

m~,is the mean atomic mass of the atoms in the unit cell and °D the Debye temperature the iscrystal. It 0D the thirdofterm less than is1 per seencent thatoffortheT~’ 0.4 term. With this approximaleading tion, Fj(q, uc)2 may be written 1 2 = exp(—2W(K)) F1{q, ic) w 1(q) r 2 X exp(iic’ r~)K’e~(q) ~—~=

.

(11.27)

3. The structure of the program FYFRE 3.1. General outline The program FYFRE is written in FORTRAN IV and consists of a calling program MAIN and nine sub-

F.~’.Hansen

/ FYFRE

225

routines. They are EIGEN, PSEPOL, STRU, RECIP, MATSC, AMATR, MATMPY, WRIMC and WRIM. Besides these subroutines standard functions and subroutines provided by FORTRAN IV are used. In the

each atom in the origin unit cell are also calculated. Temperature effects on the structure factors through the Debye Wailer factor may be taken into account, and for a given 1k0 I the partial differential cross sec-

following a brief description of the various subroutines listed above is given. Whenever possible reference to appropriate equations in this part and part I is given. The input formats are described in a separate section and the output should be self explanatory and an example is given in the test run of the program. The program uses 6 units, namely the units 1, 2, 3, 4, 5 and 8, which may be defined to be of any convenient type in the JCL cards, which are also included. The D-matrices are read from units 2 and 3 while the force matrices are read from unit 4 in the first calculation with the program, and in subsequent calculations the force matrices may be read from either unit 4 or unit 8. See the MAIN program for details. On unit 1 input data for the force adjuster program may be generated. See below. Control data for themay calculations read gram be run in are three d from unit 5. The pro1fferent modes: 1. Calculation of the normal modes of a crystal or a molecule for a specific set of force constants. 2. Generation of an “analytic” expression of the elements of the dynamical matrix in terms of the various force constants. In this mode one force constant at a time is updated from the value zero to one, while all other force constants are set to zero, and the dynamical matrix is printed without a calculation of the eigenvalues. Elements of the dynamical matrix whose values exceed a certain threshold value, ALIMIT, are also prmted and written on unit 1 for later use in the force constant adjuster program described in part III of this series. This mode is particularly valuable in the analysis of the importance of the various force constants at different points in the Brillouin zone, and if symmetry coordinates are used to block the dynamical matrix, the iterative determination of the force constants on the basis of a set of frequencies is facilitated. In particular the CPU time is reduced and the number of input and output operations are reduced. 3. Calculation of the normal modes and the dynamical structure factors for a given set of force constants and for given ic vectors. The dynamical structure factors are calculated both for phonon annihilation and creation processes. The polarization factors for

tions are also calculated at the given temperature. 3.2. The MAIN program Card types 1—13 are read from the main program and card types 14—15 are read from subroutine STRU. All the input parameters are printed for a check. The two sequences of D-matrices and the sequence of F-matrices are read from the MAIN program. The sequence in (11.3) is read from unit 2, the sequence in (11.4) is read from unit 3, while the sequence of force matrices (11.5) is read from unit 4. The calculation ofM(q) is performed in four steps. First we calculate (see eq. (11.7)) 11(q),

~

D~qexp(iq ‘Rnq)

=

(11.28)

D

for the coordinates generated by the combined action of the ith and the jthsymmetry operation on the set of basis coordinates. The D-matrices for this computation are read from unit 2, where all matrices for the coordinates belonging to the origin unit cell (Rn,, = o) have been stored by the program FYCOOR. Secondly, we calculate (see eq. (11.1)) ..

-.

11

[D

iCT ij,ik

F0~,

(q)j

(11.29)

and thirdly, (see eq. (11.7)) D~exp(iq) (R~ + R~)) = D~(q),

~

‘

q 1

q

q

(11.30)

p

is calculated. The D-matrices for this calculation are read from unit 3, and finally we calculate M~’,~/‘ (q) = [D~’(q)]CT F~’,i~kDik(q) ‘

-

,

(11.31)

‘

where M~’~~(q) is the contribution to M(q) from the coordinates generated from the basis coordinates by the combined action of the ith and the /th symmetry operation and belonging to the origin unit sell, and the coordinates generated by the combined action of the ith and the kth symmetry operation on the basis coordinates and belonging to cell n~,with the force

226

F. Y. Hansen

matrix F~ For a givenj, k, n~the force matrices Fg~hk, where I runs over all symmetry operations, are af’most identical,but as discussed in the write-up of FYCOOR, they may show individual differences. This process is repeated for all interaction types j, k, 11,,, k~,l,, considered. If we have considered NR inter~.

action types and if there is NSYM symmetry operations of the group, the procedure has to be repeated (NR X NSYM) times, and we find M(q)=

~

/ FYFRE fail to be exactly Hermitian. To form a Hermitian matrix and to allow the use of an eigenvalue-eigenvector subroutine for real symmetric matrices, the following operations are performed on m~M(q): Re~[m~M(q)],1} ~ [Re~[m’M(q)]11} +

(II.33b)

IIm[m~M(q)]11I = ~ IIm([m~M(q)]1j}

M~~~(q).

(11.32) -

+

In the 11 0 0 0 interaction type, which is always the first interaction type considered, the force matrices are symmetric, so terms like f,SC,.C5 and fP5CSC,. are included, where C,. is the rth coordinate andf,.3 the coupling constant of coordinate r and s. For all other interaction p’pes the force matrices are not symmetnc, so in order to obtain an equivalent expression to the one 11 0 conjugate 0 0 interaction type, weencountered have to add for the the complex transpose of M(q) for these interaction types to M(q), so the complete M(q) matrix may be written M(q) =

Re { [m’ M(q)]11}],

D M~

0(q),+~

[jLJ1~: ~k(q)

+ (J141/: ik)CT]

i!~1k

‘II 33 a,~

p

M(q) is now multiplied from the left with the inverse mass-matrix according to eq. (11.10), and finally the dynamical matrix, m’M(q), is blocked by the symmetry coordinate matrix, A (11.18). The A-matrix naturally depends on the direction of the wave vector and, in cases where the crystal has screw axes or glide planes in such a direction, it may also depend explicitly on IqI. In such cases it is therefore necessary to add a small subroutine AMATR, which calculates A for each wave vector q. It is generally easy to

set up the relevant subroutine using the method described in ref. 3. There are various possibilities for the output. All updated force matrices may be printed for each wave vector q or only for the first q or the print may be completely suppressed. The A-matrix is printed for each q, while the dynamical matrix may be printed or the print may be suppressed. Details are found in sect. 4. In practice, the computed matrix rn_i M(q) may -

Im~[m_iM(q)]~~}I

(11.34)

.

The sign of the imaginary part is left unchanged. Next write m_iM(q)

=

Re(m’M(q)) + i Im(m’M(q)). (11.35)

Then the expanded real symmetric matrix 1 M(q)) Im(m’ M(q)) [Re(m LIm(m_i M(q)) Re(m~M(q)) —

is formed, whose dimensions are twice the dimensions of the dynamical matrix. It is now easy to show that if A is an eigenvalue of m_iM(q) of multiplicity n then A is an eigenvalue of the expanded dynamical matrix is an eigenvector of multiplicity of the2n, dynamical and if (x1matrix + iy1, correspondx3,.+ iy~,.) ...,

ing to A then (x1, x3,.,y1, ...,y3,.) is an eigenvector of the expanded dynamical matrix and conversely. To eliminate small terms in the expanded dynamical matrix that occur due to round-off errors, the following procedure is applied: 1. Suppose i and! are fixed. 2. Let V represent the element with the largest absolute value in row i. 3. Let U represent the element with the largest abso...,

lute value in column

I.

4. Put W = Max(V, U). 5. If the i, jth element of the expanded dynamical matrix is less than ALIMIT’ W, then this element is set equal to zero. (In the present program ALIMIT = 10~). The procedure is applied to all elements of the expanded dynamical matrix. If the masses are in atomic mass units, force constants in mihidynes/A, then the printed raw eigen-

F.Y.Hansen/FYFRE values are converted by the program to various units as follows:

with P1 =

_____________

w = ~,/raweigenvalues X 1302.77 cm_i =

_____________

w ‘sjraw eigenvalues X0.161526 eV, w = ‘Vraw eigenvalues X 39.08 THz (CPS’

227

~,/x?+y~,

(11.37)

tan1 y1/x1,

0 ~ 0~<360.

(11.38)

Let Pk be the radial component of the first non-zero entry of e so

10i2)

When a raw eigenvalue is negative, a zero is printed for the frequency in the various units. In that case, some of the parameters input to the program may be erroneous or not physically meaningful. The eigenvectors are returned by subroutine EIGEN in cartesian form,but the main program prints eigenvectors in pseudo polar form, which is generated in subroutine PSEPOL and described there. Due to the expansion of the dynamical matrix each eigenvalue and eigenvector occurs twice in the output. That is, if the dimension of the dynamical matrix is (3r X 3r), then 6r eigenvalues and eigenvectors are printed, and they are pairwise The we do not print every secondidentical. eigenvalue andreason eigenvector is that the order of eigenvectors for degenerate eigenvalues may change from one problem to the other.

——

e

,—

—

,

e—~S1...S3,.~e

s,

or 0, 0 0’ 0 S=

Pk~0

03r

Let

(11.39)

.

Pk+1~ °k+1 —

—

Ok

0k

P3r,

0~q~ = 0k+I

—

°k

1

=

1, 2,

...,

n

—

k.

3.3. Subroutine EIGEN

If

This is a standard subroutine for the calculation of eigenvalues and eigenvectors of a real symmetri matnix. It is written in double precision. It may also be used to calculate the eigenvalues and eigenvectors of a complex Hermitian matrix, when the Hermitian matrix is “expanded” as described under MAIN.

~ <0 then 0~+~ = 0~+~ + 360° and if 0’A~.z~ 180° then 0~+j= 0’k~~~I180° —

and Pk+z~ The (ps, O~)pairs are printed with the matrix print subroutine WRIMC. “Ghoast” numbers may occur if the eigenvector does not fill a whole line in the output format determined by the WRIMC subroutine. See the write-up of this subroutine for details. Pk+I =

3.4. Subroutine PSEPOL This subroutine converts the eigenvectors returned from subroutine EIGEN from cartesian form to pseudo polar form. The pseudo polar form is obtained as follows: Let an eigenvector be

3.5.SubroutineSTRU where

The purpose of this subroutine is to calculate the structure factor and related quantities for the inelas-

S

tic phonon scattering of neutrons. Card types 13—15

1 = X1 +

~ ni e =

(p,, 0~),

are read from this subroutine the structure tons are calculated for as manyand reciprocal latticefacvec(11.36)

tors

t

as wanted for a given q wave vector.

F. Y. Hansen / FYFRE

228

In order to compute the various dot products in 2 (11.19) it is necessary to the expression for Fj(q, ic) consider the coordinate systems used to represent the different vectors. In the MAIN program the eigenvectors e’(q) are calculated and represented in terms of the coordinates (S1and S3r), in terms of the symmetry reciprocal lattice vectors, r in qterms of the central cartesian coordinate system (see part I). Here ...

(11.19), F1(q, g)2 ~(w + w1(q)) and for the creation process (11.22), F 2 ~(w w 1(q, ic) 1(q)) are calculated without consideration of the influence of temperature through the Debye Wailer factor. For w~(q) in THz b1 in 10i4 m1 by in atomic mass units and K 2 is m, given in A~,F1(q,x) 2 — 1 b 1374~°24w(q)i=i 1 .

.

we have chosen to represent q and in the central cartesian coordinate system and ther eigenvectors in the local cartesian coordinate systems, whose axes

—

x

exp(iic ‘r 1) ic ~(q) .

are parallel to the respective axes of the central cartesian coordinate system. The eigenvectors are 5i6r Sii (ii,

‘“,

~r

=

~

...,

S3r)

—

:

=

(~i ...

‘ ‘

‘

W(ic)

53r6r

=

72.7635

0D

may

2 (1k,

...,

k3r)

AS

=

=

(I’s, b2, b3) t ,

b11... b13 (b1,b2,b3)=(iYk)

: 1. U3i

°D

in absolute tempera-

wavevector 1k0 I of the incoming neutrons is given, then the partial differential cross sections (11.19) and (11.22) are also calculated. The N factor is not included. For the choice of the various possibilities see input formatscardtype8.

=(i,j,k)B, 3.6. Subroutine RECIP

7.

“

(11.43)

3. If temperature effects are considered and the

where bi, b2, b3 are the reciprocal lattice vectors and

7)33

so q = (ii k) Bt

_(~)4],

when ic is in A—1 and T and ture.

(b 1, b2, b3), t2

—

>< [1 + l(OD)

(11.40)

.

tl

q

(11.42)

1K12 T —

(e =

.

S

and from (11.12) e6r)

barns/eV

2. The structure factors are calculated with the may be written inclusion of the Debye Waller factor (11.26), which

3r)S,

...,

2

. . .

: ~3ri

1, Also

—

The purpose of this subroutine is to calculate the (II 41)

In subroutine RECIP the B-matrix is generated conveniently from the direct lattice vectors a1, a2, a3 represented in the central cartesian coordinate system, as these already have been used in the calculation with the program FYCOOR. With 2(11.40) andevaluated. (11.41) theIndot-products in polariare easily particular the F1(q, zationic)factors ic’ ef(q)I2 for each atom in the origin

components of the reciprocal lattice vectors bi, b2, b3 along the cartesian i,j and K axes from the X, Y and Z components of the direct lattice axes ai, a2, a3. The calculations are based on the relations

b1

a2 x a3 2ir ai a2f~.a3 a
=

2ir

b 2

unit cell of the crystal are calculated and printed. There are now three different possibilities for the calculations and the print of the results: 1. The structure factor for the annihilation process

=

3

1

ai a2 X a3 ‘

=

2~ a1 X a2 a1 a2 X a3 ‘

F.Y.Hansen/FYFRE

so a’ b•

=

2ith-•

The components of the reciprocal lattice vectors are put into the B-matrix (11.41). 3.7. Subroutine MATSC This subroutine is designed to scan the Hermitian dynamical matrix, blocked by the A-matrix, for nonzero elements. The subroutine is used when the program is run in mode 2 and a non zero element is, by definition, an element whose numerical value is larger than the parameter ALIMIT (see card type 3). Only the diagonal elements and upper triangle elements are scanned as the matrix is Hermitian. Two cases of importance may occur: 1) No force constants are known, so all have to be determined in a subsequent iteration. It is therefore necessary to know the analytic form of the blocked dynamical matrix in terms of each force constant separately, and therefore all force constants are updated, one at a time, from zero to one at each q considered. The non-zero complex elements and their positions in the matrix (kth row and lth column) are printed together with the dummy value of the force constant, which have been updated to one. For identification the components of q in reciprocal space are also printed. At the same time the results are written on unit 1 for later use as input to the force adjuster program described in part III of this series. The format on unit I is:

Columns

229

If there is n non-zero elements of the dynamical matrix for F(i)= 1, then the output format above is generated n times with the same q components and dummy force constant, and the relevant matrix elements on unit 1. 2) If some of the force constants are known and therefore do not have to be determined in a subsequent iteration, it is not necessary to determine the analytic form of the blocked dynamical matrix in terms of such force constants. It is only necessary to know their contribution to the dynamical matrix, when updated to the appropriate values. This contribution is, in contradiction to the one mentioned in 1), not characterized by a certain force constant and in order to distinguish this case from the case in 1) a 0.0 is printed, where the dummy values were printed in 1), and this is accomplished by putting the parameter MAT = I (see card type 11). The blocked dynamical matrix is scanned as before and the results are printed and written on unit 1 in the same format as in 1). 3.8. Subroutine AMATR This subroutine is used to generate the symmetry coordinate matrix A, when it depends explicitly on q. The form of A depends of course in general, on the direction of q but not on the magnitude of IqI, unless the crystal in that direction has a screw axis or a glide plane as found in the non symmorphic crystals. The application of multiplier representations of the space groups [3] is in such cases very convenient

Mnemonic

For,nat

Remarks

1—10

Q(l)

F

11—20 21—30 32—4 1

Q(2) Q(3) F(i)

F F F

43—53 55—65

H(k, 1) H(k, 1)

F F

66—68 69—71

Ic

I I

Component of q along reciprocal lattice vector b1 Component of q along reciprocal lattice vector b2 Component of q along reciprocal lattice vector b3 Force constant updated from zero to one identified by the dummy value assigned by FYCOOR. Real part of the non-zero complex matrix element (k, 1) Imaginary part of the non zero complex matrix element (k, I) The row index of the matrix element above. The column index of the matrix element above.

1

230

F. Y. Hansen / FYFRE

and simple, and leads in general to A-matrices, where the (i, f)th element has the form

3.11. Operational information The dimensions of several matrices in program

A,1(q) = A~1exp(iq’ t), (11.44) and where t is a non-primitive translation vector. It is usually very simple to set up a subroutine for calculation of A for a given problem. It is understood that the Iqi independent “skeleton” of the A(q) matrix with elements A11 has been read into the AB-matrix, whose elements then have to be multiplied in a proper way with the phase-factors in (11.44) and the result stored in the A-matrix. A subroutine AMATR illustrating the mentioned features, is included for the p ‘b3 direction in trigonal selenium. See the test example for details. Important: Subroutine AMATR always has to be included in the deck but in cases where the A-matrix has no explicit q dependence, it is not used, so the structure of the subroutine is immaterial. Subroutine AMATR has to be written specifically to the particular problem studied, as it is not feasible to write a general procedure, which may be used in all cases. 3.9. Subroutine MA TMFY This is a standard subroutine for matnx multiplication and is used in subroutine EIGEN.

MAIN and the subroutines are a function of the numben of atoms in the origin unit cell or molecule, the number of basis coordinates considered, the number of updated force constants in each nun and the total number of updated force constants on unit 8. In table 1 is shown how the dimension of the various matrices depends on those parameters together with the maximum values of these in the present set-up. In table 2 is shown which matrix dimensions in which subroutines have to be changed if the maximum value of the parameters mentioned have to be changed. The JCL cards for running the program on the IBM 370/165 computer at the computing center NEUCC are included in the program test deck. From these the definitions of the various units used in the computation may be seen. Table 1 The dimensions of matrices, whose areas may have to be changed from problem to problem. Matrices with the same areas for all problems are not included. The value of the various parameters, which determine the array dimensions, are in the present set up: N4 = 10, max. number of atoms in a unit cell or a molecule. N6 (‘~3X N4) = 20, max. number ofbasis coordinates. MSTO 100, max. number of .

.

.

cumulative updates of force constants on unit 8. Max. value

ofNOCH30

3.10. Subroutines WRIM and WRIMC Subroutines WRIM and WRIMC are standard routines for input and output of real and complex matnices respectively. They are only used for printing matrices in this program. They print 6 columns in a line and at the end of each line a set of three integers are printed. The first and the second integer give respectively the row and the column indices of the first element in the line and the last integer gives the number of elements in the line. As six columns are always printed, each line may contain a maximum of 6 real elements or 3 complex elements, and if the number of elements in a line is less, the last elements are “ghost” numbers. Three alpha numeric characters at the end of each line may be used to characterize the matrix printed. It is very important for these subroutines that the parameters N4 and N6 (card type 2) are in agreement with the matrix dimensions in the program. If this is not the case the output is nonsence.

_______________________________________ Matrix

Dimension

G(20,30) X(30) XSTO(l00) YSTO(100) VAL1(60) VAL3(60) ATOM(30)

(N6, 3 X N4) (NOCH) (MSTO) (MSTO) (6 x N4) (6 X N4) (3x N4)

AMAS(30) D(60,60) VALU(60) VAL2(60) FO’3 B(60 ‘30)

(3 (6 ~ (6 x (6 X (6)( (6 ~

A(30,30)

(3 x N4, 3 X N4)

AA(30,30)

(3 x N4, 3 x N4) (3 X N4, 3 X N4) (6 X N4)

AB(30,30) V(60)

N4) N4, 6 x N4) N4) N4) N4, 3 X N4) N4 3 x N4)

231

F.Y.Hansen/FYFRE

Table 2 The tabie shows which matrices in which subroutines have to be redimensioned, when some ofthe maximum values in table 1 have to be changed from the values in the present set-up. All local defined two-dimensional matrices in EIGEN and MATMPY have the dimension (6 X N4, 6 X N4) and all local defined one dimensional matrices have the dimension (6 X N4). Matrix

MAIN

G(20,30) X(30) Y(30) XSTO(l00) YSTO(100)

N6, N4 NOCH NOCH MSTO MSTO

VAL1(60) VAL3(60) ATOM(30)

N4 N4 N4

S(10)

N4

AMAS(30)

N4

D(60,60) VALU(60) VAL2(60) F(60,30) H(60,30)

N4 N4 N4 N4 N4

B(60,30)

N4

A(30,30)

N4

AA(30,30) AB(30,30) V(60)

N4

STRU

MATMPY

EIGEN

N4

N4

N4

MATSC

AMATR

N4 N4

N4 N4

4. Input format for FYFRE Card type 1

Columns 1—72

Mnemonic

Format

Remarks

TEXT

A4

Text for identification

Mnemonic

Format

Remarks

Read from unit 5

Card type 2

Columns 1—4

5—8

IFEIG

IFSTRU

=1 : Eigenvalues are calculated corresponding to mode 1 and 3 of the program Eigenvalues are not calculated, corresponding to mode 2 of the program

I

=1 : The structure factor is calculated, corresponding to mode 3 of the program =0 : The structure factor is not calculated corresponding to mode 1 ofthe program

IF. Y. Hansen / FYFRE

232

9—12

=1 : The symmetry coordinate matrix A depends explicitly on the wave vector. Add approriate subroutine AMATR to the deck for calculation of A. =0: The A-matrix does not depend explicitly on the wave vector. The AMATR subroutine still has to be included in the deck, but it is not used in the calculation.

IFAMAT

Print of all force matrices after each 13—16

IFFOR

I

=1 :

new update of force constants Print of all force matrices only after

the first update of force constants =0: No print of updated force constants 1: Print of Hermitian dynamical matrix 17—20

IFDYN

I

=0: blocked No printby of A. Hermitian dynamical matrix blocked by A

21—24

N4

I

Max. number of atoms in the origin unit cell or in the molecule. Determines the dimension of certain two-dimensional matrices. See table 1. Equal to 10 in the present set up.

25—28

N6

I

Maximum number of basis coordinates. Determines the dimension of certain matrices. See table 1. Equal to 20 in the present set up.

29—32

MSTO

I

Max. size of the XSTO and YSTO matrices, whose elements give the cumulative update of force constants on unit 8. If the size of the matrices is exceeded a warning is printed. Equal to 100 in the present set up.

Columns

Mnemonic

Format

Remarks

1—4

N2

I

Number of atoms in the unit cell or molecule

5—8

N3

I

Number of basis coordinates

9—12

N5

I

Number of symmetry coordinates (~3X N2)

Read from unit 5

Card type 3

F.Y.Hansen/FYFRE

233

13—16

NRD

I

Number of D-matrices on each of the units 2 and 3 as generated by program FYCOOR. The number may be found in the output from that program.

17—20

NF

I

Number of F-matrices on unit 4 as generated by program FYCOOR. The number may be found in the output from that program.

2 1—24

NSYM

I

Number of symmetry operations in the group of the crystal or the molecule.

25—34

ALIMIT

F

Threshold value for zero and nonzero values of the Hermitian dynamical matrix. ALIMIT is used in subroutine MATSC. Usually = 1.Oe—04.

Columns

Mnemonic

Format

Remarks

1—10

AMAS(1)

F

Atomic mass of atom I

11—20

AMAS(4)

F

Atomic mass of atom 2

61—70

AMAS(19)

F

Atomic mass of atom 7

Read from unit 5

Card type 4

The numbering of the atoms is identical to that in program FYCOOR. Card type 4 is repeated as many times as necessary for inclusion of all atomic masses. Read from unit 5

Card type

Columns

Mnemonic

Format

Remarks

5

1—10 11—20 21—30

ATOM(i) ATOM(i + 1) ATOM(i + 2)

F F F

The X, Y, Z components of the position vector of atom i in the unit cell of the crystal or in the molecule represented in the central cartesian coordinate system (see part I) (in A)

Card type 5 is repeated N2 times Read from unit 5 Omit this card ifIFS TRU=0

234

F.Y. Hansen /FYFRE

Card type

Columns

6

Mnemonic

Format

Remarks

1—10

s(1)

F

61—70

s(7)

F

The coherent scattering length of atom 1. (l014m) The coherent scattering length of atom 7. (

m) Card type 6 is repeated as many times as necessary for inclusion of all atomic coherent scattering lengths. 10i4

Read from unit 5. Omit this card if IFSTRU 0

Card type

Columns

Mnemonic

Format

Remarks

7

1—10 11—20 2 1—30

BASE(1, 1) BASE(1, 2) BASE(1, 3)

F F F

TheX, Y, Zcomponents of the direct lattice basis vector a1 (in A)

Card type 7 is repeated for the other two lattice basis vectors a2 and a3. This card is the same as card 5 for program FYCOOR. Read from unit 5. Omit this card ifIFSTRU =0

Card type 8

Columns

Mnemonic

Format

Remarks

1—10

TK

F

Absolute temperature

11—20

DF

F

Debye temperature of the crystal (K).

1)

21—30 KO F Incident neutron wavevector jk0 I (A— If blank no temperature effect is considered. If TK * 0 and KO = 0 : The partial differential cross section is not calculated. If TK zr 0 and KO * 0: The partial differential cross section is calculated. Read from unit 5 Omit this card ifIFSTR U = 0

Card type 9

Columns

Mnemonic

Format

Remarks

1—10

ABT(l ,1)

F

Real part of the (1,1) element of AT

11—20

ABT(1 ,l)

F

Imaginary part of the (1,1) element of AT

51—60

ABT(l,3)

F

Real part of the (1 ,3) element in

AT

235

F.Y. Hansen /FYFRE

61—70

ABT(1 ,3)

F

Imaginary part of the (1,3) element of AT. Important: For convenience the symmetry coordinates are here arranged in rows in the symmetry coordinate matrix and not in columns, as in the write-up. That is AT rather than A is generated and measures are taken in the program to deal properly with that complication. Notice: The matrix AT is a complex matrix, so each card only contains three elements.

Card type 9 is repeated as many times as necessary for the inclusion of all N5 symmetry coordinates. Read from unit 5

Card type 10

Columns

Mnemonic

Format

Remarks

1—10

V(l)

F

Component of wave vector q along the reciprocal lattice vector b1

11—20

V(2)

F

Component of wave vector q along the reciprocal lattice vector b2

21—30

V(3)

F

Component of wave vector q along the -

reciprocal lattice vector b3 q is a vector in the first Bnillouin zone

Read from unit 5

Card type 11

Columns

Mnemonic

Format

Remarks

1—4

NOCH

I

The number of force constants to be updated from virtually zero to a non-zero value. (Max. value ofNOCH in the present set up = 30)

5—8

NUNIT4

I

=1 : The force matrices are read and updated from unit 4. =0 : The force matrices are read and updated from unit 8.

9—12

NSTO

I

=1: The updated force matrices are stored on unit 8. =0 : The updated force matrices are not stored and may only be used for the wave vector given in card type 10 above.

236

F.Y. Hansen /FYFRE

13—16

MAT

I

=1 : Contributions from several force constants are considered, when program is run in mode 2 (See MATSC)

O; Contribution to the dynamical matrix from only a single force constant is considered at a time when program is run in mode 2. (See description of Subroutine MATSC)

Read from unit 5

Card type 12

Columns 1—10

Mnemonic

Format

Remarks

X(1)

F

The first force constant to be updated. The

force constants are defined by the dummy values assigned in program FYCOOR.

6 1—70

X(7)

F

The 7th force constant to be updated

Card type 7 is repeated as many times as necessary to include all force constants to be updated. Omit this card ifNOCH 0 Read from unit 5

Card type 13

Columns

Mnemonic

Format

Remarks

1—10

Y(1)

F

The update value of the first force constant to be updated and defined by X(l)

6 1—70

Y(7)

F

The update value of the seventh force constant to be updated and defined by X(7)

Card type 8 is repeated as many times as necessary to include all force constants to be updated. Omit this card ifNOCH =0 Read from unit 5

Card

type

14

Columns 1—4

Omit this card ifIFSTR U = 0 Read from unit 5

Mnemonic

Format

Remarks

NTAU

I

The number of reciprocal lattice vectors r to be considered in the calculation of the structure factor

237

F.Y. Hansen /FYFRE

Card type

Columns

15

1—10

Menmonic

Format

Remarks

T(l)

F

Component of the reciprocal lattice vectort alongb

1

11—20

T(2)

F

Component of the reciprocal lattice vector t alongb2

21—30

T(3)

F

Component of the reciprocal lattice vector r along b3

Omit this card ifIFSTRU = 0 Read from unit 5 This card is repeated NTAU times

Card types 10—15 are repeated for each wave vector q considered. The D-matrices are read unformatted from the units 2 and 3 and the force matrices are read unformatted from unit 4 (or unit 8) as defined in FYCOOR.

5. Test example

Let us illustrate the calculations with the FYFRE program on trigonal selenium. We consider the problem set up in part I. The D-matnices and the force matrices generated by FYCOOR have been stored on tape ready for use in this program. The force model corresponds to the one suggested by Nakayama et al. [10]. We have shown calculations with the program both in mode 2 and 3. In the mode 3 calculation the same force constants have been used as those of Nakayama et al. [10], and in our representation the following force constants have to be updated from virtually zero value to the values given below: fi~ritr

fi~ri~r’

are included in the calculation of the structure factor and that the partial differential cross sections are calculated. After the print of all the control parameters, the first wave vector q is printed. Because this is the first wave vector the force matrices are updated from unit 4, and a list over the updated force constants is printed. As we are going to consider q vectors along the screw axes of the Se-spirals (see fig. 2 of part I), the symmetry coordinate matrix A shows explicit q dependence. The A-matrix for this direction is shown in table 3 [3]. It is seen to have the form described in (11.44), and it is very easy to set up a proper sub-

=

1.0101 X 10b0

=

1.18 rndynes/A

=

1.0202 X 10i0

=

0.9634 mdynes/A

=

1.0102 X 10_b

=

0.02784 mdynes/A

=

1.0303 X 10~0= 0.046 mdynes/A

of lineones clearly andeach which are indicate, ghosts which elements are real

=

2.0101 X 10~0= 0.045 mdynes/A

Because IFDYN = 1, the dynamical matrix is printed, and it is seen to be blocked precisely as demanded by the symmetry coordinates in table 3. That is, only cross-terms between symmetry coordinates of the same symmetry are observed. It is noticed that the print of eigenvalues includes 2 X 9 = 18

In the output, the set of control parameters is printed and because IFSTRU = 1, data for the structure facton calculation are entered and printed. The fact that T > 0 and k0 > 0 indicates that temperature effects

routine AMATR. The transpose of A at the considered q is printed. After this the updated force matnices are printed, because IFFOR = 1. All force matrices are (3 X 3) matrices and “ghost” columns are printed as mentioned in the comments to the WRIM and WRIMC subroutines, but the integers at the end

238

F.Y. Hansen /FYFRE

Table 3 The aymmetry coordinate matrix for trigonal selenium for q = Mb3, (0 ~ < 0.5). The elements of the ith column are the components of the ith symmetry coordinate represented in the local cartesian coordinate systems (11.12). The symmetries of each symmetry coordinate are indicated below each column by crosses for: q = 0 : r; q = Mb 3 0 < M < 0.5 ~ q = 0.5 b3 : A, f1 = exp(i 2q~a3/3), 12 = exp(iq’ a3/3), e = exp(i 243), e” = exp(—i 243). 0 0

0 0

0

0

0

0

—~‘J~f2

—ff2

0

~f2

—~/~f2

0

0

0

~vif2

0

0

—~‘J~fic ~fie

—~‘J~fi ~Ii

0

0

0

0

X

x

1_’i

x

r

2 r3 A1 A2 A3

x

~‘J~fi x

0 0

0

0 0

0

~~4,~J~f2e* _~f2*

0

—~~Jif2e—~f2e

0

_~sJ~f2e* 0

if2e

—~‘J~f2e 0

0

0

0

~~Jif2e*

0

I~f1~*if1e*

4~J~f2e 0

~J~f~e

0

_~fie*

_~~J~f1e* 0

0

0

~J~f1e

0

0

x

x

x X

X

X

x

x

x

x

x

0

i~1

0

0

4fie

x x

x

x

x X

X

X

X

x

eigenvalues and that each eigenvalue is duplicated, as described in MAIN. The (trace — sum) calculation should give a small number and serves as a check on the eigenvalue calculation. Also, the eigenvectors are duplicated. Each row corresponds to an eigenvector, so that the ith row is the eigenvector belonging to a normal mode with the ith eigenvalue. The jth element in the ith row gives the contribution of the jth symmetry coordinate to the ith eigenvector. Using table 3 it is therefore very easy to determine the symmetry of the different normal modes, For simplicity only one reciprocal lattice vector = 1 b3 has been considered in the calculation of the structure factors and partial differential cross-sections. In the output the representations of t, q and i~(both for the phonon annihilation (11.20) and creation (11.23) process) are shown in the reciprocal lattice vector system b1, b2, b3 and the central cartesian coordinate system i,j, k. Then for each of the 18 2norfor ma! thethepolarization ~(q)Itogether each modes atom in origin unit factors cell areI K printed

x

with the structure factor and partial differential crosssection, both for the phonon annihilation and creation processes. In cases where card type 8 is blank only the polarization factors and the structure factors without the inclusion of temperature effects are calculated. Finally an example of the output from the program in mode 2 is shown. To include all possibilities for that mode, it is assumed that force constants 1.0101 X 10’° and 1.0202 X 10—10 are well known, while the other three force constants are going to be determined in a subsequent iteration calculation with FYADJ described in part III. The same control parameters as before are printed, while all data for the structure factor determination is omitted. Then q is printed and the parameters of card type 11 shows that one force constant is going to be updated from unit 4 and that the updated force matrices are not stored on to unit 8. The the MAT parameter are going generate analytic formshows of thethat dy-we namical matrix in terms of a force constant, namely

F. Y. Hansen / FYFRE

1.0102 X 10--iO The A matrix is printed and the updated force matrix and the dynamical matrix are also printed but left out here due to limited space. Finally, all non zero elements of the dynamical matrix is printed in the way described under MATSC. The latter is also printed on unit 1 in the format for later use in the FYADJ program. When force constants 1.0101 X 10—10 and 1.0202 X 10—10 are updated, MAT is equal to 1 and a (0.0) is printed for the force constant. The order in which the data on unit I is generated is arbitrary, as the FYADJ program sorts them out in a relevant way. Acknowledgement

The author wishes to express his gratitude to Dr. H.L. McMurry, Physics department University of Idaho, Pocatello USA, for many stimulating discussions in connection with the present work and for very valuable help in checking out the calculations.

239

References

[1] J.K. Boyter, H.L. McMurry, A.W. Solbrig and R.M. Martin, MACS C, A computing program for crystal vibrations which includes long range Coulomb interactions, nr. P00360 Argonne Code Center, Argonne,Code illinois 60439. [2] F.Y. Hansen, Computer Phys. Comm. (preceding paper). [3] F.Y. Hansen, Physical Review B, 0,0 (in press). [4] F.Y. Hansen, Computer Phys. Commun. (following paper).

[5] H.L. McMurry, A.W. Solbrig and J.K. Boyter, J. Phys. Chem. Solids 28 (1967) 2359. [6] E.B. Wilson, J.C. Decius and P.C. Cross, Molecular Vibrations (McGraw-Hill Book, 1955). [7] A.W. Solbrig, Idaho Nuclear Corporation, The use of valence force potentials in calculating crystal vibrations. II. Dipole coordinate approximation for long- and intermediate range forces, IN-1409 (July 1970).

Marshall and S.W. Lovesey, Theory of thermal neutron scattering (Oxford, 1971). [9] T. Hill, Introduction to statistical thermodynamics (Addison Wesley, 1962). [10] T. Nakayama and A. Odajima, Journal of the Physical Society of Japan, 34, no. 3 (1973). [8] W.

F. Y. Hansen / FYFRE

240

OUTPUT FROM MODE 3 TEST

CALCULATICN CONTROL

OF

IFEIG= I IFSTRU i IFAMAT I IFFCR I IFDYN 1 MAX. NUMBER OF MAX. NUMBER OF XSTO— AND YSTO NUMBER OF NUMBER OF NUMBER OF NUMBER OF NUMBER OF NUMBER OF THRESHOLD

DATA

TIlE NOR’4ALMOOES

SELENIUV.

LATTICE

FACTOR

RECIPROCAL

NEUTRCN

80.00

•

NSTO~

(0=

0.0 1

FORCECONSTANTS

7,C00

•

MAT

1/A

0.250000) 0

UPDATED

(REPRESENTED

L.02020E—l0

9.63400E—01

I.03030E—I0

4.60000E—02

1.01020E—IO

2.784000—02

2.OI010E—IO MATRIX

IN RECIPROCAL

SPACE)

FORCECONSTANTS

1180000

5.77345375—01 —2.4999762E—0I —I.4433634E—01 —3.44128995—08 —4.3301660501 2.5000256E—01 0.0 0.0 0.0 5.7735264E—0I —6.58634425—06 2.88672635—01 0.0 i.4336647E—06 —5.0000006E—0l 0.0 0.0 0.0 5.7735276E—OI 2.4999452E—0I —5.44341775—01 0.0 4.330)3325—01 2.49995335—01 0.0 0.0 0.0

0.0 1.268332

1.65976C 0.0

1.010IOE—10

SYMMETRYCOOROINATE

0.0 4.953838

K

MAVEVECTOR

WAVEVECTOR.0( 0.0 NOCP$ 5 NUNIT4~ I INITIAL

SCOH 0892000 0.892000 R.R92000

K

TEMPERATURE

INCIDENT

MASS 78.959991 78.95999) 7895999I

~

0.0 0.0

3C0.0

Z —1.651299 0.0 1.631139

3.785600 0.0

*

VECTORS

B2 83

I.0000E—04

0.0

—2.185598 0.0

LATTICE

TEMPEPATURE

Y 0.0 0.8533CC —0.853300 4.371200

42 A3

3 3 9 3. 12 N

CALCULATION:

S 0.SRS3CO —0.492600 —0.492600

VECTORS

10 20 100

N5 NJ= N5 NPD~ NF GEQUP,SSYM ALIMIT=

ATOMS IN THE UNIT CELL. BASIS COORDINATES. SYMMETRY COORDINATES. 0—MATRICES. F—MATRICES. SYMMETRYOPERATIONS IN THE VALUE.

ATOM NP I 2 3

OEEYE

TRIGOSAL

ATOMS IN UNIT CELL.N4= BASIS COORDINATES.N6 MATRIX DIMENSION MSTO=

FOR THE STRUCTURE

DIRECT

CF

PARAMETERS:

00

4.50000E_02 AC TRANSPOSE)

0.0 —1.44336045—01 —2.49997385—01 0.0 —2.50002035—01 4.3301624E—01 0.0 0.0 0.0 0.0 2.88672695-01 —6.g43o739E—ob 0.0 5.00000005—01 —1.78010475—07 0.0 0.0 0.0 0.0 I.A434206E—0l 2.49994755—01 0.0 —2.4999887E—0I —4.3301362E—0I 0.0 0.0 0.0

6.8825329K—OS 4.330L654E01 —2.5000256E—0I 5.77345255—01 —2.4494762E—oi —I.4433646501 0.0 0.0 0.0 0.0 —i.4385041E—06 5.00000065—0) 5.7735264E—01 —6.5B91381E—06 2.88672635—01 0.0 0.0 0.0 0.0 —4.3301332E—01 —2.4999833E—01 5.7735276E—O1 2.4999452E01 —1.44341770—01 0.0 0.0 0.0

0.0 2.5000197E—OI —4.3301624E—01 0.0 —1.4433610E—01 —2.49997445—01 0.0 0.0 0.0 0.0 —5.0000000E—0l 7.8284978E—07 0.0 2.88672695—01 —6.94676785—06 0.0 0.0 0.0 0.0 2.49098870—01 4.33013625—01 0.0 —1.44342065—01 2.49994755—01 0.0 0.0 0.0

0.0 0.0 0.0 0.0 0.0 0.0 5.77343345—01 5.0000298E—0l 2.8867739E—O1 0.0 0.0 0.0 0.0 0.0 0.0 5.773553BE—01 I.2815326E—06 —5.7734758E—01 0.0 0.0 0.0 0.0 0.0 0.0 5.7735187E—01 —4.9999696E—01 2.8867811E—0I

0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.8867674E—01 S.0000256E—01 0.0 0.0 0.0 0.0 0.0 0.0 0.0 —5.7134764F—01 2.05639395—06 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.88678770—01 —4.99997380—0)

1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7

7 8 8 8 9 9 9

I 4 7 I 4 7 I 4 7 I 4 7 I 4 7 1 4 7 1 4 7 I 4 7 I 4 7

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT

241

F.Y. Hansen /FYFRE UPCATEO

FORCEMATR ICES:

1.180005 00 2.784000—02 I.01030E—l0 1.180005 00 2.784000—02 1.010305—10 1.181005 00 2.784005—02 1.010305—10 0.0 2.784005—02 1.0)0305—10 0.0 2.76400E02 1.01030510 0.0 2.76400502 1.010305—10 4.500005—02 0.0 2.03010E—I0 4.50000502 0.0 2.030)05—10 4.50000E02 0.0 2.03C10E10 0.0 0.0 2.030105—10 0.0 0.0 2.030105—10 0.0 0.0 2.030105—10 DYNAMICAL

2.754005—02 1010300—I0 0.0 0.0 0.0 9.634005—01 1.020300—IC 0.0 0.0 0.0 1.020300—10 4.600000—02 8.4134300I 7.491300—0) 4.RCIRCE—01 2.784000—02 1.013300—10 0.0 0.0 9.0 9.63400E—01 1.020305—10 0.0 0.0 0.0 1.020300—10 4.600300—02 0.0 0.0 0.0 2.784005—02 1.0103CF—IC 0.0 3.0 0.0 9.634000—01 1.020300—10 0.5 0.0 0.0 1.02030010 4.600005—02 0.0 0.0 ~ 2.78400E—02 1.0)0300—10 C.0 0.0 0.0 0.0 l.02030E—I0 0.0 0.0 0.0 1.020305—10 4.600000—02—6.355570—0I—6.045550—CI—4.SOIMCC—0I 2.78400E—02 10103)E—lO 0.0 0.0 0.0 0.0 1.020305—10 0.0 0.0 C.O 1.02033E—I0 4.60000E—02 0.0 0.0 0.0 2.78400F—02 1.010300—IC 0.0 0.0 0.3 0.0 1.020305—10 0.3 0.0 0.0 1.020300—10 4.600000—02 0.0 0.0 0.0 2.0I320E—I0 2.010300—10 3.0 0.0 0.0 2.020200—10 2.020300—10 0.0 0.C 0.0 2.030205—10 2.030300—10 8.4)3400—01 2.481300—0) 4.801900—01 2010205—IC 2.010300—10 0.9 0.0 0.0 2.020205—10 2.020300—IC 0.3 0.0 0.0 2.030200—10 2.O3030~—IC 0.0 0.0 3.0 2.010205—10 2.010305—10 0.0 0.0 0.0 2.023200—10 2.020300—10 0.0 0.0 0.0 2.030205—101 2.030300—10 3,0 0.0 00 2.0102C0*I0 2.010305—10 3.0 0.9 0.0 0.0 2.020300’.lO 0.0 0.0 0.0 2.030200—IC’ 2030300—l0—5.355570—0I—6.04556E—01—4.R01NCE—01 2.010200—10 2.010300—10 0.0 0.0 0.0 0.0 2.020300—IC 0,0 0.0 0.0 2.030205—10 2.03030E—10 3.0 0.0 0.0 2.010200—10 2.010300—10 0,0 0.0 0.0 0.0 2.020300—10 ‘.3 01.0’ 0.0 2.030200—10 2.030300—Ic 0.) 00 9,3

MATRIX

3.0125540E—02 8.14907255—09 5.6345016E—08 —5.67990365—07 —4.54019745—09 —4.6566)295—10 1.02766545—08 3.6068750E—09 1.04773795—09 8.14907255—09 1.71738535—02 1.02856575—07 —1.19907785—08 —1.62748625—07 —1.65209715—09 —3.64170565—09 9.78005805—08 —1.16415325—09 5.63450165—OS 1.02856575—07 5.6556165E—03 —4.65661295—09 2.9)038305—09 8.87084755—08 —1.5)339525—05 —2.3283064E—10 2.39607745—07

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

1 13 2 13 3 13 1 13 2 1] 3 II I 13 2 113 3 13 I 13 2 13 3 10 I 13 2 13 3 13 I 13 2 13 3 13 I 13 2 13 3 II I 13 2 13 3 13 I 13 2 13 3 13 1 13 2 13 3 13 1 13 2 II 3 13 I 1 1 2 13 313

F F

F 0 F: F:

F r F F F

0 F F F F F F

F F F

F F

F F ~ F F F F F

F F

F F F

BLOCKED BY Al 0.0 4.30736695—09 I.2805585E—09 —6.79303715—03 —1.82422815—07 3.83006415—08 —4.43430245—03 —7.28759915—08 6.28642745—09 —4.30736695—OS 0.0 1.0)281335—08 —1.7764978E—O7 1.2548886E—02 1.97906055—08 —7.48550525—08 8.4851198E—03 6.98491935—09 —1.28056855—09 —1.01281335—08 0.0 7.41565605—08 59953851E—Oa —5.75574495—03 4.4004992E—O8 —8.859206OE—O8 —4.05088445—03

RAW EIGENVALUES 3.373650—02 3.373650—02 3.237700—02 3.237700—02 3.206540—02 3.206540—02 1.326520—02 1,326520—02 7.155810—03 7.155870—03 3.044350—03 3.044350—03 8.825830—04 8.825830—04 3.158930—04 3.158930—04 1.918810—04 1,918810—04

—5.67990365—07 —I .1990778E—08 —4.6566129E—09 3.26660655—03 —7.93952495—08 —5.6577846E—08 —1.56713210—OS —1.30385165—08 —6.8568625E—08 —4.5401976E—09 —1.62748625—07 2.91038305—09 —7.9395249E—08 I.0256369E—02 l.0407530E—07 —5.6401404E—08 2.98196825—03 8.42846930—OR —4.6S66129F—lO —1.6520971E—09 8.87084755—08 —5.05778465—038 I.0407530E—07 8.58671965—03 —6.72644095—08 1.24564395—08 1,13469145—02

WAVE NUMBERS 2.392865 02 2.392865 02 2.344155 02 2.344155 02 2.332855 02 2.332855 02 1.500465 02 1.500465 02 1.10204E 02 1.102045 02 7.188115 01 7.188)15 01 3.870305 01 3.87O30E 01 2.315460 II 2.315465 01 1.804615 01 1.804615 01

6.793037)5—03 I.776497R5—07 —7.4156560F—08 0.0 5.355)0485—39 —3.37604435—39 8.67250950—08 0.0 —1.74622985—OS 1.82422810—07 —1.25488865—02 —5.9953891E—08 —5.3551048E—09 0.0 1.295)2055—09 —1.41290000—09 —1.90578085—01 —8.14907255—ID —3.83006415—OR —1.97906055—08 5.75574495—03 3.3760443E—0O _I.29512050_39 0.0 3.3774086F—09 1.04773790—09 1.02466800—07

ELECTRON VOLTS 2.966830—02 2.966830—02 2.9064310—02 2.906430—02 2.892420—02 2.892420—02 1.860370—02 1.860370—02 1.366390—02 1.366390—02 8.912290—03 8.912290—03 4.798660—03 4.798660—03 2.870860—03 2.870860—03 2.237480—03 2.237480—03

-

1.02786545—08 —3.84170565—09 —1,51339925—09 —I .56113215—05 —5.64014045—OR —6.72644095—08 2.91172600—03 —9,74396245—08 —I.3539102E--07 3.60887505—09 9,78005805—08 .—2.3283064E—10 —1.3O385I65—08 2.98)96820—03 I.2456439E—38 —9.74396240—06 1.R092297E—02 2.4680048E08 1.04773795—09 —116415325—09 2.39607745—07 —6.8S68625E—08 8.42846930—08 l.1346914F—02 —1.35391025—07 2.46800480—08 2.69658495—02

CPS*E12 7.174020 00 7.178020 00 7.031905 00 7.03)905 00 6.997985 00 6.987985 00 4.501035 00 4.501035 00 3.30587E 00 3.305875 00 2.15626E 00 2.10626E 00 1.161005 00 1.161005 00 6.945635—0) 6.945835—01 5.4134IE—01 5.413415—01

4.43430245—03 7.4855052E—08 —4.40049925—08 —8. 672509SF—OS 1.41290005—09 —3.37740865—09 0.0 1.74622985—09 —8.73114915—09 7.28759915—08 —R.4851198E—03 8.85920600—08 0.0 1.90578080—07 —1.04773795—09 —1.74622985—09 0.0 —1.86264510—09 —6.28642745—39 —6.98491935—09 4.05048445—03 1.7462298E—09 8.1490725F—10 —1.02466805—07 8.73114910—09 l.9625451F—09 0.0

I 1 I 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 6 8 9 9 9

I 4 7 1 4 7 1 4 7 I 4 7 0 4 7 1 4 7 I 4

7 I 4 1 1 4 7

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

M M M N N M N N M N M N N N N N N M N M N N N 8

N N N

242

F. Y. Hansen

EIGENVECTORS IN PSEUDO—POLAR FORM 0.0 0.0 0.0 3.0 2.16254830—01 3.0 0.0 0.0 2.16638375—03 0.0 2. 16704255—01 1.47102770 02 5.53575075—01 3.0 0.0 0.0 0.0 0.0 9.63614945—01 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.02638920—01 0.0 0.0 0.0 0.0 0.0 7.02706690—01 0.0 l.36362V3E04 5.58873900 00 0.0 0.0 3.08661760—01 0.0 0.0 0.0 0.0 0.0 3.0e607055—01 0.0 0.0 0.0 0.0 0.0 3.ec47243E—o4 3.0 7.4319)540—01 l.5785310E 02 0.0 0.0 0.0 0.0 7.4288231E01 0.0 2.16488235—02 0.0 0.0 0.0 0.0 0.0 2.16496815—02 0.0 0.0 0.0 0.0 0.0 266502405,OI 0.0 0.0 0.0 0.0 0.0 2.665S048E—01 0.0 0,0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.33241060—01 0.0 0.0 0.0 2.3049061004 3.3 6.33241005—01 7.0036270E 01 0.0 0.0 6.41102490—01 0.0 0.0 -3.0 0.0 3.0 6.4110249E01 0.0 0.0 0.0

C.0 0.0 —4.42244235—0) 0.0 —2.49847750—03 4.416330910—01 —2.2518027E—01 0.0 0.0 —2.2408020E—01 01.0 0.0 (9.0 4.76545810—01 01.0 01.0 4.76486685—01 I.5935I81E—04 0.0 4.64939245—0) 0.0 0.0 4.64944660—01 0.0 9.0 —4.2880257E04 S.0364661E—01 0.0 0.0 —5.040097850I 6.05581820—01 0.0 01.0 6.05531400—01 0.0 0.0 7.6340079E—01 0.0 0.0 7.6334631E31 0.0 0.0 0.0 0.0 7.42111S0001 0.0 —2.6824954E04 7.42111445—01 0.0 —7.4614507E01 0.0 0.0 —7.46145I3E-’OI 0.0

B) RECIPROCAL WAVE

LATTICE

VECTOR.TAU

VECTOR.O

KAPPA

(ANNIHILATION

KAPPA

(CREATION

NORMAL

MODE

POLARIZATION

P48

PROCES)

PROCES)

1

FACTOR

:

ATCM

0.13 0.0 8.39966130 0.0 1.60484160 5.82740480 9.00000150 0.0 0.0 13.00000300 0.0 13.0 C.0 ‘3.99463845 3.0 0.0 9.00429995 4.49093635 0.0 9.00002440 0.0 0.0 9.00086525 “.0 0.0 9.0980728E 6.7-3051510 0.0 0.0 3.97352910 8.99999850 0.0 0.0 8.99909850 3.0 0.0 8.39999850 0.0 0.01 8.999908SE 0.0 0.13 0.0 0.0 9.00000005 0.0 13.00084695 1.69036305 0.0 9.00000I5E 0.0 0.0 (3.99999690 0.0

32

01 0) 01 01

01 01 01 01 ~II 01 01 01 01 01 01 DI

01 01 02 01 01-

83

MODE NB

1

WITH

X

Y

0.0 0.0 8.9)968720 0,0 8.0386761E 5.80335390 9.00000150 0.0 0.0 9.0000000F 0.0 0.0 0.0 8.99993130 0.0 0.0 9.00213170 0.0 0.0 9.0000092F 0.0 0.0 9.0304379F 0.0 0.0 8.66015470 6.78833160 0.0 0.0 8.933461760 9.03000005 0.0 0.0 5.00000150 0.0 0.0 8.99099850 0.0 0.0 8.99999850 0.0 0.0 0.0 0.5 8.99999850 0.0 0.0 1.69036290 0.0 8.99990690 0.0 0.0 8.99999540 0.0

0.0

0,0

1.0000

0.0

0.0

1.2683

0.0

0.2500

0.0

0.0

0.3171

1

0.1

0.0

0.7500

0.0

0.0

0.9512

:

0.0

0.0

1.2500

0.0

0.0

1.5854

I

IN

7.178025 THE

FREQ’JENCY~

01 01 0) 01

01 00 01 01 01 01 01 01 0) 01 DI

01 02 01

PROCES

•

00

31

2.285285—01 2.28526E—OI 2.285265-01

3.908975—12 8ARNS~EV AT 7. 300.0 K 1.851485—12 RARNS/EV AT T~ 300.0 K

THZ

POLARIZATION FACTOR FOR ATOM * IN TIlE CREATION PROCES • POLARIZATION FACTOR FOR ATOM 2 IN THE CREATION PROCES • ARIZATION FACTOR FOR ATOM 3 IN THE CREATION PROCES • ?kUCTURE FACTOR FOR THE CREATION PROCES PARTIAL DIFFERENTIAL CROSS—SECTION FOR THE CREATION PROCES

6.348005—01 6.347955—0) 6.347955—01 6.270945—01 8ARNSFEV AT T. 7.725465—01 BARNS~EV AT T.

300.0 K 300.0 K

I I I 2 2 2 3 3 3 ~ 4

I 4 7 I 4

4

7 I 4

3 3 3 3 3 3 3 3 3 3 13 3 3 3

7

3 I

PPE POE POE

4

3

POE

7 1

0p5

7 1

3 3 3 3 3

4

3

7 1 4 7 1 4 7 I

3 3 3 3 3 3 3 3

POE POE

4

3

7 1 4 7 1 4 7 1 4 7 1 4 7 I 4 7 I 4 7 1

1 3 3 3 3 3 3 3 3 3 3 3 3 I 3 3 3 3 3 3

5 5 5 6 5 ~ 7 7 7 5 8 8 9 9 9 10 10 10 Il II 11 12 12 12 13 13 13 14 14 14 15 15 15 16 16 16 17 17 17 18

00 1HZ

ANNIHILATION

7.178025

01

z

0.0

POLARIZATION FACTOR FOR ATOM 2 IN THE ANNIHILATION PROCES • POLARIZATION FACTOR FOR ATOM 3 IN THE ANNIHILATION PROCES • STRUCTURE FACTOR FOR THE ANNIHILATION PROCES • PARTIAL DIFFERENTIAL CROSS—SECTION FOR THE ANNIHILATION PROCES NORMAL

0.0 0.0 —8.7043303E—C1 0.0 —5.7835714F04 8.70600825—01 —1.44280495—01 0.0 0.0 —1.45119895—01 0.0 0.0 0.0 5.2R39607E—0I 0.0 0.0 5.28359230—01 0.0 0.0 —8.2979470001 0.0 0.0 —8.29193390—01 0.0 0.0 —0.04944435—04 —4.40461040—01 0.0 0.0 4,40567610—01 —7.95488420—01 0.0 0.0 —7.9548877001 0.0 0.0 5.68061200—01 0.0 0.0 5.8843279E—0l 0.0 0.0 0.0 0.0 —2.19718285—01 0.0 0.0 —2.19718280—01 0.0 —1.79506360—01 0.0 0.0 —1.79596300—01 0.0

01

:

wITH FREOUENCY= FOR

/ FYFRE

7

1 4 7 I 4

I

4

POE PRO

PPE POE POE POE POE

POE POE POE ODE ROE PRO

ODE

PPE 905 POE POE POE OPE PBS POE POE POE ORE POE POE POE

PPE PP5 POE ODE POE POE POE POE POE 000 POE POE POE PPE BPS BPS POE

18

4

3

BPS

18

7

3

P05

F. Y.

Hansen / FYFRE

243

OUTPUT FROM MODE 2 TEST

OF THE

CALCULATION CONTROL

ANALYTIC FORM TIlE

DYNAMICAL

IFAMAT IFFOR~

1 1

IFOYP4.

I

MAX. NUMBER OF MAX. NUMBER OF XSTO— AND YSTO N(JMBER OF NUMBER OF NUMBER OF NUMBER OF NUMBER OF NUMBER OF THRESHOLD

ATOMS IN UNIT CELL.N4= 8ASIS COOPDINATES.N6= MATRIX DIMENSION MSTOO

5 10 100

ATOMS IN THE UNIT CELL. P42= BASIS COORDINATES. N3 SYMMETRY COOROINATES, N5 0—MATRICES. NRD F—MATRICES. NF SYMMETRYOPEPATIONS IN THF 390U°.NSYM~ VALUE. ALIMIT AlUM NR 1

18.959991

2

78.059091

3

78.059991

INITIAL

NS

+0

0.0 0

0.250000) 0

•

MAT

FORCECOF4STANTS

UPDATED

1.0102013—10 SYMMETRYCOOROINATE 5.7734537E—0I —8.9045773E 56 —2.4999762E—01 —3.4412899E—08 —2,23640675 15 —4.330)6605—01 0.0 6.5752845E—8O 0.0 5.77352645—01 2.125OCOOE 00 —6.5863442E—06 0.0 9.52295795 27 0.4336647E—06 0.0 3.1389086E—43 0.0 5.77352765—01 6.9233O71E—41 2.4950452E—OI 0.0 5.20323565 27 4.3301332E—01 0.0 9.92295760 27 0.0

0.0

O 0.0 5

0.0 0.0 0.0 0.0 0.0 0.0 C 0.0 C 0.0 5 0.0 0.0 I 0.0 3 0.0 0.0 I 0.0 5 0.0 0.0

ELEMENTS

MATRIX

1.30005—04

OF

(REPRESENTED

IN RECIPROCAL

SPACE)

FORCECONSTANTS

I.00000E

00

A(TRANSPOSE)

0.0 3.34503695—02 —1.4433604P—OI 3.0 5.5645172E—80 —2.5000203E—O1 0.0 1.57878450—80 0.0 0.0 2.3646801E 18 28867269E—01 0.0 —7.03057695 10 5.00000095—01 3.0 1.1441455E—72 0.0 0.0 6.39670545—79 —1.44342065—01 0.0 4.07)2)635 27 —2.4999897E—01 0.0 —5.76460750 10 0.0 DYNAMICAL

WAVE VECTOR.Q 0 0 0:0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0 • 0.0

:

3 -3 9 32 12 6

MASS

WAVEV5CTOR,Q=( 0.0 NOCH I NUN(T4 1

NON ZERO

MATRIX,

PARAMETERS:

:

• • • • • • • • • • • , • • • •

3.96695O5E 28 0.0 2.52641475 08 3.2345867E 5* 0.0 —9.22046675 23 —2.37132905 22 5.77343345—01 —9.6O54I13F—50 —1.93584115 20 0.0 —8.9045773E 56 —4.14763605 22 0.0 —22364067E 15 —2.37225135 22 5.77355385—01 6.61646465—80 —7.09392613E 22 0.0 0.0 —7.59396205 22 0.0 5.30651420—50 —1. 1496011E 23 5.77351875—01 2.613384050—80

—7.59645695—45 0.0 1.30656885—70 1.76324155—38 0.0 —I .4736889E—58 —4.6)446025 22 3.0 2.0673013E—49 4.5356159E—80 0.0 3,34503655—02 2.52810005—61 0.0 6.59052095—80 1.34277345—03 0.0 6.62078865—80 —1.33702005 23 0.0 0.0 —1.33801085 23 0.0 —5.31)86755—80 —3.788046)0 22 0.0 0.57445715—82

6.8B25329E—08 3.95166175—42 4.3301654E—O1 5.7734525E—01 —8.36768815 18 —2.4999762E—01 0.0 2.1250000E 00 0.0 0.0 9.9229198E 27 —1.4385041E—06 5.77352645—01 2.36452175 08 —6.589*3815—06 0.0 2.12500005 00 0.0 0.0 1.144I455E—72 —4.33013325—01 5.7735276E—0) 6.5811524E—R0 2.4999452E—OI 0.0 4.97)21635 27 0.0

0.0 2.9997163E—60 2.5000197E—01 0.0 9.92290095 27 —1.44336105—01 0.0 2.3645261E 16 0.0 0.0 2.12500000 00 5.0000000FOI 0.0 9.9229387F 27 2,8067269E—0I 0.0 2.36453055 1. 0.0 0.0 6.9233077E—4I 2.49998875—01 0.0 5.2032356F 27 —1.4434206E01 0.0 9.9228631E 27 0.0

I 1 I 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9

1 4 7 1 4 7 1 4 7 1 4 7 I 4 7 I

4

7 1 4 7 I 4 7 1 4 7

MATRIX:

2 5005—01 2:5005—0) 2.5000—01* 2.5000—0* 2.5005—01 2.5000—01 2.5005—01 2.5000—01 2.5005—01 2.5005—031 2.500E—OI 2,5005—Cl 2.500E—01 2.5005—01 2.SOOE—0) 2.5005—01 2.5000—0) 0_OI 2.500

) 1 1 1 1 I 3 I I I I 3 I 1 ) 3 3

FORCE CONSTANT 1.01320—10 1.01020—10 1.01020—10 I.OIOZF—l0 1.0102E—I0 1.01025—13 1.01025—10 1.01020—13 1.01020—10 1.51020—10 I .OIOSE—I0 1.31020—10 1.0102E—10 1.01020—10 1.0)020—13 1.01020—10 1.0102E—I0 1.01020—10

COMPLOX —4.49565—02 —2.79400—09 —6.71I2E—I0 —3.I453E—03 —1.98390—03 4.1O09E~’O.3 —2.C783E—02 —1.8626E—00 1.06540—09 —I.2582E—02 —19356E—03 1.6404E—02 —2.3413E—03 —2.67765—09 5.34910—ID —3.1453E—03 —1.98380—03 4.1(9095—03

ELERFNT 0.0 —I.3643E—03 I.1924E—04 0.0 —2.79630—13 0.01 0.0 I.6895E—0A 4.26065—03 0.0 1.001 30—09 0.0 0.0 —3.25205—03 —4.38880—03 0.0 —8.62120—IS 0.0

ROW I I 1 2 2 3 4 4 4 5 5 6

COLUMN I 2 3 2 3 3 4 5 6 5 6 6

7

7

7 7 8 4 8

8 9 8 9 9

3 3 3 3 3 3 3

3 3 3 3 3 3 3 3 3

3 3

3 3 3 3 3 3 3 3 3

AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT AT