Functionals.md

 [[_TOC_]]
 
 # Selecting the functional
-The functional must be selected at compilation stage in *predefines.h* file:  
+The functional must be selected at the compilation stage in *predefines.h* file:  
 ```c
 /**
  * Select functional:
  *  - SLDA: 
- *      for simulating unitary Fermi gas, 
- *      it sets effective mass of particles to 1.0 which assures better convergence properties,
- *      in case of time time-dependent calculations SLDA is about 2x faster than ASLDA.
+ *      For simulating a unitary Fermi gas. 
+ *      It is equivalent to selecting ASLDA functional with SLDA_FORCE_A1 option.
  *  - ASLDA:
- *      for simulating unitary Fermi gas,
- *      at qualitative level it produces results compatible with SLDA, however it is more accurate,
- *      due to presence of current terms in the functional it has worse convergence properties.
+ *      For simulating a unitary Fermi gas.
+ *      At the qualitative level, it produces results compatible with SLDA; however, it is more accurate.
+ *      Due to the presence of current terms in the functional, it has worse convergence properties.
+ *      For more info, see: https://arxiv.org/abs/1008.3933
+ *  - SLDAE:
+ *      For simulating a Fermi gas for an arbitrary value of akF.
+ *      For small and negative akF, the functional is compatible with BDG, while for large akF, it is compatible with ASLDA.
+ *      For more info, see: https://arxiv.org/abs/2201.07626
  *  - BDG:
- *      for simulating systems in BCS regime,
+ *      for simulating systems in the BCS regime,
  *      equations of motion are equivalent to Bogoliubov-de-Gennes equations,
- *      you MUST set aBdG value in input file when using this functional.
+ *  - CUSTOMEDF:
+ *      Use this option to define your custom functional.
+ *      Then you need to provide the body of functions: compute_energy_custom( ) and compute_potentials_custom( )
+ *      in problem-definition.h file
  * */
 // #define FUNCTIONAL SLDA
 #define FUNCTIONAL ASLDA
@@ -37,7 +44,7 @@ The ASLDA functional has been designed in order to capture properties of **stron
     + [1-\alpha_{\downarrow}(n_{\uparrow},n_{\downarrow})]\dfrac{\bm{j}_{\downarrow}^2}{2n_{\downarrow}}
 \end{aligned}
 ```  
-The functional is fitted to quantum Monte Carlo data, and for spin symmetric and uniform systems it provides:
+The functional is fitted to quantum Monte Carlo data, and for spin symmetric and uniform systems, it provides:
 ```math
 E/E_{\textrm{ffg}}=\xi=0.40(1),\qquad\Delta/\varepsilon_F=0.504(24).
 ```  
