Category Archives: 3rd Post 2018

A brief introduction to the theory, parameterization and application of Projector Augmented Wave (PAW) datasets to plane wave DFT calculations

by Brandon Bocklund

Introduction

Many of the observable properties that makes materials interesting are linked to the electronic structure, so it is desireable with to perform DFT calculations with the most accurate computational methods that are affordable. Calculations within DFT scale as \( \mathcal{O}(n^3) \) for n electrons [1], so many technologically important materials, including transition metals such as Fe, are computationally demanding to calculate when considering all the electrons. This is exacerbated for core electrons, since the wave functions near the nuclei tend to oscillate, requring a large number of plane waves required to describe the wave functions near the nuclear cores in extended systems. Since many key properties are influenced by the outermost valance electrons that control bonding, it is sensible to fix the core electrons and treat them approximately to make calculations of the same level of theory more affordable or to make more complex calulations affordable. Treating core electrons as fixed is known as the frozen core approximation.

In this post, different cutoff radii will be used to generate PAW datasets for Ti. Differences in the generated PAW datasets will be demonstrated by comparing the cutoff energy required to converge the energies of hcp Ti and by comparing the formation energy of TiO2.

Treating core electrons

Core and valence electrons

Typically there are two main parameters to consider when treating electrons as core or valence that effect how the wave functions will be constructed: the core radius and the number of core vs. valence electrons. The core radius is the spherical distance from the center of the atom that forms the outer bound for the approximation to the wave functions. The number of electrons in the core vs. valance is simply the partition of how many electrons of charge density an atom can contribute to bonding. For example, Ti has 22 electrons that must be partitioned between core and valance. Both of these parameters are typically chosen by the person developing the new description.

A basic description of some methods for treating core electrons in extended systems follows.

Pseudopotentials

Pseudopotentials replace the core wave functions with a pseudo wave function with the goal of defining a smooth function to represent the effective potential of the core. The pseudo wave function are constructed so that the value and the derivative of the pseudo wave function matches the value and derivative of the true wave functions at the core radius. The advantage of pseudopotentials is their formal simplicity, however the pseudo wave function prevents the recovery of the full core electron wave functions. James Goff published a post on pseudopotentials [2] that describes their construction in a detailed, introductory way.

LAPW

Another method for treating core electrons is the Linear Augmented Plane Wave (LAPW) method. The main idea of this method is that two partial waves, which can be used to describe the all-electron wave functions, are joined at the core radius and matched to the all-electron solution. Thus, the main advantage is the core wave functions can be recovered, but this method comes at a higher computational cost than the pseudopotential method. A more complete description of the pseuodpotential and LAPW methods can be found in [3].

PAW

The Projector Augmented Wave (PAW) method presented by Blöchl [4] as a formalism that generalizes and combines favorable aspects of the pseudopotential and LAPW methods through projector functions. The basis of the PAW method is the combination of partial wave functions from the isolated atom with pseudo partial waves. The PAW method defines a transformation operator that maps the pseudo wave functions onto the all-electron wave functions,
\begin{equation} \left| \Psi \right> = \mathcal{T} \left| \tilde\Psi \right> \end{equation}
for the all-electron wave function, \( \left| \Psi \right> \), the transformation operator, \( \mathcal{T} \), and the pseudo wave functions \( \left| \tilde\Psi \right> \). The transformation operator is defined as
\begin{equation} \mathcal{T} = 1+ \sum_i \left( \left| \phi \right> – \left| \tilde\phi \right> \right) \left< \tilde p_i \right|\end{equation}
where \( \left| \phi \right> \) and \( \left| \tilde\phi \right> \) are the all-electron partial waves and pseudo partial waves, respectively, and \( \left< \tilde p_i \right| \) are the projectors that are constructed to match the all-electron solution. In the PAW method, PAW datasets are generated that can recover the full core electron wave functions, like the LAPW method, but at a lower cost than the LAPW method, closer to the computational cost of pseudopotentials. More information can be found for the PAW method in [4-6] and more detailed description of the relationship between pseudopotentials, APW, and the PAW methods can be found in [4].

Calculation details

Calculations are performed using plane waves and the PBE functional implemented by Quantum Espresso [8] in the pwscf code. PAW datasets [4] are used to describe the core states. For Ti, PAW datastets are generated as part of this work with the ld1 code distributed with Quantum Espresso, based on the PS Library 1.0 parameterizations [9]. The PS Library 1.0 parameterization for O is used directly. The reference states for Ti and O are chosen as hcp Ti and the O2 molecule in a 10x10x10 fcc box. The size of the O2 box was converged to within 0.1 mRy. Detailed discussion of the convergence of the box size is out scope for this post.

The cutoff energy convergence of the different potentials are compared in the following sections below for Ti and TiO2. Oxygen cutoff energy was converged to 0.1 mRy at a wavefunction cutoff energy of 90 Ry and a kinetic energy cutoff of 720 Ry. For the calculation of formation energy, the maximum converged cutoff energy for each set of potentials is used across all structures. In all cases, Ti required a larger cutoff energy for convergence than O, so it was used for the calculations of TiO2.

The kpoint grids were converged to within 0.1 mRy. For hcp Ti, an 11x11x7 Monkhorst-Pack grid was used. For the O2 box, a 3x3x3 Monkhorst-Pack grid, and for rutile TiO2 a 4x4x7 grid. Showing the detailed kpoints convergence tests are out of scope for this study. Even though the different PAW datasets can lead to different electronic structures and kpoints convergence, a single set of kpoints for each set of structures were used across all of the different PAW datasets for simplicity.

For comparing the cutoff energy of hcp Ti, calculations with each PAW potential used the same structure, a hexagonal primitive cell with a lattice parameter of 5.8 a.u. and a c/a ratio of 1.58 (corresponding to the experimental c/a ratio). All static calculations were performed using tetrahedron integration and an energy convergence of 1e-8 Ry. For calculations of the formation energy, the cells and ionic positions for all cells (O, Ti for each potential, TiO2 for each potential) were relaxed with a force convergence criteria of 1e-9 Ry/bohr. The Brillouin zone integration for relaxations was treated with Gaussian smearing with a smearing width of 0.001 Ry.

 

Generating PAW datasets for Ti

The ld1 program was used to generate the PAW datasets. The documentation can be found online [10]. The starting point for the tests here is the input file for generating a PAW dataset for the spin polarized, scalar relativistic [12] Ti from the PS Library version 1.0.0 [9]. The input file for the PS Library starting point follows below. The settings here create a PAW dataset for Ti with Troullier-Martins pseudization [11] via Bessel functions. The the key inputs that will be changed is `rcore`, which is the matching radius in atomic units for the core charge smoothing. The bottom section describes the number of wave functions (6) that will be pseudized and which are described by the lines below. Taking the first line (3S) as the example, the first number is the principle quantum number (n=1, starting with 1 for the lowest s, 2 for the lowest p), then the angular momentum quantum number (l=0), followed by the number of occupying electrons, the energy used to pseudize the state (=0.00), the matching radius for the norm-conserving pseudopotential (=0.85), the matching radius for the ultrasoft pseudopotential (=1.30) and finally the total angular momentum (=0.0). Changing the norm-conserving or ultrasoft pseudopotential matching radii do not change the PAW dataset, since the radii are controlled by the input parameters above. Note that the energies for the unbound states (4p=4 Ry and unoccupied 3d=0.05 Ry) need to be non-zero.

&input
title='Ti',
zed=22.,
rel=1,
config='[Ar] 4s2 4p0 3d2',
iswitch=3,
dft='PBE'
/
&inputp
lpaw=.true.,
pseudotype=3,
file_pseudopw='Ti-PSL.UPF',
author='ADC',
lloc=-1,
rcloc=1.6,
which_augfun='PSQ',
rmatch_augfun_nc=.true.,
nlcc=.true.,
new_core_ps=.true.,
rcore=0.8,
tm=.true.
/
6
3S 1 0 2.00 0.00 0.85 1.30 0.0
4S 2 0 2.00 0.00 0.85 1.30 0.0
3P 2 1 6.00 0.00 0.80 1.60 0.0
4P 3 1 0.00 4.00 0.80 1.60 0.0
3D 3 2 2.00 0.00 0.80 1.50 0.0
3D 3 2 0.00 0.05 0.80 1.50 0.0

To demonstrate how the cutoff radius affects the calculation accuracy and expense (through the cutoff energy), two new PAW datasets will be generated with core radii of 0.4 a.u. (much smaller, further referred to as “small core”) and 2.0 a.u. (much larger, further referred to as “large core”). Note that these are note expected to be accurate, but instead to demonstrate how the parameters affect might affect calculations and impress the care and detail that should be given to generating proper potentials.

Energy cutoff convergence of Ti and formation energy TiO2

Figure 1a and 1b compre the relative and absolute energy convergence with respect to the wave function energy cutoff for the differently generated Ti PAW datasets. In Figure 1a, the PS Library converges at a wave function cutoff energy of 120 Ry, the large core at 110 Ry, and the small core does not converge to less than 0.1 mRy/atom even at 140 Ry. Computational demands and the instability of the PAW dataset limited the full convergence.

