Open zzpwahaha opened 4 years ago
Hi @zzpwahaha
Thank you for your question.
& 3. It was made for integration of bound states. Could it be adopted in similar form for free states? Probably. However, note that for bound states integration starts not from 0 (as in your case) but from large radius R inwards. This is because for bound states large R value of wave function is well defined (we know wave function is 0 at large R), whereas precise potential for small R is not actually that well known. That's why radialWavefunction would normally stop with integration before reaching 0 [at self.alphaC**(1 / 3.0)
] for alkali atoms. In any case I guess that this difference in initial condition and direction of integration (you start with 0 wave function, and integrate from core outwards in your example) is making difference between two codes.
This is to adjust mesh so that we have denser mesh where radial wave function changes more quickly (close to core, for R->0) and vice versa. One can use also other meshes, for example have a constant step, but then number of integration steps required would be larger and it would take bit more time to get results.
Hope this clarifies all. I would have to think how to extend the existing method correctly for free states.
Hi Nikola Thanks for your reply. I have played around with the free state wf. Here is the way I did it
%%cython -a
import numpy as np
cdef extern from "math.h":
double fabs(double x)
def Numerovc(double [:] f, double x0, double dx, double dh):
cdef int divergentPoint = 0
cdef double maxValue = 1e16
cdef double [:] x = np.zeros(len(f))
x[0]=x0
x[1]=x0 + dh*dx
cdef double h2 = dh*dh
cdef double h12 = h2/12.
cdef double w0 = x[0]*(1-h12*f[0])
cdef double w1 = x[1]*(1-h12*f[1])
cdef double w2
cdef double xi = x[1]
cdef double fi = f[1]
for i in range(2,len(f)):
w2 = 2*w1 - w0 + h2*fi*xi # here fi=f1
fi = f[i] # fi=f2
xi = w2/(1-h12*fi)
x[i] = xi
w0 = w1
w1 = w2
# check convergence to avoid instability near origin
return x
def getRadialWf_Contin(Enconti,Rlsqrt,l,j):
feffsqrt = - np.array([effRadialPotentialsqrt(l,0.5,j, i, Enconti) for i in Rlsqrt])
ursqrt = Numerovc(feffsqrt,
0.0,
-1e-15,
Rlsqrt[1]-Rlsqrt[0]) * np.sqrt(Rlsqrt)
normsqrt = integrate.simps(ursqrt**2,x=Rlsqrt**2)
return ursqrt*1/np.sqrt(abs(normsqrt))
def effRadialPotentialsqrt(l, s, j, x, stateEnergy):
r = x * x
mu = (atom.mass - cc.m_e) / atom.mass
return -3. / (4. * r) + 4. * r * (
2. * mu * (stateEnergy - atom.potential(l, s, j, r))
- l * (l + 1) / (r**2)
)
where the method is the same as in ARC in the sense that we use the same sqrt(r) discretization as well as the parametric model potential. The start point of integration is at origin. The problem with that is in some case, the wavefunction's amplitude is too large and the normalization returns a Nan.
Hope that will help you Best, Zhenpu
Hi I am trying to compute the continuum radial wavefuntion of Rb87 with the parameteric model potential https://journals.aps.org/pra/pdf/10.1103/PhysRevA.49.982. I looked up the cpp code in
arc_c_extensions.c
and have couple of questions:x = sqrt(r)
. Is that for a better percision? since it is not the regular way we write the radial wavefunction.I try a Numerov method with the radial wf
u =r * R
, like this with the corresponding paramatric potential taken from you code. Below is my code. What I found is that the solution of the two methods sometimes agree but mostly do not. For example,En = 1.5, Rl = np.linspace(1e-7,100,10000)
does not, butEn = 0.5, Rl = np.linspace(1e-7,100,10000)
agree. I am wondering what could cause that, would it be the way of re-parametrization?