1. Problem Description
In this project, we aim to predict the optimal structures of H2O and HCN using DFT. To calculate the energy of an isolated molecule with periodic supercells, the molecule is placed in a large enough cubic supercell to minimize the influence from period images. To perform the geometric optimization, we need to take all the internal degrees of freedom like bond lengths and bond angles into consideration with appropriate stopping criterions for energy and force. We start from various initial structures with all degrees of freedom included to see if they will finally converge to the same structure to ensure that we reach the optimal structures.
The Vienna Ab initio Simulation Package (VASP) is used to perform periodic DFT calculations,1-3 employing the projected augmented-wave (PAW) pseudopotentials and the generalized gradient approximation with exchange-correlation functional by Perdew, Burke, and Ernzerhof (PBE).4-6 The conjugated gradient algorithm is chosen to perform the geometric optimization in VASP.
2. H2O Structure
Before we can perform the geometric optimization of H2O, convergence tests are required firstly with respect to the size of the unit cell as well as the energy cutoff. The test for number of k-points is not required here because all the calculations for such geometric optimization of one molecule in a large box use a 1x1x1 k-point sampling. For both H2O and HCN, all three atoms lie in one plane. In order to do the convergence test, we initially put the O atom at the (0, 0, 0) and two H atoms at (1/L, 0, 0) and (-1/L, 0, 0) according to the experimental results (H-O bond length usually around 1 Å) as shown in Figure 1 (L is the side length of the unit cell).
Figure 1. H2O in a large box (H2O is located at the vertex of the cubic unit cell)
Energy cutoffs ranging from 400.00 to 600.00 eV are used and the side length of the cubic unit cell is chosen from 20.0 to 30.0 Å with an increment of 1 Å. The results of the convergence tests are summarized below in Table 1-2 and Figure 2-3.
Table 1. Convergence tests on energy cutoff for H2O
Table 2. Convergence tests on unit cell size for H2O
Figure 2. Energy vs energy cutoff for H2O Figure 3. Energy vs unit cell size for for H2O
Based on these results, we chose the energy cutoff of 500.00 eV (tolerance of 0.005 eV) and the side length of 25.0 Å for further calculation of the water molecule.
Now, we firstly vary the bond length between H and O atoms while still keeping them linearly arranged on the y-axis. The O atom remains (0, 0, 0) and the two H atoms are located at (a/L, 0, 0) and (-a/L, 0, 0). Here, L=25.0 Å and a is chosen from 0.5 to 1.3 Å with the interval of 0.1 Å. The results are summarized below in Table 3.
Table 3. Results for H2O with different initial bond lengths
From these results, we can see that the energies for the cases with initial bond lengths of 0.5 and 0.6 Å are significantly higher as the distances between atoms in the optimal geometry, which are trapped in a local minimum, are much larger than typical bond lengths as shown in Figure 4. The distances between the O atom and the H atom are 9.790 and 4.809 Å for these two cases, respectively. As two H atoms are initially too close to the O atom, this situation physically represents the enormously compressed bonds, which brings a strong repulsive force pushing the atoms away from each other. The numerical optimization method then takes an initial step that separates the two atoms by a very large distance and they can never go back since the force is much weaker now due to the very large distance. All other cases give us the same optimized structure as shown in Figure 5. The H-O bond length is 0.943 Å and the bond angle qHOH is 180.00°.
Figure 4. Separated H2O with initially short distances Figure 5. Optimized structure for H2O
In order to ensure that the structure above is the optimal structure, we modify the initial structure through changing the bond angle that does not make components of the forces vanished by symmetry. So, with the O atom remaining at (0, 0, 0), the two H atoms are at (a/L, b/L, 0) and (-a/L, b/L, 0). L here is the side length of 25.0 Å, and according to the results above, the a is fixed at 0.9 Å while b is chosen from 0.1, 0.2, and 0.3 Å to have three different initial bond angles. After optimization, we find that these three tests finally converge at the same free energy of -14.217 eV with the same optimized structure of H2O as shown in Figure 6. The H-O bond length is 0.972 Å and the bond angle qHOH is 104.50°. The energy is lower than the previous cases with bond angle of 180.00° and since three different bond angles finally converges to the same structure, this is the optimal structure for H2O. The experimental H-O bond length and bond angle is 0.958 Å and 104.48°, respectively. We see that the DFT-predicted bond length and bond angle are very similar to the experimental results with differences less than 2%.
Figure 6. Optimal structure for H2O
3. HCN Structure
Similarly, the convergence tests are required before we optimize the structure HCN. Due to the similar size of these two molecules, the convergence test on unit cell sizes is no longer required and the side length of 25.0 Å above is used. We similarly put the C atom at the (0, 0, 0), the N atom and the H atom are at (1/L, 0, 0) and (-1/L, 0, 0) as shown in Figure 7 (L is the side length of the unit cell and is fixed at 25.0 Å here).
Figure 7. HCN in a large box (HCN is located at the vertex of the cubic unit cell)
Energy cutoffs ranging from 400.00 to 600.00 eV are tested for HCN. The results are summarized below in Table 4 and Figure 8. Based on these results, we chose the energy cutoff of 500.00 eV using the same tolerance for further calculations of HCN.
Table 4. Convergence tests on energy cutoff for HCN
Figure 8. Energy vs energy cutoff for HCN
Again as we did for the H2O molecule, we vary the N-C and H-C bond lengths while keeping all of them on the y-axis. The C atom remains (0, 0, 0) and the N and H atoms are located at (a/L, 0, 0) and (-a/L, 0, 0). Here, L=25.0 Å and a is chosen from 0.7 to 1.5 Å with the interval of 0.1 Å. The results are summarized below in Table 5.
Table 5. Results for H2O with different initial bond lengths
From these results, we can see that the energies for the cases with initial bond lengths of 0.7, 0.9, and 1.0 Å are significantly higher due to the large distances between atoms as shown in Figure 9. These atoms are pushed apart from each other with the similar reason for those H2O with very short bond lengths (only one H-C bond, the N atom is far away). For the case of 0.8 Å, it even cannot converge. The optimized structures for the rest cases are totally the same as shown in Figure 10. The N-C and H-C bond lengths are 1.161 and 1.075 Å, respectively. The bond angle qHCN is 180.00°.
Figure 9. Separated HCN with initially short distances Figure 10. Optimized structure for HCN
Similarly, we then modify the initial structure by changing the bond angle to make sure that the force components are not canceled out by symmetry. So, with the C atom remaining at (0, 0, 0), the N and H are at (a/L, b/L, 0) and (-a/L, b/L, 0). L here is the side length of 25.0 Å, and according to the results above, the a is fixed at 1.2 Å while b is chosen from 0.1, 0.2, and 0.3 Å to have three different initial bond angles. After optimization, we find that all these three tests finally give the same free energy of -19.690 eV and the same optimized structure of HCN as shown in Figure 11. The N-C and H-C bond lengths are 1.161 and 1.075 Å, respectively. The bond angle qHCN is 179.95°. The energy is almost the same compared to the previous cases with bond angle of 180.00° and since three different bond angles finally converges to the same structure, this is the optimal structure for HCN. The experimental N-C and H-C bond lengths are 1.156 and 1.064 Å, and the bond angle 180.00°. We see that the DFT-predicted bond length and bond angle are very similar to the experimental results with differences less than 2%.
Figure 11. Optimal structure for HCN
4. Conclusion
To conclude, in this project, we find the optimal structures of H2O and HCN. For H2O, the DFT-predicted H-O bond length is 0.972 Å and the bond angle qHOH is 104.50° (the experimental values are 0.958 Å and 104.48). For HCN, the DFT-predicted N-C and H-C bond lengths are 1.161 and 1.075 Å while the bond angle qHCN is 179.95° (the experimental values are 1.156 and 1.064 Å, and 180.00°). We see that the DFT can accurately predict the bond lengths and bond angles for these selected molecules with differences less than 2% compared to experimental results.
Reference
- Kresse, G. and J. Hafner, Ab initio molecular dynamics for liquid metals. Physical Review B, 1993. 47(1): p. 558.
- Kresse, G. and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium. Physical Review B, 1994. 49(20): p. 14251.
- Kresse, G. and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical review B, 1996. 54(16): p. 11169.
- Kresse, G. and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B, 1999. 59(3): p. 1758.
- Blöchl, P.E., Projector augmented-wave method. Physical review B, 1994. 50(24): p. 17953.
- Perdew, J.P., K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple. Physical review letters, 1996. 77(18): p. 3865.