The cutoff energies follow the trend that would be expected. Compared to the default PS Libary PAW dataset, increasing the core radius allows for more of the wave function to be described by the core, and therefore the dataset is softer (requires smaller energy cutoffs). In addition, the small core is not converged even at 140 Ry because the smaller radius requires a harder potential with more plane waves being required to describe the wave function. Figure 1a also shows unusual behavior for the energy to not decrease variationally with the cutoff energy and instead shows that there was some osciallation in the energy with increasing energy cutoff. Note that the large core and the the PS Library default are quite close in energy, so the differences may be hard to see in both figures. Figure 1b aims to demonstrate that the energies for different pseudopotentials are different even for the same calculation settings, so they cannot be compared.

Figure 1a: Energy cutoff convergence relative to the energy calculated at a cutoff energy of 140 Ry. Energy differences are taken from 10 Ry to 140 Ry. Calculations not shown are outside the bounds of the plot.

Figure 1b: Energy convergence for cutoff energies ranging from 10 Ry to 140 Ry. Missing results at lower energy cutoff are off the scale of the plot.

As discussed previously, energies between different pseudopotentials cannot be directly compared, therefore a metric that can cancel the differences is desired. The formation energy of TiO2 can cancel the differences caused by the different treatment of Ti and give some insight into the practical implications of different treatments of the electronic structure. Using the cutoff energies found for Ti (120 Ry for PS Library, 140 Ry for small core, 110 Ry for the large core) the energies per atom of hcp Ti, O2 in a vaccuum box, and TiO2 are calculated. The formation energy of TiO2 is calculated in the usual way:
\begin{equation} E^\mathrm{form} = E_{\mathrm{TiO_2}} – (E_{\mathrm{Ti}} + E_{\mathrm{O_2}}) \end{equation}
Since in each case, the cutoff energies are convered for all calculations within each formation energy calculation, the formation energies can be safely compared even with different cutoff energies used between them.

StructurePS LibrarySmall coreLarge core
hcp Ti (Ry/atom)-183.6910-183.6974-183.6975
O (Ry/atom)-41.4776-41.4776-41.4776
TiO2 (Ry/atom)-89.1361-89.1363-89.1357
TiO2 formation (Ry/atom)-0.2541-0.2521-0.2515

The total energies and formation energies are shown in Table 1. The differences in the formation energy are larger than the convergence criteria (0.1 mRy/atom, 0.0001 Ry/atom) for the formation energies of the PS Library to the small core PAW dataset and the PS Library to the large core PAW dataset. The largest difference is 0.0026 Ry/atom for the small core and PS Library formation energies, which corresponds to above 3 kJ/mole-atom. Differences that large could impact the relative stability of TiO2 in a more comprehensive study that used inaccurate PAW datasets.

Conclusion

PAW datasets were reviewed in context of other types of frozen core approximations, then PAW datasets were generated for Ti with larger and smaller core radii than the ones parameterized in the PS Library v1.0.0. It was shown that the different core radii can affect the energy cutoff convergence, the total energy, and the formation energy of a compund, TiO2.

References

[1] Section 10.5, D.S. Sholl, J.A. Steckel, Density Functional Theory, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2009. doi:10.1002/9780470447710.
[2] J. Goff, Transferability of Cu Pseudopotentials in Cu/CuO Systems, Density Funct. Theory Pract. Course. (2019). https://sites.psu.edu/dftap/2019/03/31/transferability-of-cu-pseudopotentials-in-cucuo-systems/.
[3] D.J. Singh, L. Nordstrom, Planewaves, Pseudopotentials and the LAPW Method, Second Edi, Springer US, 2006. doi:10.1007/978-0-387-29684-5.
[4] P.E. Blöchl, Projector augmented-wave method, Phys. Rev. B. 50 (1994) 17953–17979. doi:10.1103/PhysRevB.50.17953.
[5] X. Gonze, F. Finocchi, Pseudopotentials Plane Waves-Projector Augmented Waves: A Primer, Phys. Scr. T109 (2004) 40. doi:10.1238/physica.topical.109a00040.
[6] R.M. Martin, Electronic Structure: Basic Theory and Practical Methods, Cambridge University Press, 2004. doi:10.1017/CBO9780511805769.
[7] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (1996) 3865–3868. doi:10.1103/PhysRevLett.77.3865.
[8] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A.P. Seitsonen, A. Smogunov, P. Umari, R.M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials., J. Phys. Condens. Matter. 21 (2009) 395502. doi:10.1088/0953-8984/21/39/395502.
[9] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Comput. Mater. Sci. 95 (2014) 337–350. doi:10.1016/j.commatsci.2014.07.043.
[10] http://web.mit.edu/espresso_v6.1/amd64_ubuntu1404/qe-6.1/Doc/INPUT_LD1.html
[11] N. Troullier, J.L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B. 43 (1991) 1993–2006. doi:10.1103/PhysRevB.43.1993.
[12] Koelling, D. D. & Harmon, B. N. A technique for relativistic spin-polarised calculations. J. Phys. C Solid State
Phys. 10, 3107–3114 (1977).

Exploring the normal modes of gaseous ammonia

The aim of this post is to determine and visualize the vibrational modes of an isolated ammonia molecule. After a discussion of the theory of normal modes, the frequencies will be calculated directly using DMol3 in Materials Studio. Materials Studio also provides animations of the normal modes to help with visualization of the modes.

Normal Modes

Consider a classical system with generalized coordinates \( \{q_i \} \) for \( i=1\dots n\) moving under the influence of a potential energy \( V(q_i) \). Classical mechanics tells us that the equations of motion for this system are

\begin{equation} m_i \ddot{q}_i=-\frac{\partial V}{\partial q_i} \end{equation}

where \( m_i \) is the mass associated with the generalized coordiante \( q_i \). The fixed points of this system occur when \( \ddot{q_i}=0 \) for all \(i\) or equivalently,

\begin{equation} \frac{\partial V}{\partial q_i}=0 \end{equation}

Of course, this means that the fixed points of the system occur when the system is configured such that \(V\) is at a critical point. At such a point, the system can sit at equilibrium. Denote the configuration of the generalized coordinates that \( (2) \) holds as \( \{q^0_i\} \). If we suppose that our system is in a configuration close to \( \{q^0_i\} \), we can expand the potential in a Taylor series centered at \( \{q^0_i\} \) and truncate it at second order in \( q\):

\begin{equation} V=V_0+\frac{1}{2} \sum_{jk} \frac{\partial^2 V}{\partial q_j \partial q_k} (q_j –  q^0_j)(q_k –  q^0_k) \end{equation}

where the derivatives are evaluated at \( \{q^0_i\} \). Under such circumstances, the equations of motion reduce to a simple and familiar form:

\begin{equation}  \ddot{q}_i=-\sum_j \frac{1}{m_i} \frac{\partial^2 V}{\partial q_i \partial q_j}(q_j –  q^0_j) \end{equation}

These equations of motion describe \(n\) coupled simple harmonic oscillators. It is useful to introduce the following change of coordinates:

\[x_i=\frac{1}{\sqrt{m_i}}(q_i –  q^0_i)\]

The equations of motion then become:

\begin{equation} \ddot{x}_i=-\sum_j \frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}x_j \end{equation}

The solutions of the above equation form an \(n\)-dimensional vector space, and it is possible to construct an orthogonal basis for this vector space. Using the ansatz \( x_i(t)=x_i(0) \sin(\omega t+\delta) \), where \( \delta \) is a constant that depends on initial conditions, the equation becomes:

\begin{equation} \omega^2 x_i(0)=\sum_j \frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}x_j(0) \end{equation}

If we think of \(x_i\) as vectors, this equation says that \(x_i(0)\) is an eigenvector of the \(n \times n\) symmetric matrix \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) with eigenvalue \(\omega^2\). A theorem of linear algebra guarantees that the eigenvectors of a symmetric matrix can be made pairwise orthogonal. Therefore the \(n\) eigenvectors of \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) form an orthogonal basis for the space of solutions of \((5)\). Denote these solutions as

\begin{equation} x^a_i(t)=e^a_i \sin(\omega_a t+\delta_a) ~~(a=1,\dots,n) \end{equation}

where \(a\) is a label for the eigenvalues and eigenvectors. Any solution of \((5)\) can be written as a unique linear combination of the \(x^a_i(t)\):

\begin{equation} y_i(t)= \sum_{a=1}^n C_a x^a_i(t) \end{equation}

The solutions \(x^a_i(t)\) are called normal modes of the system and the \(\omega_a\) are called normal frequencies. Equation \((8)\) highlights for us the use of normal modes: given initial conditions of the system, we can expand in terms of the normal modes to find a solution compatible with the initial conditions. It is worth remembering, however, that this approach is only valid when the coordinates are close to the fixed points of \(V\).

Determing the normal frequencies

To determine the normal frequencies, we need to know the partial derivatives of the potential energy evaluated at the critical points. This can be accomplished using DFT calculated energies for a molecular system as follows:

