Line data Source code
1 : MODULE m_umtx
2 : !*********************************************************************
3 : !* The calculation of the "U"-contribution to Hartree-Fock matrix. *
4 : !*-------------------------------------------------------------------*
5 : !* Extension to multiple U per atom type by G.M. 2017 *
6 : !*********************************************************************
7 : USE m_juDFT
8 : USE m_constants
9 : USE m_sgaunt
10 : USE m_types
11 :
12 : IMPLICIT NONE
13 :
14 : INTERFACE umtx
15 : PROCEDURE :: umtx_all, umtx_single
16 : END INTERFACE
17 :
18 : CONTAINS
19 :
20 160 : SUBROUTINE umtx_single(u_in,f0,f2,f4,f6,u)
21 :
22 : TYPE(t_utype), INTENT(IN) :: u_in
23 : REAL, INTENT(IN) :: f0,f2,f4,f6
24 : REAL, INTENT(OUT) :: u(-lmaxU_const:,-lmaxU_const:,-lmaxU_const:,-lmaxU_const:)
25 :
26 : REAL :: uAll(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,1)
27 : REAL :: f0List(1),f2List(1)
28 : REAL :: f4List(1),f6List(1)
29 :
30 160 : f0List(1) = f0; f2List(1) = f2
31 160 : f4List(1) = f4; f6List(1) = f6
32 :
33 320 : CALL umtx_all([u_in],1,f0List,f2List,f4List,f6List,uAll)
34 :
35 448160 : u = uAll(:,:,:,:,1)
36 :
37 160 : END SUBROUTINE umtx_single
38 :
39 160 : SUBROUTINE umtx_all(u_in,n_u,f0,f2,f4,f6,u)
40 :
41 : INTEGER, INTENT(IN) :: n_u
42 : TYPE(t_utype), INTENT(IN) :: u_in(:)
43 : REAL, INTENT(IN) :: f0(:),f2(:),f4(:),f6(:)
44 : REAL, INTENT(OUT) :: u(-lmaxU_const:,-lmaxU_const:,-lmaxU_const:,-lmaxU_const:,:)
45 :
46 : INTEGER, PARAMETER :: lmmaxw1=(2*lmaxU_const+2)**2
47 : REAL, PARAMETER :: tol=1e-14
48 :
49 : INTEGER :: i,j,k,l,m,mk,nfk,itype,i_u
50 : INTEGER :: m1,m2,m3,m4,lm1,lm2,lm3,lm4,kf
51 : REAL :: uk,uq,avu,avj,cgk1,cgk2
52 160 : REAL :: fk(lmaxU_const+1,n_u)
53 160 : REAL, ALLOCATABLE :: c(:,:,:)
54 :
55 : !
56 : ! transformation to Hr-units:
57 : !
58 320 : DO i_u = 1, n_u
59 160 : itype = u_in(i_u)%atomType
60 160 : l = u_in(i_u)%l
61 160 : fk(1,i_u) = f0(i_u) / hartree_to_ev_const
62 160 : fk(2,i_u) = f2(i_u) / hartree_to_ev_const
63 160 : fk(3,i_u) = f4(i_u) / hartree_to_ev_const
64 320 : IF (l.EQ.3) THEN
65 12 : fk(4,i_u) = f6(i_u) / hartree_to_ev_const
66 148 : ELSE IF (l.GT.3) THEN
67 0 : CALL juDFT_error("LDA+U for p, d or f-states!", calledby="umtx")
68 : END IF
69 : END DO
70 :
71 : !
72 : ! evaluate Gaunt parameter
73 : !
74 5908640 : ALLOCATE(c(0:2*lmaxU_const+1,lmmaxw1,lmmaxw1),source=0.0)
75 160 : CALL sgaunt(lmaxU_const,c)
76 :
77 320 : DO i_u = 1, n_u
78 160 : l = u_in(i_u)%l
79 160 : kf = 2*l
80 880 : DO m1 = -l,l
81 720 : lm1 = l*(l+1)+m1+1
82 4336 : DO m2 = -l,l
83 3456 : lm2 = l*(l+1)+m2+1
84 21696 : DO m3 = -l,l
85 17520 : lm3 = l*(l+1)+m3+1
86 114000 : DO m4 = -l,l
87 93024 : lm4 = l*(l+1)+m4+1
88 93024 : uk = 0.0e0
89 396696 : DO k=0,kf,2
90 303672 : uq = 0.e0
91 2035680 : DO mk=-k,k
92 1732008 : IF (mk.NE.m1-m3) CYCLE
93 198480 : cgk1 = c(k/2,lm1,lm3)
94 198480 : IF (ABS(cgk1).LT.tol) CYCLE
95 197304 : IF (mk.NE.m4-m2) CYCLE
96 27112 : cgk2 = c(k/2,lm4,lm2)
97 27112 : IF (ABS(cgk2).LT.tol) CYCLE
98 2035680 : uq = uq+cgk1*cgk2
99 : END DO ! mk
100 303672 : IF (ABS(uq).LT.tol) CYCLE
101 26992 : nfk=k/2+1
102 396696 : uk=uk+uq*fk(nfk,i_u)*4*pi_const/(2*k+1)
103 : END DO ! k
104 110544 : u(m1,m2,m3,m4,i_u)=uk
105 : END DO ! m4 etc.
106 : END DO
107 : END DO
108 : END DO
109 : avu=0.e0
110 : avj=0.e0
111 :
112 : DO i = -l,l
113 : DO j = -l,l
114 : avu = avu+u(i,j,i,j,i_u)
115 : avj = avj+(u(i,j,i,j,i_u)-u(i,j,j,i,i_u))
116 : END DO
117 : END DO
118 160 : avu = avu/(2*l+1)/(2*l+1)
119 160 : avj = avj/(2*l+1)/(2*l)
120 320 : avj = avu-avJ
121 : ! WRITE (oUnit,*) 'U-matr:'
122 : ! IF (l.eq.2) WRITE (oUnit,111) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
123 : ! IF (l.eq.3) WRITE (oUnit,211) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
124 : ! WRITE (oUnit,*) 'J-matr:'
125 : ! IF (l.eq.2) WRITE (oUnit,111) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
126 : ! IF (l.eq.3) WRITE (oUnit,211) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
127 : ! PRINT*,'U-av:',avu*hartree_to_ev_const
128 : ! PRINT*,'J-av:',avj*hartree_to_ev_const
129 : 111 FORMAT (5f8.4)
130 : 211 FORMAT (7f8.4)
131 : 112 FORMAT (10e20.10)
132 : !c WRITE (9,112) ((((u(i,j,k,m,i_u),i=-l,l),j=-l,l),k=-l,l),m=-l,l)
133 : END DO
134 160 : DEALLOCATE (c)
135 :
136 160 : END SUBROUTINE umtx_all
137 : END MODULE m_umtx
|
