Line data Source code
1 : !!****m* ABINIT/m_plowannier
2 : !! NAME
3 : !! m_plowannier
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon,AGerossier,ROuterovitch)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! NOTES
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_plowannier
24 :
25 :
26 : #ifndef HAVE_CRPA_OPTIM
27 : #ifdef FC_INTEL
28 : #warning "optimization of m_plowannier is deactivated on intel fortran"
29 : !DEC$ NOOPTIMIZE
30 : #endif
31 : #endif
32 :
33 : use defs_basis
34 : use m_errors
35 : use m_abicore
36 : use m_dtset
37 : use m_dtfil
38 : use defs_wvltypes
39 : use m_xmpi
40 :
41 : use defs_datatypes, only : pseudopotential_type
42 : use defs_abitypes, only : MPI_type
43 : use m_io_tools, only : open_file
44 : use m_mpinfo, only : proc_distrb_cycle
45 : use m_crystal, only : crystal_t
46 : use m_pawtab, only : pawtab_type
47 : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_get,pawcprj_free
48 : use m_pawrad, only : pawrad_type, simp_gen
49 :
50 : implicit none
51 :
52 : private
53 :
54 : public :: init_plowannier
55 : public :: copy_orbital
56 : public :: compute_coeff_plowannier
57 : public :: destroy_plowannier
58 : public :: print_plowannier
59 : public :: get_plowannier
60 : public :: fullbz_plowannier
61 : public :: initialize_operwan
62 : public :: destroy_operwan
63 : public :: zero_operwan
64 : public :: compute_oper_ks2wan
65 : public :: normalization_plowannier
66 : public :: print_operwan
67 : public :: init_operwan_realspace
68 : public :: reduce_operwan_realspace
69 : public :: destroy_operwan_realspace
70 : public :: zero_operwan_realspace
71 : public :: compute_oper_wank2realspace
72 : !!***
73 :
74 :
75 : !!****t* m_plowannier/latom_wan_type
76 : !! NAME
77 : !! latom_wan_type
78 : !!
79 : !! FUNCTION
80 : !!
81 : !!
82 : !! SOURCE
83 :
84 : type, public :: latom_wan_type
85 :
86 : integer, allocatable :: lcalc(:)
87 : ! array of the l we want to compute the psichi with
88 :
89 : end type latom_wan_type
90 : !!***
91 :
92 :
93 : !!****t* m_plowannier/projector_wan_type
94 : !! NAME
95 : !! projector_wan_type
96 : !!
97 : !! FUNCTION
98 : !!
99 : !!
100 : !! SOURCE
101 :
102 : type, public :: projector_wan_type
103 :
104 : integer, allocatable :: lproj(:)
105 : ! gives the list of the projector chosen
106 :
107 : end type projector_wan_type
108 : !!***
109 :
110 :
111 : !!****t* m_plowannier/position_wan_type
112 : !! NAME
113 : !! position_wan_type
114 : !!
115 : !! FUNCTION
116 : !!
117 : !!
118 : !! SOURCE
119 :
120 : type, public :: position_wan_type
121 :
122 : integer, allocatable :: pos(:,:)
123 : ! size (number of position,3)
124 :
125 : end type position_wan_type
126 : !!***
127 :
128 : !!****t* m_plowannier/lorbital_type
129 : !! NAME
130 : !! lorbital_type
131 : !!
132 : !! FUNCTION
133 : !!
134 : !!
135 : !! SOURCE
136 :
137 : type, public :: lorbital_type
138 :
139 : complex(dp), allocatable :: matl(:,:,:)
140 : !details for different m
141 :
142 : real(dp), allocatable :: ph0phiint(:)
143 : ! stocks the values for each projector of the l considered
144 : ! size(total number of projectors for this l)
145 :
146 : end type lorbital_type
147 : !!***
148 :
149 : !!****t* m_plowannier/orbital_type
150 : !! NAME
151 : !! orbital_type
152 : !!
153 : !! FUNCTION
154 : !!
155 : !!
156 : !! SOURCE
157 :
158 : type, public :: orbital_type
159 :
160 : type(lorbital_type), allocatable :: atom(:)
161 : ! details of the psichi coefficients for each atom
162 : ! size of number of l chosen
163 :
164 : end type orbital_type
165 : !!***
166 :
167 :
168 : !!****t* m_plowannier/lorbital2_type
169 : !! NAME
170 : !! lorbital2_type
171 : !!
172 : !! FUNCTION
173 : !!
174 : !!
175 : !! SOURCE
176 :
177 : type, public :: lorbital2_type
178 :
179 : complex(dp), allocatable :: matl(:,:,:,:,:)
180 : ! size (2l1+1,2l2+1,nspppol,nspinor,nspinor)
181 :
182 : real(dp), allocatable :: ph0phiint(:)
183 : ! size (nproj), stocks the value of ph0phiint we may want
184 :
185 :
186 : end type lorbital2_type
187 : !!***
188 :
189 : !!****t* m_plowannier/operwan_type
190 : !! NAME
191 : !! operwan_type
192 : !!
193 : !! FUNCTION
194 : !!
195 : !!
196 : !! SOURCE
197 :
198 : type, public :: operwan_type
199 :
200 : type(lorbital2_type), allocatable :: atom(:,:)
201 : ! l chosen for each on of both atoms
202 :
203 : end type operwan_type
204 : !!***
205 :
206 : !!****t* m_plowannier/atom_index_type
207 : !! NAME
208 : !! atom_index_type_type
209 : !!
210 : !! FUNCTION
211 : !!
212 : !!
213 : !! SOURCE
214 :
215 : type, public :: atom_index_type
216 :
217 : type(operwan_type), allocatable :: position(:,:)
218 : ! size (number of positions chosen for atom1, number of positions chosen for atom2)
219 :
220 : end type atom_index_type
221 : !!***
222 :
223 : !!****t* m_plowannier/operwan_realspace_type
224 : !! NAME
225 : !! operwan_realspace_type
226 : !!
227 : !! FUNCTION
228 : !!
229 : !!
230 : !! SOURCE
231 :
232 : type, public :: operwan_realspace_type
233 :
234 : type(atom_index_type), allocatable :: atom_index(:,:)
235 : ! size (number of atom, number of atom)
236 :
237 : end type operwan_realspace_type
238 : !!***
239 :
240 : !!****t* m_plowannier/plowannier_type
241 : !! NAME
242 : !! plowannier_type
243 : !!
244 : !! FUNCTION
245 : !!
246 : !!
247 : !! SOURCE
248 :
249 : type, public :: plowannier_type
250 :
251 : integer :: nkpt
252 : ! number of k points in Brillouin zone
253 :
254 : integer :: bandi_wan
255 : ! energy band minimum considered
256 :
257 : integer :: bandf_wan
258 : ! energy band maximum considered
259 :
260 : integer :: natom_wan
261 : ! number of atoms (used to compute Wannier functions)
262 :
263 : integer :: size_wan
264 : ! sum of all the m possible for every atom considered
265 :
266 : integer, allocatable :: iatom_wan(:)
267 : ! array of each atom (used to compute Wannier functions)
268 :
269 : integer, allocatable :: nbl_atom_wan(:)
270 : ! array of the number of l considered for each atom
271 :
272 : type(latom_wan_type), allocatable :: latom_wan(:)
273 : ! for each atom, it contains an array of the l we are interested in
274 :
275 : integer, allocatable :: nbproj_atom_wan(:)
276 : ! array of the number of projectors considered for each atom
277 :
278 : type(projector_wan_type), allocatable :: projector_wan(:)
279 : ! for each atom, it contains an array of the projectors we are interested in
280 :
281 : type(position_wan_type), allocatable :: nposition(:)
282 : ! array of the number of position considered for each atom
283 :
284 : integer :: nsppol
285 : ! number of polarization
286 :
287 : integer :: nspinor
288 : ! number of spinorial components
289 :
290 : type(orbital_type), allocatable :: psichi(:,:,:)
291 : ! arrays of psichi
292 :
293 : integer, allocatable :: position(:,:)
294 : ! size natom,3, gives the position of the cell for this atom (rprim coordinates)
295 :
296 : real(dp),allocatable :: kpt(:,:)
297 : ! gives the coordinates in the BZ of the kpoint
298 : ! size (3,nkpt)
299 :
300 : real(dp),allocatable :: wtk(:)
301 : !weight of each kpoint
302 :
303 : real(dp),allocatable :: acell(:)
304 : !size of the cell
305 :
306 : end type plowannier_type
307 : !!***
308 :
309 : CONTAINS !========================================================================================*
310 : !!***
311 :
312 :
313 : !!***f* m_plowannier/init_plowannier
314 : !! NAME
315 : !! init_plowannier
316 : !!
317 : !! FUNCTION
318 : !! initialize the variables useful for the computation
319 : !!
320 : !! INPUTS
321 : !! INPUTS
322 : !! plowan_bandf = max index of band for Wannier construction
323 : !! plowan_bandi = min index of band for Wannier construction
324 : !! plowan_compute = keyword to activate Wannier calculation
325 : !! plowan_iatom(plowan_natom) = index of atoms to use for Wannier
326 : !! plowan_it(plowan_nt)= index of atoms for real space calculation
327 : !! plowan_lcalc(sum_plowan_natom Plowan_nbl()) = index of l value for Wannier construction
328 : !! plowan_natom = nb of atoms for Wannier
329 : !! plowan_nbl(plowan_natom) = nb of l values for Wannier for each atoms.
330 : !! nl = nb of l values for Wannier for all atoms.
331 : !! plowan_nt = nb of atoms for real space calculation
332 : !! plowan_projcalc(sum_plowan_natom Plowan_nbl()) = index of projectors for Wannier construction
333 : !! acell_orig(3,nimage) = cell parameters
334 : !! kpt(3,nkpt) = k-points
335 : !! nkpt = nb of k-points
336 : !! nimage
337 : !! nspinor = nb of spinors
338 : !! nsppol = nb of polarization of wfc.
339 : !! wtk = weight of k-points
340 : !!
341 : !! OUTPUT
342 : !! wan : plowannier type
343 : !!
344 : !! SIDE EFFECTS
345 : !!
346 : !! NOTES
347 : !!
348 : !! SOURCE
349 :
350 :
351 0 : subroutine init_plowannier(plowan_bandf,plowan_bandi,plowan_compute,plowan_iatom,plowan_it,&
352 0 : &plowan_lcalc,plowan_natom,plowan_nbl,plowan_nt,plowan_projcalc,acell_orig,kpt,nl,nimage,nkpt,&
353 0 : &nspinor,nsppol,wtk,t2g,wan)
354 :
355 : !Arguments ----------------------------------
356 : !scalars
357 : ! type(dataset_type), intent(in) :: dtset
358 : integer,intent(in) ::plowan_bandi,plowan_bandf,plowan_natom,plowan_nt,plowan_compute
359 : integer,intent(in) ::nkpt,nsppol,nspinor,nimage,t2g,nl
360 : integer,intent(in) ::plowan_iatom(plowan_natom)
361 : integer,intent(in) ::plowan_nbl(plowan_natom)
362 : integer,intent(in) ::plowan_lcalc(nl)
363 : integer,intent(in) ::plowan_projcalc(nl)
364 : integer,intent(in) ::plowan_it(plowan_nt*3)
365 : real(dp),intent(in) :: kpt(3,nkpt)
366 : real(dp),intent(in) :: wtk(nkpt)
367 : real(dp),intent(in) :: acell_orig(3,nimage)
368 : type(plowannier_type), intent(inout) :: wan
369 :
370 : !Local --------------------------------------
371 : integer :: iatom,ikpt,ib,iband,il,iltot,it,ittot,ltemp,nn,norbtot
372 : character(len=500) :: message
373 : !************************************************************************
374 :
375 : !! generally
376 0 : wan%nkpt = nkpt
377 0 : wan%bandi_wan = plowan_bandi
378 0 : wan%bandf_wan = plowan_bandf
379 0 : wan%nsppol = nsppol
380 0 : wan%nspinor = nspinor
381 :
382 : !! for this case
383 0 : wan%natom_wan = plowan_natom
384 :
385 : !! generally
386 0 : ABI_MALLOC(wan%kpt,(3,size(kpt,2)))
387 0 : wan%kpt = kpt
388 0 : ABI_MALLOC(wan%iatom_wan,(wan%natom_wan))
389 0 : ABI_MALLOC(wan%nbl_atom_wan,(wan%natom_wan))
390 0 : wan%nbl_atom_wan = 0
391 0 : ABI_MALLOC(wan%latom_wan,(wan%natom_wan))
392 0 : ABI_MALLOC(wan%nbproj_atom_wan,(wan%natom_wan))
393 0 : wan%nbproj_atom_wan = 0
394 0 : ABI_MALLOC(wan%projector_wan,(wan%natom_wan))
395 0 : ABI_MALLOC(wan%position,(wan%natom_wan,3))
396 0 : wan%position = 0
397 0 : ABI_MALLOC(wan%wtk,(size(wtk,1)))
398 0 : wan%wtk(:) = wtk(:)
399 0 : ABI_MALLOC(wan%acell,(3))
400 0 : wan%acell(1) = acell_orig(1,1)
401 0 : wan%acell(2) = acell_orig(2,1)
402 0 : wan%acell(3) = acell_orig(3,1)
403 :
404 : ! If we want to study twice the same atom (but at different positions), use the same iatom and modify the positions below.
405 : ! In this case, the Wannier functions will be orthonormalized for one atom.
406 : ! For this particular reason, if we use the study twice the same atom at different position, each of them should have exactly the same projectors (it could be improved though by rewriting the normalization routine).
407 :
408 :
409 :
410 0 : ABI_MALLOC(wan%nposition,(wan%natom_wan))
411 : !write(std_out,*) "plowan_it", dtset%plowan_it
412 :
413 0 : iltot=0
414 0 : do iatom=1,wan%natom_wan
415 0 : wan%iatom_wan(iatom) = plowan_iatom(iatom)
416 0 : wan%nbl_atom_wan(iatom) = plowan_nbl (iatom)
417 0 : wan%nbproj_atom_wan(iatom) = plowan_nbl (iatom)
418 :
419 : ! Now we define for each atom the selected orbital moments.
420 0 : ABI_MALLOC(wan%latom_wan(iatom)%lcalc,(wan%nbl_atom_wan(iatom)))
421 0 : ABI_MALLOC(wan%projector_wan(iatom)%lproj,(wan%nbproj_atom_wan(iatom)))
422 0 : do il=1,wan%nbl_atom_wan(iatom)
423 0 : iltot=iltot+1
424 0 : wan%latom_wan(iatom)%lcalc(il)=plowan_lcalc(iltot)
425 0 : wan%projector_wan(iatom)%lproj(il)=plowan_projcalc(iltot)
426 0 : if (t2g==1 .and. plowan_lcalc(iltot)==2) then
427 0 : wan%latom_wan(iatom)%lcalc(il)=1
428 : endif
429 : enddo
430 :
431 : !For each iatom , pos is an array of two dimensions. The first one is
432 : !the number of lattice translation and the second one is ist
433 : !coordinates.
434 0 : ABI_MALLOC(wan%nposition(iatom)%pos,(plowan_nt,3))
435 0 : ittot=0
436 0 : do it=1,plowan_nt
437 0 : wan%nposition(iatom)%pos(it,1) = plowan_it(ittot+1)
438 0 : wan%nposition(iatom)%pos(it,2) = plowan_it(ittot+2)
439 0 : wan%nposition(iatom)%pos(it,3) = plowan_it(ittot+3)
440 : !write(std_out,*) "position",wan%nposition(iatom)%pos(it,:)
441 0 : ittot=ittot+3
442 : enddo
443 : enddo
444 :
445 : !!generally
446 0 : ABI_MALLOC(wan%psichi,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%natom_wan))
447 0 : do ikpt = 1,wan%nkpt
448 0 : do iband = wan%bandi_wan,wan%bandf_wan
449 0 : ib=iband-wan%bandi_wan+1
450 0 : do iatom = 1,wan%natom_wan
451 0 : ABI_MALLOC(wan%psichi(ikpt,ib,iatom)%atom,(wan%nbl_atom_wan(iatom)))
452 0 : do il = 1,wan%nbl_atom_wan(iatom)
453 0 : nn=(2*wan%latom_wan(iatom)%lcalc(il)+1)
454 0 : ABI_MALLOC(wan%psichi(ikpt,ib,iatom)%atom(il)%matl,(nn,wan%nsppol,wan%nspinor))
455 0 : wan%psichi(ikpt,ib,iatom)%atom(il)%matl = zero
456 : end do
457 : end do
458 : end do
459 : end do
460 0 : do iatom = 1,wan%natom_wan
461 0 : do il = 1,wan%nbl_atom_wan(iatom)
462 0 : ABI_MALLOC(wan%psichi(1,1,iatom)%atom(il)%ph0phiint,(10)) ! max number of proj for l =10..
463 : end do
464 : end do
465 :
466 :
467 : !sum of all the m possible
468 0 : wan%size_wan = 0
469 0 : do iatom = 1,wan%natom_wan
470 0 : do il = 1,wan%nbl_atom_wan(iatom)
471 0 : ltemp = wan%latom_wan(iatom)%lcalc(il)
472 0 : wan%size_wan = wan%size_wan + 2*ltemp + 1
473 : end do
474 : end do
475 :
476 0 : write(message,'(2a,i5,i5)') ch10,&
477 0 : & ' == Lower and upper values of the selected bands',wan%bandi_wan,wan%bandf_wan
478 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
479 0 : write(message,'(a,i10)') ' == Number of atoms ',wan%natom_wan
480 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
481 0 : write(message,'(a,9i2)') ' == Atoms selected ',(wan%iatom_wan(ltemp),ltemp=1,wan%natom_wan)
482 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
483 0 : write(message,'(a,9i2)') ' == Nb of angular momenta used for each atom ',(wan%nbl_atom_wan(ltemp),ltemp=1,wan%natom_wan)
484 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
485 0 : norbtot=0
486 0 : do iatom=1,wan%natom_wan
487 0 : write(message,'(a,i2,a,9i2)') ' == Value of the angular momenta for atom',iatom,' is : ',&
488 0 : & (wan%latom_wan(iatom)%lcalc(ltemp),ltemp=1,wan%nbl_atom_wan(iatom))
489 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
490 0 : do ltemp=1,wan%nbl_atom_wan(iatom)
491 0 : norbtot=norbtot+2*(wan%latom_wan(iatom)%lcalc(ltemp))+1
492 : enddo
493 0 : write(message,'(a,i2,a,9i2)') ' == Value of the projectors for atom',iatom,' is : ', &
494 0 : & (wan%projector_wan(iatom)%lproj(ltemp),ltemp=1,wan%nbl_atom_wan(iatom))
495 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
496 : enddo
497 0 : if(norbtot>wan%bandf_wan-wan%bandi_wan+1) then
498 0 : write(message,'(3a,2i6)') " Number of wannier functions is larger than" ,&
499 0 : &" number of Kohn Sham bands used for Wannier functions: decrease the number of Wannier functions", &
500 0 : &" or increase the number of bands ",norbtot,wan%bandf_wan-wan%bandi_wan+1
501 : !ABI_ERROR(message)
502 : endif
503 0 : if(plowan_compute==2) then
504 0 : write(message,'(3a)') ch10,' == plowan_compute=2 => off diag blocks in the k-space Wannier Hamiltonian matrix',&
505 0 : & 'is put to zero before diagonalisation'
506 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
507 : endif
508 :
509 0 : end subroutine init_plowannier
510 : !!***
511 :
512 :
513 : !!****f* m_plowannier/copy_orbital
514 : !! NAME
515 : !! copy_orbital
516 : !!
517 : !! FUNCTION
518 : !! Copy an array of orbital_type
519 : !!
520 : !! INPUTS
521 : !! lorbital1
522 : !!
523 : !! OUTPUT
524 : !! lorbital2
525 : !!
526 : !! SOURCE
527 :
528 :
529 :
530 0 : subroutine copy_orbital(orbital1,orbital2,n1,n2,n3)
531 :
532 : !Arguments----------------
533 : integer,intent(in) :: n1,n2,n3
534 : type(orbital_type), intent(in) :: orbital1(n1,n2,n3)
535 : type(orbital_type),intent(inout) :: orbital2(n1,n2,n3)
536 :
537 : !Local variable-----------
538 : integer :: n4,n5,n6,n7
539 : integer :: i,j,k,l,m,p,q
540 :
541 0 : do i = 1,n1
542 0 : do j = 1,n2
543 0 : do k = 1,n3
544 0 : n4 = size(orbital1(i,j,k)%atom,1)
545 0 : do l = 1,n4
546 0 : n5 = size(orbital1(i,j,k)%atom(l)%matl,1)
547 0 : n6 = size(orbital1(i,j,k)%atom(l)%matl,2)
548 0 : n7 = size(orbital1(i,j,k)%atom(l)%matl,3)
549 0 : do m = 1,n5
550 0 : do p = 1,n6
551 0 : do q = 1,n7
552 0 : orbital2(i,j,k)%atom(l)%matl(m,p,q) = orbital1(i,j,k)%atom(l)%matl(m,p,q)
553 : end do
554 : end do
555 : end do
556 : end do
557 : end do
558 : end do
559 : end do
560 :
561 0 : end subroutine copy_orbital
562 : !!***
563 :
564 : !!****f* m_plowannier/allocate_orbital
565 : !! NAME
566 : !! allocate_orbital
567 : !!
568 : !! FUNCTION
569 : !! allocate an array of orbital_type
570 : !!
571 : !! INPUTS
572 : !! lorbital1
573 : !!
574 : !! OUTPUT
575 : !! lorbital2
576 : !!
577 : !! SOURCE
578 :
579 :
580 :
581 0 : subroutine allocate_orbital(orbital1,orbital2,n1,n2,n3)
582 :
583 : !Arguments----------------
584 : integer,intent(in) :: n1,n2,n3
585 : type(orbital_type), intent(in) :: orbital1(n1,n2,n3)
586 : type(orbital_type),intent(inout) :: orbital2(n1,n2,n3)
587 :
588 : !Local variable-----------
589 : integer :: n4,n5,n6,n7
590 : integer :: i,j,k,l
591 :
592 0 : do i = 1,n1
593 0 : do j = 1,n2
594 0 : do k = 1,n3
595 0 : n4 = size(orbital1(i,j,k)%atom,1)
596 0 : ABI_MALLOC(orbital2(i,j,k)%atom,(n4))
597 0 : do l = 1,n4
598 0 : n5 = size(orbital1(i,j,k)%atom(l)%matl,1)
599 0 : n6 = size(orbital1(i,j,k)%atom(l)%matl,2)
600 0 : n7 = size(orbital1(i,j,k)%atom(l)%matl,3)
601 0 : ABI_MALLOC(orbital2(i,j,k)%atom(l)%matl,(n5,n6,n7))
602 : end do
603 : end do
604 : end do
605 : end do
606 :
607 0 : end subroutine allocate_orbital
608 : !!***
609 :
610 :
611 : !!****f* m_plowannier/destroy_orbital
612 : !! NAME
613 : !! destroy_orbital
614 : !!
615 : !! FUNCTION
616 : !! destroy an array of orbital_type
617 : !!
618 : !! INPUTS
619 : !! lorbital1
620 : !!
621 : !! OUTPUT
622 : !! lorbital2
623 : !!
624 : !! SOURCE
625 :
626 :
627 :
628 0 : subroutine destroy_orbital(orbital2,n1,n2,n3)
629 :
630 : !Arguments----------------
631 : integer,intent(in) :: n1,n2,n3
632 : type(orbital_type),intent(inout) :: orbital2(n1,n2,n3)
633 :
634 : !Local variable-----------
635 : integer :: n4
636 : integer :: i,j,k,l
637 :
638 0 : do i = 1,n1
639 0 : do j = 1,n2
640 0 : do k = 1,n3
641 0 : n4 = size(orbital2(i,j,k)%atom,1)
642 0 : do l = 1,n4
643 0 : ABI_FREE(orbital2(i,j,k)%atom(l)%matl)
644 : end do
645 0 : ABI_FREE(orbital2(i,j,k)%atom)
646 : end do
647 : end do
648 : end do
649 :
650 0 : end subroutine destroy_orbital
651 : !!***
652 :
653 :
654 : !!***f* m_plowannier/compute_coeff_plowannier
655 : !! NAME
656 : !! compute_coeff_plowannier
657 : !!
658 : !! FUNCTION
659 : !! Compute the coefficient
660 : !!
661 : !! INPUTS
662 : !! cryst_struc <type(crystal_t)>=crystal structure data
663 : !! -gprimd(3,3)=dimensional reciprocal space primitive translations
664 : !! -indsym(4,nsym,natom)=indirect indexing array for atom labels
665 : !! -symrec(3,3,nsym)=symmetry operations in reciprocal space
666 : !! -nsym= number of symetry operations
667 : !! cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
668 : !! and each |p_lmn> non-local projector
669 : !! dimcprj(natom) = dimension for cprj
670 : !! dtset <type(dataset_type)>=all input variables for this dataset
671 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
672 : !! fermie= Fermi energy
673 : !! mband=maximum number of bands
674 : !! mbandcprj=
675 : !! mkmem =number of k points treated by this node
676 : !! mpi_enreg=information about MPI parallelization
677 : !! nkpt=number of k points.
678 : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
679 : !! nsppol=1 for unpolarized, 2 for spin-polarized
680 : !! occ(mband*nkpt*nsppol) = occupancies of KS states.
681 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
682 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
683 : !! usecprj=
684 : !! unpaw=file number for cprj
685 : !! nbandkss
686 : !! dtfil
687 : !!
688 : !! OUTPUT
689 : !! wan%psichi: projections <Psi|chi>
690 : !!
691 : !! SIDE EFFECTS
692 : !! (only writing, printing)
693 : !!
694 : !! NOTES
695 : !!
696 : !! SOURCE
697 :
698 :
699 0 : subroutine compute_coeff_plowannier(cryst_struc,cprj,dimcprj,dtset,eigen,fermie,&
700 0 : & mpi_enreg,occ,wan,pawtab,psps,usecprj,unpaw,pawrad,dtfil)
701 :
702 :
703 : use m_hide_lapack
704 :
705 : !Arguments ------------------------------------
706 : !scalars
707 :
708 : type(plowannier_type),intent(inout) :: wan
709 : integer,intent(in) :: unpaw,usecprj
710 : real(dp),intent(in) :: fermie
711 : type(MPI_type),intent(in) :: mpi_enreg
712 : type(dataset_type),intent(in) :: dtset
713 : type(pseudopotential_type),intent(in) :: psps
714 : type(crystal_t),intent(in) :: cryst_struc
715 : !arrays
716 : integer, intent(in) :: dimcprj(cryst_struc%natom)
717 : real(dp),intent(in) :: eigen(dtset%mband*wan%nkpt*wan%nsppol)
718 : real(dp),intent(in) :: occ(dtset%mband*wan%nkpt*wan%nsppol)
719 : type(pawcprj_type), intent(in) :: cprj(cryst_struc%natom,wan%nspinor*dtset%mband*dtset%mkmem*wan%nsppol*usecprj)
720 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
721 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
722 : type(datafiles_type),intent(in) :: dtfil
723 :
724 : !Local variables-------------------------------
725 : !scalars
726 : integer :: band_index,dimpsichi,facpara
727 : integer :: iatom,iatom1,iatom2,iband,ibandc,ibg,ierr,ikpt
728 : integer :: iband1,owrunt,opt
729 : integer :: ilmn,iorder_cprj,ispinor,isppol,itypat,ilmn2
730 : integer :: lmn_size, t2g, m1_t2g, m1_t2g_mod
731 : integer :: m1,maxnproju,me,natom,nband_k,nband_k_cprj
732 : integer :: nnn,nprocband,spaceComm
733 : integer :: plowan_greendos,plowan_hybrid,plowan_inter,plowan_computegreen
734 : real(dp) :: ph0phiint_used
735 : character(len=500) :: message
736 : character(len=50) :: mat_writing2,mat_writing2_out
737 : character(len=5000) :: mat_writing,mat_writing_out
738 : integer :: l1,count,mesh_size,il,count_total,l,proj
739 : integer :: il1,il2,im1,im2,index_c,index_l,ispinor1,ispinor2,sizem,pos1,pos2
740 : real(dp) :: int_current,sum,sum2,sum3
741 :
742 : complex(dp) :: wbase,wcurrent
743 : real(dp) :: resolution, wincrease,wmax,wmin
744 : integer :: iw,dos,shift,unt,unt2,dos_unt,dos_unt2
745 : integer :: number_of_frequencies,band_struct,prtocc,prtint
746 : real(dp) :: convert
747 : complex(dp):: xsum
748 : character(len=fnlen) :: owrfile
749 :
750 : !arrays
751 : real(dp) :: chinorm
752 0 : complex(dp), allocatable :: Fff(:)
753 0 : complex(dp), allocatable :: buffer1(:)
754 : logical :: lprojchi
755 0 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
756 0 : type(operwan_type), allocatable :: operwan(:,:,:)
757 0 : type(operwan_realspace_type) :: operwan_realspace
758 0 : type(operwan_realspace_type) :: operocc
759 0 : complex(dp), allocatable :: eigenks(:,:,:,:)
760 0 : complex(dp), allocatable :: operks(:,:,:,:)
761 0 : complex(dp), allocatable :: identityks(:,:,:,:)
762 0 : real(dp), allocatable :: ff(:)
763 0 : complex(dp), allocatable :: operwansquare(:,:,:,:)
764 0 : complex(dp), allocatable :: operwansquarereal(:,:,:)
765 0 : complex(dp), allocatable :: matrix_to_diag(:,:)
766 0 : complex(dp), allocatable :: energies(:,:)
767 0 : complex(dp), allocatable :: Ffftable(:,:)
768 : character(len = 5) :: i2s,x1
769 :
770 : !To diagonalize eigenvalues
771 0 : real(dp), allocatable :: eig(:), rwork(:)
772 0 : complex(dp), allocatable :: zwork(:)
773 : integer :: lwork,info,whole_diag
774 : !complex(dp), allocatable :: densmat(:,:)
775 : !************************************************************************
776 :
777 : ! Drive the normalization of the psichis
778 :
779 0 : if (dtset%nbandkss==-1 .and.dtset%kssform==3 )then
780 0 : if (dtset%ucrpa >= 1 .or. dtset%dmft_kspectralfunc==1) then
781 0 : opt = 0
782 : else
783 0 : opt=1
784 : endif
785 : else
786 0 : opt=0
787 : end if
788 :
789 0 : if (opt==0) then
790 0 : write(message,*)ch10,"Normalization of plowannier k-point by k-point"
791 : else
792 0 : write(message,*)ch10,"Normalization of plowannier on the sum of the k-points"
793 : endif
794 :
795 0 : ABI_COMMENT(message)
796 :
797 0 : t2g=dtset%dmft_t2g
798 0 : if (t2g==1) then
799 0 : write(message,*)ch10,"Using only t2g bands in plowannier, this is under test and may lead to bugs"
800 : end if
801 :
802 0 : ABI_COMMENT(message)
803 : !opt=1
804 : ! 0 : normalization k-point by k-point (normal use of plowan)
805 : ! 1 : normalization of the sum over k-points (use with crpa old keywords)
806 :
807 :
808 : ! Internal variables (could be put one day as input variables of ABINIT).
809 0 : plowan_computegreen = 0 !
810 : ! 0 : do nothing do not compute hybri or dos
811 : ! 1 : Compute hybridization or dos (depends on following keywords)
812 : ! 2 : not tested, probably with bug included.: compute hybridization in Wannier basis
813 :
814 : ! If plowan_computegreen>0, the following data is useful
815 0 : plowan_greendos = 1 ! For the first atom, plowan_greendos is the index
816 : ! of the angular momentum in array
817 : ! wan%latom_wan(iatom)%lcalc: it is thus betwween 0 and wan%nbl_atom_wan(iatom)
818 : ! It is not the value of the angular momentum but its
819 : ! index
820 0 : plowan_hybrid = 0 ! Same convention as for greendos
821 0 : plowan_inter = 1 ! compute all interaction between all atoms all orbitals all neighbours, requires plowan_realspace=1
822 :
823 : !owrfile = trim(dtfil%filnam_ds(4))//"_operwan_realspace"
824 0 : owrfile = "__operwan_realspace__"
825 :
826 : dos = 0
827 : if(plowan_greendos>0) dos=plowan_greendos
828 : if(plowan_hybrid>0) dos=-plowan_hybrid
829 : if(plowan_hybrid>0.and.plowan_greendos>0) then
830 : write(message,*) " plowan_hybrid and plowan_greendos cannot be both >0"
831 : ABI_ERROR(message)
832 : endif
833 :
834 : ! GREEN STUDY PARAMETERS (FREQUENCIES)
835 : ! ===================================
836 : !!To choose the frequencies in the Green study
837 0 : wmin=-2.d0 ! eV
838 0 : wmax= 2.d0 ! eV
839 0 : resolution= 0.02 ! eV
840 0 : wbase = cmplx(wmin/27.2107,0.001,kind=dp) ! most negative value of frequency
841 0 : wincrease = resolution/27.2107 ! step
842 0 : number_of_frequencies = int((wmax-wmin)/resolution)
843 :
844 :
845 :
846 : ! Select the way of diagonalisation
847 : !===================================
848 0 : whole_diag = 1! 1 for diagonalization of the whole matrix, 0 for each orbital
849 0 : if(dtset%plowan_compute==2) whole_diag = 0 ! off diagonal blocks are suppressed in the hamiltonian matrix before diagonalisation
850 0 : band_struct = 1 ! 1 for plotting band struct (Wannier bands)
851 :
852 : ! Select the real space calculation of Wannier function: Interpolation
853 : ! versus Analysis
854 : !===================================
855 0 : prtocc = 0 ! occupations have no meaning for a k-point path so the default is 0
856 : if(dtset%kptopt>0.and.dtset%plowan_realspace>=1) prtocc = 1 !1 to print the occupation in real space
857 :
858 :
859 : ! Select if computation of interactions is done
860 : !===================================
861 0 : prtint = plowan_inter !1 to print sqrt(sum of interaction squared) between orbitals, we do not use the input file to do this one
862 :
863 :
864 : ! Plot KS band structure
865 : !===================================
866 : ! !data for printing KS bands
867 : ! do ikpt = 1,wan%nkpt
868 : ! print* ,'bandstruct', real(eigen(1+(ikpt-1)*30:30+(ikpt-1)*30))!, real(eigen(7471+(ikpt-1)*30:7500+(ikpt-1)*30))
869 : ! end do
870 :
871 :
872 :
873 : !DBG_ENTER("COLL")
874 : !Fake test to keep fermie as argument. REMOVE IT AS SOON AS POSSIBLE ...
875 0 : if(fermie>huge(zero))chinorm=zero
876 :
877 0 : facpara=1 !mpi_enreg%nproc
878 0 : if(abs(dtset%pawprtvol)>=3) then
879 0 : write(message,*) ch10, " number of k-points used is nkpt = ", dtset%nkpt
880 0 : call wrtout(std_out, message,'COLL')
881 0 : write(message,*) " warning: parallelised version ", dtset%nkpt
882 0 : call wrtout(std_out, message,'COLL')
883 0 : write(message,*) " weights k-points used is wtk = wtk"
884 0 : call wrtout(std_out, message,'COLL')
885 : end if
886 :
887 0 : if(usecprj==0) then
888 0 : write(message,*) " usecprj=0 : BUG in init_plowannier",usecprj
889 0 : ABI_BUG(message)
890 : end if
891 :
892 0 : if(wan%nspinor/=dtset%nspinor) then
893 0 : write(message,*) " wan%nspinor=/dtset%nspinor, init_plowannier is not working in this case",&
894 0 : & wan%nspinor,dtset%nspinor
895 0 : ABI_ERROR(message)
896 : end if
897 :
898 :
899 :
900 : !----------------------------------- MPI-------------------------------------
901 :
902 : !Init parallelism
903 0 : spaceComm=mpi_enreg%comm_cell
904 0 : if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
905 0 : me=mpi_enreg%me_kpt
906 :
907 : !----------------------------------- MPI-------------------------------------
908 :
909 :
910 : lprojchi=.false.
911 0 : lprojchi=.true.
912 0 : natom=cryst_struc%natom
913 :
914 :
915 0 : write(message,'(2a)') ch10,&
916 0 : & ' == Prepare data for projected local orbital wannier function calculation '
917 0 : call wrtout(std_out,message,'COLL')
918 0 : if(abs(dtset%pawprtvol)>=3) then
919 0 : write(message, '(a,a)' ) ch10,&
920 0 : & '---------------------------------------------------------------'
921 : ! call wrtout(ab_out,message,'COLL');call wrtout(std_out, message,'COLL')
922 0 : call wrtout(std_out, message,'COLL')
923 0 : write(message, '(a,a,a,a,a,a,a,a,a,a,a,a)' ) ch10,&
924 0 : & ' Print useful data (as a check)',ch10,&
925 0 : & ' - Overlap of KS wfc with atomic orbital inside sphere',ch10,&
926 0 : & ' - Eigenvalues',ch10,&
927 0 : & ' - Weights of k-points',ch10,&
928 0 : & ' - Number of spins ',ch10,&
929 0 : & ' - Number of states'
930 : ! call wrtout(ab_out,message,'COLL');call wrtout(std_out, message,'COLL')
931 0 : call wrtout(std_out, message,'COLL')
932 0 : write(message, '(a,a)' ) ch10,&
933 0 : & '---------------------------------------------------------------'
934 : end if
935 0 : if(dtset%nstep==0) then
936 0 : message = 'nstep should be greater than 1'
937 0 : ABI_BUG(message)
938 : end if
939 :
940 :
941 : !********************* Max Values for U terms.
942 : !maxlpawu=0
943 0 : maxnproju=0
944 0 : do iatom=1,natom
945 : if(pawtab(dtset%typat(iatom))%lpawu.ne.-1 .and. pawtab(dtset%typat(iatom))%nproju.gt.maxnproju)&
946 : & maxnproju=pawtab(dtset%typat(iatom))%nproju
947 : end do
948 : !***************** in forlb.eig
949 0 : if(me.eq.0.and.abs(dtset%pawprtvol)>=3) then
950 0 : if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) then
951 0 : ABI_ERROR(message)
952 : end if
953 0 : rewind(unt)
954 0 : write(unt,*) " Number of bands, spins, and k-point; and spin-orbit flag"
955 0 : write(unt,*) dtset%mband,wan%nsppol,wan%nkpt,wan%nspinor,wan%bandi_wan,wan%bandf_wan
956 0 : write(unt,*) " For each k-point, eigenvalues for each band"
957 0 : write(unt,*) (dtset%wtk(ikpt)*facpara,ikpt=1,wan%nkpt)
958 0 : band_index=0
959 0 : do isppol=1,wan%nsppol
960 0 : write(unt,*) " For spin"
961 0 : write(unt,*) isppol
962 0 : do ikpt=1,wan%nkpt
963 0 : nband_k=dtset%nband(ikpt+(isppol-1)*wan%nkpt)
964 0 : write(unt,*) " For k-point"
965 0 : write(unt,*) ikpt
966 0 : do iband=wan%bandi_wan,wan%bandf_wan
967 0 : write(unt, '(2i6,4x,f20.15)' ) iband-wan%bandi_wan+1,ikpt,eigen(iband+band_index)*2.d0
968 : end do
969 0 : band_index=band_index+nband_k
970 : end do
971 : end do
972 0 : close(unt)
973 : end if ! proc=me
974 :
975 :
976 :
977 : !-----------------------------------------------------------------
978 : ! Allocates, computes and stocks the ph0phiint
979 : !-----------------------------------------------------------------
980 : !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
981 0 : do iatom = 1,wan%natom_wan
982 0 : do il = 1,wan%nbl_atom_wan(iatom)
983 0 : count = 0
984 0 : count_total = 0
985 0 : l = wan%latom_wan(iatom)%lcalc(il)
986 0 : if (t2g==1) then
987 0 : l=2
988 : endif
989 0 : proj = wan%projector_wan(iatom)%lproj(il)
990 0 : itypat = dtset%typat(wan%iatom_wan(iatom))
991 0 : lmn_size = pawtab(itypat)%lmn_size
992 : ! modif proj if proj == -2 ---> usemdft has been used
993 0 : if (proj==-2) then
994 0 : do ilmn = 1,lmn_size
995 0 : if ( psps%indlmn(1,ilmn,itypat) .eq. l .and. psps%indlmn(2,ilmn,itypat) .eq. 0 .and. proj .eq. -2 ) then
996 0 : proj=psps%indlmn(5,ilmn,itypat)
997 : end if
998 : end do
999 : end if
1000 : ! check if the choice of proj is coherent with the value of l
1001 0 : do ilmn = 1,lmn_size
1002 0 : if (psps%indlmn(1,ilmn,itypat).eq. l .and. psps%indlmn(2,ilmn,itypat) .eq. 0) then
1003 0 : count_total = count_total+1 !!counts the number total of projector for this l
1004 0 : if (psps%indlmn(5,ilmn,itypat) .eq. proj) then
1005 0 : count = count+1 !!the projector chosen is in the right l
1006 : end if
1007 : end if
1008 : end do
1009 0 : if (count .eq. 0) then
1010 0 : write(message,'(a)') " The projector choice is not consistent with the orbital l"
1011 0 : ABI_ERROR(message)
1012 : else !good choice of projector
1013 0 : wan%psichi(1,1,iatom)%atom(il)%ph0phiint = zero
1014 0 : do ilmn2 = 1,lmn_size
1015 0 : if (psps%indlmn(1,ilmn2,itypat) .eq. l .and. psps%indlmn(2,ilmn2,itypat) .eq. 0) then
1016 0 : mesh_size = pawtab(itypat)%mesh_size
1017 0 : ABI_MALLOC(ff,(mesh_size))
1018 0 : ff(1:mesh_size) = pawtab(itypat)%phi(1:mesh_size,proj)*pawtab(itypat)%phi(1:mesh_size,psps%indlmn(5,ilmn2,itypat))
1019 : ! ff(1:mesh_size) = pawtab(itypat)%tphi(1:mesh_size,proj)*pawtab(itypat)%tphi(1:mesh_size,psps%indlmn(5,ilmn2,itypat))
1020 : int_current = 0
1021 0 : call simp_gen(int_current,ff,pawrad(itypat)) !call the subroutine which does the computation
1022 0 : wan%psichi(1,1,iatom)%atom(il)%ph0phiint(psps%indlmn(3,ilmn2,itypat)) = int_current !we put the values for ikpt = 1 and iband = 1
1023 0 : ABI_FREE(ff)
1024 : end if
1025 : end do
1026 : end if
1027 : end do
1028 : end do
1029 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1030 :
1031 :
1032 : !==========================================================================
1033 : !***************** Compute <Psi|Chi>=\sum_{proja} <Psi|P_a><phi_a|Chi>
1034 : !==========================================================================
1035 :
1036 : !Allocate temporary cwaveprj storage
1037 0 : ABI_MALLOC(cwaveprj,(natom,wan%nspinor))
1038 :
1039 0 : call pawcprj_alloc(cwaveprj,0,dimcprj)
1040 :
1041 0 : nprocband=(dtset%mband/dtset%mband)
1042 0 : ibg=0
1043 0 : do isppol=1,wan%nsppol
1044 0 : do ikpt=1,wan%nkpt
1045 0 : nband_k=dtset%nband(ikpt+(isppol-1)*wan%nkpt)
1046 0 : nband_k_cprj=nband_k/nprocband
1047 : !nband_k is mband for each k, so it is mband most of the time
1048 0 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
1049 0 : do iband=wan%bandi_wan,wan%bandf_wan !loop only over the bands we are interested in
1050 0 : ibandc=iband-wan%bandi_wan+1
1051 : ! Parallelization: treat only some bands
1052 0 : if (dtset%paral_kgb==1) then
1053 0 : if (mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
1054 : else
1055 0 : if (mpi_enreg%proc_distrb(ikpt,ibandc,isppol)/=me) cycle
1056 : end if
1057 0 : do ispinor=1,wan%nspinor
1058 0 : do iatom = 1,wan%natom_wan !loop over the atom chosen
1059 0 : itypat = dtset%typat(wan%iatom_wan(iatom))
1060 0 : lmn_size = pawtab(itypat)%lmn_size !retrieve the number of different lmn
1061 : call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,iband,ibg,ikpt,&
1062 : & iorder_cprj,isppol,dtset%mband,dtset%mkmem,dtset%natom,1,nband_k_cprj,&
1063 : & wan%nspinor,wan%nsppol,unpaw,mpicomm=mpi_enreg%comm_kpt,&
1064 0 : & proc_distrb=mpi_enreg%proc_distrb)
1065 0 : chinorm=1.d0
1066 0 : do l1 = 1,wan%nbl_atom_wan(iatom) !l1 is the index of the orbital
1067 0 : l = wan%latom_wan(iatom)%lcalc(l1) !l is the value of the orbital (linked to index l1)
1068 0 : if (t2g==1) then
1069 0 : l=2
1070 0 : m1_t2g=0
1071 : endif
1072 0 : do ilmn = 1,lmn_size
1073 0 : if (psps%indlmn(1,ilmn,itypat) .eq. l) then!
1074 0 : ph0phiint_used = wan%psichi(1,1,iatom)%atom(l1)%ph0phiint(psps%indlmn(3,ilmn,itypat))
1075 0 : m1 = psps%indlmn(2,ilmn,itypat)+l+1
1076 0 : if (t2g==1) then
1077 0 : if(m1==1.or.m1==2.or.m1==4) then
1078 0 : m1_t2g=m1_t2g+1
1079 0 : m1_t2g_mod=mod(m1_t2g-1,3)+1
1080 : wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1_t2g_mod,isppol,ispinor)=&
1081 : wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1_t2g_mod,isppol,ispinor)+&
1082 : cmplx(cwaveprj(wan%iatom_wan(iatom),ispinor)%cp(1,ilmn)*ph0phiint_used,cwaveprj(&
1083 0 : wan%iatom_wan(iatom),ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp)
1084 : endif
1085 : else
1086 : wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1,isppol,ispinor)=&
1087 : wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1,isppol,ispinor)+&
1088 : cmplx(cwaveprj(wan%iatom_wan(iatom),ispinor)%cp(1,ilmn)*ph0phiint_used,cwaveprj(&
1089 0 : wan%iatom_wan(iatom),ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp)
1090 : end if
1091 : end if
1092 : end do
1093 : end do
1094 : end do ! iatom
1095 : end do ! ispinor
1096 : end do !iband
1097 0 : ibg=ibg+nband_k_cprj*wan%nspinor !useful to select the right ikpt in pawcprj_get
1098 : end do !ikpt
1099 : end do ! isppol
1100 :
1101 :
1102 : !===========================================================
1103 : !************************ new gather info for MPI
1104 : !===========================================================
1105 :
1106 0 : dimpsichi=0
1107 0 : do iatom = 1,wan%natom_wan
1108 0 : do l1 = 1,wan%nbl_atom_wan(iatom)
1109 0 : dimpsichi = dimpsichi + wan%nkpt*(wan%bandf_wan-wan%bandi_wan+1)*(2*wan%latom_wan(iatom)%lcalc(l1)+1)*wan%nsppol*wan%nspinor
1110 : end do
1111 : end do
1112 0 : dimpsichi = 2*dimpsichi !for complex
1113 0 : ABI_MALLOC(buffer1,(dimpsichi))
1114 0 : buffer1 = zero
1115 0 : nnn = 0
1116 0 : do ikpt = 1,wan%nkpt
1117 0 : do ibandc = 1,wan%bandf_wan-wan%bandi_wan+1
1118 0 : do iatom=1,wan%natom_wan
1119 0 : do l1 = 1,wan%nbl_atom_wan(iatom)
1120 0 : do m1 = 1,2*wan%latom_wan(iatom)%lcalc(l1)+1
1121 0 : do isppol = 1,wan%nsppol
1122 0 : do ispinor = 1,wan%nspinor
1123 0 : nnn=nnn+1
1124 0 : buffer1(nnn)=wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1,isppol,ispinor)
1125 : end do
1126 : end do
1127 : end do
1128 : end do
1129 : end do
1130 : end do
1131 : end do
1132 0 : call xmpi_barrier(spaceComm)
1133 0 : call xmpi_sum(buffer1,spaceComm,ierr)
1134 : if (dtset%paral_kgb==1 .and. nprocband > 1) then
1135 : call xmpi_sum(buffer1,mpi_enreg%comm_band,ierr) !build sum over band processors
1136 : end if
1137 0 : call xmpi_barrier(spaceComm)
1138 0 : nnn = 0
1139 0 : do ikpt = 1,wan%nkpt
1140 0 : do ibandc = 1,wan%bandf_wan-wan%bandi_wan+1
1141 0 : do iatom = 1,wan%natom_wan
1142 0 : do l1 = 1,wan%nbl_atom_wan(iatom)
1143 0 : do m1 = 1,2*wan%latom_wan(iatom)%lcalc(l1)+1
1144 0 : do isppol = 1,wan%nsppol
1145 0 : do ispinor = 1,wan%nspinor
1146 0 : nnn=nnn+1
1147 0 : wan%psichi(ikpt,ibandc,iatom)%atom(l1)%matl(m1,isppol,ispinor)=buffer1(nnn)
1148 : end do
1149 : end do
1150 : end do
1151 : end do
1152 : end do
1153 : end do
1154 : end do
1155 0 : ABI_FREE(buffer1)
1156 :
1157 0 : call xmpi_barrier(spaceComm)
1158 :
1159 :
1160 :
1161 : !! -------------------------------------------------------------
1162 : !! COMPUTATION OF THE OCCUPATION MATRIX BEFORE NORMALIZATION
1163 : !! -------------------------------------------------------------
1164 : !! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1165 :
1166 0 : if (dtset%prtvol >= 5) then
1167 : !Inialize an empty Wannier operator
1168 0 : ABI_MALLOC(operwan,(wan%nkpt,wan%natom_wan,wan%natom_wan))
1169 0 : call initialize_operwan(wan,operwan)
1170 :
1171 : !Creation of the KS occupation operator
1172 0 : ABI_MALLOC(eigenks,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%bandf_wan-wan%bandi_wan+1,wan%nsppol))
1173 0 : ABI_MALLOC(identityks,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%bandf_wan-wan%bandi_wan+1,wan%nsppol))
1174 0 : eigenks = czero
1175 0 : identityks=czero
1176 0 : do isppol = 1,wan%nsppol
1177 0 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
1178 0 : ibandc = iband1 + wan%bandi_wan - 1
1179 0 : do ikpt = 1,wan%nkpt
1180 0 : eigenks(ikpt,iband1,iband1,isppol) = occ(((ikpt-1)*dtset%mband+ibandc+(isppol-1)*wan%nkpt*dtset%mband))
1181 : !write(6,*) 'eigenks', ikpt,iband1,isppol,((ikpt-1)*dtset%mband+ibandc+(isppol-1)*wan%nkpt*dtset%mband),occ(((ikpt-1)*dtset%mband+ibandc+(isppol-1)*wan%nkpt*dtset%mband))
1182 : end do
1183 : end do
1184 : end do
1185 :
1186 :
1187 : !compute the occupation in wannier basis and print it
1188 0 : write(message,*)char(10),&
1189 0 : &" Print the occupation levels (not normalized) for 1 atom, 1 orbitals"
1190 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1191 0 : write(message,*)" Atom =",wan%iatom_wan(1),"orbital =",wan%latom_wan(1)%lcalc(1)
1192 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1193 0 : do ikpt = 1,wan%nkpt
1194 0 : call compute_oper_ks2wan(wan,eigenks,operwan,ikpt)
1195 : enddo
1196 0 : call init_operwan_realspace(wan,operocc)
1197 0 : write(message,*)char(10)," The occupation matrix before normalization is :"
1198 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1199 0 : call compute_oper_wank2realspace(wan,operwan,operocc)
1200 0 : if (wan%nsppol ==1)then
1201 0 : write(message,*)char(10)," For one spin :"
1202 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1203 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1204 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,1,1,1))
1205 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1206 : enddo
1207 : else
1208 0 : write(message,*)char(10)," For spin up :"
1209 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1210 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1211 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,1,1,1))
1212 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1213 : enddo
1214 0 : write(message,*)char(10)," For spin down : "
1215 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1216 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1217 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,2,1,1))
1218 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1219 : enddo
1220 : endif
1221 :
1222 : endif
1223 :
1224 : !!-------------------------------------------------------------
1225 : !!NORMALIZATION
1226 : !!--------------------------------------------------------------
1227 :
1228 0 : if(dtset%prtvol>=5) then
1229 0 : do isppol = 1,wan%nsppol
1230 0 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
1231 0 : ibandc = iband1 + wan%bandi_wan - 1
1232 0 : do ikpt = 1,wan%nkpt
1233 0 : identityks(ikpt,iband1,iband1,isppol) = one
1234 : !write(6,*) "idks", identityks(ikpt,iband1,iband1,isppol)
1235 : end do
1236 : end do
1237 : end do
1238 0 : call zero_operwan(wan,operwan)
1239 0 : do ikpt = 1,wan%nkpt
1240 0 : call compute_oper_ks2wan(wan,identityks,operwan,ikpt)
1241 : enddo
1242 0 : do ikpt = 1,wan%nkpt
1243 : ! if (ikpt<=5)then
1244 0 : write(message,*)char(10)," For ikpt=",ikpt,"the normalization matrix is before normalization :"
1245 0 : call wrtout(std_out,message,'COLL'); !call wrtout(ab_out,message,'COLL')
1246 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1247 0 : write(mat_writing,'(7f20.5)')real(operwan(ikpt,1,1)%atom(1,1)%matl(im1,:,1,1,1))
1248 0 : call wrtout(std_out,mat_writing,'COLL'); !call wrtout(ab_out,mat_writing,'COLL')
1249 : enddo
1250 : !endif
1251 : end do
1252 : endif
1253 :
1254 :
1255 0 : call normalization_plowannier(wan,opt)
1256 :
1257 :
1258 0 : if (dtset%prtvol>=5) then
1259 0 : call zero_operwan(wan,operwan)
1260 0 : do ikpt = 1,wan%nkpt
1261 0 : call compute_oper_ks2wan(wan,identityks,operwan,ikpt)
1262 : enddo
1263 0 : do ikpt = 1,wan%nkpt
1264 : !if (ikpt<=5)then
1265 0 : write(message,*)char(10)," For ikpt=",ikpt,"the normalization matrix is after normalization :"
1266 0 : call wrtout(std_out,message,'COLL'); !call wrtout(ab_out,message,'COLL')
1267 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1268 0 : write(mat_writing,'(7f20.5)')real(operwan(ikpt,1,1)%atom(1,1)%matl(im1,:,1,1,1))
1269 0 : call wrtout(std_out,mat_writing,'COLL'); !call wrtout(ab_out,mat_writing,'COLL')
1270 : enddo
1271 : !endif
1272 : end do
1273 : endif
1274 : !! -------------------------------------------------------------
1275 : !! COMPUTATION OF THE OCCUPATION MATRIX AFTER NORMALIZATION
1276 : !! -------------------------------------------------------------
1277 : !! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1278 : !compute the occupation in wannier basis and print it
1279 :
1280 0 : if (dtset%prtvol >= 5) then
1281 0 : write(message,*)char(10),&
1282 0 : &" Print the occupation levels (normalized) for 1 atom, 1 orbital"
1283 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1284 0 : write(message,*)"Atom =",wan%iatom_wan(1),"orbital =",wan%latom_wan(1)%lcalc(1)
1285 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1286 0 : call zero_operwan(wan,operwan)
1287 0 : do ikpt = 1,wan%nkpt
1288 0 : call compute_oper_ks2wan(wan,eigenks,operwan,ikpt)
1289 : enddo
1290 0 : call zero_operwan_realspace(wan,operocc)
1291 0 : write(message,*)char(10)," The occupation matrix after normalization is :"
1292 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1293 0 : call compute_oper_wank2realspace(wan,operwan,operocc)
1294 0 : if (wan%nsppol ==1)then
1295 0 : write(message,*)char(10)," For one spin :"
1296 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1297 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1298 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,1,1,1))
1299 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1300 : enddo
1301 : else
1302 0 : write(message,*)char(10)," For spin up :"
1303 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1304 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1305 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,1,1,1))
1306 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1307 : enddo
1308 0 : write(message,*)char(10)," For spin down : "
1309 0 : call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL')
1310 0 : do im1=1,2*wan%latom_wan(1)%lcalc(1)+1
1311 0 : write(mat_writing,'(7f20.10)')real(operocc%atom_index(1,1)%position(1,1)%atom(1,1)%matl(im1,:,2,1,1))
1312 0 : call wrtout(std_out,mat_writing,'COLL'); call wrtout(ab_out,mat_writing,'COLL')
1313 : enddo
1314 : endif
1315 :
1316 :
1317 0 : mat_writing = ""
1318 :
1319 0 : if (me.eq.0) then
1320 : !print operwan in the real space, in a file
1321 0 : mat_writing = trim(dtfil%filnam_ds(4))//"_wannierocc"
1322 0 : convert = 1
1323 0 : call print_operwan(wan,operwan,trim(mat_writing),convert)
1324 : end if
1325 : ! destroy operators and the occupation matrix
1326 0 : ABI_FREE(eigenks)
1327 0 : ABI_FREE(identityks)
1328 0 : call destroy_operwan(wan,operwan)
1329 0 : ABI_FREE(operwan)
1330 0 : call destroy_operwan_realspace(wan,operocc)!!Destroy the occupation matrix
1331 : endif
1332 :
1333 :
1334 : !! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1335 :
1336 :
1337 :
1338 :
1339 :
1340 :
1341 :
1342 :
1343 : !! -------------------------------------------------------------
1344 : !! TO COMPUTE THE ENERGY MATRIX
1345 : !! -------------------------------------------------------------
1346 : !! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1347 :
1348 : ! Initialize an empty Wannier operator
1349 0 : ABI_MALLOC(operwan,(wan%nkpt,wan%natom_wan,wan%natom_wan))
1350 0 : call initialize_operwan(wan,operwan)
1351 :
1352 : ! Creation of the KS occupation operator
1353 0 : ABI_MALLOC(operks,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%bandf_wan-wan%bandi_wan+1,wan%nsppol))
1354 0 : operks = czero
1355 0 : do isppol = 1,wan%nsppol
1356 0 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
1357 0 : ibandc = iband1 + wan%bandi_wan - 1
1358 0 : do ikpt = 1,wan%nkpt
1359 0 : operks(ikpt,iband1,iband1,isppol) = eigen(((ikpt-1)*dtset%mband+ibandc+(isppol-1)*wan%nkpt*dtset%mband))-fermie
1360 : end do
1361 : end do
1362 : end do
1363 :
1364 : ! Compute the energy in wannier basis
1365 0 : do ikpt = 1,wan%nkpt
1366 0 : call compute_oper_ks2wan(wan,operks,operwan,ikpt)
1367 : end do
1368 : !!In operwan, energies level are stored (shifted with fermi level) in Hartree
1369 0 : ABI_FREE(operks)
1370 :
1371 :
1372 : ! check that the eigenvalues are real
1373 0 : do ikpt = 1,wan%nkpt
1374 0 : do iatom1 = 1,wan%natom_wan
1375 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1376 0 : do m1= 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1377 0 : do isppol = 1,wan%nsppol
1378 0 : if (aimag(operwan(ikpt,iatom1,iatom1)%atom(il1,il1)%matl(m1,m1,isppol,1,1)) > 1d-8) then
1379 0 : write(mat_writing,'(a)') " An eigenvalue has an imaginary part: ikpt, atom, l, m, isppol, value"
1380 0 : write(mat_writing2,'(i0,i0,i0,i0,i0,E15.6)') ikpt, iatom1, il1, im1, isppol, &
1381 0 : & operwan(ikpt,iatom1,iatom1)%atom(il1,il1)%matl(m1,m1,isppol,1,1)
1382 0 : ABI_ERROR(message)
1383 :
1384 : end if
1385 : end do
1386 : end do
1387 : end do
1388 : end do
1389 : end do
1390 :
1391 :
1392 :
1393 :
1394 :
1395 : !! -------------------------------------------------------------
1396 : !! Transform the Wannier operator in real space (in eV)
1397 : !! 1) allocate operwan_realspace
1398 : !! -------------------------------------------------------------
1399 0 : if (dtset%plowan_realspace >= 1) then
1400 0 : call init_operwan_realspace(wan,operwan_realspace)
1401 : endif
1402 :
1403 : ! ABI_MALLOC(operwan_realspace,(wan%natom_wan,wan%natom_wan))
1404 : ! do iatom1 = 1,wan%natom_wan
1405 : ! do iatom2 = 1,wan%natom_wan
1406 : ! n1=size(wan%nposition(iatom1)%pos,1)
1407 : ! n2=size(wan%nposition(iatom2)%pos,1)
1408 : ! ABI_MALLOC(operwan_realspace(iatom1,iatom2)%position,(n1,n2))
1409 : ! do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1410 : ! do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1411 : ! n1=wan%nbl_atom_wan(iatom1)
1412 : ! n2=wan%nbl_atom_wan(iatom2)
1413 : ! ABI_MALLOC(operwan_realspace(iatom1,iatom2)%position(pos1,pos2)%atom,(n1,n2))
1414 : ! do il1 = 1,wan%nbl_atom_wan(iatom1)
1415 : ! do il2 = 1,wan%nbl_atom_wan(iatom2)
1416 : ! n1=2*wan%latom_wan(iatom1)%lcalc(il1)+1
1417 : ! n2=2*wan%latom_wan(iatom2)%lcalc(il2)+1
1418 : ! ABI_MALLOC(operwan_realspace(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl,(n1,n2,wan%nsppol,1,1))
1419 : ! operwan_realspace(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl = zero
1420 : ! end do
1421 : ! end do
1422 : ! end do
1423 : ! end do
1424 : ! end do
1425 : ! end do
1426 :
1427 : !! -------------------------------------------------------------
1428 : !! Transform the Wannier operator in real space (in eV)
1429 : !! 2) compute the value in real space (only if kptopt>0 ie BZ correctly sampled)
1430 : !! -------------------------------------------------------------
1431 0 : if (dtset%plowan_realspace >= 1 .and. dtset%kptopt > 0 ) then ! interpolation and kptopt >0 : compute Wannier functions in real space.
1432 0 : call compute_oper_wank2realspace(wan,operwan,operwan_realspace)
1433 : ! do isppol = 1,wan%nsppol
1434 : ! do iatom1 = 1,wan%natom_wan
1435 : ! do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1436 : ! do il1 = 1,wan%nbl_atom_wan(iatom1)
1437 : ! do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1438 : ! do iatom2 = 1,wan%natom_wan
1439 : ! do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1440 : ! do il2 = 1,wan%nbl_atom_wan(iatom2)
1441 : ! do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1442 : ! !sum over ikpt
1443 : ! do ikpt = 1,wan%nkpt
1444 : ! operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,isppol,1,1) =&
1445 : ! operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)&
1446 : ! + real(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)&
1447 : ! * wan%wtk(ikpt) * exp( cmplx(0.0,1.0) * two_pi * ( &
1448 : ! wan%kpt(1,ikpt) * ( wan%nposition(iatom1)%pos(pos1,1) - wan%nposition(iatom2)%pos(pos2,1) )+&
1449 : ! wan%kpt(2,ikpt) * ( wan%nposition(iatom1)%pos(pos1,2) - wan%nposition(iatom2)%pos(pos2,2) )+&
1450 : ! wan%kpt(3,ikpt) * ( wan%nposition(iatom1)%pos(pos1,3) - wan%nposition(iatom2)%pos(pos2,3)))))
1451 : ! end do
1452 : ! !end of the sum
1453 : ! end do
1454 : ! end do
1455 : ! end do
1456 : ! end do
1457 : ! end do
1458 : ! end do
1459 : ! end do
1460 : ! end do
1461 : ! end do
1462 :
1463 :
1464 :
1465 : !! -------------------------------------------------------------
1466 : !! Transform the Wannier operator in real space (in eV)
1467 : !! 3) write operwan_realspace in a file unit owrunt
1468 : !! -------------------------------------------------------------
1469 :
1470 :
1471 0 : write(message,'(4a)') ch10,&
1472 0 : & ' == Write hamiltonian in real space Wannier function to file ',trim(owrfile),' =='
1473 0 : call wrtout(std_out,message,'COLL')
1474 0 : if (me.eq.0) then
1475 0 : if (open_file(owrfile, message, newunit=owrunt, form="unformatted", status="unknown", action="write") /= 0) then
1476 0 : ABI_ERROR(message)
1477 : end if
1478 0 : rewind(owrunt)
1479 0 : do isppol = 1,wan%nsppol
1480 0 : do iatom1 = 1,wan%natom_wan
1481 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1482 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1483 0 : do iatom2 = 1,wan%natom_wan
1484 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1485 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1486 0 : write(owrunt) operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(:,:,isppol,1,1)
1487 : end do
1488 : end do
1489 : end do
1490 : end do
1491 : end do
1492 : end do
1493 : end do
1494 0 : close(owrunt)
1495 : end if
1496 : end if ! plotwan_realspace>0 and kptopt>0
1497 :
1498 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1499 : ! To change until green study (use operwan_realspace instead of the print subroutine)
1500 : !! BA?
1501 :
1502 :
1503 0 : if (dtset%plowan_realspace == 1) then ! We print the matrix of energy in eV
1504 0 : if (me.eq.0) then
1505 : !print operwan in the real space, in a file
1506 0 : mat_writing = trim(dtfil%filnam_ds(4))//"_wanniereigen"
1507 0 : convert = 27.2107
1508 0 : call print_operwan(wan,operwan,trim(mat_writing),convert)
1509 : end if
1510 : end if
1511 :
1512 :
1513 :
1514 :
1515 :
1516 : !!========================================================================================
1517 : !! Computation of the interaction ( sum(l and l') t_ll' )
1518 : !! This interaction is meaningful in the real space, a transformation is made in this loop
1519 : !!========================================================================================
1520 :
1521 : if (3==4.and.prtint .eq. 1 .and. dtset%plowan_realspace == 1 ) then
1522 : if (open_file(trim(dtfil%filnam_ds(4))//'_inter',message, newunit=unt) /= 0) then
1523 : ABI_ERROR(message)
1524 : end if
1525 : write(unt,'(a,i0,a,F7.3)') "# nsppol = ",wan%nsppol," and acell = ",dtset%acell_orig(1,1)
1526 : write(unt,'(a)') "# Interaction between an orbital and another (in Hartree) : isppol iatom1 pos1 iproj1 iatom2 pos2 iproj2 value"
1527 :
1528 : !!!to compute all interactions
1529 : ! do isppol = 1,wan%nsppol
1530 : ! do iatom1 = 1,wan%natom_wan
1531 : ! do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1532 : ! do il1 = 1,wan%nbl_atom_wan(iatom1)
1533 : ! do iatom2 = 1,wan%natom_wan
1534 : ! do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1535 : ! do il2 = 1,wan%nbl_atom_wan(iatom2)
1536 : ! if (iatom1 .ne. iatom2 .or. il1 .ne. il2 .or. pos1 .ne. pos2) then !not the same orbital on the same atom
1537 : ! if (iatom1 .lt. iatom2 .or. iatom1 .eq. iatom2 .and. (pos1 .lt. pos2 .or. (pos1 .eq. pos2 .and. il1 .lt. il2))) then ! to print only once each interac tion
1538 : ! sum = 0
1539 : ! do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1540 : ! do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1541 : ! sum2 = 0
1542 : ! do ikpt = 1,wan%nkpt
1543 : ! sum2 = sum2 + abs(real(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)*wan%wtk(ikpt)*exp(cmplx(0.0,1.0)*two_pi*(wan%kpt(1,ikpt)*(wan%nposition(iatom1)%pos(pos1,1)-wan%nposition(iatom2)%pos(pos2,1))+wan%kpt(2,ikpt)*(wan%nposition(iatom1)%pos(pos1,2)-wan%nposition(iatom2)%pos(pos2,2))+wan%kpt(3,ikpt)*(wan%nposition(iatom1)%pos(pos1,3)-wan%nposition(iatom2)%pos(pos2,3))))))
1544 : ! end do
1545 : ! sum = sum + sum2
1546 : ! end do
1547 : ! end do
1548 : ! write(unt,'(i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,E15.6)') isppol," ",iatom1," ",pos1," ",wan%projector_wan(iatom1)%lproj(il1)," ",iatom2," ",pos2," ",wan%projector_wan(iatom2)%lproj(il2) ,sum
1549 : ! end if
1550 : ! end if
1551 : ! end do
1552 : ! end do
1553 : ! end do
1554 : ! end do
1555 : ! end do
1556 : ! end do
1557 : ! end do
1558 :
1559 :
1560 :
1561 : do isppol = 1,wan%nsppol
1562 : do iatom2 = 1,wan%natom_wan
1563 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1564 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1565 : sum3 = 0
1566 : do im1 = 1,7 ! for f orbitals
1567 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1568 : sum = 0
1569 : sum2 = 0
1570 : do ikpt = 1,wan%nkpt
1571 : sum2 = sum2 + real(operwan(ikpt,1,iatom2)%atom(1,il2)%matl(im1,im2,isppol,1,1)*wan%wtk(ikpt)*&
1572 : & exp(cmplx(0.0,1.0)*two_pi*( &
1573 : & wan%kpt(1,ikpt)*(wan%nposition(1)%pos(1,1)-wan%nposition(iatom2)%pos(pos2,1))+ &
1574 : & wan%kpt(2,ikpt)*(wan%nposition(1)%pos(1,2)-wan%nposition(iatom2)%pos(pos2,2))+ &
1575 : & wan%kpt(3,ikpt)*(wan%nposition(1)%pos(1,3)-wan%nposition(iatom2)%pos(pos2,3)) )))
1576 : sum = sum + real(operwan(ikpt,iatom2,iatom2)%atom(il2,il2)%matl(im2,im2,isppol,1,1)*wan%wtk(ikpt))
1577 : end do
1578 : sum2 = sum2**2
1579 : sum3 = sum3 + sum2/sum
1580 : end do
1581 : end do
1582 : write(unt,'(i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,E15.6)') isppol,&
1583 : & " ",1," ",1," ",7," ",iatom2," ",pos2," ",wan%projector_wan(iatom2)%lproj(il2) ,sum3
1584 : end do
1585 : end do
1586 : end do
1587 : end do
1588 : close(unt)
1589 : end if
1590 :
1591 :
1592 : !! transformation back in reciprocal space with a limited number of neighbors
1593 : !!==========================================================================
1594 0 : if (dtset%plowan_realspace == 2) then
1595 : !read from the file
1596 0 : write(message,'(4a)') ch10,&
1597 0 : & ' == Read hamiltonian in real space Wannier function on file ',trim(owrfile),' =='
1598 0 : call wrtout(std_out,message,'COLL')
1599 :
1600 0 : if (open_file(owrfile, message, newunit=owrunt, form="unformatted", status="old", action="read") /= 0) then
1601 0 : ABI_ERROR(message)
1602 : end if
1603 0 : rewind(owrunt)
1604 0 : do isppol = 1,wan%nsppol
1605 0 : do iatom1 = 1,wan%natom_wan
1606 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1607 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1608 0 : do iatom2 = 1,wan%natom_wan
1609 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1610 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1611 0 : read(owrunt) operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(:,:,isppol,1,1)
1612 : end do
1613 : end do
1614 : end do
1615 : end do
1616 : end do
1617 : end do
1618 : end do
1619 0 : close(owrunt)
1620 :
1621 : !-----------------------------------------------------------------
1622 : !Set the xx' interaction to 0
1623 : ! In order to do a Wannier interpolation without a given number of
1624 : ! real space terms in the wannier hamiltonian.
1625 : !-----------------------------------------------------------------
1626 : if(3==4) then
1627 : do isppol = 1,wan%nsppol
1628 : do pos2 = 2,size(wan%nposition(2)%pos,1)
1629 : operwan_realspace%atom_index(1,2)%position(1,pos2)%atom(1,1)%matl(:,:,isppol,1,1) = zero
1630 : operwan_realspace%atom_index(2,1)%position(pos2,1)%atom(1,1)%matl(:,:,isppol,1,1) = zero
1631 : ! atoms 1 and 2 are selected for removal
1632 : ! position(1,pos2): select cell 1 and all other cells pos2.
1633 : ! atom(1,1): selected index 1 of atom1 and index 1 of atom2
1634 : ! mat1(:,:,isppol,1,1): remove all ml terms.
1635 : operwan_realspace%atom_index(1,1)%position(1,pos2)%atom(1,2)%matl(:,:,isppol,1,1) = zero
1636 : end do
1637 : end do
1638 : write(message,'(2a)') ' == Block suppressed in the real space Wannier hamiltonian'
1639 : call wrtout(std_out,message,'COLL')
1640 : endif
1641 :
1642 :
1643 :
1644 : !-----------------------------------------------------------------
1645 : ! Perform Wannier transform from real space to reciprocal space
1646 : !-----------------------------------------------------------------
1647 0 : write(message,'(2a)') ' == Perform Wannier transform from real space to reciprocal space =='
1648 0 : call wrtout(std_out,message,'COLL')
1649 0 : do isppol = 1,wan%nsppol
1650 0 : do ikpt = 1,wan%nkpt
1651 0 : do iatom1 = 1,wan%natom_wan
1652 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1653 0 : do iatom2 = 1,wan%natom_wan
1654 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1655 0 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(:,:,isppol,1,1) = zero
1656 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1657 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1658 : !sum over neigbours
1659 0 : pos1 = 1
1660 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1661 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1) =&
1662 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)&
1663 : +real(operwan_realspace%atom_index(iatom1,iatom2)%&
1664 : &position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)&
1665 : * exp( - cmplx(0.0,1.0) * two_pi * ( &
1666 : wan%kpt(1,ikpt) * ( wan%nposition(iatom1)%pos(pos1,1) - wan%nposition(iatom2)%pos(pos2,1) )+&
1667 : wan%kpt(2,ikpt) * ( wan%nposition(iatom1)%pos(pos1,2) - wan%nposition(iatom2)%pos(pos2,2) )+&
1668 0 : wan%kpt(3,ikpt) * ( wan%nposition(iatom1)%pos(pos1,3) - wan%nposition(iatom2)%pos(pos2,3)))))
1669 : end do
1670 : !end do
1671 : !end of the sum
1672 : end do
1673 : end do
1674 : end do
1675 : end do
1676 : end do
1677 : end do
1678 : end do
1679 : end do
1680 :
1681 : ! plowan_realspace==2 (Wannier interpolation)
1682 :
1683 : !! -------------------------------------------------------------
1684 : !! End Transform the Wannier operator in real space (in eV)
1685 : !! n) deallocate operwan_realspace
1686 : !! ------------------------------------------------------------
1687 : ! do iatom1 = 1,wan%natom_wan
1688 : ! do iatom2 = 1,wan%natom_wan
1689 : ! do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
1690 : ! do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
1691 : ! do il1 = 1,wan%nbl_atom_wan(iatom1)
1692 : ! do il2 = 1,wan%nbl_atom_wan(iatom2)
1693 : ! ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl)
1694 : ! end do
1695 : ! end do
1696 : ! ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom)
1697 : ! end do
1698 : ! end do
1699 : ! ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position)
1700 : ! end do
1701 : ! end do
1702 : ! ABI_FREE(operwan_realspace%atom_index)
1703 : endif
1704 0 : if (dtset%plowan_realspace >= 1) then
1705 0 : call destroy_operwan_realspace(wan,operwan_realspace)
1706 : endif
1707 :
1708 : ! ----------------------------------------------------------------------------------------
1709 : ! Here each block of the hamiltonian matrix in Wannier basis is diagonalized separately
1710 : ! ----------------------------------------------------------------------------------------
1711 :
1712 : ! off diagonal blocks are suppressed in the hamiltonian matrix before diagonalisation
1713 : ! whole_diag=-1
1714 0 : if (whole_diag .eq. 0) then
1715 0 : write(message,'(a,i5,a)')' plowan_compute =',dtset%plowan_compute,&
1716 0 : & ' Off diag blocks are suppressed in the Wannier hamiltonian before diagonalisation'
1717 0 : call wrtout(std_out,message,'COLL')
1718 : ! !To diagonalize the block matrix for each orbital
1719 0 : do iatom1 = 1,wan%natom_wan
1720 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1721 0 : do isppol = 1,wan%nsppol
1722 0 : count = 2*wan%latom_wan(iatom1)%lcalc(il1)+1
1723 0 : ABI_MALLOC(matrix_to_diag,(count,count))
1724 0 : ABI_MALLOC(eig,(count))
1725 0 : ABI_MALLOC(rwork,(3*count-2))
1726 0 : lwork = 65*count !Value to optimize the diagonalization
1727 0 : ABI_MALLOC(zwork,(lwork))
1728 0 : do ikpt = 1,wan%nkpt
1729 0 : matrix_to_diag(:,:) = operwan(ikpt,iatom1,iatom1)%atom(il1,il1)%matl(:,:,isppol,1,1)
1730 0 : call zheev('v','u',count,matrix_to_diag,count,eig,zwork,lwork,rwork,info)
1731 0 : if (info .eq. 0) then !!Correct diagonalization
1732 0 : matrix_to_diag = zero
1733 0 : do im1 = 1,count
1734 0 : matrix_to_diag(im1,im1) = eig(im1)
1735 : end do
1736 0 : operwan(ikpt,iatom1,iatom1)%atom(il1,il1)%matl(:,:,isppol,1,1) = matrix_to_diag(:,:)
1737 : else
1738 0 : write(message,'(a)') "Error in the normalization of the Wannier eigenvalues" ! BA?
1739 0 : ABI_ERROR(message)
1740 : end if
1741 : end do
1742 0 : ABI_FREE(matrix_to_diag)
1743 0 : ABI_FREE(eig)
1744 0 : ABI_FREE(rwork)
1745 0 : ABI_FREE(zwork)
1746 : end do
1747 : end do
1748 : end do
1749 : end if
1750 :
1751 :
1752 :
1753 : ! ----------------------------------------------------------------------------------------
1754 : ! Here the hamiltonian matrix in Wannier basis is diagonalized completely
1755 : ! ----------------------------------------------------------------------------------------
1756 :
1757 : if (whole_diag .eq. 1) then
1758 : ! To diagonalize the whole matrix
1759 0 : count = wan%nspinor*wan%size_wan
1760 0 : ABI_MALLOC(matrix_to_diag,(count,count))
1761 0 : ABI_MALLOC(eig,(count))
1762 0 : ABI_MALLOC(rwork,(3*count-2))
1763 0 : lwork = 65*count ! Value to optimize speed of the diagonalization
1764 0 : ABI_MALLOC(zwork,(lwork))
1765 : !First, write operwan matrix in an inversible matrix
1766 0 : do isppol = 1,wan%nsppol
1767 0 : do ikpt = 1,wan%nkpt
1768 0 : matrix_to_diag = czero
1769 0 : index_l = 0
1770 0 : do iatom1 = 1,wan%natom_wan
1771 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1772 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1773 0 : index_l = index_l + 1 ! the line changes
1774 0 : index_c = 1 ! index_c is set to one each time the line changes
1775 0 : do iatom2 = 1,wan%natom_wan
1776 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1777 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1778 0 : matrix_to_diag(index_l,index_c) = operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)
1779 0 : index_c = index_c + 1
1780 : end do !im2
1781 : end do !il2
1782 : end do ! iatom2 (the line changes)
1783 : end do ! im1
1784 : end do ! il1
1785 : end do !iatom1
1786 :
1787 : !Then, invert the matrix
1788 0 : call zheev('v','u',count,matrix_to_diag,count,eig,zwork,lwork,rwork,info)
1789 0 : if (info .eq. 0) then ! Correct diagonalization
1790 0 : matrix_to_diag = czero
1791 0 : do im1 = 1,count
1792 0 : matrix_to_diag(im1,im1) = eig(im1)
1793 : end do
1794 : else
1795 0 : write(message,'(a)') "Error in the normalization of the Wannier eigenvalues" ! BA?
1796 0 : ABI_ERROR(message)
1797 : end if
1798 : !Finally, we write the value diagonalized back into operwan
1799 0 : index_l = 0
1800 0 : do iatom1 = 1,wan%natom_wan
1801 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
1802 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
1803 0 : index_l = index_l + 1
1804 0 : index_c = 1
1805 0 : do iatom2 = 1,wan%natom_wan
1806 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
1807 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
1808 0 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1) = matrix_to_diag(index_l,index_c)
1809 0 : index_c = index_c + 1
1810 : end do
1811 : end do
1812 : end do
1813 : end do
1814 : end do
1815 : end do
1816 : !ikpt/isppol
1817 : end do
1818 : end do
1819 0 : ABI_FREE(matrix_to_diag)
1820 0 : ABI_FREE(eig)
1821 0 : ABI_FREE(rwork)
1822 0 : ABI_FREE(zwork)
1823 : end if
1824 :
1825 : !------------------------------------------------------------------------
1826 : ! we write the band structure of the atoms in the Wannier basis (in eV)
1827 : ! Warning: In the case of the diagonalisation of the whole hamltonien
1828 : ! the band structure is separated in as many files as atoms.
1829 : !------------------------------------------------------------------------
1830 : if (band_struct .eq. 1 ) then
1831 :
1832 0 : write(message,'(3a)') ch10,&
1833 0 : & ' == For each k-point of the path, gives the eigenvalues (in eV) of the Hamiltonian in the Wannier basis'
1834 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
1835 0 : write(message,'(2a,f13.4,a)') ch10,&
1836 0 : & ' (The band structure is shifted by fermie =',fermie*Ha_eV,' eV )'
1837 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
1838 :
1839 0 : write(message,'(a,i10)') ' == Number of atoms ',wan%natom_wan
1840 0 : do iatom1 = 1,wan%natom_wan
1841 0 : i2s = '(I0)' ! trick to add the atom number
1842 0 : write(x1,i2s) iatom1 ! at the end of the filename
1843 0 : if (wan%nsppol .eq. 1 .and. me.eq.0 ) then
1844 0 : if (open_file(trim(dtfil%filnam_ds(4))//"_BANDSTRUCT"//trim(x1),message,newunit=unt) /= 0) then
1845 0 : ABI_ERROR(message)
1846 : end if
1847 0 : write(unt,'(a,i0)') "#Wannier band structure for the atom ",iatom1
1848 0 : write(message,'(2a,i2)') ch10," Wannier band structure for atom ",iatom1
1849 0 : call wrtout(std_out,message,'COLL') ; call wrtout(ab_out,message,'COLL')
1850 0 : do ikpt = 1,wan%nkpt
1851 0 : mat_writing = "#ikpt ="
1852 0 : write(mat_writing2,'(i0)') ikpt
1853 0 : mat_writing = trim(mat_writing)//" "//trim(mat_writing2)
1854 0 : write(unt,*) trim(mat_writing)
1855 0 : write(mat_writing,'(i0)') ikpt
1856 0 : write(mat_writing_out,'(i0)') ikpt
1857 0 : do l1 = 1,wan%nbl_atom_wan(iatom1)
1858 0 : do m1 = 1,2*wan%latom_wan(iatom1)%lcalc(l1)+1
1859 0 : write (mat_writing2,'(F25.7)') 27.2107*real(operwan(ikpt,iatom1,iatom1)%atom(l1,l1)%matl(m1,m1,1,1,1))
1860 0 : write (mat_writing2_out,'(F12.3)') 27.2107*real(operwan(ikpt,iatom1,iatom1)%atom(l1,l1)%matl(m1,m1,1,1,1))
1861 0 : mat_writing = trim(mat_writing)//trim(mat_writing2)
1862 0 : mat_writing_out = trim(mat_writing_out)//trim(mat_writing2_out)
1863 : end do
1864 : end do
1865 0 : write(unt,*) trim(mat_writing)
1866 0 : write(ab_out,*) trim(mat_writing_out)
1867 0 : write(std_out,*) trim(mat_writing_out)
1868 : end do
1869 0 : close(unt)
1870 0 : else if (wan%nsppol .eq. 2) then
1871 0 : if (open_file(trim(dtfil%filnam_ds(4))//"_BANDSTRUCTUP"//trim(x1),message,newunit=unt) /= 0) then
1872 0 : ABI_ERROR(message)
1873 : end if
1874 0 : if (open_file(trim(dtfil%filnam_ds(4))//"_BANDSTRUCTDN"//trim(x1),message,newunit=unt2) /= 0) then
1875 0 : ABI_ERROR(message)
1876 : end if
1877 0 : write(unt,'(a,i0,a)') "#Wannier band structure for the atom ",iatom1, " polarization up"
1878 0 : write(unt2,'(a,i0,a)') "#Wannier band structure for the atom ",iatom1, " polarization down"
1879 0 : do isppol = 1,wan%nsppol
1880 0 : do ikpt = 1,wan%nkpt
1881 0 : mat_writing = "#ikpt ="
1882 0 : write(mat_writing2,'(i0)') ikpt
1883 0 : mat_writing = trim(mat_writing)//" "//trim(mat_writing2)
1884 0 : if (isppol .eq. 1) then
1885 0 : write(unt,*) trim(mat_writing)
1886 : else
1887 0 : write(unt2,*) trim(mat_writing)
1888 : end if
1889 0 : write(mat_writing,'(i0)') ikpt
1890 0 : do l1 = 1,wan%nbl_atom_wan(iatom1)
1891 0 : do m1 = 1,2*wan%latom_wan(iatom1)%lcalc(l1)+1
1892 0 : write (mat_writing2,'(F12.7)') 27.2107*real(operwan(ikpt,iatom1,iatom1)%atom(l1,l1)%matl(m1,m1,isppol,1,1))
1893 0 : mat_writing = trim(mat_writing)//trim(mat_writing2)
1894 : end do
1895 : end do
1896 0 : if (isppol .eq. 1) then
1897 0 : write(unt,*) trim(mat_writing)
1898 : else
1899 0 : write(unt2,*) trim(mat_writing)
1900 : end if
1901 0 : write(ab_out,*) trim(mat_writing)
1902 0 : write(std_out,*) trim(mat_writing)
1903 : end do
1904 : end do
1905 0 : close(unt)
1906 0 : close(unt2)
1907 : end if
1908 : end do
1909 : end if
1910 :
1911 :
1912 : !! ----------------------------------------------
1913 : !! GREEN STUDY
1914 : !! ----------------------------------------------
1915 :
1916 : !!Here we only study the first atom
1917 :
1918 : if (plowan_computegreen .eq. 1 .and. me.eq.0 ) then
1919 : if (dos .ge. 1) then ! compute partial DOS for l=dos
1920 : if (open_file(trim(dtfil%filnam_ds(4))//"_dosfromgreen",message,newunit=dos_unt) /= 0) then
1921 : ABI_ERROR(message)
1922 : end if
1923 : write(dos_unt,'(a)') "#DOS function for the first atom computed with the Green function"
1924 : write(dos_unt,'(a,i0,a,i0)') "# l = ",wan%latom_wan(1)%lcalc(dos)," bands ; nsppol = ",wan%nsppol
1925 : write(dos_unt,'(a)') "# frequency, DOS, isppol"
1926 :
1927 : if(wan%nsppol>=2) then
1928 : if (open_file(trim(dtfil%filnam_ds(4))//"_dosfromgreen_b",message,newunit=dos_unt2) /= 0) then
1929 : ABI_ERROR(message)
1930 : end if
1931 : write(dos_unt2,'(a)') "#DOS function for the first atom computed with the Green function"
1932 : write(dos_unt2,'(a,i0,a,i0)') "# l = ",wan%latom_wan(1)%lcalc(dos)," bands ; nsppol = ",wan%nsppol
1933 : write(dos_unt2,'(a)') "# frequency, DOS, isppol"
1934 : endif
1935 : end if
1936 :
1937 : if (dos .le. -1) then ! compute Hybri for l=|dos|
1938 : if (open_file(trim(dtfil%filnam_ds(4))//"_hybridization",message,newunit=dos_unt) /= 0) then
1939 : ABI_ERROR(message)
1940 : end if
1941 : write(dos_unt,'(a)') "#Hybridization obtained from the green function for the bands selected"
1942 : write(dos_unt,'(a,i0,a,i0)') "# l = ",wan%latom_wan(1)%lcalc(abs(dos))," bands ; nsppol = ",wan%nsppol
1943 : write(dos_unt,'(a)') "#isppol, frequency, F(m=0), F(m=1), F(m=2) ..."
1944 : if(wan%nsppol>=2) then
1945 : if (open_file(trim(dtfil%filnam_ds(4))//"_hybridization_b",message, newunit=dos_unt2) /= 0) then
1946 : ABI_ERROR(message)
1947 : end if
1948 : write(dos_unt2,'(a)') "#Hybridization obtained from the green function for the bands selected"
1949 : write(dos_unt2,'(a,i0,a,i0)') "# l = ",wan%latom_wan(1)%lcalc(abs(dos))," bands ; nsppol = ",wan%nsppol
1950 : write(dos_unt2,'(a)') "#isppol, frequency, F(m=0), F(m=1), F(m=2) ..."
1951 : endif
1952 : end if
1953 :
1954 : ! Method 0
1955 : sizem = 2*wan%latom_wan(1)%lcalc(abs(dos))+1 !number of m for the l orbital we want
1956 : ABI_MALLOC(energies,(sizem,wan%nsppol))
1957 : energies = czero
1958 : ABI_MALLOC(Fff,(2))
1959 : Fff = czero
1960 :
1961 : !We put the eigenenergies of the first atom in an array to use them later
1962 : do isppol = 1,wan%nsppol
1963 : do ikpt = 1,wan%nkpt
1964 : do m1 = 1,sizem
1965 : energies(m1,isppol) = energies(m1,isppol) +&
1966 : & real(operwan(ikpt,1,1)%atom(abs(dos),abs(dos))%matl(m1,m1,isppol,1,1))*wan%wtk(ikpt)
1967 : end do
1968 : end do
1969 : end do
1970 :
1971 :
1972 : shift = 0
1973 : do il1 = 1,abs(dos)-1
1974 : shift = shift + 2*wan%latom_wan(1)%lcalc(il1)+1!shift for the right choice of indices for the lorbital chosen
1975 : end do
1976 :
1977 : !We destroy the operwan which was used to compute energies before
1978 :
1979 : call destroy_operwan(wan,operwan)
1980 : ABI_FREE(operwan)
1981 :
1982 : !-----------------------------------------------------------
1983 : ! Loop over the frequencies to compute DOS or Hybridization
1984 : !-----------------------------------------------------------
1985 : do iw = 1,number_of_frequencies
1986 : ABI_MALLOC(operwan,(wan%nkpt,wan%natom_wan,wan%natom_wan))
1987 : call initialize_operwan(wan,operwan)
1988 : !!creation of the Green operator
1989 : wcurrent = wbase + (iw-1)*wincrease
1990 : ! if (allocated(operks)) then
1991 : ! ABI_FREE(operks)
1992 : ! endif
1993 : ABI_MALLOC(operks,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%bandf_wan-wan%bandi_wan+1,wan%nsppol))
1994 : operks = czero
1995 :
1996 : ! Fill diagonal elements to have DFT Green's function.
1997 : !------------------------------------------------------
1998 : do isppol = 1,wan%nsppol
1999 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
2000 : ibandc = iband1 + wan%bandi_wan - 1
2001 : do ikpt = 1,wan%nkpt
2002 : operks(ikpt,iband1,iband1,isppol) = &
2003 : & 1d0/(wcurrent-eigen(((ikpt-1)*dtset%mband+ibandc+(isppol-1)*wan%nkpt*dtset%mband))+fermie) ! 1/(w-E(kv))
2004 : end do
2005 : end do
2006 : end do
2007 :
2008 : ! Compute Green's function in wannier basis in recip space.
2009 : !----------------------------------------------------------
2010 : do ikpt = 1,wan%nkpt
2011 : call compute_oper_ks2wan(wan,operks,operwan,ikpt) !in reciprocal space
2012 : end do
2013 : ABI_FREE(operks)
2014 :
2015 :
2016 : ! Transform the operwan into a better shape for inversion
2017 : !----------------------------------------------------------
2018 : ABI_MALLOC(operwansquare,(wan%nkpt,wan%nsppol,wan%nspinor*wan%size_wan,wan%nspinor*wan%size_wan))
2019 : operwansquare = czero
2020 : do ikpt = 1,wan%nkpt
2021 : do isppol = 1,wan%nsppol
2022 : do ispinor1 = 1,wan%nspinor
2023 : do ispinor2 = 1,wan%nspinor
2024 : index_l = 0 !index_l is set to 0 at the beginning
2025 : do iatom1 = 1,wan%natom_wan
2026 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2027 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
2028 : index_l = index_l + 1 !the line changes
2029 : index_c = 1 !counter_c is set to one each time the line changes
2030 : do iatom2 = 1,wan%natom_wan
2031 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2032 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
2033 : operwansquare(ikpt,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+wan%size_wan*(ispinor2-1)) = &
2034 : & operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)
2035 : index_c = index_c + 1
2036 : end do !im2
2037 : end do !il2
2038 : end do !iatom2 (the line changes)
2039 : end do !im1
2040 : end do !il1
2041 : end do !iatom1
2042 : end do
2043 : end do
2044 : end do
2045 : end do
2046 :
2047 : ABI_MALLOC(operwansquarereal,(wan%nsppol,size(operwansquare,3),size(operwansquare,4)))
2048 : operwansquarereal = czero
2049 :
2050 : ! Transformation in the real space (T=T'=0) : compute the local quantities
2051 : !-------------------------------------------------------------------------
2052 : do isppol = 1,wan%nsppol
2053 : do index_l = 1,size(operwansquare,3)
2054 : do index_c = 1,size(operwansquare,4)
2055 : do ikpt = 1,wan%nkpt
2056 : operwansquarereal(isppol,index_l,index_c) = operwansquarereal(isppol,index_l,index_c) + &
2057 : & operwansquare(ikpt,isppol,index_l,index_c)*wan%wtk(ikpt)
2058 : end do
2059 : end do
2060 : end do
2061 : end do
2062 : ABI_FREE(operwansquare)
2063 :
2064 : if (dos .ge. 1) then ! either we compute the DOS (-imaginary part/Pi of the green function in the Wannier basis)
2065 : ! Compute the dos
2066 : !-------------------------------------------------------------------------
2067 : sum = 0
2068 : do isppol = 1,wan%nsppol
2069 : do m1 = 1,2*wan%latom_wan(1)%lcalc(dos)+1
2070 : sum = sum - aimag(operwansquarereal(isppol,shift+m1,shift+m1))
2071 : end do
2072 : if(isppol==1) write(dos_unt,'(F8.3,E15.6)') real(27.2101*wcurrent),sum/(27.2107*3.14159)
2073 : if(isppol==2) write(dos_unt2,'(F8.3,E15.6)') real(27.2101*wcurrent),sum/(27.2107*3.14159)
2074 : end do
2075 :
2076 : else ! either we compute the F part (the residual part) of the invert of the green functions
2077 : ! Compute the hybridization
2078 : !-------------------------------------------------------------------------
2079 :
2080 : do isppol = 1,wan%nsppol
2081 : ABI_MALLOC(matrix_to_diag,(sizem,sizem))
2082 : matrix_to_diag = operwansquarereal(isppol,shift+1:shift+sizem,shift+1:shift+sizem)
2083 : ! attention a isppol ci dessus
2084 : call xginv(matrix_to_diag,sizem)
2085 : operwansquarereal(isppol,shift+1:shift+sizem,shift+1:shift+sizem) = matrix_to_diag
2086 : ABI_FREE(matrix_to_diag)
2087 : mat_writing = ""
2088 : do m1 = 1,sizem
2089 : Fff(isppol) = wcurrent - operwansquarereal(isppol,shift+m1,shift+m1) - energies(m1,isppol)
2090 : write(mat_writing2,'(E15.6)') -aimag(27.2107*Fff(isppol))
2091 : mat_writing = trim(mat_writing)//" "//trim(mat_writing2)
2092 : end do
2093 : if(isppol==1) write(dos_unt,'(F10.3,a)') real(27.2107*wcurrent),trim(mat_writing)
2094 : if(isppol==2) write(dos_unt2,'(F10.3,a)') real(27.2107*wcurrent),trim(mat_writing)
2095 : end do
2096 : end if
2097 : call destroy_operwan(wan,operwan)
2098 : ABI_FREE(operwan)
2099 : ABI_FREE(operwansquarereal)
2100 : end do
2101 : close(dos_unt)
2102 : close(dos_unt2)
2103 : ABI_FREE(energies)
2104 : ABI_FREE(Fff)
2105 :
2106 : end if ! choice of the 1 plowan_computegreen
2107 :
2108 :
2109 :
2110 :
2111 : if (plowan_computegreen .eq. 2 .and. me.eq.0 ) then !! Not working ! not tested, not up to date with the code
2112 : !Method 1
2113 : ABI_MALLOC(energies,(7,wan%nsppol))
2114 : ABI_MALLOC(Ffftable,(7,wan%nsppol))
2115 : Ffftable = czero
2116 : energies = czero
2117 : ! Keep energies for later use
2118 : !----------------------------
2119 : do ikpt = 1,wan%nkpt
2120 : do im1 = 1,7
2121 : do isppol = 1,wan%nsppol
2122 : energies(im1,isppol) = energies(im1,isppol) + real(operwan(ikpt,1,1)%atom(1,1)%matl(im1,im1,isppol,1,1))*wan%wtk(ikpt)
2123 : end do
2124 : end do
2125 : end do
2126 : write(std_out,*) "energies", energies*Ha_eV
2127 :
2128 : ! Loop over frequency
2129 : !----------------------
2130 : do iw = 1,number_of_frequencies
2131 :
2132 : ABI_MALLOC(operwansquare,(wan%nkpt,wan%nsppol,wan%nspinor*wan%size_wan,wan%nspinor*wan%size_wan))
2133 : operwansquare = czero
2134 : wcurrent = wbase + (iw-1)*wincrease
2135 : Ffftable = czero
2136 :
2137 : ! create operwansquare
2138 : !----------------------------
2139 : do ikpt = 1,wan%nkpt
2140 : do isppol = 1,wan%nsppol
2141 : do ispinor1 = 1,wan%nspinor
2142 : do ispinor2 = 1,wan%nspinor
2143 : index_l = 0 ! index_l is set to 0 at the beginning
2144 : do iatom1 = 1,wan%natom_wan
2145 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2146 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
2147 : index_l = index_l + 1 ! the line changes
2148 : index_c = 1 ! counter_c is set to one each time the line changes
2149 : do iatom2 = 1,wan%natom_wan
2150 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2151 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
2152 : operwansquare(ikpt,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+wan%size_wan*(ispinor2-1)) &
2153 : & = operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)
2154 : index_c = index_c + 1
2155 : end do !im2
2156 : end do !il2
2157 : end do ! iatom2 (the line changes)
2158 : end do ! im1
2159 : end do ! il1
2160 : end do !iatom1
2161 : end do
2162 : end do
2163 : end do
2164 :
2165 : ! Create inverse of Green's function
2166 : !-----------------------------------
2167 : do isppol = 1,wan%nsppol
2168 : do im1 = 1,size(operwansquare,3)
2169 : do im2 = 1,size(operwansquare,4)
2170 : if (im1 .eq. im2) then
2171 : operwansquare(ikpt,isppol,im1,im2) = wcurrent-operwansquare(ikpt,isppol,im1,im2)
2172 : else
2173 : operwansquare(ikpt,isppol,im1,im2) = -operwansquare(ikpt,isppol,im1,im2)
2174 : end if
2175 : end do
2176 : end do
2177 : end do
2178 :
2179 :
2180 : ! Create Green's function
2181 : !-----------------------------------
2182 : ABI_MALLOC(matrix_to_diag,(size(operwansquare,3),size(operwansquare,3)))
2183 : do isppol = 1,wan%nsppol
2184 : matrix_to_diag = czero
2185 : matrix_to_diag = operwansquare(ikpt,isppol,:,:)
2186 : call xginv(matrix_to_diag,size(matrix_to_diag,1))
2187 : operwansquare(ikpt,isppol,:,:) = matrix_to_diag
2188 : end do
2189 : ABI_FREE(matrix_to_diag)
2190 :
2191 :
2192 : !! if (dos .eq. 0) then
2193 : !! !select f bands
2194 : !! do isppol = 1,wan%nsppol
2195 : !! if (allocated(matrix_to_diag)) ABI_FREE(matrix_to_diag)
2196 : !! ABI_MALLOC(matrix_to_diag,(7,7))
2197 : !! matrix_to_diag = czero
2198 : !! matrix_to_diag = operwansquare(ikpt,isppol,1:7,1:7)
2199 : !! call xginv(matrix_to_diag,size(matrix_to_diag,1))
2200 :
2201 : !! operwansquare(ikpt,isppol,1:7,1:7) = matrix_to_diag
2202 : !! end do
2203 : !! end if
2204 :
2205 : !! do im1 = 1,7
2206 : !! do isppol = 1,wan%nsppol
2207 : !! if (dos .eq. 0) then
2208 : !! Ffftable(im1,isppol) = Ffftable(im1,isppol) + (wcurrent - &
2209 : !! & operwansquare(ikpt,isppol,im1,im1) - energies(im1,isppol))*wan%wtk(ikpt)
2210 : !! end if
2211 : !! if (dos .eq. 1) then
2212 : !! Ffftable(im1,isppol) = Ffftable(im1,isppol) + operwansquare(ikpt,isppol,im1,im1)*wan%wtk(ikpt)
2213 : !! end if
2214 : !! end do
2215 : !! end do
2216 :
2217 : end do !!loop ikpt
2218 :
2219 : ABI_MALLOC(operwansquarereal,(wan%nsppol,size(operwansquare,3),size(operwansquare,4)))
2220 : operwansquarereal = czero
2221 :
2222 : ! Compute local Green's function
2223 : !--------------------------------------------
2224 : do isppol = 1,wan%nsppol
2225 : do index_l = 1,size(operwansquare,3)
2226 : do index_c = 1,size(operwansquare,4)
2227 : do ikpt = 1,wan%nkpt
2228 : operwansquarereal(isppol,index_l,index_c) = operwansquarereal(isppol,index_l,index_c)&
2229 : & + operwansquare(ikpt,isppol,index_l,index_c)*wan%wtk(ikpt)
2230 : end do
2231 : end do
2232 : end do
2233 : end do
2234 :
2235 : ! Inverse Local Correlated Green's function
2236 : !--------------------------------------------
2237 : if (dos < 0) then
2238 : !select f bands
2239 : do isppol = 1,wan%nsppol
2240 : ABI_MALLOC(matrix_to_diag,(7,7))
2241 : matrix_to_diag = operwansquarereal(isppol,1:7,1:7)
2242 : call xginv(matrix_to_diag,size(matrix_to_diag,1))
2243 : operwansquarereal(isppol,1:7,1:7) = matrix_to_diag
2244 : ABI_FREE(matrix_to_diag)
2245 : end do
2246 : write(268,*) 27.2107*real(wcurrent),27.2107*real(operwansquarereal(1,1,1)),27.2107*aimag(operwansquarereal(1,1,1))
2247 :
2248 : end if
2249 :
2250 : do im1 = 1,7
2251 : do isppol = 1,wan%nsppol
2252 : ! Compute hybridization
2253 : !--------------------------------------------
2254 : if (dos < 0) then
2255 : Ffftable(im1,isppol) = Ffftable(im1,isppol) + wcurrent - operwansquarereal(isppol,im1,im1) - energies(im1,isppol)
2256 : end if
2257 : if (dos > 0) then
2258 : ! Compute Dos
2259 : !--------------------------------------------
2260 : Ffftable(im1,isppol) = Ffftable(im1,isppol) + operwansquarereal(isppol,im1,im1)
2261 : end if
2262 : end do
2263 : end do
2264 : write(269,*) 27.2107*real(wcurrent),27.2107*real(Ffftable(1,1)),27.2107*aimag(Ffftable(1,1))
2265 : write(2699,*) 27.2107*real(wcurrent),27.2107*real(wcurrent - operwansquarereal(1,1,1) - energies(1,1)),&
2266 : & Ha_eV*real(energies(1,1)),Ha_eV*real(operwansquarereal(1,1,1))
2267 :
2268 : if (dos < 0) then
2269 : xsum=czero
2270 : do im1 = 1,7
2271 : if (wan%nsppol .eq. 1) then
2272 : write(std_out,*)'hybri',im1,27.2107*wcurrent,27.2107*Ffftable(im1,1)
2273 : else
2274 : do isppol = 1,wan%nsppol
2275 : write(std_out,*)'hybri',im1,isppol,27.2107*wcurrent,27.2107*Ffftable(im1,isppol)
2276 : xsum=xsum+Ffftable(im1,isppol)
2277 : end do
2278 : end if
2279 : end do
2280 : write(270,*)27.2107*real(wcurrent),27.2107*real(xsum),27.2107*aimag(xsum)
2281 : end if
2282 :
2283 : if (dos > 0) then
2284 : xsum=czero
2285 : do isppol = 1,wan%nsppol
2286 : Ffftable(1,isppol) = Ffftable(1,isppol)+Ffftable(2,isppol)+Ffftable(3,isppol)+Ffftable(4,isppol)&
2287 : & +Ffftable(5,isppol)+Ffftable(6,isppol)+Ffftable(7,isppol)
2288 : if (wan%nsppol .eq. 2) then
2289 : write(std_out,*)'green',isppol,27.2107*wcurrent,Ffftable(1,isppol)/(27.2107*3.14159)
2290 : else
2291 : write(std_out,*)'green',27.2107*wcurrent,Ffftable(1,isppol)/(27.2107*3.14159)
2292 : xsum=xsum+Ffftable(im1,isppol)
2293 : end if
2294 : write(271,*)27.2107*real(wcurrent),real(xsum)/(27.2107*3.14159),aimag(xsum)/(27.2107*3.14159)
2295 : end do
2296 : end if
2297 : ABI_FREE(operwansquare)
2298 : ABI_FREE(operwansquarereal)
2299 : end do !loop frequencies w
2300 : ABI_FREE(energies)
2301 : ABI_FREE(Fff)
2302 :
2303 : call destroy_operwan(wan,operwan)
2304 : ABI_FREE(operwan)
2305 : end if !! choice of the plowan_computegreen
2306 :
2307 : if(plowan_computegreen==0) then
2308 0 : call destroy_operwan(wan,operwan)
2309 0 : ABI_FREE(operwan)
2310 : end if !! choice of the plowan_computegreen
2311 :
2312 : !deallocate temporary cwaveprj/cprj storage
2313 0 : call pawcprj_free(cwaveprj)
2314 0 : ABI_FREE(cwaveprj)
2315 :
2316 :
2317 0 : end subroutine compute_coeff_plowannier
2318 : !!***
2319 :
2320 : !!****f* m_plowannier/print_plowannier
2321 : !! NAME
2322 : !! print_plowannier
2323 : !!
2324 : !! FUNCTION
2325 : !! print the wannier weight (psichi) on a forlb.ovlp file
2326 : !!
2327 : !! INPUTS
2328 : !! dtset%typat,wan
2329 : !!
2330 : !! OUTPUT
2331 : !!
2332 : !! SOURCE
2333 :
2334 :
2335 :
2336 0 : subroutine print_plowannier(wan)
2337 :
2338 : use m_abicore
2339 : use m_io_tools, only : open_file
2340 : use m_specialmsg, only : wrtout
2341 :
2342 : !Arguments-------------------------
2343 : type(plowannier_type),intent(in) :: wan
2344 : !Local variables-------------------
2345 : character(len=500) :: msg
2346 : integer :: unt,iatom,spin,ikpt,iband,ibandc,il,ispinor,im
2347 :
2348 : !Creation of the data.plowann file
2349 0 : if (open_file('data.plowann',msg,newunit=unt,form='formatted',status='replace') /= 0) then
2350 0 : ABI_ERROR(msg)
2351 : end if
2352 0 : rewind(unt)
2353 :
2354 0 : write(msg,'(2a)') ch10,' Print the psichi coefficients in data.plowann'
2355 0 : call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
2356 :
2357 : !Header of the file data.plowann
2358 0 : write(unt,'(a22,i2)')"Total number of atom =", wan%natom_wan
2359 0 : write(unt,*)"List of atoms", wan%iatom_wan(:)
2360 0 : write(unt,'(a7,2i4)')"Bands =",wan%bandi_wan,wan%bandf_wan
2361 0 : write(unt,'(a26,i2)')"Total number of orbitals =",sum(wan%nbl_atom_wan(:))
2362 0 : do iatom=1,wan%natom_wan
2363 0 : write(unt,'(a17,i2,a3,4i2)')"Orbitals for atom",wan%iatom_wan(iatom)," = ",wan%latom_wan(iatom)%lcalc(:)
2364 : enddo
2365 0 : write(unt,'(a16,i2)')"Number of spin =",wan%nsppol
2366 0 : write(unt,'(a19,i4)')"Number of k-points=",wan%nkpt
2367 0 : do ikpt=1,wan%nkpt
2368 0 : write(unt,'(a,2x,i4)')"ikpt =",ikpt
2369 0 : do spin=1,wan%nsppol
2370 0 : do ispinor=1,wan%nspinor
2371 0 : do iband=wan%bandi_wan,wan%bandf_wan
2372 0 : ibandc=iband-wan%bandi_wan+1
2373 0 : write(unt,'(2x,a,2x,i2,2x,i2)')"iband =",iband
2374 0 : do iatom=1,wan%natom_wan
2375 0 : do il=1,wan%nbl_atom_wan(iatom)
2376 0 : do im=1,2*wan%latom_wan(iatom)%lcalc(il)+1
2377 0 : write(unt,'(8x,3i3,2x,2f23.15)')iatom,wan%latom_wan(iatom)%lcalc(il),im,&
2378 0 : &real(wan%psichi(ikpt,ibandc,iatom)%atom(il)%matl(im,spin,ispinor))&
2379 0 : &,aimag(wan%psichi(ikpt,ibandc,iatom)%atom(il)%matl(im,spin,ispinor))
2380 : enddo!m
2381 : enddo!l
2382 : enddo!atom
2383 : enddo!band
2384 : enddo!spinor
2385 : enddo!spin
2386 : enddo!k-point
2387 0 : close(unt)
2388 0 : end subroutine print_plowannier
2389 : !!***
2390 :
2391 :
2392 : !!****f* m_plowannier/get_plowannier
2393 : !! NAME
2394 : !! get_plowannier
2395 : !!
2396 : !! FUNCTION
2397 : !! get the psichies (Wannier weights) from a data.plowann file
2398 : !!
2399 : !! INPUTS
2400 : !! wan
2401 : !!
2402 : !! OUTPUT
2403 : !! wan
2404 : !!
2405 : !! SOURCE
2406 :
2407 0 : subroutine get_plowannier(wan_in,wan_out,dtset)
2408 :
2409 : use m_abicore
2410 : use defs_abitypes
2411 : use m_io_tools, only : open_file
2412 : use m_specialmsg, only : wrtout
2413 :
2414 : !Arguments-------------------------
2415 : type(plowannier_type),intent(inout) :: wan_in
2416 : type(plowannier_type),intent(inout) :: wan_out
2417 : type(dataset_type),intent(in) :: dtset
2418 : !Local variables-------------------
2419 : character(len=500) :: msg
2420 : integer :: unt,iatom,spin,ikpt,iband,ibandc,il,ispinor,im,dummy,natom,bandi,bandf,nbl,nspin,nkpt
2421 : integer :: t2g
2422 : real(dp) ::xx,yy
2423 :
2424 0 : t2g=dtset%dmft_t2g
2425 :
2426 :
2427 : !Opening of the data.plowann file
2428 0 : if (open_file('data.plowann',msg,newunit=unt,form='formatted',status='old') /= 0) then
2429 0 : ABI_ERROR(msg)
2430 : end if
2431 0 : rewind(unt)
2432 :
2433 :
2434 : !Reading of the header of data.plowann
2435 0 : read(unt,'(a22,i2)') msg, natom
2436 0 : read(unt,*)
2437 0 : read(unt,'(a7,2i4)') msg, bandi,bandf
2438 0 : read(unt,'(a26,i2)') msg, nbl
2439 0 : do iatom=1,wan_in%natom_wan
2440 0 : read(unt,*)
2441 : enddo
2442 0 : read(unt,'(a16,i2)') msg, nspin
2443 0 : read(unt,'(a19,i4)') msg, nkpt
2444 :
2445 : !Testing the header
2446 : if (natom /= wan_in%natom_wan .OR.&
2447 0 : & nbl/= sum(wan_in%nbl_atom_wan(:)) .OR. nspin /= wan_in%nsppol .OR. nkpt/=wan_in%nkpt ) then
2448 0 : write(msg,'(a,3i3)')"Not the same atoms or bands in both datasets",natom,bandi,bandf
2449 0 : ABI_ERROR(msg)
2450 : endif
2451 :
2452 : call init_plowannier(bandf,bandi,dtset%plowan_compute,&
2453 : &dtset%plowan_iatom,dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,&
2454 : &dtset%plowan_nbl,dtset%plowan_nt,dtset%plowan_projcalc,dtset%acell_orig,&
2455 0 : &dtset%kptns,sum(dtset%plowan_nbl),dtset%nimage,dtset%nkpt,dtset%nspinor,dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wan_out)
2456 :
2457 0 : call destroy_plowannier(wan_in)
2458 :
2459 0 : write(msg,'(a)')"Reading of the Wannier weights from data.plowann"
2460 0 : call wrtout(std_out,msg,'COLL')
2461 0 : call wrtout(ab_out,msg,'COLL')
2462 : !Reading of the psichis
2463 0 : do ikpt=1,wan_out%nkpt
2464 0 : read(unt,*)
2465 0 : do spin=1,wan_out%nsppol
2466 0 : do ispinor=1,wan_out%nspinor
2467 0 : do iband=wan_out%bandi_wan,wan_out%bandf_wan
2468 0 : ibandc=iband-wan_out%bandi_wan+1
2469 0 : read(unt,*)
2470 0 : do iatom=1,wan_out%natom_wan
2471 0 : do il=1,wan_out%nbl_atom_wan(iatom)
2472 0 : do im=1,2*wan_out%latom_wan(iatom)%lcalc(il)+1
2473 0 : read(unt,'(8x,3i3,2x,2f23.15)')dummy,dummy,dummy,xx,yy
2474 0 : wan_out%psichi(ikpt,ibandc,iatom)%atom(il)%matl(im,spin,ispinor)=cmplx(xx,yy)
2475 : enddo!m
2476 : enddo!l
2477 : enddo!atom
2478 : enddo!band
2479 : enddo!spinor
2480 : enddo!spin
2481 : enddo!k-point
2482 0 : close(unt)
2483 0 : end subroutine get_plowannier
2484 : !!***
2485 :
2486 :
2487 : !!****f* m_plowannier/fullbz_plowannier
2488 : !! NAME
2489 : !! fullbz_plowannier
2490 : !!
2491 : !! FUNCTION
2492 : !! Reconstruct the pischis on the full BZ
2493 : !!
2494 : !! INPUTS
2495 : !! dtset,kmesh,cryst,wanibz
2496 : !!
2497 : !! OUTPUT
2498 : !! wanbz
2499 : !!
2500 : !! SOURCE
2501 :
2502 0 : subroutine fullbz_plowannier(dtset,kmesh,cryst,pawang,wanibz,wanbz)
2503 :
2504 : use m_abicore
2505 : use m_specialmsg, only : wrtout
2506 : use defs_abitypes
2507 : use m_bz_mesh, only : kmesh_t
2508 : use m_crystal, only : crystal_t
2509 : use m_pawang, only : pawang_type
2510 :
2511 : !Arguments-------------------------
2512 : type(plowannier_type),intent(inout) :: wanibz
2513 : type(plowannier_type),intent(out) :: wanbz
2514 : type(dataset_type),intent(in) :: dtset
2515 : type(kmesh_t),intent(in) :: kmesh
2516 : type(crystal_t),intent(in) :: cryst
2517 : type(pawang_type),intent(in) :: pawang
2518 : !Local variables----------------------
2519 : character(len=500) :: msg
2520 : integer :: sym,iatom,spin,ik_bz,iband,ibandc,il,ispinor,im,ik_ibz,isym,itim
2521 : integer :: at_indx,indx, iat,m1,m2,l
2522 0 : real(dp) :: kbz(3),wtk(kmesh%nbz)
2523 : !*****************************************************************************************
2524 :
2525 0 : wtk=one
2526 0 : sym=kmesh%nbz/kmesh%nibz
2527 : call init_plowannier(wanibz%bandf_wan,wanibz%bandi_wan,dtset%plowan_compute,&
2528 : &dtset%plowan_iatom,dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,&
2529 : &dtset%plowan_nbl,dtset%plowan_nt,dtset%plowan_projcalc,dtset%acell_orig,&
2530 0 : &kmesh%bz,sum(dtset%plowan_nbl),dtset%nimage,kmesh%nbz,dtset%nspinor,dtset%nsppol,wtk,dtset%dmft_t2g,wanbz)
2531 :
2532 0 : write(msg,'(a)')" Reconstruction of the full Brillouin Zone using data.plowann in the IBZ"
2533 0 : call wrtout(std_out,msg,'COLL');call wrtout(ab_out,msg,'COLL')
2534 0 : if (cryst%nsym==1) then
2535 0 : do ik_bz=1,kmesh%nbz
2536 0 : do iband=wanbz%bandi_wan,wanbz%bandf_wan
2537 0 : ibandc=iband-wanbz%bandi_wan+1
2538 0 : do iatom=1,wanbz%natom_wan
2539 0 : do il=1,wanbz%nbl_atom_wan(iatom)
2540 0 : do im=1,2*wanbz%latom_wan(iatom)%lcalc(il)+1
2541 0 : do spin=1,wanbz%nsppol
2542 0 : do ispinor=1,wanbz%nspinor
2543 0 : if (kmesh%tabi(ik_bz)==1) then
2544 : wanbz%psichi(ik_bz,ibandc,iatom)%atom(il)%matl(im,spin,ispinor)=&
2545 0 : &wanibz%psichi(kmesh%tab(ik_bz),ibandc,iatom)%atom(il)%matl(im,spin,ispinor)
2546 0 : else if (kmesh%tabi(ik_bz)==-1) then
2547 : wanbz%psichi(ik_bz,ibandc,iatom)%atom(il)%matl(im,spin,ispinor)=&
2548 0 : &conjg(wanibz%psichi(kmesh%tab(ik_bz),ibandc,iatom)%atom(il)%matl(im,spin,ispinor))
2549 : endif
2550 : enddo
2551 : enddo
2552 : enddo
2553 : enddo
2554 : enddo
2555 : enddo
2556 : enddo
2557 0 : else if (cryst%nsym>1) then
2558 0 : do ik_bz=1,kmesh%nbz
2559 0 : call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym,itim)
2560 0 : do iatom=1,wanbz%natom_wan
2561 0 : indx=cryst%indsym(4,isym,wanibz%iatom_wan(iatom))
2562 : !Link beetween full list and wan list of atom
2563 0 : do iat=1,wanbz%natom_wan
2564 0 : if (indx==wanbz%iatom_wan(iat))then
2565 0 : at_indx=iat
2566 : end if
2567 : end do
2568 : !
2569 0 : do spin=1,wanibz%nsppol
2570 0 : do ispinor=1,wanibz%nspinor
2571 0 : do il=1,wanibz%nbl_atom_wan(iatom)
2572 0 : l=wanibz%latom_wan(iatom)%lcalc(il)
2573 0 : do m1=1,2*l+1
2574 0 : do m2=1,2*l+1
2575 0 : do iband=wanibz%bandi_wan,wanibz%bandf_wan
2576 0 : ibandc=iband-wanibz%bandi_wan+1
2577 : wanbz%psichi(ik_bz,ibandc,iatom)%atom(il)%matl(m1,spin,ispinor)=&
2578 : &wanbz%psichi(ik_bz,ibandc,iatom)%atom(il)%matl(m1,spin,ispinor)+&
2579 : &wanibz%psichi(ik_ibz,ibandc,at_indx)%atom(il)%matl(m2,spin,ispinor)&
2580 0 : &*pawang%zarot(m2,m1,l+1,isym)
2581 : end do
2582 : enddo
2583 : enddo
2584 : enddo
2585 : enddo
2586 : enddo
2587 : enddo
2588 : enddo
2589 : end if
2590 0 : call destroy_plowannier(wanibz)
2591 0 : end subroutine fullbz_plowannier
2592 : !!***
2593 :
2594 : !!****f* m_plowannier/destroy_plowannier
2595 : !! NAME
2596 : !! destroy_plowannier
2597 : !!
2598 : !! FUNCTION
2599 : !! deallocate variables
2600 : !!
2601 : !! INPUTS
2602 : !! wan
2603 : !!
2604 : !! OUTPUT
2605 : !!
2606 : !! SOURCE
2607 :
2608 :
2609 0 : subroutine destroy_plowannier(wan)
2610 :
2611 : !Arguments-------------------------------------
2612 : type(plowannier_type), intent(inout) :: wan
2613 : !Local variables-------------------------------
2614 : integer :: iatom,ikpt,iband,il
2615 :
2616 0 : do iatom=1,wan%natom_wan
2617 0 : ABI_FREE(wan%latom_wan(iatom)%lcalc)
2618 0 : ABI_FREE(wan%projector_wan(iatom)%lproj)
2619 0 : ABI_FREE(wan%nposition(iatom)%pos)
2620 : enddo
2621 0 : do iatom = 1,wan%natom_wan
2622 0 : do il = 1,wan%nbl_atom_wan(iatom)
2623 0 : ABI_FREE(wan%psichi(1,1,iatom)%atom(il)%ph0phiint)
2624 : end do
2625 : end do
2626 0 : do ikpt = 1,wan%nkpt
2627 0 : do iband = wan%bandi_wan,wan%bandf_wan
2628 0 : do iatom = 1,wan%natom_wan
2629 0 : do il = 1,wan%nbl_atom_wan(iatom)
2630 0 : ABI_FREE(wan%psichi(ikpt,iband-wan%bandi_wan+1,iatom)%atom(il)%matl)
2631 : end do
2632 0 : ABI_FREE(wan%psichi(ikpt,iband-wan%bandi_wan+1,iatom)%atom)
2633 : end do
2634 : end do
2635 : end do
2636 :
2637 0 : ABI_SFREE(wan%kpt)
2638 0 : ABI_SFREE(wan%iatom_wan)
2639 0 : ABI_SFREE(wan%nbl_atom_wan)
2640 0 : ABI_SFREE(wan%latom_wan)
2641 0 : ABI_SFREE(wan%nbproj_atom_wan)
2642 0 : ABI_SFREE(wan%projector_wan)
2643 0 : ABI_SFREE(wan%position)
2644 0 : ABI_SFREE(wan%wtk)
2645 0 : ABI_SFREE(wan%acell)
2646 0 : ABI_SFREE(wan%nposition)
2647 0 : ABI_SFREE(wan%psichi)
2648 :
2649 0 : end subroutine destroy_plowannier
2650 : !!***
2651 :
2652 :
2653 :
2654 : !!****f* m_plowannier/initialize_operwan
2655 : !! NAME
2656 : !! initialize_operwan
2657 : !!
2658 : !! FUNCTION
2659 : !! initialize operwan
2660 : !!
2661 : !! INPUTS
2662 : !! wan
2663 : !!
2664 : !! OUTPUT
2665 : !! operwan
2666 : !!
2667 : !! SOURCE
2668 :
2669 0 : subroutine initialize_operwan(wan,operwan)
2670 :
2671 : !Arguments----------------------------------
2672 : type(plowannier_type), intent(in) :: wan
2673 : type(operwan_type), intent(inout) :: operwan(wan%nkpt,wan%natom_wan,wan%natom_wan)
2674 :
2675 : !Local variables----------------------------
2676 : integer :: ikpt,iatom1,iatom2,il1,il2,n1,n2
2677 :
2678 0 : do ikpt = 1,wan%nkpt
2679 0 : do iatom1 = 1,wan%natom_wan
2680 0 : do iatom2 = 1,wan%natom_wan
2681 0 : ABI_MALLOC(operwan(ikpt,iatom1,iatom2)%atom,(wan%nbl_atom_wan(iatom1),wan%nbl_atom_wan(iatom2)))
2682 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2683 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2684 0 : n1=2*wan%latom_wan(iatom1)%lcalc(il1)+1
2685 0 : n2=2*wan%latom_wan(iatom2)%lcalc(il2)+1
2686 0 : ABI_MALLOC(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl,(n1,n2,wan%nsppol,wan%nspinor,wan%nspinor))
2687 0 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl = zero
2688 : end do
2689 : end do
2690 : end do
2691 : end do
2692 : end do
2693 :
2694 0 : end subroutine initialize_operwan
2695 : !!***
2696 :
2697 :
2698 : !!****f* m_plowannier/destroy_operwan
2699 : !! NAME
2700 : !! destroy_operwan
2701 : !!
2702 : !! FUNCTION
2703 : !! destroy operwan
2704 : !!
2705 : !! INPUTS
2706 : !! wan
2707 : !!
2708 : !! OUTPUT
2709 : !! operwan
2710 : !!
2711 : !! SOURCE
2712 :
2713 0 : subroutine destroy_operwan(wan,operwan)
2714 :
2715 : !Arguments----------------------------------
2716 : type(plowannier_type), intent(in) :: wan
2717 : type(operwan_type), intent(inout) :: operwan(wan%nkpt,wan%natom_wan,wan%natom_wan)
2718 :
2719 : !Local variables----------------------------
2720 : integer :: ikpt,iatom1,iatom2,il1,il2
2721 :
2722 :
2723 0 : do ikpt = 1,wan%nkpt
2724 0 : do iatom1 = 1,wan%natom_wan
2725 0 : do iatom2 = 1,wan%natom_wan
2726 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2727 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2728 0 : ABI_FREE(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl)
2729 : end do
2730 : end do
2731 0 : ABI_FREE(operwan(ikpt,iatom1,iatom2)%atom)
2732 : end do
2733 : end do
2734 : end do
2735 0 : end subroutine destroy_operwan
2736 : !!***
2737 :
2738 : !!****f* m_plowannier/zero_operwan
2739 : !! NAME
2740 : !! zero_operwan
2741 : !!
2742 : !! FUNCTION
2743 : !! zero operwan
2744 : !!
2745 : !! INPUTS
2746 : !! wan
2747 : !!
2748 : !! OUTPUT
2749 : !! operwan
2750 : !!
2751 : !! SOURCE
2752 :
2753 0 : subroutine zero_operwan(wan,operwan)
2754 :
2755 : !Arguments----------------------------------
2756 : type(plowannier_type), intent(in) :: wan
2757 : type(operwan_type), intent(inout) :: operwan(wan%nkpt,wan%natom_wan,wan%natom_wan)
2758 :
2759 : !Local variables----------------------------
2760 : integer :: ikpt, iatom1, iatom2, il1, il2, isppol, ispinor1, ispinor2, im1, im2
2761 :
2762 :
2763 0 : do ikpt = 1,wan%nkpt
2764 0 : do isppol = 1,wan%nsppol
2765 0 : do iatom1 = 1,wan%natom_wan
2766 0 : do iatom2 = 1,wan%natom_wan
2767 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2768 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2769 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
2770 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
2771 0 : do ispinor1 = 1,wan%nspinor
2772 0 : do ispinor2 = 1,wan%nspinor
2773 0 : operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)=czero
2774 : end do
2775 : end do
2776 : end do
2777 : end do
2778 : end do
2779 : end do
2780 : end do
2781 : end do
2782 : end do
2783 : end do
2784 :
2785 0 : end subroutine zero_operwan
2786 : !!***
2787 :
2788 : !!****f* m_plowannier/compute_oper_ks2wan
2789 : !! NAME
2790 : !! compute_oper_ks2wan
2791 : !!
2792 : !! FUNCTION
2793 : !! transform ks operator into wan one
2794 : !!
2795 : !! INPUTS
2796 : !! wan,operks,option
2797 : !!
2798 : !! OUTPUT
2799 : !! if option = ikpt, gives the wan operator in reciprocal space (for each k)
2800 : !!
2801 : !! SOURCE
2802 :
2803 0 : subroutine compute_oper_ks2wan(wan,operks,operwan,option)
2804 :
2805 : !Arguments--------------------------
2806 : type(plowannier_type), intent(in) :: wan
2807 : type(operwan_type), intent(inout) :: operwan(:,:,:)
2808 : complex(dp), intent(in) :: operks(:,:,:,:)
2809 : integer, intent(in) :: option
2810 :
2811 : !Local variables--------------------
2812 : integer :: iatom1, iatom2, il1, il2, isppol, ispinor1, ispinor2, iband1, iband2, im1, im2
2813 :
2814 : ! ----------------------------------
2815 : !Transformation KS2WAN
2816 : ! ----------------------------------
2817 :
2818 :
2819 : !!operation on reciprocal space
2820 0 : do iatom1 = 1,wan%natom_wan
2821 0 : do iatom2 = 1,wan%natom_wan
2822 0 : do isppol = 1,wan%nsppol
2823 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2824 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2825 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
2826 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
2827 0 : do ispinor1 = 1,wan%nspinor
2828 0 : do ispinor2 = 1,wan%nspinor
2829 : !!sum over the bands
2830 0 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
2831 0 : do iband2 = 1,wan%bandf_wan-wan%bandi_wan+1
2832 : operwan(option,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)=&
2833 : operwan(option,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)&
2834 : + conjg(wan%psichi(option,iband2,iatom2)%atom(il2)%matl(im2,isppol,ispinor2))&
2835 0 : *operks(option,iband1,iband2,isppol)*wan%psichi(option,iband1,iatom1)%atom(il1)%matl(im1,isppol,ispinor1)
2836 : end do
2837 : end do
2838 : ! write(6,*) "operwan",im1,im2,operwan(option,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)
2839 : end do
2840 : end do
2841 : end do
2842 : end do
2843 : end do
2844 : end do
2845 : end do
2846 : end do
2847 : end do
2848 :
2849 0 : end subroutine compute_oper_ks2wan
2850 : !!***
2851 :
2852 :
2853 :
2854 : !!****f* m_plowannier/normalization_plowannier
2855 : !! NAME
2856 : !! normalization_plowannier
2857 : !!
2858 : !! FUNCTION
2859 : !! Use compute_oper_ks2wan to calculate overlap and do the normalization for the wan%psichi coefficients
2860 : !!
2861 : !! INPUTS
2862 : !! wan, opt (=0 normalize k-point by k-point; =1 normalize the sum over k)
2863 : !!
2864 : !! OUTPUT
2865 : !! wan itself is modified
2866 : !!
2867 : !! SOURCE
2868 :
2869 :
2870 0 : subroutine normalization_plowannier(wan,opt)
2871 :
2872 : use m_matrix, only : invsqrt_matrix
2873 :
2874 : !Arguments------------------
2875 : type(plowannier_type), intent(inout) :: wan
2876 : integer, intent(in) :: opt
2877 : !Local----------------------
2878 : complex(dp), allocatable :: operks(:,:,:,:)
2879 0 : type(operwan_type), allocatable :: operwan(:,:,:)
2880 0 : complex(dp), allocatable :: operwansquare(:,:,:,:)
2881 0 : complex(dp), allocatable :: tmp_operwansquare(:,:)
2882 : integer :: ikpt, iband, iband1, iband2, isppol, ispinor1, ispinor2, iatom1,nb_zeros_tot
2883 : integer :: iatom2, il1, il2, im1, im2, index_c, index_l, n1,n2,n3, nkpt,nb_of_zeros
2884 0 : type(orbital_type), allocatable :: psichinormalized(:,:,:)
2885 : !character(len = 50) :: mat_writing2
2886 : !character(len = 5000) :: mat_writing
2887 : character(len = 500) :: message
2888 :
2889 : !Initialize nkpt (wan%nkpt if opt=0, 1 if opt=1)
2890 0 : if (opt==1) then
2891 : nkpt=1
2892 : else
2893 0 : nkpt=wan%nkpt
2894 : end if
2895 :
2896 : !First, creation of the ks identity operator
2897 0 : ABI_MALLOC(operks,(wan%nkpt,wan%bandf_wan-wan%bandi_wan+1,wan%bandf_wan-wan%bandi_wan+1,wan%nsppol))
2898 0 : operks = czero
2899 0 : do iband1 = 1,wan%bandf_wan-wan%bandi_wan+1
2900 0 : do iband2 = 1,wan%bandf_wan-wan%bandi_wan+1
2901 0 : if (iband1.eq.iband2) then
2902 0 : do ikpt = 1,wan%nkpt
2903 0 : do isppol= 1,wan%nsppol
2904 0 : operks(ikpt,iband1,iband2,isppol) = cone
2905 : end do
2906 : end do
2907 : end if
2908 : end do
2909 : end do
2910 :
2911 :
2912 : !Allocation of operwan
2913 0 : ABI_MALLOC(operwan,(wan%nkpt,wan%natom_wan,wan%natom_wan))
2914 0 : call initialize_operwan(wan,operwan)
2915 :
2916 :
2917 :
2918 : !Computation of the overlap
2919 0 : do ikpt = 1,wan%nkpt
2920 0 : call compute_oper_ks2wan(wan,operks,operwan,ikpt)
2921 : end do
2922 :
2923 :
2924 :
2925 : !transform the operwan into an inversible matrix
2926 : !!operwansquare is the overlap square matrix (wan%size_wan * wan%size_wan)
2927 0 : ABI_MALLOC(operwansquare,(wan%nkpt,wan%nsppol,wan%nspinor*wan%size_wan,wan%nspinor*wan%size_wan))
2928 :
2929 0 : operwansquare = czero
2930 :
2931 0 : n1=size(wan%psichi,1)
2932 0 : n2=size(wan%psichi,2)
2933 0 : n3=size(wan%psichi,3)
2934 0 : ABI_MALLOC(psichinormalized,(n1,n2,n3))
2935 0 : call allocate_orbital(wan%psichi,psichinormalized,n1,n2,n3)
2936 0 : call copy_orbital(wan%psichi,psichinormalized,n1,n2,n3)
2937 :
2938 0 : do isppol = 1,wan%nsppol
2939 0 : do ispinor1 = 1,wan%nspinor
2940 0 : do ispinor2 = 1,wan%nspinor
2941 0 : index_l = 0 ! index_l is set to 0 at the beginning
2942 0 : do iatom1 = 1,wan%natom_wan
2943 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
2944 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
2945 0 : index_l = index_l + 1 ! the line changes
2946 0 : index_c = 1 ! counter_c is set to one each time the line changes
2947 0 : do iatom2 = 1,wan%natom_wan
2948 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
2949 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
2950 0 : do ikpt = 1,wan%nkpt
2951 0 : if (opt==1) then ! operwansquare is a sum of operwan over k points
2952 : operwansquare(1,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+wan%size_wan*(ispinor2-1)) &
2953 : & =operwansquare(1,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+wan%size_wan*(ispinor2-1))+ &
2954 0 : & (operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)*wan%wtk(ikpt))
2955 : else !operwansquare is a matrix with a k-point dimension
2956 : operwansquare(ikpt,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+wan%size_wan*(ispinor2-1)) &
2957 0 : &= operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)
2958 : end if
2959 : end do
2960 0 : index_c = index_c + 1
2961 : end do !im2
2962 : end do !il2
2963 : end do ! iatom2 (the line changes)
2964 : end do ! im1
2965 : end do ! il1
2966 : end do !iatom1
2967 : end do
2968 : end do
2969 : end do
2970 :
2971 :
2972 :
2973 :
2974 : !!Write the overlap matrix for ikpt = 1 in a nice shape
2975 : ! do isppol = 1,wan%nsppol
2976 : ! do il1 = 1,size(operwansquare,3) !! dummy variable without any meaning
2977 : ! write(mat_writing,'(a,i0,i0)') 'Overlap matrix before orthonormalization 1 ',isppol,il1
2978 : ! do il2 = 1,size(operwansquare,4)
2979 : ! write(mat_writing2,'(F10.6)') real(operwansquare(1,isppol,il1,il2))
2980 : ! mat_writing = trim(mat_writing)//trim(mat_writing2)
2981 : ! end do
2982 : ! print*,trim(mat_writing)
2983 : ! end do
2984 : ! end do
2985 :
2986 :
2987 :
2988 : !take the square root inverse of operwansquare for normalization purposes
2989 0 : nb_zeros_tot=0
2990 0 : ABI_MALLOC(tmp_operwansquare,(wan%nspinor*wan%size_wan,wan%nspinor*wan%size_wan))
2991 0 : do isppol = 1,wan%nsppol
2992 0 : do ikpt = 1,nkpt
2993 0 : write(std_out,*)"ikpt = ", ikpt
2994 0 : tmp_operwansquare(:,:)=operwansquare(ikpt,isppol,:,:)
2995 0 : call invsqrt_matrix(tmp_operwansquare,wan%nspinor*wan%size_wan,nb_of_zeros)
2996 0 : operwansquare(ikpt,isppol,:,:)=tmp_operwansquare(:,:)
2997 0 : nb_zeros_tot=nb_zeros_tot+nb_of_zeros
2998 : end do
2999 : end do
3000 0 : ABI_FREE(tmp_operwansquare)
3001 :
3002 0 : do ikpt = 1,wan%nkpt
3003 0 : do iband = 1,wan%bandf_wan-wan%bandi_wan+1
3004 0 : do iatom1 = 1,wan%natom_wan
3005 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3006 0 : psichinormalized(ikpt,iband,iatom1)%atom(il1)%matl = czero
3007 : end do
3008 : end do
3009 : end do
3010 : end do
3011 :
3012 :
3013 :
3014 : ! compute the new psichi normalized
3015 0 : do isppol = 1,wan%nsppol
3016 0 : do ispinor1 = 1,wan%nspinor
3017 0 : do iband = 1,wan%bandf_wan-wan%bandi_wan+1
3018 0 : index_l = 0
3019 0 : do iatom1 = 1,wan%natom_wan
3020 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3021 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3022 : ! sum
3023 0 : do ispinor2 = 1,wan%nspinor
3024 0 : index_l = index_l + 1 ! the line changes
3025 0 : index_c = 1 ! when the line changes, index_c is set to 1
3026 0 : do iatom2 = 1,wan%natom_wan
3027 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3028 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3029 0 : do ikpt = 1,wan%nkpt
3030 0 : if (opt==1)then ! all the psichi are normalized with the same operwansquare
3031 : psichinormalized(ikpt,iband,iatom1)%atom(il1)%matl(im1,isppol,ispinor1) =&
3032 : & psichinormalized(ikpt,iband,iatom1)%atom(il1)%matl(im1,isppol,ispinor1) +&
3033 : & wan%psichi(ikpt,iband,iatom2)%atom(il2)%matl(im2,isppol,ispinor2)*&
3034 : & operwansquare(1,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+&
3035 0 : & wan%size_wan*(ispinor2-1))
3036 : else ! each psichi is normalized with his own operwansquare
3037 : psichinormalized(ikpt,iband,iatom1)%atom(il1)%matl(im1,isppol,ispinor1) =&
3038 : & psichinormalized(ikpt,iband,iatom1)%atom(il1)%matl(im1,isppol,ispinor1) +&
3039 : & wan%psichi(ikpt,iband,iatom2)%atom(il2)%matl(im2,isppol,ispinor2)*&
3040 : & operwansquare(ikpt,isppol,index_l+wan%size_wan*(ispinor1-1),index_c+&
3041 0 : & wan%size_wan*(ispinor2-1))
3042 : end if
3043 : end do
3044 0 : index_c = index_c + 1
3045 : end do !im2
3046 : end do !il2
3047 : end do !iatom2
3048 : end do ! ispinor2
3049 : end do !im1
3050 : end do ! il1
3051 : end do ! iatom1
3052 : end do ! iband
3053 : end do ! ispinor1
3054 : end do !isppol
3055 : ! copy the new psichi normalized
3056 0 : call copy_orbital(psichinormalized,wan%psichi,n1,n2,n3)
3057 0 : call destroy_orbital(psichinormalized,n1,n2,n3)
3058 0 : ABI_FREE(psichinormalized)
3059 :
3060 0 : call destroy_operwan(wan,operwan)
3061 0 : ABI_FREE(operwan)
3062 :
3063 :
3064 :
3065 :
3066 :
3067 :
3068 :
3069 :
3070 : !!
3071 : !! !-------------------------------------------------------------
3072 : !! !check if the new norm is one
3073 :
3074 0 : ABI_MALLOC(operwan,(wan%nkpt,wan%natom_wan,wan%natom_wan))
3075 0 : call initialize_operwan(wan,operwan)
3076 0 : do ikpt = 1,wan%nkpt
3077 0 : call compute_oper_ks2wan(wan,operks,operwan,ikpt)
3078 : end do
3079 :
3080 :
3081 0 : do isppol = 1,wan%nsppol
3082 0 : do ikpt = 1,wan%nkpt
3083 0 : do iatom1 = 1,wan%natom_wan
3084 0 : do iatom2 = 1,wan%natom_wan
3085 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3086 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3087 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3088 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3089 0 : do ispinor1 = 1,wan%nspinor
3090 0 : do ispinor2 = 1,wan%nspinor
3091 0 : if (opt==0 .and. nb_zeros_tot==0) then
3092 0 : if (iatom1.eq.iatom2 .and. il1.eq.il2 .and. im1.eq.im2 .and. ispinor1.eq.ispinor2) then
3093 0 : if (abs(cmplx(1.0,0.0,dp)-&
3094 : &operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%&
3095 : &matl(im1,im2,isppol,ispinor1,ispinor2)) > 1d-8) then
3096 0 : write(message,'(a,i0,a,F18.11)') 'Normalization error for ikpt =',ikpt,&
3097 0 : &' on diag, value = ',&
3098 0 : &abs(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2))
3099 0 : ABI_ERROR(message)
3100 : end if
3101 : else
3102 0 : if (abs(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)) > 1d-8) then
3103 0 : write(message,'(a,i0,a,F10.3)') 'Normalization error for ikpt =',ikpt,&
3104 0 : &' not on diag, value = ',&
3105 0 : &abs(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2))
3106 0 : ABI_ERROR(message)
3107 : end if
3108 : end if
3109 : end if
3110 : end do
3111 : end do
3112 : end do
3113 : end do
3114 : end do
3115 : end do
3116 : end do
3117 : end do
3118 : end do
3119 : end do
3120 0 : if (opt==0 .and. nb_zeros_tot/=0) then
3121 0 : write(message,'(a,i2,a)')"The matrix inversion detects ",nb_zeros_tot,&
3122 0 : " zero(s) on the diagonals. Take results with caution or modify nkpt and/or bands for plowan"
3123 0 : ABI_COMMENT(message)
3124 : end if
3125 :
3126 : !!Uncomment to print the overlap matrix in the log file (for ikpt = 1)
3127 : ! do isppol = 1,wan%nsppol
3128 : ! do ispinor1 = 1,wan%nspinor
3129 : ! do ispinor2 = 1,wan%nspinor
3130 : ! index_l = 0 ! index_l is set to 0 at the beginning
3131 : ! do iatom1 = 1,wan%natom_wan
3132 : ! do il1 = 1,wan%nbl_atom_wan(iatom1)
3133 : ! do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3134 : ! index_l = index_l + 1 ! the line changes
3135 : ! index_c = 1 ! counter_c is set to one each time the line changes
3136 : ! do iatom2 = 1,wan%natom_wan
3137 : ! do il2 = 1,wan%nbl_atom_wan(iatom2)
3138 : ! do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3139 : ! do ikpt = 1,wan%nkpt
3140 : ! operwansquare(ikpt,isppol,index_l+wan%size_wan*(ispinor1-1),&
3141 : ! &index_c+wan%size_wan*(ispinor2-1)) = operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)
3142 : ! end do
3143 : ! index_c = index_c + 1
3144 : ! end do !im2
3145 : ! end do !il2
3146 : ! end do ! iatom2 (the line changes)
3147 : ! end do ! im1
3148 : ! end do ! il1
3149 : ! end do !iatom1
3150 : ! end do
3151 : ! end do
3152 : ! end do
3153 :
3154 : !!Print the overlap matrix in a nice shape
3155 : ! do isppol = 1,wan%nsppol
3156 : ! do il1 = 1,size(operwansquare,3) !! dummy variable without any meaning
3157 : ! write(mat_writing,'(a,i0,i0)') 'Overlap matrix after orthonormalization ',isppol,il1
3158 : ! do il2 = 1,size(operwansquare,4)
3159 : ! write(mat_writing2,'(F10.6)') real(operwansquare(9,isppol,il1,il2))
3160 : ! mat_writing = trim(mat_writing)//trim(mat_writing2)
3161 : ! end do
3162 : ! print*,trim(mat_writing)
3163 : ! end do
3164 : ! end do
3165 :
3166 :
3167 :
3168 :
3169 : !! !----------------------------------------------------------------
3170 0 : ABI_FREE(operwansquare)
3171 0 : ABI_FREE(operks)
3172 0 : call destroy_operwan(wan,operwan)
3173 0 : ABI_FREE(operwan)
3174 :
3175 :
3176 0 : end subroutine normalization_plowannier
3177 :
3178 : !!***
3179 :
3180 :
3181 :
3182 : !!****f* m_plowannier/print_operwan
3183 : !! NAME
3184 : !! print_operwan
3185 : !!
3186 : !! FUNCTION
3187 : !! Print the Wannier operator (real space) in a latex file
3188 : !!
3189 : !! INPUTS
3190 : !! wan, operwan, name
3191 : !!
3192 : !! OUTPUT
3193 : !!
3194 : !!
3195 : !! SOURCE
3196 :
3197 :
3198 0 : subroutine print_operwan(wan,operwan,name,convert)
3199 :
3200 : !Arguments----------------------------------
3201 : type(operwan_type),intent(in) :: operwan(:,:,:)
3202 : type(plowannier_type), intent(in) :: wan
3203 : character(len=*), intent(in) :: name
3204 : real(dp), intent(in) :: convert
3205 :
3206 : !Local variables----------------------------
3207 : integer :: iatom1,iatom2,pos1,pos2,il1,il2,im1,im2,isppol,ikpt,unt
3208 : real(dp) :: sum
3209 : character(len = 500) :: str1,str2,msg
3210 :
3211 0 : if (open_file(name, msg, newunit=unt) /= 0) then
3212 0 : ABI_ERROR(msg)
3213 : end if
3214 :
3215 0 : write(unt,'(a)') '\documentclass[11pt,a4paper,landscape]{article}'
3216 0 : write(unt,'(a)') '\usepackage[T1]{fontenc}'
3217 0 : write(unt,'(a)') '\usepackage{geometry,tabularx,graphicx}'
3218 0 : write(unt,'(a)') '\geometry{left=0.5cm,right=0.5cm}'
3219 :
3220 0 : write(unt,'(a)') '\begin{document}'
3221 0 : write(unt,'(a)') '\noindent'
3222 :
3223 : ! write(unt,'(a,i0,a,F7.3,a)') "% ",wan%natom_wan," atom in a ",wan%acell(1)," cell"
3224 : ! write(unt,'(a)') "% atom isppol proj"
3225 :
3226 :
3227 0 : do isppol = 1,wan%nsppol
3228 0 : write(unt,'(a)') '\begin{figure}'
3229 0 : write(unt,'(a)') '\resizebox{\linewidth}{!}{%'
3230 0 : write(unt,'(a)') '$ \left('
3231 :
3232 :
3233 0 : write(str1,'(a)') '\begin{array}{'
3234 0 : do iatom1 = 1,wan%natom_wan
3235 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
3236 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3237 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3238 0 : str1 = trim(str1)//"c"
3239 : end do
3240 0 : if (iatom1 .ne. wan%natom_wan .or. il1 .ne. wan%nbl_atom_wan(iatom1) .or. pos1 .ne. size(wan%nposition(iatom1)%pos,1)) then
3241 0 : str1 = trim(str1)//'|'
3242 : end if
3243 : end do
3244 : end do
3245 : end do
3246 0 : str1 = trim(str1)//'}'
3247 0 : write(unt,'(a)') str1
3248 :
3249 :
3250 0 : do iatom1 = 1,wan%natom_wan
3251 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
3252 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3253 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3254 0 : write(str1,'(a)') ""
3255 0 : do iatom2 = 1,wan%natom_wan
3256 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
3257 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3258 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3259 0 : sum = 0
3260 0 : do ikpt = 1,wan%nkpt
3261 : sum = sum + convert*real(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,1,1)*wan%wtk(ikpt)*&
3262 : & exp(cmplx(0.0,1.0)*two_pi*(wan%kpt(1,ikpt)*( &
3263 : & wan%nposition(iatom1)%pos(pos1,1)-wan%nposition(iatom2)%pos(pos2,1))+ &
3264 : & wan%kpt(2,ikpt)*(wan%nposition(iatom1)%pos(pos1,2)-wan%nposition(iatom2)%pos(pos2,2))+ &
3265 0 : & wan%kpt(3,ikpt)*(wan%nposition(iatom1)%pos(pos1,3)-wan%nposition(iatom2)%pos(pos2,3)))))
3266 : end do
3267 0 : write(str2,'(F10.6)') real(sum)
3268 0 : if ( len_trim(str1) .ge. 2) then
3269 0 : str1 = trim(str1)//"&"//trim(str2)
3270 : else
3271 0 : str1 = trim(str2)
3272 : end if
3273 : if (iatom2 .eq. wan%natom_wan .and. il2 .eq. wan%nbl_atom_wan(iatom2) .and. im2 &
3274 0 : & .eq. 2*wan%latom_wan(iatom2)%lcalc(il2)+1 .and. pos2 .eq. size(wan%nposition(iatom2)%pos,1)) then
3275 0 : str1 = trim(str1)//'\\'
3276 : end if
3277 : end do
3278 : end do
3279 : end do
3280 : end do
3281 0 : write(unt,'(a)') trim(str1)
3282 : end do
3283 0 : if (iatom1 .ne. wan%natom_wan .or. il1 .ne. wan%nbl_atom_wan(iatom1) .or. pos1 .ne. size(wan%nposition(iatom1)%pos,1)) then
3284 0 : write(unt,'(a)') '\hline'
3285 : end if
3286 : end do
3287 : end do
3288 : end do
3289 :
3290 :
3291 :
3292 0 : write(unt,'(a)') '\end{array} \right) $ }'
3293 :
3294 0 : if (name(len_trim(name)-2:len(trim(name))) == 'gen') then
3295 0 : write(unt,'(a,i0,a,F7.3,a)') "\caption{Energy matrix in real space for isppol = ", &
3296 0 : & isppol," in a ",wan%acell(1), " a.u. cell}"
3297 : end if
3298 :
3299 0 : if (name(len_trim(name)-2:len(trim(name))) == 'occ') then
3300 0 : write(unt,'(a,i0,a,F7.3,a)') "\caption{Occupation matrix in real space for isppol = ",&
3301 0 : & isppol," in a ",wan%acell(1), " a.u. cell}"
3302 : end if
3303 :
3304 :
3305 :
3306 0 : write(unt,'(a)') '\end{figure}'
3307 0 : write(unt,'(a)') '\end{document}'
3308 :
3309 : end do
3310 0 : close(unt)
3311 :
3312 0 : end subroutine print_operwan
3313 : !!***
3314 :
3315 :
3316 : !!****f* m_plowannier/init_operwan_realspace
3317 : !! NAME
3318 : !! init_operwan_realspace
3319 : !!
3320 : !! FUNCTION
3321 : !! Initialize an operwan_realspace type variable
3322 : !!
3323 : !! INPUTS
3324 : !! wan, operwan_realspace
3325 : !!
3326 : !! OUTPUT
3327 : !! operwan_realspace
3328 : !!
3329 : !! SOURCE
3330 0 : subroutine init_operwan_realspace(wan,oprs)
3331 :
3332 : !Arguments----------------------------------
3333 : type(operwan_realspace_type),intent(inout) :: oprs
3334 : type(plowannier_type), intent(in) :: wan
3335 :
3336 : !Local variables----------------------------
3337 : integer :: i1,i2,n1,n2,p1,p2,l1,l2,sp,pi
3338 :
3339 : !variable names is shorten to achieve not too long line lenght
3340 0 : sp=wan%nsppol
3341 0 : pi=wan%nspinor
3342 0 : ABI_MALLOC(oprs%atom_index,(wan%natom_wan,wan%natom_wan))
3343 0 : do i1 = 1,wan%natom_wan
3344 0 : do i2 = 1,wan%natom_wan
3345 0 : n1=size(wan%nposition(i1)%pos,1)
3346 0 : n2=size(wan%nposition(i2)%pos,1)
3347 0 : ABI_MALLOC(oprs%atom_index(i1,i2)%position,(n1,n2))
3348 0 : do p1 = 1,size(wan%nposition(i1)%pos,1)
3349 0 : do p2 = 1,size(wan%nposition(i2)%pos,1)
3350 0 : n1=wan%nbl_atom_wan(i1)
3351 0 : n2=wan%nbl_atom_wan(i2)
3352 0 : ABI_MALLOC(oprs%atom_index(i1,i2)%position(p1,p2)%atom,(n1,n2))
3353 0 : do l1 = 1,wan%nbl_atom_wan(i1)
3354 0 : do l2 = 1,wan%nbl_atom_wan(i2)
3355 0 : n1=2*wan%latom_wan(i1)%lcalc(l1)+1
3356 0 : n2=2*wan%latom_wan(i2)%lcalc(l2)+1
3357 0 : ABI_MALLOC(oprs%atom_index(i1,i2)%position(p1,p2)%atom(l1,l2)%matl,(n1,n2,sp,pi,pi))
3358 0 : oprs%atom_index(i1,i2)%position(p1,p2)%atom(l1,l2)%matl = czero
3359 : end do
3360 : end do
3361 : end do
3362 : end do
3363 : end do
3364 : end do
3365 :
3366 0 : end subroutine init_operwan_realspace
3367 : !!***
3368 :
3369 : !!****f* m_plowannier/reduce_operwan_realspace
3370 : !! NAME
3371 : !! reduce_operwan_realspace
3372 : !!
3373 : !! FUNCTION
3374 : !! reduce a table of operwan_realspace type
3375 : !!
3376 : !! INPUTS
3377 : !! wan,rhot1,npwx,nibz,comm,nbz,nsppol
3378 : !!
3379 : !! OUTPUT
3380 : !! rhot1
3381 : !!
3382 : !! SOURCE
3383 :
3384 0 : subroutine reduce_operwan_realspace(wan,rhot1,npwx,nibz,comm,nbz,nsppol)
3385 :
3386 :
3387 : use m_xmpi, only : xmpi_barrier,xmpi_sum
3388 : !Arguments---------------------------------------
3389 : type(plowannier_type),intent(in) :: wan
3390 : integer, intent(in) :: npwx,nibz,comm,nbz,nsppol
3391 : type(operwan_realspace_type),target,intent(inout) :: rhot1(npwx,nibz)
3392 : !Local variables----------------------------------
3393 : complex(dp),allocatable :: buffer(:)
3394 : integer :: dim,pwx,ibz, spin, ispinor1, ispinor2, iatom1, iatom2, pos1, pos2
3395 : integer :: il1, il2, im1, im2, nnn, ierr
3396 0 : complex(dp),pointer :: oper_ptr(:,:,:,:,:)
3397 :
3398 :
3399 0 : dim=0
3400 0 : do pwx=1,npwx
3401 0 : do ibz=1,nibz
3402 0 : do spin=1,wan%nsppol
3403 0 : do ispinor1=1,wan%nspinor
3404 0 : do ispinor2=1,wan%nspinor
3405 0 : do iatom1=1,wan%natom_wan
3406 0 : do iatom2=1,wan%natom_wan
3407 0 : do pos1=1,size(wan%nposition(iatom1)%pos,1)
3408 0 : do pos2=1,size(wan%nposition(iatom2)%pos,1)
3409 0 : do il1=1,wan%nbl_atom_wan(iatom1)
3410 0 : do il2=1,wan%nbl_atom_wan(iatom2)
3411 0 : do im1=1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3412 0 : do im2=1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3413 0 : dim=dim+1
3414 : enddo!im2
3415 : enddo!im1
3416 : enddo!il2
3417 : enddo!il1
3418 : enddo!pos2
3419 : enddo!pos1
3420 : enddo!iatom2
3421 : enddo!iatom1
3422 : enddo!ispinor2
3423 : enddo!ispinor1
3424 : enddo!spin
3425 : enddo!ibz
3426 : enddo!pwx
3427 0 : ABI_MALLOC(buffer,(dim))
3428 0 : nnn=0
3429 0 : do pwx=1,npwx
3430 0 : do ibz=1,nibz
3431 0 : do iatom1=1,wan%natom_wan
3432 0 : do iatom2=1,wan%natom_wan
3433 0 : do pos1=1,size(wan%nposition(iatom1)%pos,1)
3434 0 : do pos2=1,size(wan%nposition(iatom2)%pos,1)
3435 0 : do il1=1,wan%nbl_atom_wan(iatom1)
3436 0 : do il2=1,wan%nbl_atom_wan(iatom2)
3437 0 : oper_ptr=>rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl
3438 0 : do im1=1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3439 0 : do im2=1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3440 0 : do spin=1,wan%nsppol
3441 0 : do ispinor1=1,wan%nspinor
3442 0 : do ispinor2=1,wan%nspinor
3443 0 : nnn=nnn+1
3444 0 : buffer(nnn)=oper_ptr(im1,im2,spin,ispinor1,ispinor2)
3445 : enddo!ispinor2
3446 : enddo!ispinor1
3447 : enddo!spin
3448 : enddo!im2
3449 : enddo!im1
3450 : enddo!il2
3451 : enddo!il1
3452 : enddo!pos2
3453 : enddo!pos1
3454 : enddo!iatom2
3455 : enddo!iatom1
3456 : enddo!ibz
3457 : enddo!pwx
3458 0 : call xmpi_barrier(comm)
3459 0 : call xmpi_sum(buffer,comm,ierr)
3460 0 : call xmpi_barrier(comm)
3461 0 : buffer=buffer/nbz/nsppol
3462 : nnn=0
3463 0 : do pwx=1,npwx
3464 0 : do ibz=1,nibz
3465 0 : do iatom1=1,wan%natom_wan
3466 0 : do iatom2=1,wan%natom_wan
3467 0 : do pos1=1,size(wan%nposition(iatom1)%pos,1)
3468 0 : do pos2=1,size(wan%nposition(iatom2)%pos,1)
3469 0 : do il1=1,wan%nbl_atom_wan(iatom1)
3470 0 : do il2=1,wan%nbl_atom_wan(iatom2)
3471 0 : oper_ptr=>rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl
3472 0 : do im1=1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3473 0 : do im2=1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3474 0 : do spin=1,wan%nsppol
3475 0 : do ispinor1=1,wan%nspinor
3476 0 : do ispinor2=1,wan%nspinor
3477 0 : nnn=nnn+1
3478 0 : oper_ptr(im1,im2,spin,ispinor1,ispinor2)=buffer(nnn)
3479 : enddo!im2
3480 : enddo!im1
3481 : enddo!il2
3482 : enddo!il1
3483 : enddo!pos2
3484 : enddo!pos1
3485 : enddo!iatom2
3486 : enddo!iatom1
3487 : enddo!ispinor2
3488 : enddo!ispinor1
3489 : enddo!spin
3490 : enddo!ibz
3491 : enddo!pwx
3492 0 : ABI_FREE(buffer)
3493 :
3494 :
3495 0 : end subroutine reduce_operwan_realspace
3496 : !!***
3497 :
3498 : !!****f* m_plowannier/destroy_operwan_realspace
3499 : !! NAME
3500 : !! destroy_operwan_realspace
3501 : !!
3502 : !! FUNCTION
3503 : !! Destroy an operwan_realspace type variable
3504 : !!
3505 : !! INPUTS
3506 : !! wan, operwan_realspace
3507 : !!
3508 : !! OUTPUT
3509 : !! operwan_realspace
3510 : !!
3511 : !! SOURCE
3512 0 : subroutine destroy_operwan_realspace(wan,operwan_realspace)
3513 :
3514 : !Arguments----------------------------------
3515 : type(operwan_realspace_type),intent(inout) :: operwan_realspace
3516 : type(plowannier_type), intent(in) :: wan
3517 :
3518 : !Local variables----------------------------
3519 : integer :: iatom1,iatom2,pos1,pos2,il1,il2
3520 :
3521 :
3522 0 : do iatom1 = 1,wan%natom_wan
3523 0 : do iatom2 = 1,wan%natom_wan
3524 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
3525 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
3526 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3527 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3528 0 : ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl)
3529 : end do
3530 : end do
3531 0 : ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom)
3532 : end do
3533 : end do
3534 0 : ABI_FREE(operwan_realspace%atom_index(iatom1,iatom2)%position)
3535 : end do
3536 : end do
3537 0 : ABI_FREE(operwan_realspace%atom_index)
3538 :
3539 0 : end subroutine destroy_operwan_realspace
3540 : !!***
3541 :
3542 :
3543 : !!****f* m_plowannier/zero_operwan_realspace
3544 : !! NAME
3545 : !! zero_operwan_realspace
3546 : !!
3547 : !! FUNCTION
3548 : !! Set an operwan_realspace to zero
3549 : !!
3550 : !! INPUTS
3551 : !! wan, operwan_realspace
3552 : !!
3553 : !! OUTPUT
3554 : !! operwan_realspace
3555 : !!
3556 : !! SOURCE
3557 0 : subroutine zero_operwan_realspace(wan,operwan_realspace)
3558 :
3559 : !Arguments----------------------------------
3560 : type(operwan_realspace_type),intent(inout) :: operwan_realspace
3561 : type(plowannier_type), intent(in) :: wan
3562 :
3563 : !Local variables----------------------------
3564 : integer :: isppol,iatom1,iatom2,pos1,pos2,il1,il2,im1,im2,ispinor1,ispinor2
3565 :
3566 :
3567 0 : do isppol = 1,wan%nsppol
3568 0 : do ispinor1=1,wan%nspinor
3569 0 : do ispinor2=1,wan%nspinor
3570 0 : do iatom1 = 1,wan%natom_wan
3571 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
3572 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3573 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3574 0 : do iatom2 = 1,wan%natom_wan
3575 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
3576 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3577 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3578 : operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%&
3579 0 : &atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)=czero
3580 : enddo
3581 : enddo
3582 : enddo
3583 : enddo
3584 : enddo
3585 : enddo
3586 : enddo
3587 : enddo
3588 : enddo
3589 : enddo
3590 : enddo
3591 0 : end subroutine zero_operwan_realspace
3592 : !!***
3593 :
3594 :
3595 :
3596 :
3597 : !!****f* m_plowannier/compute_oper_wank2realspace
3598 : !! NAME
3599 : !! compute_operwan_wanK2realspace
3600 : !!
3601 : !! FUNCTION
3602 : !! Compute an operator from WannierK space to real space
3603 : !!
3604 : !! INPUTS
3605 : !! wan,operwan, operwan_realspace
3606 : !!
3607 : !! OUTPUT
3608 : !! operwan_realspace
3609 : !!
3610 : !! SOURCE
3611 0 : subroutine compute_oper_wank2realspace(wan,operwan,operwan_realspace)
3612 :
3613 : !Arguments----------------------------------
3614 : type(operwan_realspace_type),intent(inout) :: operwan_realspace
3615 : type(operwan_type),intent(in) :: operwan(:,:,:)
3616 : type(plowannier_type), intent(in) :: wan
3617 :
3618 :
3619 : !Local variables----------------------------
3620 : integer :: isppol,iatom1,pos1,il1,im1,iatom2,pos2,il2,im2,ikpt,ispinor1,ispinor2
3621 :
3622 :
3623 0 : do isppol = 1,wan%nsppol
3624 0 : do ispinor1=1,wan%nspinor
3625 0 : do ispinor2=1,wan%nspinor
3626 0 : do iatom1 = 1,wan%natom_wan
3627 0 : do pos1 = 1,size(wan%nposition(iatom1)%pos,1)
3628 0 : do il1 = 1,wan%nbl_atom_wan(iatom1)
3629 0 : do im1 = 1,2*wan%latom_wan(iatom1)%lcalc(il1)+1
3630 0 : do iatom2 = 1,wan%natom_wan
3631 0 : do pos2 = 1,size(wan%nposition(iatom2)%pos,1)
3632 0 : do il2 = 1,wan%nbl_atom_wan(iatom2)
3633 0 : do im2 = 1,2*wan%latom_wan(iatom2)%lcalc(il2)+1
3634 : !sum over ikpt
3635 0 : do ikpt = 1,wan%nkpt
3636 : operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%&
3637 : &atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2) =&
3638 : operwan_realspace%atom_index(iatom1,iatom2)%position(pos1,pos2)%&
3639 : &atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)&
3640 : + real(operwan(ikpt,iatom1,iatom2)%atom(il1,il2)%matl(im1,im2,isppol,ispinor1,ispinor2)&
3641 : * wan%wtk(ikpt) * exp( cmplx(0.0,1.0) * two_pi * ( &
3642 : wan%kpt(1,ikpt) * ( wan%nposition(iatom1)%pos(pos1,1) - wan%nposition(iatom2)%pos(pos2,1) )+&
3643 : wan%kpt(2,ikpt) * ( wan%nposition(iatom1)%pos(pos1,2) - wan%nposition(iatom2)%pos(pos2,2) )+&
3644 0 : wan%kpt(3,ikpt) * ( wan%nposition(iatom1)%pos(pos1,3) - wan%nposition(iatom2)%pos(pos2,3)))))
3645 : end do
3646 : !end of the sum
3647 : end do
3648 : enddo
3649 : enddo
3650 : end do
3651 : end do
3652 : end do
3653 : end do
3654 : end do
3655 : end do
3656 : end do
3657 : end do
3658 0 : end subroutine compute_oper_wank2realspace
3659 : !!***
3660 0 : END MODULE m_plowannier
3661 : !!***
|