\begin{equation} \frac{\partial^2 V}{\partial q_i \partial q_j} \approx \frac{E(q^0_i+\delta q_i,q^0_j+\delta q_j)-E(q^0_i+\delta q_i,q^0_j-\delta q_j)-E(q^0_i-\delta q_i,q^0_j+\delta q_j)+E(q^0_i-\delta q_i,q^0_j-\delta q_j)}{4(\delta q_i)(\delta q_j)} \end{equation}

This formula is exact in the limit \(\delta q_i,\delta q_j \rightarrow 0\). After determining the derivatives of \(V\) we can diagonalize the matrix \(\frac{1}{\sqrt{m_i m_j}} \frac{\partial^2 V}{\partial q_i \partial q_j}\) to determine its eigenvalues. Note that nothing restricts the sign of these eigenvalues and in fact only \(\omega_a^2>0\) correspond to oscillatory motion. Materials Studios has a function to compute these positive eigenvalues from a geometry optimization calculation. These are the modes of ammonia we will focus on in this post.

Modes of ammonia

Figure 1: A relaxed ammonia molecule NH3. The atoms are arranged in a trigonal pyramidal shape.

We performed a DMol3 geometry optimization calculation with vibrational analysis in Materials Studio to determine the oscillatory frequencies of an isolated ammonia molecule. The functional used was GGA-PBE. The energy convergence tolerance was \(10~\mu Ha\), the max. force was .002 Hartrees per Angstrom, and the max. displacement was .005 Angstroms. For this system, \(n=12\) and there are 6 positive eigenvalues of interest. The results, along with the comparisons with experimental values, are shown below:

ModeDFT Result (1/cm)
N-H wagging1064.55
H-N-H scissoring (1)1631.18
H-N-H scissoring (2)1631.18
N-H symmetric stretch3390.54
N-H asymmetric stretch (1)3524.06
N-H asymmetric stretch (2)3524.06

Materials Studio also gives animations of these modes. Click on any of the images below to view an animation of the normal modes.

Figure 2: N-H wagging

Figure 3: H-N-H scissoring (1)

Figure 4: H-N-H scissoring (2)

Figure 5: N-H symmetric stretch

 

Figure 6: N-H asymmetric stretch (1)

Figure 7: N-H asymmetric stretch (2)

References

[1] PBE functional: Perdew Burke Ernzerhof: Phys. Rev. Lett. 77, 3865 (1996)

[2]  The Generation and Use of Delocalized Internal Coordinates in Geometry optimization;
Baker Kessi Delley: J. Chem. Phys., 105, 192 (1996)
Andzelm Fitzgerald King-Smith: Chem. Phys. Lett., 335, 321 (2001)

[3] Parallel solution of partial symmetric eigenvalue problems from electronic structure calculations:
Auckenthaler, Blum, Bungartz, et al.: Parallel Computing 37, 783 (2011)

Considerations for DFT Frequency Calculations

Using DFT to Calculate Vibrational Modes/Frequencies

In this post we explore the effect of numerical method as well as displacement magnitude on the results of vibrational mode calculations using DMol3.

The vibrational modes of a system or collection of atoms/molecules are determined from the Hessian Matrix, H, which has elements [1]:

\begin{equation}H_{ij}=\frac{1}{\sqrt{m_im_j}}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}}\end{equation}

In this notation qi are the coordinates of atom i. The eigenvectors of H are the vibrational modes themselves (displacements of atoms along the normal mode) while the vibrational frequencies are given as the square roots of the eigenvalues.

Numerical Methods

Most DFT packages (DMol3 included) cannot calculate an analytic function for the energy, E with respect to all of the degrees of freedom (coordinates). Thus, numerical approximations to the derivatives in the above expression need to be used.

DMol3 offers two options for the numerical method used in estimating such derivatives; these are the Single-Point Difference (Forward Difference) Method and the Two-Point Difference (Central Difference) Method.

The Single-Point expression is [1]:

\begin{equation}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}} \approx \frac{\frac{\partial{}}{\partial{q_i}}E(q_i+\delta)-\frac{\partial{}}{\partial{q_i}}E(q_i)}{\delta}\end{equation}

The Two-Point expression is [1]:

\begin{equation}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}} \approx \frac{\frac{\partial{}}{\partial{q_i}}E(q_i+\delta)-\frac{\partial{}}{\partial{q_i}}E(q_i-\delta)}{2\delta}\end{equation}

In these expressions, δ is the displacement length (magnitude) which is a calculation parameter that can be adjusted by the user.

While the documentation for DMol3 suggests the Central Difference method may reduce rounding-errors (while doubling the computational load), it gives little explanation as to any major differences in results obtained by each method.

Similarly, other than mentioning using smaller (larger) values of δ for steeper (flatter) potential energy surfaces, little insight is given as to the effect of the displacement length on the results of these vibrational mode calculations.

Thus, in this post we will use DMol3 to calculate the vibrational modes and frequencies of methane (CH4) using both the Single-Point (SP) method and the Two-Point (TP) method as well as the effect of the magnitude of displacement.

Vibrational Modes of Methane

As methane is a non-linear molecule with 5 atoms, there are 9 vibrational modes (3N-6, N=5), however; due to symmetry (Td point group for methane) there are only 4 unique vibrational frequencies. For reference, these modes are visualized below with their frequencies taken from the NIST Webbook [2].

Fig. 1: A1 Symmetry Mode, f=2917 1/cm

Fig. 2: E Symmetry Mode, f=1534 1/cm

Fig. 3: F2 Symmetric Stretch Mode, f=3019 1/cm

Fig. 4: F2 Symmetric Bend Mode, f=1306 1/cm

 

Calculation Details

All calculations (geometry optimization, frequency analysis) were performed using the DMol3 module as implemented in Materials Studio using the PBE functional for exchange and correlation with the DNP+ basis set [3-6]. SCF cycles were converged to 1.0-5 eV and the forces on the atoms were converged to 0.001 eV/Å. Core electrons were treated explicitly (all electron).

Results & Discussion

Before calculating the vibrational modes, we must first geometry optimize the methane molecule. This was done in DMol3 using the DNP+ basis set (for non-periodic systems). Below are images of the initial configuration as well as the optimized structure.

Fig. 5 Initial methane configuration

Fig. 6 Optimized methane structure

