Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_wannier
8 : USE m_juDFT
9 : #ifdef CPP_MPI
10 : use mpi
11 : #endif
12 : CONTAINS
13 4 : SUBROUTINE wannier(&
14 : fmpi,input,kpts,sym,atoms,stars,vacuum,sphhar ,&
15 : wann,noco,nococonv,cell,enpara,banddos,sliceplot,vTot,results,&
16 4 : eig_idList,l_real,nkpt)
17 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
18 : ! Makes necessary for the construction of the wannier functions
19 : ! (W90: Yates, Mostofi, Marzari, Souza, Vanderbilt '06 f90 code)
20 : ! ab initio preliminaries: constructs the overlaps of the periodic
21 : ! parts of the wavefunctions and the projections of the
22 : ! wavefunctions
23 : ! onto a set of starting wfs, i.e. atomic-like orbitals.
24 : ! YM 06
25 : ! Mmn(k,b) = <u_{nk}|u_{m(k+b)}>, u being a periodic part
26 : ! of the wavefunction psi_nk
27 : ! A_mn^k = <psi_mk|g_n>, where g_n is a trial orbital
28 : ! which are written into the files 'WF1.mmn' and 'WF1.amn'
29 : ! Marzari Vanderbilt PRB 56,12847(1997)
30 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
31 : ! Parallelization, Optionals, Symmetry, Noco&Soc:
32 : ! Frank Freimuth
33 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
34 : ! The routine reads the bkpts file, which contains the following
35 : ! information:
36 : ! 1st line: nntot (INT) - number of the nearest neighbors for
37 : ! each k-point in the MP mesh
38 : ! 2-nkpts*nntot lines containing 5 integers i1,i2,i3,i4,i5:
39 : ! i1 - the number of the k-point in the kpts file
40 : ! i2 - number of the k-point, which is a periodic image of
41 : ! k+b in the 1st BZ
42 : ! i3-i5 - coordinates of the G-vector, connecting k-point
43 : ! i2 with the actual k+b k-point
44 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
45 : ! In general, the number of bands for each k-poin t is
46 : ! different: N(k), and also differs from the number of bands
47 : ! we are interested in: N (for instance 5 d-bands of Cu among
48 : ! the 6 s- and d-bands). While matrices Mmn(k) are square
49 : ! for each k-point, matrices Mmn(k,b) can be made so after
50 : ! defining the maximum number of bands max(N(k)).
51 : ! The matrix Amn is non-diagonal by default (N(k)*N).
52 : !ccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
53 : ! Total number of wannier functions: nwfs
54 : ! sliceplot%e1s,sliceplot%e2s: lower and upper boundaries of the energy window:
55 : ! Needed for sorting by number and sorting by energy.
56 : ! Not needed for sorting by index.
57 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
58 : ! Extension to case of higher-dimensional Wannier functions
59 : ! according to the formalism in PRB 91, 184413 (2015)
60 : ! Jan-Philipp Hanke
61 : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
62 : USE m_types
63 : USE m_constants
64 : USE m_wann_mmnk_symm
65 : USE m_wann_rw_eig
66 : USE m_abcof
67 : USE m_radfun
68 : USE m_radflo
69 : USE m_cdnread
70 : USE m_wann_mmk0_vac
71 : USE m_wann_mmkb_vac
72 : USE m_wann_updown
73 : USE m_wann_mmk0_sph
74 : USE m_wann_ujugaunt
75 : USE m_wann_mmkb_sph
76 : USE m_wann_projmethod
77 : USE m_wann_amn
78 : USE m_wann_abinv
79 : USE m_wann_kptsrotate
80 : USE m_wann_plot
81 : ! USE m_wann_read_inp !Call wann_read_inp in fleur_init
82 : USE m_wann_plot_symm
83 : USE m_wann_mmkb_int
84 : USE m_wann_postproc
85 : USE m_wann_write_mmnk
86 : USE m_wann_write_amn
87 : USE m_wann_write_nabla
88 : USE m_vsoc
89 : USE m_wann_write_matrix4
90 : USE m_wann_write_matrix5
91 : USE m_wann_orbcomp
92 : USE m_wann_anglmom
93 : #ifdef CPP_TOPO
94 : USE m_wann_surfcurr
95 : USE m_wann_surfcurr_int2
96 : USE m_wann_nabla
97 : USE m_wann_nabla_vac
98 : USE m_wann_soc_to_mom
99 : #endif
100 : USE m_wann_gwf_tools, ONLY : get_index_kq, gwf_plottemplate
101 : USE m_wann_gwf_commat
102 : USE m_wann_gwf_anglmom
103 : USE m_wann_write_mmnk2
104 : USE m_wann_uHu
105 : USE m_wann_uHu_dmi
106 : USE m_eig66_io
107 :
108 : IMPLICIT NONE
109 :
110 : #ifdef CPP_MPI
111 : INTEGER ierr(3)
112 : INTEGER cpu_index
113 : INTEGER stt(MPI_STATUS_SIZE)
114 : #endif
115 :
116 :
117 : TYPE(t_mpi), INTENT(IN) :: fmpi
118 : TYPE(t_input), INTENT(IN) :: input
119 : TYPE(t_kpts), INTENT(IN) :: kpts
120 : TYPE(t_sym), INTENT(IN) :: sym
121 : TYPE(t_atoms), INTENT(IN) :: atoms
122 : TYPE(t_stars), INTENT(IN) :: stars
123 : TYPE(t_vacuum), INTENT(IN) :: vacuum
124 : TYPE(t_sphhar), INTENT(IN) :: sphhar
125 :
126 : TYPE(t_noco), INTENT(IN) :: noco
127 : TYPE(t_nococonv), INTENT(IN) :: nococonv
128 : TYPE(t_cell), INTENT(IN) :: cell
129 : TYPE(t_enpara), INTENT(IN) :: enpara
130 : TYPE(t_banddos), INTENT(IN) :: banddos
131 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
132 : TYPE(t_potden), INTENT(IN) :: vTot
133 : TYPE(t_results), INTENT(IN) :: results
134 : TYPE(t_wann), INTENT(INOUT) :: wann
135 :
136 : LOGICAL, INTENT (in) :: l_real
137 : INTEGER, INTENT (in) :: nkpt
138 : INTEGER, INTENT (IN) :: eig_idList(:)!wann%nparampts)
139 :
140 : !ccccccccccccccccc local variables cccccccccccccccccccc
141 : INTEGER :: lmd,n,nmat,iter,ikpt,ikpt_b, pc
142 : INTEGER :: addnoco,loplod,addnoco2,igvm2,eig_id
143 : INTEGER :: noccbd,noccbd_b,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
144 : INTEGER :: jsp_start,jsp_end,nrec,nrec1,nrec_b,nbands,nbands_b
145 : INTEGER :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
146 : INTEGER :: i1,i2,i3,in,lda
147 4 : INTEGER :: n_bands(0:input%neig),nslibd,nslibd_b
148 : CHARACTER(len=8) :: dop,iop,name(10)
149 : REAL :: wronk,phase
150 : COMPLEX :: c_phase
151 4 : REAL :: eig(input%neig),eig_b(input%neig)
152 : REAL :: efermi
153 : LOGICAL :: l_p0,l_bkpts,l_proj,l_amn,l_mmn
154 : !!! energy window boundaries
155 4 : INTEGER, ALLOCATABLE :: innerEig_idList(:)
156 4 : REAL, ALLOCATABLE :: we(:),we_b(:)
157 :
158 4 : REAL, ALLOCATABLE :: eigg(:)
159 4 : REAL kpoints(nkpt)
160 : !!! a and b coeff. constructed for each k-point
161 4 : COMPLEX, ALLOCATABLE :: acof(:,:,:),acof_b(:,:,:)
162 4 : COMPLEX, ALLOCATABLE :: bcof(:,:,:),bcof_b(:,:,:)
163 4 : COMPLEX, ALLOCATABLE :: ccof(:,:,:,:),ccof_b(:,:,:,:)
164 : !!! the parameters for the number of wfs
165 : INTEGER :: nwfs
166 : !!! the potential in the spheres and the vacuum
167 4 : REAL, ALLOCATABLE :: vr(:,:,:),vz(:,:,:)
168 : !!! auxiliary potentials
169 4 : COMPLEX, ALLOCATABLE :: vpw(:,:)
170 : !!! bkpts data
171 : INTEGER nntot,ikpt_help
172 4 : INTEGER, ALLOCATABLE :: gb(:,:,:),bpt(:,:)
173 : !!! radial wavefunctions in the muffin-tins and more ...
174 4 : REAL, ALLOCATABLE :: flo(:,:,:,:,:),vso(:,:,:)
175 4 : REAL, ALLOCATABLE :: ff(:,:,:,:,:),gg(:,:,:,:,:)
176 :
177 4 : REAL :: uuilon(atoms%nlod,atoms%ntype)
178 4 : REAL :: duilon(atoms%nlod,atoms%ntype)
179 4 : REAL :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
180 : !!! the Mmn matrices
181 4 : COMPLEX, ALLOCATABLE :: mmnk(:,:,:,:),mmn(:,:,:)
182 4 : COMPLEX, ALLOCATABLE :: amn(:,:,:),nablamat(:,:,:,:)
183 4 : COMPLEX, ALLOCATABLE :: soctomom(:,:,:,:)
184 4 : COMPLEX, ALLOCATABLE :: surfcurr(:,:,:,:)
185 4 : COMPLEX, ALLOCATABLE :: socmmn(:,:,:)
186 : COMPLEX, ALLOCATABLE :: a(:)
187 4 : COMPLEX, ALLOCATABLE :: psiw(:,:,:)
188 4 : COMPLEX, ALLOCATABLE :: anglmom(:,:,:,:)
189 4 : COMPLEX, ALLOCATABLE :: orbcomp(:,:,:,:,:)
190 : !..wf-hamiltonian in real space (hopping in the same unit cell)
191 4 : COMPLEX, ALLOCATABLE :: hwfr(:,:),hwfr2(:,:)
192 : ! real, allocatable :: ei(:)
193 : COMPLEX, ALLOCATABLE :: work(:)
194 : REAL,ALLOCATABLE::centers(:,:,:)
195 : LOGICAL :: l_file
196 : LOGICAL :: l_amn2, l_conjugate
197 : CHARACTER(len=3) :: spin12(2)
198 : DATA spin12/'WF1' , 'WF2'/
199 : CHARACTER(len=30) :: task
200 4 : INTEGER,ALLOCATABLE::irreduc(:)
201 4 : INTEGER,ALLOCATABLE::mapkoper(:)
202 : INTEGER :: fullnkpts,kpt,kptibz,kptibz_b,j1,j2,j3,oper,oper_b,k
203 : REAL :: bkrot(3),dirfacs(3)
204 : INTEGER :: ios,kplot,kplotoper,plotoper,gfcut
205 : COMPLEX :: phasust
206 4 : INTEGER,ALLOCATABLE::pair_to_do(:,:)
207 4 : INTEGER :: ngopr1(atoms%nat)
208 4 : INTEGER,ALLOCATABLE::maptopair(:,:,:)
209 : INTEGER :: wannierspin,jspin2,jspin7,jspin2_b
210 : REAL, ALLOCATABLE :: rwork(:)
211 4 : REAL,ALLOCATABLE::kdiff(:,:)
212 4 : INTEGER,ALLOCATABLE :: shiftkpt(:,:)
213 : INTEGER :: unigrid(6),gfthick
214 4 : COMPLEX,ALLOCATABLE::ujug(:,:,:,:),ujdg(:,:,:,:)
215 4 : COMPLEX,ALLOCATABLE::djug(:,:,:,:),djdg(:,:,:,:)
216 4 : COMPLEX,ALLOCATABLE::ujulog(:,:,:,:,:)
217 4 : COMPLEX,ALLOCATABLE::djulog(:,:,:,:,:)
218 4 : COMPLEX,ALLOCATABLE::ulojug(:,:,:,:,:)
219 4 : COMPLEX,ALLOCATABLE::ulojdg(:,:,:,:,:)
220 4 : COMPLEX,ALLOCATABLE::ulojulog(:,:,:,:,:,:)
221 : INTEGER :: n_start,n_end,mlotot,mlolotot,err
222 : INTEGER :: mlot_d,mlolot_d,ilo,dir,length
223 : CHARACTER(len=2) :: spin012(0:2)
224 : DATA spin012/' ', '.1', '.2'/
225 : CHARACTER(len=6) :: filename
226 : REAL :: arg,hescale
227 : COMPLEX :: nsfactor,nsfactor_b,VALUE
228 : REAL :: b1(3),b2(3)
229 : REAL,PARAMETER :: bohrtocm=0.529177e-8
230 : REAL,PARAMETER :: condquant=7.7480917e-5
231 : INTEGER :: npotmatfile,ig3,maxvac,irec,imz,ivac,ipot
232 : LOGICAL :: l_orbcompinp
233 : INTEGER :: num_angl,nbasfcn,nbasfcn_b,noconbasfcn,noconbasfcn_b
234 4 : COMPLEX,ALLOCATABLE :: vxy(:,:,:)
235 :
236 :
237 : !---->gwf
238 :
239 : ! FURTHER VARIABLES
240 : REAL :: qpt_i(3),qptb_i(3)
241 4 : REAL :: alph_i(atoms%ntype),alphb_i(atoms%ntype)
242 4 : REAL :: beta_i(atoms%ntype),betab_i(atoms%ntype)
243 : REAL :: theta_i, thetab_i, phi_i, phib_i
244 : REAL :: dalph,db1,db2,coph,siph
245 4 : REAL :: zero_taual(3,atoms%nat),bqpt(3)
246 4 : REAL :: eig_qb(input%neig)
247 :
248 8 : REAL,ALLOCATABLE :: qdiff(:,:), we_qb(:)
249 : REAL,ALLOCATABLE :: energies(:,:,:)
250 4 : REAL,ALLOCATABLE :: zero_qdiff(:,:)
251 :
252 :
253 4 : INTEGER,ALLOCATABLE :: irreduc_q(:),mapqoper(:)
254 4 : INTEGER,ALLOCATABLE :: shiftqpt(:,:),pair_to_do_q(:,:)
255 4 : INTEGER,ALLOCATABLE :: maptopair_q(:,:,:)
256 4 : INTEGER,ALLOCATABLE :: gb_q(:,:,:),bpt_q(:,:)
257 :
258 : INTEGER :: nntot_q = 1
259 : INTEGER :: fullnqpts = 1
260 : INTEGER :: funit_start = 5000
261 : INTEGER :: qptibz, qptibz_b, oper_q, oper_qb
262 : INTEGER :: qpt,iqpt_help, iqpt, iqpt_b
263 : INTEGER :: nbands_qb, nmat_qb, nslibd_qb, noccbd_qb
264 : INTEGER :: sign_q = 1,band_help
265 : INTEGER :: doublespin,jspin_b,jspin3,jspin4,jspin5
266 : INTEGER :: doublespin_max,nrec5
267 : INTEGER :: count_i,count_j
268 : INTEGER :: n1,n2,ii,jj,lmplmd
269 :
270 : COMPLEX :: interchi,vacchi,amnchi
271 : COMPLEX :: phasfac,phasfac2
272 :
273 4 : COMPLEX,ALLOCATABLE :: chi(:)
274 4 : COMPLEX,ALLOCATABLE :: acof_qb(:,:,:)
275 4 : COMPLEX,ALLOCATABLE :: bcof_qb(:,:,:)
276 4 : COMPLEX,ALLOCATABLE :: ccof_qb(:,:,:,:)
277 4 : COMPLEX,ALLOCATABLE :: mmnk_q(:,:,:,:)
278 4 : COMPLEX,ALLOCATABLE :: m_int(:,:,:,:)
279 4 : COMPLEX,ALLOCATABLE :: m_sph(:,:,:,:)
280 4 : COMPLEX,ALLOCATABLE :: m_vac(:,:,:,:)
281 4 : COMPLEX,ALLOCATABLE :: ujug_q(:,:,:,:),ujdg_q(:,:,:,:)
282 4 : COMPLEX,ALLOCATABLE :: djug_q(:,:,:,:),djdg_q(:,:,:,:)
283 4 : COMPLEX,ALLOCATABLE :: ujulog_q(:,:,:,:,:)
284 4 : COMPLEX,ALLOCATABLE :: djulog_q(:,:,:,:,:)
285 4 : COMPLEX,ALLOCATABLE :: ulojug_q(:,:,:,:,:)
286 4 : COMPLEX,ALLOCATABLE :: ulojdg_q(:,:,:,:,:)
287 4 : COMPLEX,ALLOCATABLE :: ulojulog_q(:,:,:,:,:,:)
288 :
289 : CHARACTER(len=30) fname
290 :
291 : LOGICAL :: l_bqpts,l_gwf,l_nochi
292 :
293 4 : TYPE(t_usdus) :: usdus
294 4 : TYPE(t_mat) :: zMat, zzMat, zMat_b, zMat_qb
295 4 : TYPE(t_lapw) :: lapw, lapw_b, lapw_qb
296 68 : TYPE(t_wann) :: wannTemp
297 :
298 4 : eig_id = eig_idList(1)
299 :
300 : !----<gwf
301 :
302 4 : lmplmd = (atoms%lmaxd*(atoms%lmaxd+2)* (atoms%lmaxd*(atoms%lmaxd+2)+3))/2
303 :
304 :
305 : !-----initializations
306 12 : ngopr1(:)=1
307 36 : zero_taual = 0.0
308 :
309 4 : hescale=sqrt(tpi_const*condquant/bohrtocm/cell%omtil)
310 :
311 4 : CALL timestart("Wannier total")
312 :
313 4 : l_p0 = .FALSE.
314 4 : IF (fmpi%irank.EQ.0) l_p0 = .TRUE.
315 :
316 4 : lmd = atoms%lmaxd*(atoms%lmaxd+2)
317 :
318 : !!! should be changed in case the windows are really used
319 4 : nkpts = nkpt
320 :
321 : ! do we have to construct GWF ?
322 : l_gwf = .FALSE.
323 4 : l_gwf = wann%l_sgwf.OR.wann%l_socgwf
324 :
325 4 : l_nochi = .FALSE.
326 4 : INQUIRE(file='l_nochi',exist=l_nochi)
327 4 : IF(l_gwf.AND.l_p0) WRITE(*,*)'disable chi trafo: ',l_nochi
328 :
329 4 : IF(l_gwf.AND.l_p0) CALL gwf_plottemplate()
330 4 : ALLOCATE( chi(atoms%ntype) )
331 :
332 : !-----read the input file to determine what to do
333 4 : wann%atomlist_num=atoms%nat
334 4 : wann%oc_num_orbs=atoms%nat
335 :
336 : ! CALL wann_read_inp(input,l_p0,wann) !Call wann_read_inp in fleur_init
337 :
338 :
339 : !-----input file for orbital decomposition
340 4 : IF(wann%l_orbcomp.OR.wann%l_orbcomprs)THEN
341 0 : INQUIRE(file='orbcomp_inp',exist=l_orbcompinp)
342 0 : IF(l_orbcompinp)THEN
343 0 : OPEN(159,file='orbcomp_inp')
344 0 : READ(159,*)wann%oc_num_orbs,wann%l_oc_f
345 0 : ALLOCATE(wann%oc_orbs(wann%oc_num_orbs))
346 0 : DO n=1,wann%oc_num_orbs
347 0 : READ(159,*)wann%oc_orbs(n)
348 : ENDDO
349 0 : CLOSE(159)
350 : ELSE !default is all atoms including f
351 0 : wann%oc_num_orbs=atoms%nat
352 0 : wann%l_oc_f=.TRUE.
353 0 : ALLOCATE(wann%oc_orbs(wann%oc_num_orbs))
354 0 : DO n=1,wann%oc_num_orbs
355 4 : wann%oc_orbs(n)=n
356 : ENDDO
357 : ENDIF
358 : ENDIF
359 :
360 4 : IF(wann%l_updown)THEN
361 : CALL wann_updown(&
362 : wann,fmpi,input,kpts,sym,atoms,stars,vacuum,sphhar, &
363 : noco,nococonv,cell,vTot,&
364 : enpara,eig_idList(1),l_real,&
365 : fmpi%mpi_comm,atoms%l_dulo,noco%l_noco,noco%l_ss,&
366 : atoms%lmaxd,atoms%ntype,input%neig,atoms%nat,sym%nop,&
367 : lapw%dim_nvd(),input%jspins,atoms%llod,&
368 : atoms%nlod,atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
369 : atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
370 : atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
371 : atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
372 : nococonv%beta,nococonv%qss,& ! TODO: adapt if needed&
373 : stars%sk2,stars%phi2 ,fmpi%irank,fmpi%isize,&
374 : stars%ng3,&
375 : vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,sphhar%nlhd,stars%ng3,&
376 : vacuum%nvac,sym%invs,sym%invs2,input%film,sphhar%nlh,&
377 : atoms%jri,sphhar%ntypsd,sym%ntypsy,input%jspins,nkpt,&
378 : atoms%dx,stars%ng2,atoms%rmsh,sliceplot%e1s,sliceplot%e2s,&
379 : atoms%ulo_der,stars%ustep,stars%ig,stars%mx1,&
380 : stars%mx2,stars%mx3,stars%rgphs,&
381 : sliceplot%slice,sliceplot%kk,sliceplot%nnne,cell%z1,&
382 : lapw%dim_nv2d(),vacuum%nmzxy,vacuum%nmz,vacuum%delz,&
383 : stars%ig2,cell%area,sym%tau,atoms%zatom,stars%ng2,sym%nop2,&
384 : cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
385 : sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
386 : sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
387 0 : nococonv%theta,nococonv%phi)
388 0 : if(wann%l_stopupdown)then
389 0 : DO pc = 1, wann%nparampts
390 0 : CALL close_eig(eig_idList(pc))
391 : END DO
392 :
393 0 : CALL juDFT_end("updown done",fmpi%irank)
394 : endif
395 : ENDIF
396 :
397 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
398 : ! modern theory of orbital magnetization from Wannier functions
399 : ! Jan-Philipp Hanke
400 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
401 4 : IF(wann%l_matrixuHu)THEN
402 0 : wannTemp = wann
403 : CALL wann_uHu(&
404 : stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
405 : banddos ,noco,nococonv,cell,vTot,wannTemp,enpara,eig_idList,&
406 : l_real,atoms%l_dulo,noco%l_noco,noco%l_ss,atoms%lmaxd,&
407 : atoms%ntype,input%neig,atoms%nat,sym%nop,lapw%dim_nvd(),&
408 : input%jspins,atoms%llod,atoms%nlod,&
409 : atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
410 : atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
411 : atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
412 : atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
413 : nococonv%beta,nococonv%qss,stars%sk2,stars%phi2,&
414 : fmpi%irank,&
415 : fmpi%isize,stars%ng3,vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,&
416 : sphhar%nlhd,stars%ng3,vacuum%nvac,sym%invs,sym%invs2,&
417 : input%film,sphhar%nlh,atoms%jri,sphhar%ntypsd,sym%ntypsy,&
418 : input%jspins,nkpt,atoms%dx,stars%ng2,atoms%rmsh,&
419 : sliceplot%e1s,sliceplot%e2s,atoms%ulo_der,stars%ustep,&
420 : stars%ig,stars%mx1,stars%mx2,stars%mx3,&
421 : stars%rgphs,sliceplot%slice,&
422 : sliceplot%kk,sliceplot%nnne,cell%z1,lapw%dim_nv2d(),&
423 : vacuum%nmzxy,vacuum%nmz,vacuum%delz,stars%ig2,&
424 : cell%area,sym%tau,atoms%zatom,stars%ng2,stars%kv2,sym%nop2,&
425 : cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
426 : sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
427 : sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
428 : nococonv%theta,nococonv%phi,&
429 : wann%l_ms,wann%l_sgwf,wann%l_socgwf,wann%aux_latt_const,&
430 : wann%param_file,wann%param_vec,wann%nparampts,&
431 0 : wann%param_alpha,wann%l_dim)
432 :
433 :
434 0 : if(wann%l_stopuhu) then
435 :
436 0 : DO pc = 1, wann%nparampts
437 0 : CALL close_eig(eig_idList(pc))
438 : END DO
439 :
440 0 : CALL juDFT_end("wann_uHu done",fmpi%irank)
441 : endif
442 : ENDIF
443 :
444 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
445 : ! modern theory of DMI from higher-dimensional Wannier functions
446 : ! Jan-Philipp Hanke
447 : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
448 4 : IF(wann%l_matrixuHu_dmi)THEN
449 0 : wannTemp = wann
450 : CALL wann_uHu_dmi(&
451 : stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
452 : banddos ,noco,nococonv,cell,vTot,wannTemp,enpara,eig_idList,&
453 : l_real,atoms%l_dulo,noco%l_noco,noco%l_ss,atoms%lmaxd,&
454 : atoms%ntype,input%neig,atoms%nat,sym%nop,lapw%dim_nvd(),&
455 : input%jspins,lapw%dim_nbasfcn(),atoms%llod,atoms%nlod,&
456 : atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
457 : atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
458 : atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
459 : atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
460 : nococonv%beta,nococonv%qss,stars%sk2,stars%phi2,&
461 : fmpi%irank,&
462 : fmpi%isize,stars%ng3,vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,&
463 : sphhar%nlhd,stars%ng3,vacuum%nvac,sym%invs,sym%invs2,&
464 : input%film,sphhar%nlh,atoms%jri,sphhar%ntypsd,sym%ntypsy,&
465 : input%jspins,nkpt,atoms%dx,stars%ng2,atoms%rmsh,&
466 : sliceplot%e1s,sliceplot%e2s,atoms%ulo_der,stars%ustep,&
467 : stars%ig,stars%mx1,stars%mx2,stars%mx3,&
468 : stars%rgphs,sliceplot%slice,&
469 : sliceplot%kk,sliceplot%nnne,cell%z1,lapw%dim_nv2d(),&
470 : vacuum%nmzxy,vacuum%nmz,vacuum%delz,stars%ig2,&
471 : cell%area,sym%tau,atoms%zatom,stars%ng2,stars%kv2,sym%nop2,&
472 : cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
473 : sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
474 : sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
475 : nococonv%theta,nococonv%phi,&
476 : wann%l_ms,wann%l_sgwf,wann%l_socgwf,wann%aux_latt_const,&
477 : wann%param_file,wann%param_vec,wann%nparampts,&
478 0 : wann%param_alpha,wann%l_dim,l_nochi)
479 :
480 0 : DO pc = 1, wann%nparampts
481 0 : CALL close_eig(eig_idList(pc))
482 : END DO
483 :
484 0 : CALL juDFT_end("wann_uHu dmi done",fmpi%irank)
485 : ENDIF
486 :
487 4 : IF(wann%l_byenergy.AND.wann%l_byindex) CALL juDFT_error&
488 0 : ("byenergy.and.byindex",calledby ="wannier")
489 4 : IF(wann%l_byenergy.AND.wann%l_bynumber) CALL juDFT_error&
490 0 : ("byenergy.and.bynumber",calledby ="wannier")
491 4 : IF(wann%l_bynumber.AND.wann%l_byindex) CALL juDFT_error&
492 0 : ("bynumber.and.byindex",calledby ="wannier")
493 4 : IF(.NOT.(wann%l_bynumber.OR.wann%l_byindex.OR.wann%l_byenergy))&
494 0 : CALL juDFT_error("no rule to sort bands",calledby ="wannier")
495 :
496 :
497 4 : efermi=results%ef
498 4 : IF(.NOT.wann%l_fermi)efermi=0.0
499 :
500 : #ifdef CPP_MPI
501 4 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
502 : #endif
503 :
504 : !**************************************************************
505 : ! for bzsym=.true.: determine mapping between kpts and w90kpts
506 : !**************************************************************
507 4 : IF (wann%l_bzsym) THEN
508 0 : l_file=.FALSE.
509 0 : INQUIRE(file='w90kpts',exist=l_file)
510 0 : IF(.NOT.l_file) CALL juDFT_error&
511 0 : ("w90kpts not found, needed if bzsym",calledby ="wannier")
512 0 : OPEN(412,file='w90kpts',form='formatted')
513 0 : READ(412,*)fullnkpts
514 0 : CLOSE(412)
515 0 : IF(l_p0)PRINT*,"fullnkpts=",fullnkpts
516 0 : IF(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"&
517 0 : ,calledby ="wannier")
518 0 : ALLOCATE(irreduc(fullnkpts),mapkoper(fullnkpts))
519 0 : ALLOCATE(shiftkpt(3,fullnkpts))
520 0 : l_file=.FALSE.
521 0 : INQUIRE(file='kptsmap',exist=l_file)
522 0 : IF(.NOT.l_file) CALL juDFT_error&
523 0 : ("kptsmap not found, needed if bzsym",calledby ="wannier")
524 0 : OPEN(713,file='kptsmap')
525 0 : DO i=1,fullnkpts
526 0 : READ(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
527 0 : IF(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby ="wannier")
528 0 : IF(l_p0)PRINT*,i,irreduc(i),mapkoper(i)
529 : ENDDO
530 0 : CLOSE(713)
531 0 : IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error&
532 0 : ("max(irreduc(:))/=nkpts",calledby ="wannier")
533 : ELSE
534 4 : fullnkpts=nkpts
535 4 : ALLOCATE(irreduc(fullnkpts),mapkoper(fullnkpts))
536 4 : ALLOCATE(shiftkpt(3,fullnkpts))
537 : ENDIF
538 :
539 :
540 4 : IF(l_gwf) fullnqpts = wann%nparampts
541 :
542 :
543 4 : nrec = 0
544 4 : IF(l_p0)THEN
545 2 : WRITE (*,*) 'fermi energy:',efermi
546 2 : WRITE (*,*) 'emin,emax=',sliceplot%e1s,sliceplot%e2s
547 : ENDIF
548 :
549 : IF((.NOT.wann%l_matrixmmn).AND.(.NOT.wann%l_wann_plot).AND.&
550 : (.NOT.wann%l_matrixamn).AND.(.NOT.wann%l_projmethod).AND.&
551 : (.NOT.wann%l_bestproj).AND.(.NOT.wann%l_nabla).AND.&
552 : (.NOT.wann%l_mmn0).AND.(.NOT.wann%l_surfcurr).AND.&
553 : (.NOT.wann%l_offdiposop).AND.(.NOT.wann%l_anglmom).AND.&
554 4 : (.NOT.wann%l_orbcomp).AND.(.NOT.wann%l_perturb) .AND.&
555 : (.NOT.wann%l_finishgwf) ) GOTO 1911
556 :
557 : !**********************************************************
558 : !cccccccccccccc read in the bkpts file ccccccccccccccccc
559 : !**********************************************************
560 2 : IF (wann%l_matrixmmn) THEN ! for Omega functional minimization
561 2 : l_bkpts = .FALSE.
562 2 : INQUIRE (file='bkpts',exist=l_bkpts)
563 2 : IF (.NOT.l_bkpts) CALL juDFT_error("need bkpts for matrixmmn"&
564 0 : ,calledby ="wannier")
565 2 : OPEN (202,file='bkpts',form='formatted',status='old')
566 2 : REWIND (202)
567 2 : READ (202,'(i4)') nntot
568 2 : IF(l_p0)THEN
569 1 : WRITE (*,*) 'nntot=',nntot
570 1 : WRITE(*,*) 'fullnkpts=',fullnkpts
571 1 : WRITE(*,*) 'nkpts=',nkpts
572 : ENDIF
573 6 : ALLOCATE ( gb(1:3,1:nntot,1:fullnkpts),bpt(1:nntot,1:fullnkpts))
574 18 : DO ikpt=1,fullnkpts
575 146 : DO nn=1,nntot
576 : READ (202,'(2i6,3x,3i4)')&
577 128 : ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
578 128 : IF (ikpt/=ikpt_help) CALL juDFT_error("ikpt.ne.ikpt_help"&
579 0 : ,calledby ="wannier")
580 128 : IF (bpt(nn,ikpt)>fullnkpts) CALL juDFT_error("bpt.gt.fullnkpts"&
581 16 : ,calledby ="wannier")
582 : ENDDO
583 : ENDDO
584 2 : CLOSE (202)
585 2 : ALLOCATE(kdiff(3,nntot))
586 : ENDIF
587 :
588 : !**********************************************************
589 : !cccccccccccccc read in the bqpts file ccccccccccccccccc
590 : !**********************************************************
591 2 : IF ((wann%l_matrixmmn).AND.(l_gwf.OR.wann%l_ms)) THEN
592 0 : l_bqpts = .FALSE.
593 0 : INQUIRE (file='bqpts',exist=l_bqpts)
594 0 : IF (.NOT.l_bqpts) CALL juDFT_error("need bqpts for matrixmmn"&
595 0 : ,calledby ="wannier")
596 0 : OPEN (202,file='bqpts',form='formatted',status='old')
597 0 : REWIND (202)
598 0 : READ (202,'(i4)') nntot_q
599 0 : IF(l_p0)THEN
600 0 : WRITE (*,*) 'nntot_q=',nntot_q
601 0 : WRITE(*,*) 'fullnqpts=',fullnqpts
602 : ENDIF
603 : ALLOCATE ( gb_q(1:3,1:nntot_q,1:fullnqpts),&
604 0 : bpt_q(1:nntot_q,1:fullnqpts))
605 0 : DO iqpt=1,fullnqpts
606 0 : DO nn=1,nntot_q
607 : READ (202,'(2i6,3x,3i4)')&
608 0 : iqpt_help,bpt_q(nn,iqpt),(gb_q(i,nn,iqpt),i=1,3)
609 0 : IF (iqpt/=iqpt_help) CALL juDFT_error("iqpt.ne.iqpt_help"&
610 0 : ,calledby ="wannier")
611 0 : IF (bpt_q(nn,iqpt)>fullnqpts)&
612 0 : CALL juDFT_error("bpt_q.gt.fullnqpts",calledby ="wannier")
613 : ENDDO
614 : ENDDO
615 0 : CLOSE (202)
616 0 : ALLOCATE(qdiff(3,nntot_q))
617 0 : ALLOCATE(zero_qdiff(3,nntot_q))
618 2 : zero_qdiff=0.0
619 : ENDIF
620 :
621 :
622 : ! when treating gen. WF for spin spirals, the Brillouin zone
623 : ! of q-points is twice as large compared to k-BZ. Thus,
624 : ! the G-vectors connecting neighbors across the boundary
625 : ! need to be doubled
626 2 : IF(wann%l_sgwf) gb_q = 2*gb_q
627 2 : IF(wann%l_socgwf) gb_q = 2*gb_q
628 :
629 2 : IF(wann%l_finishgwf) GOTO 9110
630 : !********************************************************
631 : ! find symmetry-related elements in mmkb
632 : !********************************************************
633 2 : IF(wann%l_matrixmmn)THEN
634 6 : ALLOCATE(maptopair(3,fullnkpts,nntot))
635 2 : ALLOCATE(pair_to_do(fullnkpts,nntot))
636 : CALL wann_mmnk_symm(input,kpts,&
637 : fullnkpts,nntot,bpt,gb,wann%l_bzsym,&
638 : irreduc,mapkoper,l_p0,input%film,sym%nop,sym%invtab,sym%mrot,&
639 : sym%tau,&
640 2 : pair_to_do,maptopair,kdiff,.FALSE.,wann%param_file)
641 : ENDIF
642 :
643 : ! do the same for q-points to construct GWFs
644 2 : IF(wann%l_matrixmmn.AND.l_gwf)THEN
645 0 : ALLOCATE(maptopair_q(3,fullnqpts,nntot_q))
646 0 : ALLOCATE(pair_to_do_q(fullnqpts,nntot_q))
647 : CALL wann_mmnk_symm(input,kpts,&
648 : fullnqpts,nntot_q,bpt_q,gb_q,wann%l_bzsym,&
649 : irreduc_q,mapqoper,l_p0,.FALSE.,1,sym%invtab(1),&
650 : sym%mrot(:,:,1),sym%tau,&
651 0 : pair_to_do_q,maptopair_q,qdiff,.TRUE.,wann%param_file)
652 : ENDIF
653 :
654 :
655 : !*********************************************************
656 : !ccccccccccccccc initialize the potential cccccccccccc
657 : !*********************************************************
658 :
659 6 : ALLOCATE ( vr(atoms%jmtd,atoms%ntype,input%jspins) )
660 4 : ALLOCATE ( vso(atoms%jmtd,atoms%nat,2) )
661 2 : ALLOCATE ( vz(vacuum%nmzd,2,4) )
662 :
663 4026 : vz = 0.0
664 :
665 2 : IF(input%film)THEN
666 0 : vz(:,:,:SIZE(vTot%vac,4)) = REAL(vTot%vac(:,1,:,:))
667 0 : IF (SIZE(vTot%vac,4)==3) vz(:,:,4) = AIMAG(vTot%vac(:,1,:,3))
668 : END IF
669 :
670 4 : DO jspin = 1,input%jspins
671 6 : DO n = 1, atoms%ntype
672 946 : DO j = 1,atoms%jri(n)
673 944 : vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
674 : ENDDO
675 : ENDDO
676 : ENDDO
677 :
678 2 : IF(wann%l_soctomom)THEN
679 0 : CALL vsoc(input,atoms,vr,enpara%el0,.TRUE., vso)
680 : ENDIF
681 :
682 2 : IF(noco%l_noco.AND.input%film)THEN
683 0 : npotmatfile=25
684 0 : ALLOCATE(vpw(stars%ng3,1))
685 0 : ALLOCATE( vxy(vacuum%nmzxyd,stars%ng2-1,2) )
686 :
687 : OPEN (npotmatfile,FILE='potmat',FORM='unformatted',&
688 0 : STATUS='old')
689 0 : READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
690 0 : READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
691 0 : READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
692 0 : maxvac=2
693 0 : DO ivac = 1,maxvac
694 : !---> if the two vacuua are equivalent, the potential file has to
695 : !---> be backspaced, because the potential is the same at both
696 : !---> surfaces of the film
697 0 : IF ((ivac.EQ.2) .AND. (vacuum%nvac.EQ.1)) THEN
698 0 : DO irec = 1,4
699 0 : BACKSPACE (npotmatfile)
700 : ENDDO
701 : ENDIF
702 : !---> load the non-warping part of the potential
703 : READ (npotmatfile)&
704 0 : ((vz(imz,ivac,ipot),imz=1,vacuum%nmzd),ipot=1,4)
705 :
706 0 : DO ipot = 1,3
707 0 : READ (npotmatfile)((vxy(imz,igvm2,ivac),&
708 0 : imz=1,vacuum%nmzxy),igvm2=1,stars%ng2-1)
709 : ENDDO
710 : ENDDO
711 0 : CLOSE (npotmatfile)
712 : ENDIF
713 :
714 : !ccccccccccccccc end of the potential part ccccccccccc
715 2 : wannierspin=input%jspins
716 2 : IF(noco%l_soc) wannierspin=2
717 :
718 :
719 8 : ALLOCATE ( ff(atoms%ntype,atoms%jmtd,2,0:atoms%lmaxd,2) )
720 4 : ALLOCATE ( gg(atoms%ntype,atoms%jmtd,2,0:atoms%lmaxd,2) )
721 4 : ALLOCATE ( usdus%us(0:atoms%lmaxd,atoms%ntype,2) )
722 2 : ALLOCATE ( usdus%uds(0:atoms%lmaxd,atoms%ntype,2) )
723 2 : ALLOCATE ( usdus%dus(0:atoms%lmaxd,atoms%ntype,2) )
724 2 : ALLOCATE ( usdus%duds(0:atoms%lmaxd,atoms%ntype,2) )
725 2 : ALLOCATE ( usdus%ddn(0:atoms%lmaxd,atoms%ntype,2) )
726 4 : ALLOCATE ( usdus%ulos(atoms%nlod,atoms%ntype,2) )
727 2 : ALLOCATE ( usdus%dulos(atoms%nlod,atoms%ntype,2) )
728 2 : ALLOCATE ( usdus%uulon(atoms%nlod,atoms%ntype,2) )
729 2 : ALLOCATE ( usdus%dulon(atoms%nlod,atoms%ntype,2) )
730 8 : ALLOCATE ( usdus%uloulopn(atoms%nlod,atoms%nlod,atoms%ntype,2) )
731 :
732 2 : IF(l_gwf.AND..NOT.(wann%l_wann_plot)) THEN
733 : doublespin_max=4!2
734 : ELSE
735 2 : doublespin_max=wannierspin
736 : ENDIF
737 :
738 : !*****************************************************************c
739 : ! START Q LOOP c
740 : ! standard functionality of code for fullnqpts = nntot_q = 1 c
741 : ! and wann%l_ms = wann%l_sgwf = wann%l_socgwf = F c
742 : !*****************************************************************c
743 4 : DO iqpt = 1,fullnqpts ! loop by q-points starts
744 :
745 2 : ALLOCATE(innerEig_idList(nntot_q))
746 :
747 2 : qptibz=iqpt
748 2 : IF(wann%l_bzsym .AND. l_gwf) qptibz=irreduc_q(iqpt)
749 : IF(wann%l_bzsym .AND. l_gwf) oper_q=mapqoper(iqpt)
750 :
751 8 : qpt_i = nococonv%qss
752 4 : alph_i = nococonv%alph
753 4 : beta_i = nococonv%beta
754 2 : theta_i = nococonv%theta
755 2 : phi_i = nococonv%phi
756 2 : IF(wann%l_sgwf.OR.wann%l_ms) THEN
757 0 : qpt_i(:) = wann%param_vec(:,qptibz)
758 0 : alph_i(:) = wann%param_alpha(:,qptibz)
759 2 : ELSEIF(wann%l_socgwf) THEN
760 0 : IF(wann%l_dim(2)) phi_i = tpi_const*wann%param_vec(2,qptibz)
761 0 : IF(wann%l_dim(3)) theta_i = tpi_const*wann%param_vec(3,qptibz)
762 : ENDIF
763 :
764 2 : IF (l_gwf) THEN
765 0 : IF(wann%l_matrixmmn)THEN
766 0 : DO iqpt_b=1,nntot_q
767 :
768 0 : innerEig_idList(iqpt_b) = eig_idList(bpt_q(iqpt_b,iqpt))
769 :
770 : ! WRITE(fending,'("_",i4.4)')bpt_q(iqpt_b,iqpt)
771 : ! innerEig_idList(iqpt_b)=open_eig(fmpi%mpi_comm,
772 : ! + lapw%dim_nbasfcn(),input%neig,
773 : ! + nkpts,wannierspin,atoms%lmaxd,
774 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
775 : ! + noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
776 : ! + mpi%n_size,filename=trim(fstart)//fending,
777 : ! + layers=banddos%layers,nstars=banddos%nstars,
778 : ! + ncored=29,nsld=atoms%nat,
779 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
780 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
781 :
782 : ENDDO
783 : ENDIF
784 :
785 0 : eig_id = eig_idList(qptibz)
786 :
787 : ! WRITE(fending,'("_",i4.4)')qptibz
788 : ! eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
789 : ! + nkpts,wannierspin,atoms%lmaxd,
790 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
791 : ! + noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
792 : ! + mpi%n_size,filename=trim(fstart)//fending,
793 : ! + layers=banddos%layers,nstars=banddos%nstars,
794 : ! + ncored=29,nsld=atoms%nat,
795 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
796 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
797 :
798 2 : ELSEIF(wann%l_ms) THEN
799 :
800 0 : eig_id = eig_idList(qptibz)
801 :
802 : ! WRITE(fending,'("_",i4.4)')qptibz
803 : ! eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
804 : ! + nkpts,wannierspin,atoms%lmaxd,
805 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
806 : ! + noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
807 : ! + mpi%n_size,filename=trim(fstart)//fending,
808 : ! + layers=banddos%layers,nstars=banddos%nstars,
809 : ! + ncored=29,nsld=atoms%nat,
810 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
811 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
812 :
813 : ENDIF ! l_gwf.or.wann%l_ms
814 2 : nrec=0
815 2 : nrec_b=0
816 :
817 :
818 : !****************************************************
819 : ! cycle by spins starts!
820 : !****************************************************
821 4 : DO doublespin=1,doublespin_max ! cycle by spins
822 :
823 2 : jspin=MOD(doublespin+1,2)+1
824 2 : jspin_b=jspin
825 2 : IF(doublespin.EQ.3) jspin_b=2
826 2 : IF(doublespin.EQ.4) jspin_b=1
827 :
828 2 : nrec_b = nrec
829 :
830 2 : IF(.NOT.noco%l_noco) THEN
831 2 : nrec = (jspin-1)*nkpts
832 2 : nrec_b = (jspin_b-1)*nkpts
833 : ENDIF
834 :
835 : ! spin-dependent sign of the q-dependent phase
836 : ! in the generalized Bloch theorem
837 : ! -1: spin up, +1: spin down
838 2 : sign_q = -sign_q
839 :
840 : !...read number of bands and wannier functions from file proj
841 :
842 : !..reading the proj.1 / proj.2 / proj file
843 2 : l_proj=.FALSE.
844 4 : DO j=jspin,0,-1
845 4 : INQUIRE(file=TRIM('proj'//spin012(j)),exist=l_proj)
846 4 : IF(l_proj)THEN
847 2 : filename='proj'//spin012(j)
848 2 : EXIT
849 : ENDIF
850 : ENDDO
851 :
852 2 : IF(l_proj)THEN
853 2 : OPEN (203,file=TRIM(filename),status='old')
854 2 : REWIND (203)
855 2 : READ (203,*) nwfs,numbands
856 2 : REWIND (203)
857 2 : CLOSE (203)
858 : ELSEIF(wann%l_projmethod.OR.wann%l_bestproj&
859 0 : .OR.wann%l_matrixamn)THEN
860 : CALL juDFT_error("no proj/proj.1/proj.2"&
861 0 : ,calledby ="wannier")
862 : ENDIF
863 :
864 :
865 2 : jspin2=jspin
866 2 : IF(noco%l_soc .AND. input%jspins.EQ.1)jspin2=1
867 2 : jspin2_b=jspin_b
868 2 : IF(noco%l_soc .AND. input%jspins.EQ.1)jspin2_b=1
869 :
870 2 : jsp_start = jspin ; jsp_end = jspin
871 :
872 : !ccccccccccc read in the eigenvalues and vectors cccccc
873 2 : WRITE(*,*)'wannierspin',wannierspin
874 4 : DO jspin5=1,wannierspin!1!2
875 : ! jspin5=jspin
876 2 : jsp_start=jspin5; jsp_end=jspin5
877 2 : nrec5=0
878 : IF(.NOT.noco%l_noco) nrec5 = (jspin5-1)*nkpts
879 :
880 : CALL cdn_read0(eig_id,fmpi%irank,fmpi%isize,jspin5,input%jspins, &!wannierspin instead of DIMENSION%jspd?&
881 4 : noco%l_noco, n_bands,n_size)
882 :
883 : ENDDO
884 : !.. now we want to define the maximum number of the bands by all kpts
885 2 : nbnd = 0
886 :
887 2 : i_rec = 0 ; n_rank = 0
888 : !*************************************************************
889 : !..writing down the eig.1 and/or eig.2 files
890 :
891 : !..write individual files if multi-spiral mode wann%l_ms=T
892 : !*************************************************************
893 2 : IF(l_p0)THEN
894 : CALL wann_write_eig(&
895 : fmpi,cell,noco,nococonv,input,kpts,sym,atoms, &
896 : eig_id,l_real,&
897 : atoms%ntype,&
898 : lapw%dim_nvd(),wannierspin,&
899 : fmpi%isize,jspin,&
900 : noco%l_ss,noco%l_noco,nrec,fullnkpts,&
901 : wann%l_bzsym,wann%l_byindex,wann%l_bynumber,&
902 : wann%l_byenergy,&
903 : irreduc ,wann%band_min(jspin),&
904 : wann%band_max(jspin),&
905 : numbands,&
906 : sliceplot%e1s,sliceplot%e2s,efermi,.FALSE.,nkpts,&
907 1 : nbnd,kpoints,l_gwf,iqpt)
908 :
909 : ENDIF!l_p0
910 :
911 : ! nbnd is calculated for process zero and is sent here to the others
912 : #ifdef CPP_MPI
913 2 : IF(l_p0)THEN
914 2 : DO cpu_index=1,fmpi%isize-1
915 2 : CALL MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,fmpi%mpi_comm,ierr(1))
916 : ENDDO
917 : ELSE
918 1 : CALL MPI_RECV(nbnd,1,MPI_INTEGER,0,1,fmpi%mpi_comm,stt,ierr(1))
919 : ENDIF
920 : #endif
921 :
922 2 : PRINT*,"process: ",fmpi%irank," nbnd= ",nbnd
923 : !##################################################################
924 2 : IF(wann%l_mmn0)THEN
925 6 : ALLOCATE ( mmn(nbnd,nbnd,fullnkpts) )
926 1170 : mmn(:,:,:) = CMPLX(0.,0.)
927 2 : IF((noco%l_soc.OR.noco%l_noco) .AND. (doublespin.EQ.1))&
928 0 : ALLOCATE(socmmn(nbnd,nbnd,fullnkpts) )
929 : ENDIF
930 2 : IF(wann%l_nabla)THEN
931 0 : ALLOCATE ( nablamat(3,nbnd,nbnd,fullnkpts) )
932 0 : nablamat = CMPLX(0.,0.)
933 : ENDIF
934 :
935 2 : IF(wann%l_soctomom)THEN
936 0 : ALLOCATE ( soctomom(3,nbnd,nbnd,fullnkpts) )
937 0 : soctomom = CMPLX(0.,0.)
938 : ENDIF
939 :
940 2 : IF(wann%l_surfcurr)THEN
941 0 : ALLOCATE ( surfcurr(3,nbnd,nbnd,fullnkpts) )
942 0 : surfcurr = CMPLX(0.,0.)
943 : ENDIF
944 :
945 2 : IF(wann%l_anglmom)THEN
946 0 : IF(.NOT.ALLOCATED(anglmom))THEN
947 0 : ALLOCATE ( anglmom(3,nbnd,nbnd,fullnkpts) )
948 0 : anglmom=CMPLX(0.,0.)
949 : ENDIF
950 : ENDIF
951 :
952 2 : IF(wann%l_orbcomp)THEN
953 0 : IF(ALLOCATED(orbcomp))DEALLOCATE(orbcomp)
954 0 : IF(wann%l_oc_f)THEN
955 0 : ALLOCATE(orbcomp(16,wann%oc_num_orbs,nbnd,nbnd,fullnkpts))
956 : ELSE
957 0 : ALLOCATE(orbcomp(9,wann%oc_num_orbs,nbnd,nbnd,fullnkpts))
958 : ENDIF
959 0 : orbcomp=CMPLX(0.,0.)
960 : ENDIF
961 :
962 : !write (*,*) 'nwfs=',nwfs
963 2 : IF(wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn)THEN
964 2 : IF(.NOT.ALLOCATED(amn))THEN
965 6 : ALLOCATE ( amn(nbnd,nwfs,fullnkpts) )
966 1170 : amn(:,:,:) = CMPLX(0.,0.)
967 : ENDIF
968 : ENDIF
969 :
970 2 : IF (wann%l_projmethod.OR.wann%l_bestproj) THEN
971 0 : ALLOCATE ( psiw(nbnd,nwfs,fullnkpts) )
972 0 : psiw(:,:,:) = CMPLX(0.,0.)
973 0 : IF(.NOT.ALLOCATED(hwfr))THEN
974 0 : ALLOCATE ( hwfr(nwfs,nwfs) )
975 0 : hwfr(:,:) = CMPLX(0.,0.)
976 : ENDIF
977 : ENDIF
978 :
979 :
980 2 : IF (wann%l_matrixmmn) THEN
981 2 : IF(.NOT.ALLOCATED(mmnk))THEN
982 8 : ALLOCATE ( mmnk(nbnd,nbnd,nntot,fullnkpts) )
983 9362 : mmnk = (0.,0.)
984 : ENDIF
985 : ENDIF
986 :
987 2 : IF(wann%l_matrixmmn)THEN
988 2 : IF(.NOT.ALLOCATED(mmnk_q).AND.l_gwf)THEN
989 0 : ALLOCATE ( mmnk_q (nbnd,nbnd,nntot_q,fullnkpts) )
990 0 : mmnk_q = (0.,0.)
991 :
992 : ! allocate ( m_int(nbnd,nbnd,nntot_q,fullnkpts) )
993 : ! allocate ( m_sph(nbnd,nbnd,nntot_q,fullnkpts) )
994 : ! allocate ( m_vac(nbnd,nbnd,nntot_q,fullnkpts) )
995 : ! m_int = cmplx(0.,0.)
996 : ! m_sph = cmplx(0.,0.)
997 : ! m_vac = cmplx(0.,0.)
998 : ENDIF
999 : ENDIF
1000 :
1001 :
1002 8 : ALLOCATE ( flo(atoms%ntype,atoms%jmtd,2,atoms%nlod,2) )
1003 :
1004 4 : DO jspin4=1,wannierspin!2
1005 2 : jspin3=jspin4
1006 2 : IF(input%jspins.EQ.1) jspin3=1
1007 2 : na = 1
1008 6 : DO n = 1,atoms%ntype
1009 16 : DO l = 0,atoms%lmax(n)
1010 : !...compute the l-dependent, k-independent radial MT- basis functions
1011 :
1012 : CALL radfun(&
1013 : l,n,jspin4,enpara%el0(l,n,jspin3),vr(1,n,jspin3),atoms,&
1014 : ff(n,:,:,l,jspin4),gg(n,:,:,l,jspin4),usdus,&
1015 16 : nodeu,noded,wronk)
1016 :
1017 : ENDDO
1018 : !...and the local orbital radial functions
1019 4 : DO ilo = 1, atoms%nlo(n)
1020 :
1021 : CALL radflo(&
1022 : atoms,n,jspin4,enpara%ello0(:,:,jspin3),vr(1,n,jspin3),&
1023 : ff(n,1:,1:,0:,jspin4),gg(n,1:,1:,0:,jspin4),fmpi,&
1024 2 : usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin4))
1025 :
1026 : ENDDO
1027 : ! na = na + atoms%neq(n)
1028 : ENDDO
1029 : ENDDO!jspin3
1030 : !****************************************************************
1031 : ! calculate the k-independent uju*gaunt-matrix needed for
1032 : ! mmnmatrix
1033 : !****************************************************************
1034 : ! TODO: make this more efficient (i.e., compute ujugaunt only once
1035 : ! and not for all q-points).
1036 2 : IF(wann%l_matrixmmn)THEN
1037 0 : ALLOCATE(ujug(0:lmd,0:lmd,&
1038 8 : 1:atoms%ntype,1:nntot))
1039 0 : ALLOCATE(ujdg(0:lmd,0:lmd,&
1040 8 : 1:atoms%ntype,1:nntot))
1041 0 : ALLOCATE(djug(0:lmd,0:lmd,&
1042 6 : 1:atoms%ntype,1:nntot))
1043 0 : ALLOCATE(djdg(0:lmd,0:lmd,&
1044 6 : 1:atoms%ntype,1:nntot))
1045 0 : ALLOCATE(ujulog(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1046 10 : 1:atoms%ntype,1:nntot))
1047 0 : ALLOCATE(djulog(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1048 8 : 1:atoms%ntype,1:nntot))
1049 0 : ALLOCATE(ulojug(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1050 8 : 1:atoms%ntype,1:nntot))
1051 0 : ALLOCATE(ulojdg(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1052 8 : 1:atoms%ntype,1:nntot))
1053 0 : ALLOCATE(ulojulog(1:atoms%nlod,-atoms%llod:atoms%llod,&
1054 : 1:atoms%nlod,-atoms%llod:atoms%llod,&
1055 12 : 1:atoms%ntype,1:nntot))
1056 :
1057 : CALL wann_ujugaunt(&
1058 : atoms%llod,nntot,kdiff,atoms%lmax,atoms%ntype,&
1059 : atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
1060 : atoms%llo,flo(:,:,:,:,jspin),&
1061 : flo(:,:,:,:,jspin),&
1062 : ff(:,:,:,:,jspin),&
1063 : ff(:,:,:,:,jspin),&
1064 : gg(:,:,:,:,jspin),&
1065 : gg(:,:,:,:,jspin),atoms%jri,atoms%rmsh,atoms%dx,&
1066 : atoms%jmtd,atoms%lmaxd,lmd,&
1067 : ujug,ujdg,djug,djdg,&
1068 2 : ujulog,djulog,ulojug,ulojdg,ulojulog,.FALSE.,1)
1069 :
1070 : ! compute integrals of radial solution, according energy derivatives,
1071 : ! the spherical Bessel function and the Gaunt coefficients in order
1072 : ! to account for the overlap of the lattice periodic parts at
1073 : ! neighboring q-points
1074 2 : IF(l_gwf)THEN
1075 0 : ALLOCATE(ujug_q(0:lmd,0:lmd,&
1076 0 : 1:atoms%ntype,1:nntot_q))
1077 0 : ALLOCATE(ujdg_q(0:lmd,0:lmd,&
1078 0 : 1:atoms%ntype,1:nntot_q))
1079 0 : ALLOCATE(djug_q(0:lmd,0:lmd,&
1080 0 : 1:atoms%ntype,1:nntot_q))
1081 0 : ALLOCATE(djdg_q(0:lmd,0:lmd,&
1082 0 : 1:atoms%ntype,1:nntot_q))
1083 0 : ALLOCATE(ujulog_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1084 0 : 1:atoms%ntype,1:nntot_q))
1085 0 : ALLOCATE(djulog_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1086 0 : 1:atoms%ntype,1:nntot_q))
1087 0 : ALLOCATE(ulojug_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1088 0 : 1:atoms%ntype,1:nntot_q))
1089 0 : ALLOCATE(ulojdg_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
1090 0 : 1:atoms%ntype,1:nntot_q))
1091 0 : ALLOCATE(ulojulog_q(1:atoms%nlod,-atoms%llod:atoms%llod,&
1092 : 1:atoms%nlod,-atoms%llod:atoms%llod,&
1093 0 : 1:atoms%ntype,1:nntot_q))
1094 :
1095 : ! we need G(q+b)/2 as argument for the sph. Bessel func.
1096 : ! and additionally a spin-dependent sign (-/+ 1)^{lpp}
1097 0 : IF(wann%l_sgwf) CALL wann_ujugaunt(&
1098 : atoms%llod,nntot_q,qdiff/2.0,atoms%lmax,atoms%ntype,&
1099 : atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
1100 : atoms%llo,flo(:,:,:,:,jspin),&
1101 : flo(:,:,:,:,jspin_b),&
1102 : ff(:,:,:,:,jspin),&
1103 : ff(:,:,:,:,jspin_b),&
1104 : gg(:,:,:,:,jspin),&
1105 : gg(:,:,:,:,jspin_b),atoms%jri,atoms%rmsh,atoms%dx,&
1106 : atoms%jmtd,atoms%lmaxd,lmd,&
1107 : ujug_q,ujdg_q,djug_q,djdg_q,&
1108 : ujulog_q,djulog_q,ulojug_q,ulojdg_q,ulojulog_q,.TRUE.,&
1109 0 : sign_q)
1110 :
1111 0 : IF(wann%l_socgwf) CALL wann_ujugaunt(&
1112 : atoms%llod,nntot_q,zero_qdiff,atoms%lmax,atoms%ntype,&
1113 : atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
1114 : atoms%llo,flo(:,:,:,:,jspin),&
1115 : flo(:,:,:,:,jspin_b),&
1116 : ff(:,:,:,:,jspin),&
1117 : ff(:,:,:,:,jspin_b),&
1118 : gg(:,:,:,:,jspin),&
1119 : gg(:,:,:,:,jspin_b),atoms%jri,atoms%rmsh,atoms%dx,&
1120 : atoms%jmtd,&
1121 : atoms%lmaxd,lmd,ujug_q,ujdg_q,djug_q,djdg_q,&
1122 : ujulog_q,djulog_q,ulojug_q,ulojdg_q,ulojulog_q,&
1123 0 : .FALSE.,1)
1124 :
1125 : ENDIF ! l_gwf
1126 :
1127 : ENDIF !l_matrixmmn
1128 : ! Replace the following by calls to zmat%init further below
1129 : ! zzMat%l_real = l_real
1130 : ! zzMat%matsize1 = lapw%dim_nbasfcn()
1131 : ! zzMat%matsize2 = input%neig
1132 : ! IF(l_real) THEN
1133 : ! IF(.NOT.ALLOCATED(zzMat%data_r))&
1134 : ! ALLOCATE (zzMat%data_r(zzMat%matsize1,zzMat%matsize2))
1135 : ! ELSE
1136 : ! IF(.NOT.ALLOCATED(zzMat%data_c))&
1137 : ! ALLOCATE (zzMat%data_c(zzMat%matsize1,zzMat%matsize2))
1138 : ! END IF
1139 :
1140 : ! zMat%l_real = zzMat%l_real
1141 : ! zMat%matsize1 = zzMat%matsize1
1142 : ! zMat%matsize2 = zzMat%matsize2
1143 : ! IF (zzMat%l_real) THEN
1144 : ! IF(.NOT.ALLOCATED(zMat%data_r))&
1145 : ! ALLOCATE (zMat%data_r(zMat%matsize1,zMat%matsize2))
1146 : ! zMat%data_r = 0.0
1147 : ! ELSE
1148 : ! IF(.NOT.ALLOCATED(zMat%data_c))&
1149 : ! ALLOCATE (zMat%data_c(zMat%matsize1,zMat%matsize2))
1150 : ! zMat%data_c = CMPLX(0.0,0.0)
1151 : ! END IF
1152 :
1153 : ! zMat_b%l_real = zzMat%l_real
1154 : ! zMat_b%matsize1 = zzMat%matsize1
1155 : ! zMat_b%matsize2 = zzMat%matsize2
1156 : ! IF (zzMat%l_real) THEN
1157 : ! IF(.NOT.ALLOCATED(zMat_b%data_r))&
1158 : ! ALLOCATE (zMat_b%data_r(zMat_b%matsize1,zMat_b%matsize2))
1159 : ! zMat_b%data_r = 0.0
1160 : ! ELSE
1161 : ! IF(.NOT.ALLOCATED(zMat_b%data_c))&
1162 : ! ALLOCATE (zMat_b%data_c(zMat_b%matsize1,zMat_b%matsize2))
1163 : ! zMat_b%data_c = CMPLX(0.0,0.0)
1164 : ! END IF
1165 :
1166 2 : i_rec = 0 ; n_rank = 0
1167 :
1168 : !****************************************************************
1169 : !.. loop by kpoints starts! each may be a separate task
1170 : !****************************************************************
1171 18 : DO ikpt = wann%ikptstart,fullnkpts ! loop by k-points starts
1172 16 : kptibz=ikpt
1173 16 : IF(wann%l_bzsym) kptibz=irreduc(ikpt)
1174 16 : IF(wann%l_bzsym) oper=mapkoper(ikpt)
1175 :
1176 16 : i_rec = i_rec + 1
1177 18 : IF (MOD(i_rec-1,fmpi%isize).EQ.fmpi%irank) THEN
1178 :
1179 8 : ALLOCATE ( we(input%neig),eigg(input%neig) )
1180 :
1181 8 : n_start=1
1182 8 : n_end=input%neig
1183 :
1184 :
1185 : ! read information of diagonalization for fixed q-point iqpt
1186 : ! stored in the eig file on unit 66. the lattice respectively
1187 : ! plane-wave vectors G(k,q) are saved in (k1,k2,k3).
1188 :
1189 8 : CALL lapw%init(input,noco,nococonv,kpts,atoms,sym,kptibz,cell,fmpi)
1190 :
1191 8 : nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
1192 8 : CALL zzMat%init(l_real,nbasfcn,input%neig)
1193 8 : CALL zMat%init(l_real,nbasfcn,input%neig)
1194 :
1195 :
1196 : CALL cdn_read(&
1197 : eig_id,&
1198 : lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd?&
1199 : kptibz,jspin,lapw%dim_nbasfcn(),&
1200 : noco%l_ss,noco%l_noco,input%neig,n_start,n_end,&
1201 8 : nbands,eigg,zzMat)
1202 :
1203 :
1204 8 : nslibd = 0
1205 :
1206 : !...we work only within the energy window
1207 :
1208 72 : eig(:) = 0.
1209 :
1210 : ! print*,"bands used:"
1211 :
1212 72 : DO i = 1,nbands
1213 : IF ((eigg(i).GE.sliceplot%e1s.AND.nslibd.LT.numbands.AND.&
1214 : wann%l_bynumber).OR.&
1215 : (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
1216 64 : wann%l_byenergy).OR.(i.GE.wann%band_min(jspin2).AND.&
1217 8 : (i.LE.wann%band_max(jspin2)).AND.wann%l_byindex))THEN
1218 :
1219 : ! print*,i
1220 64 : nslibd = nslibd + 1
1221 64 : eig(nslibd) = eigg(i)
1222 64 : we(nslibd) = we(i)
1223 64 : IF(zzMat%l_real) THEN
1224 6056 : zMat%data_r(:,nslibd) = zzMat%data_r(:,i)
1225 : ELSE
1226 0 : zMat%data_c(:,nslibd) = zzMat%data_c(:,i)
1227 : END IF
1228 : ENDIF
1229 : ENDDO
1230 :
1231 : !***********************************************************
1232 : ! rotate the wavefunction
1233 : !***********************************************************
1234 8 : IF (wann%l_bzsym.AND.oper.NE.1) THEN !rotate bkpt
1235 : call wann_kptsrotate(&
1236 : atoms, &
1237 : sym%invsat,&
1238 : noco%l_noco,noco%l_soc,&
1239 : jspin,&
1240 : oper,sym%nop,sym%mrot,lapw%dim_nvd(),&
1241 : shiftkpt(:,ikpt),&
1242 : sym%tau,&
1243 : lapw,&
1244 0 : zMat,nsfactor)
1245 : ELSE
1246 8 : nsfactor=CMPLX(1.0,0.0)
1247 : ENDIF
1248 :
1249 : !******************************************************************
1250 :
1251 : !...the overlap matrix Mmn which is computed for each k- and b-point
1252 :
1253 8 : noccbd = nslibd
1254 :
1255 0 : ALLOCATE(acof(noccbd,0:lmd,atoms%nat),&
1256 0 : bcof(noccbd,0:lmd,atoms%nat),&
1257 56 : ccof(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat))
1258 :
1259 14152 : acof(:,:,:) = CMPLX(0.,0.) ; bcof(:,:,:) = CMPLX(0.,0.)
1260 24 : ccof(:,:,:,:) = CMPLX(0.,0.)
1261 :
1262 : !...generation the A,B,C coefficients in the spheres
1263 : !...for the lapws and local orbitals, summed by the basis functions
1264 :
1265 :
1266 : CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,&
1267 8 : noco,nococonv,jspin2 ,acof,bcof,ccof,zMat)
1268 :
1269 :
1270 8 : CALL wann_abinv(atoms,sym, acof,bcof,ccof)
1271 :
1272 :
1273 8 : IF((doublespin.EQ.3).OR.(doublespin.EQ.4)) GOTO 9900
1274 :
1275 :
1276 8 : IF(wann%l_orbcomp)THEN
1277 : CALL wann_orbcomp(atoms,usdus,jspin,acof,bcof,ccof,&
1278 0 : wann%oc_num_orbs,wann%oc_orbs,wann%l_oc_f,orbcomp(:,:,:,:,ikpt))
1279 : ENDIF
1280 :
1281 8 : IF(wann%l_anglmom)THEN
1282 0 : CALL wann_anglmom(atoms,usdus,jspin, acof,bcof,ccof,anglmom(:,:,:,ikpt))
1283 : ENDIF
1284 :
1285 : #ifdef CPP_TOPO
1286 : IF(wann%l_surfcurr)THEN
1287 : ! call wann_surfcurr_int(
1288 : ! > lapw%dim_nv2d(),jspin,stars%ng3,vacuum%nmzxyd,stars%ng2,sphhar%ntypsd,
1289 : ! > atoms%ntype,atoms%lmaxd,atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,stars%ng3,vacuum%nvac,
1290 : ! > vacuum%nmz,vacuum%nmzxy,stars%ng2,sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,sym%symor,
1291 : ! > sym%invs,sym%invs2,cell%z1,vacuum%delz,sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,
1292 : ! > atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,sliceplot%kk,
1293 : ! > lapw%dim_nvd(),atoms%nlod,atoms%llod,nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,
1294 : ! > k1(:,jspin),k2(:,jspin),k3(:,jspin),evac(:,jspin),
1295 : ! > vz(:,:,jspin2),
1296 : ! > nslibd,lapw%dim_nbasfcn(),input%neig,ff,gg,flo,acof,bcof,ccof,z,
1297 : ! > surfcurr(:,:,:,ikpt))
1298 :
1299 : CALL wann_surfcurr_int2(&
1300 : lapw%dim_nv2d(),jspin ,stars%ng3,&
1301 : vacuum%nmzxyd,&
1302 : stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
1303 : atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,&
1304 : atoms%neq,stars%ng3,vacuum%nvac,vacuum%nmz,&
1305 : vacuum%nmzxy,stars%ng2,sym%nop,sym%nop2,cell%volint,&
1306 : input%film,sliceplot%slice,sym%symor,&
1307 : sym%invs,sym%invs2,cell%z1,vacuum%delz,sym%ngopr,&
1308 : sym%ntypsy,atoms%jri,atoms%pos,atoms%taual,&
1309 : atoms%zatom,atoms%rmt,atoms%lmax,sym%mrot,sym%tau,&
1310 : atoms%rmsh,sym%invtab,cell%amat,cell%bmat,cell%bbmat,&
1311 : ikpt,lapw%dim_nvd(),lapw%nv(jspin),lapw%bkpt,cell%omtil,&
1312 : lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
1313 : nslibd,lapw%dim_nbasfcn(),input%neig,z,&
1314 : dirfacs,&
1315 : surfcurr(:,:,:,ikpt))
1316 :
1317 : CALL wann_surfcurr(&
1318 : dirfacs,cell%amat,&
1319 : jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
1320 : atoms%neq,noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod,&
1321 : atoms%nlo,atoms%llo, &
1322 : acof,bcof,ccof,&
1323 : us(:,:,jspin),dus(:,:,jspin),duds(:,:,jspin),&
1324 : uds(:,:,jspin),&
1325 : ulos(:,:,jspin),dulos(:,:,jspin),&
1326 : atoms%rmt,atoms%pos, &
1327 : surfcurr(:,:,:,ikpt))
1328 : WRITE(oUnit,*)"dirfacs=",dirfacs
1329 : ENDIF
1330 :
1331 : IF(wann%l_soctomom)THEN
1332 : CALL wann_soc_to_mom(&
1333 : jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
1334 : atoms%jmtd,atoms%jri,atoms%rmsh,atoms%dx,atoms%neq,&
1335 : noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod,&
1336 : vso(:,:,1), &
1337 : ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),&
1338 : acof,bcof,ccof,&
1339 : soctomom(:,:,:,ikpt))
1340 : ENDIF
1341 :
1342 : IF(wann%l_nabla)THEN
1343 : CALL wann_nabla(&
1344 : atoms%nlo,atoms%llo,&
1345 : jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
1346 : atoms%jmtd,atoms%jri,atoms%rmsh,atoms%dx,atoms%neq,&
1347 : noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod, &
1348 : ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),&
1349 : acof,bcof,ccof,&
1350 : nablamat(:,:,:,ikpt))
1351 : IF(input%film)THEN
1352 : CALL wann_nabla_vac(&
1353 : cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
1354 : stars%mx1,stars%mx2,stars%mx3,&
1355 : stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,vacuum%delz,&
1356 : stars%ig2,cell%area,cell%bmat,cell%bbmat,enpara%evac0(:,jspin),&
1357 : lapw%bkpt,vz(:,:,jspin2),nslibd,jspin,lapw%k1,lapw%k2,lapw%k3,&
1358 : wannierspin,lapw%dim_nvd(),lapw%dim_nbasfcn(),input%neig,z,nv,&
1359 : cell%omtil,&
1360 : nablamat(:,:,:,ikpt))
1361 : ENDIF
1362 : addnoco=0
1363 : DO i = n_rank+1,lapw%nv(jspin),n_size
1364 : b1 = lapw%bkpt+lapw%gvec(:,i,jspin)
1365 : b2(1)=b1(1)*cell%bmat(1,1)+b1(2)*cell%bmat(2,1)+&
1366 : b1(3)*cell%bmat(3,1)
1367 : b2(2)=b1(1)*cell%bmat(1,2)+b1(2)*cell%bmat(2,2)+&
1368 : b1(3)*cell%bmat(3,2)
1369 : b2(3)=b1(1)*cell%bmat(1,3)+b1(2)*cell%bmat(2,3)+&
1370 : b1(3)*cell%bmat(3,3)
1371 : DO j = n_rank+1,lapw%nv(jspin),n_size
1372 : !--> determine index and phase factor
1373 : i1 = lapw%k1(j,jspin) - lapw%k1(i,jspin)
1374 : i2 = lapw%k2(j,jspin) - lapw%k2(i,jspin)
1375 : i3 = lapw%k3(j,jspin) - lapw%k3(i,jspin)
1376 : in = stars%ig(i1,i2,i3)
1377 : IF (in.EQ.0) CYCLE
1378 : phase = stars%rgphs(i1,i2,i3)
1379 : phasust = CMPLX(phase,0.0)*stars%ustep(in)
1380 :
1381 : DO m = 1,nslibd
1382 : DO n = 1,nslibd
1383 : DO dir=1,3
1384 :
1385 : #if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
1386 : VALUE=phasust*z(i+addnoco,m)*CONJG(z(j+addnoco,n))
1387 : nablamat(dir,m,n,ikpt) =&
1388 : nablamat(dir,m,n,ikpt) -&
1389 : VALUE*b2(dir)
1390 : #else
1391 : VALUE=phasust*CMPLX(z(i+addnoco,m)*z(j+addnoco,n),0.0)
1392 : nablamat(dir,m,n,ikpt) =&
1393 : nablamat(dir,m,n,ikpt) - &
1394 : VALUE*b2(dir)
1395 : #endif
1396 : ENDDO
1397 : ENDDO
1398 : ENDDO
1399 :
1400 : ENDDO
1401 : ENDDO
1402 :
1403 : ENDIF
1404 : #endif
1405 : ! goto jump no longer needed?
1406 : ! if ((.not.wann%l_matrixmmn).and.(.not.wann%l_matrixamn).and.
1407 : ! (.not.wann%l_bestproj).and.(.not.wann%l_projmethod).and.
1408 : ! (.not.wann%l_mmn0)) goto 3
1409 :
1410 :
1411 : !------mmn0-matrix
1412 8 : IF(wann%l_mmn0)THEN
1413 8 : addnoco=0
1414 8 : noconbasfcn=nbasfcn
1415 8 : IF(noco%l_noco.AND.(jspin.EQ.2))THEN
1416 0 : addnoco=lapw%nv(1)+atoms%nlotot
1417 : ! noconbasfcn=nbasfcn-addnoco
1418 : ENDIF
1419 :
1420 : !-----> interstitial contribution to mmn0-matrix
1421 :
1422 : CALL wann_mmkb_int(&
1423 : cmplx_1,addnoco,addnoco,&
1424 : lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
1425 : stars%ng3,lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
1426 : lapw%nv(jspin2),input%neig,noconbasfcn,noconbasfcn,zMat,nslibd,&
1427 : lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
1428 : lapw%nv(jspin2),zMat,nslibd,&
1429 : nbnd,&
1430 : stars%rgphs,stars%ustep,stars%ig,(/ 0,0,0 /),&
1431 8 : mmn(:,:,ikpt))
1432 :
1433 : !---> spherical contribution to mmn0-matrix
1434 :
1435 : CALL wann_mmk0_sph(&
1436 : atoms%llod,noccbd,atoms%nlod,atoms%nat,atoms%ntype,&
1437 : atoms%lmaxd,atoms%lmax,lmd,atoms%ntype,atoms%neq,&
1438 : atoms%nlo,atoms%llo,acof(1:noccbd,:,:),&
1439 : bcof(1:noccbd,:,:),ccof(:,1:noccbd,:,:),&
1440 : usdus%ddn(:,:,jspin),usdus%uulon(:,:,jspin),&
1441 : usdus%dulon(:,:,jspin),usdus%uloulopn(:,:,:,jspin),&
1442 8 : mmn(:,:,ikpt))
1443 : !---> vacuum contribution to mmn0-matrix
1444 :
1445 8 : IF (input%film ) THEN
1446 :
1447 : CALL wann_mmk0_vac(&
1448 : noco%l_noco,atoms%nlotot,qpt_i,&
1449 : cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
1450 : stars%mx1,stars%mx2,stars%mx3,&
1451 : stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,vacuum%delz,&
1452 : stars%ig2,cell%area,cell%bmat,&
1453 : cell%bbmat,enpara%evac0(:,jspin2),lapw%bkpt,vz(:,:,jspin2),&
1454 : nslibd,jspin2,lapw%k1,lapw%k2,lapw%k3,wannierspin,lapw%dim_nvd(),&
1455 : lapw%dim_nbasfcn(),input%neig,zMat,lapw%nv,cell%omtil,&
1456 0 : mmn(:,:,ikpt))
1457 : ENDIF
1458 : ENDIF !l_mmn0
1459 :
1460 :
1461 : !---> overlaps with the trial orbitals
1462 8 : IF (wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn) THEN
1463 8 : l_amn2=.FALSE.
1464 8 : amnchi = cmplx_1!cmplx(1.,0.)
1465 :
1466 : CALL wann_amn (&
1467 : amnchi,nslibd,nwfs,atoms%ntype,atoms%nlod,atoms%llod,&
1468 : atoms%llo,atoms%nlo,atoms%lmaxd,atoms%jmtd,lmd,&
1469 : atoms%neq,atoms%nat,ikpt,nbnd,&
1470 : atoms%rmsh,atoms%rmt,atoms%jri,atoms%dx,atoms%lmax,&
1471 : usdus%us(:,:,jspin),usdus%dus(:,:,jspin),&
1472 : usdus%uds(:,:,jspin),&
1473 : usdus%duds(:,:,jspin),flo(:,:,:,:,jspin),&
1474 : ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),acof,bcof,ccof,&
1475 : (noco%l_soc.OR.noco%l_noco),jspin,&
1476 8 : l_amn2,amn(:,:,ikpt))
1477 8 : IF(l_amn2)THEN
1478 : CALL wann_amn (&
1479 : amnchi,nslibd,nwfs,atoms%ntype,atoms%nlod,atoms%llod,&
1480 : atoms%llo,atoms%nlo,atoms%lmaxd,atoms%jmtd,lmd,&
1481 : atoms%neq,atoms%nat,ikpt,nbnd,&
1482 : atoms%rmsh,atoms%rmt,atoms%jri,atoms%dx,atoms%lmax,&
1483 : usdus%us(:,:,jspin),usdus%dus(:,:,jspin),&
1484 : usdus%uds(:,:,jspin),&
1485 : usdus%duds(:,:,jspin),flo(:,:,:,:,jspin),&
1486 : ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),acof,bcof,ccof,&
1487 : (noco%l_soc.OR.noco%l_noco),jspin,&
1488 0 : l_amn2,amn(:,:,ikpt),lapw%bkpt)
1489 : ENDIF
1490 : ! amn(ikpt,:,:)=amn(ikpt,:,:)*conjg(nsfactor)
1491 : ENDIF
1492 :
1493 :
1494 :
1495 : !****************************************************************
1496 : !... vanderbilt mmn matrix
1497 : !***************************************************************
1498 8 : IF (wann%l_matrixmmn .AND.&
1499 : (.NOT.wann%l_skipkov)) THEN ! vanderbilt procedure Mmn matrix
1500 8 : ALLOCATE ( we_b(input%neig) )
1501 :
1502 : !!! the cycle by the nearest neighbors (nntot) for each kpoint
1503 :
1504 72 : DO ikpt_b = 1,nntot
1505 :
1506 64 : IF(pair_to_do(ikpt,ikpt_b).EQ.0)CYCLE !save time by symmetry
1507 32 : kptibz_b=bpt(ikpt_b,ikpt)
1508 32 : IF(wann%l_bzsym) oper_b=mapkoper(kptibz_b)
1509 32 : IF (wann%l_bzsym) kptibz_b=irreduc(kptibz_b)
1510 :
1511 :
1512 :
1513 : ! now we need the wavefunctions for k_b kpoint
1514 :
1515 : ! print*,"something to do"
1516 :
1517 :
1518 32 : n_start=1
1519 32 : n_end=input%neig
1520 32 : call lapw_b%init(input,noco,nococonv,kpts,atoms,sym,kptibz_b,cell,fmpi)
1521 :
1522 :
1523 32 : nbasfcn_b = MERGE(lapw_b%nv(1)+lapw_b%nv(2)+2*atoms%nlotot,lapw_b%nv(1)+atoms%nlotot,noco%l_noco)
1524 32 : CALL zMat_b%init(l_real,nbasfcn_b,input%neig)
1525 32 : CALL zzMat%init(l_real,nbasfcn_b,input%neig)
1526 :
1527 :
1528 : CALL cdn_read(&
1529 : eig_id,&
1530 : lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd?&
1531 : kptibz_b,jspin,nbasfcn_b,&
1532 : noco%l_ss,noco%l_noco,input%neig,n_start,n_end,&
1533 32 : nbands_b,eigg,zzMat)
1534 :
1535 :
1536 32 : nslibd_b = 0
1537 :
1538 288 : eig_b(:) = 0.
1539 :
1540 288 : DO i = 1,nbands_b
1541 : IF((eigg(i).GE.sliceplot%e1s.AND.nslibd_b.LT.numbands&
1542 : .AND.wann%l_bynumber).OR.&
1543 : (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
1544 256 : wann%l_byenergy).OR.(i.GE.wann%band_min(jspin2).AND.&
1545 : (i.LE.wann%band_max(jspin2)).AND.&
1546 32 : wann%l_byindex))THEN
1547 256 : nslibd_b = nslibd_b + 1
1548 256 : eig_b(nslibd_b) = eigg(i)
1549 256 : we_b(nslibd_b) = we_b(i)
1550 256 : IF (zzMat%l_real) THEN
1551 23552 : zMat_b%data_r(:,nslibd_b) = zzMat%data_r(:,i)
1552 : ELSE
1553 0 : zMat_b%data_c(:,nslibd_b) = zzMat%data_c(:,i)
1554 : END IF
1555 : ENDIF
1556 : ENDDO
1557 :
1558 : !***********************************************************
1559 : ! Rotate the wavefunction of next neighbor.
1560 : !***********************************************************
1561 32 : IF (wann%l_bzsym .AND. (oper_b.NE.1) ) THEN
1562 : call wann_kptsrotate(&
1563 : atoms, &
1564 : sym%invsat,&
1565 : noco%l_noco,noco%l_soc,&
1566 : jspin,&
1567 : oper_b,sym%nop,sym%mrot,lapw_b%dim_nvd(),&
1568 : shiftkpt(:,bpt(ikpt_b,ikpt)),&
1569 : sym%tau,&
1570 : lapw_b,&
1571 0 : zMat_b,nsfactor_b)
1572 : ELSE
1573 32 : nsfactor_b=CMPLX(1.0,0.0)
1574 : ENDIF
1575 : ! print*,"kpt2=",bkpt_b
1576 :
1577 32 : noccbd_b = nslibd_b
1578 :
1579 : !cccc we start with the Mmn matrix ccccccccccccc
1580 :
1581 : !!! matrix elements of the interstitial overlap
1582 : !!! matrix with the weights given by products of
1583 : !!! the c-coeff. for different G-vectors, bands and k-points
1584 : !!! Mmn(k,b)(IR) = \sum(G,G')C_G^(k,n)*C_G'^(k+b,m)\theta_(G-G')
1585 :
1586 : !ccccccccccc Spherical Contributions ccccccccccccc
1587 :
1588 0 : ALLOCATE (acof_b(noccbd_b,0:lmd,atoms%nat),&
1589 0 : bcof_b(noccbd_b,0:lmd,atoms%nat),&
1590 0 : ccof_b(-atoms%llod:atoms%llod,noccbd_b,atoms%nlod,&
1591 224 : atoms%nat))
1592 :
1593 : !!! get the band-dependent k-dependent ab coeff.
1594 :
1595 : CALL abcof(input,atoms,sym,cell,lapw_b,&
1596 : noccbd_b,usdus,noco,nococonv,jspin2 ,&
1597 32 : acof_b,bcof_b,ccof_b,zMat_b)
1598 :
1599 :
1600 : CALL wann_abinv(atoms,sym,&
1601 32 : acof_b,bcof_b,ccof_b)
1602 :
1603 : !cccccccc Interstitial ccccccccccccccccccccccccc
1604 : !!! matrix elements of the interstitial overlap
1605 : !!! matrix with the weights given by products of
1606 : !!! the c-coeff. for different G-vectors, bands and k-points
1607 : !!! Mmn(k,b)(IR) = \sum(G,G')C_G^(k,n)*C_G'^(k+b,m)\theta_(G-G')
1608 : !... these overlaps are the same as in the mmk0 case, only differ
1609 : !... by the b-dependence of the C-coefficients
1610 : !cccccccccccccccccccccccccccccccccccccccccccccccccccc
1611 :
1612 32 : addnoco=0
1613 32 : addnoco2=0
1614 32 : noconbasfcn=nbasfcn
1615 32 : noconbasfcn_b=nbasfcn_b
1616 32 : IF(noco%l_noco.AND.(jspin.EQ.2))THEN
1617 0 : addnoco = lapw%nv(1) + atoms%nlotot
1618 0 : addnoco2 = lapw_b%nv(1) + atoms%nlotot
1619 : ! noconbasfcn=nbasfcn-addnoco
1620 : ! noconbasfcn_b=nbasfcn_b-addnoco2
1621 : ENDIF
1622 :
1623 : CALL wann_mmkb_int(&
1624 : cmplx_1,addnoco,addnoco2,&
1625 : lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
1626 : stars%ng3,lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
1627 : lapw%nv(jspin2),input%neig,noconbasfcn,noconbasfcn_b,zMat,nslibd,&
1628 : lapw_b%k1(:,jspin2),lapw_b%k2(:,jspin2),lapw_b%k3(:,jspin2),&
1629 : lapw_b%nv(jspin2),zMat_b,nslibd_b,&
1630 : nbnd,&
1631 : stars%rgphs,stars%ustep,stars%ig,gb(:,ikpt_b,ikpt),&
1632 32 : mmnk(:,:,ikpt_b,ikpt))
1633 :
1634 :
1635 :
1636 : !ccccccccccc Spherical Contributions ccccccccccccc
1637 :
1638 64 : chi = cmplx_1
1639 : CALL wann_mmkb_sph(&
1640 : nbnd,atoms%llod,nslibd,nslibd_b,atoms%nlod,atoms%nat,&
1641 : atoms%ntype,lmd,atoms%jmtd,atoms%taual,sym%nop,atoms%lmax,&
1642 : atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,acof,bcof,ccof,&
1643 : lapw_b%bkpt,acof_b,bcof_b,ccof_b,gb(:,ikpt_b,ikpt),lapw%bkpt,&
1644 : ujug,ujdg,&
1645 : djug,djdg,ujulog,djulog,ulojug,ulojdg,ulojulog,kdiff,&
1646 : nntot,chi,&
1647 32 : mmnk(:,:,ikpt_b,ikpt))
1648 :
1649 :
1650 :
1651 :
1652 : !...vacuum contributions
1653 32 : IF (input%film ) THEN
1654 :
1655 : CALL wann_mmkb_vac(&
1656 : cmplx_1,noco%l_noco,atoms%nlotot,qpt_i,&
1657 : nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
1658 : stars%mx1,stars%mx2,stars%mx3,&
1659 : stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
1660 : vacuum%delz,stars%ig2,cell%area,cell%bmat,&
1661 : cell%bbmat,enpara%evac0(:,jspin2),&
1662 : enpara%evac0(:,jspin2_b),&
1663 : lapw%bkpt,lapw_b%bkpt,vz(:,:,jspin2),vz(:,:,jspin2_b),&
1664 : nslibd,nslibd_b,jspin2,jspin2_b,&
1665 : lapw%k1,lapw%k2,lapw%k3,lapw_b%k1,lapw_b%k2,lapw_b%k3,&
1666 : wannierspin,lapw%dim_nvd(),&
1667 : lapw%dim_nbasfcn(),input%neig,zMat,zMat_b,&
1668 : lapw%nv,lapw_b%nv,cell%omtil,&
1669 : gb(:,ikpt_b,ikpt),&
1670 0 : mmnk(:,:,ikpt_b,ikpt))
1671 :
1672 : ENDIF
1673 :
1674 : ! mmnk(:,:,ikpt_b,ikpt)=
1675 : ! mmnk(:,:,ikpt_b,ikpt)*nsfactor*conjg(nsfactor_b)
1676 :
1677 :
1678 72 : DEALLOCATE ( acof_b,bcof_b,ccof_b )
1679 :
1680 :
1681 : enddo ! end of loop by the nearest k-neighbors
1682 :
1683 8 : DEALLOCATE ( we_b )
1684 :
1685 : ENDIF!l_matrixmmn=.true.
1686 :
1687 : 9900 CONTINUE ! jump for doublespin loop
1688 :
1689 8 : IF (wann%l_matrixmmn) THEN ! vanderbilt procedure Mmn matrix
1690 :
1691 : !*******************************************c
1692 : ! START Q-NEIGHBOR LOOP c
1693 : !*******************************************c
1694 8 : ALLOCATE ( we_qb(input%neig) )
1695 :
1696 8 : DO iqpt_b=1,nntot_q
1697 8 : IF(.NOT.l_gwf) EXIT ! old functionality
1698 :
1699 0 : qptibz_b = bpt_q(iqpt_b,iqpt)
1700 0 : IF(qptibz_b.EQ.qptibz) CYCLE ! no need to compute overlaps
1701 : ! with periodic images for now
1702 :
1703 0 : qptb_i = nococonv%qss
1704 0 : alphb_i = nococonv%alph
1705 0 : betab_i = nococonv%beta
1706 0 : thetab_i = nococonv%theta
1707 0 : phib_i = nococonv%phi
1708 0 : IF(wann%l_sgwf) THEN
1709 0 : qptb_i(:) = wann%param_vec(:,qptibz_b)
1710 0 : alphb_i(:) = wann%param_alpha(:,qptibz_b)
1711 0 : ELSEIF(wann%l_socgwf) THEN
1712 0 : IF(wann%l_dim(2)) phib_i = tpi_const*wann%param_vec(2,qptibz_b)
1713 0 : IF(wann%l_dim(3)) thetab_i = tpi_const*wann%param_vec(3,qptibz_b)
1714 : ENDIF
1715 :
1716 : !if(pair_to_do_q(iqpt,iqpt_b).eq.0)cycle ! TODO: use symmetry
1717 0 : IF(wann%l_bzsym) oper_qb=mapqoper(qptibz_b)
1718 0 : IF (wann%l_bzsym) qptibz_b=irreduc_q(qptibz_b)
1719 :
1720 0 : n_start=1
1721 0 : n_end=input%neig
1722 :
1723 : ! read in diagonalization information from corresponding
1724 : ! eig file to q-point iqpt_b at a given k-point ikpt.
1725 : ! as a check verify that bkpt.eq.bqpt (same k).
1726 : ! moreover, the plane-wave vectors G(k,q+b) are stored
1727 : ! in (k1_qb,k2_qb,k3_qb) for later use.
1728 :
1729 0 : CALL lapw_qb%init(input,noco,nococonv,kpts,atoms,sym,kptibz,cell,fmpi)
1730 : CALL cdn_read(&
1731 : innerEig_idList(iqpt_b),&
1732 : lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd? !kptibz_b2?&
1733 : kptibz,jspin_b,lapw%dim_nbasfcn(),&
1734 : noco%l_ss,noco%l_noco,input%neig,n_start,&
1735 : n_end,&
1736 : nbands_qb,eigg, &
1737 0 : zzMat)
1738 :
1739 :
1740 : ! are we dealing with the same k-point at which
1741 : ! we want to construct A,B,C coefficients etc. ?
1742 0 : IF(ANY(bqpt.NE.lapw%bkpt)) CALL juDFT_error("bqpt.ne.bkpt",&
1743 0 : calledby="wannier")
1744 :
1745 0 : zMat_qb%l_real = zzMat%l_real
1746 0 : zMat_qb%matsize1 = zzMat%matsize1
1747 0 : zMat_qb%matsize2 = zzMat%matsize2
1748 0 : IF (zzMat%l_real) THEN
1749 0 : ALLOCATE (zMat_qb%data_r(zMat%matsize1,zMat%matsize2))
1750 0 : zMat_qb%data_r = 0.0
1751 : ELSE
1752 0 : ALLOCATE (zMat_qb%data_c(zMat%matsize1,zMat%matsize2))
1753 0 : zMat_qb%data_c = CMPLX(0.0,0.0)
1754 : END IF
1755 :
1756 0 : eig_qb(:) = 0.
1757 :
1758 0 : nslibd_qb = 0
1759 0 : DO i = 1,nbands_qb
1760 : IF((eigg(i).GE.sliceplot%e1s.AND.nslibd_qb.LT.numbands&
1761 : .AND.wann%l_bynumber).OR.&
1762 : (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
1763 0 : wann%l_byenergy).OR.(i.GE.wann%band_min(jspin).AND.&
1764 0 : (i.LE.wann%band_max(jspin)).AND.wann%l_byindex)) THEN
1765 0 : nslibd_qb = nslibd_qb + 1
1766 0 : eig_qb(nslibd_qb) = eigg(i)
1767 0 : we_qb(nslibd_qb) = we_qb(i)
1768 0 : IF (zzMat%l_real) THEN
1769 0 : zMat_qb%data_r(:,nslibd_qb) = zzMat%data_r(:,i)
1770 : ELSE
1771 0 : zMat_qb%data_c(:,nslibd_qb) = zzMat%data_c(:,i)
1772 : END IF
1773 : ENDIF
1774 : ENDDO
1775 :
1776 : ! check that eigenvectors and -values are identical if q=q+b
1777 0 : IF(iqpt.EQ.qptibz_b .AND. jspin.EQ.jspin_b) THEN
1778 0 : IF(zMat%l_real) THEN
1779 0 : IF(ANY(zMat%data_r.NE.zMat_qb%data_r)) &
1780 0 : WRITE(*,*)'z.ne.z_qb',iqpt,ikpt
1781 : ELSE
1782 0 : IF(ANY(zMat%data_c.NE.zMat_qb%data_c)) &
1783 0 : WRITE(*,*)'z.ne.z_qb',iqpt,ikpt
1784 : END IF
1785 0 : IF(ANY(eig.NE.eig_qb)) WRITE(*,*)'eig.ne.eiq_qb',iqpt,ikpt
1786 0 : IF(lapw%nv(jspin).NE.lapw_qb%nv(jspin)) WRITE(*,*)'nv!=nv_qb',iqpt,ikpt
1787 : ENDIF
1788 :
1789 : ! check that number of bands are the same at (k,q) and (k,q+b)
1790 0 : IF(nslibd.NE.nslibd_qb)&
1791 0 : WRITE(*,*)'nslibd.ne.nslibd_qb',ikpt,iqpt,iqpt_b
1792 :
1793 :
1794 0 : noccbd_qb = nslibd_qb
1795 0 : nsfactor_b=CMPLX(1.0,0.0)
1796 :
1797 :
1798 0 : ALLOCATE (acof_qb(noccbd_qb,0:lmd,atoms%nat),&
1799 0 : bcof_qb(noccbd_qb,0:lmd,atoms%nat),&
1800 0 : ccof_qb(-atoms%llod:atoms%llod,noccbd_qb,atoms%nlod,&
1801 0 : atoms%nat))
1802 :
1803 0 : acof_qb(:,:,:) = CMPLX(0.,0.) ; bcof_qb(:,:,:) = CMPLX(0.,0.)
1804 0 : ccof_qb(:,:,:,:) = CMPLX(0.,0.)
1805 :
1806 : ! construct the A,B,C coefficients of the wave function
1807 : ! at the point (k,q+b) using previously read information
1808 :
1809 : CALL abcof(input,atoms,sym,cell,lapw_qb,&
1810 : noccbd_qb,usdus,noco,nococonv,jspin_b ,&
1811 0 : acof_qb,bcof_qb,ccof_qb,zMat_qb)
1812 :
1813 : CALL wann_abinv(atoms,sym,&
1814 0 : acof_qb,bcof_qb,ccof_qb)
1815 :
1816 : ! check that A,B,C coefficients are the same if q+b = q
1817 0 : IF(l_gwf.AND.(iqpt.EQ.qptibz_b).AND.(jspin.EQ.jspin_b)) THEN
1818 0 : IF(ANY(acof_qb.NE.acof)) WRITE(*,*)'acof',iqpt,ikpt
1819 0 : IF(ANY(bcof_qb.NE.bcof)) WRITE(*,*)'bcof',iqpt,ikpt
1820 0 : IF(ANY(ccof_qb.NE.ccof)) WRITE(*,*)'ccof',iqpt,ikpt
1821 : ENDIF
1822 :
1823 0 : addnoco=0
1824 0 : addnoco2=0
1825 0 : IF(noco%l_noco.AND.(jspin.EQ.2))THEN
1826 0 : addnoco = lapw%nv(1) + atoms%nlotot
1827 : ENDIF
1828 0 : IF(noco%l_noco.AND.(jspin_b.EQ.2))THEN
1829 0 : addnoco2 = lapw_qb%nv(1) + atoms%nlotot
1830 : ENDIF
1831 :
1832 :
1833 : ! set up local->global transformation for overlaps
1834 0 : DO n=1,atoms%ntype
1835 0 : IF(wann%l_sgwf) THEN
1836 0 : dalph = alph_i(n)-alphb_i(n)
1837 0 : db1 = beta_i(n)/2.
1838 0 : db2 = betab_i(n)/2.
1839 0 : ELSEIF(wann%l_socgwf) THEN
1840 0 : dalph = phi_i-phib_i
1841 0 : db1 = theta_i/2.
1842 0 : db2 = thetab_i/2.
1843 : ENDIF
1844 : coph = COS(dalph)
1845 : siph = SIN(dalph)
1846 : phasfac = CMPLX(coph,siph)
1847 : phasfac2= CMPLX(coph,-siph)
1848 :
1849 0 : IF(l_p0 .AND. dalph.NE.0.0) THEN
1850 0 : WRITE(*,*)'WARNING: include dalph in chi trafo!'
1851 : ENDIF
1852 :
1853 0 : IF( (jspin.EQ.1) .AND. (jspin_b.EQ.1) ) THEN ! uu
1854 : ! chi(n) = cos(db1)*cos(db2)*phasfac
1855 : ! > + sin(db1)*sin(db2)*phasfac2
1856 0 : chi(n) = COS(db2-db1)
1857 0 : ELSEIF( (jspin.EQ.2) .AND. (jspin_b.EQ.2) ) THEN !dd
1858 : ! chi(n) = cos(db1)*cos(db2)*phasfac2
1859 : ! > + sin(db1)*sin(db2)*phasfac
1860 0 : chi(n) = COS(db2-db1)
1861 0 : ELSEIF( (jspin.EQ.1) .AND. (jspin_b.EQ.2) ) THEN ! ud
1862 : ! chi(n) = sin(db1)*cos(db2)*phasfac2
1863 : ! > - cos(db1)*sin(db2)*phasfac
1864 0 : chi(n) = -SIN(db2-db1)
1865 0 : ELSEIF( (jspin.EQ.2) .AND. (jspin_b.EQ.1) ) THEN ! du
1866 : ! chi(n) = cos(db1)*sin(db2)*phasfac2
1867 : ! > - sin(db1)*cos(db2)*phasfac
1868 0 : chi(n) = SIN(db2-db1)
1869 : ELSE
1870 0 : STOP 'problem setting up chi: jspin,jspin_b'
1871 : ENDIF
1872 : ENDDO
1873 0 : chi = CONJG(chi)
1874 : !chi = cmplx_1
1875 :
1876 : ! optional: disable chi transformation
1877 : ! instead of computing overlap w.r.t. global frame
1878 : ! only consider wave function overlaps in local frames
1879 0 : IF(l_nochi) THEN
1880 0 : IF(doublespin.LT.3) THEN
1881 0 : chi = CMPLX(1.0,0.0)
1882 : ELSE
1883 0 : chi = CMPLX(0.0,0.0)
1884 : ENDIF
1885 : ENDIF
1886 :
1887 0 : IF((iqpt.EQ.1).AND.(ikpt.EQ.1).AND.(iqpt_b.EQ.1)) THEN
1888 0 : WRITE(*,*)'dbs',doublespin,'chi',chi(1)
1889 : ENDIF
1890 :
1891 : ! muffin tin contribution to overlap is computed taking into account
1892 : ! the spin-dependent phase exp(+/- tau*b/2) and the ujugaunt integrals
1893 : ! calculated especially for the q-points before. then, q and q+b
1894 : ! take the role of k and k+b and the same for G(q+b) and G(k+b)
1895 0 : IF(wann%l_sgwf) CALL wann_mmkb_sph( &
1896 : nbnd,atoms%llod,nslibd,nslibd_qb,atoms%nlod,atoms%nat,&
1897 : atoms%ntype,lmd,atoms%jmtd,sign_q*atoms%taual/2.0,sym%nop,&
1898 : atoms%lmax,atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,&
1899 : acof,bcof,ccof,qptb_i,&
1900 : acof_qb,bcof_qb,ccof_qb,gb_q(:,iqpt_b,iqpt),qpt_i,&
1901 : ujug_q,ujdg_q,&
1902 : djug_q,djdg_q,ujulog_q,djulog_q,ulojug_q,ulojdg_q,&
1903 : ulojulog_q,qdiff, &
1904 : nntot_q,chi,&
1905 0 : mmnk_q(:,:,iqpt_b,ikpt))
1906 0 : IF(wann%l_socgwf) CALL wann_mmkb_sph( &
1907 : nbnd,atoms%llod,nslibd,nslibd_qb,atoms%nlod,atoms%nat,&
1908 : atoms%ntype,lmd,atoms%jmtd,&
1909 : zero_taual,sym%nop,atoms%lmax, &
1910 : atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,acof,bcof,ccof,&
1911 : (/ 0.0, phib_i/tpi_const, thetab_i/tpi_const /),&
1912 : acof_qb,bcof_qb,ccof_qb,gb_q(:,iqpt_b,iqpt),&
1913 : (/ 0.0, phi_i/tpi_const, theta_i/tpi_const /),&
1914 : ujug_q,ujdg_q,&
1915 : djug_q,djdg_q,ujulog_q,djulog_q,ulojug_q,ulojdg_q,&
1916 : ulojulog_q,qdiff, &
1917 : nntot_q,chi,&
1918 0 : mmnk_q(:,:,iqpt_b,ikpt))
1919 :
1920 :
1921 : IF(((doublespin.NE.3).AND.(doublespin.NE.4))&
1922 0 : .OR.(.NOT.noco%l_noco)) THEN
1923 :
1924 0 : IF(.NOT.noco%l_noco)THEN
1925 0 : interchi=chi(1)
1926 0 : vacchi=chi(1)
1927 : ELSE
1928 0 : interchi=cmplx_1
1929 0 : vacchi=cmplx_1
1930 : ENDIF
1931 :
1932 : ! interstitial contribution to overlap is computed using
1933 : ! (-/+ 1)*G(q+b)/2 as G-vector connecting neighbors
1934 : ! and lattice vectors G(k,q) (k1...) and G(k,q+b) (k1_qb...)
1935 0 : IF(wann%l_sgwf) CALL wann_mmkb_int(&
1936 : interchi,addnoco,addnoco2,&
1937 : lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
1938 : stars%ng3,lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
1939 : lapw%nv(jspin),input%neig,nbasfcn,nbasfcn_b,zMat,nslibd,&
1940 : lapw_qb%k1(:,jspin_b),lapw_qb%k2(:,jspin_b),lapw_qb%k3(:,jspin_b),&
1941 : lapw_qb%nv(jspin_b),zMat_qb,nslibd_qb,&
1942 : nbnd,&
1943 : stars%rgphs,stars%ustep,stars%ig,&
1944 : sign_q*gb_q(:,iqpt_b,iqpt)/2, &
1945 0 : mmnk_q(:,:,iqpt_b,ikpt))
1946 0 : IF(wann%l_socgwf) CALL wann_mmkb_int(&
1947 : interchi,addnoco,addnoco2,&
1948 : lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
1949 : stars%ng3,lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
1950 : lapw%nv(jspin),input%neig,nbasfcn,nbasfcn_b,zMat,nslibd,&
1951 : lapw_qb%k1(:,jspin_b),lapw_qb%k2(:,jspin_b),lapw_qb%k3(:,jspin_b),&
1952 : lapw_qb%nv(jspin_b),zMat_qb,nslibd_qb,&
1953 : nbnd,&
1954 : stars%rgphs,stars%ustep,stars%ig,(/ 0, 0, 0 /), &
1955 0 : mmnk_q(:,:,iqpt_b,ikpt))!m_int(:,:,iqpt_b,ikpt))
1956 :
1957 :
1958 : ! vacuum contribution in film calculation
1959 0 : IF (input%film ) THEN
1960 0 : IF(wann%l_sgwf) CALL wann_mmkb_vac(&
1961 : vacchi,noco%l_noco,atoms%nlotot,sign_q*2.*lapw%bkpt,&
1962 : nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
1963 : stars%mx1,stars%mx2,stars%mx3,&
1964 : stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
1965 : vacuum%delz,stars%ig2,cell%area,cell%bmat,&
1966 : cell%bbmat,enpara%evac0(:,jspin),enpara%evac0(:,jspin_b),&
1967 : sign_q*qpt_i/2.,&
1968 : sign_q*qptb_i/2.,&
1969 : vz(:,:,jspin2),vz(:,:,jspin2_b),&
1970 : nslibd,nslibd_qb,jspin,jspin_b,&
1971 : lapw%k1,lapw%k2,lapw%k3,lapw_qb%k1,lapw_qb%k2,lapw_qb%k3,&
1972 : wannierspin,lapw%dim_nvd(),&
1973 : lapw%dim_nbasfcn(),input%neig,zMat,zMat_qb,lapw%nv,&
1974 : lapw_qb%nv,cell%omtil,&
1975 : sign_q*gb_q(:,iqpt_b,iqpt)/2,&
1976 0 : mmnk_q(:,:,iqpt_b,ikpt))
1977 0 : IF(wann%l_socgwf) CALL wann_mmkb_vac(&
1978 : vacchi,noco%l_noco,atoms%nlotot,qpt_i,&
1979 : nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
1980 : stars%mx1,stars%mx2,stars%mx3,&
1981 : stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
1982 : vacuum%delz,stars%ig2,cell%area,cell%bmat,&
1983 : cell%bbmat,enpara%evac0(:,jspin),enpara%evac0(:,jspin_b),&
1984 : bqpt,bqpt,&
1985 : vz(:,:,jspin2),vz(:,:,jspin2_b),&
1986 : nslibd,nslibd_qb,jspin,jspin_b,&
1987 : lapw%k1,lapw%k2,lapw%k3,lapw_qb%k1,lapw_qb%k2,lapw_qb%k3,&
1988 : wannierspin,lapw%dim_nvd(),&
1989 : lapw%dim_nbasfcn(),input%neig,zMat,zMat_qb,lapw%nv,&
1990 : lapw_qb%nv,cell%omtil,&
1991 : (/ 0, 0, 0 /),&
1992 0 : mmnk_q(:,:,iqpt_b,ikpt))
1993 :
1994 : ! vacuum contribution in one-dimensional chain calculation where
1995 : ! q-point plays role of k-point with proper prefactor of (-/+ 1/2).
1996 : ! moreover, the k-point ikpt is treated like qss in the subroutine
1997 : ! such that a correction for sign and factor has to appear as well.
1998 : ! lattice vectors G(k,q) (k1...) and G(k,q+b) (k1_qb...) need to
1999 : ! be provided.
2000 : ENDIF!film resp. odi
2001 :
2002 :
2003 : ENDIF!doublespin
2004 :
2005 8 : DEALLOCATE ( acof_qb,bcof_qb,ccof_qb )
2006 :
2007 :
2008 : ENDDO !iqpt_b, q-neighbors
2009 : !**************************************************c
2010 : ! END Q-NEIGHBOR LOOP c
2011 : !**************************************************c
2012 :
2013 8 : DEALLOCATE ( we_qb )
2014 :
2015 : ENDIF ! if wann%l_matrixmmn = true
2016 :
2017 8 : IF(.NOT.wann%l_bzsym)oper=0
2018 8 : IF(.NOT.wann%l_plot_symm.OR.oper.EQ.1)THEN
2019 8 : IF (wann%l_wann_plot .AND. &
2020 : (doublespin.EQ.1 .OR. doublespin.EQ.2)) THEN
2021 :
2022 0 : addnoco=0
2023 0 : IF(noco%l_noco.AND.(jspin.EQ.2))THEN
2024 0 : addnoco=lapw%nv(1)+atoms%nlotot
2025 : ENDIF
2026 :
2027 0 : IF (sliceplot%slice) THEN
2028 0 : IF (ikpt.EQ.sliceplot%kk) THEN
2029 :
2030 0 : WRITE (oUnit,*) 'nnne=',sliceplot%nnne
2031 0 : WRITE (oUnit,*) 'eig(nnne)=',eig(sliceplot%nnne)
2032 0 : WRITE (oUnit,*) 'we(nnne)=',we(sliceplot%nnne)
2033 :
2034 : CALL wann_plot(&
2035 : vacuum,stars,cell,atoms,&
2036 : lapw%dim_nv2d(),jspin,stars%ng3,&
2037 : vacuum%nmzxyd,&
2038 : stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
2039 : atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,&
2040 : stars%ng3,vacuum%nvac,vacuum%nmz,vacuum%nmzxy,stars%ng2,&
2041 : sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,&
2042 : sym%symor,sym%invs,sym%invs2,cell%z1,vacuum%delz,&
2043 : sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,&
2044 : atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,&
2045 : cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,&
2046 : sliceplot%kk,lapw%dim_nvd(),atoms%nlod,atoms%llod,&
2047 : lapw%nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,&
2048 : lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),enpara%evac0(:,jspin),&
2049 : vz(:,:,jspin2),&
2050 : nslibd,lapw%dim_nbasfcn(),input%neig,&
2051 : ff(:,:,:,:,jspin),&
2052 : gg(:,:,:,:,jspin),flo,acof,bcof,ccof,zMat,&
2053 : stars%mx1,stars%mx2,stars%mx3,stars%ig,stars%ig2,&
2054 : stars%sk2,stars%phi2,&
2055 : noco%l_noco,noco%l_ss,qpt_i,&
2056 0 : addnoco,get_index_kq(ikpt,iqpt,fullnkpts),wann%l_sgwf)
2057 :
2058 : ENDIF
2059 : ELSE ! not sliceplot%slice
2060 :
2061 : CALL wann_plot(&
2062 : vacuum,stars,cell,atoms,&
2063 : lapw%dim_nv2d(),jspin ,stars%ng3,&
2064 : vacuum%nmzxyd,&
2065 : stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
2066 : atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,&
2067 : stars%ng3,vacuum%nvac,vacuum%nmz,vacuum%nmzxy,stars%ng2,&
2068 : sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,&
2069 : sym%symor,sym%invs,sym%invs2,cell%z1,vacuum%delz,&
2070 : sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,&
2071 : atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,&
2072 : cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,&
2073 : sliceplot%kk,lapw%dim_nvd(),atoms%nlod,atoms%llod,&
2074 : lapw%nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,&
2075 : lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),enpara%evac0(:,jspin),&
2076 : vz(:,:,jspin2),&
2077 : nslibd,lapw%dim_nbasfcn(),input%neig,&
2078 : ff(:,:,:,:,jspin),&
2079 : gg(:,:,:,:,jspin),flo,acof,bcof,ccof,zMat,&
2080 : stars%mx1,stars%mx2,stars%mx3,stars%ig,stars%ig2,&
2081 : stars%sk2,stars%phi2,&
2082 : noco%l_noco,noco%l_ss,qpt_i,&
2083 0 : addnoco,get_index_kq(ikpt,iqpt,fullnkpts),wann%l_sgwf)
2084 :
2085 0 : IF(wann%l_plot_symm.AND.wann%l_bzsym)THEN
2086 0 : DO kplot=1,fullnkpts
2087 0 : IF(irreduc(kplot).EQ.kptibz)THEN
2088 0 : plotoper=mapkoper(kplot)
2089 0 : IF(plotoper.LT.0)THEN
2090 0 : plotoper=-plotoper
2091 0 : l_conjugate=.TRUE.
2092 : ELSE
2093 0 : l_conjugate=.FALSE.
2094 : ENDIF
2095 0 : kplotoper=sym%invtab(plotoper)
2096 : CALL wann_plot_symm(jspin,sym%mrot(:,:,kplotoper),ikpt,&
2097 0 : kplot,l_conjugate)
2098 : ENDIF
2099 : ENDDO
2100 : ENDIF
2101 :
2102 : ENDIF
2103 : ENDIF
2104 : ENDIF !wann%l_plot_symm
2105 8 : DEALLOCATE ( acof,bcof,ccof,we,eigg )
2106 :
2107 8 : IF(wann%l_projmethod.OR.wann%l_bestproj)THEN
2108 : CALL wann_projmethod(&
2109 : fullnkpts,&
2110 : wann%l_projmethod,wann%l_bestproj,&
2111 : ikpt,nwfs,nslibd,amn,eig,&
2112 0 : psiw,hwfr)
2113 : ENDIF ! projmethod
2114 :
2115 : ENDIF ! loop by processors
2116 :
2117 : ENDDO ! end of cycle by the k-points
2118 :
2119 :
2120 2 : IF(wann%l_matrixmmn)&
2121 0 : DEALLOCATE(ujug,ujdg,djug,djdg,&
2122 2 : ujulog,djulog,ulojulog,ulojug,ulojdg)
2123 2 : IF(wann%l_matrixmmn.AND.l_gwf)&
2124 0 : DEALLOCATE(ujug_q,ujdg_q,djug_q,&
2125 0 : djdg_q,ujulog_q,djulog_q,ulojulog_q,ulojug_q,&
2126 0 : ulojdg_q)
2127 :
2128 : #ifdef CPP_MPI
2129 2 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
2130 : #endif
2131 :
2132 :
2133 : !******************************************************
2134 : ! Write down the projections.
2135 : !******************************************************
2136 : 5 CONTINUE
2137 :
2138 2 : IF(doublespin.EQ.3 .OR. doublespin.EQ.4) GOTO 912
2139 :
2140 2 : IF(wann%l_nabla)THEN
2141 :
2142 0 : nablamat=nablamat*hescale
2143 : CALL wann_write_nabla(&
2144 : fmpi%mpi_comm,l_p0,spin12(jspin)//'.nabl',&
2145 : 'Matrix elements of nabla operator',&
2146 : nbnd,fullnkpts,nbnd,&
2147 : fmpi%irank,fmpi%isize,wann%l_unformatted,&
2148 0 : nablamat)
2149 : ENDIF
2150 :
2151 2 : IF(wann%l_soctomom)THEN
2152 0 : soctomom=soctomom*hescale
2153 : CALL wann_write_nabla(&
2154 : fmpi%mpi_comm,l_p0,spin12(jspin)//'.stm',&
2155 : 'Matrix elements of stm operator',&
2156 : nbnd,fullnkpts,nbnd,&
2157 : fmpi%irank,fmpi%isize,wann%l_unformatted,&
2158 0 : soctomom)
2159 : ENDIF
2160 :
2161 2 : IF(wann%l_surfcurr)THEN
2162 0 : surfcurr=surfcurr*hescale
2163 : CALL wann_write_nabla(&
2164 : fmpi%mpi_comm,l_p0,spin12(jspin)//'.surfcurr',&
2165 : 'Surface currents',&
2166 : nbnd,fullnkpts,nbnd,&
2167 : fmpi%irank,fmpi%isize,wann%l_unformatted,&
2168 0 : surfcurr)
2169 : ENDIF
2170 :
2171 2 : IF((noco%l_soc.OR.noco%l_noco).AND.wann%l_mmn0)THEN
2172 : CALL wann_write_amn(&
2173 : fmpi%mpi_comm,&
2174 : l_p0,spin12(jspin)//'.socmmn0',&
2175 : 'Overlaps of the wavefunct. at the same kpoint',&
2176 : nbnd,fullnkpts,nbnd,&
2177 : fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
2178 0 : mmn,wann%mmn0fmt==2)
2179 : ENDIF !noco%l_soc and l_mmn0
2180 :
2181 2 : IF(wann%l_orbcomp)THEN
2182 0 : num_angl=9
2183 0 : IF(wann%l_oc_f)num_angl=16
2184 : CALL wann_write_matrix5(&
2185 : fmpi%mpi_comm,l_p0,spin12(jspin)//'.orbcomp',&
2186 : 'angular components',&
2187 : nbnd,nbnd,&
2188 : num_angl,wann%oc_num_orbs,fullnkpts,&
2189 : fmpi%irank,fmpi%isize,wann%l_unformatted,&
2190 0 : orbcomp)
2191 : ENDIF
2192 :
2193 2 : IF((noco%l_soc.OR.noco%l_noco) .AND. (doublespin.EQ.1)) THEN
2194 0 : IF(wann%l_mmn0) socmmn(:,:,:)=mmn(:,:,:)
2195 : GOTO 912
2196 : ENDIF
2197 :
2198 2 : IF(noco%l_soc.OR.noco%l_noco) THEN
2199 0 : jspin2=1
2200 0 : IF(wann%l_mmn0) mmn(:,:,:)=socmmn(:,:,:)+mmn(:,:,:)
2201 0 : IF(wann%l_mmn0) DEALLOCATE(socmmn)
2202 : ENDIF
2203 :
2204 2 : IF (wann%l_matrixamn)THEN
2205 : CALL wann_write_amn(&
2206 : fmpi%mpi_comm,&
2207 : l_p0,spin12(jspin2)//'.amn',&
2208 : 'Overlaps of the wavefunct. with the trial orbitals',&
2209 : nbnd,fullnkpts,nwfs,&
2210 : fmpi%irank,fmpi%isize,.FALSE.,.FALSE.,&
2211 2 : amn(:,:,:),wann%matrixamnfmt==2)
2212 : ENDIF !wann%l_matrixamn
2213 :
2214 2 : IF(wann%l_anglmom)THEN
2215 : CALL wann_write_matrix4(&
2216 : fmpi%mpi_comm,&
2217 : l_p0,spin12(jspin2)//'.anglmom',&
2218 : 'Matrix elements of angular momentum',&
2219 : nbnd,nbnd,3,fullnkpts,&
2220 : fmpi%irank,fmpi%isize,&
2221 0 : anglmom)
2222 : ENDIF
2223 :
2224 2 : IF (l_proj) THEN
2225 : !**************************************************************
2226 : ! for projmethod: write down WF1.umn
2227 : !*************************************************************
2228 2 : IF((wann%l_projmethod.OR.wann%l_bestproj))THEN
2229 : CALL wann_write_amn(&
2230 : fmpi%mpi_comm,l_p0,spin12(jspin2)//'.umn',&
2231 : 'transformation to first guess Wannier functions',&
2232 : nbnd,fullnkpts,nwfs,&
2233 : fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
2234 0 : psiw,.FALSE.)
2235 : #ifdef CPP_MPI
2236 0 : ALLOCATE( hwfr2(SIZE(hwfr,1),SIZE(hwfr,2)) )
2237 0 : length=nwfs*nwfs
2238 : CALL MPI_REDUCE(&
2239 : hwfr,hwfr2,length,&
2240 : MPI_DOUBLE_COMPLEX,MPI_SUM,0,&
2241 0 : fmpi%mpi_comm,ierr(1))
2242 0 : hwfr=hwfr2
2243 0 : DEALLOCATE(hwfr2)
2244 : #endif
2245 : !********************************************************
2246 : ! projmethod: hamiltonian matrix in real space
2247 : !********************************************************
2248 0 : IF(l_p0)THEN
2249 0 : WRITE (oUnit,*) 'the hamiltonian matrix in real space:'
2250 0 : DO i = 1,nwfs
2251 0 : DO j = 1,nwfs
2252 0 : WRITE (oUnit,*) ' WFs:',i,'and',j
2253 0 : WRITE (oUnit,*) ' matrix element:',hwfr(i,j)
2254 : ENDDO
2255 : ENDDO
2256 : ENDIF !l_p0
2257 0 : DEALLOCATE(hwfr)
2258 : ENDIF !wann%l_projmethod or wann%l_bestproj
2259 : ENDIF !l_proj
2260 :
2261 : !*********************************************************
2262 : !.....write down the mmn0 matrix
2263 : !*********************************************************
2264 :
2265 2 : IF(wann%l_mmn0)THEN
2266 : CALL wann_write_amn(&
2267 : fmpi%mpi_comm,&
2268 : l_p0,spin12(jspin2)//'.mmn0',&
2269 : 'Overlaps of the wavefunct. at the same kpoint',&
2270 : nbnd,fullnkpts,nbnd,&
2271 : fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
2272 2 : mmn,wann%mmn0fmt==2)
2273 : ENDIF !wann%l_mmn0
2274 :
2275 :
2276 : !*****************************************************
2277 : !.....write down the matrix M^{k,b}_{mn}
2278 : !*****************************************************
2279 :
2280 : ! 912 continue
2281 2 : IF(wann%l_matrixmmn.AND.(.NOT.wann%l_skipkov))THEN
2282 : CALL wann_write_mmnk(&
2283 : fmpi%mpi_comm,jspin2,l_p0,fullnkpts,nntot,wann,&
2284 : maptopair,pair_to_do,nbnd,bpt,gb,&
2285 : fmpi%isize,fmpi%irank," ",&
2286 2 : mmnk,wann%matrixmmnfmt==2) ! wann%l_unformatted)
2287 : ENDIF !wann%l_matrixmmn
2288 :
2289 : 912 CONTINUE
2290 :
2291 2 : IF(l_gwf .AND. wann%l_matrixmmn) THEN
2292 : ! mmnk_q = mmnk_q + m_sph+m_int+m_vac
2293 0 : IF(doublespin.EQ.doublespin_max) THEN
2294 0 : WRITE(fname,'("param_",i4.4,".mmn")')iqpt
2295 : CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
2296 : nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
2297 : fmpi%isize,fmpi%irank,fname,mmnk_q,&
2298 0 : wann%l_unformatted)
2299 : ENDIF
2300 :
2301 : IF(.FALSE.) THEN
2302 : WRITE(fname,'("param_",i4.4,"_",i1,".mmn")')iqpt,doublespin
2303 : CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
2304 : nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
2305 : fmpi%isize,fmpi%irank,fname,&
2306 : m_sph+m_int+m_vac,wann%l_unformatted)
2307 :
2308 : WRITE(fname,'("param_",i4.4,"_",i1,"_int.mmn")')iqpt,doublespin
2309 : CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
2310 : nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
2311 : fmpi%isize,fmpi%irank,fname,m_int,&
2312 : wann%l_unformatted)
2313 : WRITE(fname,'("param_",i4.4,"_",i1,"_sph.mmn")')iqpt,doublespin
2314 : CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
2315 : nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
2316 : fmpi%isize,fmpi%irank,fname,m_sph,&
2317 : wann%l_unformatted)
2318 : WRITE(fname,'("param_",i4.4,"_",i1,"_vac.mmn")')iqpt,doublespin
2319 : CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
2320 : nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
2321 : fmpi%isize,fmpi%irank,fname,m_vac,&
2322 : wann%l_unformatted)
2323 : ENDIF
2324 : ! m_int = cmplx(0.,0.)
2325 : ! m_sph = cmplx(0.,0.)
2326 : ! m_vac = cmplx(0.,0.)
2327 :
2328 : ENDIF
2329 :
2330 :
2331 2 : IF( ALLOCATED (nablamat) ) DEALLOCATE( nablamat )
2332 2 : IF( ALLOCATED (soctomom) ) DEALLOCATE( soctomom )
2333 :
2334 2 : IF( ALLOCATED (surfcurr) ) DEALLOCATE( surfcurr )
2335 2 : IF( ALLOCATED( mmn ) ) DEALLOCATE(mmn)
2336 2 : IF( ALLOCATED( amn ) ) THEN
2337 2 : IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE(amn)
2338 : ENDIF
2339 2 : IF ( ALLOCATED (psiw) ) DEALLOCATE ( psiw )
2340 2 : IF (wann%l_matrixmmn) THEN
2341 2 : IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE (mmnk)
2342 : ENDIF
2343 2 : IF (wann%l_anglmom .AND. ALLOCATED(anglmom))THEN
2344 0 : IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE (anglmom)
2345 : ENDIF
2346 2 : DEALLOCATE(flo)
2347 :
2348 4 : IF(.NOT.noco%l_noco)nrec=nrec+nkpts
2349 :
2350 : !#ifdef CPP_MPI
2351 : ! call MPI_BARRIER(fmpi%mpi_comm,ierr(1))
2352 : !#endif
2353 :
2354 : ENDDO ! end of cycle by spins
2355 :
2356 :
2357 : #ifdef CPP_MPI
2358 2 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
2359 : #endif
2360 :
2361 : ! close eig files
2362 2 : IF (l_gwf) THEN
2363 : ! CALL close_eig(eig_id)
2364 0 : IF(wann%l_matrixmmn)THEN
2365 0 : DO iqpt_b=1,nntot_q
2366 : ! CALL close_eig(innerEig_idList(iqpt_b))
2367 : ENDDO
2368 : ENDIF
2369 : ENDIF
2370 :
2371 2 : IF (wann%l_matrixmmn.AND.ALLOCATED(mmnk))THEN
2372 0 : DEALLOCATE ( mmnk )
2373 : ENDIF
2374 :
2375 2 : IF(ALLOCATED(mmnk_q)) DEALLOCATE(mmnk_q)
2376 : IF(ALLOCATED(m_int)) DEALLOCATE(m_int)
2377 : IF(ALLOCATED(m_sph)) DEALLOCATE(m_sph)
2378 : IF(ALLOCATED(m_vac)) DEALLOCATE(m_vac)
2379 :
2380 : IF ((wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn)&
2381 2 : .AND.ALLOCATED(amn))THEN
2382 0 : DEALLOCATE ( amn )
2383 : ENDIF
2384 :
2385 2 : IF(wann%l_anglmom.AND.ALLOCATED(anglmom))THEN
2386 0 : DEALLOCATE ( anglmom )
2387 : ENDIF
2388 :
2389 : IF ((wann%l_projmethod.OR.wann%l_bestproj)&
2390 2 : .AND.ALLOCATED(hwfr)) THEN
2391 0 : DEALLOCATE ( hwfr )
2392 : ENDIF
2393 :
2394 4 : DEALLOCATE(innerEig_idList)
2395 :
2396 : ENDDO ! iqpt, q-points
2397 : !************************************************c
2398 : ! END Q LOOP c
2399 : !************************************************c
2400 :
2401 2 : IF(ALLOCATED(pair_to_do))DEALLOCATE(pair_to_do,maptopair)
2402 :
2403 2 : DEALLOCATE ( vr,vz)
2404 2 : DEALLOCATE ( ff,gg )
2405 2 : IF (wann%l_bzsym) DEALLOCATE(irreduc,mapkoper)
2406 2 : IF (wann%l_bzsym.AND.l_gwf) DEALLOCATE(irreduc_q,mapqoper)
2407 2 : IF(ALLOCATED(pair_to_do_q))&
2408 0 : DEALLOCATE(pair_to_do_q,maptopair_q)
2409 :
2410 2 : IF (ALLOCATED(kdiff)) DEALLOCATE ( kdiff )
2411 2 : IF (ALLOCATED(qdiff)) DEALLOCATE(qdiff,zero_qdiff)
2412 :
2413 :
2414 : 9110 CONTINUE ! jump for l_finishgwf=T
2415 :
2416 : ! correct for previously introduced factor of 2 in the
2417 : ! G-vectors connecting neighbors across the BZ boundary
2418 2 : IF(wann%l_sgwf) gb_q = gb_q/2
2419 2 : IF(wann%l_socgwf) gb_q = gb_q/2
2420 :
2421 : ! set up input files for wannier90 --> HDWFs
2422 : IF(l_p0.AND.l_gwf.AND.(wann%l_matrixmmn.OR.wann%l_matrixamn)&
2423 2 : .AND.(.NOT.wann%l_skipkov)) THEN
2424 : CALL wann_gwf_commat(fullnkpts,nntot,bpt,fullnqpts,&
2425 : nntot_q,bpt_q,gb,gb_q,wann%aux_latt_const,&
2426 : wann%l_unformatted,wann%l_matrixamn,&
2427 : wann%l_matrixmmn,wann%l_dim,&
2428 0 : wann%nparampts,wann%param_vec/2.0)
2429 : ENDIF
2430 :
2431 2 : IF(l_p0.AND.l_gwf.AND.wann%l_anglmom) THEN
2432 0 : CALL wann_gwf_anglmom(fullnkpts,fullnqpts,wann%l_unformatted)
2433 : ENDIF
2434 :
2435 2 : IF (wann%l_matrixmmn) THEN
2436 2 : DEALLOCATE (gb,bpt)
2437 2 : IF(l_gwf) DEALLOCATE (gb_q,bpt_q)
2438 : ENDIF
2439 :
2440 :
2441 : 1911 CONTINUE
2442 :
2443 4 : IF(ALLOCATED(chi)) DEALLOCATE(chi)
2444 :
2445 4 : CALL timeStop("Wannier total")
2446 : CALL wann_postproc(&
2447 : stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
2448 : lapw ,noco,nococonv,cell,vTot,enpara,sliceplot,eig_id,l_real,& !eig_id is used here after closing the files?!&
2449 4 : wann,fullnkpts,l_proj,results%ef,wann%l_sgwf,fullnqpts)
2450 :
2451 : #ifdef CPP_MPI
2452 4 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
2453 : #endif
2454 :
2455 4 : IF(.NOT.wann%l_ldauwan) THEN
2456 4 : DO pc = 1, wann%nparampts
2457 4 : CALL close_eig(eig_idList(pc))
2458 : END DO
2459 :
2460 4 : CALL juDFT_end("wannier good",fmpi%irank)
2461 : END IF
2462 :
2463 0 : END SUBROUTINE wannier
2464 :
2465 :
2466 : END MODULE m_wannier
|