@@ -49,29 +56,30 @@ In the case of calculations for trapped system term $`\frac{\bm{j}_{\sigma}^2}{2
 * $`\bm{j}_{\sigma}\rightarrow 0`$,
 * $`\frac{\bm{j}_{\sigma}^2}{2n_{\sigma}}\rightarrow 0`$.  
 
-However, the division of very small numbers is numerically not stable operation. For this reason, we introduce the stabilization procedure:  
+However, the division of very small numbers is numerically unstable. For this reason, we introduce the stabilization procedure:  
 ```math
 \dfrac{\bm{j}_{\sigma}^2}{2n_{\sigma}}\longrightarrow f_{\textrm{reg.}}(n_{\sigma})\dfrac{\bm{j}_{\sigma}^2}{2n_{\sigma}}
 ```
-Role of $`f_{\textrm{reg.}}`$ is to exclude from computation regions of small density. The function is controlled by two parameters in `predefines.h` file:
+The role of $`f_{\textrm{reg.}}`$ is to exclude from computation regions of small density. The function is controlled by two parameters in `predefines.h` file:
 ```c
 /**
- * Meaningful only in case of ASLDA.
- * Parameters defining stabilization procedure of ASLDA functional. 
- * For regions with density smaller than SLDA_STABILIZATION_EXCLUDE_BELOW_DENISTY 
- * contribution from current term j^2/2n is assumed to be zero. 
- * For regions with density above SLDA_STABILIZATION_RETAIN_ABOVE_DENSITY 
+ * Meaningful only in the case of ASLDA and SLDAE.
+ * Parameters defining the stabilization procedure of the ASLDA functional. 
+ * For regions with density smaller than SLDA_STABILIZATION_EXCLUDE_BELOW_DENISTY, 
+ * contribution from the current term j^2/2n is assumed to be zero. 
+ * For regions with density above SLDA_STABILIZATION_RETAIN_ABOVE_DENSITY, 
  * the contribution is assumed to be intact by the stabilization procedure. 
+ * For more info see: Wiki -> Stabilization of ASLDA functional
  * */
 #define SLDA_STABILIZATION_RETAIN_ABOVE_DENSITY  1.0e-5
 #define SLDA_STABILIZATION_EXCLUDE_BELOW_DENISTY 1.0e-7
 ```
 Their meaning is presented in the figure below. 
 ![stabilization](uploads/740d294dbfb30a203cabec0b140fa326/stabilization.png)   
-A smooth transition between densities `SLDA_STABILIZATION_EXCLUDE_BELOW_DENISTY` and `SLDA_STABILIZATION_RETAIN_ABOVE_DENSITY` is introduced to avoid discontinuities for quantities, that may lead to divergences when computing derivatives.
+A smooth transition between densities `SLDA_STABILIZATION_EXCLUDE_BELOW_DENISTY` and `SLDA_STABILIZATION_RETAIN_ABOVE_DENSITY` is introduced to avoid discontinuities for quantities that may lead to divergences when computing derivatives.
 
 ## Chemical potential
-The spin-symmetric and uniform unitary Fermi is expected to be scale-invariant where the following relations for the total energy and the chemical potential satisfy: $`E=\xi E_\textrm{ffg}`$ and $`\mu=\xi \varepsilon_{\textrm{F}}`$, where $`\xi\approx 0.4`$ is the Bertch parameter. However, in numerical realization, the relation for the chemical potential is satisfied only in the limit of very low densities. It is demonstrated in the table below, where results for the uniform solution obtained on lattice `128x128x128`  with `DX=1` is presented, and $`k_{\textrm{F}}=\sqrt{2\varepsilon_{\textrm{F}}}=(6\pi^2 n_{\uparrow})^{1/3}`$.  
+The spin-symmetric and uniform unitary Fermi is expected to be scale-invariant, where the following relations for the total energy and the chemical potential satisfy: $`E=\xi E_\textrm{ffg}`$ and $`\mu=\xi \varepsilon_{\textrm{F}}`$, where $`\xi\approx 0.4`$ is the Bertch parameter. However, in numerical realization, the relation for the chemical potential is satisfied only in the limit of very low densities. It is demonstrated in the table below, where results for the uniform solution obtained on lattice `128x128x128`  with `DX=1` is presented, and $`k_{\textrm{F}}=\sqrt{2\varepsilon_{\textrm{F}}}=(6\pi^2 n_{\uparrow})^{1/3}`$.  
 
 |$`k_{\textrm{F}}`$| $`E/E_\textrm{ffg}`$ | $`\mu/\varepsilon_{\textrm{F}}`$|
 | -----------------|----------------------|---------------------------------|
@@ -81,10 +89,10 @@ The spin-symmetric and uniform unitary Fermi is expected to be scale-invariant w
 | 0.4 | 0.397 | 0.409 | 
 | 0.2 | 0.397 | 0.403 | 
 
-It is seen that the total energy $`E/E_\textrm{ffg}=\xi`$ satisfies the relation all the time, while the chemical potential $`\mu/\varepsilon_{\textrm{F}}\rightarrow\xi`$ only in low-density limit. For raw data see [energy-vs-kF.txt](uploads/221811409645103fcb1a2b3a1e506b94/energy-vs-kF.txt)
+It is seen that the total energy $`E/E_\textrm{ffg}=\xi`$ satisfies the relation all the time, while the chemical potential $`\mu/\varepsilon_{\textrm{F}}\rightarrow\xi`$ only in the low-density limit. For raw data see [energy-vs-kF.txt](uploads/221811409645103fcb1a2b3a1e506b94/energy-vs-kF.txt)
  
 # SLDA - superfluid local density approximation
-The term in ASLDA functional that depends on the currents $`\bm{j}_{\sigma}`$ introduces a significant cost to the computation. This term is responsible for maintaining Gallilean invariance of the ASLDA theory. The effective mass was found to be consistent with the bare mass to within 10% for a large range of polarizations,  Therefore, in many applications, one can set $`\alpha_{\sigma}=1`$ without losing qualitative features of ASLDA theory. Under this assumption we obtain SLDA functional:
+The term in ASLDA functional that depends on the currents $`\bm{j}_{\sigma}`$ introduces a significant cost to the computation. This term is responsible for maintaining the Galilean invariance of the ASLDA theory. The effective mass was found to be consistent with the bare mass to within 10% for a large range of polarizations. Therefore, in many applications, one can set $`\alpha_{\sigma}=1`$ without losing qualitative features of ASLDA theory. Under this assumption, we obtain the SLDA functional:
 ```math
 \mathcal{E}_{\textrm{edf}}(n,\nu,\ldots)  =  \frac{\tau_{\uparrow}}{2} + 
     \frac{\tau_{\downarrow}}{2} 