Below are plots of the various vibrational frequencies determined via the two methods for various δ (frequency as a function of 1/ δ is used for clarity). The frequencies (f#) are grouped according to symmetry (e.g. frequencies 1, 2, and 3 should be the same so they appear on the same plot). The choice of “which frequencies are identical” is  The displacement values used were 0.005, 0.0075, 0.01, 0.015, 0.02, 0.025, 0.05, 0.1, 0.25, and 0.5 Å.

Fig. 7: Vibrational Frequencies for Single Point Method various displacements.

Fig. 8: Vibrational Frequencies for Two Point Method at various displacements.

In the table below, the frequency for each mode in each method are listed. The frequencies were determined by averaging the frequencies for which the change in the frequency with displacement was ~10 cm-1 (the level of precision in acceptable experimental results). In the case of the both the SP and TP methods, all of the frequencies were converged for δ < 0.02 Å, though individual modes may have converged more quickly.

Table 1: Vibrational Frequencies for SP and TP Methods
ModeSP Frequency (1/cm)TP Frequency (1/cm)
11311.741315.76
21326.901315.76
31326.901315.76
41531.281535.83
51549.901535.83
62990.582996.03
73108.943116.29
83112.343116.29
93112.363116.29

First, it should be noted that while both the SP and TP predict 9 vibrational modes, the single-point method does not quite capture the reality of the situation. As mentioned earlier, due to belonging to the Tpoint group, there should be 4 unique modes. Using the numbering above, frequencies 1, 2, and 3 should be identical as should the pair 4, 5 and the trio of frequencies 7, 8, and 9. These modes are equivalent for all calculations done using the two-point method, as is observed by the overlapping of the points in Fig. 8. The SP method eventually converges to have approximately equal frequencies for the symmetric modes though using larger displacement values increases this slight inaccuracy.

In these cases, the single-point method can either predict a larger or smaller frequency than the two-point method with the former having some variation in the convergence patterns. This highlights how drastically calculation parameters can affect results and conclusions drawn from such calculations.

From the data above, we notice that for either the SP or TP method, as the displacement magnitude is lowered, the frequency of each mode eventually “converges” to a stable value. This is to be expected as smaller displacements are closer to the equilibrium bond lengths and the harmonic approximation is valid, while at large separations there are sever deviations from the simple harmonic approximation.

It is worth noting that had we used looser convergence criteria (or even a smaller basis set) we might expect the behavior of the frequency for very small displacements to be different. This is due to the fact that eventually the displacement size, and thus the changes in energy, would be comparable to or less than the machine precision possible using the “looser” calculation parameters. That is to say, using a smaller basis set or less demanding criteria for the energy and force convergence, the energy differences calculated numerically using the SP and TP methods would be non-significant and could produce unphysical results.

Given these arguments, our calculations using TP and SP methods are still in decent agreement, however, the above exercise is insightful in that it shows that convergence criteria for frequency calculations can have a large impact on the results, especially when using larger displacements or for larger molecules having more vibrational modes.

The most important conclusion we can draw from this work is that the single-point and two-point methods differ mostly in their accuracy of predicting symmetric modes, with the TP method outperforming the SP method. This suggests that when calculating the vibrational modes and frequencies of high-symmetry molecules, the TP method is better suited.

Finally, we can see that the calculated results for both are in decent agreement with the experimental values but with errors around 5%- for sufficiently small displacements. This may serve as a reminder that computation and simulation can help corroborate experimental findings but generally cannot predict exact values.

References

[1]     Accelrys Software Inc., Materials Studio Release Notes, Release 7.0, San Diego: Accelrys                      Software Inc., 2014.

[2]     Shimanouchi, T., Tables of Molecular Vibrational Frequencies Consolidated Volume I,                      National Bureau of Standards, 1972, 1-160.

[3]     B. Delley, J. Chem. Phys. 92, 508 (1990).

[4]     
B. Delley, J. Chem. Phys. 113, 7756 (2000).

[5]     
J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made                        simple. Phys. Rev. Lett., 77:3865, 1996.

[6]     J. P. Perdew, K. Burke, and M. Ernzerhof. Erratum: Generalized gradient approximation                 made simple. Phys. Rev. Lett., 78:1396, 1997.

 

Impact of the Zero Point Energy Correction to Hydrogen Adsorption Sites on the Cu (111) Surface

Purpose of Calculation

In this post, we calculate the zero point energy (ZPE) correction to an adsorbate atom on a surface. The ZPE contributes to the energy of the system, and can affect which binding site on a surface corresponds to the minimum energy arrangement. The ZPE arises from the fact that the adsorbed atom is free to oscillate while still bound to the surface, behaving as a quantum mechanical oscillator. One of the most important features of such an oscillator is, of course, that its lowest energy mode is not a zero energy mode, but a finite energy mode. So, when calculating a system with an adsorbate, we this finite energy (the ZPE) becomes a correction to our calculation. Since the correction energy is the energy of a harmonic oscillator, we expect that the energy should be largely determined by our spring constant (the “stiffness” of our bond” and the mass of the bonded atom. Therefore, the ZPE correction should be most important to atoms such as hydrogen. To demonstrate this, we calculate the ZPE correction for Hydrogen adsorbed at two sites on the Cu (111) surface: the hcp site and the fcc site.

Figure 1: A supercell (2×2) of Hydrogen (red) adsorbed on the (111) surface of copper

Theory of the Zero Point Energy

The ZPE for a quantum mechanical oscillator is given by the sum of the frequencies of the vibrational modes of the system, multiplied by constants. These vibrational modes can be obtained from the eigenvalues of the mass-weighted Hessian matrix of the system. The Hessian matrix is composed of double partial derivatives of the system energy with respect to the displacement of atoms in the system. The energy we calculate for each atomic displacement is in eV, and energy in general relates to time through energy equals mass times distance squared, divided by time squared. The Hessian matrix removes the distance dependence of the Hessian elements. If we then divide out the mass dependence, the elements are in units of inverse seconds squared. If we then take the square root of the eigenvalues of this system, the resulting quantity will be in units of inverse time, as our frequencies should be. In our system, since we are assuming the copper surface’s movement during vibrations will be negligible compared to the hydrogen atom’s movement, the hydrogen atom should be displaced along the x, y, and z direction with an appropriate step size to calculate the second derivatives of energies with respect the movements.

Numerically, we calculate the derivatives using the finite difference approximation. The diagonal elements can be calculated using the first equation below, and the off-diagonal elements can be calculated using the second.[3]

\begin{equation} H_{ii} \cong \frac{E(\delta x_{i})+E(-\delta x_{i})-2E_{0}}{\delta x_{i}^2} \end{equation}
\begin{equation} H_{ij} \cong \frac{E(\delta x_{i}, \delta x_{j})+E(-\delta x_{i},-\delta x_{j})-E(\delta x_{i}, -\delta x_{j})-E(-\delta x_{i}, \delta x_{j})}{4\delta x_{i}\delta x_{j}} \end{equation}

 

Delta x represents the displacement of the hydrogen atom along direction i or j, and E-naught is the energy when hydrogen is not displaced.

Now, we can diagonalize the Hessian matrix, obtaining the eigenvalues, from which we can directly calculate the vibrational frequencies, using [3]:

\begin{equation} \nu_{i} = \frac{\sqrt{\lambda_{i}}}{2\pi} \end{equation}

where lambda represents the indexed eigenvalue, and nu is the corresponding frequency. These eigenvalues are of the mass-weighted Hessian matrix, which can be obtained be dividing the Hessian matrix by the square root of the mass of the atoms that are being displaced. In our case, this means we should divide by the square root of the mass of hydrogen. Each frequency then contributes an energy

\begin{equation} E_{ZPE} = \frac{h\nu_{i}}{2} \end{equation}

to the ZPE, where h is Planck’s constant.

Calculation Methodology

The ZPE correction to hydrogen on the copper (111) surface is calculated at two sites. For these calculations, we cleave the (111) surface from bulk copper, and then form a super 2×2 supercell, so that both sites are present with all nearest neighbors for interactions. We include three layers of Cu atoms. The lattice parameters are 5.11 angstroms, 5.11 angstroms, 14.17 angstroms (a, b, c). The calculation cell contains 12 unique copper atoms and one unique hydrogen atom. We perform the calculations using Dmol3, which uses a localized basis set [1]. All calculations are performed with the Perdew-Burke-Ernzerhof (PBE) [2] exchange-correlation functional, a Generalized Gradient Approximation (GGA) functional.The on-the-fly basis set is calculated by generating a radial function calculated by solving the density equations on the fly at three hundred points on a logarithmic radial mesh for each atom. This radial function is then multiplied by appropriate spherical harmonics for each atom. We use the DNP basis set method for this calculation. All electrons are included in the Dmol3 calculations. Our SCF threshold is 10^-6 Ha.

We select a 4x4x2 k-point mesh, which test calculations show to be converged to within 0.025 Ha. The cutoff for the localized basis set is set to 4.9 angstroms, and ten angstroms of vacuum space is included in the supercell. We select 0.05 angstroms as our displacement distance. This displacement is selected after calculating the diagonal elements for the Hessian matrix for displacements of 0.01, 0.05 and 0.1 angstroms. 0.01 angstrom displacements did not provide significant changes in the calculated energy relative to the level to which energy was converged. At 0.05 angstroms, there is a significant change in calculated energy, so we elect to use 0.05 angstroms, so that our displacement does not become so large that the finite difference approximation of the derivative is no longer a good approximation.

Results

Calculations were performed for both the hcp and fcc sites with 0.05 angstroms of displacement in each possible combination of directions. Results are shown Table 1.

Displacements (Angstroms)Energy (Ha)
XYZhcpfcc
000-19685.479636-19685.485998
-0.0500-19685.478748-19685.485904
0.0500-19685.480315-19685.485877
0-0.050-19685.480315-19685.485869
00.050-19685.478748-19685.485911
00-0.05-19685.479381-19685.485882
000.05-19685.479387-19685.485759
-0.05-0.050-19685.478283-19685.485797
-0.050-0.05-19685.478244-19685.485762
0-0.05-0.05-19685.479899-19685.485717
-0.050.050-19685.4791136-19685.485798
-0.0500.05-19685.47871-19685.485685
0-0.050.05-19685.477772-19685.485658
0.05-0.050-19685.479407-19685.485774
0.050-0.05-19685.480249-19685.485728
00.05-0.05-19685.479387-19685.485771
0.050.050-19685.479321-19685.485763
0.0500.05-19685.479899-19685.485665
00.050.05-19685.477165-19685.485691

From these values, we can then calculate the appropriate Hessian matrices for each site. As a visual example, the hcp Hessian is

\begin{bmatrix} 8.36\times10^{18} & 9.166\times10^{18} & 8.16\times10^{18} \\ 9.166\times10^{18} & 8.36\times10^{18} & 9.5\times10^{17} \\ 8.16\times10^{18} & 9.5\times10^{17} & 2.016\times10^{19}\end{bmatrix}

Here, all the elements are in units of the square root of mass multiplied by inverse seconds square. Dividing all the elements by square root of the mass of hydrogen and diagonalizing, we obtain the eigenvalues: 5.1258×10^31, 3.1687×10^32 and 6.361×10^32 for the hcp site. For the fcc site, the eigenvalues are 2.0753×10^32, 2.1371×10^32 and 3.4746×10^32. The total ZPE correction for each of these eigenvalues is 0.607 eV for the hcp site and 0.576 eV for the fcc site.

These ZPE corrections are significantly large. Our convergence for the energy calculation, 0.0025 Ha corresponds to approximately 0.068 eV. The corrections are almost an order of magnitude larger than our energy convergence, meaning they represent a correction that should absolutely be included in our calculation. Also of note are that the hcp and fcc corrections differ but hundredths of an eV. When we minimized the energies for each location of the hydrogen atom before displacement, the hcp position was lower than the fcc site by approximately 0.15 eV. As a result, this correction does not change the preferred site from hcp to fcc, but it is still a significant contribution, and illustrates how the ZPE correction can be vital to the proper treatment of systems. With a more stringent convergence of the energy in the calculation, we would be able to better resolve the exact ZPE correction, but at this time, the computational resources were not sufficient to do so in a time-efficient manner.

In addition, possible other sources of error include contributions from copper atoms that Dmol3’s localized basis set do not properly capture, isolation of the hydrogen atom and other similar concerns. In an actual adsorption reaction, there will not be a single hydrogen atom with no other nearby hydrogen on the surface. Instead, something like half the sites might be occupied, or some other, more complex situation. This calculation does not reflect such concerns.

References

[1]  Delley, J., Fast Calculation of Electrostatics in Crystals and Large Molecules; Phys. Chem. 100, 6107 (1996)

[2]  John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple”, Phys. Rev. Lett. 77, 3865 – Published 28 October 1996; Erratum Phys. Rev. Lett. 78, 1396 (1997)

[3] D. Sholl and J. Steckel, Density Functional Theory: A Practical Introduction. (Wiley 2009)

Hydrogen atoms adsorbed on Cu(111) surface

In this project, the absorption of H atom in the Cu (111) surface is analyzed and calculated. The H atom absorbed on Cu(111) can bind in two distinct threefold sites, the fcc sites and hcp sites shown in Figure 1. The energy difference between these two sites are calculated using Castep and the vibrational frequencies of H in each site is analyzed by calculating the elements of the Hessian matrix using finite-difference approximations. Then based on the vibrational frequencies, the zero-point energy can be obtained.

Figure 1. The absorption of H atom in the Cu (111) surface. the grey atoms are the H atoms and the red atoms are the Cu atoms. The H atoms are in the hcp sites (a) and fcc sites (b).

1. The slab model

As shown in Figure 1, an asymmetric slab model of 4 layers are used for the calculation. During the calculation, the two bottommost layers are fixed at bulk positions and all remaining atoms allowed to relax at z direction to obtain the optimized structures.

2. The optimized structures

For both two sites, the geometry optimization in the Castep is used to obtain the optimized structures, which are shown in Figure 2. We can also get the energy difference between these two structures. for hcp structure, the energy is -6738.743eV. For fcc structure, the energy is -6738.755eV. So the energy difference is 0.012eV. And since the energy of hcp structure is larger than the fcc structure, the fcc structure is more stable than the hcp structure.

Figure 2. The optimized structures for H atoms at the hcp sites (a) and fcc sites (b). The distances are 1.729 Ångstrom in a and 1.733 Ångstrom in b.

The calculation parameters for geometry optimization are listed here:
Atomic and pseudo atomic structure for H: 1s1
Atomic structure for Cu: 1s2 2s2 2p6 3s2 3p6 3d10 4s1
pseudo atomic structure for Cu: 3d10 4s1
Functional: GGA PBE
Cut off energy: 408.2eV
k-point set: 9*9*1
Energy convergence tolerance: 1.0e-005eV/atom
Force convergence tolerance: 0.03eV/Å
Stress convergence tolerance: 0.05GPa
Displacement convergence tolerance: 0.001Å

3. The vibrational frequencies

To obtain the vibrational frequencies, we need to get the Hessian matrix first. Using finite-difference approximations, the elements in the Hessian matrix can be obtained by the equation:

Table 1 shows the calculated energy and the Hessian matrix. Then the frequencies can be obtained by the equation:

λi is the eigenvalue of Hessian matrix. So we can get the frequencies shown in Table 2. The vibrational frequencies can be seen from the wave numbers of the vibrations. All the energy shown in these two table have 3  significant digits after the decimal place because we can already see the energy differences. All the Hessian matrix are shown in SI.

Table 1 The calculated energy and the Hessian matrix for hcp structure (a) and fcc structure (b).

Table 2 The eigenvalues of Hessian matrix, the vibrational frequencies and the zero-point energy for hcp structure (a) and fcc structure (b).

4. The zero-point energy correction

Since we have get the vibrational frequencies, the zero-point energy can be calculated by the equation E=∑(hv)/2.So as shown in Table 2, the total zero-point energy for hcp structure is 0.185eV and the total zero-point energy for fcc structure is 0.192eV. Then after considering the zero-point energy, the energy difference between fcc and hcp structures is 0.0046eV.

5. Conclusion and Discussion

In this project, the absorption of H atoms binding in two distinct threefold sites, the fcc sites and hcp sites on Cu (111) surface is calculated. Without considering the zero-point energy, the energy difference for these two structures is 0.012eV. With considering the zero-point energy, the energy difference for these two structures is 0.0046eV.

In this project, I assume that the H vibrations were uncoupled from the metal atoms in the surface and calculate the vibrational frequencies only moving the H atom. But this assumption is too simple and the vibrations should be more complicated in a real system. For example, we can consider the coupling between H atom and neighboring Cu atoms, which would be an interesting project for the future.

6. Reference

[1] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, M. C. Payne, “First principles methods using CASTEP”, Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005)

Energy differences between hydrogen atoms absorbed on the Cu (111) fcc and hcp sites with(out) zero point energy included

1. Project Description

This project aims to calculate the vibrational frequencies and zero point energy (ZPE) of one hydrogen absorbed on face centered cubic (fcc) site and hexagonal close packed (hcp) site of Cu (111) surface. Energy calculation was performed using density functional theory (DFT) with CASTEP[1]. Exchange -correlation was treated in the generalized gradient approximation (GGA), as parameterized by Perdew, Burke, and Ernzerhof (PBE) [2]. The following valence electron configurations were employed in the potentials: 1s1 for H, and 3d10 4s1 for Cu.  Psedupotentials were used in place of other electrons, with core radii of 1.996Å  for Cu, and 0.599Å for H. The Brillouin zone integration was performed using a 9*9*1 Monkhorst-Pack k-point grid, and energy cutoff is 450eV. The SCF tolerance is 2.0e-6 eV/atom, and self-consistent dipole correction along the Z-axis is applied to the asymmetric slab.

2. Geometry optimization and energy calculation of two initial configurations

We will calculate the energies of hydrogen atoms adsorbed on two different threefold sites, the face centered cubic (fcc) site and hexagonal close packed (hcp) site, on Cu (111) surface. For each initial fcc and hcp configuration, we built four layers of Cu atom in the 10 Å slab and place one H atom near the surface (Figure 1&2).  Then we would optimize geometry with bottom three layers of Cu fixed, and both first layers of Cu and hydrogen moved freely. After geometry optimization of the structure,  the energy is calculated with the most stable relative Cu and hydrogen positions.

Figure 1. Top-down view of hydrogen atoms absorbed on FCC (Left) and HCP (Right) site of Cu (111) surface

Figure 2. Side view of hydrogen atoms absorbed on FCC (Left) and HCP (Right) site of Cu (111) surface after surface relaxation

3. Hessian Matrix

We aim to calculate the Hessian Matrix by deriving energies of various configurations with slight hydrogen displacement. Sholl and Steckel (2009) suggest a good finite-difference displacements range from 0.03-0.1Å [3]. In this project, we use finite-difference with 0.06 Å along x, y, and z direction, and energies will be calculated in 12 different configurations. The minimum energy for HCP and FCC configurations are -6738.805279 eV and -6738.808379 eV; and the atom position with minimum energy is marked with δx=δy=δz=0 (Table 1).

Table 1. Energies of displacements in HCP (Upper) and FCC (Lower) configurations, with a finite-difference of 0.06 Å

The Hassian matrix is calculated by second derivative of energies with displacement [3], for diagonal element, the Hessian matrix is calculated by:

\begin{equation} H_{ii}=(\frac{\partial^2 E}{\partial x_{i}^2}) \simeq \frac{E(x+\delta x_{i})+E(x-\delta x_{i})-2E_{0}(x)}{\delta x_{i}^2} \end{equation}

For off-diagonal element, the Hessian matrix is calculated by:

\begin{equation} H_{xy}=(\frac{\partial^2 E}{\partial x \partial y}) \simeq \frac{E(x+\delta x, y+\delta y)+E(x-\delta x, y-\delta y)-E(x-\delta x, y+\delta y)-E(x+\delta x, y-\delta y)}{4\delta x\delta y} \end{equation}

Thus the Hessian matrix for HCP and FCC configurations are:

4. Vibrational frequencies and zero point energy (ZPE) calculation

We can calculate the vibrational frequencies and zero point energy from the mass-weighted Hessian Matrix A, where the elements in the matrix A are [3]:

\begin{equation} A_{ij} = H_{ij} / m_i \end{equation}

The eigenvectors e and eigenvalues λ of matrix A satisfy Ae = λe. And vibrational frequencies νare obtained from the eigenvalues λi of the matrix product A[3],

\begin{equation} \nu_i = \sqrt{\lambda}/2\pi \end{equation}

Thus the eigenvalues and vibrational frequencies of the HCP and FCC mass-weighted Hessian matrix are (Table 2):

Table 2. The calculated eigenvalue and frequencies in HCP and FCC configurations

and zero point energy (ZPE) is obtained by [3]:

\begin{equation} E_{ZPE} = \Sigma_i \frac{ h \nu_i}{2} \end{equation}

Thus the zero point energy of HCP and FCC configurations are 0.175347eV and 0.623748eV.

5. Conclusion

The energies are summarized in Table 3. The energy differences between HCP and FCC configurations without ZPE are 0.003100eV, with ZPE are -0.445301eV. That’s to say, Without including ZPE, hydrogen on the FCC site has a lower energy than HCP sites (0.003100eV), which means hydrogen in FCC site is more stable. With ZPE included, hydrogen on the FCC site has a higher energy than HCP (0.445301eV), which means HCP site is more forved. Therefore, light atoms, such as hydrogen, should consider the influence of the zero point energy.

Table 3. A summary of energy with(out) ZPC when hydrogen absorbed on HCP and FCC sites.

 

Reference

[1]. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, M. C. Payne, “First principles methods using CASTEP”, Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005)

