1!
2! Copyright (C) 2004-2010 Quantum ESPRESSO 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!
8SUBROUTINE lschps (mode, z, eps, grid, nin, n, l, e, v, u, nstop)
9  !
10  ! integrates radial pauli-type scalar-relativistic equation
11  ! on a logarithmic grid
12  ! modified routine to be used in finding norm-conserving
13  ! pseudopotential
14  !
15  ! on input:
16  !   mode = 1 find energy and wavefunction of bound states,
17  !            scalar-relativistic (all-electron)
18  !   mode = 2 find energy and wavefunction of bound state,
19  !            nonrelativistic (pseudopotentials)
20  !   mode = 3 fixed-energy calculation, for logarithmic derivatives
21  !   mode = 4 find energy which produces a specified logarithmic
22  !            derivative (nonrelativistic, pseudopotentials)
23  !   mode = 5 is for pseudopotential to produce wavefunction beyond
24  !            radius used for pseudopotential construction
25  !   z    = atomic number
26  !   eps  = convergence factor: eiganvalue is considered converged if
27  !          the correction to eigenvalue is smaller in magnitude than
28  !          eps times the magnitude of the current guess
29  !   grid = structure containing radial grid information
30  !   l, n = main and angular quantum numbers
31  !   e    = starting estimate of the energy (mode=1,2)
32  !          fixed energy at which the wavefctn is calculated (mode=3,4)
33  !   v(i) = self-consistent potential
34  !   nin  = integration up to r(nin) (mode=3,4,5)
35  !
36  ! on output:
37  !   e    = final energy (mode=1,2)
38  !   u(i) = radial wavefunction (defined as the radial part of the wavefct
39  !          multiplied by r)
40  !   nstop= 0 if regular termination, 1 otherwise
41  !   nin  = last grid point for which the wavefct is calculated (mode=1,2)
42  !
43  USE kinds, ONLY : DP
44  USE radial_grids, ONLY: radial_grid_type
45  USE ld1inc, ONLY : cau_fact
46  IMPLICIT NONE
47  !
48  ! I/O variables
49  !
50  INTEGER, INTENT (in) :: mode, n, l
51  real(DP), INTENT(in) :: z, eps
52  TYPE (radial_grid_type), INTENT(in) :: grid
53  real(DP), INTENT(in) :: v(grid%mesh)
54  INTEGER, INTENT(inout) :: nin
55  real(DP), INTENT(inout) :: e
56  INTEGER, INTENT(out) :: nstop
57  real (DP), INTENT(out) :: u(grid%mesh)
58  !
59  ! local variables
60  !
61  INTEGER, PARAMETER :: maxter=60
62  real(DP), EXTERNAL:: aei, aeo, aii, aio
63  ! arrays  used as work space
64  real(DP),ALLOCATABLE :: up(:),upp(:),cf(:),dv(:),fr(:),frp(:)
65  real(DP):: al, als, cn
66  real(DP):: de, emax, emin
67  real(DP):: fss, gamma, ro, sc
68  real(DP):: sls, sn, uld, uout,  upin, upout
69  real(DP):: xkap
70  INTEGER:: i, it, mmax, n_it, node, mch, ierr
71  !
72  !
73  nstop=0
74  al   = grid%dx
75  mmax = grid%mesh
76
77  ALLOCATE(up(mmax), stat=ierr)
78  ALLOCATE(upp(mmax), stat=ierr)
79  ALLOCATE(cf(mmax), stat=ierr)
80  ALLOCATE(dv(mmax), stat=ierr)
81  ALLOCATE(fr(mmax), stat=ierr)
82  ALLOCATE(frp(mmax), stat=ierr)
83
84  uld=0.0_dp
85  !
86  !
87  IF(mode == 1 .or. mode == 3) THEN
88     !     relativistic calculation
89     !     fss=(1.0_dp/137.036_dp)**2
90     fss=(1.0_dp/cau_fact)**2
91     IF(l == 0) THEN
92        gamma=sqrt(1.0_dp-fss*z**2)
93     ELSE
94        gamma=(l*sqrt(l**2-fss*z**2) + &
95             (l+1)*sqrt((l+1)**2-fss*z**2))/(2*l+1)
96     ENDIF
97  ELSE
98     !     non-relativistic calculation
99     fss=1.0e-20_dp
100     gamma=l+1
101  ENDIF
102  !
103  sls=l*(l+1)
104  !
105  ! emin, emax = estimated bounds for e
106  !
107  IF(mode == 1 .or. mode == 2) THEN
108     emax=v(mmax)+sls/grid%r(mmax)**2
109     emin=0.0_dp
110     DO i=1,mmax
111        emin=min(emin,v(i)+sls/grid%r(i)**2)
112     ENDDO
113     IF(e > emax) e=1.25_dp*emax
114     IF(e < emin) e=0.75_dp*emin
115     IF(e > emax) e=0.5_dp*(emax+emin)
116  ELSEIF(mode == 4) THEN
117     emax=e + 10.0_dp
118     emin=e - 10.0_dp
119  ENDIF
120  !
121  DO i=1,4
122     u(i)=0.0_dp
123     up(i)=0.0_dp
124     upp(i)=0.0_dp
125  ENDDO
126  als=al**2
127  !
128  ! calculate dv/dr for darwin correction
129  !
130  CALL derv (mmax, al, grid%r, v, dv )
131  !
132  !     starting of loop on energy for bound state
133  !
134  DO n_it = 1, maxter
135     !
136     ! coefficient array for u in differential eq.
137     DO i=1,mmax
138        cf(i)=als*(sls + (v(i)-e)*grid%r(i)**2)
139     ENDDO
140     !
141     ! find classical turning point for matching
142     !
143     IF(mode == 1 .or. mode == 2) THEN
144        DO i=mmax,2,-1
145           IF(cf(i-1) <= 0.0_dp .and. cf(i) > 0.0_dp) THEN
146              mch=i
147              GOTO 40
148           ENDIF
149        ENDDO
150        !PRINT '('' warning: wfc '',2i2,'' no turning point'')', n, l
151        e=0.0_dp
152        DO i=1,mmax
153           u (i)=0.0_dp
154        ENDDO
155        nstop=1
156        GOTO 999
157     ELSE
158        mch=nin
159     ENDIF
16040   CONTINUE
161
162     !  relativistic coefficient arrays for u (fr) and up (frp).
163     DO i=1,mmax
164        fr(i)=als*(grid%r(i)**2)*0.25_dp*(-fss*(v(i)-e)**2 + &
165             fss*dv(i)/ (grid%r(i)*(1.0_dp+0.25_dp*fss*(e-v(i)))))
166        frp(i)=-al*grid%r(i)*0.25_dp*fss*dv(i)/(1.0_dp+0.25_dp*fss*(e-v(i)))
167     ENDDO
168     !
169     ! start wavefunction with series
170     !
171     DO i=1,4
172        u(i)=grid%r(i)**gamma
173        up(i)=al*gamma*grid%r(i)**gamma
174        upp(i)=(al+frp(i))*up(i)+(cf(i)+fr(i))*u(i)
175     ENDDO
176     !
177     ! outward integration using predictor once, corrector
178     ! twice
179     node=0
180     !
181     DO i=4,mch-1
182        u(i+1)=u(i)+aeo(up,i)
183        up(i+1)=up(i)+aeo(upp,i)
184        DO it=1,2
185           upp(i+1)=(al+frp(i+1))*up(i+1)+(cf(i+1)+fr(i+1))*u(i+1)
186           up(i+1)=up(i)+aio(upp,i)
187           u(i+1)=u(i)+aio(up,i)
188        ENDDO
189        IF(u(i+1)*u(i) <= 0.0_dp) node=node+1
190     ENDDO
191     !
192     uout=u(mch)
193     upout=up(mch)
194     !
195     IF(node-n+l+1 == 0 .or. mode == 3 .or. mode == 5) THEN
196        !
197        IF(mode == 1 .or. mode == 2) THEN
198           !
199           ! start inward integration at 10*classical turning
200           ! point with simple exponential
201           nin=mch+2.3_dp/al
202           IF(nin+4 > mmax) nin=mmax-4
203           xkap=sqrt(sls/grid%r(nin)**2 + 2.0_dp*(v(nin)-e))
204           !
205           DO i=nin,nin+4
206              u(i)=exp(-xkap*(grid%r(i)-grid%r(nin)))
207              up(i)=-grid%r(i)*al*xkap*u(i)
208              upp(i)=(al+frp(i))*up(i)+(cf(i)+fr(i))*u(i)
209           ENDDO
210           !
211           ! integrate inward
212           !
213           DO i=nin,mch+1,-1
214              u(i-1)=u(i)+aei(up,i)
215              up(i-1)=up(i)+aei(upp,i)
216              DO it=1,2
217                 upp(i-1)=(al+frp(i-1))*up(i-1)+(cf(i-1)+fr(i-1))*u(i-1)
218                 up(i-1)=up(i)+aii(upp,i)
219                 u(i-1)=u(i)+aii(up,i)
220              ENDDO
221           ENDDO
222           !
223           ! scale outside wf for continuity
224           sc=uout/u(mch)
225           !
226           DO i=mch,nin
227              up(i)=sc*up(i)
228              u (i)=sc*u (i)
229           ENDDO
230           !
231           upin=up(mch)
232           !
233        ELSE
234           !
235           upin=uld*uout
236           !
237        ENDIF
238        !
239        ! perform normalization sum
240        !
241        ro=grid%r(1)*exp(-0.5_dp*grid%dx)
242        sn=ro**(2.0_dp*gamma+1.0_dp)/(2.0_dp*gamma+1.0_dp)
243        !
244        DO i=1,nin-3
245           sn=sn+al*grid%r(i)*u(i)**2
246        ENDDO
247        !
248        sn=sn + al*(23.0_dp*grid%r(nin-2)*u(nin-2)**2 &
249             + 28.0_dp*grid%r(nin-1)*u(nin-1)**2 &
250             +  9.0_dp*grid%r(nin  )*u(nin  )**2)/24.0_dp
251        !
252        ! normalize u
253        cn=1.0_dp/sqrt(sn)
254        uout=cn*uout
255        upout=cn*upout
256        upin=cn*upin
257        !
258        DO i=1,nin
259           up(i)=cn*up(i)
260           u(i)=cn*u(i)
261        ENDDO
262        DO i=nin+1,mmax
263           u(i)=0.0_dp
264        ENDDO
265        !
266        ! exit for fixed-energy calculation
267        !
268        IF(mode == 3 .or. mode == 5) GOTO 999
269
270        ! perturbation theory for energy shift
271        de=uout*(upout-upin)/(al*grid%r(mch))
272        !
273        ! convergence test and possible exit
274        !
275        IF ( abs(de) < max(abs(e),0.2_dp)*eps) GOTO 999
276        !
277        IF(de > 0.0_dp) THEN
278           emin=e
279        ELSE
280           emax=e
281        ENDIF
282        e=e+de
283        IF(e > emax .or. e < emin) e=0.5_dp*(emax+emin)
284        !
285     ELSEIF(node-n+l+1 < 0) THEN
286        ! too few nodes
287        emin=e
288        e=0.5_dp*(emin+emax)
289
290     ELSE
291        ! too many nodes
292        emax=e
293        e=0.5_dp*(emin+emax)
294     ENDIF
295  ENDDO
296
297  !PRINT '('' warning: wfc '',2i2,'' not converged'')', n, l
298  u=0.0_dp
299  nstop=1
300  !
301  ! deallocate arrays and exit
302  !
303999 CONTINUE
304  DEALLOCATE(frp)
305  DEALLOCATE(fr)
306  DEALLOCATE(dv)
307  DEALLOCATE(cf)
308  DEALLOCATE(upp)
309  DEALLOCATE(up)
310  RETURN
311
312END SUBROUTINE lschps
313!
314!----------------------------------------------------------------
315SUBROUTINE lschps_meta (mode, z, eps, grid, nin, n, l, e, v, vtau, &
316                        u, nstop)
317  !----------------------------------------------------------------
318  !
319  ! Meta-GGA version of lschps
320  ! vtau is the meta-GGA potential
321  !
322  USE kinds, ONLY : DP
323  USE radial_grids, ONLY: radial_grid_type
324  USE ld1inc, ONLY : cau_fact
325  IMPLICIT NONE
326  !
327  ! I/O variables
328  !
329  INTEGER, INTENT (in) :: mode, n, l
330  real(DP), INTENT(in) :: z, eps
331  TYPE (radial_grid_type), INTENT(in) :: grid
332  real(DP), INTENT(in) :: v(grid%mesh), vtau(grid%mesh)
333  INTEGER, INTENT(inout) :: nin
334  real(DP), INTENT(inout) :: e
335  INTEGER, INTENT(out) :: nstop
336  real (DP), INTENT(out) :: u(grid%mesh)
337  !
338  ! local variables
339  !
340  INTEGER, PARAMETER :: maxter=60
341  real(DP), EXTERNAL:: aei, aeo, aii, aio, d2u, int_0_inf_dr
342  ! arrays used as work space
343  real(DP), ALLOCATABLE :: up(:),upp(:),cf(:),dv(:),fr(:),frp(:),dvtau(:)
344  real(DP):: al, als, cn
345  real(DP):: de, emax, emin
346  real(DP):: fss, gamma, ro, sc
347  real(DP):: sls, sn, uld, uout,  upin, upout
348  real(DP):: xkap, mm, mvt, fac1, fact2, alpha
349  INTEGER :: i, it, mmax, n_it, node, mch, ierr
350  LOGICAL ::  rel, find_e, fixed_e
351  !
352  !
353  IF (mode<1.or.mode>5) STOP 'lschps: wrong mode'
354  !
355  nstop=0
356  al   = grid%dx
357  mmax = grid%mesh
358  !
359  ALLOCATE(up(mmax), stat=ierr)
360  ALLOCATE(upp(mmax), stat=ierr)
361  ALLOCATE(cf(mmax), stat=ierr)
362  ALLOCATE(dv(mmax), stat=ierr)
363  ALLOCATE(fr(mmax), stat=ierr)
364  ALLOCATE(frp(mmax), stat=ierr)
365  ALLOCATE(dvtau(mmax), stat=ierr)
366  !
367  uld=0.0_dp
368  !
369  rel    = ( mode == 1 .or. mode == 3 )
370  find_e = ( mode == 1 .or. mode == 2 )
371  fixed_e= ( mode == 3 .or. mode == 5 )
372  !
373  ! fss  : square of the fine structure constant
374  ! gamma: u(r)=r^gamma is the asymptotic behavior of u(r) for r->0
375  !
376  IF (rel) THEN
377     ! relativistic calculation
378     !    fss=(1.0_dp/137.036_dp)**2
379     fss=(1.0_dp/cau_fact)**2
380     IF(l == 0) THEN
381        gamma=sqrt(1.0_dp-fss*z**2)
382     ELSE
383        gamma=(l*sqrt(l**2-fss*z**2) + &
384              (l+1)*sqrt((l+1)**2-fss*z**2))/(2*l+1)
385     ENDIF
386  ELSE
387     ! nonrelativistic calculation
388     fss=0.0_dp
389     gamma=l+1
390  ENDIF
391  sls=l*(l+1)
392  !
393  ! emin, emax = estimated bounds for e
394  !
395  IF(find_e) THEN
396     emax=v(mmax)+sls/grid%r(mmax)**2
397     emin=0.0_dp
398     DO i=1,mmax
399        emin=min(emin,v(i)+sls/grid%r(i)**2)
400     ENDDO
401     IF(e > emax) e=1.25_dp*emax
402     IF(e < emin) e=0.75_dp*emin
403     IF(e > emax) e=0.5_dp*(emax+emin)
404  ELSEIF(mode == 4) THEN
405     emax=e + 10.0_dp
406     emin=e - 10.0_dp
407  ENDIF
408  !
409  DO i=1,mmax
410     u(i)=0.0_dp
411     up(i)=0.0_dp
412     upp(i)=0.0_dp
413  ENDDO
414  als=al**2
415  !
416  ! calculate dv = dv/dr for darwin correction
417  ! calculate dvtau = d tau / dr
418  !
419  CALL derV(mmax,al,grid%r,v,dv)
420  CALL derV(mmax,al,grid%r,vtau,dvtau)
421  !
422  !     starting of loop on energy for bound state
423  !
424  DO n_it = 1, maxter
425     !
426     !     find classical turning point for matching
427     !
428     DO i=1,mmax
429        cf(i)=(sls + (v(i)-e)*grid%r(i)**2)
430     ENDDO
431     IF(find_e) THEN
432        DO i=mmax,2,-1
433           IF(cf(i-1) <= 0.0_dp .and. cf(i) > 0.0_dp) THEN
434              mch=i
435              GOTO 40
436           ENDIF
437        ENDDO
438        !PRINT '('' warning: wfc '',2i2,'' no turning point'')', n, l
439        e=0.0
440        DO i=1,mmax
441           u (i)=0.0
442           ! up(i)=0.0
443        ENDDO
444        nstop=1
445        GOTO 999
446     ELSE
447        mch=nin
448     ENDIF
44940   CONTINUE
450     !
451     ! coefficient array for u in differential eq.
452     !
453     DO i=1,mmax
454        mm=1.0_dp+0.25_dp*fss*(e-v(i))
455        mvt=1.0_dp/(1.0_dp+0.5_dp*vtau(i))
456        fac1=grid%r(i)*fss*0.25_dp*dv(i)*mvt/mm
457        !
458        cf(i)=als*(sls+mvt*(mm*(v(i)-e)*grid%r(i)**2                     &
459             &           +grid%r(i)*0.5_dp*dvtau(i)))
460        !
461        !            cf(i)=als*(sls+mvt*(mm*(v(i)-e)*grid%r(i)**2))
462        !
463        !     relativistic coefficient arrays for u (fr) and up (frp).
464        !
465        fr(i)=als*fac1
466        frp(i)=-al*(0.5*mvt*grid%r(i)*dvtau(i)+fac1)
467        !            frp(i)=-al*fac1
468        !
469     ENDDO
470     !
471     !     start wavefunction with series
472     !
473     fac1=1.0_dp+0.5_dp*vtau(1)
474     alpha=-z/(l+1.0_dp)/fac1*grid%r(1)
475     fact2=-l/(2.0_dp*(l+1.0_dp))
476     DO i=1,4
477        u(i) = grid%r(i)**gamma*exp(alpha)*(1+0.5_dp*vtau(i))**fact2
478        up(i)=al*(gamma+(fact2-z/(l+1.0_dp))*grid%r(i)/fac1)*u(i)
479        upp(i)=d2u(al,cf(i),fr(i),frp(i),u(i),up(i))
480        fac1=(1.0_dp+0.5*vtau(i+1))
481        alpha = alpha - z/(l+1.0_dp)/fac1*(grid%r(i+1)-grid%r(i))
482     ENDDO
483     !
484     !     outward integration using predictor once, corrector twice
485     !
486     node=0
487     DO i=4,mch-1
488        u (i+1)=u (i)+aeo(up,i)
489        up(i+1)=up(i)+aeo(upp,i)
490        DO it=1,2
491           upp(i+1)=d2u(al,cf(i+1),fr(i+1),frp(i+1),u(i+1),up(i+1))
492           up (i+1)=up(i)+aio(upp,i)
493           u  (i+1)=u (i)+aio(up,i)
494        ENDDO
495        IF(u(i+1)*u(i) <= 0.0_dp) node=node+1
496     ENDDO
497     !
498     uout=u(mch)
499     upout=up(mch)
500     !
501     IF(node-n+l+1 == 0 .or. fixed_e) THEN
502        !
503        IF (find_e) THEN
504           !
505           ! good number of nodes: start inward integration
506           ! at 10*classical turning point with simple exponential
507           !
508           nin=mch+2.3_dp/al
509           IF(nin+4 > mmax ) nin=mmax-4
510           !
511           mm=1.0_dp+0.25_dp*fss*(e-v(nin))
512           fac1=0.5_dp*dvtau(nin)/(1.0_dp+mm*0.5_dp*vtau(nin))
513           fact2=0.5_dp*dvtau(nin)/grid%r(nin)
514           alpha=1.0_dp+mm*0.5_dp*vtau(nin)
515           xkap=sqrt(fac1**2.0_dp + 4.0_dp*(sls/grid%r(nin)**2 + &
516                   (v(nin)-e+fact2)/alpha))
517           !
518           xkap=(-fac1+xkap)/2.0_dp
519           DO i=nin,nin+4
520              u (i) =exp(-xkap*(grid%r(i)-grid%r(nin)))
521              up(i) =-grid%r(i)*al*xkap*u(i)
522              upp(i)=d2u(al,cf(i),fr(i),frp(i),u(i),up(i))
523           ENDDO
524           !
525           ! integrate inward
526           !
527           DO i=nin,mch+1,-1
528              u (i-1)=u (i)+aei(up,i)
529              up(i-1)=up(i)+aei(upp,i)
530              DO it=1,2
531                 upp(i-1)=d2u(al,cf(i-1),fr(i-1),frp(i-1),u(i-1),up(i-1))
532                 up(i-1)=up(i)+aii(upp,i)
533                 u (i-1)=u (i)+aii(up,i)
534              ENDDO
535           ENDDO
536           !
537           ! scale outside wf for continuity
538           !
539           sc=uout/u(mch)
540           DO i=mch,nin
541              up(i)=sc*up(i)
542              u (i)=sc*u (i)
543           ENDDO
544           !
545           upin=up(mch)
546           !
547        ELSE
548           !
549           ! this is used only in fixed-logarithmic-derivative calculation
550           !
551           upin=uld*uout
552           !
553        ENDIF
554        !
555        ! normalization: upp is used to store u^2
556        !
557        DO i=1,nin
558           upp(i)=u(i)**2
559        ENDDO
560        !
561        ! integral over dr (includes series expansion for r->0) ...
562        ! assumes u(r)=r^(l+1) behaviour at the origin
563        !
564        sn=int_0_inf_dr ( upp, grid, nin-3, 2*l+2 )
565        !
566        ! ...extrapolation for the asymptotic behaviour (maybe useless...)
567        !
568        sn=sn + al*(23.0_dp*grid%r(nin-2)*upp(nin-2)     &
569                  + 28.0_dp*grid%r(nin-1)*upp(nin-1)     &
570                  +  9.0_dp*grid%r(nin  )*upp(nin  ))/24.0_dp
571        !
572        ! ...normalize u
573        !
574        cn=1.0_dp/sqrt(sn)
575        uout=cn*uout
576        upout=cn*upout
577        upin=cn*upin
578        !
579        DO i=1,nin
580           up(i)=cn*up(i)
581           u(i)=cn*u(i)
582        ENDDO
583        DO i=nin+1,mmax
584           u(i) =0.0_dp
585           up(i)=0.0_dp
586        ENDDO
587        !
588        !
589        ! exit for fixed-energy calculation
590        !
591        IF (fixed_e)  GOTO 999
592        !
593        ! perturbation theory for energy shift
594        de=uout*(upout-upin)/(al*grid%r(mch))
595        !
596        ! convergence test and possible exit
597        !
598        IF ( abs(de) < max(abs(e),0.2_dp)*eps) GOTO 999
599        !
600        IF(de > 0.0_dp) THEN
601           emin=e
602        ELSE
603           emax=e
604        ENDIF
605        e=e+de
606        IF(e > emax .or. e < emin) e=0.5_dp*(emax+emin)
607        !
608     ELSEIF(node-n+l+1 < 0) THEN
609        !
610        ! too few nodes
611        !
612        emin=e
613        e=0.5_dp*(emin+emax)
614        !
615     ELSE
616        !
617        ! too many nodes
618        !
619        emax=e
620        e=0.5_dp*(emin+emax)
621        !
622     ENDIF
623     !
624     ! loop back to converge e
625     !
626  ENDDO
627  !
628  !PRINT '('' warning: wfc '',2i2,'' not converged'')', n, l
629  nstop=1
630  u=0.0_dp
631  nstop=1
632  !
633  ! deallocate arrays and exit
634  !
635999 CONTINUE
636
637  DEALLOCATE(frp)
638  DEALLOCATE(fr)
639  DEALLOCATE(dv)
640  DEALLOCATE(cf)
641  DEALLOCATE(upp)
642  DEALLOCATE(up)
643  !do i=1,mmax
644  !   up(i)=up(i)/grid%r(i)/al
645  !enddo
646  RETURN
647  !
648END SUBROUTINE lschps_meta
649!
650!---------------------------------------------------------------
651FUNCTION d2u(al,cf,fr,frp,u,up)
652  !---------------------------------------------------------------
653  ! second derivative from radial KS equation
654  USE kinds, ONLY : DP
655  IMPLICIT NONE
656  real(dp):: d2u
657  real(dp), INTENT(in):: al,cf,fr,frp,u,up
658  !
659  d2u = (al+frp)*up + (cf+fr)*u
660  RETURN
661END FUNCTION d2u
662
663!-----------------------------------------------------------
664FUNCTION estimatealpha(mmax,u,up,al,r)
665  !-----------------------------------------------------------
666  USE kinds, ONLY : DP
667  IMPLICIT NONE
668  INTEGER, INTENT(in) :: mmax
669  real(dp) :: estimatealpha
670  real(dp), INTENT(in) :: u(mmax),up(mmax),al,r(mmax)
671  INTEGER i,istart,iend
672  estimatealpha = 0.0_dp
673  istart=5
674  iend=100
675  DO i=istart,iend
676     IF(u(i) > 1.0d-8) THEN
677        estimatealpha=estimatealpha+(1.0_dp-up(i)/u(i)/al)/r(i)
678     ENDIF
679  ENDDO
680  estimatealpha=estimatealpha/(iend-istart+1)
681  RETURN
682END FUNCTION estimatealpha
683
684FUNCTION aei(y,j)
685  !
686  USE kinds, ONLY : DP
687  IMPLICIT NONE
688  INTEGER j
689  real(DP):: y(j+3), aei
690  !
691  aei=-(4.16666666667e-2_dp)*(55.0_dp*y(j)-59.0_dp*y(j+1) &
692       +37.0_dp*y(j+2)-9.0_dp*y(j+3))
693  RETURN
694END FUNCTION aei
695!
696! adams extrapolation and interpolation formulas for
697! outward and inward integration, abramowitz and stegun, p. 896
698FUNCTION aeo(y,j)
699  !
700  USE kinds, ONLY : DP
701  IMPLICIT NONE
702  INTEGER:: j
703  real(DP):: y(j), aeo
704  !
705  aeo=(4.16666666667e-2_dp)*(55.0_dp*y(j)-59.0_dp*y(j-1) &
706       +37.0_dp*y(j-2)-9.0_dp*y(j-3))
707  RETURN
708END FUNCTION aeo
709!
710FUNCTION aii(y,j)
711  !
712  USE kinds, ONLY : DP
713  IMPLICIT NONE
714  INTEGER:: j
715  real(DP) :: y(j+2), aii
716  !
717  aii=-(4.16666666667e-2_dp)*(9.0_dp*y(j-1)+19.0_dp*y(j) &
718       -5.0_dp*y(j+1)+y(j+2))
719  RETURN
720END FUNCTION aii
721!
722FUNCTION aio(y,j)
723  !
724  USE kinds, ONLY : DP
725  IMPLICIT NONE
726  INTEGER :: j
727  real(DP):: y(j+1), aio
728  !
729  aio=(4.16666666667e-2_dp)*(9.0_dp*y(j+1)+19.0_dp*y(j) &
730       -5.0_dp*y(j-1)+y(j-2))
731  RETURN
732END FUNCTION aio
733
734SUBROUTINE derV (mmax,al,r,v,dv)
735  ! dv = dv/dr
736  USE kinds, ONLY : dp
737  IMPLICIT NONE
738  INTEGER, INTENT(in)  :: mmax
739  REAL(dp), INTENT(in) :: al, r(mmax), v(mmax)
740  REAL(dp), INTENT(out):: dv(mmax)
741  !
742  INTEGER :: i
743  !
744  dv(1)=(-50.0_dp*v(1)+96.0_dp*v(2)-72.0_dp*v(3)+32.0_dp*v(4) &
745       -6.0_dp*v(5))/(24.0_dp*al*r(1))
746  dv(2)=(-6.0_dp*v(1)-20.0_dp*v(2)+36.0_dp*v(3)-12.0_dp*v(4) &
747       +2.0_dp*v(5))/(24.0_dp*al*r(2))
748  !
749  DO i=3,mmax-2
750     dv(i)=(2.0_dp*v(i-2)-16.0_dp*v(i-1)+16.0_dp*v(i+1) &
751          -2.0_dp*v(i+2))/(24.0_dp*al*r(i))
752  ENDDO
753  !
754  dv(mmax-1)=( 3.0_dp*v(mmax)+10.0_dp*v(mmax-1)-18.0_dp*v(mmax-2)+ &
755       6.0_dp*v(mmax-3)-v(mmax-4))/(12.0_dp*al*r(mmax-1))
756  dv(mmax)=( 25.0_dp*v(mmax)-48.0_dp*v(mmax-1)+36.0_dp*v(mmax-2)-&
757       16.0_dp*v(mmax-3)+3.0_dp*v(mmax-4))/(12.0_dp*al*r(mmax))
758  !
759  RETURN
760  !
761END SUBROUTINE derV
762