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