Purpose of Calculation
Geometric optimization of a unit cell is useful for finding the minimum energy configuration of a molecule or solid without needing to manually change atom positions. Through numeric minimization algorithms, the proper lattice parameters and atom positions can be obtained. For molecules, this can also find the minimum energy bond angles and bond lengths. However, the different contexts between solids and molecules can present very different challenges. These calculations will compare the convergence and speed of the geometry optimization calculations of water and hydrogen cyanide (HCN), using two different basis sets: Castep’s [1] plane wave basis set and DMol3’s calculated-on-the-fly localized basis set [2].
Figure 1: Optimized geometry of a water molecule
Figure 2: Optimized geometry of a hydrogen cyanide molecule in the H-C-N isomer
Figure 3: Optimized geometry of a hydrogen cyanide molecule in the H-N-C isomer
Calculation Methodology
The bond angle and bond lengths of water and hydrogen cyanide are calculated using two different basis sets.
All calculations are performed with the Perdew-Burke-Ernzerhof (PBE) [3] exchange-correlation functional, a Generalized Gradient Approximation (GGA) functional. The Koelling-Harmon relativistic treatment was used for atomic solutions [4] in the Castep calculations for the plane wave basis set. 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.
Convergence Calculations
Calculations to check the convergence of the minimum energy output with respect to the energy cutoff, the k-point mesh and the unit cell size were performed for the Castep calculations. It is necessary to check the convergence with respect to the unit cell size because the plane wave basis fills the entire cell. Due to the computational assumption of repeating cells, each unit cell must be large enough so that each molecule does not interact with those in adjacent cells. Convergence tests must be used to idea what is a sufficiently large cell. These tests are shown in figures 4-9 for both water and hydrogen cyanide.
Fig. 4 Convergence of Castep with respect to Cutoff Energy for water
Fig. 5 Convergence of Castep with respect to k-point number for water
Fig. 6 Convergence of Castep with respect to unit cell size for water
Fig. 7 Convergence of Castep with respect to Cutoff Energy for hydrogen cyanide
Fig. 8 Convergence of Castep with respect to k-point number for hydrogen cyanide
Fig. 9 Convergence of Castep with respect to unit cell size for hydrogen cyanide
For DMol3, we do not need to converge the unit cell size. The generated basis sets are highly localized, and we instead test the convergence of the maximum radius of the radial functions, referred to as the cutoff radius. The cutoff radius must be large enough that each atom “sees” the other atoms during all of the geometry optimization, but can still be smaller enough that the basis does not extend as far as the unit cell extends in the Castep calculations. As a result, we check the convergence of energy calculations with respect to cutoff radius and k-point mesh. These convergence tests are shown in figures 10-13 for both water and hydrogen cyanide.
Fig. 10 Convergence of DMol3 with respect to k-point number for water
Fig. 11 Convergence of DMol3 with respect to cutoff radius for water
Fig. 12 Convergence of DMol3 with respect to k-point number for hydrogen cyanide
Fig. 13 Convergence of DMol3 with respect to cutoff radius for hydrogen cyanide
All parameters were converged to within 0.01 eV
Calculation Results
The key metric for comparing the two sets of results is the computational time for both the geometry optimization and the convergence tests. Given an arbitrarily long time, either basis set can give an equally accurate and precise calculation of the bond angles and bond lengths. The important comparison is how quickly the computations can be performed using each basis set to within reasonable agreement with each other.
Both basis sets agree that there is one minimum for the water molecule: a planar bent structure, with bond angle of approximately 104 degrees and a bond length of approximately 0.97 angstroms. Both basis sets also report two isomers of hydrogen cyanide, the H-C-N configuration, and the C-N-H configuration, both with a linear orientation and bond angle of approximately 180 degrees. Disagreement on the exact bond angles and lengths is minor, but does occur. Exact results are included in table one. By better converging the minimum energy, a geometry optimization with a stricter cutoff could be run in order to lessen such disagreement. This, however, is not the purpose of this example case.
| Bond | Castep Bond Length (Angstroms) | Castep Bond Angle (Degrees) | DMol3 Bond Length (Angstroms) | DMol3 Bond Angle (Degrees) |
|---|---|---|---|---|
| H-O | 0.970 | 0.971 | ||
| H-O-H | 103.708 | 103.708 | ||
| H-C (H-C-N Isomer) | 1.075 | 1.076 | ||
| C-N (H-C-N Isomer) | 1.159 | 1.164 | ||
| H-C-N (H-C-N Isomer) | 175.507 | 179.406 | ||
| H-N (H-N-C Isomer) | 1.006 | 1.005 | ||
| N-C (H-N-C Isomer) | 1.176 | 1.181 | ||
| H-N-C (H-N-C Isomer) | 179.970 | 179.959 |
The total computation time for both convergence studies and the geometry optimizations for Castep is 8486 seconds (141.43 minutes). For DMol3, the total time is 102 seconds, providing a similar degree of accuracy. This result is expected, due to the more localized basis set being favorable for describing the more local behavior of a molecule confined to a small portion of a cell. The size of the disparity between the two methods is very noteworthy, however, and shows how helpful it can be to choose the correct basis set for a give task.
References
[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] Delley, J., Fast Calculation of Electrostatics in Crystals and Large Molecules; Phys. Chem. 100, 6107 (1996)
[3] 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)
[4] D D Koelling and B N Harmon 1977 J. Phys. C: Solid State Phys. 10 3107
