Author Archives: Mark Bronson

ab initio Phase Diagram for LiH and Li

Leave a reply

Introduction

In this project the phase diagram for LiH and Li will constructed.  To do this three single point energy calculations must be preformed, one for the LiH crystal, the Li crystal and H_2 molecule. Before the single point energies could be calculated, the geometry of the crystals and molecule were optimized.

The calculations were performed with CASTEP [2] and used the GGA XC-functional PBE [1], OTFG ultrasoft pseudopotentials, the Koelling-Hamon relativistic treatment, a k-point grid of 11x11x11 for the  crystals and 2x2x2 for H_2 for the and a energy cutoff of 600eV. The SCF tolerance was set to 10^{-7} eV/atom with a convergence window of three steps, For the geometry optimizations the energy convergence was set to 10^{-6} eV/atom, the max force was set to 0.01 eV/A, the max stress was set to 0.02 GPa, and the max displacement was set to 0.0005A. For the two crystals a full cell optimization was used with hard comprehensibility.

The k-point grid used converged the energy of the Li and LiH crystals to within 0.005eV and energy cutoff used converged the energy of the Li and LiH crystals to within 0.02eV.

Pseudo atomic configuration Li: 1s2 2s1

Pseudo atomic configuration H: 1s1

Optimization

The geometries were optimized using the settings above and the results are shown in Figure 1, 2 and 3

Figure 1: Primitive bcc lattice of Li, a = 2.983Å

Figure 2: Primitive NaCl type lattice of LiH, a = 2.835Å

Figure 3: hydrogen molecule with bond length 0.752Å in a 15Å sided cube of vacuum.

The Grand Partition Function

To find the equilibrium point between LiH and Li the grand partition functions,\Omega [3] of the two were set to be equal and the chemical potential of atomic hydrogen, \mu_H, was solved for

\displaystyle \Omega_{Li} = E_{Li} - \Omega^M_S      (eq 1)

\displaystyle \Omega_{LiH} = E_{LiH} - \mu _{H} N_{H,LiH} - \Omega^{M}_{S}      (eq 2)

where E_{Li} and E_{LIH} are the internal energies of Li and LiH crystals respectively,  N_{H,LIH} is the number of hydrogen atoms in the LiH crystal and \Omega^{M}_{S} is an additive constant that is the same for all materials.  When those two equations are combined the equilibrium chemical potential of atomic hydrogen can be expressed as,

\displaystyle \mu_H = \frac{ E_{LiH} - E_{Li} } { N_{H,LiH} }      (eq 3)

In addition the, chemical potential of atomic hydrogen can be expressed in terms of the chemical potential of molecular hydrogen \mu_{H_2} ,

\mu_H = \frac{1} {2} \mu_{H_2}      (eq 4)

Assuming that the hydrogen behaves like an ideal gas one can write the chemical potential of molecular hydrogen as,

\mu_{H_2} = E^{Total}_{H_2} + \hat{\mu}_{H_2}(T, p^o) +kTln(p/p^o)     (eq 5)

where E^{Total}_{H_2} is the total energy of an isolated H_2 molecule at T=0K, \hat{\mu}_{H_2}(T, p^o) is the difference in the chemical potential of H_2 between T=0K and the temperature of interest at the reference pressure, p^o is the reference pressure, p is the pressure of the system, and T is the temperature of the system.

Chemical Potential Difference

The values for the difference in the chemical potential of H_2 between T=0K and the temperature of interest at the reference pressure can be evaluated using data tabulated in the NIST-JANAF Thermochemical Tables and the following equation

\hat{\mu}_{H_2}(T, p^o) = [H^o (T) - H^o (T_r) ] - TS (T) - [H^o (0) - H^o (T_r)]      (eq 6)

Where T_r = 298.15K;  [H^o (T) - H^o (T_r) ], the difference in the enthalpy at temperature T and the enthalpy at the reference temperature;  and S(T), the entropy at temperature T, are values given in the reference table. Selected values are listed in Table 1, along with the calculated value of \hat{\mu}_{H_2}(T, p^o) .

Table 1: Values from NIST-JANAF Thermochemical Tables and the calculated difference in chemical potential.

Results

Single point calculations were performed on the optimized geometries for the Li crystal, the LiH crystal and the molecular hydrogen, the results of which are listed in Table 2.  Using equations 3 and 4 the value of the chemical potential of atomic hydrogen and molecular hydrogen were found, the resulting values can be found in Table 2.

Table 2: Results of single point energy calculations for the three systems.

 

As we know the value of \mu_{H_2}, E^{Total}_{H_2} and \hat{\mu}_ {H_2}(T, p^o) equation 5 can be rewritten to solve for p.

\displaystyle p(T) = p^0 e^{ \bigg( \frac{\mu_{H_2} - E^{Total}_{H_2} - \hat{\mu}_{H_2}(T, p^o) } { kT } \bigg) }     (eq 7)

