• Energy
• Densities
• Potentials

Energy

Energy is computed from the formula:

E=\int \mathcal{E}_{\textrm{edf}}(n,\nu,\ldots)\,d^3r+\sum_{\sigma}\int V_{\sigma}^{\textrm{(ext)}}(r)n_{\sigma}(r)\,d^3r\\-\frac{1}{2}\int\left(\Delta^{\textrm{(ext)}}(r)\nu^*(r)+\textrm{h.c.}\right)d^3r-\sum_{\sigma}\int \vec{v}_{\sigma}^{\textrm{(ext)}}(r)\cdot\vec{j}_{\sigma}(r)\,d^3r

The intrinsic energy is assumed to have a generic structure:

E_{\textrm{edf}} = \int \mathcal{E}_{\textrm{edf}}\,d^3r = \int\left(\mathcal{E}_{\textrm{kin}} + \mathcal{E}_{\textrm{pot}} + \mathcal{E}_{\textrm{pair}} + \mathcal{E}_{\textrm{curr}}\right)d^3r

where:

• E_kin:

E_{\textrm{kin}} = \int\mathcal{E}_{\textrm{kin}}\,d^3r = \int \sum_{\sigma=\{\uparrow,\downarrow\}}\frac{\alpha_{\sigma}}{2}\left(\tau_{\sigma}-\frac{\vec{j}_{\sigma}^2}{n_{\sigma}}\right)d^3r

• E_pot:

E_{\textrm{pot}} = \int\mathcal{E}_{\textrm{pot}}\,d^3r =\int D(n_{\uparrow},n_{\downarrow}, \nu, \ldots)\,d^3r

• E_pair:

E_{\textrm{pair}} = \int\mathcal{E}_{\textrm{pair}}\,d^3r =\int g(n_{\uparrow},n_{\downarrow})\nu^{\dagger}\nu\,d^3r

• E_curr:

E_{\textrm{curr}} = \int\mathcal{E}_{\textrm{curr}}\,d^3r =\int \sum_{\sigma=\{\uparrow,\downarrow\}}\frac{\vec{j}_{\sigma}^2}{2n_{\sigma}}\,d^3r

Contributions to energies from the external potentials are:

• E_potext:

E_{\textrm{pot.ext}} = \int \sum_{\sigma=\{\uparrow,\downarrow\}} V_{\sigma}^{\textrm{(ext)}}(r)n_{\sigma}(r)\,d^3r

• E_pairext:

E_{\textrm{pair.ext}} = -\frac{1}{2}\int\left(\Delta^{\textrm{(ext)}}(r)\nu^*(r)+\textrm{h.c.}\right)d^3r

• E_velext:

E_{\textrm{vel.ext}} = -\int \sum_{\sigma=\{\uparrow,\downarrow\}}\vec{v}_{\sigma}^{\textrm{(ext)}}(r)\cdot\vec{j}_{\sigma}(r)\,d^3r

The total energy E_tot is computed as the sum of all these contributions.

Densities

Densities are computed according to formulas:

• nu:
  \nu(r) = \frac{1}{2}\sum_{|E_n|<E_c} u_{n,\uparrow}(r)v_{n,\downarrow}^{*}(r)( f_{\beta}(-E_n)-f_{\beta}(E_n) )
• rho_a:
  n_{\uparrow}(r)= \sum_{|E_n|<E_c}|u_{n,\uparrow}(r)|^2 f_{\beta}(E_n)
• rho_b:
  n_{\downarrow}(r) = \sum_{|E_n|<E_c}|v_{n,\downarrow}(r)|^2 f_{\beta}(-E_n)
• tau_a:
  \tau_{\uparrow}(r) = \sum_{|E_n|<E_c}|\nabla u_{n,\uparrow}(r)|^2 f_{\beta}(E_n)
• tau_b:
  \tau_{\downarrow}(r) = \sum_{|E_n|<E_c}|\nabla v_{n,\downarrow}(r)|^2 f_{\beta}(-E_n)