@@ -94,58 +102,59 @@ The term in ASLDA functional that depends on the currents $`\bm{j}_{\sigma}`$ in
     g(n_{\uparrow},n_{\downarrow})\nu^{\dagger}\nu 
 ```  
 
-*Note 1*: SLDA functional exhibits much better convergence properties than ASLDA, i.e `st-wslda` codes typically converge in a significantly smaller number of iterations.  
+*Note 1*: SLDA functional exhibits much better convergence properties than ASLDA, i.e, `st-wslda` codes typically converge in a significantly smaller number of iterations.  
 
-*Note 2*: SLDA functional provides by a factor about 2x better performance in case of `td-wslda` codes than ASLDA variant.  
+*Note 2*: SLDA functional provides a factor of about 2x better performance in the case of `td-wslda` codes than the ASLDA variant.  
 
-Below we provide a plot showing relative energy change (in %) between simulations of SLDA and ASLDA that compute the energy of N fermions (both even and odd particle numbers) in the spin-symmetric unitary Fermi gas trapped in an isotropic harmonic oscillator.  For more info see [supplemental material of this paper](https://arxiv.org/abs/1306.4266).  
+Below, we provide a plot showing the relative energy change (in %) between simulations of SLDA and ASLDA that compute the energy of N fermions (both even and odd particle numbers) in the spin-symmetric unitary Fermi gas trapped in an isotropic harmonic oscillator.  For more info, see [supplemental material of this paper](https://arxiv.org/abs/1306.4266).  
 ![slda_vs_aslda](uploads/073efe244d1ab67bed1314af5cf0fc21/slda_vs_aslda.png)
 
 # BdG - Bogoliubov de-Gennes functional
-The BdG functional is equivalent to  Bogoliubov de-Gennes mean-field approximation. This approximation is valid in BCS regime, under assumption $`-ak_F<1`$. The functional has form:
+The BdG functional is equivalent to the Bogoliubov de-Gennes mean-field approximation. This approximation is valid in the BCS regime, under the assumption $`-ak_F<1`$. The functional has form:
 ```math
 \mathcal{E}_{\textrm{edf}}(n,\nu,\ldots)  =  \alpha_\uparrow\frac{\tau_{\uparrow}}{2} + 
     \alpha_\downarrow\frac{\tau_{\downarrow}}{2} 
     + g\nu^{\dagger}\nu 
 ```  
-The (bare) coupling constant $`g`$ is related to scattering length by the formula $`g=4\pi\hbar^2a/m_r`$, where $`m_r`$ is twice of reduced mass, i.e.: $`\frac{1}{m_r}=\frac{1}{2}(\alpha_\uparrow+\alpha_\downarrow)`$. Default values of effective masses are $`\alpha_\uparrow=\alpha_\downarrow=1`$. See [mass imbalanced gas in a harmonic trap](https://gitlab.fizyka.pw.edu.pl/gabrielw/wslda/-/wikis/st-wslda-examples#example-3-mass-imbalanced-gas-in-a-harmonic-trap) to learn how to modify effective masses. 
+The (bare) coupling constant $`g`$ is related to scattering length by the formula $`g=4\pi\hbar^2a/m_r`$, where $`m_r`$ is twice of reduced mass, i.e.: $`\frac{1}{m_r}=\frac{1}{2}(\alpha_\uparrow+\alpha_\downarrow)`$. Default values of effective masses are $`\alpha_\uparrow=\alpha_\downarrow=1`$. See [mass imbalanced gas in a harmonic trap](https://gitlab.fizyka.pw.edu.pl/wtools/wslda/-/wikis/st-wslda-examples-archival#example-3-mass-imbalanced-gas-in-a-harmonic-trap) to learn how to modify effective masses. 
 
-Below we provide predictions of BdG functional for total energy and paring gap as a function of coupling constant $`-ak_F`$ computed for uniform and spin-symmetric system. Results are compared with predictions of BCS theory, and as expected the agreement is observed for $`-ak_F<1`$.  
+Below we provide predictions of BdG functional for total energy and pairing gap as a function of coupling constant $`-ak_F`$ computed for a uniform and spin-symmetric system. Results are compared with predictions of BCS theory, and as expected, the agreement is observed for $`-ak_F<1`$.  
 ![BdG-E](uploads/a3e8e87c80efd2c4a31d2ea50dec6191/BdG-E.png)  
 ![BdG-D](uploads/7529bf3d47e71594ef42d2b2b3a51865/BdG-D.png) 
 
 # SLDAE -  SLDA Extended
-SLDAE is extension of SLDA-type functional to arbitrary value of $`\lambda=|ak_F|`$. The functional has generic form:
+SLDAE is an extension of the SLDA-type functional to an arbitrary value of $`\lambda=|ak_F|`$. The functional has a generic form:
 ```math
 	\mathcal{E}  = A_\lambda \frac{\tau}{2}
 	        + \frac{3}{5}B_\lambda n \varepsilon_F
 	        + \frac{C_\lambda }{n^{1/3}} |\nu|^2
                 + (1 - A_\lambda) \frac{j^2}{2n}
 ``` 
-where functional parameters $`A_\lambda`$, $`B_\lambda`$ and $`C_\lambda`$ are constructed in such a way to reproduce properly quasiparticle properties in the whole regime of $`\lambda`$, and also satisfy constraints in limits $`\lambda\rightarrow 0`$ and $`\lambda\rightarrow \infty`$. Precisely, the SLDAE is constructed in such a way to reproduce data for the:
+where functional parameters $`A_\lambda`$, $`B_\lambda`$, and $`C_\lambda`$ are constructed in such a way to reproduce quasiparticle properties in the whole regime of $`\lambda`$ properly, and also satisfy constraints in limits $`\lambda\rightarrow 0`$ and $`\lambda\rightarrow \infty`$. Precisely, the SLDAE is constructed in such a way to reproduce data for the:
 * the ground-state energy per unit volume $`E = 3n\xi_\lambda \varepsilon_F/5`$, and thus corresponding chemical potential $`\mu/\varepsilon_F = \zeta_\lambda`$ which must verify  the thermodynamic relationship $`\zeta_\lambda = \xi_\lambda + (\lambda/5) \xi_\lambda^\prime`$,
 * the pairing gap function $`\Delta_\lambda=\eta_\lambda \varepsilon_F`$,
 * the effective mass of particle $`\alpha_\lambda=m/m^\star`$.
 The plot below shows $`\lambda=|ak_F|`$ dependence of these quantities in the SLDAE functional. 
 ![sldae_qp_properties](uploads/bc7018319141055f9f032d7526a2dc9c/sldae_qp_properties.png)
-For more details see https://arxiv.org/abs/2201.07626.  
+For more details, see https://arxiv.org/abs/2201.07626.  
 **NOTE**: the functional has been constructed for spin-symmetric systems, $`N_a=N_b`$.
 
 ## SLDAE with $`m/m^\star=1`$
-The last term of SLDAE functional (depending on currents $`j`$) introduces a significant cost to the computation process. In many cases, it is sufficient to proceed with the functional that neglects corrections related to the effective mass (actually the effective mass is not known with high precision). The W-SLDA Toolkit provides a variant of the SLDAE functional where the constraint $`\alpha_\lambda = m/m^\star=1`$ is imposed. In order to activate it you need to use flag (in `predefines.h`):
+The last term of SLDAE functional (depending on currents $`j`$) introduces a significant cost to the computation process. In many cases, it is sufficient to proceed with the functional that neglects corrections related to the effective mass (actually the effective mass is not known with high precision). The W-SLDA Toolkit provides a variant of the SLDAE functional where the constraint $`\alpha_\lambda = m/m^\star=1`$ is imposed. In order to activate it, you need to use a flag (in `predefines.h`):
 ```c
 /**
- * Meaningful only in case SLDAE.
- * Sets effective mass to be equal, and speeds-up computation (approximately by a factor of two)
+ * Sets the effective mass to be equal to 1.
+ * Computation converges faster, but it also decreases the accuracy of the functional.
+ * If functional is SLDA or BDG, this option is activated automatically.
  * */