Since the relationship of temperature to pressure is known at the equilibrium point between LiH and Li, the phase diagram can be constructed and is presented in Figure 4.

Figure 4: Phase diagram of Li and LiH

References

1. Perdew, J. P.; Burke, K.; Ernzerhof, M. Physical Review Letters 1996, 77 (18), 3865–3868.

2. Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. I. J.; Refson, K.; Payne, M. C. Zeitschrift für Kristallographie – Crystalline Materials 2005, 220 (5/6).

3. Sholl, D.; Steckel, J. A.; Sholl. Density Functional Theory: a Practical Introduction; Wiley: Somerset, 2011.

Lattice Parameter Prediction of Hafnium

Leave a reply

For this problem we want to predict the lattice parameters of Hafnium and to compare them with the experimental values, given that Hafnium is observed to be a hcp metal, see Figure 1, with a ratio of c to a of 1.58.

Figure 1: Structure of Hf lattice with a =3.1 and c/a=1.58

Methods

To find the optimal values for the lattice parameters, energy calculations were be preformed at differing values of lattice parameter a and c. To turn this into a pseudo one parameter optimization the ratio of c to a was be held fixed while a was varied, then this was repeated for multiple c/a ratios.

The exchange-correlation functional used was PBE, the pseudopotential used was OTFG ultrasoft and the relativistic treatment was that that of Koelling-Harmon.

Convergence tests with respect to the number of k-points and the energy cutoff energy were preformed and the results of which can be found in the Appendix. Using those results the energy cutoff was set to 435.4eV, and the k-point grid was set to 10x10x6. Also, the SCF tolerance was set to 1.0e-6eV/atom with a convergence window of three steps, all other parameters where kept to the quality fine preset for CASTEP.

Material Studio was used to create the crystal of Hafnium as a hcp metal of group P63/MMC and change the lattice parameters, while CASTEP was used to perform the calculations.

Atomic electron configuration for Hf is:
1s2 2s2 2p6 3s2 3p6 3d10 4s2 4p6 4d10 4f14 5s2 5p6 5d2 6s2

And the pseudo atomic electron configuration for Hf used was:
4f14 5s2 5p6 5d2 6s2

Optimization at c/a=1.58

As a starting point, the c/a ratio of was set to 1.58 and the lattice parameter a was set to the range from 2.5 Å to 4.0 Å in increments of 0.1 Å to determine the range of a where the optimal value for a should lie. The results of those calculations can be found in Table 1 and Figure 2.

Table 1: Results of the calculations for c/a=1.58

Figure 2: Cohesive Energy per atom versus Lattice Parameter a, for c/a =1.58

From Table 1 and Figure 2 it can be seen that the minimum cohesive energy occurs at a = 3.2 Å, a volume per atom of 22.419 Å3.

Optimization of c/a

To optimize the ratio of c to a, energy calculations were preformed for a c to a ratio of 1.56 and 1.60 for lattice parameter a ranging from 3.0 Å to 3.4 Å, with the resulting cohesive energy versus volume plots are shown in Figure 3.

 

Figure 3: Cohesive energy per atom versus Volume per atom plots for the different c/a ratios tested.

From Figure 3 we can see that the cohesive energy was at it’s minimum with a c to a ratio of 1.58.

To obtain the optimum value for a, and by extension c, a quadratic function was fit to the plot of the cohesive energy versus lattice parameter a for the c to a ratio of 1.58. The minimum point of the function was used to find the optimal value of a.

Figure 4: Plot of Cohesive Energy per atom vs lattice parameter a with a c/a=1.58.

From the quadratic fit the optimal values of a and c were found to be 3.20 Å and 5.06 Å, respectively.

From experiment the lattice parameters a and c of Hafnium were found to be 3.1964 Å and 5.0511 Å, respectively.  This corresponds to a relative error for a and c of  0.2% and 0.2%, respectively.

Appendix

Validation of the Cutoff Energy

To validate that the cutoff energy was low enough to obtain the level of precision we want from the calculations we can look at the convergence of the free energy with respect to the cutoff energy. While this calculation was done at the start of this project, it was redone for an a of 3.2Å  and a c/a ratio of 1.58 as it will serve as a stronger validation for the results presented.

Figure 5: Free Energy vs. Cutoff Energy for a=3.2, c/a=1.58 and k-points: 10x10x6

From Figure 5, one can see that the free energy for the cutoff used, 435.4eV, was within 0.05eV of the converged value.

Validation of k-points

To validate that the number of k-points was high enough for the level of precision we want from the calculations, we can look at the convergence of the free energy with respect to the number of k-points.

Figure 6: Free Energy vs. Number of k-points, a=3.2Å, c/a=1.58 and Ecut=480eV

Figure 6 shows that the k-point values used, 10x10x6, was within 0.001eV of the converged value.

