Accurate ionic forces and geometry optimization in linear-scaling density-functional theory with local orbitals

Linear scaling methods for density-functional theory (DFT) simulations are formulated in terms of localized orbitals in real space, rather than the delocalized eigenstates of conventional approaches. In local-orbital methods, relative to conventional DFT, desirable properties can be lost to some extent, such as the translational invariance of the total energy of a system with respect to small displacements and the smoothness of the potential-energy surface. This has repercussions for calculating accurate ionic forces and geometries. In this work we present results from onetep, our linear scaling method based on localized orbitals in real space. The use of psinc functions for the underlying basis set and on-the-fly optimization of the localized orbitals results in smooth potential-energy surfaces that are consistent with ionic forces calculated using the Hellmann-Feynman theorem. This enables accurate geometry optimization to be performed. Results for surface reconstructions in silicon are presented, along with three example systems demonstrating the performance of a quasi-Newton geometry optimization algorithm: an organic zwitterion, a point defect in an ionic crystal, and a semiconductor nanostructure.

Figure 1

Total energy as calculated by onetep and castep at a fixed geometry as a function of cutoff energy $E_c$ and NGWF truncation radius $R_{\varphi }$, for the CO${}_2$ molecule (left) and an ${\mathrm{H}}_2$O dimer (right). Convergence properties are dominated by the oxygen pseudopotential and so are broadly the same for the two systems. The results for 4.0-Å spheres to 6.0-Å spheres agree to high accuracy, particularly at a cutoff energy of 900 eV or above, indicating convergence at around 4.0 Å.

Figure 2

Convergence properties of a single force component calculated using onetep and castep at a fixed, off-equilibrium geometry as a function of cutoff energy $E_c$ and NGWF truncation radius $R_{\varphi }$, for the CO${}_2$ molecule (left) and an ${\mathrm{H}}_2$O dimer (right). Force shown is the force acting on an O atom along the bond axis for CO${}_2$ (hence a very strong force), and the force on the hydrogen atom not involved in the dimer-dimer $h$ bond (hence a very weak force, near equilibrium). Forces converge with both $R_{\varphi }$ and $E_c$, but care must be taken to accurately converge with respect to both quantities simultaneously.

Figure 3

Convergence properties of the RMS error in the onetep and castep forces at a fixed, off-equilibrium geometry as a function of cutoff energy $E_c$ and NGWF truncation radius $R_{\varphi }$, for the CHAPS molecule (see text). The quantity plotted is the RMS difference between the calculated result for a given $E_c$ and $R_{\varphi }$ and the calculated result for $E_c=1200$ eV in castep, which is taken to be the converged result. We see that with respect to both $R_{\varphi }$ and $E_c$ it is possible to obtain systematic convergence to the plane-wave result.

Figure 4

Total energy as a function of C-O bond-length for a carbon dioxide molecule subjected to a symmetric stretch. castep (plus symbols); onetep (open diamonds). The curve shows a polynomial fit $E(x)$ to the data points. On the scale of this plot the data are indistinguishable. Inset: the Hellmann-Feynman force on each oxygen atom. The curve shows the numerical force $F(x)\equiv -\frac{1}{2}\frac{dE}{dx}$.

Figure 5

Interaction potential of a water dimer. castep (plus symbols); onetep (open diamonds). The curves show polynomial fits $E(x)$ to the data points. Inset: the Hellmann-Feynman force on each water molecule. The curves show the numerical force $F(x)\equiv -\frac{dE}{dx}$ for castep (dashed line) and onetep (solid line).

Figure 6

Schematic of the reconstructed Si(001) surface. Bonds indicated by dashed lines do not lie in the plane of the diagram. The surface atoms (indicated by white circles) pair up to form dimers, that then buckle out of the plane of the surface. For the purposes of this work the bonds have been labeled alphabetically and the buckling angle denoted by $\theta$.

Figure 7

Systems for which geometry optimization was performed with the BFGS algorithm in onetep for illustration of convergence behavior. Left: CHAPS molecule (100 atoms) Center: Al vacancy in $2\times 2\times 1$ supercell of $\alpha -$alumina (119 atoms) Right: H-terminated GaAs nanocrystal (204 atoms).

Figure 8

Convergence behavior of the maximum force $\left|F_{\max }\right|$ on any atom as the BFGS algorithm proceeds for the typical systems shown in Fig. 7. CHAPS Zwitterion (open circles); Al vacancy in alumina (filled circles); GaAs nanocrystal (squares). Dashed line shows convergence threshold of 0.11 eV/Å. Simultaneous convergence of total energy and the maximum displacement were required, with $\left|d{\mathbf{R}}_{\max }\right|=0.002$ Å  for 2 successive iterations was also required. From the starting coordinates (see text) convergence was achieved in 29, 21, and 29 iterations, respectively.

Figure 9

Scaling of total computational time for SCF total-energy calculation for a series of DNA molecules of increasing length (squares), compared to scaling of total computational time for forces calculation (triangles), and the three force components present: Ewald, local potential, and nonlocal potential (diamonds, empty circles, and filled circles, respectively), for the same simulations.