umtx.f90 4.01 KB
Newer Older
1 2 3 4
MODULE m_umtx
  USE m_juDFT
  !*********************************************************************
  !* The calculation of the "U"-contribution to Hartree-Fock matrix.   *
Daniel Wortmann's avatar
Daniel Wortmann committed
5 6
  !*-------------------------------------------------------------------*
  !* Extension to multiple U per atom type by G.M. 2017                *
7 8
  !*********************************************************************
CONTAINS
9
  SUBROUTINE umtx(atoms,f0,f2,f4,f6,&
Daniel Wortmann's avatar
Daniel Wortmann committed
10
                  u)
11

Daniel Wortmann's avatar
Daniel Wortmann committed
12
    USE m_constants
13
    USE m_sgaunt
Daniel Wortmann's avatar
Daniel Wortmann committed
14
    USE m_types
15 16
    IMPLICIT NONE

17 18 19 20 21 22
    TYPE(t_atoms), INTENT(IN)  :: atoms
    REAL,          INTENT(IN)  :: f0(atoms%n_u),f2(atoms%n_u),f4(atoms%n_u),f6(atoms%n_u)
    REAL,          INTENT(OUT) :: u(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
                                    -lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)

    INTEGER, PARAMETER         :: lmaxw=3,lmmaxw1=(2*lmaxw+2)**2
23

Daniel Wortmann's avatar
Daniel Wortmann committed
24 25
    INTEGER i,j,k,l,m,mk,nfk,itype,i_u
    INTEGER m1,m2,m3,m4,lm1,lm2,lm3,lm4,kf
26
    REAL    uk,uq,avu,avj,cgk1,cgk2,tol
27
    REAL    fk(lmaxU_const+1,atoms%n_u)
28 29 30 31 32 33
    REAL,   ALLOCATABLE :: c(:,:,:)
    !
    tol = 1.0e-14
    !
    ! transformation to Hr-units:
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
34 35 36 37 38 39 40 41 42 43 44 45
    DO i_u = 1, atoms%n_u
       itype = atoms%lda_u(i_u)%atomType
       l = atoms%lda_u(i_u)%l
       fk(1,i_u) = f0(i_u) / hartree_to_ev_const
       fk(2,i_u) = f2(i_u) / hartree_to_ev_const
       fk(3,i_u) = f4(i_u) / hartree_to_ev_const
       IF (l.EQ.3) THEN
          fk(4,i_u) = f6(i_u) / hartree_to_ev_const
       ELSE IF (l.GT.3) THEN
          CALL juDFT_error("LDA+U for p, d or f-states!", calledby="umtx")
       END IF
    END DO
46 47 48
    !
    ! evaluate Gaunt parameter
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
49
    ALLOCATE(c(0:2*lmaxw+1,lmmaxw1,lmmaxw1))
50 51 52 53
    DO k = 1,lmmaxw1
       DO j = 1,lmmaxw1
          DO i = 0,2*lmaxw+1
             c(i,j,k) = 0.0
Daniel Wortmann's avatar
Daniel Wortmann committed
54 55 56 57
          END DO
       END DO
    END DO

58
    CALL sgaunt(lmaxw,lmmaxw1,lmaxU_const,c)
Daniel Wortmann's avatar
Daniel Wortmann committed
59 60 61

    DO i_u = 1, atoms%n_u                     !!! over U parameters
       l = atoms%lda_u(i_u)%l
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
       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
Daniel Wortmann's avatar
Daniel Wortmann committed
82
                      END DO                   ! mk
83 84
                      IF (ABS(uq).LT.tol) CYCLE
                      nfk=k/2+1
Daniel Wortmann's avatar
Daniel Wortmann committed
85 86 87 88 89 90 91
                      uk=uk+uq*fk(nfk,i_u)*4*pi_const/(2*k+1)
                   END DO                     ! k
                   u(m1,m2,m3,m4,i_u)=uk
                END DO                       ! m4 etc.
             END DO
          END DO
       END DO
92 93 94 95 96
       avu=0.e0
       avj=0.e0

       DO i = -l,l
          DO j = -l,l
Daniel Wortmann's avatar
Daniel Wortmann committed
97 98 99 100
             avu = avu+u(i,j,i,j,i_u)
             avj = avj+(u(i,j,i,j,i_u)-u(i,j,j,i,i_u))
          END DO
       END DO
101 102 103 104
       avu = avu/(2*l+1)/(2*l+1)
       avj = avj/(2*l+1)/(2*l)
       avj = avu-avJ
       !        WRITE (6,*) 'U-matr:'
Daniel Wortmann's avatar
Daniel Wortmann committed
105 106
       !        IF (l.eq.2) WRITE (6,111) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
       !        IF (l.eq.3) WRITE (6,211) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
107
       !        WRITE (6,*) 'J-matr:'
Daniel Wortmann's avatar
Daniel Wortmann committed
108 109
       !        IF (l.eq.2) WRITE (6,111) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
       !        IF (l.eq.3) WRITE (6,211) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
Daniel Wortmann's avatar
Daniel Wortmann committed
110 111
       !         PRINT*,'U-av:',avu*hartree_to_ev_const
       !         PRINT*,'J-av:',avj*hartree_to_ev_const
112 113 114
111    FORMAT (5f8.4)
211    FORMAT (7f8.4)
112    FORMAT (10e20.10)
Daniel Wortmann's avatar
Daniel Wortmann committed
115 116
       !c         WRITE (9,112) ((((u(i,j,k,m,i_u),i=-l,l),j=-l,l),k=-l,l),m=-l,l)
    END DO
117 118 119 120
    DEALLOCATE (c)

  END SUBROUTINE umtx
END MODULE m_umtx