1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8!*********************************************************************** 9! This file arw.F contains modified versions of routines written by 10! A.R.Williams around 1985 (possibly based on previous work). 11! They are used for the solution of the radial Schrodinger's and 12! Poisson's equations in the atomic program. It also contains some 13! routines written by J.M.Soler in collaboration with A.R.Williams 14! and based on algorithms developed by ARW. 15! J. M. Soler, Nov. 1998 and April 2003 16!*********************************************************************** 17! The routines contained in this file are: 18! SUBROUTINE EGOFV 19! SUBROUTINE YOFE 20! SUBROUTINE NRMLZG 21! SUBROUTINE BCORGN 22! SUBROUTINE BCRMAX 23! SUBROUTINE NUMIN 24! SUBROUTINE NUMOUT 25! SUBROUTINE VHRTRE 26! SUBROUTINE QVLOFZ ---> moved to periodic_table.f 27! SUBROUTINE LMXOFZ ---> moved to periodic_table.f 28! SUBROUTINE CNFIG ---> moved to periodic_table.f 29!*********************************************************************** 30 31 SUBROUTINE EGOFV(H,S,NR,E,G,Y,L,Z,A,B,RMAX,NPRIN,NNODE,DGDR) 32!*********************************************************************** 33! Finds the radial atomic wavefunction and energy for a given potential, 34! number of nodes, and logarithmic derivative at the boundary. 35! Input: 36! real*8 H(NR) : Radial hamiltonian 37! real*8 S(NR) : Metric function. H and S are defined for a radial 38! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 39! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 40! For a logarithmic mesh (see below), they are 41! S(j) = (A*r(j))^2 42! H(j) = S(j)*V(r(j)) + A^2/4 43! where V(r) is the total effective radial potential 44! in Rydbergs, including the kinetic term l*(l+1)/r^2 45! integer NR : Number of radial points (including r(1)=0) 46! integer L : Angular momentum 47! real*8 Z : Atomic number 48! real*8 A,B : Log-mesh parameters: 49! r(j) = B*(exp(A*(j-1)) - 1), j=1,2,...,NR 50! real*8 RMAX : Maximum radius. Must be equal to r(NR) 51! integer NPRIN : Principal quantum number 52! integer NNODE : Required number of nodes for wavefunction 53! real*8 DGDR : Required logarithmic derivative dlogG/dr at RMAX 54! Output: 55! real*8 E : Solution energy 56! real*8 G(NR) : Radial wavefunction 57! Auxiliary: 58! real*8 Y(NR) : Intermediate array used in the Numerov method 59!*********************************************************************** 60 use parallel, only: ionode 61 use sys, only: die 62 use precision, only: dp 63 64 IMPLICIT NONE 65 INTEGER :: NR, L, NPRIN, NNODE 66 real(dp) :: H(NR), S(NR), E, G(NR), Y(NR), 67 & Z, A, B, RMAX, DGDR 68 69 INTEGER, PARAMETER :: MAXITER = 40 ! Max. iterations 70 real(dp) ,PARAMETER :: TOL = 1.D-5 ! Convergence tol. 71 INTEGER :: I, NCOR, NITER, NN, NN1, NN2 72 real(dp) :: DE, DE12, E1, E2, T 73 74 NCOR=NPRIN-L-1 75 NN=NNODE 76 NN1=NNODE 77 NN2=NNODE-1 78 E1=E 79 E2=E 80! Initial bracketing energy range 81 DE12=0.5D0 82 DE=0.D0 83 DO NITER = 1,MAXITER+1 84! New energy estimate from YOFE 85 E=E+DE 86 IF (E.GT.E1 .AND. E.LT.E2 .AND. 87 & NN.GE.NNODE-1 .AND. NN.LE.NNODE) THEN 88! New energy and number of nodes are within desired ranges 89 IF (ABS(DE).LT.TOL) EXIT ! NITER loop converged 90 ELSE 91! Use bisection to force energy within range 92 E=0.5D0*(E1+E2) 93 END IF 94! Find wavefunction Y for energy E, and new energy estimate E+DE 95 CALL YOFE(E, DE, DGDR, RMAX, H, S, Y, NR, L, NCOR, NN, Z, A, B) 96 IF (NN.LT.NNODE) THEN 97! Too few nodes. Increase lower energy bound 98 E1=E 99 NN1=NN 100 IF (NN2.LT.NNODE) THEN 101! Energy not yet bracketed. Increase higher energy bound 102 DE12=DE12*2.D0 103 E2=E1+DE12 104 END IF 105 ELSE 106! Too many nodes. Decrease higher energy bound 107 E2=E 108 NN2=NN 109 IF (NN1.GE.NNODE) THEN 110! Energy not yet bracketed. Decrease lower energy bound 111 DE12=DE12*2.D0 112 E1=E2-DE12 113 END IF 114 END IF 115 END DO ! NITER 116 117! No-convergence error exception 118 IF (NITER.GT.MAXITER) THEN 119 IF (IOnode) WRITE(6,'(A,/,A,F3.0,2(A,I2),2(A,F12.5))') 120 & ' EGOFV: ERROR: Too many iterations. Stopping.', 121 & ' Z=',Z,' L=',L,' NNODE=',NNODE,' E=',E,' DE=',DE 122 call die() 123 END IF 124 125! Find true waveftn G from auxiliary function Y and normalize 126 G(1) = 0.D0 127 DO I=2,NR 128 T=H(I)-E*S(I) 129 G(I)=Y(I)/(1.D0-T/12.D0) 130 END DO 131 CALL NRMLZG(G,S,NR) 132 END SUBROUTINE EGOFV 133 134 135 SUBROUTINE YOFE( E, DE, DGDR, RMAX, H, S, Y, NR, L, 136 & NCOR, NNODE, Z, A, B ) 137!*********************************************************************** 138! Integrates the radial Scroedinger equation for a given energy 139! Input: 140! real*8 E : Energy of the required wavefunction 141! real*8 DGDR : Required logarithmic derivative dlogG/dr at RMAX 142! real*8 RMAX : Maximum radius. Must be equal to r(NR) 143! real*8 H(NR) : Radial hamiltonian 144! real*8 S(NR) : Metric function. H and S are defined for a radial 145! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 146! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 147! For a logarithmic mesh (see below), they are 148! S(j) = (A*r(j))^2 149! H(j) = S(j)*V(r(j)) + A^2/4 150! where V(r) is the total effective radial potential 151! in Rydbergs, including the kinetic term l*(l+1)/r^2 152! integer NR : Number of radial points (including r(1)=0) 153! integer L : Angular momentum 154! integer NCOR : Number of core states of same L below the given one 155! integer NNODE : Required number of nodes for wavefunction 156! real*8 Z : Atomic number 157! real*8 A,B : Log-mesh parameters: 158! r(j) = B*(exp(A*(j-1)) - 1), j=1,2,...,NR 159! Output: 160! real*8 DE : Estimate of energy change required to eliminate kink 161! real*8 Y(NR) : Numerov function, related to G above by 162! Y(j)=G(j)*(12-H(j)+E*S(j))/12 163! integer NNODE : Actual number of nodes of wavefunction 164!*********************************************************************** 165 use precision, only: dp 166 IMPLICIT NONE 167 INTEGER, intent(in) :: NR, L, NCOR 168 INTEGER, intent(out) :: NNODE 169 real(dp) , intent(in) :: E, DGDR, RMAX, H(NR), S(NR), 170 & Z, A, B 171 real(dp) , intent(out) :: DE, Y(NR) 172 173 real(dp) ,PARAMETER:: TMAX =1.D0 ! Max negative kin engy 174 real(dp) ,PARAMETER:: DLGMAX=1.D3 ! Max log deriv at Rmax 175 INTEGER :: KNK, NR0, NNDIN 176 real(dp) :: GIN, GOUT, GSG, GSGIN, GSGOUT, RATIO, 177 & T, XIN, XOUT, Y2, YN, ZDR 178 179 ! Find where the negative kinetic energy H-ES becomes too large 180 ! and make Y=0 beyond that point 181 DO NR0 = NR,1,-1 182 IF ( H(NR0)-E*S(NR0) .LT. TMAX ) EXIT 183 Y(NR0)=0.D0 184 END DO ! NR0 185 186 ! Find boundary condition at origin 187 ZDR = Z*A*B 188 CALL BCORGN(E,H,S,L,ZDR,Y2) 189 190 ! Integrate Schroedinger equation outwards 191 KNK=NR0 ! A new value for the kink position KNK will be returned 192 CALL NUMOUT(E, H, S, Y, NCOR, KNK, NNODE, Y2, GOUT, GSGOUT, XOUT) 193 194 ! Find if kinetic energy is sufficiently non negative to use 195 ! Numerov at Rmax. Otherwise use zero-value boundary condition 196 IF (NR0.EQ.NR .AND. ABS(DGDR).LE.DLGMAX) THEN 197 CALL BCRMAX(E,DGDR,RMAX,H,S,NR0,YN,A,B) 198 ELSE 199 YN=0.D0 200 END IF 201 202 ! Integrate Schroedinger equation inwards 203 CALL NUMIN(E,H,S,Y,NR0,NNDIN,YN,GIN,GSGIN,XIN,KNK) 204 205 ! Make wavefunction continuous at R(KNK) 206 RATIO = GOUT/GIN 207 Y(KNK:NR0) = Y(KNK:NR0)*RATIO 208 XIN = XIN*RATIO 209 GSGIN = GSGIN*RATIO**2 210 211 ! Estimate the energy change required to eliminate kink 212 GSG=GSGOUT+GSGIN 213 T=H(KNK)-E*S(KNK) 214 DE=GOUT*(XOUT+XIN+T*GOUT)/GSG 215 216 ! Find the number of nodes 217 NNODE=NNODE+NNDIN 218 IF (DE.LT.0.D0) NNODE=NNODE+1 219 END SUBROUTINE YOFE 220 221 222 223 SUBROUTINE NRMLZG(G,S,N) 224!*********************************************************************** 225! Normalizes the radial wavefunction, using Simpson's rule for the norm 226! Input: 227! real*8 G(N) : Wavefunction to be normalized. G is related to the 228! true radial wavefunction psi by 229! G(j) = (dr/dj)^(3/2) * r(j) * psi(r(j)) 230! real*8 S(N) : Metric function defined for a logarithmic mesh 231! r(j) = B*(exp(A*(j-1)) - 1), j=1,2,...,N 232! as S(j) = (dr/dj)^2 = (A*r(j))^2 233! integer N : Number of radial points (including r(1)=0) 234! Output: 235! real*8 G(N) : Normalized wavefunction 236!*********************************************************************** 237 use alloc, only: re_alloc, de_alloc 238 use precision, only: dp 239 IMPLICIT NONE 240 INTEGER :: N 241 real(dp) :: G(N), S(N) 242 243 INTEGER :: I 244 real(dp) :: NORM, SRNRM 245 real(dp), POINTER :: gaux(:) 246 247 nullify(gaux) 248 call re_alloc( gaux, 1, N, name='gaux', routine='NRMLZG' ) 249 gaux = g*g 250 251 call integrator(gaux,s,n,norm) 252 253 call de_alloc( gaux, name='gaux', routine='NRMLZG' ) 254 255 ! Normalize wavefunction 256 SRNRM = SQRT(NORM) 257 DO I=1,N 258 G(I) = G(I)/SRNRM 259 END DO 260 261 END SUBROUTINE NRMLZG 262 263 264 SUBROUTINE BCORGN(E,H,S,L,ZDR,Y2) 265!************************************************************************ 266! Finds the boundary condition at the origin (actually the wavefunction 267! at the second point) to integrate the radial Schroedinger equation. 268! Input: 269! real*8 E : Energy of the required wavefunction 270! real*8 H(3) : Radial hamiltonian (only the first 3 points are used) 271! real*8 S(3) : Metric function. H and S are defined for a radial 272! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 273! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 274! For a logarithmic mesh (see below), they are 275! S(j) = (A*r(j))^2 276! H(j) = S(j)*V(r(j)) + A^2/4 277! where V(r) is the total effective radial potential 278! in Rydbergs, including the kinetic term l*(l+1)/r^2 279! integer L : Angular momentum 280! real*8 ZDR : Atomic charge times (dr/dj) at r=0 281! Output: 282! real*8 Y2 : Value at j=2 of the Numerov function, related to G 283! above by Y(j)=G(j)*(12-H(j)+E*S(j))/12 284! Notes: 285! - Local variable D(j) is the inverse of the diagonal of the 286! tri-diagonal Numerov matrix 287! - For L=0, G(j) vanishes at the origin (j=1), but its first and second 288! derivatives are finite, making Y(1) finite 289! - For L=1, G and G' vanish at the origin, but G'' and Y are finite 290! - For L>1, G, G', G'', and Y all vanish at the origin 291!************************************************************************ 292 use precision, only: dp 293 IMPLICIT NONE 294 INTEGER :: L 295 real(dp) :: E, H(3), S(3), ZDR, Y2 296 297 real(dp) :: C0, C1, C2, D2, T2, T3 298 299 T2=H(2)-E*S(2) 300 D2=-((24.D0+10.D0*T2)/(12.D0-T2)) 301 IF (L<2) THEN 302 IF (L==0) THEN 303 C0=ZDR/6.D0 304 C0=C0/(1.D0-0.75D0*ZDR) 305 ELSE 306 C0=1.D0/12.D0 307 C0=(-C0)*8.D0/3.D0 308 END IF 309 C1=C0*(12.D0+13.D0*T2)/(12.D0-T2) 310 T3=H(3)-E*S(3) 311 C2=(-5.D-1)*C0*(24.D0-T3)/(12.D0-T3) 312 D2=(D2-C1)/(1.D0-C2) 313 END IF 314 Y2=(-1.D0)/D2 315 316 END SUBROUTINE BCORGN 317 318 319 320 SUBROUTINE BCRMAX(E,DGDR,RMAX,H,S,N,YN,A,B) 321!************************************************************************ 322! Finds the boundary condition at Rmax (actually the wavefunction at the 323! last point) to integrate the radial Schroedinger equation inwards. 324! Input: 325! real*8 E : Energy of the required wavefunction 326! real*8 DGDR : Required logarithmic derivative dlogG/dr at RMAX 327! real*8 RMAX : Maximum radius. Must be equal to r(N) 328! real*8 H(N+1): Radial hamiltonian (only the last 3 points are used) 329! real*8 S(N+1): Metric function. H and S are defined for a radial 330! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 331! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 332! For a logarithmic mesh (see below), they are 333! S(j) = (A*r(j))^2 334! H(j) = S(j)*V(r(j)) + A^2/4 335! where V(r) is the total effective radial potential 336! in Rydbergs, including the kinetic term l*(l+1)/r^2 337! integer N : Number of radial points (including r(1)=0) 338! real*8 A,B : Log-mesh parameters: 339! r(j) = B*(exp(A*(j-1)) - 1), j=1,2,...,NR 340! Output: 341! real*8 YN : Value at j=N of the Numerov function, related to G 342! above by Y(j)=G(j)*(12-H(j)+E*S(j))/12 343!************************************************************************ 344 use precision, only: dp 345 IMPLICIT NONE 346 INTEGER :: N 347 real(dp) :: E, DGDR, RMAX, H(N+1), S(N+1), YN, A, B 348 349 real(dp) :: BETA, C1, C2, C3, DG, DN, TN, TNM1, TNP1 350 351 TNM1=H(N-1)-E*S(N-1) 352 TN =H(N )-E*S(N ) 353 TNP1=H(N+1)-E*S(N+1) 354 BETA=1.D0+B/RMAX 355 DG=A*BETA*(DGDR+1.D0-0.5D0/BETA) 356 357 C2=24.D0*DG/(12.D0-TN) 358 DN=-(24.D0+10.D0*TN)/(12.D0-TN) 359 360 C1= (12.D0-2.D0*TNM1)/(12.D0-TNM1) 361 C3=-(12.D0-2.D0*TNP1)/(12.D0-TNP1) 362 YN=-(1.D0-C1/C3)/(DN-C2/C3) 363 364 END SUBROUTINE BCRMAX 365 366 367 368 SUBROUTINE NUMIN(E,H,S,Y,NR,NNODE,YN,G,GSG,DY,KNK) 369!*********************************************************************** 370! Integrates Schroedinger's equation inwards, using Numerov's method 371! Input: 372! real*8 E : Energy of the required wavefunction 373! real*8 H(NR) : Radial hamiltonian 374! real*8 S(NR) : Metric function. H and S are defined for a radial 375! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 376! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 377! For a logarithmic mesh (see below), they are 378! S(j) = (A*r(j))^2 379! H(j) = S(j)*V(r(j)) + A^2/4 380! where V(r) is the total effective radial potential 381! in Rydbergs, including the kinetic term l*(l+1)/r^2 382! integer NR : Number of radial points (including r(1)=0) 383! real*8 YN : Value at j=N of the Numerov function, defined below 384! integer KNK : Fixes the 'kink point' r(KNK) up to which the 385! inward integration is required 386! Output: 387! real*8 Y(NR) : Numerov function related to G(j) by 388! Y(j)=(1-T(j)/12)*G(j), with T(j)=H(j)-E*S(j) 389! integer NNODE : Number of nodes of Y(j) between KNK and NR 390! real*8 G : Value of G, defined above, at r(KNK) 391! real*8 GSG : Probability G*S*G=(4*pi*r**2)*psi**2 at r(KNK) 392! real*8 DY : Value of Y(j)-Y(j-1) at j=KNKout 393! integer KNK : Made equal to NR-2 if the input value is larger 394! Algorithm: see routine NUMOUT 395!*********************************************************************** 396 use precision, only: dp 397 IMPLICIT NONE 398 INTEGER :: NR, NNODE, KNK 399 real(dp) :: E, H(NR), S(NR), Y(NR), YN, G, GSG, DY 400 401 INTEGER :: I 402 real(dp) :: T 403 404 Y(NR)=YN 405 T=H(NR)-E*S(NR) 406 G=Y(NR)/(1.D0-T/12.D0) 407 GSG=G*S(NR)*G 408 I=NR-1 409 Y(I)=1.D0 410 T=H(I)-E*S(I) 411 G=Y(I)/(1.D0-T/12.D0) 412 GSG=GSG+G*S(I)*G 413 DY=Y(I)-Y(NR) 414 NNODE=0 415 KNK = MIN(KNK,NR-2) 416 DO I = NR-2,KNK,-1 417 DY=DY+T*G 418 Y(I)=Y(I+1)+DY 419 IF( Y(I)*Y(I+1) .LT. 0.D0) NNODE=NNODE+1 420 T=H(I)-E*S(I) 421 G=Y(I)/(1.D0-T/12.D0) 422 GSG=GSG+G*S(I)*G 423 END DO 424 425 END SUBROUTINE NUMIN 426 427 428 SUBROUTINE NUMOUT( E, H, S, Y, NCOR, KNK, NNODE, Y2, G, GSG, DY ) 429!*********************************************************************** 430! Integrates Schroedinger's equation outwards, using Numerov's method, 431! up to the first maximum of psi after NCOR nodes 432! Input: 433! real*8 E : Energy of the required wavefunction 434! real*8 H(KNK) : Radial hamiltonian 435! real*8 S(KNK) : Metric function. H and S are defined for a radial 436! variable mesh so that d2G/dj2=(H(j)-E*S(j))*G(j) 437! where G(j)=(dr/dj)^(3/2)*r(j)*psi(r(j)) 438! For a logarithmic mesh (see below), they are 439! S(j) = (A*r(j))^2 440! H(j) = S(j)*V(r(j)) + A^2/4 441! where V(r) is the total effective radial potential 442! in Rydbergs, including the kinetic term l*(l+1)/r^2 443! integer NCOR : Number of core states of same L below the given one 444! integer KNK : Number of radial points (including r(1)=0) 445! real*8 Y2 : Value at j=N of the Numerov function, defined below 446! Output: 447! real*8 Y(KNK) : Numerov function related to G(j) by 448! Y(j)=(1-T(j)/12)*G(j), with T(j)=H(j)-E*S(j) 449! integer NNODE : One plus the number of nodes up to output KNK 450! real*8 G : Value of G, defined above, at r(KNKout) 451! real*8 GSG : Probability G*S*G=(4*pi*r**2)*psi**2 at r(KNKout) 452! real*8 DY : Value of Y(j)-Y(j-1) at j=KNKout 453! integer KNK : Position of first maximum of Y(j) after NCOR nodes, 454! or KNK=KNKin-4 (whatever is lower) 455! Algorithm: 456! The Numerov equation is 457! (g(j+1) - 2*g(j) + g(j-1)) / dx2 = 458! (t(j+1)*g(j+1) + 10*t(j)*g(j) + t(j-1)*g(j-1)) / 12 459! or, for dx=1, 460! (1-t(j+1)/12)*g(j+1) = (2+10*t(j)/12)*g(j) - ((1-t(j-1)/12)*g(j-1) 461! Then, defining y(j)=(1-t(j)/12)*g(j) and after some simple algebra 462! dy(j) = y(j+1)-y(j) = t(j)*g(j) + y(j)-y(j-1) = t(j)*g(j) + dy(j-1) 463!*********************************************************************** 464 use precision, only: dp 465 IMPLICIT NONE 466! Input/Output variables 467 INTEGER, intent(in) :: NCOR 468 INTEGER, intent(inout) :: KNK 469 INTEGER, intent(out) :: NNODE 470 real(dp) , intent(in) :: E, H(KNK), S(KNK), Y2 471 real(dp) , intent(out) :: Y(KNK), G, GSG, DY 472! Local variables 473 INTEGER :: I 474 real(dp) :: DYL, T 475 476 Y(1) = 0.D0 477 Y(2) = Y2 478 T = H(2)-E*S(2) 479 G = Y(2)/(1.D0-T/12.D0) 480 GSG = G*S(2)*G 481 Y(3) = 1.D0 482 T = H(3)-E*S(3) 483 G = Y(3)/(1.D0-T/12.D0) 484 GSG = GSG+G*S(3)*G 485 DY = Y(3)-Y(2) 486 NNODE = 0 487 DO I= 4, KNK-4 488 DYL = DY 489 DY = DY+T*G 490 Y(I) = Y(I-1)+DY 491 IF( Y(I)*Y(I-1) .LT. 0.D0) NNODE = NNODE + 1 492 T = H(I)-E*S(I) 493 G = Y(I)/(1.D0-T/12.D0) 494 GSG = GSG+G*S(I)*G 495 ! End loop if |Y(j)| has a maximum after NCOR nodes 496 IF (NNODE.GE.NCOR .AND. DYL*DY.LE.0.D0) EXIT 497 END DO 498 KNK = MIN(I,KNK-4) 499 500 END SUBROUTINE NUMOUT 501 502 503 504 SUBROUTINE VHRTRE(R2RHO,V,R,DRDI,SRDRDI,NR,A) 505!*********************************************************************** 506! Finds the Hartree potential created by a radial electron density, 507! using Numerov's method to integrate the radial Poisson equation: 508! d2(r*V)/dr2 = -4*pi*rho*r = -(4*pi*r2*rho)/r 509! Input: 510! real*8 R2RHO(NR) : 4*pi*r**2*rho, with rho the electron density 511! real*8 R(NR) : Logarithmic radial mesh R(i)=B*(exp(A*(i-1)-1) 512! real*8 DRDI(NR) : dr/di at the mesh points 513! real*8 SRDRDI(NR): sqrt(dr/di) at the mesh points 514! integer NR : Number of radial mesh points, including r(1)=0 515! real*8 A : The parameter A in r(i)=B*(exp(A*(i-1)-1) 516! Output: 517! real*8 V(NR) : Electrostatic potential created by rho, in Ryd 518! The constants of integration are fixed so that 519! V=finite at the origin and V(NR)=Q/R(NR), 520! where Q is the integral of rho up to R(NR) 521! Algorithm: see routine NUMOUT 522!*********************************************************************** 523 use precision, only: dp 524 use alloc, only: re_alloc, de_alloc 525 IMPLICIT NONE 526 INTEGER :: NR 527 real(dp) :: R2RHO(NR),V(NR),R(NR),DRDI(NR),SRDRDI(NR),A 528 529 INTEGER :: IR 530 real(dp) :: A2BY4, BETA, DV, DY, DZ, Q, QBYY, QPARTC, QT, 531 . T, V0, Y, YBYQ 532 real(dp), POINTER :: gaux(:) 533 534 535 ! Find some constants 536 A2BY4=A*A/4.D0 537 YBYQ=1.D0-A*A/48.D0 538 QBYY=1.D0/YBYQ 539 540 ! Use Simpson's rule to find the total charge QT, and the 541 ! potential at the origin V0: 542 ! QT = Int(4*pi*r**2*rho*dr) = Int((4*pi*r**2*rho)*(dr/di)*di) 543 ! V0 = Int(4*pi*r*rho*dr) = Int((4*pi*r**2*rho)/r*(dr/di)*di) 544 545 call integrator(r2rho(2),drdi(2),nr-1,QT) 546 547 nullify(gaux) 548 call re_alloc( gaux, 2, NR, name='gaux', routine='VHRTRE' ) 549 550 gaux = r2rho(2:nr)/r(2:nr) 551 call integrator( gaux, drdi(2), nr-1, V0 ) 552 553 call de_alloc( gaux, name='gaux', routine='VHRTRE' ) 554 555 ! Fix V(1) and V(2) to start Numerov integration. To find a 556 ! particular solution of the inhomog. eqn, V(2) is fixed by 557 ! setting rV(2)=0. Notice that V=finite => rV=0 at r=0 558 V(1)=2.D0*V0 ! Factor 2 because we use Rydbergs 559 T=SRDRDI(2)/R(2) 560 BETA=DRDI(2)*T*R2RHO(2) 561 DY=0.D0 562 Y=0.D0 563 Q=(Y-BETA/12.D0)*QBYY 564 V(2)=2.D0*T*Q 565 566 ! Integrate Poisson's equation outwards, using Numerov's method 567 DO IR = 3,NR 568 DY=DY+A2BY4*Q-BETA 569 Y=Y+DY 570 T=SRDRDI(IR)/R(IR) 571 BETA=T*DRDI(IR)*R2RHO(IR) 572 Q=(Y-BETA/12.D0)*QBYY 573 V(IR)=2.D0*T*Q 574 END DO 575 576 ! Add a solution (finite at r=0) of the homogeneous equation 577 ! d2(r*V)/dr2=0 => rV=const*r => V=const, to ensure that 578 ! V(NR)=Q/R(NR). Notice that V(1) is set independently 579 QPARTC = R(NR)*V(NR)/2.D0 580 DZ=QT-QPARTC 581 DV=2.D0*DZ/R(NR) 582 DO IR=2,NR 583 V(IR)=V(IR)+DV 584 ENDDO 585 586 END SUBROUTINE VHRTRE 587 588 SUBROUTINE INTEGRATOR(F,S,NP,VAL) 589!*********************************************************************** 590! Integrates a radial function tabulated on a logarithmic grid, 591! using a generalized Simpson's rule valid for both even and odd 592! number of points. Note that the "h" is 1 as the reparametrization 593! involves a mapping of integers to reals. 594! 595! Alberto Garcia, Dec. 2006, based on code by Art Williams. 596! 597! Input: 598! real*8 F(NP) : Function to be integrated. 599! real*8 S(NP) : Metric function defined for a logarithmic mesh 600! r(j) = B*(exp(A*(j-1)) - 1), j=1,2,...,NP 601! as S(j) = (dr/dj)^2 = (A*r(j))^2 602! integer NP : Number of radial points (including r(1)=0) 603! Output: 604! real*8 VAL : Value of the integral 605!*********************************************************************** 606 use precision, only: dp 607 IMPLICIT NONE 608 INTEGER :: NP 609 real(dp) :: F(NP), S(NP) 610 611 INTEGER :: I, N 612 real(dp) :: VAL 613 614 IF (MOD(NP,2).EQ.1) THEN 615 N = NP ! ODD 616 ELSE 617 IF (NP .EQ. 2) THEN 618 ! Special case of trapezoidal rule 619 VAL = 0.5D0 * (F(1)*S(1) + F(2)*S(2)) 620 RETURN 621 ENDIF 622 N = NP - 3 ! EVEN: TAKE A FINAL FOUR-POINT INTERVAL 623 ENDIF 624! 625! STANDARD EXTENDED SIMPSON RULE WITH ALTERNATING 4/3 AND 2/3 FACTORS 626! FOR THE SECTION MADE UP OF PAIRS OF INTERVALS (THREE-POINT SEGMENTS) 627! 628 VAL = 0.D0 629 DO I = 2,N-1,2 630 VAL=VAL + F(I)*S(I) 631 END DO 632 VAL = VAL * 2.D0 633 DO I = 3,N-2,2 634 VAL=VAL + F(I)*S(I) 635 END DO 636 VAL = VAL * 2.D0 637 DO I = 1,N,N-1 ! first and last points at 1/3 638 VAL=VAL + F(I)*S(I) 639 END DO 640 VAL = VAL/3.D0 641 642 IF (MOD(NP,2).EQ.0) THEN ! EVEN 643! ADD THE CONTRIBUTION OF THE 644! FINAL FOUR-POINT SEGMENT USING SIMPSON'S 3/8 RULE 645! (SEE NUMERICAL RECIPES (FORTRAN), P. 105) 646 I = NP - 3 647 VAL = VAL + 648 $ (3.D0/8.D0) * ( (F(I)*S(I) + F(I+3)*S(I+3)) + 649 $ 3.D0 * (F(I+1)*S(I+1) + F(I+2)*S(I+2)) ) 650 ENDIF 651 END SUBROUTINE INTEGRATOR 652