CODE 1D DPFT Momentum Space ================================== Thomas Fermi Approximation -------------------------- .. math:: E_1^{TF}[T-\mu]=g\int {drdp\over (2\pi)^3}[ T(p)+V_{\text{ext}}(r)-\mu ]\eta(\mu -T-V_{\text{ext}}) \\ \Rightarrow \rho={\delta E_1^{TF} \over \delta T}= g\int {dr \over (2\pi)^3} 1\cdot \eta(\mu -T-V_{\text{ext}})\\ = g{1 \over (2\pi)^3} {4\over 3}\pi R^3={g \over 6\pi^3}R^3 \\ \text{Where }R \text{ is the radius of the sphere such that } \mu -T-V_{\text{ext}}>0 - For harmonic oscillator in 1D, and ignore the various constants, we have .. math:: V=x^2\Rightarrow \mu -T- x^2>0 \Rightarrow x < \sqrt{\mu-T}\equiv R\\ \Rightarrow \rho= g{1 \over 2\pi} 2R ={g \over \pi} \sqrt{\mu-T} Effective Kinetic Envergy :math:`T` ----------------------------------- .. math:: T(p)={\delta \over \delta \rho}E_{\text{int}} +T_{\text{kin}}\\ \text{where } T_{\text{kin}} = {1\over 2}p^2 The Interaction --------------- - I don’t know what :math:`E_{\text{int}}` should look like at the moment, let just assume the same form as in [1D DPFT In Python] so that .. math:: T_{\text{int}} = { \delta E_{\text{int}} \over \delta \rho}= w \rho Why momentum space ? -------------------- - As can be seen above, in position space, we have fixed expression for density :math:`n=\sqrt{\alpha^{-1}[\mu - V_{\text{ext}}]}`, this expression does not change, only :math:`V_{\text{ext}}` changes - In momentum space - We start with :math:`\mu -T- V_{\text{ext}}(r)>0` , then :math:`R=V_{\text{ext}}^{\text{Inv}}(\mu-T)` and :math:`\rho = {g\over 6\pi^3}R^3` - The position space and momentum space are only equivalent, at least computationally, for harmonic oscillators! And this is a result of the fundamental relation :math:`T={1\over 2}p^2` - Therefore, depends on :math:`V_{\text{ext}}`, some problem may be easier to solve in momentum space! Simple DPFT in Momentum space (Hamonic Oscillator 1D) ----------------------------------------------------- .. code:: python import numpy as np import matplotlib.pyplot as plt class SimpleDPFT: def __init__(self,args): self.Tkin = args["Tkin"] self.maxIteration = args["maxIteration"] self.mu = args["mu"] self.g = args["g"] self.N = args["N"] self.dp = args["dp"] self.learningRate = args["learningRate"] self.gradientMax = args["gradientMax"] self.w = args["w"] self.p = args["p"] self.nChangeMax = args["nChangeMax"] def selfConsistentLoop(self): self.T = self.Tkin self.mu = self.enforceN(self.mu) self.n = self.getDensity(self.mu) for i in range(self.maxIteration): self.T = self.Tkin + 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.N): muNmin = muNmin/2 while(self.getN(muNmax)self.N): muNmax = middlePoint else: muNmin = middlePoint if(self.getN(muNmin)/self.getN(muNmax)>(1-1e-3)): break return muNmin def getN(self,mu): return np.sum(self.getDensity(mu)*self.dp) # ================================================= # ================= Usage ========================= # ================================================= p = np.linspace(-10, 10, num=100) dp = p[1]-p[0] args = { "Tkin": p**2/2, "maxIteration": 4000, "mu": 1, "g": 1, "N": 20, "dp": dp, "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)*dp,"/ 20") plt.plot(p,n,label="w="+str(w*5)) plt.legend() plt.show() .. parsed-literal:: Particle Number Check: 19.65459900651429 / 20 Particle Number Check: 19.560823172965314 / 20 Particle Number Check: 19.494538379097083 / 20 .. image:: imgs/1D-DPFT-In-Momentum-Space/output_6_1.webp Simple DPFT in Momentum space (Coulomb potential 1D) ------------------------------------------------------ :math:`{Z\over T-\mu}>r, \text{ Then }\rho = {g\over \pi}{Z\over T-\mu}` **Gradient descent doesn’t work for Coulomb potential, so we use Bisection** .. code:: python import numpy as np import matplotlib.pyplot as plt class SimpleDPFT: def __init__(self,args): self.p = args["p"] self.g = args["g"] self.N = args["N"] self.w = args["w"] self.mixRatio = args["mixRatio"] self.maxIteration=args["maxIteration"] self.Tkin = self.p**2/2 self.dp = self.p[1] - self.p[0] def selfConsistentLoop(self): self.T = self.Tkin self.mu = self.enforceN() self.n = self.getDensity(self.mu) for i in range(self.maxIteration): self.T = self.T + self.w * self.n self.mu = self.enforceN() nOld = self.n self.n = (1-self.mixRatio)*self.n \ + self.mixRatio *self.getDensity(self.mu) ratio = self.n/nOld if(np.sum(ratio)/ratio.shape[0] < (1+1e-3) \ and np.sum(ratio)/ratio.shape[0] > (1-1e-3)): break return self.n def enforceN(self): # mu must be negative # if abs(mu)*2, TminusMu increase, N decrease muNmin = muNmax = -1 while(self.getN(muNmin) > self.N): muNmin = muNmin*2 while(self.getN(muNmax) < self.N): muNmax = muNmax/2 if(muNmax==0): print("Max N you can get with such w:", self.getN(muNmax)) assert False, "Your N is not achievable since w is too large!" for i in range(self.maxIteration): middlePoint = (muNmin+muNmax)/2 if(self.getN(middlePoint)>self.N): muNmax = middlePoint else: muNmin = middlePoint if(self.getN(muNmin)/self.getN(muNmax)>(1-1e-3)): break return muNmin def getN(self,mu): return np.sum(self.getDensity(mu)*self.dp) def getDensity(self,mu): return self.g/np.pi/(self.T-mu) # ================================================= # ================= Usage ========================= # ================================================= p = np.linspace(-1,1,num=100) N = 0.1 args = { "p":p, "g":1, "N":N, "mixRatio":0.5, "maxIteration":1000, } wList = [0,20,40,80] for w in wList: args["w"] = w simpleDPFT = SimpleDPFT(args) n = simpleDPFT.selfConsistentLoop() print("Particle Number Check:",np.sum(n)*(p[1]-p[0]),"/",N) plt.plot(p,n,label="w="+str(w)) plt.legend() plt.show() .. parsed-literal:: Particle Number Check: 0.09995211570740857 / 0.1 Particle Number Check: 0.0999575484394774 / 0.1 Particle Number Check: 0.09996421750569122 / 0.1 Particle Number Check: 0.09994714485865916 / 0.1 .. image:: imgs/1D-DPFT-In-Momentum-Space/output_8_1.webp