LCOV - code coverage report
Current view: top level - core - setcor.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 101 156 64.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_setcor
       2             :   USE m_juDFT
       3             : CONTAINS
       4         660 :   SUBROUTINE setcor(itype,jspins,atoms,input,bmu,nst,kappa,nprnc,occ)
       5             :     !
       6             :     !     *****************************************************
       7             :     !     sets the values of kappa and occupation numbers of
       8             :     !     the neutral atoms.
       9             :     !         following code by m. weinert  february 1982
      10             :     !     *****************************************************
      11             : 
      12             :     USE m_types
      13             :     IMPLICIT NONE
      14             :     TYPE(t_atoms),INTENT(IN)   :: atoms
      15             :     TYPE(t_input),INTENT(IN)   :: input
      16             :     !
      17             :     !     .. Scalar Arguments ..
      18             :     INTEGER,INTENT (IN) :: itype,jspins
      19             :     REAL,INTENT (INOUT) :: bmu
      20             :     !     ..
      21             :     INTEGER,INTENT (OUT) :: nst
      22             :     INTEGER,INTENT (OUT) :: kappa(:),nprnc(:)
      23             :     REAL,INTENT (OUT)    :: occ(:,:)
      24             :     !     ..
      25             :     !     .. Local Scalars ..
      26             :     INTEGER iz,jz,jz0,k,n,m,i,tempInt
      27             :     INTEGER k_h(2),n_h(2)
      28             :     REAL fj,l,bmu_l,o_h(2), fac(2),tempReal
      29             :     LOGICAL l_clf
      30             :     CHARACTER(len=13) :: fname
      31             :     !     ..
      32             : 
      33         660 :     l_clf = .FALSE.  
      34         660 :     WRITE(fname,"('corelevels.',i2.2)") NINT(atoms%zatom(itype))
      35         660 :     INQUIRE (file=fname, exist=l_clf)
      36             : 
      37         660 :     IF (l_clf.AND..NOT.input%l_inpXML) THEN
      38           0 :        OPEN (61,file=fname,form='formatted')
      39           0 :        READ (61,'(i3)') nst
      40           0 :        IF (bmu.LT.0.001) bmu = 999.
      41           0 :        IF (nst > size(kappa))  CALL juDFT_error("corelevels: nst > nstd" ,calledby ="setcor")
      42           0 :        DO n = 1, nst
      43           0 :           fac(1) = 1.0 ; fac(2) = 1.0
      44           0 :           READ (61,'(4i3)') nprnc(n),kappa(n),n_h(1),n_h(2)
      45           0 :           IF ( n_h(1) < 0 )  fac(1) = -0.5
      46           0 :           IF ( n_h(2) < 0 )  fac(2) = -0.5
      47           0 :           IF (jspins.EQ.1) THEN
      48           0 :              occ(n,1) = fac(1) * n_h(1) + fac(2) * n_h(2) 
      49             :           ELSE
      50           0 :              occ(n,1) = fac(1) * n_h(1) ; occ(n,2) = fac(2) * n_h(2)
      51             :           ENDIF
      52             :           !         write(*,*) nprnc(n),kappa(n),occ(n,1), occ(n,2)
      53             :        ENDDO
      54           0 :        CLOSE (61)
      55           0 :        RETURN
      56             :     ENDIF
      57             : 
      58         660 :     jz0 = atoms%zatom(itype) + 0.01e0
      59         660 :     jz = jz0
      60         660 :     k = 0
      61        5218 :     DO n = 1,7 !maximal main quantuum number =7
      62        2609 :        IF (jz.LE.0) EXIT
      63             :        !--->    s states
      64        2609 :        k = k + 1
      65        2609 :        nprnc(k) = n
      66        2609 :        kappa(k) = -1
      67        2609 :        occ(k,1) = 2
      68        2609 :        jz = jz - 2
      69        2609 :        IF (jz.LE.0) EXIT
      70             :        !--->    p states
      71        2032 :        IF (n.EQ.1) CYCLE
      72        1372 :        k = k + 1
      73        1372 :        nprnc(k) = n
      74        1372 :        kappa(k) = 1
      75        1372 :        occ(k,1) = 2
      76        1372 :        jz = jz - 2
      77        1372 :        IF (jz.LE.0) EXIT
      78        1324 :        k = k + 1
      79        1324 :        nprnc(k) = n
      80        1324 :        kappa(k) = -2
      81        1324 :        occ(k,1) = 4
      82        1324 :        jz = jz - 4
      83        1324 :        IF (jz.LE.0) EXIT
      84             :        !--->    d functions
      85        1289 :        iz = 0
      86        1289 :        IF (n.EQ.3 .AND. jz0.GT.20) iz = MIN(jz0-20,4)
      87        1289 :        IF (n.EQ.4 .AND. jz0.GT.38) iz = MIN(jz0-38,4)
      88        1289 :        IF (n.EQ.4 .AND. jz0.EQ.41) iz = 4
      89        1289 :        IF (n.EQ.5 .AND. jz0.GT.70) iz = MIN(jz0-70,4)
      90        1289 :        IF (n.EQ.5 .AND. (jz0.EQ.57.OR.jz0.EQ.64)) iz = 1
      91        1289 :        IF (n.EQ.6 .AND. jz0.GT.88) iz = 1
      92        1289 :        IF (n.EQ.6 .AND. jz0.EQ.90) iz = 2
      93        1289 :        IF (iz.NE.0) THEN 
      94         645 :           k = k + 1
      95         645 :           nprnc(k) = n
      96         645 :           kappa(k) = 2
      97         645 :           occ(k,1) = iz
      98         645 :           jz = jz - iz
      99         645 :           IF ((n==6).AND.(iz.GE.4)) CYCLE
     100         645 :           IF (iz.GE.4 .AND. .NOT.(n.EQ.4 .AND. jz0.EQ.41) .AND. .NOT. (n.EQ.5 .AND. jz0.EQ.74)) THEN
     101         635 :              iz = 1
     102         635 :              IF (n.EQ.3 .AND. jz0.GT.25) iz = MIN(jz0-24,6)
     103         635 :              IF (n.EQ.3 .AND. jz0.EQ.29) iz = 6
     104         635 :              IF (n.EQ.4 .AND. jz0.GT.43) iz = jz0 - 41
     105         635 :              IF (n.EQ.4 .AND. jz0.GT.45) iz = 6
     106         635 :              IF (n.EQ.5 .AND. jz0.GT.75) iz = jz0 - 74
     107         635 :              IF (n.EQ.5 .AND. jz0.GT.77) iz = 6
     108         635 :              k = k + 1
     109         635 :              nprnc(k) = n
     110         635 :              kappa(k) = -3
     111         635 :              occ(k,1) = iz
     112         635 :              jz = jz - iz
     113             :              
     114             :           ENDIF
     115             :        ENDIF
     116             :        !--->    f states
     117        1289 :        IF (n==4) THEN
     118             :           !+gu  IF (jz0.LE.57) GO TO 50
     119          10 :           IF (jz0.LE.62) CYCLE
     120          10 :           k = k + 1
     121          10 :           nprnc(k) = n
     122          10 :           kappa(k) = 3
     123          10 :           iz = 6
     124             :           IF (jz0.LT.62) THEN
     125             :              iz = jz0 - 56
     126             :              occ(k,1) = iz
     127             :              jz = jz - iz
     128             :              CYCLE
     129             :           ENDIF
     130          10 :           occ(k,1) = iz
     131          10 :           jz = jz - iz
     132          10 :           iz = 8
     133          10 :           k = k + 1
     134          10 :           nprnc(k) = n
     135          10 :           kappa(k) = -4
     136          10 :           IF (jz0.LT.70) THEN
     137           0 :              iz = jz0 - 62
     138           0 :              IF (jz0.EQ.64) iz = 1
     139           0 :              occ(k,1) = iz
     140           0 :              jz = jz - iz
     141           0 :              CYCLE
     142             :           ENDIF
     143          10 :           occ(k,1) = iz
     144          10 :           jz = jz - iz
     145             :        ENDIF
     146        1289 :        IF (n.NE.5) CYCLE
     147          10 :        IF (jz0.LE.90) CYCLE
     148           0 :        k = k + 1
     149           0 :        nprnc(k) = n
     150           0 :        kappa(k) = 3
     151           0 :        iz = jz0 - 89
     152           0 :        occ(k,1) = iz
     153        2609 :        jz = jz - iz
     154             :     ENDDO
     155         660 :     nst = k
     156         660 :     IF (k.GE.1) occ(k,1) = occ(k,1) + jz
     157             :     !
     158             :     ! add magnetic moments
     159             :     !
     160         660 :     IF (jspins.EQ.2) THEN
     161         515 :        bmu_l = bmu
     162        5665 :        DO k = 1,nst
     163        5150 :           occ(k,jspins) = occ(k,1)/2.0 
     164        5665 :           occ(k,1) = occ(k,jspins)  
     165             :        ENDDO
     166        1569 :        kloop:DO k = nst,1,-1
     167        1042 :           fj = iabs(kappa(k)) - 0.5e0
     168        1042 :           l = fj + 0.5e0*isign(1,kappa(k)) + 0.01e0
     169             :           ! polarize (d,f) only
     170        1042 :           IF (l.GT.1.99) THEN
     171         527 :              IF (2*occ(k,1).GE.ABS(bmu_l)) THEN
     172         515 :                 occ(k,1) = occ(k,1) + bmu_l/2.
     173         515 :                 occ(k,jspins) = occ(k,jspins) - bmu_l/2.
     174         515 :                 EXIT kloop
     175             :              ELSE
     176          12 :                 IF (bmu_l.GT.0) THEN
     177          12 :                    occ(k,1) = 2.0*occ(k,1)
     178          12 :                    occ(k,jspins) = 0.0
     179          12 :                    bmu_l = bmu_l - occ(k,1)
     180             :                 ELSE
     181           0 :                    occ(k,jspins) = 2.0*occ(k,jspins)
     182           0 :                    occ(k,1) = 0.0
     183           0 :                    bmu_l = bmu_l + occ(k,jspins)
     184             :                 ENDIF
     185             :              ENDIF
     186             :           ENDIF
     187             :        ENDDO kloop
     188             :     ENDIF
     189             :     !
     190         660 :     IF (atoms%zatom(itype).EQ.65) THEN
     191           0 :        k_h(1) = kappa(15) ; n_h(1) = nprnc(15) ; o_h(1) = occ(15,1)
     192           0 :        k_h(2) = kappa(16) ; n_h(2) = nprnc(16) ; o_h(2) = occ(16,1)
     193           0 :        kappa(15)= kappa(17) ; nprnc(15)=nprnc(17) ; occ(15,1)=occ(17,1)
     194           0 :        kappa(16)= kappa(18) ; nprnc(16)=nprnc(18) ; occ(16,1)=occ(18,1)
     195           0 :        kappa(17)= kappa(19) ; nprnc(17)=nprnc(19) ; occ(17,1)=occ(19,1)
     196           0 :        kappa(18) = k_h(1) ; nprnc(18) =  n_h(1)  ; occ(18,1)= o_h(1)
     197           0 :        kappa(19) = k_h(2) ; nprnc(19) =  n_h(2)  ; occ(19,1)= o_h(2)
     198             : 
     199           0 :        IF (jspins.EQ.2) THEN
     200           0 :           o_h(1) = occ(15,jspins) ; o_h(2) = occ(16,jspins)
     201           0 :           occ(15,jspins) = occ(17,jspins) 
     202           0 :           occ(16,jspins) = occ(18,jspins) 
     203           0 :           occ(17,jspins) = occ(19,jspins) 
     204           0 :           occ(18,jspins) = o_h(1) 
     205           0 :           occ(19,jspins) = o_h(2) 
     206             :        ENDIF
     207             :     ENDIF
     208             : 
     209             :     ! modify default electron configuration according to explicitely provided setting in inp.xml
     210         660 :     IF(input%l_inpXML) THEN
     211         470 :        nst = max(nst,atoms%numStatesProvided(itype))
     212         470 :        IF (atoms%numStatesProvided(itype).NE.0) THEN
     213           0 :           IF (bmu.LT.0.001) bmu = 999.0
     214             :        END IF
     215         470 :        DO n = 1, atoms%numStatesProvided(itype)
     216           0 :           IF((nprnc(n).NE.atoms%coreStateNprnc(n,itype)).OR.(kappa(n).NE.atoms%coreStateKappa(n,itype))) THEN
     217           0 :              m = 0
     218           0 :              DO m = n, nst
     219           0 :                 IF((nprnc(m).EQ.atoms%coreStateNprnc(n,itype)).AND.(kappa(m).EQ.atoms%coreStateKappa(n,itype))) THEN
     220             :                    EXIT
     221             :                 END IF
     222             :              END DO
     223           0 :              DO i = m-1, n, -1
     224           0 :                 nprnc(i+1) = nprnc(i)
     225           0 :                 kappa(i+1) = kappa(i)
     226           0 :                 occ(i+1,:) = occ(i,:)
     227             :              END DO
     228             :           END IF
     229           0 :           nprnc(n) = atoms%coreStateNprnc(n,itype)
     230           0 :           kappa(n) = atoms%coreStateKappa(n,itype)
     231         470 :           IF (jspins.EQ.1) THEN
     232           0 :              occ(n,1) = atoms%coreStateOccs(n,1,itype) + atoms%coreStateOccs(n,2,itype)
     233             :           ELSE
     234           0 :              occ(n,1) = atoms%coreStateOccs(n,1,itype)
     235           0 :              occ(n,2) = atoms%coreStateOccs(n,2,itype)
     236             :           END IF
     237             :        END DO
     238             :     ELSE
     239         190 :        IF (atoms%zatom(itype)>92.01e0)  CALL juDFT_error(" z > 92",calledby ="setcor")
     240             : 
     241             :     END IF
     242             : 
     243             :   END SUBROUTINE setcor
     244             :       END MODULE m_setcor

Generated by: LCOV version 1.13