DFTpy

This problem set gives a (broken) solver for the 1D Hohnberg-Kohn equations using the Thomas-Fermi approximation. The numerical methods are all there - fix the science and turn it in.

<source lang="python">

1. !/usr/bin/env python
3. QMII, USF Spring 2014
4. Homework 3 - due Thurs., Apr 10, 2014:
6. Solve the 1D Hohnberg-Kohn equations using the Thomas-Fermi approximation
7. to the energy functional, E = int v_ext n dx + T[n] + U[n]
9. Find solutions over L = 1nm in an external potential
10. given by a uniform electric field (in the +x direction) of 0.1 V / nm.
12. - Look at all 4 solutions (lm = -0.05, -0.1, -0.15, and -0.2 eV/nm^2)
13. - Your final result should use M=100 points.
15. 1. Fix the units in the code to make the length and energy scales
16. consistently in nanometers and eV / nm^2.
18. 2. Find the energy and multiplier, lm for the conditions above
19. (hint: loop over a range of lm from -1 to -10)
20. How are the two mathematically related?
22. 3. Describe the electron density that results.
24. 4. What does the solution suggest about modeling bonding with TF theory?
26. 5. What energy term is problematic when there is only 1 electron (and why)?

from numpy import * from scipy.optimize import fmin, fsolve

M = 100 eps0 = 8.854e-12 # <- check units!

1. length scale

dx = 1.0/M # nm

x = arange(M)*dx G = -abs(x[:,newaxis] - x[newaxis,:]) / eps0

Me = 0.91e-30 tf = 3*Planck**2/(10*Me)\

       *(3/(8*pi))**(2/3.0) \
       * 1.0 # TODO: more units...

v_ext = arange(M)*0.01 # <- check units!

1. v_ext = G[:,M/2]

def T(n):

   return sum(tf*n**(2/3.0) * n)*dx

def U(n):

   return dot(n, dot(G, n))*0.5*dx**2

def E(n):

   n = abs(n)
   return T(n) + U(n) + dx*sum(v_ext*n) - dx*lm*sum(n)

1. check units on these as well. Are the dx-es right?

def dT(n):

   return tf*(5/3.0)*n**(2/3.0)

def dU(n):

   return dot(G, n)*dx

def dE(n):

   n = abs(n)
   return dT(n) + dU(n) + v_ext - lm

def dE_2(m):

   n = cumsum(m)
   n += 10.0/M - sum(n)
   return dT(n) + dU(n) + v_ext

for lm in -0.05*arange(4):

   n0 = zeros(M)+1.0
   n = fsolve(dE, n0)
   #n = fmin(E, n0)
   n = abs(n)
   print n
   print sum(n)
   print E(n)

1. print G
2. print v_ext

</source>