Line data Source code
1 : MODULE m_add_selfen
2 :
3 : USE m_types
4 : USE m_types_selfen
5 : USE m_constants
6 : USE m_juDFT
7 :
8 : IMPLICIT NONE
9 :
10 : CONTAINS
11 :
12 0 : SUBROUTINE add_selfen(g0,selfen,gfinp,input,atoms,noco,nococonv,occDFT,g,mmpMat)
13 :
14 : !Calculates the interacting Green's function for the mt-sphere with
15 : !
16 : ! (G)^-1 = (G_0)^-1 - mu 1 + V_dc - selfen + V_U
17 : !
18 : !The term mu * unity is there to ensure that the number of particles
19 : !doesn't change and is determined by a two-step process
20 : !The occupation as a function of mu has a peak in the region where
21 : !something is inside the energy interval between e_bot and e_fermi
22 : !To determine where we have the same number of particles we first
23 : !search for the maximum occupation
24 : !Then the desired chemical potential is found with the bisection method
25 : !to the right of the maximum
26 : !TODO: Parallelization (OMP over chemical potentials in first loop??)
27 :
28 : TYPE(t_greensf), INTENT(IN) :: g0
29 : TYPE(t_gfinp), INTENT(IN) :: gfinp
30 : TYPE(t_input), INTENT(IN) :: input
31 : TYPE(t_atoms), INTENT(IN) :: atoms
32 : TYPE(t_noco), INTENT(IN) :: noco
33 : TYPE(t_nococonv), INTENT(IN) :: nococonv
34 : REAL, INTENT(IN) :: occDFT(:)
35 : TYPE(t_selfen), INTENT(INOUT) :: selfen
36 : TYPE(t_greensf), INTENT(INOUT) :: g
37 : COMPLEX, INTENT(INOUT) :: mmpMat(-lmaxU_const:,-lmaxU_const:,:)
38 :
39 : INTEGER :: l,ispin,m,mp,iMatch
40 : REAL :: mu_a,mu_b,mu_step
41 : REAL :: mu,nocc,nTarget,muMax,nMax
42 : LOGICAL :: l_fullMatch,l_invalidElements
43 :
44 : !Are we matching the spin polarized self-energy with one chemical potential
45 0 : l_fullMatch = SIZE(selfen%muMatch).EQ.1 .AND. input%jspins.EQ.2
46 :
47 0 : l = g0%elem%l
48 :
49 : !Search for the maximum of occupation
50 0 : DO iMatch = 1, SIZE(selfen%muMatch)
51 :
52 : !Target occupation
53 0 : nTarget = MERGE(SUM(occDFT(:)),occDFT(iMatch),l_fullMatch)
54 :
55 : !Interval where we expect the correct mu
56 0 : mu_a = -2.0
57 0 : mu_b = 1.5
58 0 : mu_step = 0.01
59 :
60 : !----------------------------------------------------
61 : ! Scan the given Interval for the maximum Occupation
62 : !----------------------------------------------------
63 0 : mu = mu_a
64 0 : muMax = 0.0
65 0 : nMax = 0.0
66 0 : DO WHILE(mu.LE.mu_b)
67 :
68 0 : mu = mu + mu_step
69 :
70 : CALL getOccupationMtx(g0,gfinp,input,atoms,noco,nococonv,selfen,mu,l_fullMatch,iMatch,&
71 0 : g,mmpMat,nocc,l_invalidElements)
72 :
73 0 : IF(nocc.GT.nMax) THEN
74 0 : muMax = mu
75 0 : nMax = nocc
76 : ENDIF
77 :
78 : #ifdef CPP_DEBUG
79 : WRITE(1337,'(2f15.8)') mu,nocc
80 : #endif
81 :
82 : ENDDO
83 :
84 : !Sanity check for the maximum occupation
85 0 : IF(nMax-2*(2*l+1).GT.1.0) THEN
86 : !These oscillations seem to emerge when the lorentzian smoothing is done inadequately
87 : CALL juDFT_error("Something went wrong with the addition of the selfenergy: nMax>>2*(2*l+1)",&
88 0 : calledby="add_selfen")
89 0 : ELSE IF(nMax-nTarget.LT.-0.1) THEN
90 : CALL juDFT_error("Something went wrong with the addition of the selfenergy: nMax<nTarget",&
91 0 : calledby="add_selfen")
92 : ENDIF
93 :
94 : !------------------------------------------------------------------
95 : ! Find the matching chemical potential on the right of the maximum
96 : !------------------------------------------------------------------
97 : !Set up the interval for the bisection method (mu_max,mu_b)
98 : mu_a = muMax
99 0 : DO WHILE (ABS(nocc-nTarget).GT.1e-8.AND.ABS((mu_b - mu_a)/2.0).GT.1e-8)
100 0 : mu = (mu_a + mu_b)/2.0
101 :
102 : CALL getOccupationMtx(g0,gfinp,input,atoms,noco,nococonv,selfen,mu,l_fullMatch,iMatch,&
103 0 : g,mmpMat,nocc,l_invalidElements)
104 :
105 0 : IF((nocc - nTarget).GT.0.0) THEN
106 : !The occupation is to big --> choose the right interval
107 : mu_a = mu
108 0 : ELSE IF((nocc - nTarget).LT.0.0) THEN
109 : !The occupation is to small --> choose the left interval
110 0 : mu_b = mu
111 : ENDIF
112 : ENDDO
113 0 : selfen%muMatch(iMatch) = mu
114 : !----------------------------------------------------
115 : ! Check if the final mmpMat contains invalid elements
116 : !----------------------------------------------------
117 0 : IF(l_invalidElements) THEN
118 0 : CALL juDFT_error("Invalid Element/occupation in final density matrix",calledby="add_selfen")
119 : ENDIF
120 : ENDDO
121 :
122 : !Test throw out elements smaller than 1e-4
123 0 : DO ispin = 1, MERGE(3,input%jspins,gfinp%l_mperp)
124 0 : DO m = -l, l
125 0 : DO mp =-l, l
126 0 : IF(ABS(mmpMat(m,mp,ispin)).LT.1e-4) mmpMat(m,mp,ispin) = cmplx_0
127 : ENDDO
128 : ENDDO
129 : ENDDO
130 :
131 0 : END SUBROUTINE add_selfen
132 :
133 0 : SUBROUTINE getOccupationMtx(g0,gfinp,input,atoms,noco,nococonv,selfen,mu,l_fullMatch,iMatch,&
134 0 : g,mmpMat,nocc,l_invalidElements)
135 :
136 :
137 : TYPE(t_greensf), INTENT(IN) :: g0 !DFT Green's Function
138 : TYPE(t_gfinp), INTENT(IN) :: gfinp
139 : TYPE(t_input), INTENT(IN) :: input
140 : TYPE(t_atoms), INTENT(IN) :: atoms
141 : TYPE(t_noco), INTENT(IN) :: noco
142 : TYPE(t_nococonv), INTENT(IN) :: nococonv
143 : TYPE(t_selfen), INTENT(IN) :: selfen !Atomic self-energy (with removed LDA+U potential)
144 : REAL, INTENT(IN) :: mu !chemical potential shift
145 : LOGICAL, INTENT(IN) :: l_fullMatch !Are spins matched individually?
146 : INTEGER, INTENT(IN) :: iMatch !Index for current matching
147 : TYPE(t_greensf), INTENT(INOUT) :: g !Impurity Green's Function
148 : COMPLEX, INTENT(INOUT) :: mmpMat(-lmaxU_const:,-lmaxU_const:,:) !Occupation matrix of g
149 : REAL, INTENT(INOUT) :: nocc !trace of the occupation matrix
150 : LOGICAL, INTENT(INOUT) :: l_invalidElements !Are there invalid elements in the resulting Occupation matrix
151 :
152 : INTEGER :: l,ns,matsize,start,end
153 : INTEGER :: ipm,iz,ispin,m
154 :
155 0 : TYPE(t_mat) :: vmat,gmat
156 :
157 0 : l = g0%elem%l
158 0 : ns = 2*l+1
159 0 : matsize = ns*MERGE(2,1,l_fullMatch)
160 0 : CALL vmat%init(.false.,matsize,matsize)
161 0 : CALL gmat%init(.false.,matsize,matsize)
162 0 : IF(iMatch>1) CALL g%reset()
163 :
164 : !Select the correct section from the selfenergy
165 0 : start = MERGE(1,1+(iMatch-1)*ns,l_fullMatch)
166 0 : end = MERGE(2*ns,iMatch*ns,l_fullMatch)
167 0 : DO ipm = 1, 2
168 0 : DO iz = 1, g0%contour%nz
169 : !Read selfenergy
170 0 : vmat%data_c = selfen%data(start:end,start:end,iz,ipm)
171 0 : IF(.NOT.gfinp%l_mperp.AND.l_fullMatch) THEN
172 : !Dismiss spin-off-diagonal elements
173 0 : vmat%data_c(1:ns,ns+1:2*ns) = cmplx_0
174 0 : vmat%data_c(ns+1:2*ns,1:ns) = cmplx_0
175 : ENDIF
176 :
177 : !Read in the DFT-Green's Function at the energy point
178 0 : IF(l_fullMatch) THEN
179 0 : CALL g0%getFullMatrix(atoms,iz,ipm.EQ.2,gmat)
180 : ELSE
181 0 : CALL g0%get(atoms,iz,ipm.EQ.2,iMatch,gmat)
182 : ENDIF
183 :
184 : !----------------------------------------------------
185 : !Solve the Dyson equation at the current energy point
186 : !----------------------------------------------------
187 0 : CALL add_pot(gmat,vmat,mu)
188 :
189 : !Set the Impurity-Green's Function at the energy point
190 0 : IF(l_fullMatch) THEN
191 0 : CALL g%set(iz,ipm.EQ.2,gmat)
192 : ELSE
193 0 : CALL g%set(iz,ipm.EQ.2,gmat,spin=iMatch)
194 : ENDIF
195 :
196 : ENDDO
197 : ENDDO
198 :
199 : !Get the occupation matrix
200 0 : IF(l_fullMatch) THEN
201 0 : mmpMat = g%occupationMatrix(gfinp,input,atoms,noco,nococonv,check=.TRUE.,occError=l_invalidElements)
202 : ELSE
203 0 : mmpMat(:,:,iMatch) = g%occupationMatrix(iMatch,gfinp,input,atoms,noco,nococonv,check=.TRUE.,occError=l_invalidElements)
204 : ENDIF
205 :
206 : !Compute the trace
207 0 : nocc = 0.0
208 0 : DO ispin = MERGE(1,iMatch,l_fullMatch), MERGE(input%jspins,iMatch,l_fullMatch)
209 0 : DO m = -l, l
210 0 : nocc = nocc + REAL(mmpMat(m,m,ispin))
211 : ENDDO
212 : ENDDO
213 :
214 :
215 0 : END SUBROUTINE getOccupationMtx
216 :
217 0 : SUBROUTINE add_pot(gmat,vmat,mu)
218 :
219 : TYPE(t_mat), INTENT(INOUT) :: gmat
220 : TYPE(t_mat), INTENT(IN) :: vmat
221 : REAL, INTENT(IN) :: mu
222 :
223 : INTEGER i
224 :
225 0 : IF(vmat%matsize1.NE.gmat%matsize1.OR.vmat%matsize2.NE.gmat%matsize2) &
226 0 : CALL juDFT_error("vmat & gmat dimension do not match",hint="This is a bug in FLEUR, please report",calledby="add_pot")
227 :
228 0 : CALL gmat%inverse()
229 0 : gmat%data_c = gmat%data_c - vmat%data_c
230 0 : DO i = 1, gmat%matsize1
231 0 : gmat%data_c(i,i) = gmat%data_c(i,i) - mu
232 : ENDDO
233 0 : CALL gmat%inverse()
234 :
235 0 : END SUBROUTINE add_pot
236 :
237 : END MODULE m_add_selfen
|
