CODE 1D DPFT
=======================

:Date: 8 Oct 2019

Note
------

The :code:`.rst` file for this page is generated using Jupyter Notebook, from the :code:`.ipynb` file: `1D DPFT In Python <https://github.com/tesla-cat/Physics/blob/master/ipynb/1D_DPFT_In_Python.ipynb>`__


Thomas Fermi Kinetic Energy
---------------------------------------------

.. math::

   E_{\text{kin}}=\int dx\; {1\over 3}\alpha \bar{n}^3 \\
   \Rightarrow {\delta E_{\text{kin}} \over \delta n}=\alpha n^2 

Where :math:`\alpha \equiv {\pi^2\hbar^2\over 2mg^2}`

And :math:`g` is spin multiplicity. For spin-polarized fermi gas (For
example electron gas with all the electrons spin-up), :math:`g=1` from
Pauli principle

Without Interaction
---------------------------------------------

.. math::

   \begin{cases}
    \delta\left(E_{\text{kin}}+E_{\text{ext}}\right)=0 \\
   {\delta E_{\text{ext}} \over \delta n}= V_{\text{ext}} - \mu
    \end{cases}
   \\ \Rightarrow n_{\text{nonInt}}[\mu-V_{\text{ext}}]=\sqrt{\alpha^{-1}[\mu-V_{\text{ext}}]_+}

is a **Functional** of :math:`\mu-V_{\text{ext}}`

Where :math:`\mu` is the **chemical potential**

The Interaction
---------------------------------------------

.. math::

   E_{\text{int}} = w{1\over 2}\int{dr\; n^2(r)} \\
   \Rightarrow V_{\text{int}}={\delta E_{\text{int}} \over \delta n}=wn

-  Where :math:`w` is strength of interaction

-  Then we have

.. math::


   \delta\left(E_{\text{kin}}+E_{\text{ext}}+E_{\text{int}} \right)=0
   \\ \Rightarrow n_{\text{int}}(x,\mu)=\left(-{w\over 2\alpha}+\sqrt{ ( {w\over 2\alpha})^2 +{\mu-V_{\text{ext}} \over \alpha}  }\right) \eta(\mu-V_{\text{ext}})

Self Consistent Loop
---------------------------------------------

**What is the idea?**


-  We will find :math:`n` **with interaction**
-  However, **we will NOT use the expression** :math:`n_{\text{int}}`
   above
-  Instead, **we use the expression** :math:`n_{\text{nonInt}}`:

   1. Specifically we will use :math:`n=\sqrt{\alpha^{-1}[\mu-V]_+}`
   2. We initialize :math:`V` with :math:`V_{\text{ext}}`
   3. Then we find the chemical potential :math:`\mu` using **gradient
      descent**, the loss function would be :math:`N-\int n(x) \;dx`, we
      do this to make sure the total particle number is constant
   4. Then we calculate :math:`n_{\text{modified}}` based on :math:`\mu`
      and :math:`V` using above formula
   5. Then we update :math:`n` by mixing the **old** :math:`n` with
      :math:`n_{\text{modified}}`, that is
      :math:`n=(1-p)n+p\; n_{\text{modified}}`, the idea is similar to
      **gradient descent**, the :math:`p` here can be understood as
      **learning rate**
   6. Using the **new** :math:`n` we can calculate
      :math:`V_{\text{int}}=wn`, then we **add the interaction energy
      to** :math:`V_{\text{ext}}` as our modified :math:`V`, that is
      :math:`V=V_{\text{ext}}+V_{\text{int}}`
   7. Then we can go back to :math:`3` until :math:`n` converges

Simple implementation in Python
---------------------------------------------

.. code:: python

    import numpy as np
    import matplotlib.pyplot as plt
    
    class SimpleDPFT:
      def __init__(self,args):
        self.Vext = args["Vext"]
        self.maxIteration = args["maxIteration"]
        self.mu = args["mu"]
        self.alpha = args["alpha"]
        self.N = args["N"]
        self.dx = args["dx"]
        self.learningRate = args["learningRate"]
        self.gradientMax = args["gradientMax"]
        self.w = args["w"]
        self.p = args["p"]
        self.nChangeMax = args["nChangeMax"]
      
      def selfConsistentLoop(self):
        self.V = self.Vext
        self.mu = self.enforceN(self.mu)
        self.n = self.getDensity(self.mu)
        for i in range(self.maxIteration):
          self.V = self.Vext + self.w * self.n
          self.mu = self.enforceN(self.mu)
          nOld = self.n
          self.n = (1-self.p)*self.n + self.p*self.getDensity(self.mu)
          if(np.sum(abs(self.n-nOld))<self.nChangeMax):
            break
        return self.n
    
      def enforceN(self, mu):
        # gradient descent
        for i in range(self.maxIteration):
          gradient = (self.cost(mu+1)-self.cost(mu))
          mu = mu - self.learningRate * gradient
          if(abs(gradient)<self.gradientMax):
            break
        return mu
      def cost(self, mu):
        return abs(self.N - np.sum(self.getDensity(mu)*self.dx) )
    
      def getDensity(self,mu):
        muMinusV = mu - self.V
        muMinusV[ muMinusV<0 ] = 0
        return np.sqrt(1/self.alpha * muMinusV)
    
    # =================================================
    # ================= Usage =========================
    # =================================================
    x = np.linspace(-10, 10, num=100)
    dx = x[1]-x[0]
    args = {
        "Vext": x**2,
        "maxIteration": 4000,
        "mu": 1,
        "alpha": np.pi**2/2,
        "N": 20,
        "dx": dx,
        "learningRate": 1e-1,
        "gradientMax": 1e-6,
        "w": 5,
        "p": 0.5,
        "nChangeMax": 1e-6,
    }
    
    for w in range(3):
      args["w"] = w*5
      simpleDPFT = SimpleDPFT(args)
      n = simpleDPFT.selfConsistentLoop()
      print("Particle Number Check:", np.sum(n)*dx,"/ 20")
      plt.plot(x,n,label="w="+str(w*5))
    plt.legend()
    plt.show()



.. parsed-literal::

    Particle Number Check: 19.651993996020014 / 20
    Particle Number Check: 19.524746349214976 / 20
    Particle Number Check: 19.440909250694254 / 20



.. image:: imgs/1D-DPFT-In-Python/output_5_1.webp