CASTEP is primarily a solid state code with periodic boundary conditions, and is not necessarily the first choice for performing vibrational spectroscopy calculations on molecules. Nevertheless it can sometimes be convenient to use it in this mode (for example, if it is desired to compare a molecular crystal with an isolated molecule at exactly the same level of theory).
With certain limitations this can be done in CASTEP using the “molecule in large box” method. The idea is to place the molecule in a vacuum by constructing a large supercell. If the supercell is sufficiently large that the periodic images of the molecule do not interact then the frequencies and eigenvectors in the molecular limit can be recovered. Because the zero charge density in a large volume must still be represented using a plane-wave basis sat, such calculations can become expensive, but it is nevertheless essential to perform careful cell size convergence tests.
A molecule in free space possesses rotational as well as translational symmetry, which results in three (or two for a linear molecule) free librational modes of zero frequency. A periodic CASTEP calculation does not possess this symmetry and the librational mode frequencies will have nonzero values, either real and positive or imaginary. If a molecular electric field calculation is desired the parameter
may be specified with the keyword value of MOLECULE or LINEAR_MOLECULE. This omits the 6 (or 5) lowest modes from the infra-red polarizability calculation, assuming them to be the librations.
As with crystalline systems it is desirable to maximise the use of symmetry by optimally orienting the molecule. In this case an additional criterion applies - the simulation cell vectors should be organised for maximal compatibility with the molecular symmetry. Any incompatible symmetry operations will not be found if using the SYMMETRY_GENERATE keyword. For example a tetrahedral molecule is best modelled in a cubic cell with the 3-fold axes oriented along the cube diagonals. By contrast a hexagonal molecule such as benzene is better modelled using a hexagonal simulation cell with the molecule commensurately oriented.
A limitation of this method is that CASTEP implements only the crystallographic point group symmetry operations, and does not include all of the molecular point groups. It is not always possible to take full advantage of molecular symmetry, and 5-fold rotation axes and icosahedral groups can not be represented. In such cases some degeneracies will be lifted and only approximately satisfied, and the group theoretical analysis of the eigenvectors will not be correct for the molecular point group.
Finally, it is sometimes desirable to compute the vibrational spectrum of an ion rather than a neutral molecule. This adds a number of complications of different degrees of severity. First is the well-known result that semi-local DFT (LDA and GGA) severely underbinds anions, and the description of the molecular states in some systems may be substantially incorrect. Second, the introduction of a non-zero molecular charge in a periodic cell necessitates the addition of a compensating charge to avoid a divergent Coulomb energy. CASTEP implicitly adds a uniform charge distribution which integrates over the cell to the negative of the sum of ionic and electronic charge in the system. This model while effective, gives rise to additional electrostatic terms which decrease only as 1∕L3. Consequently it is impossible to completely converge such a calculation with respect to cell size. It is beyond the scope of this document to describe the techniques which may be used to recover the “infinite” volume limit. The reader is referred to ref [17] and the extensive literature which cites this paper for further reading. An example of a successful approach involving an extrapolation of frequencies of a molecular anion to the infinite cell limit is given in Parker et al. [18].
A powerful technique in experimental spectroscopy is isotopic substitution where the sample is modified by the substitution of a common isotope of some element by a heavier or lighter one. This can provide valuable site-specific spectroscopic information, particularly if the substitution can be performed in a site-specific manner. Force constants and dynamical matrices do not depend on nuclear mass, and are therefore one set of computations will suffice irrespective of isotopic substitution. The phonons tool (section 3.4) allows recomputation of the frequencies for different isotopic substitutions without recomputing dynamical matrices.
The simplest case is complete substitution. The new .cell file should contain a block specifying the new mass, e.g.
Running the phonon code will compute the vibrational frequencies and eigenvectors for a completely deuterated sample.
Slightly more complicated is the case where only one or a few of a number of sites is to be substituted. In that case the above prescription is inadequate as it would change all of the sites containing the species in question. It is therefore necessary to use CASTEP’s sub-species labelling capabilities. For example the original cell containing a methane molecule might contain
The continuation cell for input to phonons should contain
The new atom label may contain any alphanumeric extension following the colon up to a maximum of 8 characters.
This defines a new system with one hydrogen atom of a methane molecule replaced by deuterium. Notice that this has broken the initial symmetry of the cell. The phonons program generates a new system of reduced symmetry with three distinct atom types, and copies the dynamical matrix data from the original methane system from the .check file over to the corresponding atoms in the new system. It then diagonalises the dynamical matrix, applies whatever post-processing is specified and writes the frequencies and eigenvectors as usual.
This straightforward approach fails to take into account isotopic shifts in bond lengths and geometry, and is therefore an approximate one. Since isotopic shifts in bond lengths depend in turn on the vibrational frequency, within the quasiharmionic approximation a self-consistent approach under which bond lengths and frequencies are simultaneously adjusted to self-consistency would be required. This would also require new DFPT electronic calculations at each stage and is beyond the scope of this utility.
It is sometimes desirable to compute the vibrational frequencies of a restricted region of an ab initio model. Consider the case of a molecule adsorbed on a surface. If only the frequencies of the molecular vibrations are needed it would be desirable to compute a subset of the modes were this possible. While it is not practical (as it would require prior knowledge of the eigenvectors of the full calculation), an alternative approach is to apply constraints to certain atoms.
CASTEP implements a technique known as constrained lattice dynamics, also known as the partial hessian method. Nominated atoms are assumed to be “frozen”, and the corresponding entries of the dynamical matrix are set to zero. The model is effectively one whereby the atoms in the region of the system deemed “irrelevant” are assigned a mass of infinity. It is not necessary to perform any computations for perturbation of these atoms, and a considerable saving of computational effort may be achieved.
Frozen atoms are specified a CASTEP .cell file using the same syntax as applies to geometry operations. The block
constrains silicon atoms numbered 1-2 and C atom number 4 to be fixed. These do not move during geometry optimisation, and their perturbations are not considered during a lattice dynamics calculation. In fact the constrained lattice dynamics method does not make full use of the generality of CASTEP’s linear constraints block (see for example the tutorial on MD at http://www.castep.org) but only identifies atoms which are fully constrained not to move. As in the example above, there should be a single line for each atom which creates a uniquely numbered constraint. This should contain a “1” in all of the x, y, z positions.
Except in rather specialised geometries the presence of fixed atom constraints is incompatible with most symmetry operations, and therefore symmetry should usually be turned off during a constrained lattice dynamics calculation. There is also an incompatibility with the acoustic sum rule (section 2.3.4) as constraining the atoms breaks the translational invariance of the Hamiltonian. Acoustic sum rule correction is therefore disabled automatically if constraints are present.