-#define SLDAE_FORCE_A1
+#define SLDA_FORCE_A1
 ```
 ## Custom SLDAE
-The W-SLDA Toolkit allows for creating custom SLDAE-like functionals. For more info see [Custom SLDAE functional](Custom SLDAE functional).
+The W-SLDA Toolkit allows the creation of custom SLDAE-like functionals. For more info, see [Custom SLDAE functional](Custom-SLDAE-functional).
 
 # Regularization scheme of pairing field
 By default, W-SLDA implements a regularization scheme of the pairing field known as Superfluid Local Density Approximation (SLDA). Precisely, the implemented formulas are (9.86)-(9.87) from paper [arXiv:1008.3933](https://arxiv.org/abs/1008.3933).  
 
-*Note 1*: The regularization scheme is only valid for 3D systems. If you want to execute pure 2D or 1D calculations, you need to provide your regularization procedure. For more info see [Strict 2D or 1D modes](Strict 2D or 1D mode).  
+*Note 1*: The regularization scheme is only valid for 3D systems. If you want to execute pure 2D or 1D calculations, you need to provide your regularization procedure. For more info, see [Strict 2D or 1D modes](Strict 2D or 1D mode).  
 
-*Note 2*: The regularization uses average effective mass $`\alpha_{+}`$ and average chemical potential $`\mu_{+}`$. Thus, the regularization is valid for spin-imbalanced and mass-imbalanced systems as well. 
+*Note 2*: The regularization uses average effective mass $`\alpha_{+}`$ and average chemical potential $`\mu_{+}`$. Thus, the regularization is valid for both spin- and mass-imbalanced systems.