Cohesive Energy

The DFT program used, CASTEL does not report the cohesive energy, it reports the energy of the free atom and then the free energy of the crystal.  To correct for that we define the cohesive energy as the free energy of the crystal divided by the number of atoms in the crystal’s cell minus the energy of the free atom.

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

Leave a reply

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

Formation Energy of bcc CuPd Alloy

Leave a reply

For this problem we want to find the formation energy of bcc CuPd alloy from fcc Cu and fcc Pd. Then use that to justify why CuPd is the favored low temperature crystal structure of Pd and Cu when they are mixed with this stoichiometry.

Methods

The optimal lattice parameter for fcc Cu, fcc Pd and bcc CuPd was be found, then using the optimal lattice parameter the cohesive energy of each crystal was calculated.

The exchange-correlation functional used was PBE, the pseudopotential used was OTFG ultrasoft and the relativistic treatment was that that of Koelling-Harmon. The pseudo atomic configuration of the pseudo potential was 3d10 4s1 for the copper atom and 4s2 4p6 4d10 5s0.5 for the palladium atom.

Convergence tests with respect to the number of k-points and the energy cutoff energy were preformed for each of the three crystals and the results of which can be found in the following section. To get a convergence of 0.01eV in the cohesive energy a energy cutoff of 4000eV and a k-point grid of 10x10x10 were used. Also, the SCF tolerance was set to 1.0e-6eV/atom with a convergence window of three steps, all other parameters where kept to the quality fine preset for CASTEP.

For the calculation the space group for Cu and Pd was Fm-3m and Pm-3m for the CuPd alloy.

Convergence

k-points

For each crystal the k-point grid was varied and the energy calculated, the results shown in Figures 1,2 and 3.

Figure 1: fcc Cu Cohesive Energy vs. Number of k-points, a=3.615Å and Ecut=408.2eV

Figure 2: fcc Pd Cohesive Energy vs. Number of k-points, a=3.891Å and Ecut=560.6eV

Figure 3: bcc CuPc Cohesive Energy vs. Number of k-points, a=3.012Å and Ecut=560.6eV

From these three graphs one can see that a k-point grid of 10x10x10 will converge the cohesive energy of all the crystals to at least 0.01ev

Energy Cutoff

For each crystal the energy cutoff was varied and the energy calculated, the results shown in Figures 4, 5 and 6. The k-point grid was selected to be 11x11x11 for each of the crystals to ensure that the convergence of the k-point grid did not effect the convergence of the energy cutoff.

Figure 4: fcc Cu Cohesive Energy vs. Energy Cutoff, a=3.615Å and k-points 11x11x11

Figure 5: fcc Pd Cohesive Energy vs. Energy Cutoff, a=3.891Å and k-points 11x11x11

Figure 6: bcc CuPd Cohesive Energy vs. Energy Cutoff, a=3.012Å and k-points 11x11x11

From these three graphs we can see that an energy cutoff of 5000eV converged the cohesive energies to at least 0.02eV, the limiting case being the CuPd alloy in Figure 6.

Optimization

For each of the three crystals, the lattice parameter was optimized in CASTEP, using the exchange correlation functional PBE, a k-point grid of 10x10x10, a cutoff energy of 5000eV, the BFGS algorithm for geometry optimization with full cell optimization. For the convergence tolerance for the energy was 1e-5 eV per atom, for the max force was 0.03 eV per Å, for the max stress was 0.05 GPa, and for the max displacement was 0.001 Å.

For both the fcc Cu and fcc Pd crystal CASTEP changed the unit cell to reduce the computational cost, in the new cell a, b and c are still equal, but the angles α, β, γ, changed from 90° to 60°.

The results of the optimizations are listed in Table 1.

Table 1: Lattice parameters after optimization of the crystals unit cell
CrystalLattice Parameter
fcc Cu2.566
fcc Pd2.786
bcc CuPd3.013

Results

For each of the three crystals the energy of the crystal was taken from the last cycle of the geometry optimization and can be found in Table 2.

Table 2: Cohesive Energies of optimized crystal structures, the Final Free Energy is per unit cell
Final free energy (E-TS) (eV)Energy of Free Cu atom (eV)Energy of Free Pd atom (eV)Cohesive Energy (eV)
Cu-1680.93-1677.14-3.79
Pd-3493.26-3490.46-2.79
CuPd-5174.43-1677.14-3490.46-6.83

From the cohesive energies in Table 2, the formation energy of the bcc CuPc alloy was found to be -0.24eV.

From the literature we know that the bcc CuPd alloy is the favored low temperature crystal structure of Pd and Cu when they are mixed with this stoichiometry, this would mean that the CuPd alloy is more thermodynamically  stable than the two separate mono-atomic crystals, leading to a formation energy that is negative. This matches with the results of the calculations.