etabinit.F90 4.15 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
MODULE m_etabinit
  USE m_juDFT
  !     *******************************************************
  !     *****   set up etab via old core program          *****
  !     *******************************************************
  !     modified to run with core-levels as provided by setcor
  !     ntab & ltab transport this info to core.F        gb`02
  !------------------------------------------------------------
CONTAINS
  SUBROUTINE etabinit(atoms,DIMENSION, vr,&
       etab,ntab,ltab,nkmust)

    USE m_constants, ONLY : c_light
    USE m_setcor
    USE m_differ
    USE m_types
    IMPLICIT NONE
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_atoms),INTENT(IN)       :: atoms
    !
    !     .. Scalar Arguments ..
    !     ..
    !     .. Array Arguments ..
    REAL   , INTENT (IN) :: vr(atoms%jmtd,atoms%ntypd) 
    REAL   , INTENT (OUT):: etab(100,atoms%ntypd)
    INTEGER, INTENT (OUT):: ntab(100,atoms%ntype),ltab(100,atoms%ntype)
    INTEGER, INTENT (OUT):: nkmust(atoms%ntypd)
    !     ..
    !     .. Local Scalars ..
    REAL  c,d,dxx,e,fj,fl,fn,rn,rnot,t2 ,z,t1,rr,weight
    REAL  bmu
    INTEGER i,ic,iksh,ilshell,j,jatom,korb,l, nst,ncmsh ,nshell,ipos,ierr
    !     ..
    !     .. Local Arrays ..
    INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
    REAL eig(DIMENSION%nstd),occ(DIMENSION%nstd,1),vrd(DIMENSION%msh),a(DIMENSION%msh),b(DIMENSION%msh)
    !     ..
    !
    c = c_light(1.0)
    !
    WRITE (6,FMT=8020)
    !
    ncmsh = DIMENSION%msh
    !     ---> set up densities
    DO  jatom = 1,atoms%ntype
       z = atoms%zatom(jatom)
       rn = atoms%rmt(jatom)
       dxx = atoms%dx(jatom)
       bmu = 0.0
       CALL setcor(jatom,1,atoms,bmu,nst,kappa,nprnc,occ)
       rnot = atoms%rmsh(1,jatom)
       d = EXP(atoms%dx(jatom))
       rn = rnot* (d** (ncmsh-1))
       WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
       WRITE (16,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
       DO j = 1,atoms%jri(jatom)
          vrd(j) = vr(j,jatom)
       ENDDO
#ifdef CPP_CORE
       !--->    linear extension of the potential with slope t1 / a.u.
       t1=0.125
       t2  = vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1
       rr = atoms%rmt(jatom)
       d = EXP(atoms%dx(jatom))
#else
       t2 = vrd(atoms%jri(jatom))/ (atoms%jri(jatom)-DIMENSION%msh)
#endif
       IF (atoms%jri(jatom).LT.DIMENSION%msh) THEN
          DO i = atoms%jri(jatom) + 1,DIMENSION%msh
#ifdef CPP_CORE
             rr = d*rr
             vrd(i) = rr*( t2 + rr*t1 )
#else
             vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
#endif
          ENDDO
       END IF

       nst = atoms%ncst(jatom)
       DO  korb = 1,nst
          fn = nprnc(korb)
          fj = iabs(kappa(korb)) - .5e0
          weight = 2*fj + 1.e0
          fl = fj + (.5e0)*isign(1,kappa(korb))
          e = -2* (z/ (fn+fl))**2
          CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,DIMENSION%msh,vrd,&
               e, a,b,ierr)
          IF (ierr/=0)  CALL juDFT_error("error in core-levels",calledby="etabinit")
          WRITE (6,FMT=8010) fn,fl,fj,e,weight
          WRITE (16,FMT=8010) fn,fl,fj,e,weight
          eig(korb) = e
       ENDDO
       ic = 0
       DO korb = 1,nst
          fn = nprnc(korb)
          fj = iabs(kappa(korb)) - .5e0
          weight = 2*fj + 1.e0
          fl = fj + (.5e0)*isign(1,kappa(korb))
          DO i = 1, INT(weight)
             ic = ic + 1
             IF (kappa(korb).GT.0) THEN
                ipos = ic + 1 + i 
             ELSEIF (kappa(korb).LT.-1) THEN
                ipos = ic - 2*(iabs(kappa(korb))-1) + MAX(i-2,0)
             ELSE
                ipos = ic
             ENDIF
             etab(ipos,jatom) = eig(korb)
             ntab(ipos,jatom) = NINT(fn)
             ltab(ipos,jatom) = NINT(fl)
          ENDDO
       ENDDO
       nkmust(jatom) = ic

       DO i=1,nkmust(jatom)
          WRITE(6,'(f12.6,2i3)') etab(i,jatom),ntab(i,jatom), ltab(i,jatom)
       ENDDO

    ENDDO
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f8.6,5x,&
                'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x, 'weight')
8010 FORMAT (12x,2f5.0,f6.1,f10.4,f10.0)
8020 FORMAT (/,/,12x,'core e.v. initialization')

  END SUBROUTINE etabinit
END MODULE m_etabinit