from scipy import * from scipy import weave from scipy import integrate from scipy import optimize from scipy import interpolate from pylab import * from scipy.misc import derivative code_Numerov=""" // void Numerov(const container& F, int Nmax, double dh, container& Solution) double dx = dh; double h2 = dx*dx; double h12 = h2/12; double w0 = (1-h12*F(0))*Solution(0); double Fx = F(1); double w1 = (1-h12*Fx)*Solution(1); double Phi = Solution(1); double w2; for (int i=2; i=(n_max-l): break u0=u1 dE = dE/1.091 Eb += Ebl E0 = Ebl[0][1] Eb.sort(comp) return Eb def ChargeDensity(Eb,R,Z,Z_screened): rho = zeros(len(R),dtype=float) Nelec=0 for i,(l,Ene) in enumerate(Eb): u = ComputeSchrod(Ene,R,l,Z_screened,True) dN = 2*(2*l+1) if Z >= Nelec+dN: ferm=1. else: ferm= (Z-Nelec)/(dN+0.0) Nelec += dN*ferm rho += u**2*(dN*ferm) / (4*pi*R**2) if Nelec>=Z: break print 'Adding ', dN*ferm, 'electrons of l=', l, 'Ene=', Ene return rho def rhsPotential(u, r, rhoSpline): """Right hand side (derivative) for the Poisson equation to computer V_H gives an array d/dr(U,U') == (U',U'') """ return (u[1], -8.*pi*r*rhoSpline(r)) def HartreeU(R,rho): dh=(R[-1]-R[0])/(len(R)-1.) Nmax = len(R) U = zeros(len(R),dtype=float) weave.inline(code_Hart, ['U', 'Nmax', 'R', 'rho'],type_converters=weave.converters.blitz, compiler = 'gcc') Solution = zeros(len(R),dtype=float) Solution[1]=5*dh weave.inline(code_NumerovU, ['U', 'Nmax', 'dh', 'Solution'],type_converters=weave.converters.blitz, compiler = 'gcc') return Solution Z_screened=1. Z=2 n_max=2 #5 lmax=1 #3 #R0 = linspace(1e-7,90,2**12+1) # For high precision R0 = linspace(1e-7,20,2**10+1) R=R0[::-1] Eb = FindBoundStates(Z_screened,R,n_max,lmax) rho = ChargeDensity(Eb,R,Z,Z_screened) # Computing U by odeint rhoSpline = interpolate.UnivariateSpline(R[::-1],rho[::-1],s=0) U = integrate.odeint(rhsPotential, [0.0, 5.], R0, args=(rhoSpline,))[:,0] # Computing U using Numerov U2 = HartreeU(R0,rho[::-1]) plot(R0, U) plot(R0, U2) show()