h2.py (1786B)
1 #!/usr/bin/python3.9 2 3 4 from dolfin import * 5 6 import numpy as np 7 import os 8 9 10 #verwenden atomare einheiten. hier werden alle in der schrodingergleicchung notigen naturkonstannten 1 und die langeneinheit wird der bohrradius 11 12 def system(i): 13 14 d = i/100.0 #d in mesh units 15 dh = d/2.0 16 17 #define mesh and function space 18 mesh = UnitCubeMesh(20, 20, 20) #1 mesh unit = bohrradius = 0.529 Angstrom 19 V = FunctionSpace(mesh, 'CG', 1) 20 21 # Potential for the Hydrogen Atom 22 # { 1/r if r > 0 23 # V = { 24 # { 0 else 25 26 formula = 'sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])>0 ?\ 27 -(1/sqrt((x[0]-' + str(dh) + ')*(x[0]-' + str(dh) + ')+x[1]*x[1]+x[2]*x[2]) + 1/sqrt((x[0]+' + str(dh) + ')*(x[0]+' + str(dh) + ')+x[1]*x[1]+x[2]*x[2]) - 1.0/' + str(d) + ') : 0' 28 29 30 31 pot = Expression(formula,degree=3) 32 33 34 # Boundary 0 everywhere 35 def boundary(x, on_boundary): 36 return on_boundary 37 bc = DirichletBC(V, 0, boundary) 38 39 u = TrialFunction(V) 40 v = TestFunction(V) 41 42 # Assemble system 43 a = (1/2*inner(grad(u), grad(v)) + pot*u*v)*dx 44 45 # Define Matrices 46 A = PETScMatrix() 47 assemble(a, tensor=A) 48 49 # Apply boundary on the system 50 bc.apply(A) 51 52 # Create eigensolver 53 eigensolver = SLEPcEigenSolver(A) 54 eigensolver.parameters["spectrum"] = "smallest magnitude" 55 56 values = 1 57 eigensolver.solve(values) 58 59 u = Function(V) 60 61 62 for j in range(values): 63 E, E_c, R, R_c = eigensolver.get_eigenpair(j) 64 65 print('E_0 = ', E) 66 67 file = open("energy.dat","a") 68 file.write(str(d)) 69 file.write(" ") 70 file.write(str(E)) 71 file.write("\n") 72 file.close 73 74 def main(): 75 for i in range(100, 300, 2): 76 system(i) 77 78 if __name__ == "__main__": 79 main()