1! 2! Copyright (C) 2002 Vanderbilt group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! This routine has been modified in order to be compatible with the 9! ld1 code. The numerical algorithm is unchanged. 10! ADC Nov 2003 11! 12!------------------------------------------------------------------------- 13! 14subroutine dirsol(idim1,mesh,ncur,lcur,jcur,it,e0,thresh,grid,snl,ruae,nstop) 15! 16! subroutine to compute solutions to the full dirac equation 17! 18! the dirac equation in rydberg units reads: 19! 20! df(r) k alpha 21! ----- = + - f(r) - ( e - v(r) ) ----- g(r) 22! dr r 2 23! 24! dg(r) k 4 alpha 25! ----- = - - g(r) + ( e - v(r) + -------- ) ----- f(r) 26! dr r alpha**2 2 27! 28! where 29! alpha is the fine structure constant 30! f(r) is r*minor component 31! g(r) is r*major component 32! k is quantum number 33! k = - (l+1) if j = l+0.5 34! k = + l if j = l-0.5 35! IMPORTANT: on output, snl(:,1) contains the MAJOR component 36! snl(:,2) contains the MINOR component 37! 38!---------------------------------------------------------------------------- 39! 40use io_global, only : stdout 41use kinds, only : DP 42use radial_grids, only: radial_grid_type 43use ld1inc, only : cau_fact, zed 44implicit none 45integer :: idim1 46type(radial_grid_type),intent(in)::grid 47real(DP) :: ruae(idim1), & ! the all electron potential 48 snl(idim1,2) ! the wavefunction 49 50real(DP) :: e0, & ! the starting energy eigenvalue 51 jcur, & ! the j of the state 52 thresh 53 54integer :: mesh, & ! the dimension of the mesh 55 it, & ! the iteration 56 ncur, & ! the n of the state 57 lcur ! the l of the state 58 59real(DP) :: tbya, abyt, & 60 emin, emax, & 61 zz(idim1,2,2), & 62 tolinf,alpha2,alpha, & 63 yy(idim1,2), & 64 vzero, & 65 f0,f1,f2,g0,g1,g2, & 66 gout, gpout, & 67 gin, gpin, & 68 factor,gamma0, & 69 ecur, & 70 decur, decurp 71real(DP) :: f(idim1), int_0_inf_dr 72 73integer :: itmax, & ! maximum number of iterations 74 iter, & ! current iteration 75 ir, & ! counter 76 ig, & ! auxiliary 77 kcur, & ! current k 78 nctp, & ! index of the classical turning point 79 nodes, & ! the number of nodes 80 ninf, & ! practical infinite 81 nstop ! 0 if all ok, 1 otherwise 82! 83! r o u t i n e i n i t i a l i s a t i o n 84if (mesh.ne.grid%mesh) call errore('dirsol','mesh dimension is not as expected',1) 85nstop=0 86! 87! set the maximum number of iterations for improving wavefunctions 88! 89itmax = 100 90! 91! set ( 2 / fine structure constant ) 92!tbya = 2.0_DP * 137.04_DP 93tbya = 2.0_DP * cau_fact 94 95! set ( fine structure constant / 2 ) 96abyt = 1.0_DP / tbya 97 98if (jcur.eq.lcur+0.5_DP) then 99 kcur = - ( lcur + 1 ) 100else 101 kcur = lcur 102endif 103! 104! set initial upper and lower bounds for the eigen value 105emin = - 1.0e10_DP 106emax = 1.0_DP 107ecur=e0 108! 109do iter = 1,itmax 110 yy = 0.0_DP 111! 112! define the zz array 113! =================== 114! 115 if ( iter .eq. 1 ) then 116 do ir = 1,mesh 117 zz(ir,1,1) = grid%rab(ir) * DBLE(kcur) / grid%r(ir) 118 zz(ir,2,2) = - zz(ir,1,1) 119 enddo 120 endif 121 do ir = 1,mesh 122 zz(ir,1,2) = - grid%rab(ir) * ( ecur - ruae(ir) ) * abyt 123 zz(ir,2,1) = - zz(ir,1,2) + grid%rab(ir) * tbya 124 enddo 125! 126! ============================================== 127! classical turning point and practical infinity 128! ============================================== 129! 130 do nctp = mesh,10,-1 131 if ( zz(nctp,1,2) .lt. 0.0_DP ) goto 240 132 enddo 133 call errore('dirsol', 'no classical turning point found', 1) 134! 135! jump point out of classical turning point loop 136240 continue 137 138 if ( nctp .gt. mesh - 10 ) then 139! write(stdout,*) 'State nlk=', ncur, lcur, kcur, nctp, mesh 140! write(stdout,*) 'ecur, ecurmax=', ecur, ruae(mesh-10) 141 write(stdout,*) 'classical turning point too close to mesh',ncur,lcur,kcur 142 snl=0.0_DP 143 e0=e0-0.2_DP 144 nstop=1 145 goto 700 146 endif 147! 148 tolinf = log(thresh) ** 2 149 do ninf = nctp+10,mesh 150 alpha2 = (ruae(ninf)-ecur) * (grid%r(ninf) - grid%r(nctp))**2 151 if ( alpha2 .gt. tolinf ) goto 260 152 enddo 153! 154! jump point out of practical infinity loop 155260 continue 156! 157 if (ninf.gt.mesh) ninf=mesh 158! 159! =========================================================== 160! analytic start up of minor and major components from origin 161! =========================================================== 162! 163! with finite nucleus so potential constant at origin we have 164! 165! f(r) = sum_n f_n r ** ( ig + n ) 166! g(r) = sum_n g_n r ** ( ig + n ) 167! 168! with 169! 170! f_n+1 = - (ecur-v(0)) * abyt * g_n / ( ig - kcur + n + 1 ) 171! g_n+1 = (ecur-v(0)+tbya**2 ) * abyt * f_n / ( ig + kcur + n + 1) 172! 173! if kcur > 0 ig = + kcur , f_0 = 1 , g_0 = 0 174! if kcur < 0 ig = - kcur , f_0 = 0 , g_1 = 1 175! 176! Uncomment the following instruction if you need this boundary conditions 177! 178! vzero = ruae(1) 179! 180! set f0 and g0 181! if ( kcur .lt. 0 ) then 182! ig = - kcur 183! f0 = 0 184! g0 = 1 185! else 186! ig = kcur 187! f0 = 1 188! g0 = 0 189! endif 190! 191! f1 = - (ecur-vzero) * abyt * g0 / DBLE( ig - kcur + 1 ) 192! g1 = (ecur-vzero+tbya**2) * abyt * f0 / DBLE( ig + kcur + 1 ) 193! f2 = - (ecur-vzero) * abyt * g1 / DBLE( ig - kcur + 2 ) 194! g2 = (ecur-vzero+tbya**2) * abyt * f1 / DBLE( ig + kcur + 2 ) 195! 196! 197! do ir = 1,5 198! yy(ir,1) = grid%r(ir)**ig * ( f0 + grid%r(ir) * ( f1 + grid%r(ir) * f2 ) ) 199! yy(ir,2) = grid%r(ir)**ig * ( g0 + grid%r(ir) * ( g1 + grid%r(ir) * g2 ) ) 200! enddo 201 202! 203! The following boundary conditions are for the Coulomb nuclear potential 204! ADC 05/10/2007 205! 206 vzero = ruae(1)+zed*2.0_DP/grid%r(1) 207! 208 gamma0=sqrt(kcur**2-4.0_DP*(abyt*zed)**2) 209 if ( kcur .lt. 0 ) then 210 ig = - kcur 211 f0 = (kcur+gamma0)/(2.0_DP*abyt*zed) 212 g0 = 1.0_DP 213 else 214 ig = kcur 215 f0 = 1.0_DP 216 g0 = (kcur-gamma0)/(2.0_DP*abyt*zed) 217 endif 218! 219 f1 = - (ecur-vzero) * abyt * g0 / ( gamma0 - kcur + 1.0_DP ) 220 g1 = (ecur-vzero+tbya**2) * abyt * f0 / ( gamma0 + kcur + 1.0_DP ) 221 f2 = - (ecur-vzero) * abyt * g1 / ( gamma0 - kcur + 2.0_DP ) 222 g2 = (ecur-vzero+tbya**2) * abyt * f1 / ( gamma0 + kcur + 2.0_DP ) 223! 224! 225 do ir = 1,5 226 yy(ir,1) = grid%r(ir)**gamma0*(f0+grid%r(ir)*(f1+grid%r(ir)*f2)) 227 yy(ir,2) = grid%r(ir)**gamma0*(g0+grid%r(ir)*(g1+grid%r(ir)*g2)) 228 enddo 229 230! =========================== 231! outward integration to nctp 232! =========================== 233! 234! fifth order predictor corrector integration routine 235 call cfdsol(zz,yy,6,nctp,idim1) 236! 237! save major component and its gradient at nctp 238 gout = yy(nctp,2) 239 gpout = zz(nctp,2,1)*yy(nctp,1) + zz(nctp,2,2)*yy(nctp,2) 240 gpout = gpout / grid%rab(nctp) 241! 242! ============================================== 243! start up of wavefunction at practical infinity 244! ============================================== 245! 246 do ir = ninf,ninf-4,-1 247 alpha = sqrt( ruae(ir) - ecur ) 248 yy(ir,2) = exp ( - alpha * ( grid%r(ir) - grid%r(nctp) ) ) 249 yy(ir,1) = ( DBLE(kcur)/grid%r(ir) - alpha ) * yy(ir,2)*tbya / & 250 & ( ecur - ruae(ir) + tbya ** 2 ) 251 enddo 252! 253! ========================== 254! inward integration to nctp 255! ========================== 256! 257! fifth order predictor corrector integration routine 258 call cfdsol(zz,yy,ninf-5,nctp,idim1) 259! 260! save major component and its gradient at nctp 261 gin = yy(nctp,2) 262 gpin = zz(nctp,2,1)*yy(nctp,1) + zz(nctp,2,2)*yy(nctp,2) 263 gpin = gpin / grid%rab(nctp) 264! 265! 266! =============================================== 267! rescale tail to make major component continuous 268! =============================================== 269! 270 factor = gout / gin 271 do ir = nctp,ninf 272 yy(ir,1) = factor * yy(ir,1) 273 yy(ir,2) = factor * yy(ir,2) 274 enddo 275! 276 gpin = gpin * factor 277! 278! ================================= 279! check that the number of nodes ok 280! ================================= 281! 282! count the number of nodes in major component 283 call nodeno(yy(1,2),1,ninf,nodes,idim1) 284 285 if ( nodes .lt. ncur - lcur - 1 ) then 286! energy is too low 287 emin = ecur 288! write(stdout,*) 'energy too low' 289! write(stdout,'(i5,3f12.5,2i5)') & 290! & iter,emin,ecur,emax,nodes,ncur-lcur-1 291 if ( ecur * 0.9_DP .gt. emax ) then 292 ecur = 0.5_DP * ecur + 0.5_DP * emax 293 else 294 ecur = 0.9_DP * ecur 295 endif 296 goto 370 297 endif 298! 299 if ( nodes .gt. ncur - lcur - 1 ) then 300! energy is too high 301 emax = ecur 302! 303! write(stdout,*) 'energy too high' 304! write(stdout,'(i5,3f12.5,2i5)') & 305! & iter,emin,ecur,emax,nodes,ncur-lcur-1 306 if ( ecur * 1.1_DP .lt. emin ) then 307 ecur = 0.5_DP * ecur + 0.5_DP * emin 308 else 309 ecur = 1.1_DP * ecur 310 endif 311 goto 370 312 endif 313! 314! 315! ======================================================= 316! find normalisation of wavefunction 317! ======================================================= 318! 319 do ir = 1,ninf 320 f(ir) = (yy(ir,1)**2 + yy(ir,2)**2) 321 enddo 322 factor=int_0_inf_dr(f,grid,ninf,2*ig) 323! 324! 325! ========================================= 326! variational improvement of the eigenvalue 327! ========================================= 328! 329 decur = gout * ( gpout - gpin ) / factor 330! 331! to prevent convergence problems: 332! do not allow decur to exceed 20% of | ecur | 333! do not allow decur to exceed 70% of distance to emin or emax 334 if (decur.gt.0.0_DP) then 335 emin=ecur 336 decurp=min(decur,-0.2_DP*ecur,0.7_DP*(emax-ecur)) 337 else 338 emax=ecur 339 decurp=-min(-decur,-0.2_DP*ecur,0.7_DP*(ecur-emin)) 340 endif 341! 342! write(stdout,'(i5,3f12.5,1p2e12.4)') & 343! & iter,emin,ecur,emax,decur,decurp 344! 345! test to see whether eigenvalue converged 346 if ( abs(decur) .lt. thresh ) goto 400 347 348 ecur = ecur + decurp 349! 350! jump point from node check 351370 continue 352! 353! ======================================================= 354! check that the iterative loop is not about to terminate 355! ======================================================= 356! 357 if ( iter .eq. itmax ) then 358! eigenfunction has not converged in allowed number of iterations 359 360! write(stdout,999) it,ncur,lcur,jcur,e0,ecur 361!999 format('iter',i4,' state',i4,i4,f4.1,' could not be converged.',/, & 362! & ' starting energy for calculation was',f10.5, & 363! & ' and end value =',f10.5) 364 write(stdout,*) 'state nlj',ncur,lcur,jcur, ' not converged' 365 snl=0.0_DP 366 nstop=1 367 goto 700 368 endif 369! 370! close iterative loop 371enddo 372! 373! jump point on successful convergence of eigenvalue 374400 continue 375! 376! normalize the wavefunction and exit. Note also that in this routine 377! yy(.,1) is the small component 378! yy(.,2) is the large component 379! 380! 381snl=0.0_DP 382do ir=1,mesh 383 snl(ir,1)=yy(ir,2)/sqrt(factor) 384 snl(ir,2)=yy(ir,1)/sqrt(factor) 385enddo 386e0=ecur 387700 continue 388return 389end subroutine dirsol 390