[2]. J. P. Perdew, K. Burke and M. Ernzerhof, “Generalized gradient approximation made simple”. Phys. Rev. Lett. 77, 3865 (1996).

[3]. Sholl, D. S., & Steckel, J. A. (2009). What is density functional theory? (pp. 1-33). John Wiley & Sons, Inc.

 

Energy and frequency calculation for hydrogen on Cu(111)

Introduction

Hydrogen atoms on Cu(111) can bind in two distinct threefold sites, HCP and FCC. This post is aimed at showing energy difference between these two sites without and with zero point energy correction by using DFT calculation.

For harmonic oscillator, in classical mechanics, the energy for equilibrium state is zero kinetic energy plus potential energy E0; But in quantum mechanics, due to Heisenburg uncertainty principle, the energy changes to: \begin{equation} E = E_0 + \Sigma_i \frac{h \nu_i}{2} \end{equation}

E0 is still the potential energy and ν is the vibrational frequency of the oscillator. Zero point energy is defined as: \begin{equation} E_{ZPE} = \Sigma_i \frac{ h \nu_i}{2} \end{equation}

Normally, zero point energy is not trivial for light atoms, like H.
In order to get zero point energy, vibration frequencies of H on Cu surface are needed and these frequencies are calculated from mass-weighed Hessian Matrix A, which has a relation with Hessian matrix H:
\begin{equation} A_{ij} = H_{ij} / m_i \end{equation}
where mi is mass of the atom associated with ith coordinate. In this post m =1.00794*1.66054*10^-27 kg. The eigenvalues of matrix A λ could give vibration frequencies by: \begin{equation} \nu_i = \sqrt{\lambda}/2\pi \end{equation}

