1! 2! Copyright (C) 2004 PWSCF 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#undef DEBUG 9!--------------------------------------------------------------- 10subroutine c6_tfvw (mesh, zed, grid, rho_input) 11 !-------------------------------------------------------------------- 12 ! 13 use kinds, only : DP 14 use constants, only : e2, pi, fpi, BOHR_RADIUS_ANGS 15 use ld1inc, only : lsd, nwf, oc, ll 16 use radial_grids, only: radial_grid_type 17 ! 18 implicit none 19 ! 20 ! I/O variables 21 ! 22 type(radial_grid_type), intent(in):: grid 23 integer mesh 24 real (kind=8) :: rho_input(mesh) 25 real (kind=8) :: zed, rho(mesh) 26 ! 27 !real (kind=8) :: vw_lambda=9.0_dp ! describes exactly the low-q limit 28 !real (kind=8) :: vw_lambda=3.0_dp ! interpolate between high and low q's 29 real (kind=8) :: vw_lambda=1.0_dp ! describes exactly the high-q limit 30 ! 31 ! local variables 32 ! 33 logical :: csi, l_add_tf_term 34 real (kind=8) :: error, error2, e, charge, beta, u, alpha, dalpha, c6, du1, & 35 du2, factor, ze2, thresh 36 real (kind=8), allocatable :: veff(:), y(:), yy(:) 37 real (kind=8), allocatable :: dvpot(:), dvscf(:), drho(:), dvhx(:), dvxc(:), pp(:) 38 complex (kind=8), allocatable :: dy(:), drho_old(:) 39 integer i, iter, n, l, ly, iu, Nu, Nc, counter, nstop, mesh_save 40 41 allocate ( veff(mesh),y(mesh),yy(mesh)) 42 allocate ( dvpot(mesh),dvscf(mesh),drho(mesh),dvhx(mesh),dvxc(mesh),pp(mesh) ) 43 allocate ( dy(mesh), drho_old(mesh) ) 44 ! 45 write(6,'(/,/,/,5x,10(''-''),'' Compute C6 from polarizability with TFvW approx.'',10(''-''),/)') 46 if (vw_lambda.ne.1.d0) write(6,*) " value of vw_lambda ", vw_lambda 47 ! 48 if (mesh.ne.grid%mesh) call errore('c6_tfwv',' mesh dimension is not as expected',1) 49 do i = 1, mesh 50 rho(i) = rho_input(i) / (fpi*grid%r(i)**2) 51 end do 52 ! 53 counter = 1 54 do i = 1, mesh 55 if (rho(i) .gt. 1.0d-30) counter = counter + 1 56 enddo 57 mesh_save = mesh 58 mesh = counter 59#ifdef DEBUG 60 write (*,*) "mesh ", mesh 61#endif 62 ! 63 if (lsd .ne. 0) call errore ('c6_tfvw', 'implemented only for non-magnetic ions', lsd) 64 csi = .true. 65 do i = 1, nwf 66 csi = csi .and. ( ((ll(i).eq.0) .and. (oc(i).eq.2 )) .or. & 67 ((ll(i).eq.1) .and. (oc(i).eq.6 )) .or. & 68 ((ll(i).eq.2) .and. (oc(i).eq.10)) .or. & 69 ((ll(i).eq.3) .and. (oc(i).eq.14)) ) 70 enddo 71 if (.not. csi) call errore ('c6_tfvw', 'implemented only for closed-shell ions', 1) 72! rho = 0.d0 73! open (7,file='rho.out',status='unknown',form='formatted') 74! do i=1,mesh 75! read (7,'(P5E20.12)') r(i), rho(i), y(i), y(i), y(i) 76! write (6,'(P5E15.6)') r(i), rho(i) 77! end do 78! close (7) 79! 80! compute unperturbed effective potential 81! 82 call veff_of_rho(mesh,grid%dx,grid%r,grid%r2,rho,y,veff) 83#ifdef DEBUG 84 write (*,*) "veff(1:3)" 85 write (*,*) veff(1:3) 86 write (*,*) "veff(mesh-5:mesh)" 87 write (*,*) veff(mesh-5:mesh) 88#endif 89! 90! check that veff and y are what we think 91! 92 n = 1 93 l = 0 94 e = -1.d-7 95 charge=0.d0 96 ze2 = - zed * e2 97 thresh = 1.d-14 98 99 do i=1,mesh 100 charge = charge + rho(i) * fpi * grid%r2(i) * grid%r(i) * grid%dx 101 end do 102! call solve_scheq(n,l,e,mesh,dx,r,sqr,r2,veff,zed,yy) 103 call ascheq (n, l, e, mesh, grid, veff, ze2, thresh, yy, nstop) 104 105 error = 0.d0 106 do i=1,mesh 107 error = error + (y(i)-yy(i)/grid%sqr(i)*sqrt(charge))**2 * grid%r2(i) * grid%dx 108 end do 109 110#ifdef DEBUG 111 write (*,*) "ascheq called with mesh" 112 write (*,*) "nstop", nstop, e 113 write (*,*) "error ",error 114 write (*,*) "y(1:3)" 115 write (*,*) y(1:3) 116 write (*,*) "y(mesh-2:mesh)" 117 write (*,*) y(mesh-2:mesh) 118 write (*,*) "yy(1:3)" 119 write (*,*) yy(1:3) 120 write (*,*) "yy(mesh-2:mesh)" 121 write (*,*) yy(mesh-2:mesh) 122 write (*,*) grid%sqr(1:3) 123 write (*,*) "sqrt(charge)", sqrt(charge) 124#endif 125 if (error > 1.d-8) then 126 call errore('c6_tfvw','auxiliary funtions veff(r) and y(r) are inaccurate',1) 127 end if 128! 129! initialize external perturbation (electric field) 130! 131 call init_dpot(grid%r,mesh,dvpot) 132! 133! derivative of xc-potential 134! 135 call dvxc_dn(mesh, rho, dvxc) 136! 137!write(*,'(1PE20.12)')sum(abs(dvxc)) 138!stop 139 write(6,'(5x,''Frequency dependent polarizability is written into freq-pol.dat'',/)') 140 141 c6 = 0.0d0 142 alpha = 0.0d0 143 144 open(1, file = 'freq-pol.dat') 145 write (1,'(15x," u alpha(angstrom) alpha(a.u.) ",/)') 146 ! 147 Nu = 230 148 Nc = 50 149 du1 = 0.1d0 150 du2 = 0.25d0 151 u = -du1 152 ! 153 do iu=0,Nu 154 ! 155 if (iu .le. 50) then 156 u = u + du1 157 else 158 u = u + du2 159 endif 160 ! 161 if (iu.eq.0) then 162 do i=1,mesh 163 dvscf(i) = dvpot(i) 164 drho_old(i) = 0.d0 165 end do 166 end if 167 beta = 0.05 168 dalpha = 1.0d+99 169 alpha = 0.d0 170 counter = 0 171 do while (dalpha > 1.d-9) 172 counter = counter + 1 173 ! 174 ! solve Sternheimer equation for the auxiliary wavefunction 175 ! 176 l = 1 177 ly = 0 178 call sternheimer(u*vw_lambda,l,ly,mesh,grid%dx,grid%r,grid%sqr,grid%r2,veff,zed,y,dvscf,dy) 179 dy = dy*vw_lambda 180 ! compute drho of r 181 ! 182 call drho_of_r(mesh, grid%dx, grid%r, grid%r2, y, dy, drho) 183#ifdef DEGUG 184 write (*,*) "========================" 185 write (*,*) "drho(1:3)" 186 write (*,*) drho(1:3) 187 write (*,*) "drho(20:22)" 188 write (*,*) drho(20:22) 189 write (*,*) "drho(40:42)" 190 write (*,*) drho(40:42) 191 write (*,*) "drho(mesh-2:mesh)" 192 write (*,*) drho(mesh-2:mesh) 193#endif 194 ! 195 ! compute dv of drho (including the TF term) 196 ! 197 l_add_tf_term = .true. 198 call dv_of_drho(mesh, grid%dx, grid%r,grid%r2,rho,drho,dvhx,dvxc,pp, l_add_tf_term) 199 200#ifdef DEGUG 201 write (*,*) "========================" 202 write (*,*) "dvhx(1:3)" 203 write (*,*) dvhx(1:3) 204 write (*,*) "dvhx(20:22)" 205 write (*,*) dvhx(20:22) 206 write (*,*) "dvhx(40:42)" 207 write (*,*) dvhx(40:42) 208 write (*,*) "dvhx(mesh-2:mesh)" 209 write (*,*) dvhx(mesh-2:mesh) 210 write (*,*) "========================" 211 write (*,*) "pp(1:3)" 212 write (*,*) pp(1:3) 213 write (*,*) "pp(20:22)" 214 write (*,*) pp(20:22) 215 write (*,*) "pp(40:42)" 216 write (*,*) pp(40:42) 217 write (*,*) "pp(mesh-2:mesh)" 218 write (*,*) pp(mesh-2:mesh) 219#endif 220 ! 221 ! mix 222 ! 223 error = 0.d0 224 error2 = 0.d0 225 do i=1,mesh 226 dvscf(i) = dvscf(i) + beta * (dvpot(i)+dvhx(i) -dvscf(i)) 227 error = error + abs (drho(i) -drho_old(i)) 228 error2 = error2 + abs (drho(i) -drho_old(i))* grid%r(i) * grid%dx 229 drho_old(i) = drho(i) 230 end do 231 dalpha = abs(alpha + pp(mesh)) 232 alpha = -pp(mesh) 233! write (*,'(4e16.6)') alpha, dalpha, error, error2 234 end do 235 236 write (1,'(17x, f8.4, 3x, 1PE14.6, 9x, 1PE14.6)') u, pp(mesh)*BOHR_RADIUS_ANGS**3, pp(mesh) 237 if (iu .eq. 0) & 238 write (6,'(5x, "Static polarizability: ", f10.5, " (in A^3) --->", f10.5,& 239 & " (in e^2a0^3)")') pp(mesh)*BOHR_RADIUS_ANGS**3, pp(mesh) 240 241 if (iu .eq. 0) factor = 0.5d0 * du1 242 if (iu .gt. 0 .and. iu .lt. Nc) factor = du1 243 if (iu .eq. Nc) factor = 0.5d0 * ( du1 + du2) 244 if (iu .gt. Nc .and. iu .lt. Nu) factor = du2 245 if (iu .eq. Nu) factor = 0.5d0 * du2 246 c6 = c6 + factor*alpha*alpha 247 248 end do 249 250 c6 = c6 * 3.d0 / pi 251 252 write (*,'(/, 5x, a, f12.6)') "C6 coefficient in units [e2*a0**5]", c6/e2 253 ! 254 write(6,'(/,5x,20(''-''),'' End of C6 calculation '',20(''-''),/)') 255 256 deallocate ( dy ) 257 deallocate ( veff, y, yy ) 258 deallocate ( dvpot, dvscf, drho, dvhx, pp ) 259 260 mesh = mesh_save 261 return 262end subroutine 263 264!-------------------------------------------------------------------- 265subroutine veff_of_rho(mesh,dx,r,r2,rho,y,veff) 266 !-------------------------------------------------------------------- 267 ! compute unperturbed auxiliary wavefunction y and 268 ! the corresponding effective potential veff 269 ! 270 use constants, only : e2, pi, fpi 271 implicit none 272 ! 273 ! I/O variables 274 ! 275 integer mesh 276 real (kind=8) :: dx, r(mesh), r2(mesh), rho(mesh), veff(mesh), y(mesh) 277 ! 278 ! local variables 279 ! 280 real (kind=8), allocatable :: vold(:) 281 real (kind=8) :: dx2, error 282 integer i, iter, k 283 284! 285! compute auxiliary wavefunction y 286! 287 do i=1,mesh 288 y(i) = sqrt(rho(i)*r(i)*fpi) 289 end do 290! 291! compute effective potential veff 292! 293 allocate (vold(mesh)) 294 295 do i=1,mesh 296 vold(i) = 0.d0 297 end do 298 dx2= dx*dx 299 error = 1.d0 300 k = 0 301 do while (error > 1.d-9) 302 ! 303 k=k+1 304 ! 305 do i=2,mesh-1 306 veff(i) = ( y(i+1)/y(i) + y(i-1)/y(i) -2.d0 )/dx2 & 307 - ( vold(i+1)*y(i+1)/y(i) + vold(i-1)*y(i-1)/y(i) -2.d0*vold(i) )/12.d0 308 end do 309 veff(1) = veff(2) + (veff(3)-veff(2))*(r(1)-r(2))/(r(3)-r(2)) 310 veff(mesh) = (y(mesh-1)/y(mesh) -2.d0 )/dx2 & 311 - (vold(mesh-1)*y(mesh-1)/y(mesh) -2.d0*vold(mesh) )/12.d0 312! 313! the routine that integrates the Sh.Eq. requires that v(mesh) is an upper bound 314! 315 veff(mesh) = max(veff(mesh),veff(mesh-1)) 316! 317 error = 0.d0 318 do i=1,mesh 319 error = error + abs( veff(i) - vold(i) ) 320 vold(i) = veff(i) 321 end do 322 error = error / mesh 323! write (*,*) 'iteration # ', k, error 324 end do 325 326 deallocate (vold) 327 ! 328 do i=1,mesh 329 veff(i) = (veff(i) -0.25d0)/r2(i) 330 end do 331 ! 332! open (7,file='veff.out',status='unknown',form='formatted') 333! write (7,*) "# r(i), rho(i), y(i), veff(i), veff(i)*r(i) " 334! do i=0,mesh 335! write (7,'(P6E15.6)') r(i), rho(i), y(i), veff(i), veff(i)*r(i) 336! end do 337! close (7) 338 ! 339 return 340end subroutine 341! 342!-------------------------------------------------------------------- 343subroutine dv_of_drho(mesh,dx,r,r2,rho,drho,dvhx,dvxc,pp, l_add_tf_term) 344 !-------------------------------------------------------------------- 345 use constants, only : e2, pi, fpi 346! use flags, only : HartreeFock, rpa 347 implicit none 348 ! 349 ! I/O variables 350 ! 351 logical l_add_tf_term 352 integer mesh 353 real (kind=8) :: dx, r(mesh), r2(mesh) 354 real (kind=8) :: rho(mesh), drho(mesh), dvhx(mesh), pp(mesh), dvxc(mesh) 355 ! 356 ! local variables 357 ! 358 real (kind=8) :: dr3, kf2, charge 359 real (kind=8), allocatable :: qq(:) 360 integer i 361 362 allocate (qq(mesh)) 363 364 do i=1,mesh 365 dr3 = fpi * r2(i) * r(i) * dx 366 pp(i) = drho(i) * r(i) * dr3 /3.d0 367 qq(i) = drho(i) / r2(i) * dr3 /3.d0 368 end do 369 do i=2,mesh 370 pp(i) = pp(i) + pp(i-1) 371 end do 372! write (*,*) "pp in dv_of_drho" 373! write (*,*) pp(1:6) 374! write (*,*) pp(mesh-5:mesh) 375! write (*,*) "qq -prima in dv_of_drho" 376! write (*,*) qq(1:6) 377! write (*,*) qq(mesh-5:mesh) 378 do i=mesh-1,1,-1 379 qq(i) = qq(i) + qq(i+1) 380 end do 381! write (*,*) "qq -dopo in dv_of_drho" 382! write (*,*) qq(1:6) 383! write (*,*) "r2 in dv_of_drho" 384! write (*,*) r2(1:6) 385! write (*,*) "r in dv_of_drho" 386! write (*,*) r(1:6) 387 do i=1,mesh 388 dvhx(i) = e2 * ( pp(i) / r2(i) + qq(i) * r(i) ) ! Hartree term 389 end do 390 391! write (*,*) "Hartree in dv_of_drho" 392! write (*,*) dvhx(1:6) 393! write (*,*) dvhx(mesh-5:mesh) 394! add TF term 395 if (l_add_tf_term) then 396 do i=1,mesh 397 kf2 = ( 3.d0*pi*pi*rho(i) )**(2.d0/3.d0) 398 dvhx(i) = dvhx(i) + e2/3.d0* kf2 / rho(i) * drho(i) 399 end do 400 end if 401! 402!add xc term 403! 404! write (*,*) "dvxc in dv_of_drho" 405! write (*,*) dvxc(1:6) 406 do i=1,mesh 407 dvhx(i) = dvhx(i) + dvxc(i) * drho(i) 408 end do 409! write (*,*) "Hartree + dvxc in dv_of_drho" 410! write (*,*) dvhx(1:3) 411! write (*,*) dvhx(mesh-2:mesh) 412! 413 deallocate (qq) 414 415 return 416end subroutine 417!---------------------------------------------------------------------- 418subroutine dvxc_dn(mesh, rho, dvxc) 419 !------------------------------------------------------------------- 420 ! compute the derivative of xc-potential w.r.t local density. 421 ! some routine in PH and flibs will be called 422 ! 423 use funct, only : dft_is_gradient 424 ! 425 implicit none 426 ! 427 ! I/O variables 428 ! 429 integer :: mesh 430 real(kind=8) :: rho(mesh), dvxc(mesh) 431 ! 432 ! local variables 433 ! 434 integer :: i 435 ! 436 ! 437 ! 438 if ( dft_is_gradient() ) & 439 call errore ('dvxc_dn', 'gradient correction to dvxc not yet implemented', 1) 440 ! 441 ! LDA only 442 ! 443 CALL dmxc_lda( mesh, rho, dvxc ) 444 ! 445 return 446 ! 447end subroutine dvxc_dn 448!-------------------------------------------------------------------- 449subroutine drho_of_r(mesh, dx, r, r2, y, dy, drho) 450 !-------------------------------------------------------------------- 451 ! compute the first order variation of the density from 452 ! the zeroth and first order auxiliary wavefunctions y and dy 453 ! 454 use constants, only : e2, pi, fpi 455 implicit none 456 ! 457 ! I/O vaiables 458 ! 459 integer mesh 460 real (kind=8) :: dx, r(mesh), r2(mesh), y(mesh), drho(mesh) 461 complex (kind=8) :: dy(mesh) 462 ! local variables 463 integer i 464 do i=1,mesh 465 drho(i) = 2.d0 * y(i) * real(dy(i)) * r(i) / (fpi*r2(i)) 466 end do 467 468 return 469end subroutine 470!-------------------------------------------------------------------- 471subroutine init_dpot(r,mesh,dvpot) 472 !-------------------------------------------------------------------- 473 ! 474 ! initialize external potential 475 ! 476 use constants, only : e2 477 implicit none 478 ! 479 ! I/O variables 480 ! 481 integer mesh 482 real (kind=8) :: r(mesh), dvpot(mesh) 483 ! 484 ! local variables 485 ! 486 integer i 487 488 do i =1,mesh 489 dvpot (i) = - e2*r(i) 490 end do 491 492 return 493end subroutine 494 495!-------------------------------------------------------------------- 496subroutine sternheimer(u, l, ll, mesh, dx, r, sqr, r2, vpot, zed, y, dvpot, dy) 497 !-------------------------------------------------------------------- 498 ! 499 ! solve the sternheimer equation for imaginary frequency 500 ! in radial coordinates by Numerov method 501 ! 502 implicit none 503 ! 504 ! I/O variables 505 ! 506 integer mesh, l, ll 507 real (kind=8) :: u, dx, zed 508 real (kind=8) :: r(mesh), sqr(mesh), r2(mesh) 509 real (kind=8) :: vpot(mesh), y(mesh), dvpot(mesh) 510 complex (kind=8) :: dy(mesh) 511 ! 512 ! local variables 513 ! 514 integer i, icl, j 515 real (kind=8) :: ddx12, sqlhf, x2l2 516 complex (kind=8) :: gg, aa, bb, fac, e 517 complex (kind=8), allocatable :: f(:), g(:), yy(:) 518 519 allocate ( f(mesh), g(mesh), yy(mesh) ) 520 521 ddx12=dx*dx/12.d0 522 sqlhf = (l+0.5d0)**2 523 x2l2 = 2*l+2 524 e = cmplx (0.d0, u, kind=8) 525 ! 526 ! set up the f-function and determine the position of its last 527 ! change of sign 528 ! f < 0 (approximatively) means classically allowed region 529 ! f > 0 " " " forbidden " 530 ! 531 icl = 2 532! f(0) = ddx12 *( sqlhf + r2(0) * (vpot(0)-e) ) 533 f(1) = ddx12 *( r2(1) * (vpot(1)-e) ) 534 do i=2,mesh 535! f(i) = ddx12 * ( sqlhf + r2(i) *(vpot(i)-e) ) 536 f(i) = ddx12 * ( r2(i) *(vpot(i)-e) ) 537 if( real(f(i)) .ne. sign(real(f(i)),real(f(i-1))) & 538 .and. real(f(i)).gt.0d0 ) icl=i 539 end do 540 ! write (*,*) icl 541 542 do i=1,mesh 543 f(i) = ddx12 * ( sqlhf + r2(i) *(vpot(i)-e) ) 544 end do 545 546 do i=1,mesh 547 f(i)=1.0d0-f(i) 548 g(i)= ddx12 * r2(i) * dvpot(i) * y(i) 549 end do 550 ! step 1) determine a solution to the homogeneous equation that 551 ! satisfies proper boundary conditions at 0 and infty 552 ! NB it will NOT satisfy the homogeneous equation at icl 553 ! 554 ! determination of the wave-function in the first two points 555 ! 556 yy(1) = r(1)**(l+1) *(1.d0 - 2.d0*zed*r(1)/x2l2) / sqr(1) 557 yy(2) = r(2)**(l+1) *(1.d0 - 2.d0*zed*r(2)/x2l2) / sqr(2) 558 ! 559 ! outward integration 560 ! 561 do i =2, icl-1 562 yy(i+1)=( (12.d0-10.d0*f(i))*yy(i)-f(i-1)*yy(i-1) )/f(i+1) 563 end do 564 ! 565 ! rescale to 1 at icl 566 ! 567 fac = 1.d0/yy(icl) 568 do i =1, icl 569 yy(i)= yy(i) * fac 570 end do 571 ! 572 ! determination of the wave-function in the last two points 573 ! assuming y(mesh+1) = 0 and y(mesh) = dx 574 ! 575 yy(mesh) = dx 576 yy(mesh-1) = (12.d0-10.d0*f(mesh))*yy(mesh)/f(mesh-1) 577 ! 578 ! inward integration 579 ! 580 do i = mesh-1,icl+1,-1 581 yy(i-1)=( (12.d0-10.d0*f(i))*yy(i)-f(i+1)*yy(i+1) )/f(i-1) 582 if (abs(yy(i-1)).gt.1.d6) then 583 fac = 1.d0/yy(i-1) 584 do j=i-1,mesh 585 yy(j) = yy(j) * fac 586 end do 587 end if 588 end do 589 ! 590 ! rescale to 1 at icl 591 ! 592 fac = 1.d0 / yy(icl) 593 do i = icl, mesh 594 yy(i)= yy(i) * fac 595 end do 596 ! 597 ! step 2) determine a solution to the inhomogeneous equation that 598 ! satisfies proper boundary conditions at 0 and infty 599 ! and matches in value at icl 600 ! NB it will NOT satisfy the inhomogeneous equation at icl 601 ! 602 ! determination of the wave-function in the first two points 603 ! 604 dy(1) = r(1)**(l+1) *(1.d0 - 2.d0*zed*r(1)/x2l2) / sqr(1) 605 dy(2) = r(2)**(l+1) *(1.d0 - 2.d0*zed*r(2)/x2l2) / sqr(2) 606 ! 607 ! outward integration 608 ! 609 do i =2, icl-1 610 gg = g(i+1) + 10.d0*g(i) + g(i-1) 611 dy(i+1)=( (12.d0-10.d0*f(i))*dy(i)-f(i-1)*dy(i-1) + gg )/f(i+1) 612 end do 613 ! 614 ! choose the solution that goes to 1 in icl 615 ! 616 fac = dy(icl) - 1.d0 617 do i = 1,icl 618 dy(i) = dy(i) - fac * yy(i) 619 end do 620 ! 621 ! determination of the wave-function in the last two points 622 ! assuming y(mesh+1) = 0 and y(mesh) = dx 623 ! 624 dy(mesh) = dx 625 dy(mesh-1) = (12.d0-10.d0*f(mesh))*yy(mesh)/f(mesh-1) 626 ! 627 ! inward integration 628 ! 629 do i = mesh-1,icl+1,-1 630 gg = g(i+1) + 10.d0*g(i) + g(i-1) 631 dy(i-1)=( (12.d0-10.d0*f(i))*dy(i)-f(i+1)*dy(i+1) + gg )/f(i-1) 632 if (abs(dy(i-1)).gt.1.d6) then 633 if (abs(yy(i-1)).gt. 1.d-12) then 634 fac = ( dy(i-1) - 1.d0 ) / yy(i-1) 635 do j = i-1, mesh 636 dy(j) = dy(j) - fac * yy(j) 637 end do 638 else 639 do j = i-1, mesh 640 dy(j) = 0.d0 641 end do 642 end if 643 end if 644 end do 645 ! 646 ! choose the solution that goes to 1 in icl 647 ! 648 fac = dy(icl) - 1.d0 649 do i = icl, mesh 650 dy(i) = dy(i) - fac * yy(i) 651 end do 652 ! 653 ! step 3) chose the proper combination of dy and yy so that the 654 ! inhomogeneous equation is satisfied also in icl 655 i = icl 656 gg = g(i+1) + 10.d0*g(i) + g(i-1) 657 aa= (12.d0-10.d0*f(i))*dy(i)-f(i+1)*dy(i+1)-f(i-1)*dy(i-1) + gg 658 bb= (12.d0-10.d0*f(i))*yy(i)-f(i+1)*yy(i+1)-f(i-1)*yy(i-1) 659 ! write (*,*) aa, bb 660 fac = aa/bb 661 do i=1,mesh 662 dy(i) = dy(i) - fac * yy(i) 663 end do 664 deallocate ( f, g, yy ) 665 666 return 667 668end subroutine 669 670