Using forces to accelerate first-principles anharmonic vibrational calculations

High-level vibrational calculations have been used to investigate anharmonicity in a wide variety of materials using density functional theory methods. We have developed an efficient approach for describing strongly anharmonic systems using a vibrational self-consistent field method. By far the most computationally expensive part of the calculations is the mapping of an accurate Born-Oppenheimer (BO) energy surface within the region of interest. Here we present an improved method which reduces the computational cost of the mapping. In this approach we use data from a set of energy calculations for different vibrational distortions of the materials and the corresponding forces on the atoms. Results using both energies and forces are presented for the test cases of the hydrogen molecule, solid hydrogen under high pressure including mapping of two-dimensional subspaces of the BO surface, and the bcc phases of the metals Li and Zr. The use of force data speeds up the anharmonic calculations by up to 40%.

Figure 1

Convergence of anharmonic correction to the energy at 0 K for ${\mathrm{H}}_2,\mathrm{\Delta }{E}_{\text{anh}}=E_{\text{anh}}-E_{\text{har}}$, with respect to the number of mapping points used per mapping direction, for the basic VSCF method, as well as the VSCF+f method, fitting with both cubic and quintic splines. The inset shows the results at low numbers of mapping points on a different energy scale.

Figure 2

Views of the structures of solid hydrogen considered at 100 GPa. Both of the $Cmca$ structures are shown looking along the $x$ axis, while the $C2/c$-24 structure is shown looking along the $y$ axis. As these are layered structures, atoms in inequivalent layers are denoted by different colors. Green, purple, translucent green, and translucent purple denote the first, second, third, and fourth inequivalent layers, respectively.

Figure 3

Convergence of the anharmonic correction to the energy at 0 K, $\mathrm{\Delta }{E}_{\text{anh}}$, of the Cmca-4, Cmca-12, and $C2/c$-24 phases of solid hydrogen at 100 GPa, with respect to the number of mapping points used per mapping direction. The convergence of the basic VSCF method as well as the VSCF+f method of fitting with a quintic spline are shown. The inset shows the results at low numbers of mapping points on a different energy scale.

Figure 4

Results of anharmonic vibrational calculations for the Cmca-4 structure of solid hydrogen including selected 2-D subspaces of the BO surface. The left-hand column shows BO surfaces mapped in the labeled subspace, where “puc” stands for “per unit cell,” while the right-hand column shows the convergence of the correction to the vibrational energy due to the relevant 2-D term with respect to the number of mapping points used. (a) Subspace corresponding to directions 4 and 5. (b) Subspace corresponding to directions 4 and 7.

Figure 5

Convergence with respect to the number of mapping points used per mapping direction of contribution of 1-D terms corresponding to the mapping directions defined by soft modes to the total anharmonic vibrational energy of the bcc phase of zirconium at 0 K. The convergence of the basic VSCF method as well as that of the VSCF+f method using both a cubic spline and a quintic spline are shown.

Figure 6

The shape of the BO surface as mapped along one of the soft modes of zirconium; “puc” stands for “per unit cell.”

Figure 7

Convergence of the temperature at which the bcc phase of zirconium becomes dynamically stable with respect to the number of mapping points used per mapping direction. Results using the basic VSCF method and the VSCF+f method using both a cubic spline and a quintic spline are shown. The VSCF+f cubic spline results are hidden behind the VSCF+f quintic spline results.