Harold T. Stokes,
Department of Physics and Astronomy, Brigham
Young University, Provo, Utah, 84602, USA, branton_campbell@byu.edu

Larry L. Boyer
Complex Systems Theory Branch, Naval Research Laboratory, Washington, D.C.,
boyer@dave.nrl.navy.mil

This program calculates phonon frequencies and eigenmodes in a
crystal, using the method of frozen phonons. The calculations are
done in two steps. First the atomic displacements in each frozen
phonon are calculated by frozsl_init, taking full advantage of the
symmetry of the crystal. The user must then calculate the energy of
each of these frozen phonons, using whatever program he wishes. Then,
with these energies as input, frozsl calculates the phonon frequencies
and eigenmodes. For more details about frozen phonons and symmetry,
see H. T. Stokes, *Ferroelectrics* **164**, 183--188
(1995).

**Space-group settings.** Frozsl uses the space-group settings in
*International Tables of Crystallography*. For monoclinic
crystals, we use the setting with unique axis *b*, cell choice 1.
For trigonal crystals, we use the setting with hexagonal axes. For
crystals with two different origin choices, we use the setting with
origin choice 2.

**Symbols for k points and irreducible representations (irreps).**
Frozsl uses the convention of Miller and Love, *Tables of Irreducible
Representations of Space Groups and Co-Representations of Magnetic
Space Groups.*

(1) title line. This line is copied to the output, but is otherwise ignored by the program.

(2) space group symmetry of the crystal (1-230).

(3) lattice parameters of the crystal: a,b,c,alpha,beta,gamma, where a,b,c are the lengths of the lattice vectors defining the conventional unit cell, alpha is the angle between b and c, beta is the angle between a and c, and gamma is the angle between a and b. Distances are in bohrs, and angles are in degrees.

(4) number of displacement amplitudes for each frozen phonon. One is usually sufficient.

(5) displacement amplitudes in bohrs.

(6) number of terms in polynomial to which you wish to fit the energies for the different amplitudes. One is usually sufficient.

(7) power of each term in the polynomial. The power is 2 in the harmonic approximation.

(8) number of Wyckoff positions occupied by atoms in the crystal.

(9) for each Wyckoff position (each on a separate line),

(a) the
symbol for the atom

(b) the mass of the atom in atomic mass
units,

(c) the symbol for the Wyckoff position

(d) the values
for x,y,z if needed (These are *not* the atomic coordinates. These
are the dimensionless structural parameters for the Wyckoff position,
as defined in *International Tables for Crystallography*.)

(10) number of k vectors to be considered.

(11) for each k vector (each on a separate line),

(a) the symbol
of the k vector, using the notation of Miller and Love.

(b)
parameters *a,b,c*, if needed. For example, if the k vector lies
on a line of symmetry in the first Brillouin zone, then a value for
*a* must be given to specify the position on the line.
Similarly, values for *a,b* must be given to specify the position
of the k vector in a plane of symmetry, and values for *a,b,c*
must be given to specify the position of the k vector at a general
point.

(c) The LO mode at k=0 will be treated if in addition the
direction of the k vector is given (see below for more details).

(12) number of specified irreps, on a separate line following (11) above. If zero is entered, then every irrep for that k vector is considered.

(13) for each specified irrep (each on a separate line),

(a) the symbol for the irrep, using the notation of Miller and Love,

(b) the order parameter (optional), giving each of the n components, where
n is the dimension of the irrep. If the order parameter is not given,
the program finds the optimum choice.

(1) title, as read from the input.

(2) space group symmetry of the crystal, as read from the input.

(3) lattice parameters, as read from the input.

(4) number of displacement amplitudes for each frozen phonon, as read from input.

(5) displacement amplitudes, as read from the input.

(6) number of terms in fitting polynomial, as read from the input.

(7) power of each term in the polynomial, as read from the input.

(8) number of Wyckoff positions, as read from the input.

(9) each Wyckoff position, along with its symbol and atomic mass, and, if given as input, the effective charge tensor for the first atom in each Wyckoff position.

(10) number of k vectors to be considered, as read from input.

For each k vector:

(11) k vector symbol and coordinates. These are given in terms of the lattice vectors of the conventional reciprocal lattice.

(12) irreps for each k vector.

For each irrep:

(13) order parameter used by frozsl. This defines a subgroup for which the space group symmetry, lattice parameters, and Wyckoff positions of the atoms are also given.

(14) for irreps at k=0, an indication of whether the modes are longitudinal or transverse. If longitudinal, the direction of the k vector is given.

(15) for irreps at k=0, an indication of whether acoustic modes are present.

(16) vectors defining the superlattice of the phonon modes. Components are given in both dimensionless units and in cartesian coordinates (bohrs).

(17) number of atoms in the unit cell of the superlattice.

(18) type of each atom. These are numbered according to the Wyckoff positions read from the input. Atoms in the first Wyckoff position are type 1, atoms in the second Wyckoff position are type 2, etc.

(19) position of each atom. Components are given in both dimensionless units and in cartesian coordinates (bohrs).

(20) for longitudinal modes at k=0, the Born effective charge tensor for each atom.

(21) number of symmetry vectors. Each phonon mode is a linear combination of these symmetry vectors.

(22) atomic displacements in each symmetry vector and the mass of the atoms displaced. Components are given in (1) cartesian coordinates (bohrs) and (2) structural parameters for the Wyckoff positions in the subgroup.

(23) number of frozen phonon modes

For each frozen mode:

(24) coefficients in the linear combination of symmetry vectors in frozen mode.

(25) for each amplitude, the value of the amplitude and atomic positions for the frozen mode with the amplitude. These are the positions you should use when calculating the energy of the frozen phonon. Atomic positions are given in (1) dimensionless coordinates, (2) Cartesian coordinates (bohrs), and (3) structural parameters for the Wyckoff positions in the subgroup.

Besides the input to frozsl_init as described above, an energy of each frozen phonon must be entered in the box provided. The units must be in millihartrees (1 hartree = 2 rydbergs = 27.196 eV). If calculating longitudinal modes at k=0, then the high-frequency dielectric tensor and the Born effective charge tensor for the first atom of each Wyckoff position must also be entered in the box provided. (See below for more details.)

(1) title, as read from the input.

(2) each irrep symbol.

(3) for irreps at k=0, an indication of whether the mode is longitudinal or transverse.

(4) for each amplitude of each frozen phonon mode, the value of the amplitude and the energy.

(5) each eigenmode obtained from the diagonalization of the dynamical matrix, along with the frequency and eigenvector. The eigenvectors have been adjusted so that they are coefficients in the linear combination of the symmetry vectors listed in the output of frozsl_init. Phonon modes which are unstable have imaginary frequencies.

The accuracy of the frequencies and eigenvectors can be improved by iterating the calculation. An iteration of frozsl_init takes the eigenvectors from the output of frozsl, and uses them in calculating new frozen phonons. The energy of each of these frozen phonons must be calculated and entered in the box, following the energies for the old frozen phonons. The iteration of frozsl then calculates the eigenmode frequencies and eigenvectors from these new frozen mode energies.

By default, frozsl only treats the transverse optical (TO) phonons at
k=0. The longitudinal optical (LO) modes at k=0 require special
treatment because of the electric field associated with
long-wavelength longitudinal modes. To specify the calculation of LO
modes, enter the direction of the wave following the k-vector symbol GM
on input line (11) to frozsl_init. The direction is specified in
cartesian coordinates using the same coordinate system that
frozsl_init uses. For example, GM 1 0 0 would specify an LO mode in
the *x* direction.

The calculation of the LO mode requires the high-frequency dielectric tensor for the crystal and the Born effective charge tensor for each atom. In the box provided, enter first the high-frequency dielectric tensor and then, for the first atom of each Wyckoff position, the Born effective charge tensor. The tensors should be given in the Cartesian coordinate system used by frozsl_init. Each of the 9 elements of each tensor should be entered consecutively. If the dielectric tensor is not known, you may enter a unit tensor: 1 0 0 0 1 0 0 0 1. If the Born effective charge tensor is not known, you may enter a unit tensor multiplied by the charge. For example, for an atom of charge 2, enter 2 0 0 0 2 0 0 0 2.