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