Constructing Hessian Matrix requires us to do second-order differentiation for energy to coordinate: \begin{equation} H_{ij} = (\frac{\partial E}{\partial x_i \partial x_j})_{x=0} \end{equation}

which could be obtained from finite difference approximation. The method used in this post is the simple one which involving three points in one differentiation:

\begin{equation} H_{ij} \simeq \frac{E(\delta x_i, \delta x_j) – 2 E_0 + E(- \delta x_i, – \delta x_j)}{\delta x_i \delta x_j} \end{equation}

The model built in this post is H on the Cu(111) surface with three layers, vacuum above the surface is 10 Å. 2*2 supercell is adopted in case the hydrogen atoms in periodic structure are too close.

Before individual energy calculations, geometry optimization is conducted first. The hydrogen atom is put near the HCP site (or FCC site). Bottom two layers are fixed and the first layer and hydrogen are allowed to relax. After geometry optimization is finished, the distance between H atom and nearest three Cu atom could be checked. Fig. 1, 2 and 3 show the diagrams of the model. Fig. 1 shows that after geometry optimization, the distances between H atom and nearest three Cu atoms are close. So does Fig. 2. Fig.3 shows the side-view of the model.

 

Fig. 1 diagram for HCP site after geometry optimization

Fig. 2 diagram for FCC site after geometry optimization

 

Fig. 3 side-view of ‘H on Cu(111)’ model

 

Based on the geometry optimized, H atom is moved to three directions (x, y, z) according to elements in Hessian Matrix and energy calculation is conducted after each movement. The move step is 0.04Å.

Functional: GGA PBE[1]; k points: 5*5*1 ; cutoff energy: 408.2 eV; SCF tolerence: 2.0*10^-6 eV/atom; Pseudopotentials: OTFG ultrasoft, for Cu: inner shell is 1s2 2s2 2p6 3s2 3p6, outer shell is 3d10 4s1, for H: outer shell is 1s1; dipole correction is included since the model is not symmetric.

(‘CASTEP’ calculation module[2] is adopted for all calculations)

Calculation results

HCP site

Table 1 shows the calculation results at HCP site. δx means movement of H atom in x direction. δy and δz have similar meaning. δx=δy=δz=0 means the result of geometry optimization. ‘1’ means one movement which is 0.04Å.

Table 1.

δx δy δz Energy (eV)
0 0 0  -2.01835196E+004
1 0 0  -2.01834837E+004
1 1 0  -2.01834811E+004
1 0 1  -2.01835046E+004
0 1 0  -2.01834837E+004
1 1 0  -2.01834811E+004
0 1 1  -2.01835050E+004
0 0 1  -2.01835070E+004
1 0 1  -2.01835046E+004
0 1 1  -2.01835050E+004
-1 0 0  -2.01834834E+004
-1 -1 0  -2.01834815E+004
-1 0 -1  -2.01834479E+004
0 -1 0  -2.01834837E+004
-1 -1 0  -2.01834815E+004
0 -1 -1  -2.01834479E+004
0 0 -1  -2.01834499E+004
-1 0 -1  -2.01834479E+004
0 -1 -1  -2.01834479E+004

Mass-weighted Hessian matrix calculated based on data in Table 1:

[  4.31362677e+29          2.29142726e+29             2.59356062e+29]

[  2.29142726e+29          4.29567826e+29             2.58159494e+29]

[  2.59356062e+29          2.58159494e+29             4.92387633e+29]

Eigenvalues of the matrix are 9.51368753 e+29, 2.01335272 e+29, 2.00614111 e+29.

Vibration frequencies of H atom at HCP site (s^-1) are 1.5524*10^14, 0.7141*10^14, 0.7129*10^14.

Zero point energy for H at HCP site is 0.6161 eV.

 

FCC site

Table 2 shows the calculation results at FCC site.

Table 2.

δx δy δz Energy (eV)
0 0 0  -2.01835358E+004
1 0 0  -2.01834780E+004
1 1 0  -2.01834759E+004
1 0 1  -2.01835094E+004
0 1 0  -2.01834777E+004
1 1 0  -2.01834759E+004
0 1 1  -2.01835091E+004
0 0 1  -2.01835112E+004
1 0 1  -2.01835094E+004
0 1 1  -2.01835091E+004
-1 0 0  -2.01834783E+004
-1 -1 0  -2.01834765E+004
-1 0 -1  -2.01834259E+004
0 -1 0  -2.01834784E+004
-1 -1 0  -2.01834765E+004
0 -1 -1  -2.01834261E+004
0 0 -1  -2.01834280E+004
-1 0 -1  -2.01834259E+004
0 -1 -1  -2.01834261E+004

Mass-weighted Hessian matrix calculated based on data in Table 2:

[  6.89821314e+29          3.56577192e+29             4.07730464e+29]

[  3.56577192e+29          6.91017881e+29             4.08029606e+29]

[  4.07730464e+29          4.08029606e+29             7.92127857e+29]

Eigenvalues of the matrix are 1.51030097 e+30, 3.33850851 e+29, 3.28815230 e+29.

Vibration frequencies of H atom at FCC site (s^-1) are 1.9559*10^14, 0.9196*10^14, 0.9126*10^14

Zero point energy for H at FCC site is 0.7833 eV.

Conclusion:

Table 3 shows the energies of H at two three-hold sites on Cu(111) without and with zero point energy.

Table 3

Energy (eV) HCP FCC E(HCP) – E(FCC)
Without ZPE -20183.5196 -20183.5358 0.0162
With ZPE -20182.9034 -20182.7525 -0.1509

When zero point energy is not included, the H atom on FCC site has lower energy. But when zero point energy is included, the H atom on HCP site has lower energy. Zero point energy is important for light atoms, like H and should be considered for systems involving light atoms.

Reference:

  1. J. P. Perdew, K. Burke, and M. Ernzerhof Phys. Rev. Lett. 77, 3865 (1996)
  2. First principles methods using CASTEP” Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005) S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, M. C. Payne
  3. R. Ryberg, Carbon-Monoxide Adsorbed on Cu(100) Studied by Infrared-Spectroscopy, Surf. Sci. 114 (1982), 627 641.

Vibrations of Hydrogen on a Copper(1 1 1 ) Surface

In this project a hydrogen atom will be adsorbed to a Cu(1 1 1) surface and the vibration modes of the hydrogen atom will be calculated.  There are two sites where the hydrogen atom can be stably adsorbed, the fcc site and hcp site, as such both sites will be considered. 

For these calculations the functional PBE was use, along with OTFG ultrasoft pseudopotentials, the Koelling-Hamon relativistic treatment, a k-point grid of 11x11x1 and a energy cutoff of 600eV. The dipole correction was not used for these calculations as the corrections were smaller that the level of precision guaranteed by the k-point grid and energy cut selection.

Pseudo atomic configuration Cu: 3d10 4s1

