umtx.f90 3.94 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
Daniel Wortmann's avatar
Daniel Wortmann committed
9 10
  SUBROUTINE umtx(atoms,lmaxb,f0,f2,f4,f6,&
                  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 17
    IMPLICIT NONE

    INTEGER, PARAMETER   :: lmaxw=3,lmmaxw1=(2*lmaxw+2)**2
Daniel Wortmann's avatar
Daniel Wortmann committed
18 19 20 21
    TYPE(t_atoms),  INTENT(IN)  :: atoms
    INTEGER, INTENT (IN) :: lmaxb
    REAL,    INTENT (IN) :: f0(atoms%n_u),f2(atoms%n_u),f4(atoms%n_u),f6(atoms%n_u)
    REAL,    INTENT (OUT) :: u(-lmaxb:lmaxb,-lmaxb:lmaxb,-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u)
22

Daniel Wortmann's avatar
Daniel Wortmann committed
23 24
    INTEGER i,j,k,l,m,mk,nfk,itype,i_u
    INTEGER m1,m2,m3,m4,lm1,lm2,lm3,lm4,kf
25
    REAL    uk,uq,avu,avj,cgk1,cgk2,tol
Daniel Wortmann's avatar
Daniel Wortmann committed
26
    REAL    fk(lmaxb+1,atoms%n_u)
27 28 29 30 31 32
    REAL,   ALLOCATABLE :: c(:,:,:)
    !
    tol = 1.0e-14
    !
    ! transformation to Hr-units:
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
33 34 35 36 37 38 39 40 41 42 43 44
    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
45 46 47
    !
    ! evaluate Gaunt parameter
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
48
    ALLOCATE(c(0:2*lmaxw+1,lmmaxw1,lmmaxw1))
49 50 51 52
    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
53 54 55 56 57 58 59 60
          END DO
       END DO
    END DO

    CALL sgaunt(lmaxw,lmmaxw1,lmaxb,c)

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

       DO i = -l,l
          DO j = -l,l
Daniel Wortmann's avatar
Daniel Wortmann committed
96 97 98 99
             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
100 101 102 103
       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
104 105
       !        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)
106
       !        WRITE (6,*) 'J-matr:'
Daniel Wortmann's avatar
Daniel Wortmann committed
107 108
       !        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
109 110
       !         PRINT*,'U-av:',avu*hartree_to_ev_const
       !         PRINT*,'J-av:',avj*hartree_to_ev_const
111 112 113
111    FORMAT (5f8.4)
211    FORMAT (7f8.4)
112    FORMAT (10e20.10)
Daniel Wortmann's avatar
Daniel Wortmann committed
114 115
       !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
116 117 118 119
    DEALLOCATE (c)

  END SUBROUTINE umtx
END MODULE m_umtx