1!
2! Copyright (C) 2002-2005 FPMD-CPV groups
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!-----------------------------------------------------------------------------!
9!   This module is basad on a similar module from CP2K
10!-----------------------------------------------------------------------------!
11
12    MODULE splines
13
14! routines for handling splines
15!   allocate_spline: allocates x and y vectors for splines
16!   init_spline:   generate table for spline (allocate spl%y2)
17!   spline:        return value of spline for given abscissa (optional:also y1)
18!   spline_1:      return value of 1. derivative of spline for given abscissa
19!   spline_int:    return value of integral on given interval of spline
20!   kill_spline:   destructor ( spl%x,y und/oder spl%y2)
21
22! NB: splines are always "natural splines", i.e. values of first
23!     derivative at the end-points cannot be specified
24!-----------------------------------------------------------------------------!
25      USE kinds, ONLY : DP
26
27      IMPLICIT NONE
28
29      PRIVATE
30      PUBLIC :: spline_data, allocate_spline, init_spline, spline, spline_1, &
31        spline_int, kill_spline, splineh, splinedx, splintdx, nullify_spline
32
33      TYPE spline_data
34        REAL (DP), POINTER :: x(:)  ! array containing x values
35        REAL (DP), POINTER :: y(:)  ! array containing y values
36                                     ! y(i) is the function value corresponding
37                                     ! to x(i) in the interpolation table
38        REAL (DP), POINTER :: y2(:) ! second derivative of interpolating function
39        INTEGER :: n       ! number of element in the interpolation table
40        INTEGER :: pos
41        REAL (DP) :: h, invh, h26, h16
42        REAL (DP) :: xmin, xmax  ! ... added by Carlo Cavazzoni
43      END TYPE spline_data
44
45!-----------------------------------------------------------------------------!
46
47    CONTAINS
48
49!-----------------------------------------------------------------------------!
50
51      SUBROUTINE nullify_spline( spl )
52        TYPE (spline_data), INTENT (INOUT) :: spl
53        NULLIFY( spl%x )
54        NULLIFY( spl%y )
55        NULLIFY( spl%y2 )
56        spl%n = 0
57        spl%pos = 0
58        spl%h = 0.0d0
59        spl%invh = 0.0d0
60        spl%h26 = 0.0d0
61        spl%h16 = 0.0d0
62        spl%xmin = 0.0d0
63        spl%xmax = 0.0d0
64        RETURN
65      END SUBROUTINE nullify_spline
66
67      SUBROUTINE allocate_spline( spl, nn, xmin, xmax )
68
69        IMPLICIT NONE
70
71        TYPE (spline_data), INTENT (INOUT) :: spl
72        INTEGER, INTENT (IN) :: nn
73        REAL(DP), INTENT (IN), OPTIONAL :: xmin, xmax
74
75        INTEGER err
76
77        IF( PRESENT( xmin ) .AND. .NOT. PRESENT( xmax ) ) &
78          CALL errore(' allocate_spline ', ' wrong number of arguments ', 1 )
79
80        spl%n = nn
81
82        IF ( associated(spl%x) ) THEN
83          DEALLOCATE (spl%x,STAT=err)
84          IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%x',1)
85          NULLIFY (spl%x)
86        END IF
87
88        ! note that spl%x is not allocated if we use a regular x grid
89
90        IF( PRESENT( xmin ) .AND. PRESENT( xmax ) ) THEN
91          IF( xmin >= xmax ) &
92            CALL errore(' allocate_spline ', ' wrong interval ', 1)
93          spl%xmin = xmin
94          spl%xmax = xmax
95          spl%h    = ( xmax - xmin ) / DBLE( nn - 1 )
96          spl%invh = 1.0d0 / spl%h
97        ELSE
98          spl%xmin = 0
99          spl%xmax = 0
100          ALLOCATE (spl%x(1:nn),STAT=err)
101          IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%x',1)
102        END IF
103
104        IF (associated(spl%y)) THEN
105          DEALLOCATE (spl%y,STAT=err)
106          IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%y',1)
107          NULLIFY (spl%y)
108        END IF
109
110        ALLOCATE (spl%y(1:nn),STAT=err)
111        IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%y',1)
112
113        IF (associated(spl%y2)) THEN
114          DEALLOCATE (spl%y2,STAT=err)
115          IF (err/=0) CALL errore(' allocate_spline ','could not deallocate spl%y2',1)
116          NULLIFY (spl%y2)
117        END IF
118
119        ALLOCATE (spl%y2(1:nn),STAT=err)
120        IF (err/=0) CALL errore(' allocate_spline ','could not allocate spl%y2',1)
121
122        RETURN
123      END SUBROUTINE allocate_spline
124
125
126!-----------------------------------------------------------------------
127
128      SUBROUTINE init_spline( spl, endpt, y1a, y1b )
129
130!   endpt: 's': regular spacing
131!          'l': left; 'r': right; 'b': both = specify 1-deriv for each endpoints
132
133        IMPLICIT NONE
134        TYPE (spline_data), INTENT (INOUT) :: spl
135        CHARACTER (len=*), INTENT (IN), OPTIONAL :: endpt
136        REAL (DP), INTENT (IN), OPTIONAL :: y1a, y1b  ! end point derivative
137        INTEGER :: err, i, k, n
138        REAL (DP) :: p, qn, sig, un, y1l, y1r, dyp, dym, dxp, dxm, dxpm
139        REAL (DP), POINTER :: ww(:)
140        CHARACTER (len=8) :: ep
141        LOGICAL :: reg, lep, rep
142
143        !  shortcat for regular mesh without table of x values
144
145        IF( .NOT. ASSOCIATED( spl%x ) ) THEN
146          CALL splinedx( spl%xmin, spl%xmax, spl%y(:), spl%n, 0.0d0, 0.0d0, spl%y2(:) )
147          RETURN
148        END IF
149
150        !  Find out if y first derivative is given at endpoints
151
152        IF ( .NOT. present(endpt) ) THEN
153          ep = ' '
154        ELSE
155          ep = endpt
156        END IF
157        reg = ( scan(ep,'sS') > 0 )
158        lep = ( scan(ep,'lL') > 0 ) .OR. ( scan(ep,'bB') > 0 )
159        rep = ( scan(ep,'rR') > 0 ) .OR. ( scan(ep,'bB') > 0 )
160
161        !  check input parameter consistency
162
163        IF ( ( lep .OR. rep ) .AND. .NOT. present(y1a) ) &
164          CALL errore( 'init_spline', 'first deriv. at end-point missing', 1 )
165        IF ( lep .AND. rep .AND. .NOT. present(y1b) ) &
166          CALL errore( 'init_spline', 'first deriv. at end-point missing', 1 )
167
168        !  define endpoints derivative
169
170        IF ( lep ) y1l = y1a
171        IF ( rep .AND. .NOT. lep ) y1r = y1a
172        IF ( rep .AND. lep ) y1r = y1b
173
174        spl%pos = 1
175        ALLOCATE ( ww( 1 : spl%n ), STAT = err )
176        IF (err/=0) CALL errore('init_spline','could not allocate ww',1)
177
178        n = spl % n
179
180        IF ( lep ) THEN
181          spl%y2(1) = -0.5d0
182          dxp  = spl%x(2) - spl%x(1)
183          dyp  = spl%y(2) - spl%y(1)
184          ww(1) = ( 3.0d0 / dxp ) * ( dyp / dxp - y1l )
185        ELSE
186          spl%y2(1) = 0
187          ww(1) = 0.d0
188        END IF
189
190        DO i = 2, n - 1
191
192          dxp  = spl%x(i+1) - spl%x(i)
193          dxm  = spl%x(i)   - spl%x(i-1)
194          dxpm = spl%x(i+1) - spl%x(i-1)
195
196          sig = dxm / dxpm
197          p = sig * spl%y2(i-1) + 2.0d0
198          spl%y2(i) = ( sig - 1.0d0 ) / p
199
200          dyp  = spl%y(i+1) - spl%y(i)
201          dym  = spl%y(i)   - spl%y(i-1)
202
203          ww(i) = ( 6.0d0 * ( dyp / dxp  - dym / dxm ) / dxpm  - sig * ww(i-1) ) / p
204
205        END DO
206
207        IF ( rep ) THEN
208          qn = 0.5d0
209          dxm = spl%x(n) - spl%x(n-1)
210          dym = spl%y(n) - spl%y(n-1)
211          un = ( 3.0d0 / dxm ) * ( y1r - dym / dxm )
212        ELSE
213          qn = 0
214          un = 0
215        END IF
216
217        spl % y2(n) = ( un - qn * ww(n-1) ) / ( qn * spl%y2(n-1) + 1.0d0 )
218
219        DO k = n - 1, 1, -1
220          spl % y2(k) = spl%y2(k) * spl%y2(k+1) + ww(k)
221        END DO
222
223        DEALLOCATE ( ww, STAT = err )
224        IF (err/=0) CALL errore('init_spline','could not deallocate ww',1)
225
226        IF ( reg ) THEN
227          spl%h = ( spl%x(n) - spl%x(1) ) / ( n - 1.0d0 )
228          spl%h16 = spl%h / 6.0d0
229          spl%h26 = spl%h**2 / 6.0d0
230          spl%invh = 1 / spl%h
231        ELSE
232          spl%h = 0.0d0
233          spl%invh = 0.0d0
234        END IF
235
236        RETURN
237      END SUBROUTINE init_spline
238
239!-----------------------------------------------------------------------
240
241      FUNCTION interv( spl, xx )
242
243        IMPLICIT NONE
244
245        TYPE (spline_data), INTENT (IN) :: spl
246        REAL (DP), INTENT (IN) :: xx
247        INTEGER :: interv
248        INTEGER :: khi, klo, i, p, n, k
249
250        !   if we have a regular mesh use a quick position search
251
252        IF ( spl%h /= 0 ) THEN
253          i = ( xx - spl%x(1) ) * spl%invh + 1
254          IF ( i < 1 .OR. i > spl%n ) &
255            CALL errore('interv', 'illegal x-value passed to interv',1)
256          interv = i
257          RETURN
258        END IF
259
260        p = spl%pos
261        IF ( p >= spl%n .OR. p <= 1 ) p = spl%n / 2
262        i = 0
263        n = spl%n
264
265        !   check if interval is close to previous interval
266
267        IF ( xx < spl%x(p+1) ) THEN
268          IF ( xx >= spl%x(p) ) THEN
269            i = spl%pos
270          ELSE IF ( p > 1 .AND. xx >= spl%x(p-1) ) THEN
271            i = p - 1
272          ELSE
273            klo = 1
274            khi = p + 1
275          END IF
276        ELSE IF ( (p + 2) <= n .AND. xx < spl%x(p+2) ) THEN
277          i = p + 1
278        ELSE
279          klo = p + 1
280          khi = n
281        END IF
282
283        !   perform binary search
284
285        IF ( i == 0 ) THEN
286          IF ( xx < spl%x(1) .OR. xx > spl%x(n) ) &
287            CALL errore('interv', 'xx value out of spline-range',1)
288          DO WHILE ( (khi - klo) > 1 )
289            k = ( khi + klo ) / 2
290            IF ( spl%x(k) > xx ) THEN
291              khi = k
292            ELSE
293              klo = k
294            END IF
295          END DO
296          i = klo
297        END IF
298
299        interv = i
300        RETURN
301      END FUNCTION interv
302
303
304!-----------------------------------------------------------------------
305      FUNCTION spline( spl, xx, y1 )
306
307        IMPLICIT NONE
308
309        TYPE (spline_data), INTENT (INOUT) :: spl
310        REAL (DP), INTENT (IN) :: xx
311        REAL (DP), INTENT (OUT), OPTIONAL :: y1
312        REAL (DP) :: spline
313
314        INTEGER :: khi, klo
315        REAL (DP) :: a, b, h, invh, ylo, yhi, y2lo, y2hi, a3ma, b3mb
316
317        !  shortcat for regular mesh without table of x values
318
319        IF( .NOT. ASSOCIATED( spl%x ) ) THEN
320          IF( PRESENT( y1 ) ) &
321            CALL errore(' spline ', ' y1 without x table not implemented ', 1 )
322          CALL splintdx( spl%xmin, spl%xmax, spl%y, spl%y2, spl%n, xx, a )
323          spline = a
324          RETURN
325        END IF
326
327        spl%pos = interv( spl, xx )
328        klo = spl%pos
329        khi = spl%pos + 1
330
331        IF ( spl%h /= 0 ) THEN
332          h    = spl%h
333          invh = spl%invh
334        ELSE
335          h    = spl%x( khi ) - spl%x( klo )
336          invh = spl%invh
337          IF ( h == 0.0d0 ) &
338            CALL errore('spline','bad spl%x input',1)
339        END IF
340
341        a = ( spl%x( khi ) - xx ) * invh
342        b = 1 - a
343        a3ma = a**3 - a
344        b3mb = b**3 - b
345        ylo  = spl%y( klo )
346        yhi  = spl%y( khi )
347        y2lo = spl%y2( klo )
348        y2hi = spl%y2( khi )
349        spline = a * ylo + b * yhi + ( a3ma * y2lo + b3mb * y2hi ) * ( h**2 ) / 6.0d0
350
351        IF ( present( y1 ) ) then
352          y1 = ( yhi - ylo ) * invh + &
353          ( ( 1.0d0 - 3 * a**2 ) * y2lo + ( 3 * b**2 - 1.0d0 ) * y2hi ) * h / 6.0d0
354        END IF
355
356        RETURN
357      END FUNCTION spline
358
359!-----------------------------------------------------------------------
360
361      FUNCTION splineh(spl,xx,y1)
362        IMPLICIT NONE
363        TYPE (spline_data), INTENT (IN) :: spl
364        REAL (DP), INTENT (IN) :: xx
365        REAL (DP), INTENT (OUT) :: y1
366        REAL (DP) :: splineh
367
368        INTEGER :: khi, klo, i
369        REAL (DP) :: a, b, h, invh, t, ylo, yhi, y2lo, y2hi, d, d0
370
371! fast spline for pair potentials without checks
372        h = spl%h
373        invh = spl%invh
374        d=xx-spl%x(1); i=INT(d*spl%invh); d0=DBLE(i)*h; i=i+1
375        i = (xx-spl%x(1))*invh + 1
376
377        a = (spl%x(i+1)-xx)*invh
378        b = 1 - a
379        t = -a*b
380!    b=(d-d0)*invh; a=1-b; t=-a*b
381        ylo = spl%y(i)
382        yhi = spl%y(i+1)
383        y2lo = spl%y2(i)
384        y2hi = spl%y2(i+1)
385        splineh = a*ylo + b*yhi + ((a+1)*y2lo+(b+1)*y2hi)*t*spl%h26
386        y1 = (yhi-ylo)*invh + ((1.d0-3*a*a)*y2lo+(3*b*b-1.d0)*y2hi)*spl%h16
387
388      END FUNCTION splineh
389!-----------------------------------------------------------------------
390      FUNCTION spline_1(spl,xx)
391        IMPLICIT NONE
392        TYPE (spline_data), INTENT (INOUT) :: spl
393        REAL (DP), INTENT (IN) :: xx
394        REAL (DP) :: spline_1
395
396        INTEGER :: khi, klo
397        REAL (DP) :: a, b, h
398
399        spl%pos = interv(spl,xx)
400        klo = spl%pos
401        khi = spl%pos + 1
402
403        h = spl%x(khi) - spl%x(klo)
404        IF (h==0.d0) CALL errore('spline','bad spl%x input',1)
405        a = (spl%x(khi)-xx)/h
406        b = 1 - a
407        spline_1 = (spl%y(khi)-spl%y(klo))/h + ((1.d0-3*a**2)*spl%y2(klo)+(3*b** &
408          2-1.d0)*spl%y2(khi))*h/6.d0
409
410        RETURN
411      END FUNCTION spline_1
412
413
414
415!-----------------------------------------------------------------------
416      FUNCTION stamm(spl,p,x)
417        IMPLICIT NONE
418        TYPE (spline_data), INTENT (IN) :: spl
419        INTEGER, INTENT (IN) :: p
420        REAL (DP), INTENT (IN) :: x
421        REAL (DP) :: stamm
422        REAL (DP) :: a, b, aa, bb, h
423
424        h = spl%x(p+1) - spl%x(p)
425        b = (x-spl%x(p))/h
426        a = 1 - b
427        aa = a**2
428        bb = b**2
429        stamm = 0.5d0*h*(bb*spl%y(p+1)-aa*spl%y(p)) + h**3/12.d0*(aa*(1-0.5d0*aa)* &
430          spl%y2(p)-bb*(1-0.5d0*bb)*spl%y2(p+1))
431
432        RETURN
433      END FUNCTION stamm
434
435
436
437!-----------------------------------------------------------------------
438      FUNCTION spline_int(spl,x0,x1)
439        IMPLICIT NONE
440        TYPE (spline_data), INTENT (INOUT) :: spl
441        REAL (DP), INTENT (IN) :: x0, x1
442        REAL (DP) :: spline_int
443
444        INTEGER :: j, pa, pb
445        REAL (DP) :: h, vorz, xa, xb, i1, i2
446
447        vorz = 1
448        xa = min(x0,x1)
449        xb = max(x0,x1)
450        IF (x0>x1) vorz = -1
451        IF (xa<spl%x(1) .OR. xb>spl%x(spl%n)) CALL errore('spline_int', &
452          'illegal integration range',1)
453
454        pa = interv(spl,xa)
455        pb = interv(spl,xb)
456
457        IF (pa==pb) THEN
458          spline_int = vorz*(stamm(spl,pa,xb)-stamm(spl,pa,xa))
459          RETURN
460        END IF
461
462        i1 = 0
463        i2 = 0
464        DO j = pa + 1, pb - 1
465          h = spl%x(j+1) - spl%x(j)
466          i1 = i1 + h*(spl%y(j)+spl%y(j+1))
467          i2 = i2 + h**3*(spl%y2(j)+spl%y2(j+1))
468        END DO
469        h = spl%x(pa+1) - spl%x(pa)
470        i1 = i1 + h*spl%y(pa+1)
471        i2 = i2 + h**3*spl%y2(pa+1)
472        h = spl%x(pb+1) - spl%x(pb)
473        i1 = i1 + h*spl%y(pb)
474        i2 = i2 + h**3*spl%y2(pb)
475
476        spline_int = vorz*(i1/2.-i2/24.d0+stamm(spl,pb,xb)-stamm(spl,pa,xa))
477
478        RETURN
479      END FUNCTION spline_int
480
481!-----------------------------------------------------------------------
482      SUBROUTINE kill_spline(spl,what)
483!   deallocate splines
484!   what=='a' or not present: deallocate all (spl%x, spl%y, spl%y2)
485!   what=='d': deallocate only data (spl%x, spl%y)
486!   what=='2': deallocate only table of 2. derivatives (spl%y2)
487        IMPLICIT NONE
488        TYPE (spline_data), INTENT (INOUT) :: spl
489        CHARACTER, INTENT (IN), OPTIONAL :: what
490        CHARACTER :: w
491        INTEGER :: err
492
493        w = 'a'
494        IF (present(what)) w = what
495        SELECT CASE (w)
496        CASE ('d','D')
497          IF (associated(spl%x)) THEN
498            DEALLOCATE (spl%x,STAT=err)
499            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%x',1)
500            NULLIFY (spl%x)
501          END IF
502          IF (associated(spl%y)) THEN
503            DEALLOCATE (spl%y,STAT=err)
504            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y',1)
505            NULLIFY (spl%y)
506          END IF
507        CASE ('2')
508          IF (associated(spl%y2)) THEN
509            DEALLOCATE (spl%y2,STAT=err)
510            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y2',1)
511            NULLIFY (spl%y2)
512          END IF
513        CASE ('a','A')
514          IF (associated(spl%x)) THEN
515            DEALLOCATE (spl%x,STAT=err)
516            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%x',1)
517            NULLIFY (spl%x)
518          END IF
519          IF (associated(spl%y)) THEN
520            DEALLOCATE (spl%y,STAT=err)
521            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y',1)
522            NULLIFY (spl%y)
523          END IF
524          IF (associated(spl%y2)) THEN
525            DEALLOCATE (spl%y2,STAT=err)
526            IF (err/=0) CALL errore('kill_spline', 'could not deallocate spl%y2',1)
527            NULLIFY (spl%y2)
528          END IF
529        END SELECT
530
531        RETURN
532      END SUBROUTINE kill_spline
533
534!=-----------------------------------------------------------------------=!
535! Subroutines: splinedx, splintdx
536! added for compatibility with SISSA code
537! Carlo Cavazzoni 15-03-2000
538!=-----------------------------------------------------------------------=!
539
540      SUBROUTINE splinedx(xmin,xmax,y,n,yp1,ypn,y2)
541        USE kinds
542        IMPLICIT NONE
543        INTEGER, INTENT(IN)  :: n
544        REAL(DP),  INTENT(IN)  :: yp1,ypn,xmin,xmax,y(:)
545        REAL(DP),  INTENT(OUT) :: y2(:)
546        INTEGER :: i, k
547        REAL(DP)  :: p, qn, sig, un, dx
548        REAL(DP)  :: u(n)
549
550        dx = (xmax-xmin)/DBLE(n-1)
551        if ( yp1 .gt. 0.99d30 ) then
552          y2(1)=0.d0
553          u(1)=0.0d0
554        else
555          y2(1)=-0.5d0
556          u(1)=(3.d0/dx) * ( (y(2)-y(1)) / dx - yp1 )
557        endif
558        do i=2,n-1
559          sig=0.5d0
560          p=sig*y2(i-1)+2.d0
561          y2(i)=(sig-1.d0)/p
562          u(i) = (6.0d0 * ( (y(i+1)-y(i))/ dx - (y(i)-y(i-1))/ dx ) &
563                        / (2.0d0*dx) - sig * u(i-1) ) / p
564        end do
565        if ( ypn .gt. 0.99d30 ) then
566          qn=0.d0
567          un=0.d0
568        else
569          qn=0.5d0
570          un= ( 3.d0 / dx ) * ( ypn - (y(n)-y(n-1)) / dx )
571        endif
572        y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0d0)
573        do k=n-1,1,-1
574          y2(k)=y2(k)*y2(k+1)+u(k)
575        end do
576        return
577      END SUBROUTINE splinedx
578
579
580      SUBROUTINE splintdx(xmin,xmax,ya,y2a,n,x,y)
581        USE kinds
582        IMPLICIT NONE
583        INTEGER, INTENT(IN)  :: n
584        REAL(DP),  INTENT(IN)  :: x,xmin,xmax,ya(:),y2a(:)
585        REAL(DP),  INTENT(OUT) :: y
586        INTEGER :: khi,klo
587        REAL(DP)  :: a,b,h,dx,xhi,xlo
588        dx  = (xmax-xmin)/DBLE(n-1)
589        klo = INT(x/dx)
590        khi = klo+1
591        IF(klo.LT.1) THEN
592          CALL errore(' splintdx ',' klo less than one ',klo)
593        END IF
594        IF(khi.GT.n) THEN
595          CALL errore(' splintdx ',' khi grether than N ',khi)
596        END IF
597        xlo  = xmin + DBLE(klo-1) * dx
598        xhi  = xmin + DBLE(khi-1) * dx
599
600        a   = (xhi-x)/dx
601        b   = (x-xlo)/dx
602        y   = a*ya(klo) + b*ya(khi) + &
603              ( (a*a*a-a)*y2a(klo) + (b*b*b-b)*y2a(khi) ) * (dx*dx)/6.0d0
604        RETURN
605      END SUBROUTINE splintdx
606
607!-----------------------------------------------------------------------
608
609      SUBROUTINE nr_spline( x, y, n, yp1, ypn, y2 )
610      INTEGER :: n
611      REAL(DP) :: yp1, ypn, x(n), y(n), y2(n)
612      INTEGER :: i, k
613      REAL(DP) :: p, qn, sig, un
614      REAL(DP) :: u( n )
615      if ( yp1 .gt. 0.99d30 ) then
616        y2(1)=0.d0
617        u(1)=0.d0
618      else
619        y2(1)=-0.5d0
620        u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
621      endif
622      do i=2,n-1
623        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
624        p=sig*y2(i-1)+2.d0
625        y2(i)=(sig-1.d0)/p
626        u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / &
627             (x(i)-x(i-1))) / (x(i+1)-x(i-1))-sig*u(i-1))/p
628      end do
629      if ( ypn .gt. 0.99d30 ) then
630        qn=0.d0
631        un=0.d0
632      else
633        qn=0.5d0
634        un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
635      endif
636      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
637      do k=n-1,1,-1
638        y2(k)=y2(k)*y2(k+1)+u(k)
639      end do
640      return
641      END SUBROUTINE nr_spline
642
643
644      SUBROUTINE nr_splint( xa, ya, y2a, n, x, y )
645      INTEGER :: n
646      REAL(DP) :: x,y,xa(n),y2a(n),ya(n)
647      INTEGER :: k,khi,klo
648      REAL(DP) :: a,b,h
649      klo=1
650      khi=n
6511     if (khi-klo.gt.1) then
652        k=(khi+klo)/2
653        if(xa(k).gt.x)then
654          khi=k
655        else
656          klo=k
657        endif
658      goto 1
659      endif
660      h=xa(khi)-xa(klo)
661      if (h.eq.0.) &
662        CALL errore(' nr_splint ', 'bad xa input in splint' , 1 )
663      a=(xa(khi)-x)/h
664      b=(x-xa(klo))/h
665      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))* &
666        (h**2)/6.d0
667      return
668      END SUBROUTINE nr_splint
669
670
671      SUBROUTINE nr_splie2( x1a, x2a, ya, m, n, y2a )
672      INTEGER :: m, n
673      REAL(DP) :: x1a(m), x2a(n), y2a(m,n), ya(m,n)
674      INTEGER :: j,k
675      REAL(DP) :: y2tmp(n), ytmp(n)
676      do j = 1, m
677        do k = 1, n
678          ytmp(k) = ya(j,k)
679        end do
680        call nr_spline( x2a, ytmp, n, 1.d30, 1.d30, y2tmp )
681        do k = 1, n
682          y2a(j,k) = y2tmp(k)
683        end do
684      end do
685      return
686      END SUBROUTINE nr_splie2
687
688
689      SUBROUTINE nr_splin2( x1a, x2a, ya, y2a, m, n, x1, x2, y )
690      INTEGER :: m, n
691      REAL(DP) :: x1, x2, y, x1a(m), x2a(n), y2a(m,n), ya(m,n)
692      INTEGER :: j, k
693      REAL(DP) :: y2tmp( MAX(n,m) ), ytmp( n ), yytmp( MAX(n,m) )
694      do j = 1, m
695        do k = 1, n
696          ytmp(k)  = ya(j,k)
697          y2tmp(k) = y2a(j,k)
698        end do
699        call nr_splint( x2a, ytmp, y2tmp, n, x2, yytmp(j) )
700      end do
701      call nr_spline( x1a, yytmp, m, 1.d30, 1.d30, y2tmp )
702      call nr_splint( x1a, yytmp, y2tmp, m, x1, y )
703      return
704      END SUBROUTINE nr_splin2
705
706!-----------------------------------------------------------------------
707
708    END MODULE splines
709