• j_a_x, j_a_y, j_a_z:
  \vec{j}_{\uparrow}(r) = -\sum_{|E_n|<E_c} \textrm{Im}[u_{n,\uparrow}(r)\nabla u_{n,\uparrow}^*(r)]f_{\beta}(E_n)
• j_b_x, j_b_y, j_b_z:
  \vec{j}_{\downarrow}(r) = \sum_{|E_n|<E_c} \textrm{Im}[v_{n,\downarrow}(r)\nabla v_{n,\downarrow}^*(r)]f_{\beta}(-E_n)

In these formulas E_{n} denotes quasi-particle energy, and E_c is the energy cut-off scale. Fermi distribution function f_{\beta}(E)=1/(\exp(\beta E)+1) is introduced to model temperature T=1/\beta effects.

Densities are accessible to users through structure wslda_density:

typedef struct
{
    int nx;             
    int ny;             /// for 1d code it is set to 1
    int nz;             /// for 1d and 2d code it is set to 1
    int datadim;        /// dimensionality of data
    int blocklength;    /// number of elements in potential array = nx*ny*nz
    
    // densities
    double complex *nu;
    double *rho_a;
    double *rho_b;
    double *tau_a;
    double *tau_b;
    double *j_a_x;
    double *j_a_y;
    double *j_a_z;
    double *j_b_x;
    double *j_b_y;
    double *j_b_z;  
} wslda_density; 

Potentials

Minimization of the functional with respect to quasiparticle orbitals provides Bogoliubov-de Gennes type equations. Its general form is:

\begin{pmatrix}h_{\uparrow}(r)  &  \Delta(r)+ \Delta_{\textrm{ext}}(r) \\\Delta^*(r)+ \Delta^*_{\textrm{ext}}(r) &  -h^*_{\downarrow}(r)\end{pmatrix} \begin{pmatrix}u_{n\uparrow}(r) \\ v_{n\downarrow}(r)\end{pmatrix}= E_n\begin{pmatrix}u_{n\uparrow}(r) \\ v_{n\downarrow}(r)\end{pmatrix}

where single particle hamiltonian is given by

h_{\sigma} = -\dfrac{1}{2}\vec{\nabla}\alpha_{\sigma}(r)\vec{\nabla} + V_{\sigma}(r)-\left(\mu_{\sigma}-V_{\sigma}^{\textrm{(ext)}}(r)\right)-\dfrac{i}{2}\left\lbrace \vec{A}_{\sigma}(r)-\vec{v}_{\sigma}^{\textrm{(ext)}}(r),\vec{\nabla} \right\rbrace

Potentials entering the hamiltonian are:

• alpha_a and alpha_b:
  \alpha_{\sigma} = 2\dfrac{\delta\mathcal{E}_{\textrm{edf}}}{\delta \tau_{\sigma}} -- effective mass,
• V_a and V_b:
  V_{\sigma}=\dfrac{\delta\mathcal{E}_{\textrm{edf}}}{\delta n_{\sigma}} -- mean-field potential,
• A_a_x, A_a_y, A_a_z, A_b_x, A_b_z, and A_b_z:
  \vec{A}_{\sigma}=\dfrac{\delta\mathcal{E}_{\textrm{edf}}}{\delta \vec{j}_{\sigma}} -- current potential,
• delta:
  \Delta(r)=-\dfrac{\delta\mathcal{E}_{\textrm{edf}}}{\delta \nu^*} -- paring potential.

Potentials are accessible for user through structure wslda_potential:

typedef struct
{
    int nx;             
    int ny;             /// for 1d code it is set to 1
    int nz;             /// for 1d and 2d code it is set to 1
    int datadim;        /// dimensionality of data
    int blocklength;    /// number of elements in potential array = nx*ny*nz
    
    // potentials
    double complex *delta;
    double *alpha_a; /// effective mass, spin-a
    double *alpha_b; /// effective mass, spin-b
    double *V_a;
    double *V_b;
    double *A_a_x;
    double *A_a_y;
    double *A_a_z;
    double *A_b_x;
    double *A_b_y;
    double *A_b_z;
} wslda_potential;