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