Line data Source code
1 : MODULE m_ssomat
2 : USE m_judft
3 : IMPLICIT NONE
4 : CONTAINS
5 0 : SUBROUTINE ssomat(seigvso,h_so,theta,phi,eig_id,atoms,kpts,sym,&
6 0 : cell,noco,nococonv, input,fmpi, enpara,v,results,ef )
7 : USE m_types_nococonv
8 : USE m_types_mat
9 : USE m_types_setup
10 : USE m_types_mpi
11 : USE m_types_enpara
12 : USE m_types_potden
13 : USE m_types_misc
14 : USE m_types_kpts
15 : USE m_types_tlmplm
16 : USE m_types_usdus
17 : USE m_types_lapw
18 : USE m_constants
19 : USE m_eig66_io
20 : USE m_spnorb
21 : USE m_abcof
22 : USE m_fermifct
23 : #ifdef CPP_MPI
24 : USE mpi
25 : #endif
26 : IMPLICIT NONE
27 :
28 : TYPE(t_mpi),INTENT(IN) :: fmpi
29 :
30 :
31 : TYPE(t_input),INTENT(IN) :: input
32 : TYPE(t_noco),INTENT(IN) :: noco
33 : TYPE(t_nococonv),INTENT(IN) :: nococonv
34 : TYPE(t_sym),INTENT(IN) :: sym
35 : TYPE(t_cell),INTENT(IN) :: cell
36 : TYPE(t_kpts),INTENT(IN) :: kpts
37 : TYPE(t_atoms),INTENT(IN) :: atoms
38 : TYPE(t_enpara),INTENT(IN) :: enpara
39 : TYPE(t_potden),INTENT(IN) :: v
40 : TYPE(t_results),INTENT(IN) :: results
41 : INTEGER,INTENT(IN) :: eig_id
42 : REAL,INTENT(in) :: theta(:),phi(:) ! more than a single angle at once...
43 : REAL,INTENT(IN) :: ef(:) !Multiple Fermi energies (bandfillings)
44 : REAL,INTENT(OUT) :: seigvso(:,0:)
45 : REAL,INTENT(OUT) :: h_so(0:,:,:)
46 : ! ..
47 : ! .. Locals ..
48 : #ifdef CPP_MPI
49 : INTEGER:: ierr
50 : #endif
51 : INTEGER :: neigf=1 !not full-matrix
52 : INTEGER :: ilo,js,jsloc,nk,n,l ,lm,band,nr,ne,nat,m
53 : INTEGER :: na,nef
54 : REAL :: r1,r2
55 : COMPLEX :: c1,c2
56 :
57 0 : COMPLEX, ALLOCATABLE :: matel(:,:,:)
58 0 : REAL, ALLOCATABLE :: eig_shift(:,:,:,:)
59 0 : Real, allocatable :: w_iks(:)
60 0 : COMPLEX, ALLOCATABLE :: acof(:,:,:,:,:), bcof(:,:,:,:,:)
61 0 : COMPLEX, ALLOCATABLE :: ccof(:,:,:,:,:,:)
62 0 : COMPLEX,ALLOCATABLE :: soangl(:,:,:,:,:,:,:)
63 :
64 0 : TYPE(t_rsoc) :: rsoc
65 0 : TYPE(t_mat) :: zmat
66 0 : TYPE(t_usdus):: usdus
67 0 : TYPE(t_lapw) :: lapw
68 :
69 0 : IF (ANY(atoms%neq/=1)) CALL judft_error('(spin spiral + soc) does not work'//&
70 0 : ' properly for more than one atom per type!',calledby="ssomat")
71 :
72 :
73 :
74 : ! needed directly for calculating matrix elements
75 0 : seigvso=0.0
76 0 : ALLOCATE(eig_shift(input%neig,0:atoms%ntype,kpts%nkpt,SIZE(theta)));eig_shift=0.0
77 : ALLOCATE( acof(input%neig,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat,2,2),&
78 0 : bcof(input%neig,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat,2,2) )
79 0 : ALLOCATE( ccof(-atoms%llod:atoms%llod,input%neig,atoms%nlod,atoms%nat,2,2) )
80 :
81 0 : ALLOCATE( matel(neigf,input%neig,0:atoms%ntype) )
82 :
83 :
84 :
85 0 : CALL usdus%init(atoms,2)
86 :
87 :
88 : !Calculate radial and angular matrix elements of SOC
89 : !many directions of SOC at once...
90 0 : CALL spnorb(atoms,noco,nococonv,input,fmpi, enpara, v%mt, usdus, rsoc,.FALSE.)
91 :
92 : ALLOCATE(soangl(atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,&
93 0 : atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,SIZE(theta)))
94 0 : soangl=0.0
95 0 : DO nr=1,SIZE(theta)
96 0 : CALL spnorb_angles(atoms,fmpi,theta(nr),phi(nr),soangl(:,:,:,:,:,:,nr))
97 : ENDDO
98 :
99 0 : DO nk=fmpi%irank+1,kpts%nkpt,fmpi%isize
100 0 : CALL lapw%init(input,noco,nococonv, kpts,atoms,sym,nk,cell)
101 0 : zMat%matsize1=lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot
102 0 : zmat%matsize2=input%neig
103 0 : zmat%l_real=.FALSE.
104 0 : IF (ALLOCATED(zmat%data_c)) DEALLOCATE(zmat%data_c)
105 0 : ALLOCATE(zmat%data_c(zMat%matsize1,zmat%matsize2))
106 0 : CALL read_eig(eig_id,nk,1,neig=ne,eig=eig_shift(:,0,nk,1),zmat=zmat)
107 0 : DO jsloc= 1,2
108 0 : eig_shift(:,0,nk,1)=0.0 !not needed
109 : CALL abcof(input,atoms,sym, cell,lapw,ne,usdus,noco,nococonv,jsloc , &
110 0 : acof(:,:,:,jsloc,1),bcof(:,:,:,jsloc,1),ccof(:,:,:,:,jsloc,1),zMat)
111 : ENDDO
112 :
113 : ! rotate abcof into global spin coordinate frame
114 : nat= 0
115 0 : DO n= 1,atoms%ntype
116 0 : DO na= 1,atoms%neq(n)
117 0 : nat= nat+1
118 0 : r1= nococonv%alph(n)
119 0 : r2= nococonv%beta(n)
120 0 : DO lm= 0,atoms%lmaxd*(atoms%lmaxd+2)
121 0 : DO band= 1,input%neig
122 0 : c1= acof(band,lm,nat,1,1)
123 0 : c2= acof(band,lm,nat,2,1)
124 0 : acof(band,lm,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c1
125 0 : acof(band,lm,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX(-SIN(r2/2.),0.) *c2
126 0 : acof(band,lm,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX(+SIN(r2/2.),0.) *c1
127 0 : acof(band,lm,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c2
128 0 : c1= bcof(band,lm,nat,1,1)
129 0 : c2= bcof(band,lm,nat,2,1)
130 0 : bcof(band,lm,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c1
131 0 : bcof(band,lm,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX(-SIN(r2/2.),0.) *c2
132 0 : bcof(band,lm,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX(+SIN(r2/2.),0.) *c1
133 0 : bcof(band,lm,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c2
134 : ENDDO ! band
135 : ENDDO ! lm
136 0 : DO ilo = 1,atoms%nlo(n)
137 0 : l = atoms%llo(ilo,n)
138 0 : DO band= 1,input%neig
139 0 : DO m = -l, l
140 0 : c1= ccof(m,band,ilo,nat,1,1)
141 0 : c2= ccof(m,band,ilo,nat,2,1)
142 0 : ccof(m,band,ilo,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.))*CMPLX( COS(r2/2.),0.)*c1
143 0 : ccof(m,band,ilo,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.))*CMPLX(-SIN(r2/2.),0.)*c2
144 0 : ccof(m,band,ilo,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.))*CMPLX(+SIN(r2/2.),0.)*c1
145 0 : ccof(m,band,ilo,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.))*CMPLX( COS(r2/2.),0.)*c2
146 : ENDDO
147 : ENDDO
148 : ENDDO
149 : ENDDO
150 : ENDDO
151 0 : DO nr=1,size(theta) !loop over angles
152 : ! matrix elements within k
153 : CALL ssomatel(neigf,input,atoms, noco, &
154 : soangl(:,:,:,:,:,:,nr),rsoc%rsopp(:,:,:,:),rsoc%rsoppd(:,:,:,:),&
155 : rsoc%rsopdp(:,:,:,:),rsoc%rsopdpd(:,:,:,:),rsoc%rsoplop(:,:,:,:), &
156 : rsoc%rsoplopd(:,:,:,:),rsoc%rsopdplo(:,:,:,:),rsoc%rsopplo(:,:,:,:),&
157 : rsoc%rsoploplop(:,:,:,:,:),&
158 : .TRUE.,&
159 : acof,bcof, ccof,&
160 : acof,bcof, ccof,&
161 0 : matel )
162 0 : eig_shift(:,0:,nk,nr)=matel(1,:,0:)
163 : ENDDO
164 : ENDDO
165 :
166 : !Collect data from distributed k-loop
167 : #ifdef CPP_MPI
168 0 : IF (fmpi%irank==0) THEN
169 0 : CALL MPI_REDUCE(MPI_IN_PLACE,eig_shift,SIZE(eig_shift),MPI_DOUBLE_PRECISION,MPI_SUM,0,fmpi%mpi_comm,ierr)
170 : ELSE
171 0 : CALL MPI_REDUCE(eig_shift,eig_shift,SIZE(eig_shift),MPI_DOUBLE_PRECISION,MPI_SUM,0,fmpi%mpi_comm,ierr)
172 : ENDIF
173 : #endif
174 0 : h_so=0.0
175 0 : IF (fmpi%irank==0) THEN
176 : !Sum all shift using weights
177 0 : DO nr=1,SIZE(theta)
178 0 : DO nk=1,kpts%nkpt
179 0 : DO nef=1,size(ef)
180 0 : w_iks=kpts%wtkpt(nk)*fermifct(results%eig(:,nk,1),ef(nef),input%tkb)
181 : !for first angle, also add unmodified eigenvalue sum
182 0 : if (nr==1) seigvso(nef,0)=seigvso(nef,nr)+dot_PRODUCT(w_iks,eig_shift(:,0,nk,nr)+results%eig(:,nk,1))
183 0 : seigvso(nef,nr)=seigvso(nef,nr)+dot_PRODUCT(w_iks,eig_shift(:,0,nk,nr)+results%eig(:,nk,1))
184 0 : DO n=0,atoms%ntype
185 0 : H_so(n,nef,nr)=H_so(n,nef,nr)+dot_PRODUCT(w_iks,eig_shift(:,n,nk,nr))
186 : enddo
187 : enddo
188 : ENDDO
189 : ENDDO
190 : !seigvso= results%seigv+seigvso !now included in sum above
191 : ENDIF
192 0 : END SUBROUTINE ssomat
193 :
194 : ! ==================================================================== !
195 :
196 0 : SUBROUTINE ssomatel(neigf,input,atoms, noco,&
197 0 : soangl,rsopp,rsoppd,rsopdp,rsopdpd,rsoplop,&
198 0 : rsoplopd,rsopdplo,rsopplo,rsoploplop,&
199 : diag, &
200 0 : acof1,bcof1,ccof1,acof2,bcof2,ccof2,&
201 0 : matel )
202 : USE m_types
203 : IMPLICIT NONE
204 : TYPE(t_input),INTENT(IN) :: input
205 : TYPE(t_noco),INTENT(IN) :: noco
206 : TYPE(t_atoms),INTENT(IN) :: atoms
207 :
208 : LOGICAL, INTENT(IN) :: diag
209 : INTEGER, INTENT(IN) :: neigf
210 : REAL, INTENT(IN) :: &
211 : rsopp(:,:,:,:), rsoppd(:,:,:,:),&
212 : rsopdp(:,:,:,:), rsopdpd(:,:,:,:), &
213 : rsoplop(:,:,:,:),rsoplopd(:,:,:,:),&
214 : rsopdplo(:,:,:,:),rsopplo(:,:,:,:),&
215 : rsoploplop(:,:,:,:,:)
216 : COMPLEX, INTENT(IN) :: &
217 : soangl(:,-atoms%lmaxd:,:,:,-atoms%lmaxd:,:), &
218 : acof1(:,0:,:,:,:), &
219 : bcof1(:,0:,:,:,:),&
220 : ccof1(-atoms%llod:,:,:,:,:,:),&
221 : acof2(:,0:,:,:,:), &
222 : bcof2(:,0:,:,:,:),&
223 : ccof2(-atoms%llod:,:,:,:,:,:)
224 :
225 : Complex, INTENT(OUT) :: matel(neigf,input%neig,0:atoms%ntype)
226 :
227 : INTEGER :: band1,band2,bandf, n ,na, l,m1,m2,lm1,lm2,&
228 : jsloc1,jsloc2, js1,js2,jsnumber,ilo,ilop,nat
229 0 : COMPLEX, ALLOCATABLE :: sa(:,:),sb(:,:),sc(:,:,:),ral(:,:,:)
230 0 : COMPLEX, ALLOCATABLE :: ra(:,:),rb(:,:),rc(:,:,:),rbl(:,:,:)
231 :
232 : ! with the following nesting of loops the calculation of the
233 : ! matrix-elements is of order
234 : ! natall*lmd*neigd*(lmd+neigd) ; note that lmd+neigd << lmd*neigd
235 :
236 0 : matel(:,:,:)= CMPLX(0.,0.)
237 0 : ALLOCATE ( sa(2,0:atoms%lmaxd*(atoms%lmaxd+2)),sb(2,0:atoms%lmaxd*(atoms%lmaxd+2)),ra(2,0:atoms%lmaxd*(atoms%lmaxd+2)),rb(2,0:atoms%lmaxd*(atoms%lmaxd+2)) )
238 0 : ALLOCATE ( sc(2,-atoms%llod:atoms%llod,atoms%nlod),rc(2,-atoms%llod:atoms%llod,atoms%nlod) )
239 0 : ALLOCATE ( ral(2,-atoms%llod:atoms%llod,atoms%nlod),rbl(2,-atoms%llod:atoms%llod,atoms%nlod) )
240 :
241 : ! within one k-point loop over global spin
242 0 : IF (diag) THEN
243 : jsnumber= 2
244 : ELSE
245 0 : jsnumber= 1
246 : ENDIF
247 0 : DO js2= 1,jsnumber
248 0 : IF (diag) THEN
249 0 : js1= js2
250 : ELSE
251 : js1= 2
252 : ENDIF
253 :
254 : ! loop over MT
255 0 : na= 0
256 0 : DO n= 1,atoms%ntype
257 0 : DO nat= 1,atoms%neq(n)
258 0 : na= na+1
259 :
260 0 : DO band2= 1,input%neig ! loop over eigenstates 2
261 :
262 0 : DO l= 1,atoms%lmax(n) ! loop over l
263 0 : DO m1= -l,l ! loop over m1
264 0 : lm1= l*(l+1) + m1
265 :
266 0 : DO jsloc2= 1,2
267 0 : sa(jsloc2,lm1) = CMPLX(0.,0.)
268 0 : sb(jsloc2,lm1) = CMPLX(0.,0.)
269 0 : DO m2= -l,l
270 0 : lm2= l*(l+1) + m2
271 :
272 : sa(jsloc2,lm1)= sa(jsloc2,lm1) + &
273 : CONJG(acof2(band2,lm2,na,jsloc2,js2))&
274 0 : * soangl(l,m2,js2,l,m1,js1)
275 : sb(jsloc2,lm1)= sb(jsloc2,lm1) + &
276 : CONJG(bcof2(band2,lm2,na,jsloc2,js2))&
277 0 : * soangl(l,m2,js2,l,m1,js1)
278 :
279 : ENDDO ! m2
280 : ENDDO ! jsloc2
281 :
282 : ENDDO ! m1
283 : ENDDO ! l
284 :
285 0 : DO ilo = 1, atoms%nlo(n) ! LO-part
286 0 : l = atoms%llo(ilo,n)
287 0 : DO m1 = -l, l
288 0 : DO jsloc2= 1,2
289 0 : sc(jsloc2,m1,ilo) = CMPLX(0.,0.)
290 0 : IF (l==0) CYCLE
291 0 : DO m2= -l, l
292 : sc(jsloc2,m1,ilo) = sc(jsloc2,m1,ilo) +&
293 : CONJG(ccof2(m2,band2,ilo,na,jsloc2,js2))&
294 0 : * soangl(l,m2,js2,l,m1,js1)
295 : ENDDO
296 : ENDDO
297 : ENDDO
298 : ENDDO ! ilo
299 :
300 0 : DO l= 1,atoms%lmax(n) ! loop over l
301 0 : DO m1= -l,l ! loop over m1
302 0 : lm1= l*(l+1) + m1
303 :
304 0 : DO jsloc1= 1,2
305 0 : ra(jsloc1,lm1)= CMPLX(0.,0.)
306 0 : rb(jsloc1,lm1)= CMPLX(0.,0.)
307 0 : DO jsloc2= 1,2
308 : ra(jsloc1,lm1)= ra(jsloc1,lm1) + &
309 : sa(jsloc2,lm1) * rsopp(n,l,jsloc1,jsloc2) &
310 0 : + sb(jsloc2,lm1) * rsoppd(n,l,jsloc1,jsloc2)
311 : rb(jsloc1,lm1)= rb(jsloc1,lm1) +&
312 : sa(jsloc2,lm1) * rsopdp(n,l,jsloc1,jsloc2)&
313 0 : + sb(jsloc2,lm1) * rsopdpd(n,l,jsloc1,jsloc2)
314 : ENDDO ! jsloc2
315 : ENDDO ! jsloc1
316 :
317 : ENDDO ! m1
318 : ENDDO ! l
319 :
320 0 : DO ilo = 1, atoms%nlo(n) ! LO-part
321 0 : l = atoms%llo(ilo,n)
322 0 : DO m1 = -l, l
323 0 : lm1= l*(l+1) + m1
324 0 : DO jsloc1= 1,2
325 0 : ral(jsloc1,m1,ilo) = CMPLX(0.,0.)
326 0 : rbl(jsloc1,m1,ilo) = CMPLX(0.,0.)
327 0 : rc(jsloc1,m1,ilo) = CMPLX(0.,0.)
328 0 : DO jsloc2= 1,2
329 : ral(jsloc1,m1,ilo) = ral(jsloc1,m1,ilo) +&
330 0 : sc(jsloc2,m1,ilo) * rsopplo(n,ilo,jsloc1,jsloc2)
331 : rbl(jsloc1,m1,ilo) = rbl(jsloc1,m1,ilo) +&
332 0 : sc(jsloc2,m1,ilo) * rsopdplo(n,ilo,jsloc1,jsloc2)
333 : rc(jsloc1,m1,ilo) = rc(jsloc1,m1,ilo) +&
334 : sa(jsloc2,lm1) * rsoplop(n,ilo,jsloc1,jsloc2)&
335 0 : + sb(jsloc2,lm1) * rsoplopd(n,ilo,jsloc1,jsloc2)
336 : ENDDO
337 : ENDDO
338 : ENDDO
339 : ENDDO ! ilo
340 :
341 0 : DO l= 1,atoms%lmax(n) ! loop over l
342 0 : DO m1= -l,l ! loop over m1
343 0 : lm1= l*(l+1) + m1
344 :
345 0 : DO jsloc1= 1,2
346 0 : DO bandf= 1,neigf
347 0 : IF (neigf==input%neig) THEN
348 : band1= bandf
349 : ELSE
350 0 : band1= band2
351 : ENDIF
352 : matel(bandf,band2,n)= matel(bandf,band2,n) +&
353 : acof1(band1,lm1,na,jsloc1,js1)*ra(jsloc1,lm1) &
354 0 : + bcof1(band1,lm1,na,jsloc1,js1)*rb(jsloc1,lm1)
355 : ENDDO ! band1
356 : ENDDO ! jsloc1
357 :
358 : ENDDO ! m1,lm1
359 : ENDDO ! l
360 :
361 0 : DO ilo = 1, atoms%nlo(n) ! LO-part
362 0 : l = atoms%llo(ilo,n)
363 0 : IF (l==0) CYCLE
364 0 : DO m1 = -l, l
365 0 : lm1= l*(l+1) + m1
366 :
367 0 : DO jsloc1= 1,2
368 0 : DO bandf= 1,neigf
369 0 : IF (neigf==input%neig) THEN
370 : band1= bandf
371 : ELSE
372 0 : band1= band2
373 : ENDIF
374 : matel(bandf,band2,n)= matel(bandf,band2,n) +&
375 : ccof1(m1,band1,ilo,na,jsloc1,js1)*rc(jsloc1,m1,ilo)&
376 : + acof1(band1,lm1,na,jsloc1,js1)*ral(jsloc1,m1,ilo)&
377 0 : + bcof1(band1,lm1,na,jsloc1,js1)*rbl(jsloc1,m1,ilo)
378 : ENDDO ! band1
379 : ENDDO ! jsloc1
380 :
381 0 : DO ilop = 1,atoms%nlo(n)
382 0 : IF (atoms%llo(ilop,n).EQ.l) THEN
383 0 : DO jsloc1= 1,2
384 0 : DO bandf= 1,neigf
385 0 : IF (neigf==input%neig) THEN
386 : band1= bandf
387 : ELSE
388 0 : band1= band2
389 : ENDIF
390 0 : DO jsloc2= 1,2
391 : matel(bandf,band2,n)= matel(bandf,band2,n) +&
392 : ccof1(m1,band1,ilo,na,jsloc1,js1)*&
393 : rsoploplop(n,ilo,ilop,jsloc1,jsloc2)*&
394 0 : sc(jsloc2,m1,ilop)
395 : ENDDO ! jsloc2
396 : ENDDO ! band1
397 : ENDDO ! jsloc1
398 : ENDIF
399 : ENDDO ! ilop
400 :
401 : ENDDO ! m1
402 : ENDDO ! ilo
403 :
404 : ENDDO ! band2
405 : ENDDO ! nat,na
406 : ENDDO ! n
407 : ENDDO ! js2,js1
408 :
409 0 : DO n= 1,atoms%ntype
410 0 : DO band2= 1,input%neig
411 0 : DO bandf= 1,neigf
412 0 : matel(bandf,band2,0)= matel(bandf,band2,0) + matel(bandf,band2,n)
413 : ENDDO
414 : ENDDO
415 : ENDDO
416 :
417 0 : IF (diag) THEN
418 0 : DO n= 1,atoms%ntype
419 0 : DO band2= 1,input%neig
420 0 : IF (neigf==input%neig) THEN
421 0 : bandf= band2
422 : ELSE
423 0 : bandf= 1
424 : ENDIF
425 0 : IF (ABS(AIMAG(matel(bandf,band2,n)))>1.e-12) THEN
426 0 : PRINT *,bandf,band2,n,AIMAG(matel(bandf,band2,n))
427 0 : CALL judft_error('Stop in ssomatel: diagonal matrix element not real')
428 : ENDIF
429 : ENDDO
430 : ENDDO
431 : ENDIF
432 :
433 0 : DEALLOCATE ( sa,sb,ra,rb )
434 :
435 0 : END SUBROUTINE ssomatel
436 : END MODULE m_ssomat
|
