LCOV - code coverage report
Current view: top level - core - core_.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 208 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       MODULE m_core
       2             : c...............................................................core
       3             :       CONTAINS
       4           0 :       SUBROUTINE core(
       5           0 :      >                mrad,vt,bt,zz,stval,dx,nlshell,nqntab,lqntab,jtop,
       6             :      X                ectab,
       7           0 :      <                rhochr,rhospn)
       8             :       USE m_felim
       9             :       USE m_findlim
      10             :       USE m_cnodes
      11             :       USE m_coredir
      12             :       USE m_ccsdnt
      13             : c
      14             :       IMPLICIT NONE
      15             : C     ..
      16             : C     .. Scalar Arguments ..
      17             :       INTEGER, INTENT (IN) :: mrad,jtop,nlshell
      18             :       REAL,    INTENT (IN) :: dx,stval,zz
      19             : C     ..
      20             : C     .. Array Arguments ..
      21             :       INTEGER, INTENT (IN) :: lqntab(15),nqntab(15)
      22             :       REAL,    INTENT (IN) :: bt(mrad),vt(mrad)
      23             :       REAL,    INTENT (OUT) :: rhochr(mrad),rhospn(mrad)
      24             :       REAL,    INTENT (INOUT) :: ectab(100)
      25             : C     ..
      26             : C     .. Local Scalars ..
      27             :       REAL bhf,bsum,btot,dvstep,ec,elim,qcor,spcor,tolvar,vz,
      28             :      +     xmj,ec_sv
      29             :       INTEGER i,ic,ic1,ic2,ie,iflag,ii,ilshell,imin,ir,ish,istart,iter,
      30             :      +        itermax,iv,j,jv,kap1,kap2,l,lll,muem05,n,ncor,nmatch,node,
      31             :      +        nqn,nrmax,nsh,nsol,nval,nvar,nzero,s,t
      32             :       LOGICAL ferro
      33             : C     ..
      34             : C     .. Local Arrays ..
      35           0 :       REAL bb(mrad),bhff(15),dedv(4,4),dv(4),dvde(4,4),err(4),errnew(4),
      36           0 :      +     fc(2,2,mrad),fck(2,2,mrad),gc(2,2,mrad),gck(2,2,mrad),
      37           0 :      +     piw(2,2),pow(2,2),qiw(2,2),qow(2,2),rc(mrad),rc2(mrad),
      38           0 :      +     rchr(mrad),rspn(mrad),var(4),varnew(4),vv(mrad)
      39             :       INTEGER kap(2)
      40             :       CHARACTER txtl(0:3)
      41             :       CHARACTER*3 txtk(4)
      42             : C     ..
      43             : C     .. External Functions ..
      44             :       REAL rsimp
      45             :       EXTERNAL rsimp
      46             : C     ..
      47             : C     .. External Subroutines ..
      48             :       EXTERNAL cfnorm,coreerr,corehff,nwrfst,rinvgj
      49             : C     ..
      50             : C     .. Intrinsic Functions ..
      51             :       INTRINSIC abs,exp,iabs,int,sign
      52             : C     ..
      53             : C     .. Data statements ..
      54             : 
      55             :       DATA txtl/'s','p','d','f'/
      56             :       DATA txtk/'1/2','3/2','5/2','7/2'/
      57             : !      DATA nqntab/1,2,2,3,3,3,4,4,4,4,5,5,5,6,6/
      58             : !      DATA lqntab/0,0,1,0,1,2,0,1,2,3,0,1,2,0,1/
      59             :       DATA tolvar/1.0E-12/
      60             :       DATA dvstep/0.010/
      61             :       DATA itermax/150/
      62             : C     ..
      63             : c
      64           0 :       nrmax = mrad
      65           0 :       DO n = 1,nrmax
      66           0 :          rc(n) = exp(stval+ (n-1)*dx)
      67           0 :          rc2(n) = rc(n)**2
      68             :       END DO
      69             : c
      70           0 :       DO n = 1,nrmax
      71           0 :          rhochr(n) = 0.00
      72           0 :          rhospn(n) = 0.00
      73             :       END DO
      74             : c
      75           0 :       bsum = 0.0
      76           0 :       DO n = 1,jtop
      77           0 :          vv(n) = vt(n)
      78           0 :          bb(n) = bt(n)
      79           0 :          bsum = bsum + abs(bb(n))
      80             :       END DO
      81             : C
      82             : 
      83           0 :       ncor = 0
      84           0 :       DO 10 ilshell = 1,nlshell
      85           0 :          ncor = ncor + 2* (2*lqntab(ilshell)+1)
      86           0 :    10 CONTINUE
      87           0 :       nval = int(zz) - ncor
      88             : 
      89           0 :       IF (bsum.GT.1.0E-8) THEN
      90           0 :          ferro = .true.
      91             :       ELSE
      92           0 :          ferro = .false.
      93             :          WRITE (6,FMT=
      94             :      +'(''  PARAMAGNETIC CASE'',/,
      95           0 :      +   ''  *****************'')')
      96             :       END IF
      97             : 
      98             :       WRITE (6,FMT=
      99             :      +'(/,''  ATOMIC NUMBER  : '',F6.2,/,
     100             :      +    ''  CORE ELECTRONS : '',I5,/,
     101           0 :      +    ''  VAL. ELECTRONS : '',I5)') zz,ncor,nval
     102             :       WRITE (6,FMT=
     103             :      +'(  /,
     104             :      +    ''  MESH    RC(1)  : '',F12.6,//10x,
     105             :      +    ''MUE  KAP ITER    ENERGY '',
     106           0 :      +    ''(Hart)'', " HF (KG) ")') rc(1)
     107             : 
     108             : 
     109             : C                   ---------------------------------------
     110             : C--->               INITIALIZE QUANTUM NUMBERS  NQN  AND  L
     111             : C                   ---------------------------------------
     112           0 :       btot = 0.0
     113           0 :       ic = 0
     114             : C ---> nl  - loop
     115           0 :       DO 190 ilshell = 1,nlshell
     116           0 :          nqn = nqntab(ilshell)
     117           0 :          l = lqntab(ilshell)
     118           0 :          nsh = 2* (2*l+1)
     119           0 :          bhff(ilshell) = 0.0
     120           0 :          ish = 0
     121             : C ---> \mu  - loop
     122           0 :          DO 180 muem05 = -l - 1,+l
     123           0 :             xmj = muem05 + 0.5
     124           0 :             kap1 = -l - 1
     125           0 :             kap2 = l
     126           0 :             kap(1) = kap1
     127           0 :             kap(2) = kap2
     128             : c
     129           0 :             lll = l* (l+1)
     130           0 :             IF (abs(xmj).GT.l) THEN
     131           0 :                nsol = 1
     132             :             ELSE
     133           0 :                nsol = 2
     134             :             END IF
     135             : c
     136           0 :             nvar = 2*nsol
     137             : c
     138             : c----> s - loop : solutions for each (nl,\mu)
     139             : c
     140           0 :             DO 170 s = 1,nsol
     141             : c
     142           0 :                DO i = 1,2
     143           0 :                   DO j = 1,2
     144           0 :                      DO n = 1,mrad
     145           0 :                         gc(i,j,n) = 0.0
     146           0 :                         gck(i,j,n) = 0.0
     147           0 :                         fc(i,j,n) = 0.0
     148           0 :                         fck(i,j,n) = 0.0
     149             :                      END DO
     150             :                   END DO
     151             :                END DO
     152             : c
     153           0 :                ic = ic + 1
     154           0 :                ish = ish + 1
     155           0 :                t = 3 - s
     156           0 :                ec = ectab(ic)
     157           0 :                write(*,*) ic,ec
     158             : c+gu
     159           0 :                IF (ic > 1) THEN
     160           0 :                  IF (ectab(ic-1) > 0.1) THEN
     161           0 :                    vv(:) = vv(:) + 2*ec_sv
     162           0 :                    bb(:) = bb(:) + 2*ec_sv
     163             :                  ENDIF
     164             :                ENDIF
     165           0 :                IF (ec.GT.0.1) THEN !  GOTO 170
     166           0 :                   vv(:) = vv(:) - 2*ectab(ic)
     167           0 :                   bb(:) = bb(:) - 2*ectab(ic)
     168           0 :                   ec_sv = ectab(ic)
     169           0 :                   ec = -ectab(ic)
     170             :                ELSE
     171             :                  ec_sv = 0.0
     172             :                ENDIF
     173           0 :                istart = 1
     174             : c
     175           0 :                IF (ish.GT.1) GO TO 30
     176             : c
     177             : c energy limits : elim and emax
     178             : c
     179             :                CALL felim(
     180             :      >                    mrad,lll,zz,nqn,vv,rc,
     181           0 :      <                    elim)
     182           0 :                write(*,*) 'elim',elim
     183           0 :    20          IF (ec.LE.elim) ec = elim*0.70
     184             : c
     185             : c turning point and physical infinity : nmatch and nzero
     186             : c
     187             :                CALL findlim(
     188             :      >                      mrad,lll,ec,vv,rc,
     189           0 :      <                      nmatch,nzero)
     190           0 :                write(*,*) 'nmatch,nzero',nmatch,nzero
     191             :    30          CONTINUE
     192             : c
     193             : c nodes for large component : if IFLAG=0 number of nodes is OK
     194             : c
     195           0 :                iflag = 0
     196             :                CALL cnodes(mrad,iflag,s,ec,l,xmj,nqn,vv,bb,rc,dx,
     197           0 :      +                     nmatch,nzero,gc,fc,pow,qow,piw,qiw,node)
     198           0 :                IF (iflag.NE.0) GO TO 20
     199             : c
     200             : c setup of Newton-Raphson algorithm
     201             : c
     202             :                CALL nwrfst(mrad,nsol,s,t,nmatch,nzero,ferro,ec,rc,
     203           0 :      +                     pow,piw,gc,err,var,dv,varnew,errnew)
     204             : c
     205           0 :                CALL coreerr(err,var,s,nsol,pow,qow,piw,qiw)
     206           0 :                DO 40 iv = 1,nvar
     207           0 :                   dv(iv) = var(iv)
     208           0 :    40          CONTINUE
     209             : c
     210           0 :                iter = 0
     211             : C---> iterations start
     212           0 :    50          iter = iter + 1
     213             : C                         ----------------------------------
     214             : C                         CHECK WHETHER NUMBER OF NODES O.K.
     215             : C                         ----------------------------------
     216           0 :                IF (iter.GT.1) THEN
     217           0 :                   node = 0
     218           0 :                   DO 60 n = 2, (nmatch-1)
     219           0 :                      IF (gc(s,s,n)*gc(s,s,n-1).LT.0.0) node = node + 1
     220           0 :    60             CONTINUE
     221             : c
     222           0 :                   IF (node.NE. (nqn-l-1)) THEN
     223           0 :                      IF (node.GT. (nqn-l-1)) THEN
     224           0 :                         ec = 1.2*ec
     225           0 :                         write (*,*) 'up',node,ec
     226             :                      ELSE
     227           0 :                         ec = 0.8*ec
     228           0 :                         write (*,*) 'do',node,ec
     229             :                      END IF
     230           0 :                      istart = istart + 1
     231           0 :                      IF (istart.LT.20) GO TO 20
     232             :                   END IF
     233             :                END IF
     234             : c                    - atomic energy search -
     235           0 :                DO 90 iv = 2,nvar
     236           0 :                   DO 70 jv = 1,nvar
     237           0 :                      varnew(jv) = var(jv)
     238           0 :    70             CONTINUE
     239           0 :                   varnew(iv) = var(iv) + dv(iv)*dvstep
     240             : c
     241           0 :                   IF (abs(var(iv)).EQ.0.00) THEN
     242           0 :                      IF (ferro) THEN
     243           0 :                         WRITE (6,FMT='(A,I3,A)') ' VAR ',iv,
     244           0 :      +                    ' = 0 ??????!!!!!'
     245             :                      END IF
     246             :                   ELSE
     247           0 :                      IF (abs(dv(iv)/var(iv)).LT.
     248             :      +                   tolvar) varnew(iv) = var(iv)*
     249           0 :      +                   (1.00+sign(dvstep*tolvar,dv(iv)))
     250             :                   END IF
     251             : c
     252           0 :                   CALL coreerr(errnew,varnew,s,nsol,pow,qow,piw,qiw)
     253             : c
     254           0 :                   DO 80 ie = 1,nvar
     255           0 :                      IF ((errnew(ie)-err(ie)).EQ.0.00) THEN
     256           0 :                         dedv(ie,iv) = 0.00
     257           0 :                         IF ((ie.EQ.iv) .AND. .NOT.ferro) THEN
     258           0 :                            dedv(ie,iv) = 1.00
     259             :                         END IF
     260             :                      ELSE
     261             :                         dedv(ie,iv) = (errnew(ie)-err(ie))/
     262           0 :      +                                (varnew(iv)-var(iv))
     263             :                      END IF
     264           0 :    80             CONTINUE
     265           0 :    90          CONTINUE
     266             : c
     267           0 :                DO 100 jv = 1,nvar
     268           0 :                   varnew(jv) = var(jv)
     269           0 :   100          CONTINUE
     270           0 :                varnew(1) = var(1) + dv(1)*dvstep
     271           0 :                IF (abs(dv(1)/var(1)).LT.tolvar) varnew(1) = var(1)*
     272           0 :      +             (1.00+sign(dvstep*tolvar,dv(1)))
     273             : c
     274             :                CALL coredir(mrad,varnew(1),l,xmj,1,vv,bb,rc,dx,
     275           0 :      +                      nmatch,nzero,gc,fc,pow,qow,piw,qiw)
     276             :                CALL coredir(mrad,varnew(1),l,xmj,2,vv,bb,rc,dx,
     277           0 :      +                      nmatch,nzero,gc,fc,pow,qow,piw,qiw)
     278             : c
     279           0 :                CALL coreerr(errnew,varnew,s,nsol,pow,qow,piw,qiw)
     280             : c
     281           0 :                DO 110 ie = 1,nvar
     282           0 :                   dedv(ie,1) = (errnew(ie)-err(ie))/ (varnew(1)-var(1))
     283           0 :   110          CONTINUE
     284             : c
     285           0 :                CALL rinvgj(dvde,dedv,4,nvar)
     286             : c
     287           0 :                DO 130 iv = 1,nvar
     288           0 :                   dv(iv) = 0.00
     289           0 :                   DO 120 ie = 1,nvar
     290           0 :                      dv(iv) = dv(iv) + dvde(iv,ie)*err(ie)
     291           0 :   120             CONTINUE
     292           0 :                   var(iv) = var(iv) - dv(iv)
     293           0 :   130          CONTINUE
     294             : c
     295           0 :                IF (var(1).GT.0.00) THEN
     296           0 :                   WRITE (6,FMT=*) ' WARNING FROM <CORE> E=',var(1)
     297           0 :                   var(1) = -0.20
     298             :                END IF
     299             : c
     300             :                CALL coredir(mrad,var(1),l,xmj,1,vv,bb,rc,dx,
     301           0 :      +                      nmatch,nzero,gc,fc,pow,qow,piw,qiw)
     302             :                CALL coredir(mrad,var(1),l,xmj,2,vv,bb,rc,dx,
     303           0 :      +                      nmatch,nzero,gc,fc,pow,qow,piw,qiw)
     304             : c
     305           0 :                CALL coreerr(err,var,s,nsol,pow,qow,piw,qiw)
     306             : c
     307           0 :                ec = var(1)
     308             : c
     309             : C---> check relative change in parameters
     310             : C---> parameters 3 and 4 = 0 for paramagnetic case !
     311             : c
     312           0 :                IF (iter.LT.itermax) THEN
     313           0 :                   IF ((nsol.EQ.2) .AND. .NOT.ferro) THEN
     314           0 :                      DO iv = 3,4
     315           0 :                         err(iv) = 0.00
     316           0 :                         errnew(iv) = 0.00
     317           0 :                         var(iv) = 0.00
     318           0 :                         varnew(iv) = 0.00
     319           0 :                         dv(iv) = 0.00
     320             :                      END DO
     321             :                   END IF
     322           0 :                   DO 140 iv = 1,nvar
     323           0 :                      IF (abs(var(iv)).EQ.0.00) THEN
     324           0 :                         IF (ferro) THEN
     325           0 :                            WRITE (6,FMT='(A,I3,A)') ' VAR ',iv,
     326           0 :      +                       ' = 0 ??????!!!!!'
     327             :                         END IF
     328             :                      ELSE
     329           0 :                         IF (abs(dv(iv)/var(iv)).GT.tolvar) GO TO 50
     330             :                      END IF
     331           0 :   140             CONTINUE
     332             :                ELSE
     333             :                   WRITE (6,FMT=
     334             :      +'('' ITERATION NOT CONVERGED AFTER'',I3,'' STEPS !'',/,
     335             :      + '' PARAMETERS:'',4E18.10,/,'' LAST CORR.:'',4E18.10,/,
     336           0 :      +  '' LAST ERROR:'',4E18.10)') itermax, (var(iv),iv=1,4),
     337           0 :      +              (dv(iv),iv=1,4), (err(ie),ie=1,4)
     338             :                END IF
     339             : c
     340             : C---> end of iterations
     341             : c
     342             :                CALL cfnorm(mrad,s,t,nsol,nmatch,jtop,var,gc,fc,rc,rc2,
     343           0 :      +                     dx,gck,fck)
     344             : C                --------------------------
     345             : C--->             CALCULATE  CHARGE DENSITY
     346             : C                --------------------------
     347             :                CALL ccsdnt(mrad,s,jtop,nsol,l,xmj,kap1,kap2,gck,fck,
     348           0 :      +                     rc2,rchr,rspn)
     349           0 :                DO ir = 1,jtop
     350           0 :                   rhochr(ir) = rhochr(ir) + rchr(ir)
     351           0 :                   rhospn(ir) = rhospn(ir) + rspn(ir)
     352             :                END DO
     353             : c                -----------------
     354             : c--->             HYPERFINE FIELD
     355             : c                -----------------
     356             :                CALL corehff(mrad,kap1,kap2,xmj,s,nsol,bhf,gck,fck,
     357           0 :      +                      rc,dx,jtop)
     358           0 :                bhff(ilshell) = bhff(ilshell) + bhf
     359             : c                 - final info. -
     360           0 :                ectab(ic) = ec + 2*ec_sv
     361             : c
     362           0 :                IF (ish.LT.nsh) THEN
     363           0 :                   WRITE (6,FMT=8000) nqn,txtl(l),txtk(iabs(kap(s))),
     364           0 :      +              (2*muem05+1),kap(s),iter,ec/2.
     365             :                ELSE
     366           0 :                   WRITE (6,FMT=8000) nqn,txtl(l),txtk(iabs(kap(s))),
     367           0 :      +              (2*muem05+1),kap(s),iter,ec/2.
     368             : 
     369             : C              ----------------------------
     370             : C--->          CHECK CONSISTENCY OF RESULTS
     371             : C              ----------------------------
     372           0 :                   IF (l.NE.0) THEN
     373           0 :                      ic1 = ic - nsh + 1
     374           0 :                      ic2 = ic
     375           0 :                      IF (ectab(ic2).GE.ectab(ic1)) THEN
     376             :                         imin = ic1
     377             :                         vz = +1.00
     378             :                      ELSE
     379           0 :                         imin = ic2
     380           0 :                         vz = -1.00
     381             :                      END IF
     382           0 :                      iflag = 0
     383           0 :                      ii = 1
     384           0 :                      DO 150 i = ic1 + 1,ic2,2
     385           0 :                         IF (vz* (ectab(i)-ectab(i-ii)).LT.0.0) iflag = 1
     386           0 :                         ii = 2
     387           0 :   150                CONTINUE
     388           0 :                      IF (ectab(ic1+2).GT.ectab(imin)) iflag = 1
     389           0 :                      DO 160 i = ic1 + 4,ic2 - 1,2
     390           0 :                         IF (ectab(i).GT.ectab(imin)) iflag = 1
     391           0 :                         IF (vz* (ectab(i)-ectab(i-ii)).GT.0.0) iflag = 1
     392           0 :   160                CONTINUE
     393             : c
     394           0 :                      IF (ferro .AND. (iflag.EQ.1)) WRITE (6,FMT=
     395             :      +'('' >>> CHECK DATA '',                               '' E(KAP,MJ)
     396             :      + SHOULD BE MONOTONOUS AND '',                         '' E(+L,MJ)
     397           0 :      +< E(-L-1,MJ) '',//)')
     398             :                   END IF
     399             :                END IF
     400             : c
     401           0 :   170       CONTINUE
     402             : c----> s - loop end
     403           0 :   180    CONTINUE
     404             : C----> \mu  - loop end
     405           0 :          WRITE (6,FMT='(I4,A1,20X,F16.3)') nqn,txtl(l),bhff(ilshell)
     406           0 :          btot = bhff(ilshell) + btot
     407           0 :   190 CONTINUE
     408             : C----> nl  - loop end
     409           0 :       WRITE (6,FMT='(''HFTOT:'',20X,F16.3)') btot
     410             : c
     411           0 :       qcor = rsimp(mrad,rhochr,rc,jtop,dx)
     412           0 :       WRITE (6,FMT=8020) 'charge',qcor
     413           0 :       spcor = rsimp(mrad,rhospn,rc,jtop,dx)
     414           0 :       WRITE (6,FMT=8020) 'spin',spcor
     415             : c
     416             : C----> that's all
     417             : c
     418             : c
     419             :  8000 FORMAT (I4,A1,A3,I3,'/2',2I4,1X,F14.7,F16.3,:,F16.3,/)
     420             :  8010 FORMAT (/,'  NQN=',I2,'  L=',I2,'  KAP=',I2,'  MJ=',I2,
     421             :      +       '/2    IC=',I3,'  ISH=',I2,/,' E(',A5,')=',F15.5,/,
     422             :      +       ' NMATCH  =',I5,'    R=',F10.5,/,' NZERO   =',I5,'    R=',
     423             :      +       F10.5,/,' NODES   =',I5,'  RAT=',E11.4)
     424             :  8020 FORMAT (' integrated ',A,':',F12.8)
     425             : c
     426           0 :       RETURN
     427             :       END SUBROUTINE core
     428             :       END MODULE m_core

Generated by: LCOV version 1.13