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