Line data Source code
1 : MODULE m_local_Hamiltonian
2 : USE m_judft
3 : IMPLICIT NONE
4 : PRIVATE
5 : PUBLIC:: local_ham
6 : !*********************************************************************
7 : ! sets up the local Hamiltonian, i.e. the Hamiltonian in the
8 : ! l',m',l,m,u- basis which is independent from k!
9 : ! shifts this local Hamiltonian to make it positive definite
10 : ! and does a cholesky decomposition
11 : !*********************************************************************
12 : CONTAINS
13 748 : SUBROUTINE local_ham(sphhar,atoms,sym,noco,nococonv,enpara,&
14 : fmpi,v,vx,inden,input,hub1inp,hub1data,td,ud,alpha_hybrid,l_dfptmod)
15 : !! Should probably be called tlmplm_postprocess or something, as there is a
16 : !! lot more happening than only the Cholesky decomposition.
17 :
18 : ! Add onto the t_{L'L}^{\mu} matrices from tlmplm some contributions from
19 : ! DFT+U, DFT+HIA, DFT+OPC, constraint fields and the diagonal E_{l} terms
20 : ! etc. In the diagonal case, Cholesky decompose the nonspherical part of
21 : ! the Hamiltonian by shifting the diagonal part upwards until the matrix
22 : ! is positive-definite.
23 : USE m_spnorb
24 : USE m_tlmplm
25 : USE m_types
26 :
27 : TYPE(t_mpi), INTENT(IN) :: fmpi
28 : TYPE(t_noco), INTENT(IN) :: noco
29 : TYPE(t_nococonv), INTENT(IN) :: nococonv
30 : TYPE(t_input), INTENT(IN) :: input
31 : TYPE(t_sphhar), INTENT(IN) :: sphhar
32 : TYPE(t_atoms), INTENT(IN) :: atoms
33 : TYPE(t_sym), INTENT(IN) :: sym
34 : TYPE(t_enpara), INTENT(IN) :: enpara
35 : TYPE(t_hub1inp), INTENT(IN) :: hub1inp
36 : TYPE(t_hub1data), INTENT(INOUT) :: hub1data
37 : TYPE(t_potden), INTENT(IN) :: v,vx,inden
38 : TYPE(t_tlmplm), INTENT(INOUT) :: td
39 : TYPE(t_usdus), INTENT(INOUT) :: ud
40 :
41 : ! Scalar Arguments
42 :
43 : REAL, INTENT(IN) :: alpha_hybrid
44 :
45 : LOGICAL, INTENT(IN),OPTIONAL :: l_dfptmod
46 :
47 : ! Local Scalars
48 : INTEGER :: l,lm,j1,j2,jsp
49 : INTEGER :: n,m,s
50 : COMPLEX :: one
51 :
52 748 : CALL timestart("local_hamiltonian")
53 2150 : CALL td%init(atoms,input%jspins,(noco%l_noco.AND.noco%l_soc.AND..NOT.noco%l_ss).OR.any(noco%l_constrained))
54 :
55 4578 : DO jsp=1,MERGE(4,input%jspins,any(noco%l_unrestrictMT).OR.any(noco%l_spinoffd_ldau))
56 :
57 4312 : DO j1=merge(jsp,1,jsp<3),merge(jsp,2,jsp<3)
58 1232 : j2 = MERGE(j1,3-j1,jsp<3)
59 1232 : one = MERGE(CMPLX(1.,0.),CMPLX(0.,1.),jsp<4)
60 1232 : one = MERGE(CONJG(one),one,j1<j2)
61 :
62 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(l,m,lm,s)&
63 1232 : !$OMP SHARED(atoms,sphhar,sym,enpara,nococonv,j1,j2,jsp,fmpi,v,vx,input,hub1inp,hub1data,td,ud,alpha_hybrid,one,l_dfptmod,noco)
64 : DO n = 1,atoms%ntype
65 : IF (j1==j2.OR.noco%l_unrestrictMT(n)) THEN
66 : CALL tlmplm(n,sphhar,atoms,sym,enpara,nococonv,j1,j2,jsp,fmpi,v,vx,input,hub1inp,hub1data,td,ud,alpha_hybrid,one,PRESENT(l_dfptmod))
67 : END IF
68 : !Copy local hamiltonian for non_spherical setup
69 : call restrict_to_lnonsph(td%h_loc(:,:,n,j1,j2),td%h_loc2(n),td%h_loc2_nonsph(n),td%h_loc_nonsph(:,:,n,j1,j2))
70 : ! Now add diagonal contributions to the matrices:
71 : IF (jsp<3) THEN
72 : DO l = 0,atoms%lmax(n)
73 : DO m = -l,l
74 : lm = l*(l+1) + m
75 : s = td%h_loc2(n)
76 : td%h_loc(lm,lm,n,jsp,jsp) = td%h_loc(lm,lm,n,jsp,jsp) + enpara%el0(l,n,jsp)
77 : td%h_loc(lm,lm+s,n,jsp,jsp) = td%h_loc(lm,lm+s,n,jsp,jsp) + 0.5 ! Symmetrized from 1.0
78 : td%h_loc(lm+s,lm,n,jsp,jsp) = td%h_loc(lm+s,lm,n,jsp,jsp) + 0.5 ! Symmetrized from 0.0
79 : td%h_loc(lm+s,lm+s,n,jsp,jsp) = td%h_loc(lm+s,lm+s,n,jsp,jsp) + enpara%el0(l,n,jsp)*ud%ddn(l,n,jsp)
80 : END DO
81 : END DO
82 : END IF
83 : ! Store Matrices needed for LOs
84 : if (atoms%nlotot>0) call restrict_to_lnonsph(td%h_loc(:,:,n,j1,j2),td%h_loc2(n),td%h_loc2_nonsph(n),td%h_loc_LO(:,:,n,j1,j2))
85 : ENDDO
86 : !$OMP end parallel do
87 : !Add LDA+U
88 1232 : call add_ldaU(fmpi,inden,jsp,atoms,v,ud,input,hub1data,enpara,td%h_loc_nonsph,td%h_loc2_nonsph,j1,j2,.false.)
89 : !Add LDA+U also to LO part
90 1232 : if (atoms%nlotot>0) call add_ldaU(fmpi,inden,jsp,atoms,v,ud,input,hub1data,enpara,td%h_loc_LO,td%h_loc2_nonsph,j1,j2,.true.)
91 : ! Create Cholesky decomposition of local hamiltonian
92 : ! For DFPT, do not decompose!
93 1232 : IF (jsp<3.AND..NOT.PRESENT(l_dfptmod)) THEN
94 1144 : call cholesky_decompose(td%h_loc_nonsph(:,:,:,j1,j2),td%e_shift(:,jsp),atoms,ud,jsp)
95 : endif
96 :
97 4558 : IF (any(noco%l_constrained)) CALL tlmplm_constrained(atoms,v,enpara,input,hub1data,ud,nococonv,td)
98 : END DO
99 : END DO
100 :
101 : !Setup of soc parameters for first-variation SOC
102 748 : IF (noco%l_soc.AND.noco%l_noco.AND..NOT.noco%l_ss) THEN
103 64 : CALL spnorb(atoms,noco,nococonv,input,fmpi,enpara,v%mt,ud,td%rsoc,.FALSE.,hub1inp,hub1data)
104 : END IF
105 748 : CALL timestop("local_hamiltonian")
106 748 : END SUBROUTINE
107 :
108 0 : SUBROUTINE tlmplm_constrained(atoms,v,enpara,input,hub1data,ud,nococonv,td)
109 : USE m_radovlp
110 : USE m_types
111 : USE m_tlmplm
112 : TYPE(t_input), INTENT(IN) :: input
113 : TYPE(t_atoms), INTENT(IN) :: atoms
114 : TYPE(t_enpara), INTENT(IN) :: enpara
115 : TYPE(t_potden), INTENT(IN) :: v
116 : TYPE(t_tlmplm), INTENT(INOUT) :: td
117 : TYPE(t_usdus), INTENT(INOUT) :: ud
118 : TYPE(t_nococonv), INTENT(IN) :: nococonv
119 : TYPE(t_hub1data), INTENT(IN) :: hub1data
120 :
121 0 : REAL, ALLOCATABLE :: uun21(:,:),udn21(:,:),dun21(:,:),ddn21(:,:)
122 :
123 : COMPLEX :: c
124 : INTEGER :: n,l,s
125 :
126 : ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),udn21(0:atoms%lmaxd,atoms%ntype), &
127 0 : & dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype))
128 :
129 0 : CALL rad_ovlp(atoms,ud,input,hub1data,v%mt,enpara%el0,uun21,udn21,dun21,ddn21)
130 :
131 0 : DO n = 1,atoms%ntype
132 : ! In a constraint calculation, we have to calculate the local spin
133 : ! off-diagonal contributions
134 :
135 0 : s=atoms%lnonsph(n)+1
136 : ! First ispin=1,jspin=2 case
137 0 : DO l=0, atoms%lnonsph(n)
138 0 : c = (-0.5)*CMPLX(nococonv%b_con(1,n),nococonv%b_con(2,n))
139 0 : td%h_off(l ,l ,n,1,2) =td%h_off(l ,l ,n,1,2) + uun21(l,n)*c
140 0 : td%h_off(l ,l+s,n,1,2) =td%h_off(l ,l+s,n,1,2) + udn21(l,n)*c
141 0 : td%h_off(l+s,l ,n,1,2) =td%h_off(l+s,l ,n,1,2) + dun21(l,n)*c
142 0 : td%h_off(l+s,l+s,n,1,2) =td%h_off(l+s,l+s,n,1,2) + ddn21(l,n)*c
143 : END DO
144 :
145 : ! Then ispin=2,jspin=1 case
146 0 : DO l = 0, atoms%lnonsph(n)
147 0 : c = (-0.5)*CMPLX(nococonv%b_con(1,n),-nococonv%b_con(2,n))
148 0 : td%h_off(l ,l ,n,2,1) =td%h_off(l ,l ,n,2,1) + uun21(l,n)*c
149 0 : td%h_off(l ,l+s,n,2,1) =td%h_off(l ,l+s,n,2,1) + udn21(l,n)*c
150 0 : td%h_off(l+s,l ,n,2,1) =td%h_off(l+s,l ,n,2,1) + dun21(l,n)*c
151 0 : td%h_off(l+s,l+s,n,2,1) =td%h_off(l+s,l+s,n,2,1) + ddn21(l,n)*c
152 : END DO
153 : END DO
154 0 : END SUBROUTINE tlmplm_constrained
155 :
156 2024 : SUBROUTINE add_ldaU(fmpi,inden,jsp,atoms,v,ud,input,hub1data,enpara,mat,mat_half,j1,j2,l_lomatrix)
157 : ! Include contribution from LDA+U and LDA+HIA (latter are behind LDA+U contributions)
158 : USE m_radovlp
159 : USE m_opc_setup
160 : TYPE(t_mpi), INTENT(IN) :: fmpi
161 : TYPE(t_input), INTENT(IN) :: input
162 : TYPE(t_atoms), INTENT(IN) :: atoms
163 : TYPE(t_enpara), INTENT(IN) :: enpara
164 : TYPE(t_hub1data), INTENT(IN) :: hub1data
165 : TYPE(t_usdus), INTENT(INOUT) :: ud
166 : TYPE(t_potden), INTENT(IN) :: v,inden
167 : INTEGER,INTENT(IN) :: jsp,j1,j2
168 : COMPLEX,INTENT(INOUT) :: mat(:,:,:,:,:)
169 : INTEGER,INTENT(IN) :: mat_half(:)
170 : LOGICAL,INTENT(IN) :: l_lomatrix
171 :
172 : INTEGER :: i_u,n,s,l,lp,m,lm,mp,lmp,i_opc
173 2024 : REAL, ALLOCATABLE :: uun21(:,:),udn21(:,:),dun21(:,:),ddn21(:,:)
174 2024 : REAL, ALLOCATABLE :: opc_corrections(:)
175 :
176 2024 : IF(atoms%n_u+atoms%n_hia+atoms%n_opc==0) return !No LDA+U
177 :
178 274 : IF (j1.ne.j2) THEN
179 : !Calculate overlap integrals for the off-diagonal LDA+U contributions
180 0 : ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),source=0.0)
181 0 : ALLOCATE(udn21(0:atoms%lmaxd,atoms%ntype),source=0.0)
182 0 : ALLOCATE(dun21(0:atoms%lmaxd,atoms%ntype),source=0.0)
183 0 : ALLOCATE(ddn21(0:atoms%lmaxd,atoms%ntype),source=0.0)
184 0 : CALL rad_ovlp(atoms,ud,input,hub1data,v%mt,enpara%el0, uun21,udn21,dun21,ddn21)
185 : ELSE
186 274 : IF (atoms%n_opc > 0) THEN
187 : CALL opc_setup(input,atoms,fmpi,v,inden,jsp,opc_corrections)
188 : END IF
189 :
190 : END IF
191 :
192 : !!$OMP PARALLEL DO DEFAULT(NONE) private(n,s,l,lp,m,lm,mp,lmp)&
193 : !!$OMP shared(atoms,mat_half,l_lomatrix,v,mat,uun21,udn21,dun21,ddn21,j1,j2,ud,jsp)
194 650 : DO i_u=1,atoms%n_u+atoms%n_hia
195 376 : if (l_lomatrix.and..not.atoms%lda_u(i_u)%use_lo) cycle
196 256 : n=atoms%lda_u(i_u)%atomType
197 256 : s=mat_half(n)
198 : ! Found a "U" for this atom type
199 256 : l = atoms%lda_u(i_u)%l
200 256 : lp = atoms%lda_u(i_u)%l
201 :
202 1722 : DO m = -l,l
203 1192 : lm = l* (l+1) + m +1 !indexing from 1
204 7456 : DO mp = -lp,lp
205 5888 : lmp = lp*(lp+1) + mp +1 !indexing from 1
206 7080 : IF (j1==j2) THEN
207 5888 : mat(lm ,lmp ,n,j1,j2) = mat(lm ,lmp ,n,j1,j2) + v%mmpMat(m,mp,i_u,jsp)
208 5888 : mat(lm+s,lmp+s,n,j1,j2) = mat(lm+s,lmp+s,n,j1,j2) + v%mmpMat(m,mp,i_u,jsp) * ud%ddn(lp,n,jsp)
209 0 : ELSE IF(j1>j2) THEN
210 0 : mat(lm ,lmp ,n,j1,j2) = mat(lm ,lmp ,n,j1,j2) + v%mmpMat(m,mp,i_u,3) * uun21(l,n)
211 0 : mat(lm ,lmp+s,n,j1,j2) = mat(lm ,lmp+s,n,j1,j2) + v%mmpMat(m,mp,i_u,3) * udn21(l,n)
212 0 : mat(lm+s,lmp ,n,j1,j2) = mat(lm+s,lmp ,n,j1,j2) + v%mmpMat(m,mp,i_u,3) * dun21(l,n)
213 0 : mat(lm+s,lmp+s,n,j1,j2) = mat(lm+s,lmp+s,n,j1,j2) + v%mmpMat(m,mp,i_u,3) * ddn21(l,n)
214 : ELSE
215 : ! For this part of the Hamiltonian we need to perform Hermitian conjugation on mmpMat
216 0 : mat(lm ,lmp ,n,j1,j2) = mat(lm ,lmp ,n,j1,j2) + conjg(v%mmpMat(mp,m,i_u,3)) * uun21(l,n)
217 0 : mat(lm ,lmp+s,n,j1,j2) = mat(lm ,lmp+s,n,j1,j2) + conjg(v%mmpMat(mp,m,i_u,3)) * udn21(l,n)
218 0 : mat(lm+s,lmp ,n,j1,j2) = mat(lm+s,lmp ,n,j1,j2) + conjg(v%mmpMat(mp,m,i_u,3)) * dun21(l,n)
219 0 : mat(lm+s,lmp+s,n,j1,j2) = mat(lm+s,lmp+s,n,j1,j2) + conjg(v%mmpMat(mp,m,i_u,3)) * ddn21(l,n)
220 : END IF
221 : END DO
222 : END DO
223 : END DO
224 : !!$OMP end parallel do
225 466 : DO i_opc=1,atoms%n_opc
226 192 : n=atoms%lda_opc(i_opc)%atomType
227 192 : s=mat_half(n)
228 : ! Found an "OPC" for this atom type
229 192 : l=atoms%lda_opc(i_opc)%l
230 1426 : DO m = -l,l
231 960 : lm = l*(l+1) + m +1 !indexing from 1
232 960 : mat(lm ,lm ,n,j1,j2) = mat(lm ,lm ,n,j1,j2) + opc_corrections(i_opc) * m
233 1152 : mat(lm+s,lm+s,n,j1,j2) = mat(lm+s,lm+s,n,j1,j2) + opc_corrections(i_opc) * m * ud%ddn(l,n,jsp)
234 : END DO
235 : END DO
236 274 : END SUBROUTINE
237 :
238 1144 : subroutine cholesky_decompose(matrix,e_shift,atoms,ud,jsp)
239 : USE m_types
240 : TYPE(t_atoms), INTENT(IN) :: atoms
241 : TYPE(t_usdus),INTENT(IN) :: ud
242 : COMPLEX,INTENT(INOUT) :: matrix(0:,0:,:)
243 : REAL, INTENT(OUT) :: e_shift(:)
244 : INTEGER,INTENT(IN) :: jsp
245 :
246 : REAL, PARAMETER :: e_shift_min=0.5
247 : REAL, PARAMETER :: e_shift_max=65.0
248 :
249 : INTEGER :: n,info,s,l,lp,lmp,mp
250 1144 : COMPLEX,ALLOCATABLE :: mat(:,:)
251 :
252 3146 : e_shift=e_shift_min
253 :
254 3146 : DO n=1,atoms%ntype
255 2002 : s=atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+1
256 2002 : info=1
257 4004 : cholesky_loop:DO WHILE(info.ne.0)
258 28372208 : mat=matrix(0:,0:,n)
259 : !Mat is now using a lower bound of 1!!
260 : ! Add shift onto the diagonal terms to make matrix positive definite
261 16392 : DO lp = 0,atoms%lnonsph(n)
262 123146 : DO mp = -lp,lp
263 106754 : lmp = lp* (lp+1) + mp +1
264 106754 : mat(lmp,lmp)=e_shift(n)+mat(lmp,lmp)
265 121144 : mat(lmp+s,lmp+s)=e_shift(n)*ud%ddn(lp,n,jsp)+mat(lmp+s,lmp+s)
266 : END DO
267 : END DO
268 2002 : IF (lmp.NE.s) CALL judft_error("BUG in local_Hamiltonian:cholesky")
269 :
270 : ! Perform cholesky decomposition
271 2002 : CALL zpotrf("L",2*s,mat(:,:),SIZE(mat,1),info)
272 :
273 : ! Set upper part to zero
274 215510 : DO l=1,2*s
275 12912184 : DO lp=1,l-1
276 12910182 : mat(lp,l)=0.0
277 : END DO
278 : END DO
279 :
280 4004 : IF (info.NE.0) THEN
281 0 : e_shift(n)=e_shift(n)*2.0
282 0 : IF (e_shift(n)>e_shift_max) THEN
283 0 : CALL judft_error("Potential shift at maximum")
284 : END IF
285 : END IF
286 : END DO cholesky_loop
287 28371350 : matrix(0:,0:,n)=mat
288 : ENDDO
289 1144 : END SUBROUTINE
290 :
291 :
292 3484 : subroutine restrict_to_lnonsph(mat,s2,s,mat_nonsph)
293 : COMPLEX,INTENT(IN) :: mat(0:,0:)
294 : INTEGER,INTENT(IN) :: s,s2
295 : COMPLEX,INTENT(OUT) :: mat_nonsph(0:,0:)
296 : ! Set up local hamiltonian
297 11909124 : mat_nonsph(0:s-1,0:s-1) = mat(0:s-1,0:s-1)
298 11909124 : mat_nonsph(s:s+s-1,0:s-1) = mat(s2:s+s2-1,0:s-1)
299 11909124 : mat_nonsph(0:s-1,s:s+s-1) = mat(0:s-1,s2:s+s2-1)
300 11909124 : mat_nonsph(s:s+s-1,s:s+s-1)= mat(s2:s+s2-1,s2:s+s2-1)
301 3484 : end subroutine
302 :
303 : END MODULE m_local_Hamiltonian
|
