4.6 Computation of vibrational and phonon properties

There is an older perl script and source codes for this step for the entire finite-difference calculation included with FHI-aims. To compile the tools for this, go to the build directory in which you compiled FHI-aims. From there, you can type either make vibrations, make vibrations.mpi, or make vibrations.scalapack.mpi.

This will create two files in the src/vibrations subdirectory of the FHI-aims binary directory. One is the actual diagonalization subprogram, aims.vibrations.*.x the other one the perl-wrapper numerical_vibrations.pl. This location of the Perl script (which is in the same folder as the diagonalization subprogram) should be displayed after the make command.

You will need to edit the header of the numerical_vibrations.pl script and give three key parameters: The absolute location of the FHI-aims binary directory ($AIMS_BINDIR), the aims executable to be used ($EXE), and the calling command for the executable ($CALLING_PREFIX and, if needed, $CALLING_SUFFIX).

The $CALLING_PREFIX is required for specification of parallel runs, for example

$CALLING_PREFIX = mpirun -np 2

for a two-core parallel run. In addition, some parallel environments (most notably IBM’s poe require that the number of cores is specified after the executable is called, for those instances you can set $CALLING_PREFIX. Note that the postprocessing routines for the vibrations are also parallel. Please also check that $HESSIAN_DIAGONALIZATION_EXE has been automatically set correctly.

Having compiled the vibrations script, running it is very simple. We use the same water molecule that was relaxed in the example section 4.1. This test case is contained in the directory testcases/H2O-vibrations.

We are here calculating a second derivative of the energy from the numerical change on the forces arising from small finite displacements $\delta$, in the following way:

Here, $E$ is the total energy as a function of all atomic coordinates $x_i$ ($i$=1,…,3$M$ for $M$ individual atoms), and ${𝐅}_i$ are the components of the analytical first derivatives (forces) on each atom (again, $i$=1,…,3$M$) due to a displacement of coordinate another $x_j$ by an amount plus or minus $\delta$.

Since we are interested in the Hessian in a local minimum of the potential energy surface (ideally, ${𝐅}_i(\{x_j\})$=0 for all $i$ at the local minimum geometry $\{x_j\}$), this expression is the difference of two force values that are themselves already near zero. The smaller the displacement $\delta$, the smaller the absolute magnitude of the forces to be subtracted, giving greater weight to any residual numerical noise in the calculation. The larger the displacement $\delta$, the larger does the difference become, but at the same time, we will also move out of the region of the potential energy surface that is exactly harmonic, and introduce systematic errors into the Hessian. We are thus faced with a tradeoff between choosing a value $\delta$ that is large enough to be free of any potential numerical noise, yet small enough to avoid large systematic errors due to real anharmonicities. The default value for $\delta$ is 0.0025 Å, but this value should be checked explicitly in any practical calculation.

Since the vibrational frequencies depend strongly on the accuracy of the forces, and on any near-zero residual forces on the fully relaxed geometry itself, we apply the following changes to the earlier file control.in:

• •

  We use the tight species defaults for elements H and O. Among other things, these defaults prescribe denser integration grids than light and thus lead to more accurate forces. Obviously, the previous optimum geometry from the light settings should be postrelaxed to the exact minimum for the modified tight settings, and our script does this automatically (but check the output just in case).

• •

  We increase the accuracy to which the forces are converged in the self-consistency cycle, sc_accuracy_forces, to 1E-5. However, be sure to check that the remaining settings, especially sc_accuracy_eev are accurate enough. You do not want your FHI-aims calculation to spend more than 1-2 s.c.f. cycles per geometry on actually converging the forces, because each s.c.f. iteration with forces can be up to a factor 10 more expensive than an s.c.f. iteration during which the forces are not yet checked.

• •

  We change the relax_geometry bfgs convergence criterion to 1E-4. This is, again, a rather harsh requirement, and one should verify (especially for large molecules) that the relaxation algorithm is not spending large amounts of time on failed relaxation steps near the minimum, probing any residual numerical noise of integration grids etc. rather than following an actual, smooth potential energy surface.

These new settings contain all that is necessary to get close enough to the actual minimum. Again, the convergence settings quoted above should normally be fine, but in case of doubt, check. If, for some reason, those stringent criteria can not be reached exactly, it is still preferable to obtain a slightly less converged Hessian within an amount of CPU time that is actually finite.

In order to run the calculation of the vibrations, change to the directory containing the input files control.in and geometry.in and run the numerical_vibrations.pl script in the directory. That’s it. The calculation should take no more than a few minutes to complete. It first relaxes the structure to its minimum configuration and then applies six finite displacements to each of the atoms: 18 single-point calculations in total for the water molecule. The result (using the default settings for $\delta$) can be visited in subdirectory test_default.

The output stream contains all of the important data, which is additionally saved in a number of temporary and output files. The file basic.vib.out contains the output of the normal mode analysis, while the file basic.xyz contains the actual eigenmodes (vibrations), which can be read by a number of molecular viewers.

A short note on the actual physical output basic.vib.out: The frequencies are given in cm¹, the corresponding zero point energies in eV, and their cumulative total. It is important to ensure that the first six eigenmodes are close to zero, which means that the structure in question actually corresponds to a minimum on the potential energy surface. FHI-aims does not do an a priori reduction of the translational and rotational eigenmodes, which allows an explicit check on the quality of the strucure.

In the default version (no command line options specified), the vibration script automatically assumes a jobname basic. This name, as well as the finite difference displacement $\delta$, can be changed on the command line, leading to the following complete calling form of the script:

numerical_vibrations.pl $⟨$jobname$⟩$ delta

$⟨$jobname$⟩$ is a name of choice that will be appended to all output files produced during the run. One can thus run the same job with different settings $\delta$ multiple times in the same directory, starting from the same input files.

The second parameter, delta, is the finite displacement $\delta$ used for calculating the Hessian matrix. Its default is 0.0025 Å, which has been doing a good job for all the cases we are aware of, but still, please verify that it works for your particular system. (See the brief discussion above).

For example, for our test case, subdirectory test_delta_0.001 contains a second example run with the following command line:

numerical_vibrations.pl test_0.001 0.001

For the default settings, our resulting frequencies look like this:

  Mode number      Frequency [cm^(-1)]   Zero point energy [eV]   IR-intensity [D^2/Ang^2]
            1             -12.66044778              -0.00078485                 1.89406419
            2              -0.15559990              -0.00000965                 0.00266735
            3              -0.00081379              -0.00000005                 0.00000000
            4               0.04885082               0.00000303                 0.00002821
            5               6.58934967               0.00040849                 0.00000000
            6               7.01167236               0.00043467                 5.41560478
            7            1592.96295435               0.09875111                 1.63341394
            8            3708.86739272               0.22992045                 0.04298854
            9            3813.74446188               0.23642200                 1.07616306

In contrast, the smaller value $\delta$=0.001 Å produces the following results:

  Mode number      Frequency [cm^(-1)]   Zero point energy [eV]   IR-intensity [D^2/Ang^2]
            1              -5.13719038              -0.00031847                 1.89407621
            2              -0.06616895              -0.00000410                 0.00323346
            3              -0.00073155              -0.00000005                 0.00000000
            4               0.01938144               0.00000120                 0.00002697
            5               2.47172391               0.00015323                 0.00000000
            6               2.70760068               0.00016785                 5.41498492
            7            1593.00852915               0.09875393                 1.63340808
            8            3708.82166272               0.22991761                 0.04299176
            9            3813.70936996               0.23641982                 1.07615430

So, in this case, the rather tight relaxation and force convergence settings do produce a set of translations and rotations that are closer to zero for $\delta$=0.001 Å than the default.

However, once again be warned that for larger molecules the same harsh convergence criteria may not be applicable, and a too small value $\delta$ can in fact inflate any residual noise. So, a small value $\delta$=0.001 Å or even less should not be blindly expected to improve numerical accuracy.

It should also be noted that the “physical” vibrational frequencies 7-9, above the six translations and rotations, remain completely unimpressed by the change of $\delta$, either way. This is generally true, and would even remain true for a yet larger value $\delta$=0.005 Å, which would also be a reasonable default choice.

Finally, a word about restarting the numerical_vibrations.pl script if it performed a part of the needed single-point force calculations but did not finish all of them, typically as a result of queue time limits on a shared computer. In short, one can simply restart the script in the same directory as before right away. The script will figure out how far it got during its last run and it will then pick up from there to continue and eventually finish the calculation of vibrational frequencies and free-energy terms.

Raman spectra: A third optional argument is available in the vibrations script, which allows for the calculation of harmonic Raman spectra. In order to obtain such spectra, one should type the following,

numerical_vibrations.pl $⟨$jobname$⟩$ delta polar  .

The polarizabilities are obtained with DFPT, and the derivatives with respect to the displacement are then calculated in a similar way as for IR spectra. All the formulae can be found in the following reference: Neugebauer et al., J. Comput. Chem. 23: 895-910, 2002. In particular, Eq. (43) is used to produce the final Raman intensities. Keep in mind that Raman spectra obtained in such fashion are much more time-consuming than IR spectra, the latter requiring only DFT. Note that Raman spectra for solids can also be calculated in this manner (only $\mathrm{\Gamma }$-point is involved). The Raman intensities that are outputted from this script correspond to combinations of the different tensor components combined to form what is commonly called “powder” spectrum.

In the following we provide the description of tags defining the output of free energies.

The calculation of the free energy of vibration and rotation (rigid rotor approximation) can be requested with this keyword. They are calculated between the temperatures Tstart (K) and Tend (K) with a total of Tpoints steps. The implemented equations read:

with N being the number of atoms in the molecule and ${\omega }_i$ depicting the vibrational frequencies.

where $I_A$, $I_B$, and $I_C$ depict the moments of inertia of the molecule.

The calculation of the free energy of translation can be requested with this keyword. It only works if vibrations free_energy is used as well. The free energy of translation is calculated within the temperature range given with vibrations free_energy and between pressures pstart (Pa) and pend (Pa) with ppoints steps. The equation reads:

with $p$ denoting the partial pressure.

A newer implementation of the same technique (finite-difference approach) as the Perl script can be achieved through a python script get$\mathrm{\_}$vibrations.py found in folder utilities, and through which the running mechanism of this script is described below. The script works for nonperiodic and periodic systems (only $\mathrm{\Gamma }$-point) as well.
python get_vibrations.py [options] <name> <mode>
The script functions in three modes:

• •

  mode 0: The script creates folders containing the aims input files for single calculations. These are the 6N calculations for the N unconstrained atoms in the geometry.in in the execution path. The atoms are displaced by +/-delta (default 0.0025Å) in x, y, and z. Moreover, the control.in is copied to all folders. The user must then explicitly run all these calculations.
  

• •

  mode 1: The script calculates normal modes based on FHI-aims output files, after calculations have been run following their creation with mode 0.
  

• •

  mode 2: The script calculates IR and/or Raman spectra based on FHI-aims output files, after calculations have been run following their creation with mode 0.
  

The available options with the discussed modes:

• •

  -d : The amount of displacement in AA along Cartesian coordinates(mode 1,2 and 3), default is 0.0025 AA.

• •

  -p : Plotting the spectra (Raman/IR) (mode 2 only).

• •

  -b : Broadening for IR and Raman spectrum in $c{m}^{-1}$, default is 5.

• •

  -w : If imaginary frequencies are present, then files of the system displaced along the geometries of the imaginary frequencies are outputted (mode 2 and 3).

• •

  -m : Printing the geometries, frequencies and vibrations in molden format.

• •

  -M : For calculating both Raman and infrared intensities (mode 2 only).

• •

  -R : For calculating Raman intensities only (mode 2 only).

• •

  -I : For calculating infrared intensities only (mode 2 only).

For further details about the calculation of the polarization of solids from the Berry phase formalism, please see Section 3.39, where the k-point grid mentioned above is better explained.

1. python get$\mathrm{\_}$vibrations.py name 0 : Creating all folders as discussed above.
2. python get$\mathrm{\_}$vibrations.py name 1 -w: Post-processing on the FHI-aims output files generated by running the calculations set up in example 1. The output files will be:

• •

  masses.name.dat : Atomic mass of atoms in a.m.u and their position as in geometry.in .

• •

  hessian.name.dat : Hessian matrix before mass weighting.

• •

  massweighted$\mathrm{\_}$Hessian.name.dat : Mass weighted hessian matrix (Dynamical matrix.)

• •

  mass$\mathrm{\_}$vec.name.dat : Diagonal matrix of 1/sqrt(mass)

• •

  eigen$\mathrm{\_}$vectors.name.dat: Normal modes written in a matrix form, each column corresponds to a normal mode arranged from lower to higher frequency.

• •

  car$\mathrm{\_}$eig$\mathrm{\_}$vec.name.dat : Eigen vectors that are not normalized (not mass weighted).

• •

  normalmodes.name.dat: Normal modes and their corresponding frequencies.

• •

  name.distorted.vibration : Geometry displaced along imaginary frequency mode (if exists).

3. python get$\mathrm{\_}$vibrations.py name 2 -R -p -m > output: Post-processing step to compute Raman spectra. The output files will be:

• •

  masses.name.dat: Atomic mass of atoms in a.m.u and their position as in geometry.in .

• •

  hessian.name.dat: Hessian matrix before mass weighting.

• •

  grad$\mathrm{\_}$dipole.name.dat: Derivative of dipole moment with respect to displacement delta.

• •

  name.Raman: Raman intensities and the corresponding frequencies.

• •

  name.molden: Molden geometries and displacements for the mode.

• •

  name.xyz: All vibrational displacement, this file can be opened using jmol to view the vibrations.

• •

  name$\mathrm{\_}$Raman$\mathrm{\_}$spectrum.pdf: Raman spectra plot

4. python get$\mathrm{\_}$vibrations.py name 2 -I -p > output:
Post-processing step to compute IR spectra. . The output files will be:

• •

  masses.name.dat: Atomic mass of atoms in a.m.u and their position as in geometry.in .

• •

  hessian.name.dat: Hessian matrix before mass weighting.

• •

  grad$\mathrm{\_}$dipole.name.dat: Derivative of dipole moment with respect to displacement delta.

• •

  name.ir: Infrared intensities and the corresponding frequencies.

• •

  name.molden: Molden geometries and displacements for the mode.

• •

  name.xyz: All vibrational displacement, this file can be opened using jmol to view the vibrations.

• •

  name$\mathrm{\_}$IR$\mathrm{\_}$spectrum.pdf: IR spectra plot

It is important to note that the atoms that will not enter in the vibrational analysis are followed with constrain$\mathrm{\_}$relaxation .true. in geometry.in file

The following lists the keywords to activate the calculation of the vibrational frequencies with DFPT.

Here is an example, the following need to be added to control.in:

DFPT vibration_reduce_memory
DFPT_mixing   0.2       #default is 0.2
DFPT_sc_accuracy_dm   1E-6  # default is 1.0d-6

Here is an example for using DFPT polarizability, the following need to be added to control.in:

DFPT polarizability
DFPT_mixing   0.5             #default is 0.2
DFPT_sc_accuracy_dm   1.0d-6  # default is 1.0d-3
dfpt_pulay_steps 6            # default is 8

Here is an example for using DFPT dielectric, the following need to be added to control.in:

ΨDFPT dielectric
ΨDFPT_mixing   0.5             #default is 0.2
ΨDFPT_sc_accuracy_dm   1.0d-6  # default is 1.0d-3

The recommended way to compute phonon properties with FHI-aims is through the Python package FHI-vibes [170]. FHI-vibes is a python package built on top of the Atomic Simulation Environment (ASE, [186]), and tightly integrates FHI-aims with phonopy [298]. With that, FHI-vibes provides a unified user interface to run and process the necessary calculations. It thus allows harmonic phonon calculations via phonopy [298], but also gives access to more advanced analysis and post-processing methods, especially for the study of anharmonic effects.

Details on how to obtain, set up, and use FHI-vibes –including tutorials on how to perform phonon calculations– can be found in the official documentation: https://vibes.fhi-berlin.mpg.de/.

Alternatively, it is also possible to use phonopy directly (without FHI-vibes) by using the interface to FHI-aims included in phonopy (Version 2.5 and higher). Documentation and examples can be found in the phonopy source code at https://github.com/phonopy/phonopy/tree/develop/example/diamond-FHI-aims (cf. README.md).

Please note: The former interface phonopy-FHI-aims is not officially supported anymore. The last FHI-aims release including phonopy-FHI-aims and respective documentation was the 210226 release.

Note that the former, native phonon infrastucture in FHI-aims is no longer needed. The previous make targets phonons and phonons.mpi still exist, but we no longer recommend them. FHI-vibes/Phonopy is, simply, more stable and general.

A third level of infrastructure for vibrations and phonons in non-periodic and periodic systems based on density functional perturbation theory (DFPT) is nearing readiness and usable to some degree. However, parts of this implementation are not yet optimized. The DFPT implementation should therefore still be considered experimental at this time. We recommend to cross-check the DFPT results always with finite-difference methods described above.

Here is an example for using DFPT phonon, the following need to be added to control.in:

DFPT phonon
DFPT_mixing   0.5             #default is 0.2
DFPT_sc_accuracy_dm   1.0d-6  # default is 1.0d-3
dfpt_pulay_steps 6            # default is 8

Here is an example for using DFPT phonon_reduce_memory, the following need to be added to control.in:

DFPT phonon_reduce_memory
DFPT_mixing   0.5             #default is 0.2
DFPT_sc_accuracy_dm   1.0d-6  # default is 1.0d-3
dfpt_pulay_steps 6            # default is 8