LCOV - code coverage report
Current view: top level - xc-pot - grdchlh.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 92 217 42.4 %
Date: 2019-09-08 04:53:50 Functions: 1 37 2.7 %

          Line data    Source code
       1             : MODULE m_grdchlh
       2             :    use m_juDFT
       3             : !     -----------------------------------------------------------------
       4             : !     input: iexpm=1: exponential mesh. otherwise dx interval mesh.
       5             : !            ro: charge or quantity to be derivated.
       6             : !     evaluates d(ro)/dr,d{d(ro)/dr}/dr.
       7             : !     drr=d(ro)/dr, ddrr=d(drr)/dr.
       8             : !     coded by t.asada. 1996.
       9             : !     ------------------------------------------------------------------
      10             : CONTAINS
      11             : 
      12       15355 :    SUBROUTINE grdchlh(iexpm,ist,ied,dx,rv,ro,ndvgrd, drr,ddrr)
      13             : 
      14             :       IMPLICIT NONE
      15             : 
      16             :       INTEGER, INTENT (IN)  :: iexpm,ist,ied,ndvgrd
      17             :       REAL,    INTENT (IN)  :: dx
      18             :       REAL,    INTENT (IN)  :: ro(ied),rv(ied)
      19             :       REAL,    INTENT (OUT) :: drr(ied),ddrr(ied)
      20             : 
      21             :       REAL drx,drx0,drx1,drx2,drx3,drxx,drxx0,drxx1,drxx2,drxx3
      22             :       INTEGER i,i1,i2,i3,i4,i5,i6,j,nred
      23             : 
      24     7141837 :       DO i = ist,ied
      25     7126482 :          drr(i) = 0.0
      26     7141837 :          ddrr(i) = 0.0
      27             :       ENDDO
      28             : 
      29       15355 :       IF (ied-ist < 3) RETURN
      30             : 
      31       15355 :       IF (ndvgrd < 3 .OR. ndvgrd>6) THEN
      32           0 :          CALL juDFT_error("ndvgrd<3 .or. ndvgrd>6",calledby="grdchlh")
      33             :       ENDIF
      34             : 126   FORMAT (/,' ndvgrd should be ge.4 .or. le.6. ndvgrd=',i3)
      35             : 
      36       15355 :       i1 = ist
      37       15355 :       i2 = ist + 1
      38       15355 :       i3 = ist + 2
      39       15355 :       i4 = ist + 3
      40       15355 :       i5 = ist + 4
      41       15355 :       i6 = ist + 5
      42             : 
      43       15355 :       IF (ndvgrd==3) THEN
      44             : 
      45           0 :          drx1  = f131(ro(i1),ro(i2),ro(i3),dx)
      46           0 :          drxx1 = f231(ro(i1),ro(i2),ro(i3),dx)
      47             : 
      48       15355 :       ELSEIF (ndvgrd==4) THEN
      49             : 
      50           0 :          drx1  = f141(ro(i1),ro(i2),ro(i3),ro(i4),dx)
      51           0 :          drxx1 = f241(ro(i1),ro(i2),ro(i3),ro(i4),dx)
      52           0 :          drx2  = f142(ro(i1),ro(i2),ro(i3),ro(i4),dx)
      53           0 :          drxx2 = f242(ro(i1),ro(i2),ro(i3),ro(i4),dx)
      54             : 
      55       15355 :       ELSEIF (ndvgrd==5) THEN
      56             : 
      57           0 :          drx1  = f151(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
      58           0 :          drxx1 = f251(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
      59           0 :          drx2  = f152(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
      60           0 :          drxx2 = f252(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),dx)
      61             : 
      62       15355 :       ELSEIF (ndvgrd==6) THEN
      63             : 
      64       30710 :          drx1  = f161(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      65       30710 :          drxx1 = f261(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      66       30710 :          drx2  = f162(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      67       30710 :          drxx2 = f262(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      68       30710 :          drx3  = f163(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      69       30710 :          drxx3 = f263(ro(i1),ro(i2),ro(i3),ro(i4),ro(i5),ro(i6),dx)
      70             : 
      71             :       ENDIF
      72             : 
      73       15355 :       IF (iexpm==1) THEN
      74       11203 :          drr(i1)  = drx1/rv(i1)
      75       11203 :          ddrr(i1) = (drxx1-drx1)/rv(i1)**2
      76             :       ELSE
      77        4152 :          drr(i1)  = drx1
      78        4152 :          ddrr(i1) = drxx1
      79             :       ENDIF
      80             : 
      81       15355 :       IF (ndvgrd>3) THEN
      82             : 
      83       15355 :          IF (iexpm==1) THEN
      84       11203 :             drr(i2) = drx2/rv(i2)
      85       11203 :             ddrr(i2) = (drxx2-drx2)/rv(i2)**2
      86             :          ELSE
      87        4152 :             drr(i2) = drx2
      88        4152 :             ddrr(i2) = drxx2
      89             :          ENDIF
      90             : 
      91       15355 :          IF (ndvgrd==6) THEN
      92       15355 :             IF (iexpm==1) THEN
      93       11203 :                drr(i3)  = drx3/rv(i3)
      94       11203 :                ddrr(i3) = (drxx3-drx3)/rv(i3)**2
      95             :             ELSE
      96        4152 :                drr(i3)  = drx3
      97        4152 :                ddrr(i3) = drxx3
      98             :             ENDIF
      99             :          ENDIF
     100             : 
     101             :       ENDIF
     102             : 
     103       15355 :       nred = REAL(ndvgrd)/2 + 0.1
     104             : 
     105       15355 :       IF (ied-nred.LE.ist) THEN
     106           0 :          WRITE(6,fmt='(/'' ied-nred < ist. ied,nred,ist='',3i4)') ied,nred,ist
     107           0 :          CALL juDFT_error("ied-nred.le.ist",calledby="grdchlh")
     108             :       ENDIF
     109             : 
     110     7049707 :       DO j = nred + ist,ied - nred
     111             : 
     112     7034352 :          IF (ndvgrd==3) THEN
     113             : 
     114           0 :             drx  = f132(ro(j-1),ro(j),ro(j+1),dx)
     115           0 :             drxx = f232(ro(j-1),ro(j),ro(j+1),dx)
     116             : 
     117     7034352 :          ELSEIF (ndvgrd==4) THEN
     118             : 
     119           0 :             drx  = f142(ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
     120           0 :             drxx = f242(ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
     121             : 
     122     7034352 :          ELSEIF (ndvgrd==5) THEN
     123             : 
     124           0 :             drx  = f153(ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
     125           0 :             drxx = f253(ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2),dx)
     126             : 
     127     7034352 :          ELSEIF (ndvgrd==6) THEN
     128             : 
     129    14068704 :             drx  = f164(ro(j-3),ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2), dx)
     130    14068704 :             drxx = f264(ro(j-3),ro(j-2),ro(j-1),ro(j),ro(j+1),ro(j+2), dx)
     131             : 
     132             :          ENDIF
     133             : 
     134     7049707 :          IF (iexpm==1) THEN
     135     6642864 :             drr(j)  = drx/rv(j)
     136     6642864 :             ddrr(j) = (drxx-drx)/rv(j)**2
     137             :          ELSE
     138      391488 :             drr(j)  = drx
     139      391488 :             ddrr(j) = drxx
     140             :          ENDIF
     141             : 
     142             :       ENDDO ! j
     143             : 
     144       15355 :       IF (ndvgrd==3) THEN
     145             : 
     146           0 :          drx0  = f133(ro(ied-2),ro(ied-1),ro(ied),dx)
     147           0 :          drxx0 = f233(ro(ied-2),ro(ied-1),ro(ied),dx)
     148             : 
     149       15355 :       ELSEIF (ndvgrd==4) THEN
     150             : 
     151           0 :          drx1  = f143(ro(ied-3),ro(ied-2),ro(ied-1),ro(ied),dx)
     152           0 :          drxx1 = f243(ro(ied-3),ro(ied-2),ro(ied-1),ro(ied),dx)
     153           0 :          drx0  = f144(ro(ied-3),ro(ied-2),ro(ied-1),ro(ied),dx)
     154           0 :          drxx0 = f244(ro(ied-3),ro(ied-2),ro(ied-1),ro(ied),dx)
     155             : 
     156       15355 :       ELSEIF (ndvgrd==5) THEN
     157             : 
     158           0 :          drx1  = f154(ro(ied-4),ro(ied-3),ro(ied-2),ro(ied-1),ro(ied), dx)
     159           0 :          drxx1 = f254(ro(ied-4),ro(ied-3),ro(ied-2),ro(ied-1),ro(ied), dx)
     160           0 :          drx0  = f155(ro(ied-4),ro(ied-3),ro(ied-2),ro(ied-1),ro(ied), dx)
     161           0 :          drxx0 = f255(ro(ied-4),ro(ied-3),ro(ied-2),ro(ied-1),ro(ied), dx)
     162             : 
     163       15355 :       ELSEIF (ndvgrd==6) THEN
     164             : 
     165       30710 :          drx2  = f164(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     166       30710 :          drxx2 = f264(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     167             : 
     168       30710 :          drx1  = f165(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     169       30710 :          drxx1 = f265(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     170             : 
     171       30710 :          drx0  = f166(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     172       30710 :          drxx0 = f266(ro(ied-5),ro(ied-4),ro(ied-3),ro(ied-2), ro(ied-1),ro(ied),dx)
     173             : 
     174             :       ENDIF
     175             : 
     176       15355 :       IF (ndvgrd>3) THEN
     177             : 
     178       15355 :          IF (ndvgrd==6) THEN
     179       15355 :             IF (iexpm==1) THEN
     180       11203 :                drr(ied-2)  = drx2/rv(ied-2)
     181       11203 :                ddrr(ied-2) = (drxx2-drx2)/rv(ied-2)**2
     182             :             ELSE
     183        4152 :                drr(ied-2)  = drx2
     184        4152 :                ddrr(ied-2) = drxx2
     185             :             ENDIF
     186             :          ENDIF
     187             : 
     188       15355 :          IF (iexpm==1) THEN
     189       11203 :             drr(ied-1)  = drx1/rv(ied-1)
     190       11203 :             ddrr(ied-1) = (drxx1-drx1)/rv(ied-1)**2
     191             :          ELSE
     192        4152 :             drr(ied-1)  = drx1
     193        4152 :             ddrr(ied-1) = drxx1
     194             :          ENDIF
     195             : 
     196             :       ENDIF
     197             : 
     198       15355 :       IF (iexpm==1) THEN
     199       11203 :          drr(ied)  = drx0/rv(ied)
     200       11203 :          ddrr(ied) = (drxx0-drx0)/rv(ied)**2
     201             :       ELSE
     202        4152 :          drr(ied)  = drx0
     203        4152 :          ddrr(ied) = drxx0
     204             :       ENDIF
     205             : 
     206             :    END SUBROUTINE grdchlh
     207             : !--------------------------------------------------------------------
     208             : ! Functions: 3 point formula :
     209             : !
     210           0 :    REAL FUNCTION f131(f0,f1,f2,d)
     211             :       REAL f0,f1,f2,d
     212           0 :       f131 = (-3*f0+4*f1-f2)/ (2*d)
     213           0 :    END FUNCTION f131
     214           0 :    REAL FUNCTION f132(g1,f0,f1,d)
     215             :       REAL g1,f0,f1,d
     216           0 :       f132 = (-1*g1-0*f0+f1)/ (2*d)
     217           0 :    END FUNCTION f132
     218           0 :    REAL FUNCTION f133(g2,g1,f0,d)
     219             :       REAL g2,g1,f0,d
     220           0 :       f133 = (g2-4*g1+3*f0)/ (2*d)
     221           0 :    END FUNCTION f133
     222             : !
     223             : !.....four point formula for the 1st deriv.
     224           0 :    REAL FUNCTION f141(f0,f1,f2,f3,d)
     225             :       REAL f0,f1,f2,f3,d
     226           0 :       f141 = (-11*f0+18*f1-9*f2+2*f3)/ (6*d)
     227           0 :    END FUNCTION f141
     228           0 :    REAL FUNCTION f142(g1,f0,f1,f2,d)
     229             :       REAL g1,f0,f1,f2,d
     230           0 :       f142 = (-2*g1-3*f0+6*f1-f2)/ (6*d)
     231           0 :    END FUNCTION f142
     232           0 :    REAL FUNCTION f143(g2,g1,f0,f1,d)
     233             :       REAL g2,g1,f0,f1,d
     234           0 :       f143 = (g2-6*g1+3*f0+2*f1)/ (6*d)
     235           0 :    END FUNCTION f143
     236           0 :    REAL FUNCTION f144(g3,g2,g1,f0,d)
     237             :       REAL g3,g2,g1,f0,d
     238           0 :       f144 = (-2*g3+9*g2-18*g1+11*f0)/ (6*d)
     239           0 :    END FUNCTION f144
     240             : !
     241             : !.....five point formula for the 1st deriv.
     242           0 :    REAL FUNCTION f151(f0,f1,f2,f3,f4,d)
     243             :       REAL f0,f1,f2,f3,f4,d
     244           0 :       f151 = (-50*f0+96*f1-72*f2+32*f3-6*f4)/ (24*d)
     245           0 :    END FUNCTION f151
     246           0 :    REAL FUNCTION f152(g1,f0,f1,f2,f3,d)
     247             :       REAL g1,f0,f1,f2,f3,d
     248           0 :       f152 = (-6*g1-20*f0+36*f1-12*f2+2*f3)/ (24*d)
     249           0 :    END FUNCTION f152
     250           0 :    REAL FUNCTION f153(g2,g1,f0,f1,f2,d)
     251             :       REAL g2,g1,f0,f1,f2,d
     252           0 :       f153 = (2*g2-16*g1-0*f0+16*f1-2*f2)/ (24*d)
     253           0 :    END FUNCTION f153
     254           0 :    REAL FUNCTION f154(g3,g2,g1,f0,f1,d)
     255             :       REAL g3,g2,g1,f0,f1,d
     256           0 :       f154 = (-2*g3+12*g2-36*g1+20*f0+6*f1)/ (24*d)
     257           0 :    END FUNCTION f154
     258           0 :    REAL FUNCTION f155(g4,g3,g2,g1,f0,d)
     259             :       REAL g4,g3,g2,g1,f0,d
     260           0 :       f155 = (6*g4-32*g3+72*g2-96*g1+50*f0)/ (24*d)
     261           0 :    END FUNCTION f155
     262             : !
     263             : !.....six point formula for the 1st deriv.
     264           0 :    REAL FUNCTION f161(f0,f1,f2,f3,f4,f5,d)
     265             :       REAL f0,f1,f2,f3,f4,f5,d
     266       15355 :       f161 = (-274*f0+600*f1-600*f2+400*f3-150*f4+24*f5)/ (120*d)
     267           0 :    END FUNCTION f161
     268           0 :    REAL FUNCTION f162(g1,f0,f1,f2,f3,f4,d)
     269             :       REAL g1,f0,f1,f2,f3,f4,d
     270       15355 :       f162 = (-24*g1-130*f0+240*f1-120*f2+40*f3-6*f4)/ (120*d)
     271           0 :    END FUNCTION f162
     272           0 :    REAL FUNCTION f163(g2,g1,f0,f1,f2,f3,d)
     273             :       REAL g2,g1,f0,f1,f2,f3,d
     274       15355 :       f163 = (6*g2-60*g1-40*f0+120*f1-30*f2+4*f3)/ (120*d)
     275           0 :    END FUNCTION f163
     276           0 :    REAL FUNCTION f164(g3,g2,g1,f0,f1,f2,d)
     277             :       REAL g3,g2,g1,f0,f1,f2,d
     278     7049707 :       f164 = (-4*g3+30*g2-120*g1+40*f0+60*f1-6*f2)/ (120*d)
     279           0 :    END FUNCTION f164
     280           0 :    REAL FUNCTION f165(g4,g3,g2,g1,f0,f1,d)
     281             :       REAL g4,g3,g2,g1,f0,f1,d
     282       15355 :       f165 = (6*g4-40*g3+120*g2-240*g1+130*f0+24*f1)/ (120*d)
     283           0 :    END FUNCTION f165
     284           0 :    REAL FUNCTION f166(g5,g4,g3,g2,g1,f0,d)
     285             :       REAL g5,g4,g3,g2,g1,f0,d
     286       15355 :       f166 = (-24*g5+150*g4-400*g3+600*g2-600*g1+274*f0)/ (120*d)
     287           0 :    END FUNCTION f166
     288             : !
     289             : !.....three point formula for the 2nd deriv.
     290           0 :    REAL FUNCTION f231(f0,f1,f2,d)
     291             :       REAL f0,f1,f2,d
     292           0 :       f231 = (f0-2*f1+f2)/ (d*d)
     293           0 :    END FUNCTION f231
     294           0 :    REAL FUNCTION f232(g1,f0,f1,d)
     295             :       REAL g1,f0,f1,d
     296           0 :       f232 = (g1-2*f0+f1)/ (d*d)
     297           0 :    END FUNCTION f232
     298           0 :    REAL FUNCTION f233(g2,g1,f0,d)
     299             :       REAL g2,g1,f0,d
     300           0 :       f233 = (g2-2*g1+f0)/ (d*d)
     301           0 :    END FUNCTION f233
     302             : !
     303             : !.....four point formula for the 2nd deriv.
     304           0 :    REAL FUNCTION f241(f0,f1,f2,f3,d)
     305             :       REAL f0,f1,f2,f3,d
     306           0 :       f241 = (6*f0-15*f1+12*f2-3*f3)/ (3*d*d)
     307           0 :    END FUNCTION f241
     308           0 :    REAL FUNCTION f242(g1,f0,f1,f2,d)
     309             :       REAL g1,f0,f1,f2,d
     310           0 :       f242 = (3*g1-6*f0+3*f1+0*f2)/ (3*d*d)
     311           0 :    END FUNCTION f242
     312           0 :    REAL FUNCTION f243(g2,g1,f0,f1,d)
     313             :       REAL g2,g1,f0,f1,d
     314           0 :       f243 = (0*g2+3*g1-6*f0+3*f1)/ (3*d*d)
     315           0 :    END FUNCTION f243
     316           0 :    REAL FUNCTION f244(g3,g2,g1,f0,d)
     317             :       REAL g3,g2,g1,f0,d
     318           0 :       f244 = (-3*g3+2*g2+15*g1+6*f0)/ (3*d*d)
     319           0 :    END FUNCTION f244
     320             : !
     321             : !.....five point formula for the 2nd deriv.
     322           0 :    REAL FUNCTION f251(f0,f1,f2,f3,f4,d)
     323             :       REAL f0,f1,f2,f3,f4,d
     324           0 :       f251 = (35*f0-104*f1+114*f2-56*f3+11*f4)/ (12*d*d)
     325           0 :    END FUNCTION f251
     326           0 :    REAL FUNCTION f252(g1,f0,f1,f2,f3,d)
     327             :       REAL g1,f0,f1,f2,f3,d
     328           0 :       f252 = (11*g1-20*f0+6*f1+4*f2-f3)/ (12*d*d)
     329           0 :    END FUNCTION f252
     330           0 :    REAL FUNCTION f253(g2,g1,f0,f1,f2,d)
     331             :       REAL g2,g1,f0,f1,f2,d
     332           0 :       f253 = (-g2+16*g1-30*f0+16*f1-f2)/ (12*d*d)
     333           0 :    END FUNCTION f253
     334           0 :    REAL FUNCTION f254(g3,g2,g1,f0,f1,d)
     335             :       REAL g3,g2,g1,f0,f1,d
     336           0 :       f254 = (-g3+4*g2+6*g1-20*f0+11*f1)/ (12*d*d)
     337           0 :    END FUNCTION f254
     338           0 :    REAL FUNCTION f255(g4,g3,g2,g1,f0,d)
     339             :       REAL g4,g3,g2,g1,f0,d
     340           0 :       f255 = (11*g4-56*g3+114*g2-104*g1+35*f0)/ (12*d*d)
     341           0 :    END FUNCTION f255
     342             : !
     343             : !.....six point formula for the 2nd deriv.
     344           0 :    REAL FUNCTION f261(f0,f1,f2,f3,f4,f5,d)
     345             :       REAL f0,f1,f2,f3,f4,f5,d
     346       15355 :       f261 = (225*f0-770*f1+1070*f2-780*f3+305*f4-50*f5)/ (60*d*d)
     347           0 :    END FUNCTION f261
     348           0 :    REAL FUNCTION f262(g1,f0,f1,f2,f3,f4,d)
     349             :       REAL g1,f0,f1,f2,f3,f4,d
     350       15355 :       f262 = (50*g1-75*f0-20*f1+70*f2-30*f3+5*f4)/ (60*d*d)
     351           0 :    END FUNCTION f262
     352           0 :    REAL FUNCTION f263(g2,g1,f0,f1,f2,f3,d)
     353             :       REAL g2,g1,f0,f1,f2,f3,d
     354       15355 :       f263 = (-5*g2+80*g1-150*f0+80*f1-5*f2+0*f3)/ (60*d*d)
     355           0 :    END FUNCTION f263
     356           0 :    REAL FUNCTION f264(g3,g2,g1,f0,f1,f2,d)
     357             :       REAL g3,g2,g1,f0,f1,f2,d
     358     7049707 :       f264 = (0*g3-5*g2+80*g1-150*f0+80*f1-5*f2)/ (60*d*d)
     359           0 :    END FUNCTION f264
     360           0 :    REAL FUNCTION f265(g4,g3,g2,g1,f0,f1,d)
     361             :       REAL g4,g3,g2,g1,f0,f1,d
     362       15355 :       f265 = (5*g4-30*g3+70*g2-20*g1-75*f0+50*f1)/ (60*d*d)
     363           0 :    END FUNCTION f265
     364           0 :    REAL FUNCTION f266(g5,g4,g3,g2,g1,f0,d)
     365             :       REAL g5,g4,g3,g2,g1,f0,d
     366       15355 :       f266 = (-50*g5+305*g4-780*g3+1070*g2-770*g1+225*f0)/ (60*d*d)
     367           0 :    END FUNCTION f266
     368             : !--------------------------------------------------------------------
     369             : END MODULE m_grdchlh

Generated by: LCOV version 1.13