# DFTpy

From Predictive Chemistry

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">

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

- 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!

- 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)

- 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)

- print G
- print v_ext

</source>