16. Potential calculation

In this section, we introduce details of potential calculations. Recall the Kohn-Sham equation

\[\left[ -\frac{1}{2}(\nabla+i\mathbf{k})^2 + v_{\text{NA}}(\mathbf{r}) + v_{\delta H}(\mathbf{r}) + v_{xc}(\mathbf{r})\right]u_{\mathbf{k}i}(\mathbf{r}) = \epsilon_{\mathbf{k}i}u_{\mathbf{k}i}(\mathbf{r})\]

where

\[-\nabla^2 v_{\delta H}(\mathbf{r}) = 4\pi(\rho(\mathbf{r}) - \rho_{NA}(\mathbf{r}))\]

It is often useful to inspect and plot various contributions to the potential.

First, we need a converged self-consistent calculation. Suppose we are interested in AlP

info.savepath = './results/alp_lcao_scf'
info.calculationType = 'self-consistent'
atom.element = [1 2]
atom.fracxyz = [0 0 0
0.25 0.25 0.25]
domain.latvec = 5.4510*...
[0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0];
domain.lowres = 0.3
units.latvec = 'A'
element(1).species = 'Al'
element(2).species = 'P'
kpoint.gridn = [4,4,4]

Then, we set up a potential calculation as follows

info.savepath = './results/alp_lcao_pot'
info.calculationType = 'potential'
rho.in = './results/alp_lcao_scf'

Let’s inspect the results

>> h5disp('results/alp_lcao_pot.h5','/potential')
HDF5 alp_lcao_pot.h5
Group '/potential'
   Dataset 'vdh'
      Size:  15625x1
      MaxSize:  15625x1
      Datatype:   H5T_IEEE_F64LE (double)
      ChunkSize:  []
      Filters:  none
      FillValue:  0.000000
      Attributes:
            'real':  1.000000
   Dataset 'vna'
      Size:  15625x1
      MaxSize:  15625x1
      Datatype:   H5T_IEEE_F64LE (double)
      ChunkSize:  []
      Filters:  none
      FillValue:  0.000000
      Attributes:
            'real':  1.000000
   Dataset 'vxc'
      Size:  15625x2
      MaxSize:  15625x2
      Datatype:   H5T_IEEE_F64LE (double)
      ChunkSize:  []
      Filters:  none
      FillValue:  0.000000
      Attributes:
            'real':  1.000000

By default, there three potentials are saved: vdh, vna and vxc.

To control which potential is calculated and saved, set option.savePotential to a list of strings specifying which potentials to compute and save. For example,

info.savepath = './results/alp_lcao_pot'
info.calculationType = 'potential'
rho.in = './results/alp_lcao_scf'
option.savePotential = {'vh','vps'}

The available options are the following.

vatom: neutral atom charge electrostatic potential.

\[-\nabla^2 v_{atom}(\mathbf{r}) = 4\pi\rho_{NA}(\mathbf{r})\]
vdh: \(\delta\)-Hartree potential.

\[-\nabla^2 v_{\delta H}(\mathbf{r}) = 4\pi(\rho(\mathbf{r}) - \rho_{NA}(\mathbf{r}))\]
veff: effective potential.

\[\begin{split}v_{\text{eff}} = v_{\text{NA}}(\mathbf{r}) + v_{\delta H}(\mathbf{r}) + v_{xc}(\mathbf{r})\\ v_{\text{eff}} = v_{\text{PS}}(\mathbf{r}) + v_{H}(\mathbf{r}) + v_{xc}(\mathbf{r})\end{split}\]
vh: Hartree potential.

\[-\nabla^2 v_{H}(\mathbf{r}) = 4\pi\rho(\mathbf{r})\]
vna: neutral atom potential.

\[v_{\text{NA}} = v_{\text{PS}}(\mathbf{r}) + v_{\text{atom}}(\mathbf{r})\]
vps: pseudo-ion potential.

\[v_{\text{PS}} = \sum\limits_{\mathbf{T}I} v_{loc,I}(\mathbf{r}-\mathbf{R}_I-\mathbf{T})\]
vxc: exchange-correlation potential.

\[v_{xc}^\sigma = \frac{\delta E_{xc}[\rho^\uparrow,\rho^\downarrow]}{\delta\rho^\sigma}\]