from scipy import constants as c from matplotlib import pyplot as plt from scipy.optimize import curve_fit from scipy.integrate import simps from math import pi import numpy as np #============================================ POTENTIALS ================================================ def pot_Volkov(r) : return(144.86*np.exp(-r*r/0.82**2) - 83.34*np.exp(-r*r/1.6**2)) def pot_MTV(r) : return(1458.05*np.exp(-r*3.11)/r - 578.09*np.exp(-1.55*r)/r) def pot_ATS3(r) : return( 1000.0*np.exp(-r*r*3.0) - 326.7*np.exp(-r*r*1.05) - 43.0*np.exp(-r*r*0.6) ) def pot_Minnesota(r) : return( 200.*np.exp(-1.487*r**2) - 178.*np.exp(-0.639*r**2) ) def pot_res(r) : return( -8.*np.exp(-0.16*r**2) + 4.*np.exp(-0.04*r**2) ) #======================================================================================================== def normalization(u) : norm = simps(np.abs(u)**2, r_grid) u_norm = u/np.sqrt(norm) return(u_norm) #======================================================================================================== def numerov_wf(E, potential) : k = np.zeros(len(r_grid)) u = np.zeros(len(r_grid)) C = 2.*m_red/hbarc**2 #C = 1./41.47 #initial conditions [for r->0; u(r)=r^(l+1)] chnage -> start from 0 u[0] = r_grid[0]**(l+1) u[1] = r_grid[1]**(l+1) #initializing k[r] = 2\mu/hbarc^2 (E - V[r] - hbarc^2 * l*(l+1)/(2.*\mu * r^2)) k[0] = 0 for i in range(1,len(r_grid)): k[i] = C * (E - potential(r_grid[i]) - l*(l+1) / (r_grid[i]**2 * C)) # recursive relation for forward Numerov algorithm for i in range(2,len(r_grid)) : u[i] = ( 2. * (1.-5.*h**2/12.*k[i-1]) * u[i-1] - (1.+h**2/12.*k[i-2]) * u[i-2] ) / (1.+h**2/12.*k[i]) u_norm = normalization(u) return(u_norm) #======================================================================================================== def bound_states(potential) : print("Bound state calculation : \n") energy_grid_points = 10 #numerical accuracy of bound state energy solution eps = 1.e-8 #maximal number of iterations itr_max = 100 itr=0 Eb_prev = Eb_min Eb_curr = Eb_max while itreps : itr+=1 Eb_grid = np.linspace(Eb_prev, Eb_curr, num=energy_grid_points, endpoint=True) for i in range(0,len(Eb_grid)) : u = numerov_wf(Eb_grid[i], potential) if i>0 and u_prev*u[-1] < 0: Eb_curr=Eb_grid[i] line = str(itr) + ". \t E \t = \t"+ "{0:.6f}".format(Eb_curr) + "\t \t u(rmax) \t = \t " + "{0:.6f}".format(u[-1]) print(line) break else : u_prev = u[len(u)-1] Eb_prev = Eb_grid[i] print("=========================================================================") return(Eb_curr,u) #======================================================================================================== def bess_j(x) : if l==0 : return(np.sin(x)/x) elif l==1 : return(np.sin(x)/x**2 - np.cos(x)/x) elif l==2 : return(3.*np.sin(x)/x**3 - 3.*np.cos(x)/x**2 - np.sin(x)/x) else : print("Orbital momentum l<3 --- > STOP") exit(0) #======================================================================================================== def bess_n(x) : if l==0 : return(-np.cos(x)/x) elif l==1 : return(- np.cos(x)/x**2 - np.sin(x)/x) elif l==2 : return(-3.*np.cos(x)/x**3 - 3.*np.sin(x)/x**2 + np.cos(x)/x) else : print("Orbital momentum l<3 --- > STOP") exit(0) #======================================================================================================== def calc_phaseshifts(u,k) : #matching is performed at two last points (rmax-h) and rmax r1 = r_grid[len(r_grid)-2] r2 = r_grid[len(r_grid)-1] u1 = u[len(r_grid)-2] u2 = u[len(r_grid)-1] f = r1 * u2 / (r2 * u1) #beta tan_delta = (f * bess_j(k*r1) - bess_j(k*r2) ) / (f * bess_n(k*r1) - bess_n(k*r2) ) return(k/tan_delta) #======================================================================================================== def scattering(potential) : print("Phaseshifts calculation : ") phaseshifts_grid = np.zeros(len(kscatt_points)) for i in range(0,len(kscatt_points)) : E = kscatt_points[i]**2/(2.*m_red)*hbarc**2 u = numerov_wf(E, potential) phaseshifts_grid[i] = calc_phaseshifts(u, kscatt_points[i]) print("Done") print("=========================================================================") return(phaseshifts_grid) #======================================================================================================== def rms_radius(u) : rms = simps(np.abs(u)**2 * r_grid**2 , r_grid) return(np.sqrt(rms)) #======================================================================================================== def ere_funct(x, a, reff, p) : return(-1/a + 0.5*reff*x + p*reff**3*x**2) #======================================================================================================== def ere_funct1(x, a, reff) : return(-1/a + 0.5*reff*x) #======================================================================================================== def ere_param(phsf) : popt, pcov = curve_fit(ere_funct, kscatt_points**2, phsf) perr = np.sqrt(np.diag(pcov)) a = popt[0] a_err = perr[0] reff = popt[1] reff_err = perr[1] p = popt[2] p_err = perr[2] return(a, a_err, reff, reff_err, p, p_err) #======================================================================================================== def bound_energy_scatt(a,reff) : return(( ( 2.*reff/a -2. +2.*np.sqrt(1-2.*reff/a) ) / reff**2 ) / (2.*m_red) * hbarc**2) #======================================================================================================== def output(Eb, u, phaseshifts_grid) : #=======================plotting bound state wave function======================= plt.subplot(1,2,1) plt.subplots_adjust(wspace=0.5) plt.plot(r_grid,u,"-",color="royalblue") #zero line zline = np.zeros((2,2)) zline[0,0] = -100. zline[0,1] = 100. zline[1,0] = 0 zline[1,1] = 0 plt.plot(zline[0,:],zline[1,:],"--",color="k") plt.xlim(0,rmax) plt.ylim(0,1.0) plt.xlabel("r [fm]") plt.ylabel("u(r)") #saving bound state wave function to file output = open("u_wf.dat","w") output.write("#r [fm] \t\t u\n") for i in range(0,rgrid_points) : line="{0:.6f}".format(r_grid[i])+"\t"+"{0:.6f}".format(u[i])+"\n" output.write(line) output.close() #=============================plotting phaseshifts============================== plt.subplot(1,2,2) plt.plot(kscatt_points**2, phaseshifts_grid,"-",color="royalblue") plt.xlabel(r'$k^2$ [${\rm fm}^{-2}$]') plt.ylabel(r'$k {\rm cot}(\delta)$ [${\rm fm}^{-1}$]') plt.xlim(0, kscatt_max**2) #plt.subplot(1,3,3) #plt.plot(kscatt_points[:]**2/2.,np.arctan(kscatt_points[:]/phaseshifts_grid[:])/pi,"-",color="royalblue") #plt.plot(kscatt_points[:]**2/2.,kscatt_points[:]/phaseshifts_grid[:],"-",color="royalblue") #plt.xlabel(r'$k^2$ [${\rm fm}^{-2}$]') #plt.ylabel(r'$k {\rm cot}(\delta)$ [${\rm fm}^{-1}$]') #plt.xlim(0, kscatt_max**2/2.) #plt.ylim(-5,5) #saving phase shifts to file output = open("phaseshifts.dat","w") output.write("#k [fm^-1] \t\t kcotd [fm^-1]\n") for i in range(0,nkscatt_points) : line="{0:.6f}".format(kscatt_points[i])+"\t"+"{0:.6f}".format(phaseshifts_grid[i])+"\n" output.write(line) output.close() #================================other properties=============================== line="Bound state properties :\n\n" #Bound state energy line+=r'E_b = '+"{0:.6f}".format(Eb)+" MeV"+"\n" #Rms radius #rms = rms_radius(u) #line+=r'\sqrt{} = '+"{0:.6f}".format(rms)+" fm" + "\n" #scattering length and effective range line+="\n\nScattering properties :\n\n" a, a_err, reff, reff_err, p, p_err = ere_param(phaseshifts_grid) line+=r'a = '+ "{0:.4f}".format(a)+r' +- '+"{0:.4f}".format(a_err)+" fm"+"\n" line+=r'r_eff = '+ "{0:.4f}".format(reff)+r' +- '+"{0:.4f}".format(reff_err)+" fm"+"\n" line+=r'P = '+ "{0:.4f}".format(p)+r' +- '+"{0:.4f}".format(p_err)+" fm"+"\n" #Bound state energy from a and reff Eb_scatt = bound_energy_scatt(a,reff) line+="\n\n"+r'Bound state energy check from (a;r_eff) :'+"\n\n" line+=r'E_b^{scatt} = '+"{0:.6f}".format(Eb_scatt)+" MeV"+"\n" print(line) print("=========================================================================") plt.show() #======================================================================================================== ##################################################### # # # NUMEROV ALGORITHM FOR BOUND STATES AND SCATTERING # # # ##################################################### #constants hbarc = 197.3269804 #particle masses [MeV] m1 = c.value("proton mass energy equivalent in MeV") m2 = c.value("neutron mass energy equivalent in MeV") #m1 = hbarc**2*2. #m2 = hbarc**2*2. m_red = m1*m2/(m1+m2) #orbital mometum (here we assume l <= 3) l = 0 #step [fm] h = 0.1 #relative distance (boundaries) [fm] and r_grid rmin = 0 rmax = 60 rgrid_points = int((rmax-rmin)/h) + 1 r_grid = np.linspace(rmin, rmax, num=rgrid_points, endpoint=True) # energy region Eb in (bound states only) [MeV] Eb_min = -5. Eb_max = 0 # relative momentum array for scattering phaseshifts (scattering states only) [fm^-1] kscatt_min = 0.00001 kscatt_max = 0.1 nkscatt_points = 1000 kscatt_points = np.linspace(kscatt_min, kscatt_max, num = nkscatt_points, endpoint=True) #################BOUND STATE CALCULATION#################### Eb, u = bound_states(pot_Volkov) ##################SCATTERING CALCULATION#################### phaseshifts_grid = scattering(pot_Volkov) ###########################OUTPUT########################### output(Eb, u, phaseshifts_grid)