Line data Source code
1 : module m_coulombPotential
2 :
3 : use m_types
4 : use m_uj2f
5 : use m_umtx
6 :
7 : implicit none
8 :
9 : contains
10 :
11 160 : subroutine coulombPotential(density, ldau, jspins, l_spinoffd, potential, interaction_energy)
12 :
13 : !This subroutine calculates the DFT+U potential matrix excluding the double counting
14 :
15 : COMPLEX, INTENT(IN) :: density(-lmaxU_const:,-lmaxU_const:,:)
16 : TYPE(t_utype), INTENT(IN) :: ldau
17 : INTEGER, INTENT(IN) :: jspins
18 : LOGICAL, INTENT(IN) :: l_spinoffd
19 : COMPLEX, INTENT(INOUT) :: potential(-lmaxU_const:,-lmaxU_const:,:)
20 : REAL, INTENT(OUT) :: interaction_energy
21 :
22 : INTEGER :: m,mp,q,p,ispin,jspin,spin_dim
23 : REAL :: spin_deg,f0,f2,f4,f6,energy_contribution,total_charge
24 : COMPLEX :: vu
25 : REAL, ALLOCATABLE :: umatrix(:,:,:,:)
26 160 : COMPLEX, ALLOCATABLE :: Vdc(:,:,:)
27 :
28 160 : spin_dim = SIZE(density,3)
29 160 : spin_deg = 1.0 / (3 - jspins)
30 :
31 160 : IF(SIZE(potential,3) /= spin_dim) CALL juDFT_error('Mismatch in dimensions between potential and density', calledby='dftUPotential')
32 :
33 160 : IF(l_spinoffd.AND.spin_dim < 3) THEN
34 0 : CALL juDFT_error('Spin offdiagonal parts missing', calledby='dftUPotential')
35 : ENDIF
36 :
37 160 : IF(.NOT.l_spinoffd) spin_dim = MIN(2,spin_dim)
38 :
39 :
40 160 : CALL uj2f(jspins,ldau,f0,f2,f4,f6)
41 :
42 : ALLOCATE ( umatrix(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
43 448160 : -lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const), source=0.0 )
44 12952 : ALLOCATE(Vdc(-lmaxU_const:lmaxU_const, -lmaxU_const:lmaxU_const, SIZE(density,3)), source=cmplx_0)
45 160 : CALL umtx(ldau,f0,f2,f4,f6,umatrix)
46 :
47 : !--------------------------------------------------------!
48 : ! s -- s' !
49 : ! V = > ( <m,p|V|m',q> - <m,p|V|q,m'> d ) n !
50 : ! m,m' -- s,s' p,q !
51 : ! p,q,s' !
52 : !--------------------------------------------------------!
53 12472 : potential = cmplx_0
54 :
55 376 : DO ispin = 1, spin_dim
56 1400 : DO m = -ldau%l,ldau%l
57 6384 : DO mp =-ldau%l,ldau%l
58 5144 : vu = cmplx_0
59 5144 : IF(ispin < 3) THEN
60 13664 : DO jspin = 1, jspins
61 8520 : IF (ispin == jspin) THEN
62 32280 : DO p = -ldau%l,ldau%l
63 181616 : DO q = -ldau%l,ldau%l
64 176472 : vu = vu + density(p, q,jspin) * ( umatrix(m,p,mp,q) - umatrix(m,p,q,mp) )
65 : END DO
66 : END DO
67 : END IF
68 13664 : IF (ispin /= jspin .OR. jspins == 1) THEN
69 32280 : DO p = -ldau%l,ldau%l
70 184992 : DO q = -ldau%l,ldau%l
71 176472 : vu = vu + umatrix(m,p,mp,q) * density(p, q,jspin)
72 : END DO
73 : END DO
74 : END IF
75 : END DO
76 : ELSE
77 0 : DO p = -ldau%l,ldau%l
78 0 : DO q = -ldau%l,ldau%l
79 0 : vu = vu - umatrix(m,p,q,mp) * density(p,q,ispin)
80 : ENDDO
81 : ENDDO
82 : ENDIF
83 6168 : potential(m,mp,ispin) = potential(m,mp,ispin) + vu
84 : END DO ! m' loop
85 : END DO ! m loop
86 : END DO ! outer spin loop
87 :
88 : !Take into account spin-degeneracy for non-spinpolarized calculations
89 12472 : potential = potential * spin_deg
90 :
91 : !----------------------------------------------------------------------+
92 : ! s !
93 : ! ee 1 --- s s 1 s 1 !
94 : ! E (n) = - > n ( V + d ( U (n - -) - J (n - -) )) !
95 : ! 2 --- m,m' m,m' m,m' 2 2 !
96 : ! m,m' !
97 : !----------------------------------------------------------------------+
98 :
99 160 : interaction_energy = 0.0
100 376 : DO ispin = 1,spin_dim
101 1400 : DO m = -ldau%l,ldau%l
102 6384 : DO mp = -ldau%l,ldau%l
103 6168 : IF(ispin < 3) THEN
104 5144 : interaction_energy = interaction_energy + REAL(potential(m,mp,ispin)*density(m,mp,ispin))
105 : ELSE
106 0 : DO p = -ldau%l,ldau%l
107 0 : DO q = -ldau%l,ldau%l
108 : interaction_energy = interaction_energy + umatrix(m,p,q,mp) *&
109 : REAL( density(m,mp,ispin)*conjg(density(q,p,ispin)) &
110 0 : + conjg(density(mp,m,ispin))*density(p,q,ispin) )
111 : ENDDO
112 : ENDDO
113 : ENDIF
114 : END DO
115 : END DO
116 : END DO
117 160 : interaction_energy = interaction_energy/2.0
118 :
119 320 : end subroutine
120 :
121 : end module m_coulombPotential
|