Pseudo atomic configuration H: 1s1

Optimization

To calculate this first the hydrogen atom was place on either the hcp or fcc lattice site of a Cu(111) surface and then the geometry of the system was optimized holding the bottom three layers fixed and allowing the top layer of the Cu surface and the hydrogen atom move freely. For the convergence, the energy convergence was set to 10-5 eV/atom, the max force was set to 0.03 eV/A, the max stress was set to 0.05 GPa, and the max displacement was set to 0.001A. The results of the optimization can be seen in Figures 1 and 2.

Figure 1: Optimized structure for hydrogen adsorbed to fcc site.

Figure 2: Optimized structure for hydrogen adsorbed to hcp site.

Construction of the Hessian Matrix

In this work, only the vibrations of the hydrogen are wanted, so an approximation that can be used is to only allow the hydrogen atom to move while constructing the Hessian. This will simplify the Hessian to a 3×3 matrix. From the following equation the Hessian matrix can be calculated.

H_{ij} = \bigg(\frac{\partial^2 E}{\partial x_i \partial x_j}\bigg) = \frac{E(\delta x_i, \delta x_j) - 2E(x_0)+ E(-\delta x_i, -\delta x_j)}{\delta x_i \delta x_j}

Using a displacement of 0.1A the Hessian matrix was computed and then divided by the mass of a hydrogen atom to get the mass-weighted Hessian matrix whose matrix elements are given by the following,

A_{ij} = H_{ij}/m_i

Then the eigenvalues and eigenvectors of the mass-weighted Hessian, A,  were calculated using the following relation.

\vec{A}\vec{e} = \lambda \vec{e}

Once the eigenvalues, \lambda_i, where found the normal mode frequencies, \nu_i, were calculated using the following relation.

\nu_i = \frac{\sqrt{\lambda_i}}{2\pi}

The results of these calculations for the hydrogen at the fcc site and hcp site are given in Tables 1 and 2 respectively.

Table 1: Normal mode frequency and Eigenvector for H atom is fcc site.
Eigenvectors()
Normal ModeFrequency (cm-1)xyz
1742.7-0.1880.769-0.611
2375.90.787-0.254-0.562
3366.7-0.588-0.586-0.558

Table 2: Normal mode frequency and Eigenvector for H atom is hcp site.

Eigenvectors()
Normal ModeFrequency (cm-1)xyz
1740.9-0.182 -0.763-0.620
2372.7 0.7880.264-0.557
3363.2-0.5890.590-0.553

Zero-point Energy Correction

The zero point energy for a quantum mechanical system is not that of the classical system, the corrected form is as the following,

E = E_0 +\sum_i \frac{h \nu_i}{2}

This correction is 0.0921eV for the the system with the hydrogen at the fcc site and 0.0915eV for the system with the hydrogen at the hcp site.

Conclusion

The difference in the energy of the hydrogen adsorbing to the fcc site and the hydrogen adsorbing to the hcp site is 0.0154eV. Taking into account the zero-point energy correction the difference shrinks to 0.0149eV

Zero point energy of H on Cu(111) hollow sites


This project aims to find zero point vibrational energy (ZPVE) for H atom adsorbed at Cu(111) hollow site through calculating Hessian matrix, which was built through second order finite difference approximation.

H adsorption optimization
The geometry optimization was performed in Material Studio with CASTEP Calculation Package [1]. The functional of Perdew Burke, and Ernzerhoff was employed [2]. Cu bulk was first optimized with medium SCF tolerance. Energy cutoff of 480 eV and 9 X 9 X 9 K points sampling were used. The pseudopotential was solved through Koelling-Harmon treatment. After Cu optimization, 111 surface was cut from bulk and H atom was place at hollow hcp site. K points of 11 X 11 X 1 was applied for slab optimization. Figure 1 shows optimized configuration, in which four layers of Cu was chosen and 20 angstrom vacuum was built above the surface to ensure no interaction between super cells.


Figure 1. Optimized H adsorption at Cu (111) hollow hcp site


Figure 2. Optimized H adsorption at Cu (111) hollow fcc site

Energy of structure with 0.05 angstrom H displacement at hcp site

deltaxdeltaydeltazE(eV)
0
00-5057.904815339
-0.0500-5057.900058684
0.0500-5057.900254389
0-0.050-5057.899918579
00.050-5057.900491436
00-0.05
-5057.902027622
000.05-5057.897504906
0.050.050-5057.895183611
0.05
-0.05
0-5057.894736028
-0.050.050-5057.896567514
-0.05-0.050-5057.895964697
0.0500.05-5057.893632073
0.050-0.05-5057.896676147
-0.0500.05-5057.893493162
-0.050-0.05-5057.896416514
00.050.05-5057.893875845
00.05-0.05-5057.896914666
0-0.050.05-5057.893352868
0-0.05-0.05-5057.896277175

Energy of structure with 0.05 angstrom H displacement at fcc site

deltaxdeltaydeltazE(eV)
000-5057.908254387
-0.0500-5057.903846876
0.0500-5057.902812323
0-0.050-5057.905355245
00.050-5057.901336034
00-0.05-5057.903434736
000.05-5057.900649874
0.050.050-5057.896923968
0.05-0.050-5057.900397672
-0.050.050-5057.895898678
-0.05-0.050-5057.900587944
0.0500.05-5057.896058663
0.050-0.05-5057.897022898
-0.0500.05-5057.896937081
-0.050-0.05-5057.898270256
00.050.05-5057.894816130
00.05-0.05-5057.895334570
0-0.050.05-5057.898164719
0-0.05-0.05-5057.900003390

Then finite difference approximation was used to calculate Hessian matrix elements. The diagonal element can be obtained as equation (1) and off-diagonal element can be obtained from equation (2).
\begin{equation} H_{ii} \cong \frac{E(\delta x_{i})+E(-\delta x_{i})-2E_{0}}{\delta x_{i}^2} \end{equation}
\begin{equation} H_{ij} \cong \frac{E(\delta x_{i}, \delta x_{j})+E(-\delta x_{i},-\delta x_{j})-E(\delta x_{i}, -\delta x_{j})-E(-\delta x_{i}, \delta x_{j})}{4\delta x_{i}\delta x_{j}} \end{equation}

After calculating eigenvalues of Hessian matrix, we can calculate vibrational frequency and corresponding zero point energies, which is shown in following table.

siteZero point energy (eV)
hcp0.1893
fcc0.2001

Conclusions
In this project, we applied harmonic approximation and use finite difference approximation to calculate vibrational frequency of H atom adsorbed at the Cu(111) hollow sites and further zero point energy. According to our calculation, the zero point energy of H on hcp site is 0.1893 eV and H on fcc site is 0.2001 eV. The difference of the zero point energies between H on fcc site and hcp site are less 0.011 eV.

Reference
[1] “First principle methods using CASTEP” Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005)
[2] Perdew, J. P; Burke, K; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865-3868

Zero Point Energy Correction of H adsorbed on Cu(111)


1. Problem Description

In this project, the zero point energy (ZPE) correction of one hydrogen atom adsorbed at two different hollow sites of Cu(111) was studied.

In order to get the zero point energy, we firstly initiated the geometry optimization to find the optimal conformation of H adsorbed on Cu(111) as a stating point, by putting the H atom at the fcc-hollow site and hcp-hollow site above the surface for about 2.0 Å. Then, to get the weighted Hessian matrix required to calculate the vibrational frequencies, the H atom was displaced with a step size of 0.05 Å or 0.10 Å along the x, y, and z direction and the corresponding energies were calculated. Different displacements were tested to confirm the correct ZPE. The weighted Hessian matrix consists of the second order derivatives of energies with respect to the displacement along different directions, multiplying the weighted factor. After that, the vibrational frequencies could be found by diagonalizing the Hessian matrix and thus, the ZPE was calculated to correct the binding energy of H at different sites.

 

2. Computational Method

The surface was cleaved from bulk Cu using Materials Studio along the <111> direction with three layers included. In order to minimize the interactions between periodic images of H atoms, we chose to construct a (2×2) supercell including 10 Å of vacuum space as shown in Figure 1. H atom was initially located at two different hollows site above the surface. During the geometry optimization, the upmost layer along with the H atoms were allowed to move while the bottom two layers were totally fixed at the bulk position. After the geometry optimization, the single point calculations were performed to calculate the energy with different displacement of the H atom.

                                 

Figure 1. Top-down view and side view of the Cu(111)-2×2 supercell

 

All the calculations were carried out using CASTEP with Materials Studio and the GGA-PBE was chosen to be the exchange-correlation functional. The pseudopotential was solved with the Koelling-Harmon atomic solver and the valence electrons were 3d10 4s1 and 1s1 for Cu and H, respectively. The energy cutoff of 408.2 eV and the k-point sampling of 6x6x2 were used. The energy cutoff and k-point sampling were tested to give a convergence within 0.02 eV. The results of convergence tests are summarized below in Table 1. The convergence criteria for force and energy were 0.02 eV/ Å and 10e-6 eV, respectively. The self-consistent dipole correction along the Z-axis was also introduced to better describe the adsorption of the H atom.

