umtx.f90 3.85 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
MODULE m_umtx
  USE m_juDFT
  !*********************************************************************
  !* The calculation of the "U"-contribution to Hartree-Fock matrix.   *
  !*********************************************************************
CONTAINS
  SUBROUTINE umtx(&
       lmaxb,ntype,n_u,lda_u,f0,f2,f4,f6,&
       u)

Daniel Wortmann's avatar
Daniel Wortmann committed
11
    USE m_constants
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
    USE m_sgaunt
    IMPLICIT NONE

    INTEGER, PARAMETER   :: lmaxw=3,lmmaxw1=(2*lmaxw+2)**2
    INTEGER, INTENT (IN) :: n_u,lmaxb,ntype
    INTEGER, INTENT (IN) :: lda_u(ntype)
    REAL,    INTENT (IN) :: f0(n_u),f2(n_u),f4(n_u),f6(n_u)
    REAL,    INTENT (OUT) :: u(-lmaxb:lmaxb,-lmaxb:lmaxb,&
         -lmaxb:lmaxb,-lmaxb:lmaxb,n_u)

    INTEGER i,j,k,l,m,mk,nfk,n,itype
    INTEGER m1,m2,m3,m4,lm1,lm2,lm3,lm4,kf,l_l(n_u)
    REAL    uk,uq,avu,avj,cgk1,cgk2,tol
    REAL    fk(lmaxb+1,n_u)
    REAL,   ALLOCATABLE :: c(:,:,:)
    !
    tol = 1.0e-14
    !
    ! transformation to Hr-units:
    !
    n = 0
    DO itype = 1,ntype
       IF (lda_u(itype).GE.0) THEN
          n = n + 1
          l_l(n) = lda_u(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
37 38 39
          fk(1,n) = f0(n) / hartree_to_ev_const
          fk(2,n) = f2(n) / hartree_to_ev_const
          fk(3,n) = f4(n) / hartree_to_ev_const
40
          IF (l_l(n).EQ.3) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
41
             fk(4,n) = f6(n) / hartree_to_ev_const
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
          ELSEIF (l_l(n).GT.3) THEN
             CALL juDFT_error("LDA+U for p, d or f-states!", calledby="umtx")
          ENDIF
       ENDIF
    ENDDO
    !
    ! evaluate Gaunt parameter
    !
    ALLOCATE( c(0:2*lmaxw+1,lmmaxw1,lmmaxw1) )
    DO k = 1,lmmaxw1
       DO j = 1,lmmaxw1
          DO i = 0,2*lmaxw+1
             c(i,j,k) = 0.0
          ENDDO
       ENDDO
    ENDDO
    CALL sgaunt(lmaxw,lmmaxw1,lmaxb,&
         c)
    !
    ! lda_u(n) is here only the 'l' for atom 'n'
    !
    DO  n = 1,n_u                      !!! over d-atoms
       l = l_l(n)
       kf = 2*l
       DO m1 = -l,l
          lm1 = l*(l+1)+m1+1
          DO m2 = -l,l
             lm2 = l*(l+1)+m2+1
             DO m3 = -l,l
                lm3 = l*(l+1)+m3+1
                DO m4 = -l,l
                   lm4 = l*(l+1)+m4+1
                   uk = 0.0e0
                   DO k=0,kf,2
                      uq = 0.e0
                      DO mk=-k,k
                         IF (mk.NE.m1-m3)  CYCLE
                         cgk1 = c(k/2,lm1,lm3)
                         IF (ABS(cgk1).LT.tol) CYCLE
                         IF (mk.NE.m4-m2)  CYCLE
                         cgk2 = c(k/2,lm4,lm2)
                         IF (ABS(cgk2).LT.tol) CYCLE
                         uq = uq+cgk1*cgk2
                      ENDDO                   ! mk
                      IF (ABS(uq).LT.tol) CYCLE
                      nfk=k/2+1
                      uk=uk+uq*fk(nfk,n)*4*pi_const/(2*k+1)
                   ENDDO                     ! k
                   u(m1,m2,m3,m4,n)=uk
                ENDDO                       ! m4 etc.
             ENDDO
          ENDDO
       ENDDO
       avu=0.e0
       avj=0.e0

       DO i = -l,l
          DO j = -l,l
             avu = avu+u(i,j,i,j,n)
             avj = avj+(u(i,j,i,j,n)-u(i,j,j,i,n))
          ENDDO
       ENDDO
       avu = avu/(2*l+1)/(2*l+1)
       avj = avj/(2*l+1)/(2*l)
       avj = avu-avJ
       !        WRITE (6,*) 'U-matr:'
       !        IF (l.eq.2) WRITE (6,111) ((u(i,j,i,j,n),i=-l,l),j=-l,l)
       !        IF (l.eq.3) WRITE (6,211) ((u(i,j,i,j,n),i=-l,l),j=-l,l)
       !        WRITE (6,*) 'J-matr:'
       !        IF (l.eq.2) WRITE (6,111) ((u(i,j,j,i,n),i=-l,l),j=-l,l)
       !        IF (l.eq.3) WRITE (6,211) ((u(i,j,j,i,n),i=-l,l),j=-l,l)
Daniel Wortmann's avatar
Daniel Wortmann committed
113 114
       !         PRINT*,'U-av:',avu*hartree_to_ev_const
       !         PRINT*,'J-av:',avj*hartree_to_ev_const
115 116 117 118 119 120 121 122 123
111    FORMAT (5f8.4)
211    FORMAT (7f8.4)
112    FORMAT (10e20.10)
       !c         WRITE (9,112) ((((u(i,j,k,m,n),i=-l,l),j=-l,l),k=-l,l),m=-l,l)
    ENDDO
    DEALLOCATE (c)

  END SUBROUTINE umtx
END MODULE m_umtx