1!*****************************************************************
2! Module for solving scalar relativistic radial equations
3!   Uses program adapted by Marc Torrent and Francois Jollet from
4!      USPS pgm of David Vanderbilt based on two coupled first order
5!      differential equations
6!      Previous version, based on second order differential equation
7!        from formalism of Shadwick, Talman, and Norman, Comp. Phys. Comm.
8!      54, 95-102 (1989)  found to be unstable
9!   09-16-06  NAWH
10!*****************************************************************
11
12!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13! This module contains the following active subroutines:
14!     Allocate_scalar_relativistic, deallocate_scalar_relativistic
15!        Azeroexpand, wfnsrinit, wfnsrasym, unboundsr, boundsr,
16!        scalarrelativisticturningpt, prepareforcfdsol
17!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
18
19#if defined HAVE_CONFIG_H
20#include "config.h"
21#endif
22
23MODULE radialsr
24
25  USE io_tools
26  USE GlobalMath
27  USE gridmod
28  USE atomdata
29
30  IMPLICIT NONE
31
32  !REAL(8), parameter :: inverse_fine_structure=137.035999139d0
33  !!! moved to globalmath
34  Real(8), private :: gamma,c1,c2,MA,MB,qq
35  Real(8), private, allocatable :: ww(:),jj(:)
36     ! jj stores (r+(alpha/2)**2*(E*r-rv) == r*M(r)
37     ! ww stores kappa*(kappa+1)/(r**2*M(r)) - (E - V(r))
38
39CONTAINS
40
41!******************************************************************
42! Allocate_scalar_relativistic
43!******************************************************************
44 Subroutine Allocate_scalar_relativistic(Grid)
45   Type(GridInfo), INTENT(IN) :: Grid
46
47   INTEGER :: n,i
48
49   n=Grid%n
50   allocate(ww(n),jj(n), stat=i)
51   if (i/=0) then
52      write(std_out,*)  'Allocate_scalar_relativistic: error in allocation ',i,n
53      stop
54   endif
55
56 End subroutine Allocate_scalar_relativistic
57
58!******************************************************************
59! Deallocate_scalar_relativistic
60!******************************************************************
61 Subroutine deallocate_scalar_relativistic
62
63    deallocate(ww,jj)
64
65 end subroutine deallocate_scalar_relativistic
66
67!*******************************************************************
68!  Subroutine Azeroexpand(Grid,Pot,l,energy)
69!      If finitenucleus==.true. assumes potential is non-singular
70!          at origin and Pot%v0 and Pot%v0p are properly set
71!      Otherwise, assumes nuclear potential is -2*Z/r
72!*******************************************************************
73 Subroutine Azeroexpand(Grid,Pot,l,energy,nr)
74   Type(GridInfo), INTENT(IN) :: Grid
75   Type(PotentialInfo), INTENT(IN) :: Pot
76   Integer, INTENT(IN) :: l
77   Real(8), INTENT(IN) :: energy
78   Integer, optional, INTENT(IN) :: nr
79
80   Integer :: i,j,k,n
81   Real(8) :: nz,xx,yy,angm,alpha2,balpha2
82   Real(8) :: Tm10,Tm11,T00,Tm21,Tm22,term
83
84   n=Grid%n
85   if (present(nr)) n=min(n,nr)
86
87!  check for possible ionic charge
88    n=Grid%n
89    qq=-Pot%rv(n)/2
90    if(qq<0.001d0) qq=0
91
92   nz=Pot%nz
93    ww=0; jj=0;
94   balpha2=inverse_fine_structure**2
95   alpha2=1.d0/balpha2
96   jj(1:n)=(Grid%r(1:n) + &
97&       0.25d0*alpha2*(energy*Grid%r(1:n)-Pot%rv(1:n)))
98   angm=l*(l+1)
99   ww(2:n)=(Pot%rv(2:n)/Grid%r(2:n)-energy) &
100&     + angm/(Grid%r(2:n)*jj(2:n))
101   ww(1)=0
102
103   if (.not.finitenucleus) then
104      gamma=sqrt(angm+1.d0-alpha2*nz**2)
105      term=1.d0+0.25d0*alpha2*(energy-Pot%v0)
106      Tm21=2*gamma+1;   Tm22=2*(2*gamma+2)
107      Tm10=nz*(2.d0+alpha2*(energy-Pot%v0))-(2*balpha2/nz)*term*(gamma-1.d0)
108      Tm11=nz*(2.d0+alpha2*(energy-Pot%v0))-(2*balpha2/nz)*term*(gamma)
109      T00=-alpha2*nz*Pot%v0p+term*(energy-Pot%v0) + &
110&        (Pot%v0p/nz+(4*balpha2**2/(nz*nz))*term**2)*(gamma-1.d0)
111      c1=-Tm10/Tm21
112      c2=-(Tm11*C1+T00)/Tm22
113      !write(std_out,*) 'Azeroexpand: ', gamma,c1,c2
114      MA=0; MB=0
115
116   else  ! version for finite nuclear size
117       gamma=l+1.d0
118       term=1.d0+0.25d0*alpha2*(energy-Pot%v0)
119       Tm21=2*l+2;      Tm22=2*(2*l+3)
120       Tm10=(0.25d0*alpha2*Pot%v0p/term)*(l)
121       Tm11=(0.25d0*alpha2*Pot%v0p/term)*(l+1)
122       T00=(energy-Pot%v0)*term+l*((0.25d0*alpha2*Pot%v0p/term)**2)
123       c1=-Tm10/Tm21
124       c2=-(Tm11*C1+T00)/Tm22
125       !write(std_out,*) 'Azeroexpand: ', gamma,c1,c2
126       MA=0; MB=0
127   endif
128
129  end subroutine Azeroexpand
130
131!*******************************************************************
132! SUBROUTINE wfnsrinit(Grid,l,wfn,lwfn,istart)
133!*******************************************************************
134  SUBROUTINE wfnsrinit(Grid,l,wfn,lwfn,istart)
135   ! returns the solution of the scalar relativistic equations near r=0
136   !  using power series expansion
137   Type(GridInfo), INTENT(IN) :: Grid
138   INTEGER, INTENT(IN) :: l
139    REAL(8),INTENT(INOUT) :: wfn(:),lwfn(:)
140    INTEGER, INTENT(OUT) :: istart
141
142    REAL(8) :: rr,M
143    INTEGER :: i,j,n
144
145    wfn=0; lwfn=0
146
147    istart=6
148    do i=1,istart
149       rr=Grid%r(i+1)
150       if (.not.finitenucleus) then
151          wfn(i+1)=1+rr*(c1+rr*c2)
152          lwfn(i+1)=(gamma-1)+rr*(c1*gamma+rr*c2*(gamma+1))
153          wfn(i+1)=wfn(i+1)*(rr**gamma)
154          lwfn(i+1)=lwfn(i+1)*(rr**gamma)/jj(i+1)
155
156       else   ! finite nucleus case
157          M=MA-MB*rr
158          wfn(i+1)=(1+rr*(c1+rr*c2))*(rr**(l+1))
159          lwfn(i+1)=(l+rr*((l+1)*c1+rr*(l+2)*c2))*(rr**(l+1))/M
160       endif
161
162    enddo
163
164  End SUBROUTINE wfnsrinit
165
166 subroutine wfnsrasym(Grid,wfn,lwfn,energy,iend)
167  ! returns the solution of the scalar relativistic equations near r=inf
168  !  using exp(-x*r) for upper component
169   Type(GridInfo), INTENT(IN) :: Grid
170    REAL(8),INTENT(INOUT) :: wfn(:),lwfn(:)
171    REAL(8), INTENT(IN) :: energy
172    INTEGER, INTENT(OUT) :: iend
173
174    REAL(8) :: rr,x,m,qx
175    INTEGER :: i,j,n
176
177    if (energy>0.d0) then
178       write(std_out,*) 'Error in wfnsrasym -- energy > 0', energy
179       stop
180    endif
181
182    wfn=0; lwfn=0;
183    n=Grid%n
184
185
186    m=1.d0+0.25d0*energy/(inverse_fine_structure**2)
187    x=sqrt(-m*energy)
188    qx=qq     !  Possible net ionic charge
189    qx=(qx/x)*(1.d0+0.5d0*energy/(inverse_fine_structure**2))
190    iend=5
191    do i=n-iend,n
192       wfn(i)=exp(-x*(Grid%r(i)-Grid%r(n-iend)))
193       if (qx>0.d0) then
194               rr=(Grid%r(i)/Grid%r(n-iend))**qx
195               wfn(i)=wfn(i)*rr
196       endif
197       lwfn(i)=-wfn(i)*(x*Grid%r(i)+(1.d0-qx))/jj(i)
198    enddo
199  end subroutine wfnsrasym
200
201  !**********************************************************************
202  !      subroutine unboundsr(Grid,Pot,nr,l,energy,wfn,nodes)
203  !  pgm to solve radial scalar relativistic equation for unbound states
204  !    at energy 'energy' and at angular momentum l
205  !
206  !    with potential rv/r, given in uniform linear or log mesh of n points
207  !   assuming p(r)=C*r**(l+1)*polynomial(r) for r==0;
208  !
209  !  nz=nuclear charge
210  !
211  !  Does not use Noumerov algorithm -- but uses coupled first-order
212  !       equations from David Vanderbilt, Marc Torrent, and Francois Jollet
213  !
214  ! also returns node == number of nodes for calculated state
215  !************************************************************************
216  SUBROUTINE unboundsr(Grid,Pot,nr,l,energy,wfn,nodes)
217    TYPE(GridInfo), INTENT(IN)  :: Grid
218    TYPE(PotentialInfo), INTENT(IN)  :: Pot
219    INTEGER, INTENT(IN) :: nr,l
220    REAL(8), INTENT(IN) :: energy
221    REAL(8), INTENT(INOUT) :: wfn(:)
222    INTEGER, INTENT(INOUT) :: nodes
223
224    INTEGER :: n,i,j,k,ierr,istart
225    REAL(8) :: scale
226    REAL(8), allocatable :: lwfn(:),zz(:,:,:),yy(:,:)
227
228    n=Grid%n
229    IF (nr > n) THEN
230       WRITE(STD_OUT,*) 'Error in unboundsr -- nr > n', nr,n
231       STOP
232    ENDIF
233
234    call Azeroexpand(Grid,Pot,l,energy,nr)
235
236    allocate(lwfn(nr),zz(2,2,nr),yy(2,nr),stat=ierr)
237       if (ierr/=0) then
238          write(std_out,*) ' allocation error in unboundsr ', nr,ierr
239          stop
240       endif
241
242    lwfn=0;zz=0;yy=0;
243
244    call wfnsrinit(Grid,l,wfn,lwfn,istart)
245    call prepareforcfdsol(Grid,1,istart,nr,wfn,lwfn,yy,zz)
246    call cfdsol(Grid,zz,yy,istart,nr)
247    call getwfnfromcfdsol(1,nr,yy,wfn)
248    nodes=countnodes(2,nr,wfn)
249    !
250    ! normalize to unity within integration range
251    !
252    !call filter(nr,wfn,machine_zero)     !NH votes to remove
253    scale=1.d0/overlap(Grid,wfn(1:nr),wfn(1:nr),1,nr)
254    scale=SIGN(SQRT(scale),wfn(nr-2))
255    wfn(1:nr)=wfn(1:nr)*scale
256
257    deallocate(lwfn,yy,zz)
258
259  END SUBROUTINE unboundsr
260
261!******************************************************************
262!  SUBROUTINE boundsr(Grid,Pot,eig,wfn,l,nroot,emin,ierr,success)
263!******************************************************************
264  SUBROUTINE boundsr(Grid,Pot,eig,wfn,l,nroot,emin,ierr,success)
265    !  pgm to solve radial scalar relativistic equation for nroot bound state
266    !    energies and wavefunctions for angular momentum l
267    !    with potential rv/r, given in uniform linear or log mesh of n points
268    !  nz=nuclear charge
269    !  emin=is estimate of lowest eigenvalue; used if nz=0
270    !     otherwise, set to the value of -(nz/(l+1))**2
271    !
272    !  It is assumed that the wavefunction has np-l-1 nodes, where
273    !    np is the principle quantum number-- np=1,2,..nroot
274    !
275    !  Does not use Noumerov algorithm -- but uses coupled first-order
276    !       equations from David Vanderbilt, Marc Torrent, and Francois Jollet
277    !
278    !  Corrections are also needed for r>n*h, depending on:
279    !         e0 (current guess of energy eigenvalue
280    !         the extrapolated value of rv == r * v
281    !
282    ! ierr=an nroot digit number indicating status of each root
283    !   a digit of 1 indicates success in converging root
284    !              2 indicates near success in converging root
285    !              9 indicates that root not found
286    !
287    ! first check how many roots expected =  ntroot (returned as argument)
288    !
289    TYPE(GridInfo), INTENT(IN) :: Grid
290    TYPE(PotentialInfo), INTENT(IN) :: Pot
291    REAL(8), INTENT(INOUT) :: eig(:),wfn(:,:)
292    INTEGER, INTENT(IN) :: l,nroot
293    INTEGER, INTENT(INOUT) :: ierr
294    REAL(8), INTENT(INOUT) :: emin
295    LOGICAL, INTENT(INOUT) :: success
296
297    REAL(8), PARAMETER :: convre=1.d-10,vlrg=1.d30
298    INTEGER, PARAMETER :: niter=1000
299
300    REAL(8), POINTER :: rv(:)
301    REAL(8), ALLOCATABLE :: p1(:),p2(:),dd(:)
302    INTEGER :: n
303    REAL(8) :: nz,h,v0,v0p
304    REAL(8) :: err,convrez,energy,zeroval
305    REAL(8) :: scale,emax,best,rout
306    REAL(8) :: arg,r,r2,veff,pppp1,rin,dele,x,rvp1,pnp1,bnp1
307    INTEGER :: iter,i,j,k,node,match,mxroot,ntroot,ir,iroot
308    INTEGER :: least,many,ifac,istart,iend
309    LOGICAL :: ok
310    REAL(8), allocatable :: lwfn(:),zz(:,:,:),yy(:,:)
311    ! integer :: icount=0
312
313    n=Grid%n
314    h=Grid%h
315    ALLOCATE(p1(n),p2(n),dd(n),stat=i)
316    IF (i/=0) THEN
317       WRITE(STD_OUT,*) ' Error in boundsr allocation ',i,n
318       STOP
319    ENDIF
320
321    success=.true.
322    allocate(lwfn(n),zz(2,2,n),yy(2,n),stat=i)
323       if (i/=0) then
324          write(std_out,*) ' allocation error in boundsr ', n,i
325          stop
326       endif
327
328
329    nz=Pot%nz
330    v0=Pot%v0
331    v0p=Pot%v0p
332    rv=>Pot%rv
333    err=n*nz*(h**4)
334    convrez=convre
335    IF (nz>0.001d0) convrez=convre*nz
336    ierr=0
337
338    WRITE(STD_OUT,*) 'z , l = ',nz,l
339    ! check how many roots expected by integration outward at
340    !   energy = 0
341    energy = 0
342    call Azeroexpand(Grid,Pot,l,energy)
343    lwfn=0;zz=0;yy=0;
344    call wfnsrinit(Grid,l,p1,lwfn,istart)
345    !
346    !start outward integration
347    call prepareforcfdsol(Grid,1,istart,n,p1,lwfn,yy,zz)
348    call cfdsol(Grid,zz,yy,istart,n)
349    call getwfnfromcfdsol(1,n,yy,p1)
350    node=countnodes(2,n,p1)
351
352    WRITE(STD_OUT,*) ' nodes at e=0  ', node
353
354    mxroot=node+1
355    ntroot=node
356    IF (mxroot.LT.nroot) THEN
357       WRITE(STD_OUT,*)'error in boundsr - for l = ',l
358       WRITE(STD_OUT,*) nroot,' states requested but only',mxroot,' possible'
359       DO ir=mxroot+1,nroot
360          ierr=ierr+9*(10**(ir-1))
361       ENDDO
362       success=.false.
363    ENDIF
364    mxroot=min0(mxroot,nroot)
365    !
366    IF (nz.EQ.0) energy=-ABS(emin)
367    IF (nz.NE.0) energy=-1.1d0*(nz/(l+1.d0))**2
368    emin=energy-err
369    emax=0.d0
370
371    DO iroot=1,mxroot
372       best=1.d10; dele=1.d10
373       energy=emin+err
374       IF (energy.LT.emin) energy=emin
375       IF (energy.GT.emax) energy=emax
376       ok=.FALSE.
377       !write(std_out,*) 'iter,iroot,energy',iter,iroot,energy
378       !write(std_out,*) 'emin,max',emin,emax
379       BigIter: DO iter=1,niter
380          !write(std_out,*) 'In iter with energy', iter,energy,niter,l,iroot
381          !  start inward integration
382          !  start integration at n
383          call Azeroexpand(Grid,Pot,l,energy)
384          ! find classical turning point
385          call ClassicalTurningPoint(Grid,Pot%rv,l,energy,match)
386          match=max(match,10); match=min(match,n-20)
387          call wfnsrasym(Grid,p2,lwfn,energy,iend)
388          call prepareforcfdsol(Grid,n-iend,n,n,p2,lwfn,yy,zz)
389          call cfdsol(Grid,zz,yy,n-iend,match)
390          call getwfnfromcfdsol(match,n,yy,p2)
391          match=match+6
392          rin=Gfirstderiv(Grid,match,p2)/p2(match)
393
394          call wfnsrinit(Grid,l,p1,lwfn,istart)
395          call prepareforcfdsol(Grid,1,istart,n,p1,lwfn,yy,zz)
396          call cfdsol(Grid,zz,yy,istart,match+6)
397          call getwfnfromcfdsol(1,match+6,yy,p1)
398          node= countnodes(2,match+6,p1)
399
400          !icount=icount+1
401          !  do i=1,match+6
402          !    write(100+icount,'(1p,2e15.7)') Grid%r(i),p1(i)
403          !  enddo
404          rout=Gfirstderiv(Grid,match,p1)/p1(match)
405          ! check whether node = (iroot-1)
406          !   not enough nodes -- raise energy
407          IF (node.LT.iroot-1) THEN
408             emin=MAX(emin,energy)-err
409             energy=emax-(emax-energy)*ranx()
410             ifac=9
411             !   too many nodes -- lower energy
412          ELSEIF (node.GT.iroot-1) THEN
413             IF (energy.LE.emin) THEN
414                ierr=ierr+9*(10**(iroot-1))
415                WRITE(STD_OUT,*) 'boundsr error -- emin too high',l,nz,emin,energy
416                do i=2,n
417                   write(999,'(1p,4e15.7)') Grid%r(i),jj(i)/Grid%r(i),ww(i),Pot%rv(i)
418                enddo
419                STOP
420             ENDIF
421             emax=MIN(emax,energy+err)
422             energy=emin+(energy-emin)*ranx()
423             !   correct number of nodes -- estimate correction
424          ELSEIF (node.EQ.iroot-1) THEN
425             DO j=1,match
426                p1(j)=p1(j)/p1(match)
427                         !write(std_out,*) 'j,p1',j,p1(j)
428             ENDDO
429             DO j=match,n
430                p1(j)=p2(j)/p2(match)
431                         !write(std_out,*) 'j,p2',j,p1(j)
432             ENDDO
433             scale=1.d0/overlap(Grid,p1,p1)
434             dele=(rout-rin)*scale
435                  !write(std_out,*) 'energy,dele,scale',energy,dele,scale
436             x=ABS(dele)
437             IF (x.LT.best) THEN
438                scale=SQRT(scale)
439                p1(1:n)=p1(1:n)*scale
440                call filter(n,p1,machine_zero)
441                wfn(1:n,iroot)=p1(1:n)
442                eig(iroot)=energy
443                !write(std_out,*) 'root',l,iroot,eig(iroot),emin,emax
444                best=x
445             ENDIF
446             IF (ABS(dele).LE.convrez) THEN
447                !write(std_out,*) 'iter with dele' , iter,dele
448                ok=.TRUE.
449                !  eigenvalue found
450                ierr=ierr+10**(iroot-1)
451                IF (iroot+1.LE.mxroot) THEN
452                   emin=energy+err
453                   emax=0
454                   energy=(emin+emax)/2
455                   IF (energy.LT.emin) energy=emin
456                   IF (energy.GT.emax) energy=emax
457                   best=1.d10
458                ENDIF
459                EXIT BigIter
460             ENDIF
461             IF (ABS(dele).GT.convrez) THEN
462                energy=energy+dele
463                ! if energy is out of range, pick random energy in correct range
464                IF (emin-energy.GT.convrez.OR.energy-emax.GT.convrez)         &
465&                    energy=emin+(emax-emin)*ranx()
466                ifac=2
467             ENDIF
468          ENDIF
469       ENDDO BigIter !iter
470       IF (.NOT.ok) THEN
471          success=.false.
472          ierr=ierr+ifac*(10**(iroot-1))
473          WRITE(STD_OUT,*) 'no convergence in boundsr',iroot,l,dele,energy
474          WRITE(STD_OUT,*) ' best guess of eig, dele = ',eig(iroot),best
475          IF (iroot.LT.mxroot) THEN
476             DO ir=iroot+1,mxroot
477                ierr=ierr+9*(10**(ir-1))
478             ENDDO
479          ENDIF
480        ! reset wfn with hydrogenic form
481        j=iroot+l+1
482        wfn(:,iroot)=0
483        x=(j)*sqrt(abs(eig(iroot)))
484        do i=2,n
485           wfn(i,iroot)=hwfn(x,j,l,Grid%r(i))
486        enddo
487       ENDIF
488    ENDDO !iroot
489
490    DEALLOCATE(p1,p2,dd,lwfn,yy,zz)
491             write(std_out,*) 'returning from boundsr -- ierr=',ierr
492  END SUBROUTINE Boundsr
493
494  subroutine prepareforcfdsol(Grid,i1,i2,n,wfn,lwfn,yy,zz)
495     Type(gridinfo), INTENT(IN) :: Grid
496     INTEGER, INTENT(IN) :: i1,i2,n
497     REAL(8), INTENT(IN) :: wfn(:),lwfn(:)
498     REAL(8), INTENT(OUT) :: yy(:,:),zz(:,:,:)
499
500     INTEGER :: i
501
502      yy=0;zz=0
503      yy(1,i1:i2)=wfn(i1:i2)
504      yy(2,i1:i2)=lwfn(i1:i2)
505
506      do  i=2,n
507       zz(1,1,i)=1.d0/Grid%r(i)
508       zz(1,2,i)=jj(i)/Grid%r(i)
509       zz(2,2,i)=-1.d0/Grid%r(i)
510       zz(2,1,i)=ww(i)
511      enddo
512
513   end subroutine prepareforcfdsol
514
515END MODULE radialsr
516