Energy Cutoff (eV) 400 408.2 410 420
Energy (eV) -20182.9745 -20183.1051 -20183.1281 -20183.1374
K-point 5x5x2 6x6x2 7x7x2 8x8x2
Energy (eV) -20183.0103 -20183.1051 -20183.0376 -20183.1224
Vacuum Space (Å) 10 11 12 13
Energy (eV) -20167.1164 -20167.1163 -20167.1157 -20167.1148

Table 1: Summary of convergence tests

 

3. ZPE at FCC-Hollow Site

As the only atom displaced is H in our system, there weighted Hessian matrix can be easily achieved by multiplying a factor to the Hessian matrix. In order to construct the Hessian matrix for vibrational frequency, the H atom should be displaced along the x, y, and z direction with an appropriate step size to calculate the second derivatives of energies with respect the movements.

For the diagonal elements in the Hessian matrix, the second derivative was calculated using three points as shown below in equation 1, while for the off-diagonal elements, it was calculated using four points as in equation 2.

\begin{equation} H_{ii} \cong \frac{E(\delta x_{i})+E(-\delta x_{i})-2E_{0}}{\delta x_{i}^2} \end{equation}
\begin{equation} H_{ij} \cong \frac{E(\delta x_{i}, \delta x_{j})+E(-\delta x_{i},-\delta x_{j})-E(\delta x_{i}, -\delta x_{j})-E(-\delta x_{i}, \delta x_{j})}{4\delta x_{i}\delta x_{j}} \end{equation}

 

Here,  dxi represents the displacement in the i direction (can be along x, y, and z direction) and  E0 is the minima energy.

In this project, the step size of 0.05 Å and 0.10 Å were both applied in order to ensure that the displacements were appropriately chosen. Due to the symmetry, there were in total 18 sets of energies of different displacements required and the results for H atom at fcc-hollow site with different step sizes are summarized below in Table 2.

Displacement Direction Energy (eV)
X Y Z Step size = 0.05 Å Step size = 0.10 Å
0 0 0 -20183.1051 -20183.1051
1 0 0 -20183.1017 -20183.0910
-1 0 0 -20183.1022 -20183.0943
1 1 0 -20183.0991 -20183.0829
-1 -1 0 -20183.0986 -20183.0773
1 -1 0 -20183.0995 -20183.0841
-1 1 0 -20183.0987 -20183.0779
1 0 1 -20183.0970 -20183.0760
-1 0 -1 -20183.0970 -20183.0704
1 0 -1 -20183.0962 -20183.0647
-1 0 1 -20183.0972 -20183.0777
0 1 0 -20183.1020 -20183.0929
0 -1 0 -20183.1022 -20183.0933
0 1 1 -20183.0970 -20183.0767
0 -1 -1 -20183.0967 -20183.0678
0 1 -1 -20183.0967 -20183.0857
0 -1 1 -20183.0974 -20183.0684
0 0 1 -20183.0996 -20183.0848
0 0 -1 -20183.1004 -20183.0776

Table 2: Summary of energies with different displacements

 

According to the equations above, we were able to calculate the elements in the Hessian matrices as well as the eigenvalues these matrices. Thus, the zero point energy can be calculated using the equation 3 and 4 as show below.

\begin{equation} \nu_{i} = \frac{\sqrt{\lambda_{i}}}{2\pi} \end{equation}
\begin{equation} E_{ZPE} = \frac{h\nu_{i}}{2} \end{equation}

 

Here,  vi is the vibrational frequency and  lambda is the eigenvalue of the Hessian matrix and h is the Planck’s constant.

For the displacements with step size of 0.05 Å, the weighted Hessian matrix is

\begin{bmatrix} 2.412\times10^{28} & 4.786\times10^{26} & -5.744\times10^{26} \\ 4.786\times10^{26} & 2.297\times10^{28} & 3.829\times10^{26} \\ -5.744\times10^{26} & 3.829\times10^{26} & 3.906\times10^{28}\end{bmatrix}

The eigenvalues are 3.90907327e+28, 2.42876582e+28, and 2.27843185e+28. Corresponding zero point energies are 0.130 eV, 0.102 eV, and 0.099 eV.

For the displacements with step size of 0.10 Å,  the weighted Hessian matrix is

\begin{bmatrix} 2.383\times10^{28} & 4.308\times10^{26} & -9.573\times10^{26} \\ 4.308\times10^{26} & 2.297\times10^{28} & 3.590\times10^{26} \\ -9.573\times10^{26} & 3.590\times10^{26} & 3.800\times10^{28}\end{bmatrix}

The eigenvalues are 3.80789093e+28, 2.39781505e+28, and 2.27653408e+28. Corresponding zero point energies are 0.128 eV, 0.102 eV, and 0.099 eV.

We can see that these two different choices bring us similar ZPE with differences less than 2%, which means that our choice of step size of the displacement is appropriate. To sum up, based on the results of 0.10 Å, the total ZPE for H atom at the fcc-hollow site is 0.329 eV. The corresponding eigenvectors are [-0.06640341, -0.91052413, -0.40808872], [0.02181974, -0.41021872, 0.91172611], and [0.99755425, -0.05163733, -0.04710732].

 

4. ZPE at HCP-Hollow Site

Similarly, we carried out the calculations by initially locating the H atom at the hcp-hollow site. Accordingly, the results for H atom at hcp-hollow site with different step sizes are summarized below in Table 3.

Displacement Direction Energy (eV)
X Y Z Step size = 0.05 Å Step size = 0.10 Å
0 0 0 -20183.0973 -20183.0973
1 0 0 -20183.0948 -20183.0874
-1 0 0 -20183.0937 -20183.0830
1 1 0 -20183.0913 -20183.0708
-1 -1 0 -20183.0912 -20183.0753
1 -1 0 -20183.0915 -20183.0718
-1 1 0 -20183.0916 -20183.0766
1 0 1 -20183.0897 -20183.0703
-1 0 -1 -20183.0881 -20183.0560
1 0 -1 -20183.0895 -20183.0631
-1 0 1 -20183.0889 -20183.0678
0 1 0 -20183.0945 -20183.0859
0 -1 0 -20183.0944 -20183.0857
0 1 1 -20183.0895 -20183.0697
0 -1 -1 -20183.0890 -20183.0607
0 1 -1 -20183.0889 -20183.0602
0 -1 1 -20183.0893 -20183.0692
0 0 1 -20183.0917 -20183.0768
0 0 -1 -20183.0925 -20183.0774

Table 3: Summary of energies with different displacements

 

For the displacements with step size of 0.05 Å, the weighted Hessian matrix is

\begin{bmatrix} 2.317\times10^{28} & 5.505\times10^{26} & 1.101\times10^{26} \\ 5.505\times10^{26} & 2.202\times10^{28} & -2.393\times10^{26} \\ 1.101\times10^{26} & -2.393\times10^{26} & 3.868\times10^{28}\end{bmatrix}

The eigenvalues of this Hessian matrix are 3.98502936e+28, 2.35385631e+28, and 2.16250165e+28. The corresponding zero point energies are 0.131 eV, 0.101 eV, and 0.097 eV.

For the displacements with step size of 0.10 Å, the weighted Hessian matrix is

\begin{bmatrix} 2.336\times10^{28} & 5.744\times10^{26} & 5.744\times10^{26} \\ 5.744\times10^{26} & 2.183\times10^{28} & 2.872\times10^{26} \\ 5.744\times10^{26} & 2.872\times10^{26} & 3.983\times10^{28}\end{bmatrix}

The eigenvalues of this Hessian matrix are 3.87576398e+28, 2.33336757e+28, and 2.17737216e+28. The corresponding zero point energies are 0.130 eV, 0.101 eV, and 0.097 eV. 

We can see that these two different choices bring us similar ZPE with differences less than 2%, which means that our choice of step size of the displacement is appropriate. To sum up, based on the results of 0.10 Å, the total ZPE for H atom at the fcc-hollow site is 0.328 eV. The corresponding eigenvectors are [-0.07002192, -0.91680843, 0.39314022], [0.01196003, -0.39485082, -0.9186674 ], and [-0.99747375, 0.05962489, -0.03861327].

 

5. Summary

The energies of Cu slab with the H atom adsorbed at different sites before and after ZPE correction are summarized below in Table 4.

Energy (eV) ZPE(eV) Energy with ZPE (eV)
FCC-Hollow -20183.1051 0.3297 -20182.7754
HCP-Hollow -20183.0973 0.3273 -20182.7701
Difference 0.0078 0.0024 0.0053

Table 4. Summary of ZPE at different sites

 

We can see that ZPE is significant for the H atom adsorbed on Cu(111) but the ZPE at different hollow sites are quite similar. The differences of the total energy with and without ZPE at different sites are both within 0.01 eV, which is very small and the adsorption energy of the H atom at these two sites are almost the same. However, such a small binding energy difference is within the range of convergence error, which means that we cannot decide which site is preferred.

Reference

  1. Clark, Stewart J., et al. “First principles methods using CASTEP.” Zeitschrift für Kristallographie-Crystalline Materials 220.5/6 (2005): 567-570.
  2. Perdew, J.P., K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple. Physical review letters, 1996. 77(18): p. 3865.