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