Introduction

In this project we calculate the vibrational normal modes of ammonia in the harmonic approximation. We performed these calculations for all electrons within the framework of DMol³⁽¹׳²⁾ using a basis set of atomic orbitals and the PBE⁽³⁾ functional for the exchange-correlation energy. The specific basis set used was DNP 3.5 (double numerical with polarization), and the convergence required of the self-consistent fields was 10⁻⁶.

Figure 1: The ammonia molecule at equilibrium, shown here, was calculated to have NH bond lengths and HNH angles respectively of 1.022Å and 106°.

Methods

In the harmonic approximation, one calculates the second derivatives of atomic energy with respect to all degrees of freedom, generating a matrix of second derivatives known as the Hessian. Weighting the entries of the Hessian according to the masses of the atoms under consideration, one obtains the mass-weighted Hessian, whose eigenvalues (in Cartesian coordinates) are the squares of the frequencies of the normal modes. A four-atom system such as ammonia in three-dimensional space has 12 degrees of freedom and hence 12 normal modes; however, three modes are zero-frequency whole-molecule translations and three are whole-molecule rotations. These are discarded.

To better understand what excitations or displacements qualify as “small” and to explore the stability of numerical differentiation for a range of step sizes, we calculated the second derivative associated with stretching an NH bond for varying finite difference steps and plotted the resulting parabolae against the DFT-calculated energies in Figure 2.

Figure 2: The energy associated with a displacement from equilibrium by 0.005Å, 0.01Å, 0.02Å, and 0.04Å was used to create the approximate potentials labeled in the plot legend by the step sizes used to calculate the effective spring constant. Notice that the effective spring constants all agree closely, and for displacement not exceeding 0.02Å agree with the DFT-calculated energies within 3% or less.

However, it would be naive to assume that the same step size is equally valid for all calculations, a point illustrated in Table 1.

Table 1: Numerical second derivatives in units of (kcal/mol)/Ų with fixed step size 0.005Å are compared here to the derivatives as calculated by Materials Studio. By “H1xx”, for example, we mean the diagonal matrix element for which both differentiations in the second derivative are taken in the x direction for atom H1. The agreement is generally good to within a few percent except for the some of the smaller results, which differ by about 15%. Moreover, notice the difference between the second derivatives in z; by symmetry these should be equal but in the Hessian calculated by Materials Studio they are not. Differentiation by hand, however, reproduced that symmetry.

As we mentioned before, the Hessian is calculated in Cartesian coordinates so that the squared normal mode frequencies will be the eigenvalues of the Hessian. In more general coordinate systems, one would have a more general matrix equation for the normal mode frequencies. Though a general matrix equation would be more difficult to solve, with judiciously-chosen coordinates we could remove the possibility of zero modes. In our example of ammonia, a natural set of coordinates would be stretching and polar bending of the three NH bonds. These generalized coordinates are sufficient to describe any relative orientation of the four atoms without allowing for simultaneous translation or rotation. Having not only six degrees of freedom but also high symmetry (as there would be no “preferred” Cartesian axes, only identical bond lengths and angles) would also greatly reduce the number of finite difference calculations necessary to generate a complete Hessian. However, for the sake of generalizability to arbitrary molecular geometries, we perform these calculations in Cartesian coordinates. A matrix representation of the translations and rotations is used to project those modes (i.e. the zero modes) out from the results.

We may still use some symmetry arguments to simplify our work. For example, suppose we choose to place N at the origin and one of the Hs (say H1) in the xz-plane. Then the other two H atoms (H2 and H3) are mirror reflections of one another, and any calculation we do for H2 can be reused for H3. Likewise, any positive y-displacement for H1 by symmetry has the same effect as an equal displacement in the opposite direction. Finally, if we choose to orient the atom’s dipole axis with the z-axis, displacing any H atom in the z-direction would have the same effect. Nonetheless, the number of calculations required to do this by hand is infeasible, and hence we use a Hessian calculated automatically by Materials Studio.

Results

Table 2: These frequencies, reported in terms of wavenumber, were obtained in this work through Materials Studio.

Figure 3: The six diagrams above corerspond to the six normal modes of ammonia described in Table 2 in order of ascending frequency or wavenumber. Where possible (i.e. diagrams 1 through 3) arrow size is used to represent the relative magnitude of the displacements; elsewhere (i.e. diagrams 4 through 6) only direction is indicated. In all diagrams, the relative directions indicate phase. For example, in diagram 4 all bonds are compressed (and later stretched) simultaneously, whereas in diagram 5 the upper-left bond stretches while the lower one compresses, and vice versa.

References

[1] B. Delley, “An all‐electron numerical method for solving the local density functional for polyatomic molecules”.J. Chem. Phys. 92, 508 (1990).
[2] B. Delley, “From molecules to solids with the DMol³ approach”. J. Chem. Phys. 113, 7756 (2000).
[3] J. P. Perdew, K. Burke and M. Ernzerhof, “Generalized gradient approximation made simple”. Phys. Rev. Lett. 77, 3865 (1996).