Forces and geometry optimization on real-space grids

\label{sec:forces} A function represented on a real-space grid is not invariant under translations as one would expect from a physical system. The potential of an atom sitting on top of a grid point might be slightly different from the potential of the same atom located between points. This implies that a rigid displacement of the system produces an artificial variation of the energy and other properties. If we plot the energy of the atom as a function of this rigid displacement, the energy shows an oscillation that gives this phenomenon the name of the “egg-box effect”.

The egg-box effect is particularly problematic for calculations where the atoms are allowed to move, for example to study the dynamics of the atoms (molecular dynamics) or to find the minimum energy configuration (geometry optimization).

In Octopus we have studied several schemes to control the egg-box effect \cite{Andrade2010thesis}. The first step is to use pseudo-potential filtering to eliminate Fourier components of the potential that cannot be represented on the grid \cite{Tafipolsky_2006}.

Additionally, we have found a formulation for the forces that reduces the spurious effect of the grid on the calculations. One term in the forces is the expectation value of the derivative of the ionic potential with respect to the ionic position \(\vec{R_{\alpha}}\), which can be evaluated as \[\label{eq:forcespot} \vec{F}_{\alpha} = \vec{F}_{\alpha}^{\mathrm{ion-ion}} -\sum_{n} \left< \varphi_{n} \left| \frac{\partial \hat{v}_{\alpha}}{\partial \vec{R}_{\alpha}} \right| \varphi_{n} \right>\ . %\int \mathrm{d}\vec{r} %\varphi^*_{n}(\vec{r})\frac{\partial v_{\alpha}(\vec{r}- %\vec{R}_\alpha)}{\partial \vec{R}_{\alpha}}\varphi_{n}(\vec{r})\ .\] (For simplicity, we consider only local potentials here, but the results are valid for non-local potentials as well.) This term can be rewritten such that it does not include the derivative of the ionic potential \(v_{\alpha}\), but the gradient of the orbitals with respect to the electronic coordinates \cite{hirose2005first}: \[\label{eq:forcesgrad} \vec{F}_\alpha = \vec{F}_\alpha^{\mathrm{ion-ion}}+\sum_{n} \left[ \left< \frac{\partial \varphi_{n}}{\partial r_{\vec{r}}} \left| \hat{v}_{\alpha} \right| \varphi_{n} \right> + \mathrm{cc.} \right]\,. %\int \mathrm{d}\vec{r} %\frac{\partial \varphi^*_{n}(\vec{r})}{\partial \vec{r}}v_{\alpha}(\vec{r}- %\vec{R}_\alpha)\varphi_{n}(\vec{r}) + \mathrm{cc.}\,.\] The first advantage of this formulation is that it is easier to implement than eq. \eqref{eq:forcespot}, as it does not require the derivatives of the potential, which can be quite complex and difficult to code, especially when relativistic corrections are included. However, the main benefit of using eq. \eqref{eq:forcesgrad} is that it is more precise when discretized on a grid, as the orbitals are smoother than the ionic potential. We illustrate this point in Fig. \ref{fig:n2forces}, where the forces obtained with the two methods are compared. While taking the derivative of the atomic potential gives forces with a considerable oscillation due to the grid, using the derivative of the orbitals gives a force that is considerably smoother.