Line data Source code
1 : !!****m* ABINIT/m_mlwfovlp
2 : !! NAME
3 : !! m_mlwfovlp
4 : !!
5 : !! FUNCTION
6 : !! Interface with Wannier90
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2005-2026 ABINIT group (BAmadon, CEspejo, FJollet, TRangel, DRH, hexu, MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_mlwfovlp
23 :
24 : use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use defs_wannier90
27 : use m_abicore
28 : use m_errors
29 : use m_atomdata
30 : use m_xmpi
31 : use m_sort
32 : #ifdef FC_NAG
33 : use f90_unix_dir
34 : #endif
35 : use netcdf
36 : use m_nctk
37 : use m_hdr
38 : use m_dtset
39 : use m_dtfil
40 : use m_krank
41 : use m_yaml
42 :
43 : use defs_datatypes, only : pseudopotential_type
44 : use defs_abitypes, only : MPI_type
45 : use m_time, only : cwtime, cwtime_report
46 : use m_io_tools, only : delete_file, get_unit, open_file
47 : use m_hide_lapack, only : matrginv, xheev
48 : use m_fstrings, only : strcat, sjoin, itoa
49 : use m_numeric_tools, only : uniformrandom, simpson_int, c2r, l2int, isdiagmat, get_diag, blocked_loop
50 : use m_special_funcs, only : besjm
51 : use m_geometry, only : xred2xcart, rotmat, wigner_seitz
52 : use m_crystal, only : crystal_t
53 : use m_fftcore, only : sphereboundary
54 : use m_ebands, only : ebands_t
55 : use m_pawang, only : pawang_type
56 : use m_pawrad, only : pawrad_type, simp_gen
57 : use m_pawtab, only : pawtab_type
58 : use m_pawcprj, only : pawcprj_type
59 : use m_paw_sphharm, only : ylm_cmplx, initylmr
60 : use m_paw_overlap, only : smatrix_pawinit
61 : use m_evdw_wannier, only : evdw_wannier
62 : use m_abstract_wf, only: abstract_wf, wann_ksetting_t, cg_cprj, wfd_wf
63 : use m_wannier_io, only: write_eigenvalues, write_Amn, compute_and_write_unk, write_Mmn
64 :
65 : implicit none
66 :
67 : private
68 : !!***
69 :
70 : public :: mlwfovlp
71 : !!***
72 :
73 : !----------------------------------------------------------------------
74 :
75 : !!****t* m_mlwfovlp/wan_t
76 : !! NAME
77 : !! wan_t
78 : !!
79 : !! FUNCTION
80 : !! This object stores the results of the Wannnierization algorithm.
81 : !! It can be constructed by reading the ABIWAN.nc file produced by
82 : !! Abinit when we call wannier90 in library mode. See mlwfovlp routine.
83 : !!
84 : !! SOURCE
85 :
86 : type,public :: wan_t
87 :
88 : integer :: spin = -1
89 : ! Spin index.
90 :
91 : integer :: nwan = -1
92 : ! Number of Wannier functions.
93 :
94 : integer :: max_nwan = -1
95 : ! Max number of Wannier functions over spins (used to dimension arrays)
96 :
97 : integer :: num_bands = -1
98 : ! Number of bands seen by wannier90 for this spin.
99 :
100 : !integer :: nbndep, ! Number of remaining bands after excluding bands in Wannierizatin step
101 : !integer :: nbndskip, ! Number of bands to be skipped in Wannierization step, leading to
102 : ! the exclusion from the original Hamiltonian
103 : integer :: nkbz = -1
104 : ! Number of k-points in the full BZ.
105 :
106 : integer :: nr_h = -1, nr_e = -1, nr_p = -1
107 : ! Number of lattice points for H, electrons, phonons.
108 :
109 : integer :: ngkpt(3) = -1
110 : ! K-mesh divisions.
111 :
112 : !integer :: nshiftk
113 : ! Number of shifts. At present only 1 shift is supported.
114 :
115 : logical :: have_disentangled
116 : ! True if disentanglement has been used.
117 :
118 : real(dp) :: spread(3) = -one
119 : ! Spread of wannier functions.
120 :
121 : type(krank_t) :: krank
122 : ! Used to find the index of the kpoint from its coordinates.
123 :
124 : integer,allocatable :: exclude_bands(:)
125 : ! FIXME: Is this still needed?
126 :
127 : integer,allocatable :: dimwin(:), winstart(:)
128 : ! (nkbz)
129 : ! Number of bands within the outer window at each k-point
130 :
131 : integer :: bmin = -1, bmax = -1
132 : ! Minimum and maximum band included in the Wannierization.
133 :
134 : integer,allocatable :: r_h(:,:), r_e(:,:), r_p(:,:)
135 : ! Lattice points for H in the Wannier representation
136 : ! (3, nr_h)
137 :
138 : integer,allocatable :: ndegen_h(:), ndegen_e(:), ndegen_p(:)
139 :
140 : real(dp),allocatable :: rmod_h(:), rmod_p(:), rmod_e(:)
141 : ! (nr_h), (nr_p)
142 : ! Lenght of the lattice points in Bohr
143 :
144 : real(dp),allocatable :: all_eigens(:,:)
145 : ! (mband, nkbz)
146 : ! All KS eigenvalues (before possible filtering done by wannier90)
147 :
148 : real(dp),allocatable :: centres(:,:)
149 : ! Wannier centers.
150 : ! (3, nwan)
151 :
152 : real(dp),allocatable :: spreads(:)
153 : ! Wannier centers.
154 : ! (nwan)
155 :
156 : real(dp),allocatable :: kbz(:,:)
157 : ! kpoints in the full BZ used by Wannier90.
158 : ! (3, nkbz)
159 :
160 : logical,allocatable :: band_in(:)
161 : ! (num_bands)
162 :
163 : logical,allocatable :: lwindow(:,:)
164 : ! (num_bands, nkbz)
165 :
166 : logical :: keep_umats
167 : ! True of u_mat_opt and u_mat should be saved in memory.
168 :
169 : complex(dp),allocatable :: u_mat_opt(:,:,:)
170 : complex(dp),allocatable :: u_mat(:,:,:)
171 :
172 : complex(dp),allocatable :: u_k(:,:,:)
173 : ! (max_dimwin, nwan, nkbz)
174 : ! total rotation matrix: the product of the optimal subspace x the rotation among the nwan Wannier functions.
175 : ! on the coarse ab-initio k-mesh
176 :
177 : complex(dp),allocatable :: hwan_r(:,:,:)
178 : ! (nr_h, nwan, nwan)
179 : ! KS Hamiltonian in the Wannier representation.
180 :
181 : integer :: my_npert = -1, my_pert_start = -1
182 : ! My number of perturbations and my initial perturbation.
183 :
184 : type(xcomm_t), pointer :: pert_comm => null()
185 : ! MPI-communicator for parallelism over perturbations.
186 :
187 : complex(dp),allocatable :: grpe_wwp(:,:,:,:,:)
188 : ! (nr_p, nr_e, nwan, nwan, my_npert))
189 : ! e-ph matrix elements in the Wannier representation.
190 : ! NB: These matrix elements are in the atomic represention and distributed inside pert_comm.
191 :
192 : contains
193 : procedure :: from_abiwan => wan_from_abiwan
194 : ! Initialize a wan_t instance from a ABIWAN.nc file
195 :
196 : procedure :: load_gwan => wan_load_gwan
197 : ! Read g(R_p, R_e) in the Wannier representation from the GWAN.nc file.
198 :
199 : procedure :: print => wan_print
200 : ! Print info on the object.
201 :
202 : procedure :: interp_ham => wan_interp_ham
203 : ! Interpolate Hamiltonian at an arbitray k-point.
204 : ! energies, and optionally diagonal velocities.
205 :
206 : procedure :: setup_eph_ws_kq => wan_setup_eph_ws_kq
207 : ! Prepare interpolation of e-ph matrix elements.
208 :
209 : procedure :: interp_eph_manyq => wan_interp_eph_manyq
210 : ! Interpolate e-ph matrix elements.
211 :
212 : procedure :: interp_eph_manyk => wan_interp_eph_manyk
213 : ! Interpolate e-ph matrix elements for many k points at fixed q.
214 :
215 : procedure :: prepare_eph_q => wan_prepare_eph_q
216 : ! Fourier transform the e-ph vertex from R_p to one q point.
217 :
218 : procedure :: interp_eph_manyk_from_q => wan_interp_eph_manyk_from_q
219 : ! Interpolate a k batch from a q-prepared e-ph vertex.
220 :
221 : procedure :: eph_kbatch_size => wan_eph_kbatch_size
222 : ! Select a k-batch size from a configurable workspace limit.
223 :
224 : procedure :: ncwrite_gwan => wan_ncwrite_gwan
225 : ! Write g in the Wannier representation to netcdf file.
226 :
227 : procedure :: free => wan_free
228 : ! Free memory.
229 :
230 : end type wan_t
231 : !!***
232 :
233 : public :: wan_interp_ebands
234 : ! Build new ebands_t object on a k-mesh via Wannier interpolation.
235 :
236 : contains
237 : !!***
238 :
239 : !!****f* m_mlwfovlp/mlwfovlp
240 : !! NAME
241 : !! mlwfovlp
242 : !!
243 : !! FUNCTION
244 : !! Routine which computes overlap M_{mn}(k,b) and projection A_{mn}(k)
245 : !! for Wannier code (www.wannier.org f90 version).
246 : !! Various file are written (wannier90.*) which can be used to run a
247 : !! separate wannier calculation with the wannier90 code.
248 : !!
249 : !! INPUTS
250 : !! crystal<crystal_t>=Info on the crystalline structure.
251 : !! ebands<ebands_t>=The object describing the band structure.
252 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
253 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
254 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
255 : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
256 : !! dtset <type(dataset_type)>=all input variables for this dataset
257 : !! dtfil <type(datafiles_type)>=variables related to files
258 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
259 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
260 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
261 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
262 : !! mband=maximum number of bands
263 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
264 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
265 : !! mgfft=maximum size of 1D FFTs
266 : !! mgfftc=maximum size of 1D FFTs (coarse grid)
267 : !! mkmem =number of k points treated by this node.
268 : !! mpi_enreg=information about MPI parallelization
269 : !! mpw=maximum dimensioned size of npw.
270 : !! natom=number of atoms in cell.
271 : !! nattyp(ntypat)= # atoms of each type.
272 : !! nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
273 : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
274 : !! nkpt=number of k points.
275 : !! npwarr(nkpt)=number of planewaves in basis at this k point
276 : !! nsppol=1 for unpolarized, 2 for spin-polarized
277 : !! ntypat=number of types of atoms in unit cell.
278 : !! occ(mband*nkpt*nsppol) Occupation number for each band (often 2) for each k point.
279 : !! prtvol=control print volume and debugging output
280 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
281 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
282 : !! ucvol=unit cell volume (bohr**3)
283 : !! xred(3,natom)=reduced dimensionless atomic coordinates
284 : !!
285 : !! OUTPUT
286 : !! (only writing, printing)
287 : !!
288 : !! SIDE EFFECTS
289 : !! (only writing, printing)
290 : !!
291 : !! NOTES
292 : !!
293 : !! SOURCE
294 :
295 19 : subroutine mlwfovlp(mywfc, crystal, ebands, hdr, atindx1, &
296 : !&cg,cprj, &
297 19 : dtset,dtfil,eigen,gprimd,kg,&
298 : mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
299 38 : nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
300 19 : pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred, exclude_bands)
301 :
302 : !Arguments ------------------------------------
303 : !scalars
304 : class(abstract_wf), pointer :: mywfc
305 : integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpw,natom,nfft,nkpt
306 : integer,intent(in) :: nsppol,ntypat,prtvol
307 : real(dp),intent(in) :: ucvol
308 : type(crystal_t),intent(in) :: crystal
309 : type(ebands_t),intent(in) :: ebands
310 : type(hdr_type),intent(in) :: hdr
311 : type(MPI_type),intent(inout) :: mpi_enreg
312 : type(dataset_type),intent(in) :: dtset
313 : type(datafiles_type),intent(in) :: dtfil
314 : type(pawang_type),intent(in) :: pawang
315 : type(pseudopotential_type),intent(in) :: psps
316 : !arrays
317 : integer,intent(in) :: atindx1(natom)
318 : integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt)
319 : ! real(dp), optional, intent(in) :: cg(2,mcg)
320 : ! type(pawcprj_type), optional, intent(in) :: cprj(natom,mcprj)
321 : real(dp),optional, intent(in) :: eigen(mband*nkpt*nsppol),gprimd(3,3),rprimd(3,3)
322 : real(dp),intent(in) :: occ(mband*nkpt*nsppol)
323 : real(dp),intent(in) :: xred(3,natom)
324 : !type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
325 : !type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
326 : type(pawrad_type),intent(in) :: pawrad(:)
327 : type(pawtab_type),intent(in) :: pawtab(:)
328 : integer, intent(inout) :: exclude_bands(:,:)
329 :
330 : !Local variables-------------------------------
331 : !scalars
332 : integer :: i, ierr, ikpt1,ikpt2,intot,isppol, iwan, lproj,lwanniersetup,mwan
333 : #if defined HAVE_WANNIER90
334 : integer :: ncid, ncerr, nrpts
335 : character(len=fnlen) :: abiwan_fname
336 38 : integer :: have_disentangled_spin(nsppol)
337 19 : integer,allocatable :: irvec_r_h(:,:),ndegen_h(:)
338 19 : real(dp),allocatable :: rmods_r_h(:)
339 : !type(wan_t) :: wan
340 : #endif
341 : integer,parameter :: master = 0
342 : integer :: nntot,num_nnmax
343 : integer :: max_num_bands,nprocs,comm,rank
344 57 : integer :: nwan(nsppol),nband_inc(nsppol),num_bands(nsppol)
345 : logical :: gamma_only,leig,lmmn,lwannierrun,spinors,test_matrix_output !,have_disentangled
346 : character(len=fnlen) :: wfnname
347 : character(len=1000) :: msg
348 38 : character(len=fnlen) :: seed_name(nsppol)
349 38 : character(len=fnlen) :: filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol)
350 38 : character(len=fnlen) :: filew90_mmn(nsppol),filew90_eig(nsppol)
351 : !arrays
352 : integer :: g1temp(3),ngkpt(3), units(2)
353 19 : integer,allocatable :: g1(:,:,:)
354 19 : integer,allocatable:: ovikp(:,:)
355 19 : integer,allocatable :: proj_l(:,:),proj_m(:,:),proj_radial(:,:)
356 19 : integer,allocatable :: proj_s_loc(:)
357 : real(dp) :: real_lattice(3,3), recip_lattice(3,3)
358 19 : real(dp),allocatable :: cm1(:,:,:,:,:,:),cm2_paw(:,:,:)
359 19 : real(dp),allocatable :: eigenvalues_w(:,:,:)
360 19 : real(dp),allocatable :: proj_site(:,:,:),proj_x(:,:,:),proj_z(:,:,:),proj_zona(:,:)
361 19 : real(dp),allocatable :: wann_centres(:,:,:),wann_spreads(:,:),xcart(:,:)
362 19 : real(dp),allocatable :: proj_s_qaxis_loc(:,:)
363 19 : complex(dp),allocatable :: A_paw(:,:,:,:)
364 19 : complex(dp),allocatable :: M_matrix(:,:,:,:,:),U_matrix(:,:,:,:)
365 19 : complex(dp),allocatable :: U_matrix_opt(:,:,:,:)
366 19 : complex(dp),pointer :: A_matrix(:,:,:,:)
367 19 : logical,allocatable :: band_in(:,:),lwindow(:,:,:)
368 19 : character(len=3),allocatable :: atom_symbols(:)
369 19 : logical,allocatable:: just_augmentation(:,:)
370 19 : type(pawcprj_type), pointer :: cprj_ptr(:, :)
371 : #ifdef HAVE_WANNIER90
372 38 : real(dp) :: spreadw(3,nsppol)
373 : #endif
374 : !************************************************************************
375 :
376 : ABI_UNUSED((/crystal%natom, ebands%nkpt, hdr%nkpt/))
377 : ABI_UNUSED(atindx1)
378 : ABI_UNUSED((/mcg, mcprj, prtvol/))
379 19 : ABI_UNUSED_A(pawang)
380 :
381 57 : units = [std_out, ab_out]
382 :
383 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384 : !1) Initialize variables and allocations
385 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
386 : !
387 : !Some initialization and checks
388 : !
389 19 : lwanniersetup=1 ! 1 is mandatory ( 0 is for debug)
390 : !to use lwanniersetup=0, one would need to define which bands to exclude.
391 19 : lwannierrun=.true. ! .false. and .true. are possible
392 19 : lmmn=.true. ! .false. and .true. are possible
393 19 : leig=.true. ! .false. and .true. are possible
394 19 : gamma_only=.false. !not yet implemented
395 19 : spinors=.false.
396 19 : if (dtset%nspinor == 2) spinors = .true.
397 : ! Raw A_mn and M_mn entries are suitable regression quantities only for
398 : ! scalar wavefunctions. Spinor eigenvectors have additional gauge freedom.
399 19 : test_matrix_output = dtset%nspinor == 1
400 :
401 : !mpi initialization
402 19 : comm=MPI_enreg%comm_cell
403 19 : nprocs=xmpi_comm_size(comm)
404 19 : rank=MPI_enreg%me_kpt
405 : !write(std_out,'("master ",i0," rank ",i0," nprocs ",i0)') master,rank,nprocs
406 :
407 : !Generate seed names for wannier90 files, and file names
408 : call mlwfovlp_seedname(dtfil%fnameabo_w90,filew90_win,filew90_wout,filew90_amn,&
409 19 : filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name)
410 :
411 19 : call hdr%vs_dtset(dtset)
412 :
413 : ! Check that the WFK contains the full Brillouin-zone mesh required by Wannier90.
414 : ! Use the arrays stored in the WFK header: dtset%istwfk and dtset%wtk may only
415 : ! contain the compact input representation and therefore need not have nkpt entries.
416 19 : if (rank == master) then
417 1299 : if (.not. (all(hdr%istwfk(1:nkpt) == 1) .and. all(hdr%wtk(1:nkpt) == hdr%wtk(1)))) then
418 0 : write(msg, '(4a)') ch10, &
419 0 : ' mlwfovlp: you should provide k-points in the full Brillouin zone ', ch10, &
420 0 : ' with explicit k-points (or kptopt=3) and istwfk 1'
421 0 : call wrtout(units, msg)
422 0 : ABI_ERROR(msg)
423 : end if
424 : end if
425 :
426 19 : if(MPI_enreg%paral_spinor==1) then
427 0 : ABI_ERROR('Parallelization over spinorial components not yet available !')
428 : end if
429 :
430 : ! MG: TODO: Why this check?
431 19 : if (psps%npsp/=psps%ntypat) then
432 0 : ABI_ERROR("Alchemical mixing not supported in mlwfovlp")
433 : end if
434 :
435 19 : if (nsppol==2) then
436 1 : write(msg, '(3a)' ) ch10,' mlwfovlp: Calculating matrices for both spin polarization ',ch10
437 1 : call wrtout(units, msg)
438 : end if
439 :
440 : ! get lattice parameters in wannier90 format
441 76 : do i=1, 3
442 228 : real_lattice(:,i)=Bohr_Ang*rprimd(i,:)
443 247 : recip_lattice(:,i)=two_pi*gprimd(i,:)/Bohr_Ang
444 : end do
445 :
446 : ! Allocations.
447 19 : num_nnmax=12 ! limit fixed for compact structure in wannier_setup.
448 76 : ABI_MALLOC(g1,(3,nkpt,num_nnmax))
449 76 : ABI_MALLOC(ovikp,(nkpt,num_nnmax))
450 57 : ABI_MALLOC(atom_symbols,(natom))
451 57 : ABI_MALLOC(xcart,(3,natom))
452 76 : ABI_MALLOC(band_in,(mband,nsppol))
453 76 : ABI_MALLOC(proj_site,(3,mband,nsppol))
454 57 : ABI_MALLOC(proj_l,(mband,nsppol))
455 57 : ABI_MALLOC(proj_m,(mband,nsppol))
456 57 : ABI_MALLOC(proj_radial,(mband,nsppol))
457 57 : ABI_MALLOC(proj_x,(3,mband,nsppol))
458 57 : ABI_MALLOC(proj_s_loc,(mband))
459 57 : ABI_MALLOC(proj_s_qaxis_loc,(3,mband))
460 57 : ABI_MALLOC(proj_z,(3,mband,nsppol))
461 76 : ABI_MALLOC(proj_zona,(mband,nsppol))
462 : ! ABI_MALLOC(exclude_bands, (mband,nsppol))
463 : !
464 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
465 : !2) Call to Wannier setup
466 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
467 : !
468 19 : nullify(A_matrix)
469 :
470 : call mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,&
471 : g1,lwanniersetup,mband,natom,nband_inc,nkpt,&
472 : nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,&
473 : proj_l,proj_m,proj_radial,proj_site,proj_s_loc, proj_s_qaxis_loc, proj_x,proj_z,proj_zona,&
474 19 : real_lattice,recip_lattice,rprimd,seed_name,spinors,xcart,xred,exclude_bands)
475 :
476 39 : do isppol=1, nsppol
477 20 : write(msg, '(6a)' ) ch10,&
478 20 : ' mlwfovlp : mlwfovlp_setup done -',ch10,&
479 40 : '- see ',trim(filew90_wout(isppol)),' for details.'
480 39 : call wrtout(units, msg)
481 : end do
482 :
483 : ! some allocations after wannier90 setup
484 39 : max_num_bands=maxval(num_bands(:))
485 39 : mwan=maxval(nwan(:))
486 95 : ABI_MALLOC(eigenvalues_w,(max_num_bands,nkpt,nsppol))
487 133 : ABI_MALLOC(M_matrix,(max_num_bands,max_num_bands,nntot,nkpt,nsppol))
488 114 : ABI_MALLOC(A_matrix,(max_num_bands,mwan,nkpt,nsppol))
489 :
490 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
491 : !3) Write Eigenvalues (file seed_name.eig)
492 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
493 : if (leig) then
494 19 : call write_eigenvalues(filew90_eig,eigen, band_in, eigenvalues_w, nsppol, nkpt, mband, dtset, rank, master )
495 : end if !leig
496 :
497 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
498 : !4) Calculate overlaps (file seed_name.mmn)
499 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
500 : !
501 : !First calculate indices and shift
502 : !
503 : !write(std_out,*) "Computes shift for cg"
504 19 : write(msg, '(a,a)' ) ch10,' mlwfovlp : compute shifts for g-points '
505 19 : call wrtout(std_out, msg)
506 : !----------------------------------------------------------------------
507 : !Compute shifts for g points (icg,iwav)
508 : !(here mband is not used, because shifts are internal variables of abinit)
509 : !----------------------------------------------------------------------
510 : !call mywfc%init(cg, cprj, dtset, dtfil, hdr, &
511 : ! & MPI_enreg, nprocs, psps, pawtab, rank)
512 :
513 : !TODO uncomment
514 : ! call mywfc%kset%set_ovikp( ovikp=ovikp, nntot=nntot, num_nnmax=num_nnmax)
515 : !
516 : !Shifts computed.
517 : !
518 : if( lmmn) then
519 : !
520 : ! In case of parallelization write out cg for all k-points
521 : !
522 : !if (nprocs > 1) then
523 : ! call write_cg_and_cprj(dtset, cg, cprj, dtfil, iwav, npwarr, mband, natom, &
524 : ! &nsppol, nkpt, MPI_enreg, rank, psps, pawtab)
525 : ! call mywfc%write_cg_and_cprj_tmpfile()
526 : !end if !MPI nprocs>1
527 : !
528 : ! End of MPI preliminarities
529 : ! Calculate PW contribution of overlaps
530 : !
531 133 : ABI_MALLOC(cm1,(2,mband,mband,nntot,nkpt,nsppol))
532 : ! this loops over spin internally
533 : call mlwfovlp_pw(mywfc,cm1,g1,kg,mband, mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
534 19 : npwarr,hdr%nspinor,nsppol,ovikp)
535 :
536 : !mlwfovlp_pw(mywfc,cm1,g1,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
537 : ! & npwarr,nspinor,nsppol,ovikp,seed_name)
538 19 : write(msg, '(a,a)' ) ch10,' mlwfovlp: PW part of overlap computed '
539 19 : call wrtout(std_out, msg)
540 : !
541 : ! compute PAW Contribution and add it to PW contribution
542 : !
543 19 : if(psps%usepaw==1) then
544 5 : write(msg, '(a,a)' ) ch10,'** smatrix_pawinit: PAW part of overlap '
545 5 : call wrtout(std_out, msg)
546 20 : ABI_MALLOC(cm2_paw,(2,mband,mband))
547 11 : do isppol=1,nsppol
548 283 : do ikpt1=1,nkpt
549 : ! MPI:cycle over k-points not treated by this node
550 272 : if (nprocs>1 ) then !sometimes we can have just one processor
551 0 : if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-rank) /=0) CYCLE
552 : end if
553 :
554 272 : write(msg, '(a,i0,a,2(i0,1x))')' processor',rank,' computes PAW part for kpt and spin',ikpt1,isppol
555 272 : call wrtout(std_out, msg)
556 :
557 3222 : do intot=1,nntot
558 2944 : ikpt2= ovikp(ikpt1,intot)
559 11776 : g1temp(:)=g1(:,ikpt1,intot)
560 : ! TODO : smatrix_pawinit: use high level wfd.
561 2944 : cprj_ptr=>mywfc%get_cprj_ptr()
562 : call smatrix_pawinit(atindx1,cm2_paw,cprj_ptr,ikpt1,ikpt2,isppol,&
563 : & g1temp,gprimd,hdr%kptns,mband,mband,mkmem,mpi_enreg,&
564 : & natom,dtset%nband,nkpt,dtset%nspinor,nsppol,dtset%ntypat,pawang,pawrad,pawtab,rprimd,&
565 2944 : & dtfil%fnametmp_cprj,dtset%typat,xred)
566 : ! cm1(:,:,:,intot,ikpt1,isppol)=four_pi*cm2_paw(:,:,:)
567 : ! write(6,*) "ikpt1=",ikpt1
568 : ! do iband=1,mband
569 : ! write(6,*) "iband=",iband
570 : ! write(6,*) "Wannier PW overlap",cm1(:,iband,iband,intot,ikpt1,isppol)
571 : ! write(6,*) "Wannier PAW overlap",four_pi*cm2_paw(:,iband,iband)
572 : ! write(6,*) "Wannier PW+PAW overlap",cm1(:,iband,iband,intot,ikpt1,isppol)+four_pi*cm2_paw(:,iband,iband)
573 : ! enddo
574 2473744 : cm1(:,:,:,intot,ikpt1,isppol)=cm1(:,:,:,intot,ikpt1,isppol)+four_pi*cm2_paw(:,:,:)
575 : end do ! intot
576 : end do ! ikpt1
577 : end do ! isppol
578 5 : ABI_FREE(cm2_paw)
579 5 : write(msg, '(a,a)' ) ch10,' mlwfovlp : PAW part of overlap computed '
580 5 : call wrtout(std_out, msg)
581 : end if ! usepaw
582 :
583 19 : call xmpi_sum(cm1,comm,ierr)
584 :
585 : call write_Mmn(filew90_mmn, band_in, cm1, ovikp, g1, M_matrix, nkpt, nsppol, nntot, mband, num_bands, msg, &
586 19 : iam_master=(rank==master), test_matrix_output=test_matrix_output)
587 38 : ABI_FREE(cm1)
588 :
589 : ! erase temporary files created for parallel runs
590 : !if (nprocs > 1) call mywfc%remove_tmpfile(prtvol)
591 : !end if !MPI nprocs>1
592 : end if !lmmn
593 :
594 : !Deallocate arrays no longer needed
595 19 : ABI_FREE(ovikp)
596 19 : ABI_FREE(g1)
597 :
598 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
599 : !5) Calculate initial projections
600 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
601 :
602 19 : if (dtset%w90iniprj/=0) then
603 :
604 : ! Set value for lproj (type of projections to be computed)
605 : ! In PAW, options 5 and 6 are not in use.
606 : ! 5 means that there will be a contribution from inside the spheres and another from the PW part
607 : ! 6 means that we take into account just the inside-spheres contribution
608 : ! 2 means that PW part will be calculated
609 19 : lproj=dtset%w90iniprj
610 19 : if(dtset%w90iniprj == 5 ) lproj=2 ! Necessary to calculate PW contribution
611 76 : ABI_MALLOC(just_augmentation,(mwan,nsppol))
612 151 : just_augmentation(:,:)=.false.
613 :
614 19 : if( psps%usepaw==1 .and. (dtset%w90iniprj==2 .or. dtset%w90iniprj>4)) then
615 4 : if (dtset%w90iniprj==6) just_augmentation(:,:)=.true.
616 4 : if (dtset%w90iniprj==5) then
617 0 : do isppol=1,nsppol
618 0 : do iwan=1,nwan(isppol)
619 : ! Trick to skip the planewave contribution for some Wannier functions
620 : ! (Not in production).
621 0 : if(proj_radial(iwan,isppol) > 4) then
622 0 : just_augmentation(iwan,isppol)=.true.
623 0 : proj_radial(iwan,isppol)=proj_radial(iwan,isppol)-3
624 0 : write(msg, '(2a,2(i0,1x))')' ','Skiping planewave contribution for iwan, ispin=',iwan,isppol
625 0 : call wrtout(std_out, msg)
626 : end if !proj_radial>4
627 : end do !iwan
628 : end do !isppol
629 : end if !w90iniprj == 5
630 : end if !paw
631 : !
632 : ! Call mlwfovlp_proj (plane waves part of projections)
633 : !
634 19 : if (dtset%w90iniprj/=6) then ! option 6 not yet in use
635 : ! call mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,&
636 : !& lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,&
637 : !& nattyp,nkpt,npwarr,&
638 : !& dtset%nspinor,nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,&
639 : !& proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,ucvol)
640 : call mlwfovlp_proj(A_matrix,band_in,mywfc,dtset,gprimd,just_augmentation,kg,&
641 : & lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,&
642 : & nattyp,nkpt,npwarr,&
643 : & dtset%nspinor,nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,&
644 19 : & proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,ucvol)
645 19 : write(msg, '(a,a,a,a)' ) ch10,&
646 19 : ' mlwfovlp: mlwfovlp_proj done -',ch10,&
647 38 : ' Projectors computed.'
648 19 : call wrtout(std_out, msg)
649 : end if !w90proj/=6
650 : !
651 : ! Calculate inside-sphere part of projections (PAW)
652 : !
653 19 : if (psps%usepaw ==1 .and. ( dtset%w90iniprj>4)) then
654 0 : ABI_MALLOC(A_paw,(max_num_bands,mwan,nkpt,nsppol))
655 : call mlwfovlp_projpaw(A_paw,band_in,mywfc,just_augmentation,max_num_bands,mband,mkmem,&
656 : & mwan,natom,dtset%nband,nkpt,&
657 : & dtset%nspinor,nsppol,dtset%ntypat,nwan,pawrad,pawtab,&
658 : & proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,&
659 0 : & rprimd,dtset%typat,xred)
660 : !
661 0 : write(msg, '(a,a,a,a)' ) ch10,&
662 0 : ' mlwfovlp: mlwfovlp_proj_paw done -',ch10,&
663 0 : ' Inside-spheres part of projectors computed.'
664 0 : call wrtout(std_out, msg)
665 : !
666 : ! Add in-sphere contribution to A_matrix
667 : !
668 : !
669 : ! w90iniprj==5. Plane waves + augmentation contributions
670 : !
671 0 : if(dtset%w90iniprj==5) A_matrix(:,:,:,:)=A_matrix(:,:,:,:)+A_paw(:,:,:,:)
672 : !
673 : ! w90iniprj==6. Just augmentation contribution
674 : !
675 0 : if(dtset%w90iniprj==6) A_matrix(:,:,:,:)=A_paw(:,:,:,:)
676 0 : ABI_FREE(A_paw)
677 : end if !usepaw==1
678 :
679 19 : ABI_FREE(just_augmentation)
680 19 : call xmpi_sum(A_matrix,comm,ierr)
681 :
682 : ! write projections to file
683 19 : if (rank==master) then
684 19 : if(dtset%w90iniprj==1) then
685 6 : call write_Amn(A_matrix, filew90_ramn, nsppol, mband, nkpt, num_bands, nwan, band_in, test_matrix_output)
686 : else
687 13 : call write_Amn(A_matrix, filew90_amn, nsppol, mband, nkpt, num_bands, nwan, band_in, test_matrix_output)
688 : end if
689 : end if
690 : end if !dtset%w90iniprj/=0
691 :
692 : ! Deallocations
693 19 : ABI_FREE(proj_site)
694 19 : ABI_FREE(proj_l)
695 19 : ABI_FREE(proj_m)
696 19 : ABI_FREE(proj_radial)
697 19 : ABI_FREE(proj_x)
698 19 : ABI_FREE(proj_z)
699 19 : ABI_FREE(proj_zona)
700 19 : ABI_FREE(proj_s_loc)
701 19 : ABI_FREE(proj_s_qaxis_loc)
702 :
703 :
704 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
705 : !6) write files for wannier function plot
706 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
707 19 : if( dtset%w90prtunk>0) then
708 : call compute_and_write_unk(wfnname, psps%usepaw, dtset%w90prtunk, &
709 : mpi_enreg, ngfft, nsppol, dtset%nspinor, &
710 : nkpt, mband, mpw, mgfftc, mkmem, nprocs, rank, npwarr, &
711 4 : band_in, dtset, kg, mywfc)
712 : end if
713 :
714 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
715 : !7) Call to Wannier90
716 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
717 : if (lwannierrun) then
718 19 : if (lwanniersetup.ne.1) ABI_ERROR("lwanniersetup.ne.1")
719 114 : ABI_MALLOC(U_matrix,(mwan,mwan,nkpt,nsppol))
720 114 : ABI_MALLOC(U_matrix_opt,(max_num_bands,mwan,nkpt,nsppol))
721 95 : ABI_MALLOC(lwindow,(max_num_bands,nkpt,nsppol))
722 76 : ABI_MALLOC(wann_centres,(3,mwan,nsppol))
723 76 : ABI_MALLOC(wann_spreads,(mwan,nsppol))
724 : ! Initialize
725 48135 : U_matrix(:,:,:,:)=czero
726 65031 : U_matrix_opt(:,:,:,:)=czero
727 7399 : lwindow(:,:,:)=.false.
728 487 : wann_centres(:,:,:)=zero
729 151 : wann_spreads(:,:)=zero
730 : ! write(std_out,*) seed_name; write(std_out,*) ngkpt
731 19 : ABI_CHECK(isdiagmat(dtset%kptrlatt), "kptrlatt should be diagonal! Please use nkgpt with one shift")
732 19 : ngkpt(1)=dtset%kptrlatt(1,1)
733 19 : ngkpt(2)=dtset%kptrlatt(2,2)
734 19 : ngkpt(3)=dtset%kptrlatt(3,3)
735 :
736 : ! write(std_out,*) nkpt
737 : ! write(std_out,*) rprimd*Bohr_Ang
738 : ! write(std_out,*) two_pi*gprimd/Bohr_Ang
739 : ! write(std_out,*) mband
740 : ! write(std_out,*) "nwan",nwan
741 : ! write(std_out,*) nntot
742 : ! write(std_out,*) natom
743 : ! write(std_out,*) atom_symbols
744 : ! write(std_out,*) xcart
745 : ! write(std_out,*) num_bands,num_bands,nntot,nkpt
746 : ! write(std_out,*) wann_spreads
747 : ! wann_spreads=2
748 : ! do i=1, nkpt
749 : ! do j=1, nntot
750 : ! write(std_out,*) i,j
751 : ! do k=1, num_bands
752 : ! do l=1, num_bands
753 : ! write(std_out,*) "m",M_matrix(l,k,j,i,1)
754 : ! enddo
755 : ! enddo
756 : ! enddo
757 : ! enddo
758 :
759 : #if defined HAVE_WANNIER90
760 39 : do isppol=1,nsppol
761 : ! when nsppol>1, master runs isppol 1 and rank==1 runs isppol 2
762 20 : if(nprocs>1 .and. isppol==1.and.rank.ne.master) cycle
763 20 : if(nprocs>1 .and. isppol==2.and.rank.ne.1) cycle
764 :
765 20 : write(msg, '(8a)' ) ch10,&
766 20 : '** mlwfovlp: calling wannier90 library subroutine wannier_run ',ch10,&
767 20 : ' Calculation is running ',ch10,&
768 40 : '- see ',trim(filew90_wout(isppol)),' for details.'
769 20 : call wrtout(units, msg)
770 :
771 20 : write(msg, '(a,i0)') ' Spin channel: ', isppol
772 20 : call wrtout(units, msg)
773 20 : write(msg, '(a,3(i0,1x))') ' Uniform k-mesh: ', ngkpt
774 20 : call wrtout(units, msg)
775 20 : write(msg, '(a,i0)') ' Number of k-points: ', nkpt
776 20 : call wrtout(units, msg)
777 20 : write(msg, '(a,i0)') ' Number of input bands: ', num_bands(isppol)
778 20 : call wrtout(units, msg)
779 20 : write(msg, '(a,i0)') ' Number of Wannier functions: ', nwan(isppol)
780 20 : call wrtout(units, msg)
781 20 : write(msg, '(a,i0)') ' Number of k-point neighbours: ', nntot
782 20 : call wrtout(units, msg)
783 20 : write(msg, '(a,es16.8)') ' Fermi energy [eV]: ', ebands%fermie * Ha_eV
784 20 : call wrtout(units, msg)
785 20 : write(msg, '(a,2(es16.8,1x))') ' Eigenvalue range (min, max) [eV]: ', &
786 18504 : minval(ebands%eig(:,:,isppol)) * Ha_eV, maxval(ebands%eig(:,:,isppol)) * Ha_eV
787 20 : call wrtout(units, msg)
788 :
789 : call wannier_run(trim(seed_name(isppol)),ngkpt,nkpt,& ! input
790 : real_lattice,recip_lattice,hdr%kptns,num_bands(isppol),& ! input
791 : nwan(isppol),nntot,natom,atom_symbols,& ! input
792 : xcart*Bohr_Ang,gamma_only,M_matrix(:,:,:,:,isppol),A_matrix(:,:,:,isppol),eigenvalues_w(:,:,isppol),& ! input
793 : U_matrix(1:nwan(isppol),1:nwan(isppol),:,isppol),& ! output
794 : U_matrix_opt(1:num_bands(isppol),1:nwan(isppol),:,isppol),& ! output
795 : lwindow_loc=lwindow(1:num_bands(isppol),:,isppol),& ! output
796 : wann_centres_loc=wann_centres(:,1:nwan(isppol),isppol),& ! output
797 164 : wann_spreads_loc=wann_spreads(1:nwan(isppol),isppol),spread_loc=spreadw(:,isppol)) ! output
798 :
799 20 : write(msg, '(a)') ' Wannier90 physical results:'
800 20 : call wrtout(units, msg)
801 20 : write(msg, '(a,3(es16.8,1x))') '- Spreads (Omega_total, Omega_I, Omega_tilde) [Ang^2]: ', &
802 40 : spreadw(:,isppol)
803 20 : call wrtout(units, msg)
804 20 : write(msg, '(a)') ' Wannier function centres [Ang] and spreads [Ang^2]:'
805 20 : call wrtout(units, msg)
806 132 : do iwan=1,nwan(isppol)
807 112 : write(msg, '(a,i0,a,3(f14.8,1x),a,f14.8)') '- WF ', iwan, ': centre = ', &
808 224 : wann_centres(:,iwan,isppol), ' spread = ', wann_spreads(iwan,isppol)
809 132 : call wrtout(units, msg)
810 : end do
811 :
812 20 : write(msg, '(7a)' ) ch10,&
813 40 : ' mlwfovlp : mlwfovlp_run completed -',ch10,'- see ',trim(filew90_wout(isppol)),' for details.',ch10
814 39 : call wrtout(units, msg)
815 : end do !isppol
816 :
817 : ! collect output of wannier90 from different processors
818 19 : call xmpi_sum(U_matrix,comm,ierr)
819 19 : call xmpi_sum(U_matrix_opt,comm,ierr)
820 19 : call xmpi_lor(lwindow,comm)
821 19 : call xmpi_sum(wann_centres,comm,ierr)
822 19 : call xmpi_sum(wann_spreads,comm,ierr)
823 :
824 : ! Output ABIWAN.nc file
825 19 : if (hdr%kptopt == 0) then
826 4 : ABI_WARNING("Output of ABIWAN.nc requires kptopt /= 0. ABIWAN.nc file won't be produced!")
827 : ! Need kptrlatt in wigner_seitz and client code needs to know the k-grid.
828 : end if
829 :
830 19 : if (rank == master .and. hdr%kptopt /= 0) then
831 15 : abiwan_fname = strcat(dtfil%filnam_ds(4), "_ABIWAN.nc")
832 15 : call wrtout(std_out, sjoin(" Saving wannier90 ouput results to:", abiwan_fname))
833 : call wigner_seitz([zero, zero, zero], [2, 2, 2], dtset%kptrlatt, crystal%rmet, &
834 15 : nrpts, irvec_r_h, ndegen_h, rmods_r_h, prtvol=prtvol)
835 : ! We know if disentanglement has been done by looking at the output values of lwindow
836 : ! Not elegant but it is the only way to avoid the parsing of the wannier input.
837 : ! In wannier_run lwindow is set to True if not disentanglement
838 31 : have_disentangled_spin = 0
839 31 : do isppol=1,nsppol
840 : !if nwan(isppol) < num_bands(isppol)
841 4342 : if (.not. all(lwindow(:,:,isppol))) have_disentangled_spin(isppol) = 1
842 : end do
843 :
844 15 : NCF_CHECK(nctk_open_create(ncid, abiwan_fname, xmpi_comm_self))
845 15 : NCF_CHECK(hdr%ncwrite(ncid, fform_from_ext("ABIWAN"), nc_define=.True.))
846 15 : NCF_CHECK(crystal%ncwrite(ncid))
847 15 : NCF_CHECK(ebands%ncwrite(ncid))
848 :
849 : ncerr = nctk_def_dims(ncid, [ &
850 : nctkdim_t("mwan", mwan), &
851 : nctkdim_t("max_num_bands", max_num_bands), &
852 : nctkdim_t("nrpts", nrpts) &
853 60 : ], defmode=.True.)
854 15 : NCF_CHECK(ncerr)
855 :
856 30 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "nntot"])
857 15 : NCF_CHECK(ncerr)
858 : !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"])
859 : !NCF_CHECK(ncerr)
860 :
861 : ncerr = nctk_def_arrays(ncid, [ &
862 : nctkarr_t("nwan", "int", "number_of_spins"), &
863 : nctkarr_t("num_bands", "int", "number_of_spins"), &
864 : nctkarr_t("band_in_int", "int", "max_number_of_states, number_of_spins"), &
865 : nctkarr_t("lwindow_int", "int", "max_num_bands, number_of_kpoints, number_of_spins"), &
866 : nctkarr_t("exclude_bands", "int", "max_number_of_states, number_of_spins"), &
867 : nctkarr_t("spread", "dp", "three, number_of_spins"), &
868 : !nctkarr_t("A_matrix", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), &
869 : nctkarr_t("irvec", "int", "three, nrpts"), &
870 : nctkarr_t("ndegen", "int", "nrpts"), &
871 : nctkarr_t("have_disentangled_spin", "int", "number_of_spins"), &
872 : nctkarr_t("U_matrix", "dp", "two, mwan, mwan, number_of_kpoints, number_of_spins"), &
873 : nctkarr_t("U_matrix_opt", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), &
874 : nctkarr_t("wann_centres", "dp", "three, mwan, number_of_spins"), &
875 : nctkarr_t("wann_spreads", "dp", "mwan, number_of_spins") &
876 210 : ])
877 15 : NCF_CHECK(ncerr)
878 :
879 : ! Write data.
880 15 : NCF_CHECK(nctk_set_datamode(ncid))
881 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nntot"), nntot))
882 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nwan"), nwan))
883 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "num_bands"), num_bands))
884 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "band_in_int"), l2int(band_in)))
885 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "lwindow_int"), l2int(lwindow)))
886 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "exclude_bands"), exclude_bands))
887 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "spread"), spreadw))
888 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "A_matrix"), c2r(A_matrix)))
889 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "irvec"), irvec_r_h))
890 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ndegen"), ndegen_h))
891 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "have_disentangled_spin"), have_disentangled_spin))
892 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix"), c2r(U_matrix)))
893 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix_opt"), c2r(U_matrix_opt)))
894 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_centres"), wann_centres))
895 15 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_spreads"), wann_spreads))
896 15 : NCF_CHECK(nf90_close(ncid))
897 :
898 15 : ABI_FREE(irvec_r_h)
899 15 : ABI_FREE(ndegen_h)
900 15 : ABI_FREE(rmods_r_h)
901 : end if
902 :
903 : ! CALL SILVESTRELLI'S APPROACH TO EVALUATE vdW INTERACTION ENERGY USING MLWF!!
904 : ! ----------------------------------------------------------------------------------------------
905 19 : if (dtset%vdw_xc==10.or.dtset%vdw_xc==11.or.dtset%vdw_xc==12.or.dtset%vdw_xc==14.and.rank==master) then
906 : ! vdw_xc==10,11,12,14 starts the vdW interaction using MLWFs
907 3 : call evaluate_vdw_with_mlwf()
908 : end if
909 :
910 : #else
911 : ABI_UNUSED(occ)
912 : #endif
913 19 : ABI_FREE(wann_centres)
914 19 : ABI_FREE(wann_spreads)
915 19 : ABI_FREE(U_matrix)
916 19 : ABI_FREE(U_matrix_opt)
917 57 : ABI_FREE(lwindow)
918 : end if !lwannierrun
919 :
920 19 : ABI_FREE(band_in)
921 19 : ABI_FREE(atom_symbols)
922 19 : ABI_FREE(xcart)
923 19 : ABI_FREE(eigenvalues_w)
924 19 : ABI_FREE(M_matrix)
925 19 : ABI_FREE(A_matrix)
926 : !ABI_FREE(exclude_bands)
927 :
928 19 : call mywfc%free()
929 57 : ABI_FREE_SCALAR(mywfc)
930 :
931 : contains
932 : !!***
933 :
934 : !--------------------------------------------------------------------------------
935 : !> van der Waals interaction energy using MLWFs
936 : !--------------------------------------------------------------------------------
937 3 : subroutine evaluate_vdw_with_mlwf()
938 : #ifdef HAVE_WANNIER90
939 : integer :: ii, jj, ikpt, iband, kk
940 : real(dp) :: corrvdw
941 : complex(dp) :: caux,caux2,caux3
942 3 : real(dp),allocatable :: csix(:,:,:,:)
943 3 : real(dp),allocatable :: occ_arr(:,:,:),occ_wan(:,:,:)
944 3 : real(dp),allocatable :: tdocc_wan(:,:)
945 :
946 3 : write(std_out,*) 'nwan(nsppol)=',ch10
947 6 : do ii=1,nsppol
948 6 : write(std_out,*) 'nsppol=',ii, 'nwan(nsppol)=',nwan(ii),ch10
949 : end do
950 3 : write(std_out,*) 'mwan=', mwan, ch10
951 :
952 15 : ABI_MALLOC(occ_arr,(mband,nkpt,isppol))
953 15 : ABI_MALLOC(occ_wan,(mwan,nkpt,nsppol))
954 12 : ABI_MALLOC(tdocc_wan,(mwan,nsppol))
955 :
956 249 : occ_arr(:,:,:)=zero
957 126 : occ_wan(:,:,:)=zero
958 18 : tdocc_wan(:,:)=zero
959 3 : jj = 0
960 6 : do isppol=1,nsppol
961 30 : do ikpt=1,nkpt
962 123 : do iband=1,num_bands(isppol)
963 96 : jj = jj + 1
964 120 : occ_arr(iband,ikpt,isppol) = occ(jj)
965 : end do
966 : end do
967 : end do
968 :
969 6 : do isppol=1,nsppol
970 30 : do ikpt=1,nkpt
971 123 : do iwan=1,nwan(isppol)
972 96 : caux=czero
973 96 : caux2=czero
974 96 : caux3=czero
975 480 : do iband=1,num_bands(isppol) !nband_inc(isppol) !nwan(isppol)
976 2016 : do ii=1,nwan(isppol)
977 1536 : caux=U_matrix(ii,iwan,ikpt,isppol)*U_matrix_opt(iband,ii,ikpt,isppol)
978 : ! DEBUG
979 : ! if(ISNAN(dble(caux))) then
980 : ! write(std_out,*) 'NaN: caux(ikpt,iwan,iband,ii):',ikpt,iwan,iband,ii,ch10
981 : ! end if
982 : ! END DEBUG
983 8064 : do kk=1,nwan(isppol)
984 6144 : caux2=conjg(U_matrix(kk,iwan,ikpt,isppol))*conjg(U_matrix_opt(iband,kk,ikpt,isppol))
985 7680 : caux3= caux3+caux*caux2*occ_arr(iband,ikpt,isppol) !take care here as exclude_bands case is not well
986 : ! DEBUG
987 : ! if(ISNAN(dble(caux2))) then
988 : ! write(std_out,*) 'NaN: caux2(ikpt,iwan,iband,kk):',ikpt,iwan,iband,kk,ch10
989 : ! end if
990 : ! if(ISNAN(dble(caux3))) then
991 : ! write(std_out,*) 'NaN: caux3(ikpt,iwan,iband,kk,jj):',ikpt,iwan,iband,kk,jj
992 : ! end if
993 : ! END DEBUG
994 : end do
995 : end do
996 : end do
997 120 : occ_wan(iwan,ikpt,isppol) = dble(caux3)
998 : ! DEBUG
999 : ! write(std_out,*) occ_wan(iwan,ikpt,isppol)
1000 : ! END DEBUG
1001 : ! end do
1002 : end do
1003 : end do
1004 : end do
1005 :
1006 3 : write(std_out,*) ch10,'MLWFs Occupation Matrix diagonal terms:',ch10
1007 :
1008 6 : do jj=1,nsppol
1009 111 : forall(iwan=1:nwan(jj)) tdocc_wan(iwan,jj) = sum(occ_wan(iwan,1:nkpt,jj)) / real(nkpt,dp)
1010 3 : write(std_out,*) 'tdocc_wan(iwan),isppol:',ch10
1011 6 : write(std_out,*) (tdocc_wan(iwan,jj),iwan=1,nwan(jj)),jj
1012 : end do
1013 :
1014 18 : ABI_MALLOC(csix,(mwan,mwan,nsppol,nsppol))
1015 :
1016 : call evdw_wannier(csix,corrvdw,mwan,natom,nsppol,nwan,tdocc_wan,dtset%vdw_nfrag,&
1017 3 : dtset%vdw_supercell,dtset%vdw_typfrag,dtset%vdw_xc,rprimd,wann_centres,wann_spreads,xcart)
1018 :
1019 3 : ABI_FREE(csix)
1020 3 : ABI_FREE(occ_arr)
1021 3 : ABI_FREE(occ_wan)
1022 3 : ABI_FREE(tdocc_wan)
1023 :
1024 : #endif
1025 3 : end subroutine evaluate_vdw_with_mlwf
1026 :
1027 : end subroutine mlwfovlp
1028 : !!***
1029 :
1030 : !!****f* m_mlwfovlp/mlwfovlp_seedname
1031 : !! NAME
1032 : !! mlwfovlp_seedname
1033 : !!
1034 : !! FUNCTION
1035 : !! Get seed name and file names of all wannier90 related files
1036 : !!
1037 : !! INPUTS
1038 : !! fname_w90=root name of file appended with _w90
1039 : !!
1040 : !! OUTPUT
1041 : !! filew90_win= main input file for Wannier90
1042 : !! filew90_wout= main output file for Wannier90
1043 : !! filew90_amn= file containing Amn matrix
1044 : !! filew90_ramn= file containing Amn matrix (random initial projections)
1045 : !! filew90_mmn= file containing Mmn matrix
1046 : !! filew90_eig= file containing eigenvalues
1047 : !! nsppol= number of spin polarizations
1048 : !! seed_name= common seed name for all wannier90 related files
1049 : !!
1050 : !! SOURCE
1051 :
1052 19 : subroutine mlwfovlp_seedname(fname_w90,filew90_win,filew90_wout,filew90_amn,&
1053 19 : & filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name)
1054 :
1055 : !Arguments ------------------------------------
1056 : integer,intent(in) :: nsppol
1057 : character(len=fnlen),intent(out) :: filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol)
1058 : character(len=fnlen),intent(out) :: filew90_mmn(nsppol),filew90_eig(nsppol),seed_name(nsppol)
1059 : character(len=fnlen),intent(in) :: fname_w90
1060 :
1061 : !Local variables-------------------------------
1062 : integer:: isppol, units(2)
1063 : character(len=fnlen) :: test_win1,test_win2,test_win3
1064 : logical :: lfile
1065 : character(len=2000) :: msg
1066 : character(len=10):: postfix
1067 : ! *************************************************************************
1068 :
1069 57 : units = [std_out, ab_out]
1070 39 : seed_name(:)=trim(fname_w90)
1071 :
1072 39 : do isppol=1,nsppol
1073 20 : if(nsppol==1) postfix='.win'
1074 20 : if(nsppol==2 .and. isppol==1) postfix='_up.win'
1075 20 : if(nsppol==2 .and. isppol==2) postfix='_down.win'
1076 :
1077 20 : filew90_win(isppol) = trim(seed_name(isppol))//trim(postfix)
1078 20 : test_win1 = filew90_win(isppol)
1079 20 : inquire(file=filew90_win(isppol),exist=lfile)
1080 :
1081 20 : if(.not.lfile) then
1082 4 : seed_name(isppol)='wannier90'
1083 4 : filew90_win(isppol)=trim(seed_name(isppol))//trim(postfix)
1084 4 : test_win2=filew90_win(isppol)
1085 4 : inquire(file=filew90_win(isppol),exist=lfile)
1086 : end if
1087 :
1088 20 : if(.not.lfile) then
1089 3 : seed_name(isppol)='w90'
1090 9 : filew90_win=trim(seed_name(isppol))//trim(postfix)
1091 3 : test_win3=filew90_win(isppol)
1092 3 : inquire(file=filew90_win(isppol),exist=lfile)
1093 : end if
1094 :
1095 39 : if(.not. lfile) then
1096 : write(msg,'(12a)')&
1097 0 : ' wannier90 interface needs one of the following input files:',ch10,&
1098 0 : ' ',trim(test_win1),ch10,&
1099 0 : ' ',trim(test_win2),ch10,&
1100 0 : ' ',trim(test_win3),ch10,&
1101 0 : ' Action: read wannier90 tutorial and/or user manual and supply proper *.win file'
1102 0 : ABI_ERROR(msg)
1103 : end if
1104 : end do !isppol
1105 :
1106 : ! Files having different names for different spin polarizations
1107 19 : if(nsppol==1) then
1108 18 : filew90_win(1) =trim(seed_name(1))//'.win'
1109 18 : filew90_wout(1)=trim(seed_name(1))//'.wout'
1110 18 : filew90_ramn(1)=trim(seed_name(1))//'random.amn'
1111 18 : filew90_amn(1) =trim(seed_name(1))//'.amn'
1112 18 : filew90_mmn(1) =trim(seed_name(1))//'.mmn'
1113 18 : filew90_eig(1) =trim(seed_name(1))//'.eig'
1114 1 : elseif(nsppol==2) then
1115 1 : filew90_win(1) =trim(seed_name(1))//'_up.win'
1116 1 : filew90_win(2) =trim(seed_name(2))//'_down.win'
1117 1 : filew90_wout(1)=trim(seed_name(1))//'_up.wout'
1118 1 : filew90_wout(2)=trim(seed_name(2))//'_down.wout'
1119 1 : filew90_ramn(1)=trim(seed_name(1))//'random_up.amn'
1120 1 : filew90_ramn(2)=trim(seed_name(2))//'random_down.amn'
1121 1 : filew90_amn(1)=trim(seed_name(1))//'_up.amn'
1122 1 : filew90_amn(2)=trim(seed_name(2))//'_down.amn'
1123 1 : filew90_mmn(1)=trim(seed_name(1))//'_up.mmn'
1124 1 : filew90_mmn(2)=trim(seed_name(2))//'_down.mmn'
1125 1 : filew90_eig(1)=trim(seed_name(1))//'_up.eig'
1126 1 : filew90_eig(2)=trim(seed_name(2))//'_down.eig'
1127 : end if
1128 :
1129 : ! change also seed_name for nsppol=2
1130 19 : if(nsppol==2) then
1131 1 : seed_name(1)=trim(seed_name(1))//'_up'
1132 1 : seed_name(2)=trim(seed_name(2))//'_down'
1133 : end if
1134 :
1135 19 : write(msg, '(a,a)' ) ch10,'---------------------------------------------------------------'
1136 19 : call wrtout(units, msg)
1137 19 : write(msg, '(5a)' ) ch10,&
1138 19 : ' Calculation of overlap and call to wannier90 library ',ch10,&
1139 38 : ' to obtain maximally localized wannier functions ',ch10
1140 19 : call wrtout(units, msg)
1141 :
1142 19 : if(nsppol==1) then
1143 : write(msg, '(23a)' ) &
1144 18 : ' - ',trim(filew90_win(1)),' is a mandatory secondary input',ch10,&
1145 18 : ' - ',trim(filew90_wout(1)),' is the output for the library',ch10,&
1146 18 : ' - ',trim(filew90_ramn(1)),' contains random projections',ch10,&
1147 18 : ' - ',trim(filew90_amn(1)),' contains projections',ch10,&
1148 18 : ' - ',trim(filew90_mmn(1)),' contains the overlap',ch10,&
1149 36 : ' - ',trim(filew90_eig(1)),' contains the eigenvalues'
1150 1 : else if(nsppol==2) then
1151 : write(msg, '(41a)' ) &
1152 1 : ' - ',trim(filew90_win(1)),&
1153 1 : ' and ',trim(filew90_win(2)),ch10,'are mandatory secondary input',ch10,&
1154 1 : ' - ',trim(filew90_wout(1)),&
1155 1 : ' and ',trim(filew90_wout(2)),ch10,' are the output for the library',ch10,&
1156 1 : ' - ',trim(filew90_ramn(1)),&
1157 1 : ' and ',trim(filew90_ramn(2)),ch10,' contain random projections',ch10,&
1158 1 : ' - ',trim(filew90_amn(1)),&
1159 1 : ' and ',trim(filew90_amn(2)),ch10,' contain projections',ch10,&
1160 1 : ' - ',trim(filew90_mmn(1)),&
1161 1 : ' and ',trim(filew90_mmn(2)),ch10,' contain the overlap',ch10,&
1162 1 : ' - ',trim(filew90_eig(1)),&
1163 2 : ' and ',trim(filew90_eig(2)),ch10,' contain the eigenvalues'
1164 : end if
1165 19 : call wrtout(units, msg)
1166 :
1167 19 : write(msg, '(a,a)' ) ch10,'---------------------------------------------------------------'
1168 19 : call wrtout(units, msg)
1169 :
1170 19 : end subroutine mlwfovlp_seedname
1171 : !!***
1172 :
1173 : !!****f* m_mlwfovlp/mlwfovlp_setup
1174 : !! NAME
1175 : !! mlwfovlp_setup
1176 : !!
1177 : !! FUNCTION
1178 : !! Routine which creates table g1 and ovikp necessary to compute
1179 : !! overlap for Wannier code (www.wannier.org f90 version).
1180 : !!
1181 : !! INPUTS
1182 : !! atom_symbols(natom)= table of symbol for each atom
1183 : !! and each |p_lmn> non-local projector
1184 : !! dtset <type(dataset_type)>=all input variables for this dataset
1185 : !! filew90_win(nsppol) secondary input files for w90
1186 : !! lwanniersetup= flag: only 1 is fully working.
1187 : !! natom =number of atoms in cell.
1188 : !! mband=maximum number of bands
1189 : !! natom=number of atoms in cell.
1190 : !! nkpt=number of k points.
1191 : !! num_bands(isppol)=number of bands actually used to construct the wannier function
1192 : !! nwan(isppol)= number of wannier fonctions (read in wannier90.win).
1193 : !! dtset <type(dataset_type)>=all input variables for this dataset
1194 : !! real_lattice(3,3)=dimensional primitive translations for real space
1195 : !! in format required by wannier90
1196 : !! recip_lattice(3,3)=dimensional primitive translations for reciprocal space
1197 : !! in format required by wannier90
1198 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1199 : !! seed_name=character string for generating wannier90 filenames
1200 : !! xcart(3,natom)=atomic coordinates in bohr
1201 : !! xred(3,natom)=reduced dimensionless atomic coordinates
1202 : !!
1203 : !! OUTPUT
1204 : !! band_in(mband,nsppol) = band to take into account for wannier calculation
1205 : !! g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b
1206 : !! from k2 in the case where k1+b does not belong to the 1st BZ.
1207 : !! nband_inc(nsppol) = # of included bands
1208 : !! nntot = number of k-point neighbour
1209 : !! ovikp(nkpt,nntot)= gives nntot value of k2 (in the BZ) for each k1 (k2=k1+b mod(G))
1210 : !!
1211 : !! SIDE EFFECTS
1212 : !! (only writing, printing)
1213 : !!
1214 : !! SOURCE
1215 :
1216 19 : subroutine mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,&
1217 19 : & g1,lwanniersetup,mband,natom,nband_inc,nkpt,&
1218 19 : & nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,&
1219 19 : & proj_l,proj_m,proj_radial,proj_site,proj_s_loc, &
1220 19 : & proj_s_qaxis_loc,proj_x,proj_z,proj_zona,&
1221 19 : & real_lattice,recip_lattice,rprimd,seed_name,spinors,xcart,xred,exclude_bands)
1222 :
1223 : !Arguments---------------------------
1224 : !scalars
1225 : integer,intent(in) :: lwanniersetup,mband,natom,nkpt,nsppol
1226 : integer,intent(in) :: num_nnmax
1227 : integer,intent(out) :: nband_inc(nsppol),nntot,num_bands(nsppol),nwan(nsppol)
1228 : logical,intent(in) :: gamma_only,spinors
1229 : type(dataset_type),intent(in) :: dtset
1230 : !arrays
1231 : integer,intent(out) :: g1(3,nkpt,num_nnmax),ovikp(nkpt,num_nnmax)
1232 : integer,intent(out) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol)
1233 : real(dp),intent(in) :: real_lattice(3,3)
1234 : real(dp),intent(in) :: recip_lattice(3,3),rprimd(3,3),xred(3,natom)
1235 : real(dp),intent(out) :: proj_site(3,mband,nsppol),proj_x(3,mband,nsppol),proj_z(3,mband,nsppol)
1236 : real(dp),intent(out) :: proj_zona(mband,nsppol),xcart(3,natom)
1237 : logical,intent(out) :: band_in(mband,nsppol)
1238 : integer,intent(out) :: exclude_bands(mband,nsppol)
1239 : character(len=3),intent(out) :: atom_symbols(natom)
1240 : character(len=fnlen),intent(in) :: seed_name(nsppol),filew90_win(nsppol)
1241 : integer, optional, intent(out) :: proj_s_loc(mband)
1242 : real(dp), optional, intent(out) :: proj_s_qaxis_loc(3,mband)
1243 :
1244 : !Local variables---------------------------
1245 : !scalars
1246 : integer :: iatom,icb,ikpt,ikpt1,intot,isppol,itypat,jj,mband_,unt
1247 : real(dp) :: znucl1
1248 : character(len=2) :: symbol
1249 : character(len=500) :: msg
1250 : character(len=fnlen) :: filew90_nnkp
1251 : type(atomdata_t) :: atom
1252 : !arrays
1253 : integer :: ngkpt(3)
1254 : ! *************************************************************************
1255 :
1256 : !^^^^^^^^^^^^^^^^read wannier90.nnkp^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1257 19 : if (lwanniersetup==0) then !this part is not coded for nsppol>1
1258 0 : isppol=1
1259 0 : filew90_nnkp=trim(seed_name(isppol))//'.nnkp'
1260 0 : if (open_file(filew90_nnkp,msg,newunit=unt,form='formatted',status='old') /= 0) then
1261 0 : ABI_ERROR(msg)
1262 : end if
1263 0 : read(unt,*)
1264 0 : read(unt,*) nntot , mband_, nwan(1)
1265 0 : write(msg, '(2a,2(i0,1x))')ch10,' mlwfovlp_setup nntot, mband, nwan ', nntot,mband_,nwan(1)
1266 0 : call wrtout(std_out,msg)
1267 0 : if (mband_ /= mband) then
1268 0 : write(msg, '(4a)' )'mband_ is not equal to mband ',ch10,'Action: check ',trim(filew90_nnkp)
1269 0 : ABI_ERROR(msg)
1270 : end if
1271 0 : if(nwan(1) > mband) then
1272 0 : write(msg, '(4a)' )'nwan > mband ',ch10,'Action: check ',trim(filew90_nnkp)
1273 0 : ABI_ERROR(msg)
1274 : end if
1275 0 : if(nwan(1) == 0) then
1276 0 : write(msg, '(4a)' )'nwan = 0 ',ch10, 'Action: check ',trim(filew90_nnkp)
1277 0 : ABI_ERROR(msg)
1278 : end if
1279 0 : do ikpt=1,nkpt
1280 0 : do intot=1,nntot
1281 : ! ikpt1: k point (ikpt=ikpt1)
1282 : ! ovikp(intot,ikpt): neighbour number intot for ikpt
1283 : ! g1(1:3,intot,ikpt): non reciprocal space vector between the 2 k-points
1284 0 : read(unt,*) ikpt1,ovikp(ikpt,intot),(g1(jj,ikpt,intot),jj=1,3)
1285 0 : if(ikpt1 /= ikpt) then
1286 0 : ABI_WARNING("ikpt1 .ne ikpt : ?")
1287 : end if
1288 : end do
1289 : end do
1290 0 : close(unt)
1291 0 : write(msg, '(3a)' )ch10,trim(filew90_nnkp),'wannier90.nnkp has been read !'
1292 0 : call wrtout(std_out,msg)
1293 :
1294 0 : ABI_ERROR(' exclude bands is not given in this case (not implemented) ')
1295 :
1296 : ! ^^^^^^^^^^^^^^^^^^^^^^^ call wannier_setup begin^^^^^^^^^^^^^^^^^^^^^^^^
1297 19 : else if (lwanniersetup==1) then
1298 39 : num_bands(:)=mband
1299 : ! num_nnmax=12 !limit fixed for compact structure in wannier_setup.
1300 7927 : ovikp=0.d0
1301 : ! "When nshiftk=1, kptrlatt is initialized as a diagonal (3x3) matrix, whose diagonal
1302 : ! elements are the three values ngkpt(1:3)"
1303 19 : ngkpt(1)=dtset%kptrlatt(1,1)
1304 19 : ngkpt(2)=dtset%kptrlatt(2,2) ! have to verify that kptrlatt is diagonal
1305 19 : ngkpt(3)=dtset%kptrlatt(3,3)
1306 19 : ABI_CHECK(isdiagmat(dtset%kptrlatt), "kptrlatt must be diagonal please use ngkpt and nshiftk 1.")
1307 54 : do iatom=1,natom
1308 35 : itypat=dtset%typat(iatom)
1309 35 : znucl1=dtset%znucl(itypat)
1310 35 : call atomdata_from_znucl(atom, znucl1)
1311 35 : symbol=trim(adjustl(atom%symbol))
1312 : !write(309,*) symbol
1313 35 : atom_symbols(iatom)=symbol
1314 194 : xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+ rprimd(:,2)*xred(2,iatom)+ rprimd(:,3)*xred(3,iatom)
1315 : end do ! iatom
1316 : ! write(std_out,*) xcart; write(std_out,*) Bohr_Ang; write(std_out,*) rprimd*Bohr_Ang
1317 : ! write(std_out,*) seed_name; write(std_out,*) ngkpt; write(std_out,*) nkpt
1318 : ! write(std_out,*) mband; write(std_out,*) natom; write(std_out,*) atom_symbols
1319 19 : write(msg, '(a,a)' )ch10,' mlwfovlp_setup: calling wannier90 library subroutine wannier_setup'
1320 19 : call wrtout(std_out,msg)
1321 :
1322 : #if defined HAVE_WANNIER90
1323 39 : nwan(:)=0
1324 39 : num_bands(:)=0
1325 :
1326 39 : do isppol=1,nsppol
1327 : #ifdef HAVE_WANNIER90_V1
1328 : call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input
1329 : & ,real_lattice,recip_lattice,dtset%kptns& !input
1330 : & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input
1331 : & ,gamma_only,spinors& !input
1332 : & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output
1333 : & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output
1334 : & ,proj_m(:,isppol),proj_radial(:,isppol)& !output
1335 : & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output
1336 : & ,proj_zona(:,isppol),exclude_bands(:,isppol)) !output
1337 : #else
1338 : !WANNIER90_V2 has the 2 optional arguments
1339 39 : if (present(proj_s_loc)) then
1340 : call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input
1341 : & ,real_lattice,recip_lattice,dtset%kptns& !input
1342 : & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input
1343 : & ,gamma_only,spinors& !input
1344 : & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output
1345 : & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output
1346 : & ,proj_m(:,isppol),proj_radial(:,isppol)& !output
1347 : & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output
1348 : & ,proj_zona(:,isppol),exclude_bands(:,isppol)& !output
1349 164 : & ,proj_s_loc,proj_s_qaxis_loc) !output
1350 : else
1351 : !no proj_s_loc provided
1352 : call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input
1353 : & ,real_lattice,recip_lattice,dtset%kptns& !input
1354 : & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input
1355 : & ,gamma_only,spinors& !input
1356 : & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output
1357 : & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output
1358 : & ,proj_m(:,isppol),proj_radial(:,isppol)& !output
1359 : & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output
1360 0 : & ,proj_zona(:,isppol),exclude_bands(:,isppol)) !output
1361 : end if
1362 : #endif
1363 : end do !isppol
1364 : ! if we do not have w90, avoid complaints about unused input variables
1365 : #else
1366 : ABI_UNUSED(gamma_only)
1367 : ABI_UNUSED(real_lattice)
1368 : ABI_UNUSED(recip_lattice)
1369 : ABI_UNUSED(spinors)
1370 : #endif
1371 :
1372 39 : do isppol=1,nsppol
1373 20 : write(std_out,*) "1", nntot,nwan(isppol)
1374 20 : write(std_out,*) "2", num_bands(isppol) ! states on which wannier functions are computed
1375 20 : write(std_out,*) "3", proj_site(:,1:nwan(isppol),isppol)
1376 20 : write(std_out,*) "4", proj_l(1:nwan(isppol),isppol)
1377 20 : write(std_out,*) "5", proj_m(1:nwan(isppol),isppol)
1378 20 : write(std_out,*) "6", proj_radial(1:nwan(isppol),isppol)
1379 20 : write(std_out,*) "7", proj_z(:,1:nwan(isppol),isppol)
1380 20 : write(std_out,*) "8", proj_x(:,1:nwan(isppol),isppol)
1381 20 : write(std_out,*) "9", proj_zona(1:nwan(isppol),isppol)
1382 39 : write(std_out,*) "10", exclude_bands(:,isppol)
1383 : end do!isppol
1384 : end if ! lwanniersetup
1385 :
1386 39 : do isppol=1,nsppol
1387 214 : band_in(:,isppol)=.true.
1388 214 : do icb=1,mband
1389 214 : if(exclude_bands(icb,isppol) /= 0) band_in(exclude_bands(icb,isppol),isppol)=.false.
1390 : end do
1391 20 : nband_inc(isppol)=0
1392 233 : do icb=1, mband
1393 214 : if (band_in(icb,isppol)) nband_inc(isppol) = nband_inc(isppol)+1
1394 : end do
1395 : end do !isppol
1396 :
1397 29 : if (any(mband.gt.num_bands(:))) then
1398 10 : write(msg, '(a,a)' )ch10,' The following bands are excluded from the calculation of wannier functions:'
1399 10 : call wrtout(std_out,msg)
1400 :
1401 20 : do isppol=1,nsppol
1402 10 : if(nsppol==2) then
1403 0 : write(msg,'("For spin",i2)')isppol
1404 0 : call wrtout(std_out,msg)
1405 : end if !nsppol
1406 32 : do jj=1,mband-num_bands(isppol),10
1407 12 : write(msg,'(10i7)') exclude_bands(jj:min(jj+9,mband-num_bands(isppol)),isppol)
1408 22 : call wrtout(std_out,msg)
1409 : end do
1410 : end do !isppol
1411 : end if
1412 :
1413 39 : do isppol=1,nsppol
1414 20 : if(nsppol==2) then
1415 2 : write(msg,'("For spin",i2)')isppol
1416 2 : call wrtout(std_out,msg)
1417 : end if !nsppol
1418 20 : write(msg, '(a,i0,3a)' )ch10,nwan(isppol),' wannier functions will be computed (see ',trim(filew90_win(isppol)),')'
1419 20 : call wrtout(std_out,msg)
1420 : ! write(std_out,*) exclude_bands(icb),band_in(icb)
1421 20 : write(msg, '(a,i0,a)' )ch10,num_bands(isppol),' bands will be used to extract wannier functions'
1422 20 : call wrtout(std_out,msg)
1423 20 : if (num_bands(isppol).lt.nwan(isppol)) then
1424 : write(msg, '(4a)' )&
1425 0 : ' number of bands is lower than the number of wannier functions',ch10,&
1426 0 : ' Action : check input file and ',trim(filew90_win(isppol))
1427 0 : ABI_ERROR(msg)
1428 20 : else if (num_bands(isppol)==nwan(isppol)) then
1429 15 : write(msg, '(4a)' )ch10,&
1430 15 : ' Number of bands is equal to the number of wannier functions',ch10,&
1431 30 : ' Disentanglement is not necessary'
1432 15 : call wrtout(std_out,msg)
1433 : else if (num_bands(isppol).gt.nwan(isppol)) then
1434 5 : write(msg, '(4a)' )ch10,&
1435 5 : ' Number of bands is larger than the number of wannier functions',ch10,&
1436 10 : ' Disentanglement will be necessary'
1437 5 : call wrtout(std_out,msg)
1438 : end if
1439 20 : write(msg, '(2x,a,a,i0,1x,a)' )ch10,' Each k-point has: ', nntot,' neighbours'
1440 39 : call wrtout(std_out,msg)
1441 : end do !isppol
1442 :
1443 19 : end subroutine mlwfovlp_setup
1444 : !!***
1445 :
1446 : !!****f* m_mlwfovlp/mlwfovlp_pw
1447 : !! NAME
1448 : !! mlwfovlp_pw
1449 : !!
1450 : !! FUNCTION
1451 : !! Routine which computes PW part of overlap M_{mn}(k,b)
1452 : !! for Wannier code (www.wannier.org f90 version).
1453 : !!
1454 : !! INPUTS
1455 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
1456 : !! g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b
1457 : !! iwav(mband,nkpt,nsppol): shift for pw components in cg.
1458 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
1459 : !! mband=maximum number of bands
1460 : !! mgfft=maximum size of 1D FFTs
1461 : !! mkmem =number of k points treated by this node.
1462 : !! mpi_enreg=information about MPI parallelization
1463 : !! mpw=maximum dimensioned size of npw.
1464 : !! nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
1465 : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
1466 : !! nkpt=number of k points.
1467 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1468 : !! nspinor=number of spinorial components of the wavefunctions
1469 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1470 : !! ovikp(nkpt,nntot)= gives nntot value of k2 (in the BZ) for each k1 (k2=k1+b mod(G))
1471 : !! seed_name= seed_name of files containing cg for all k-points to be used with MPI
1472 : !!
1473 : !! OUTPUT
1474 : !! cm1(2,mband,mband,nntot,nkpt,nsppol): overlap <u_(nk1)|u_(mk1+b)>.
1475 : !!
1476 : !! SIDE EFFECTS
1477 : !! (only writing, printing)
1478 : !!
1479 : !! SOURCE
1480 :
1481 19 : subroutine mlwfovlp_pw(mywfc,cm1,g1,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
1482 19 : npwarr,nspinor,nsppol,ovikp)
1483 :
1484 : !Arguments ------------------------------------
1485 : !scalars
1486 : integer,intent(in) :: mband,mkmem,mpw,nfft,nkpt,nntot
1487 : integer,intent(in) :: nspinor,nsppol
1488 : ! character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI
1489 : class(abstract_wf) :: mywfc
1490 : type(MPI_type),intent(in) :: mpi_enreg
1491 : !arrays
1492 : integer,intent(in) :: g1(3,nkpt,nntot),kg(3,mpw*mkmem),ngfft(18),npwarr(nkpt)
1493 : ! integer,intent(in) :: iwav(mband,nkpt,nsppol)
1494 : integer,intent(in) :: ovikp(nkpt,nntot)
1495 : ! real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
1496 : real(dp),intent(out) :: cm1(2,mband,mband,nntot,nkpt,nsppol)
1497 :
1498 : !Local variables-------------------------------
1499 : !scalars
1500 : integer :: iband1,iband2,ierr,ig,ig1,ig1b,ig2,ig2b
1501 : integer :: ig3,ig3b,igk1,igk2,ikg,ikpt,ikpt1,ikpt2,imntot,index,intot
1502 : integer :: ispinor,isppol,me,n1,n2,n3,npoint,npoint2,npw_k,npw_k2
1503 : integer :: nprocs,comm
1504 19 : integer,allocatable :: indpwk(:,:),kg_k(:,:), invpwk(:,:)
1505 : character(len=500) :: msg
1506 : logical:: lfile
1507 19 : real(dp),allocatable :: cg_read(:,:) !to be used in case of MPI
1508 : !************************************************************************
1509 :
1510 19 : write(msg, '(a,a)' ) ch10, '** mlwfovlp_pw : compute pw part of overlap'
1511 19 : call wrtout(std_out, msg)
1512 :
1513 : !initialize flags
1514 19 : lfile=.false.
1515 : !mpi initialization
1516 19 : comm=MPI_enreg%comm_cell
1517 19 : nprocs=xmpi_comm_size(comm)
1518 19 : me=MPI_enreg%me_kpt
1519 :
1520 19 : if(nprocs>1) then
1521 0 : ABI_MALLOC(cg_read,(2,nspinor*mpw*mband))
1522 : end if
1523 :
1524 :
1525 : !****************compute intermediate quantities (index, shifts) ******
1526 : !------------compute index for g points--------------------------------
1527 : !ig is a plane waves which belongs to the sphere ecut for ikpt (they
1528 : !are npwarr(ikpt))
1529 : !npoint is the position in the grid of planes waves
1530 : !(they are nfft)
1531 : !indpwk is a application ig-> npoint
1532 : !invpwk is not an application (some npoint have no ig corresponding)
1533 : !cg are ordered with npw_k !
1534 : !----------------------------------------------------------------------
1535 : !------------compute index for g points--------------------------------
1536 : !----------------------------------------------------------------------
1537 19 : write(msg, '(a,a)' ) ch10,' first compute index for g-points'
1538 19 : call wrtout(std_out, msg)
1539 :
1540 : ! Allocations
1541 57 : ABI_MALLOC(kg_k,(3,mpw))
1542 76 : ABI_MALLOC(indpwk,(nkpt,mpw))
1543 76 : ABI_MALLOC(invpwk,(nkpt,nfft))
1544 :
1545 19 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1546 7395000 : invpwk=0
1547 205756 : indpwk=0
1548 54231 : kg_k=0
1549 :
1550 38 : do isppol=1,1 !invpwk is not spin dependent so we just do it once
1551 19 : ikg=0
1552 678 : do ikpt=1,nkpt
1553 : ! MPI:cycle over k-points not treated by this node
1554 640 : if (nprocs>1 ) then !sometimes we can have just one processor
1555 : !print *, "MPI_enreg%proc_distrb(ikpt,1,isppol)", ikpt, isppol , MPI_enreg%proc_distrb(ikpt,1,isppol)
1556 0 : if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-MPI_enreg%me) /= 0) CYCLE
1557 : end if
1558 :
1559 : ! write(std_out,*)'me',me,'ikpt',ikpt,'isppol',isppol
1560 6992088 : do npoint=1,nfft
1561 6992088 : if(invpwk(ikpt,npoint)/=0 )then
1562 0 : write(std_out,*) "error0 , invpwk is overwritten"
1563 0 : write(std_out,*) ikpt,npoint
1564 0 : ABI_ERROR("Aborting now")
1565 : end if
1566 : end do
1567 640 : npw_k=npwarr(ikpt)
1568 : ! write(std_out,*) ikpt,npw_k,nfft
1569 741452 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
1570 185843 : do ig=1,npw_k
1571 185203 : if(ig.gt.mpw) then
1572 0 : write(std_out,*)"error ig",ig,"greater than mpw ",mpw
1573 0 : ABI_ERROR("Aborting now")
1574 : end if
1575 185203 : if(indpwk(ikpt,ig)/=0) then
1576 0 : write(std_out,*) "error, indpwk is overwritten"
1577 0 : write(std_out,*) ikpt,ig,indpwk(ikpt,ig)
1578 0 : ABI_ERROR("Aborting now")
1579 : end if
1580 185203 : ig1=modulo(kg_k(1,ig),n1)
1581 185203 : ig2=modulo(kg_k(2,ig),n2)
1582 185203 : ig3=modulo(kg_k(3,ig),n3)
1583 185203 : indpwk(ikpt,ig)=ig1+1+n1*(ig2+n2*ig3)
1584 185203 : npoint=indpwk(ikpt,ig)
1585 185203 : if(npoint.gt.nfft) then
1586 0 : ABI_ERROR("error npoint")
1587 : end if
1588 : ! write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
1589 185203 : if(invpwk(ikpt,npoint)/=0) then
1590 0 : write(std_out,*) "error, invpwk is overwritten"
1591 0 : write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
1592 0 : ABI_ERROR("Aborting now")
1593 : end if
1594 185843 : invpwk(ikpt,npoint)=ig
1595 : ! write(std_out,*)'ikpt,npoint,invpwk',ikpt,npoint,invpwk(ikpt,npoint)
1596 : ! if(ikpt.eq.1) write(std_out,*) "ig npoint",ig, npoint
1597 : ! write(std_out,*) "ikpt ig npoint",ikpt,ig, npoint
1598 : end do
1599 19 : ikg=ikg+npw_k
1600 :
1601 : end do !ikpt
1602 : end do !isppol
1603 : !write(std_out,*) "index for g points has been computed"
1604 :
1605 19 : call xmpi_sum(invpwk,comm,ierr)
1606 :
1607 : !----------------------------------------------------------------------
1608 : !------------test invpwk-----------------------------------------------
1609 : !----------------------------------------------------------------------
1610 : !write(std_out,*) "TEST INVPWK"
1611 : !ikpt=3
1612 : !isppol=1
1613 : !do ig=1,npwarr(ikpt)
1614 : !npoint=indpwk(ikpt,ig)
1615 : !write(std_out,*) "ig npoint ",ig, npoint
1616 : !write(std_out,*) "ig npoint inv",invpwk(ikpt,npoint),npoint
1617 : !end do
1618 : !do ig3=1,n3
1619 : !do ig2=1,n2
1620 : !do ig1=1,n1
1621 : !npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
1622 : !ig=invpwk(ikpt,npoint)
1623 : !! if(ig/=0) write(std_out,*) "ig npoint",ig, npoint
1624 : !end do
1625 : !end do
1626 : !end do
1627 :
1628 19 : ABI_FREE(kg_k)
1629 19 : ABI_FREE(indpwk)
1630 :
1631 : !***********************************************************************
1632 : !**calculate overlap M_{mn}(k,b)=<\Psi_{k,m}|e^{-ibr}|\Psi_{k+b,n}>*****
1633 : !***********************************************************************
1634 19 : write(msg, '(a,a)' ) ch10,' mlwfovlp_pw : compute overlaps '
1635 19 : call wrtout(std_out, msg)
1636 19 : write(msg, '(a,a)' ) ch10," nkpt nntot mband "
1637 19 : call wrtout(std_out, msg)
1638 19 : write(msg, '(i6,2x,i6,2x,i6,2x,i6)' ) nkpt,nntot,mband
1639 19 : call wrtout(std_out, msg)
1640 3614231 : cm1=zero
1641 19 : write(msg, '(a)' ) ' '
1642 19 : call wrtout(std_out, msg)
1643 :
1644 39 : do isppol=1,nsppol
1645 20 : imntot=0
1646 743 : do ikpt1=1,nkpt
1647 : ! MPI:cycle over k-points not treated by this node
1648 704 : if (nprocs>1) then
1649 0 : if (ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-me) /= 0) CYCLE
1650 : end if
1651 704 : write(msg, '(a,i0,a,i0,a,i0)' ) ' Processor: ',me,' computes k-point: ',ikpt1,', and spin: ',isppol
1652 704 : call wrtout(std_out, msg)
1653 :
1654 7108 : do intot=1,nntot
1655 6384 : lfile=.false. !flag to know if this kpt will be read from a file, see below
1656 6384 : imntot=imntot+1
1657 6384 : ikpt2= ovikp(ikpt1,intot)
1658 : ! write(std_out,*)'me',me,'ikpt1',ikpt1,'ikpt2',ikpt2,'intot',intot,'isppol',isppol
1659 :
1660 : !
1661 : ! MPI: if ikpt2 not found in this processor then
1662 : ! read info from an unformatted file
1663 : ! TODO: also get MPI mapping to retrieve who has this wf k-point
1664 : !
1665 : ! if (nprocs>1) then
1666 : ! if ( ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me) /=0) then
1667 : ! lfile=.true.
1668 : ! write(cg_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol
1669 : ! iunit=1000+ikpt2+ikpt2*(isppol-1)
1670 : ! npw_k2=npwarr(ikpt2)
1671 : ! open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios)
1672 : ! if(ios /= 0) then
1673 : ! write(msg,*) " mlwfovlp_pw: file",trim(cg_file), "not found"
1674 : ! ABI_ERROR(msg)
1675 : ! end if
1676 : ! !
1677 : ! do iband2=1,mband
1678 : ! do ipw=1,npw_k2*nspinor
1679 : ! index=ipw+(iband2-1)*npw_k2*nspinor
1680 : ! read(iunit) (cg_read(ii,index),ii=1,2)
1681 : ! ! if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
1682 : ! ! if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
1683 : ! end do
1684 : ! end do
1685 : ! close(iunit)
1686 : ! end if
1687 : ! end if
1688 :
1689 6384 : if(nprocs>1) then
1690 : !call mywfc%read_cg(cg_read, ikpt2)
1691 : !call mywfc%read_cg( ikpt2, isppol, cg_read)
1692 :
1693 0 : if (ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me) /= 0) then
1694 0 : lfile=.true.
1695 0 : call mywfc%load_cg(ikpt2, isppol, cg_read)
1696 : endif
1697 : end if
1698 :
1699 6384 : npw_k=npwarr(ikpt1)
1700 6384 : npw_k2=npwarr(ikpt2)
1701 133872 : do ig3=1,n3
1702 2849904 : do ig2=1,n2
1703 66270784 : do ig1=1,n1
1704 : ! write(std_out,*) isppol,ikpt1,iband1,iband2,intot
1705 63427264 : npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
1706 63427264 : if(npoint.gt.nfft) then
1707 0 : ABI_ERROR("error npoin Aborting now")
1708 : end if
1709 63427264 : ig1b=ig1+g1(1,ikpt1,intot)
1710 63427264 : ig2b=ig2+g1(2,ikpt1,intot)
1711 63427264 : ig3b=ig3+g1(3,ikpt1,intot)
1712 : ! write(std_out,*) ig1,ig2,ig3; write(std_out,*) ig1b,ig2b,ig3b
1713 63427264 : if(ig1b.lt.1) ig1b=ig1b+n1
1714 63427264 : if(ig2b.lt.1) ig2b=ig2b+n2
1715 63427264 : if(ig3b.lt.1) ig3b=ig3b+n3
1716 63427264 : if(ig1b.gt.n1) ig1b=ig1b-n1
1717 63427264 : if(ig2b.gt.n2) ig2b=ig2b-n2
1718 63427264 : if(ig3b.gt.n3) ig3b=ig3b-n3
1719 63427264 : npoint2=ig1b+(ig2b-1)*n1+(ig3b-1)*n2*n1
1720 63427264 : if(npoint2.gt.nfft) then
1721 0 : ABI_ERROR("error npoint c")
1722 : end if
1723 63427264 : igk1=invpwk(ikpt1,npoint)
1724 63427264 : igk2=invpwk(ikpt2,npoint2)
1725 :
1726 : ! if(intot==10) write(std_out,*)'Before igk1 and igk2',ikpt1,ikpt2,isppol
1727 :
1728 66144000 : if(igk1/=0.and.igk2/=0) then
1729 17453252 : do iband2=1,mband
1730 216841732 : do iband1=1,mband
1731 468974288 : do ispinor=1,nspinor
1732 253737536 : if(lfile) index=ispinor + nspinor*(igk2-1) + nspinor*npw_k2*(iband2-1) !In case of MPI, see below
1733 : ! TODO : Check if the index in the cg_elems are correct.
1734 : !
1735 : ! If MPI sometimes the info was read from an unformatted file
1736 : ! If that is the case lfile==.true.
1737 : !
1738 : ! TODO: this filter should be outside, not inside 1000 loops!!!
1739 199388480 : if(lfile) then
1740 : cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ &
1741 : & mywfc%cg_elem(1, igk1, ispinor, iband1, ikpt1, isppol) *cg_read(1,index)&
1742 0 : & + mywfc%cg_elem(2, igk1, ispinor, iband1, ikpt1, isppol)*cg_read(2,index)
1743 : cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ &
1744 : & mywfc%cg_elem(1, igk1, ispinor, iband1, ikpt1, isppol)*cg_read(2,index)&
1745 0 : &- mywfc%cg_elem(2, igk1, ispinor, iband1, ikpt1, isppol)*cg_read(1,index)
1746 : else
1747 : ! TODO: Here it is very inefficient.
1748 : ! Could be replaced with the fftbox and dotproduct.
1749 : ! cgtk_rotate. sphere.
1750 : !
1751 : cm1(1,iband1,iband2,intot,ikpt1,isppol)=&
1752 : & cm1(1,iband1,iband2,intot,ikpt1,isppol) &
1753 : & + mywfc%cg_elem(1, igk1, ispinor,iband1, ikpt1, isppol) &
1754 : & *mywfc%cg_elem(1, igk2, ispinor,iband2, ikpt2, isppol) &
1755 : & +mywfc%cg_elem(2, igk1, ispinor,iband1, ikpt1, isppol) &
1756 253737536 : & *mywfc%cg_elem(2, igk2, ispinor,iband2, ikpt2, isppol)
1757 : cm1(2,iband1,iband2,intot,ikpt1,isppol)= &
1758 : & cm1(2,iband1,iband2,intot,ikpt1,isppol) &
1759 : & + mywfc%cg_elem(1, igk1, ispinor,iband1, ikpt1, isppol) &
1760 : & *mywfc%cg_elem( 2, igk2, ispinor,iband2, ikpt2, isppol) &
1761 : & -mywfc%cg_elem( 2, igk1, ispinor,iband1, ikpt1, isppol) &
1762 253737536 : & *mywfc%cg_elem( 1, igk2, ispinor,iband2, ikpt2, isppol)
1763 : end if
1764 : end do !ispinor
1765 : end do ! iband1
1766 : end do ! iband2
1767 : end if
1768 : end do ! ig1
1769 : end do ! ig2
1770 : end do ! ig3
1771 : end do ! intot
1772 : end do ! ikpt1
1773 : end do ! isppol
1774 :
1775 19 : ABI_FREE(invpwk)
1776 19 : ABI_SFREE(cg_read)
1777 :
1778 38 : end subroutine mlwfovlp_pw
1779 : !!***
1780 :
1781 : !!****f* m_mlwfovlp/mlwfovlp_proj
1782 : !! NAME
1783 : !! mlwfovlp_proj
1784 : !!
1785 : !! FUNCTION
1786 : !! Routine which computes projection A_{mn}(k) for Wannier code (www.wannier.org f90 version).
1787 : !!
1788 : !! INPUTS
1789 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
1790 : !! cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
1791 : !! and each |p_lmn> non-local projector
1792 : !! dtset <type(dataset_type)>=all input variables for this dataset
1793 : !! filew90_win = secondary input file for wannier90 (WAS NOT USED IN v6.7.1 - so has been temporarily removed)
1794 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
1795 : !! lproj= flag 0: no projections, 1: random projections,
1796 : !! 2: projections on atomic orbitals
1797 : !! 3: projections on projectors
1798 : !! mband=maximum number of bands
1799 : !! mkmem =number of k points treated by this node.
1800 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1801 : !! mpi_enreg=information about MPI parallelization
1802 : !! mpw=maximum dimensioned size of npw.
1803 : !! natom=number of atoms in cell.
1804 : !! nattyp(ntypat)= # atoms of each type.
1805 : !! nkpt=number of k points.
1806 : !! nspinor=number of spinorial components of the wavefunctions
1807 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1808 : !! ntypat=number of types of atoms in unit cell.
1809 : !! num_bands=number of bands actually used to construct the wannier function
1810 : !! nwan= number of wannier fonctions (read in wannier90.win).
1811 : !! proj_l(mband)= angular part of the projection function (quantum number l)
1812 : !! proj_m(mband)= angular part of the projection function (quantum number m)
1813 : !! proj_radial(mband)= radial part of the projection.
1814 : !! proj_site(3,mband)= site of the projection.
1815 : !! proj_x(3,mband)= x axis for the projection.
1816 : !! proj_z(3,mband)= z axis for the projection.
1817 : !! proj_zona(mband)= extension of the radial part.
1818 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
1819 : !!
1820 : !! OUTPUT
1821 : !! A_matrix(num_bands,nwan,nkpt,nsppol)= Matrix of projections needed by wannier_run
1822 : !! ( also wannier90random.amn is written)
1823 : !!
1824 : !! SIDE EFFECTS
1825 : !! (only writing, printing)
1826 : !!
1827 : !! SOURCE
1828 :
1829 19 : subroutine mlwfovlp_proj(A_matrix,band_in,mywfc, dtset,gprimd,just_augmentation,kg,&
1830 19 : lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,nattyp,&
1831 19 : nkpt,npwarr,nspinor,&
1832 38 : nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,proj_radial,&
1833 19 : proj_site,proj_x,proj_z,proj_zona,psps,ucvol)
1834 :
1835 : !Arguments ------------------------------------
1836 : !scalars
1837 : complex(dp),parameter :: c1=(1._dp,0._dp)
1838 : integer,intent(in) :: lproj,max_num_bands,mband,mkmem,mpw,mwan,natom,nkpt,nspinor,nsppol
1839 : integer,intent(in) :: ntypat
1840 : type(MPI_type),intent(in) :: mpi_enreg
1841 : type(dataset_type),intent(in) :: dtset
1842 : type(pseudopotential_type),intent(in) :: psps
1843 : !arrays
1844 : integer ::nattyp(ntypat)
1845 : integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt),num_bands(nsppol),nwan(nsppol),proj_l(mband,nsppol)
1846 : integer,intent(in) :: proj_m(mband,nsppol)
1847 : integer,intent(inout)::proj_radial(mband,nsppol)
1848 : !real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
1849 : real(dp),intent(in) :: gprimd(3,3),proj_site(3,mband,nsppol)
1850 : real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol)
1851 : complex(dp),intent(out) :: A_matrix(max_num_bands,mwan,nkpt,nsppol)
1852 : !character(len=fnlen),intent(in) :: filew90_win(nsppol)
1853 : logical,intent(in) :: band_in(mband,nsppol)
1854 : logical,intent(in)::just_augmentation(mwan,nsppol)
1855 : !type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
1856 : !type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1857 : type(pawtab_type),intent(in) :: pawtab(:)
1858 : class(abstract_wf), intent(inout) :: mywfc
1859 :
1860 : !Local variables-------------------------------
1861 : !scalars
1862 : integer :: iatom,iatprjn,iband,iband1,iband2,ibg,icat,icg,icg_shift
1863 : integer :: idum,ikg,ikpt,ilmn,ipw,iproj
1864 : integer :: ispinor,isppol,itypat,iwan,jband,jj1,libprjn
1865 : integer :: lmn_size,natprjn,nband_k,nbprjn,npw_k, sumtmp
1866 : integer :: max_lmax,max_lmax2,mproj,nprocs,comm,rank, idx
1867 : real(dp),parameter :: qtol=2.0d-8
1868 : real(dp) :: arg,norm_error,norm_error_bar
1869 : real(dp) :: ucvol,x1,x2,xnorm,xnormb,xx,yy,zz
1870 38 : complex(dp) :: amn_tmp(nspinor)
1871 : complex(dp) :: cstr_fact
1872 : character(len=500) :: msg
1873 : !arrays
1874 38 : integer :: kg_k(3,mpw),lmax(nsppol),lmax2(nsppol),nproj(nsppol)
1875 19 : integer,allocatable :: lprjn(:),npprjn(:)
1876 : real(dp) :: kpg(3),kpt(3)
1877 19 : real(dp),allocatable :: amn(:,:,:,:,:),amn2(:,:,:,:,:,:,:)
1878 19 : real(dp),allocatable :: gsum2(:),kpg2(:),radial(:)
1879 19 : complex(dp),allocatable :: gf(:,:),gft_lm(:), ylmc_fac(:,:,:),ylmcp(:)
1880 : !Tables 3.1 & 3.2, User guide
1881 : integer,parameter :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/)
1882 : ! integer,parameter :: mtransfo(0:3,7)=&
1883 : !& reshape((/1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,-2,-1,2,1,0,0,0,-1,1,2,-2,-3,3/),(/4,7/))
1884 : !************************************************************************
1885 :
1886 : !mpi initialization
1887 19 : comm=MPI_enreg%comm_cell
1888 19 : nprocs=xmpi_comm_size(comm)
1889 19 : rank=MPI_enreg%me_kpt
1890 :
1891 : ! Check input variables
1892 19 : if ((lproj/=1).and.(lproj/=2).and.(lproj/=5)) then
1893 0 : write(msg, '(3a)' )' Value of lproj no allowed ',ch10,' Action: change lproj.'
1894 0 : ABI_ERROR(msg)
1895 : end if
1896 :
1897 19 : write(msg, '(a,a)' )ch10,'** mlwfovlp_proj: compute A_matrix of initial guess for wannier functions'
1898 19 : call wrtout(std_out,msg)
1899 :
1900 : !Initialize to 0.d0
1901 65031 : A_matrix(:,:,:,:)=cmplx(0.d0,0.d0)
1902 :
1903 : !********************* Write Random projectors
1904 19 : if(lproj==1) then
1905 6 : idum=123456
1906 : ! Compute random projections
1907 36 : ABI_MALLOC(amn,(2,mband,mwan,nkpt,nsppol))
1908 3324 : amn=zero
1909 12 : do isppol=1,nsppol
1910 60 : do ikpt=1,nkpt
1911 : ! MPI: cycle over kpts not treated by this node
1912 48 : if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
1913 : ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') ikpt,rank
1914 :
1915 310 : do iband1=1,mband
1916 256 : xnormb=0.d0
1917 1280 : do iband2=1,nwan(isppol)
1918 1024 : x1=uniformrandom(idum)
1919 1024 : x2=uniformrandom(idum)
1920 1024 : xnorm=sqrt(x1**2+x2**2)
1921 1024 : xnormb=xnormb+xnorm
1922 1024 : amn(1,iband1,iband2,ikpt,isppol)=x1
1923 1280 : amn(2,iband1,iband2,ikpt,isppol)=x2
1924 : end do
1925 1328 : do iband2=1,nwan(isppol)
1926 1024 : amn(1,iband1,iband2,ikpt,isppol)=amn(1,iband1,iband2,ikpt,isppol)/xnormb
1927 1280 : amn(2,iband1,iband2,ikpt,isppol)=amn(2,iband1,iband2,ikpt,isppol)/xnormb
1928 : end do !iband2
1929 : end do !iband1
1930 : end do !ikpt
1931 : end do !isppol
1932 12 : do isppol=1,nsppol
1933 60 : do ikpt=1,nkpt
1934 : ! MPI: cycle over kpts not treated by this node
1935 48 : if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
1936 246 : do iband2=1,nwan(isppol)
1937 : jband=0
1938 1264 : do iband1=1,mband
1939 1216 : if(band_in(iband1,isppol)) then
1940 768 : jband=jband+1
1941 768 : if(jband.gt.num_bands(isppol)) then
1942 0 : ABI_ERROR('Value of jband is above num_bands ')
1943 : end if
1944 768 : A_matrix(jband,iband2,ikpt,isppol)=cmplx(amn(1,iband1,iband2,ikpt,isppol),amn(2,iband1,iband2,ikpt,isppol))
1945 : end if
1946 : end do !iband1
1947 : end do !iband2
1948 : end do !ikpt
1949 : end do !isppol
1950 6 : ABI_FREE(amn)
1951 : end if
1952 :
1953 : !********************* Projection on atomic orbitals based on .win file
1954 19 : if( lproj==2) then !based on .win file
1955 27 : nproj(:)=nwan(:)/nspinor !if spinors, then the number of projections are
1956 27 : mproj=maxval(nproj(:))
1957 : ! half the total of wannier functions. obtain lmax and lmax2
1958 27 : lmax(:)=0
1959 27 : lmax2(:)=0
1960 :
1961 27 : do isppol=1,nsppol
1962 89 : do iproj=1,nproj(isppol)
1963 89 : lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iproj,isppol)))
1964 : end do !iproj
1965 27 : lmax2(isppol)=(lmax(isppol)+1)**2
1966 : end do !isppol
1967 : max_lmax=maxval(lmax(:))
1968 27 : max_lmax2=maxval(lmax2(:))
1969 : ! Allocate arrays
1970 65 : ABI_MALLOC(ylmc_fac,(max_lmax2,mproj,nsppol))
1971 :
1972 : ! get ylmfac, factor used for rotations and hybrid orbitals
1973 27 : do isppol=1,nsppol
1974 : !print *, "nproj", nproj(isppol)
1975 : !print *, "isppol", isppol
1976 : !print *, ylmc_fac(1:lmax2(isppol),1:nproj(isppol),isppol)
1977 : !print *, "lmax, lmax2: ", lmax(isppol),lmax2(isppol)
1978 : !print *, nproj(isppol),proj_l(:,isppol),proj_m(:,isppol),proj_x(:,:,isppol)
1979 : !print *, proj_z(:,:,isppol)
1980 : call mlwfovlp_ylmfac(ylmc_fac(1:lmax2(isppol),1:nproj(isppol),isppol),lmax(isppol),lmax2(isppol),&
1981 27 : & mband,nproj(isppol),proj_l(:,isppol),proj_m(:,isppol),proj_x(:,:,isppol),proj_z(:,:,isppol))
1982 : end do
1983 :
1984 13 : norm_error=zero
1985 13 : norm_error_bar=zero
1986 13 : icg=0
1987 :
1988 27 : do isppol=1,nsppol
1989 : ! Allocate arrays
1990 : ! this has to be done this way because the variable icg changes at the end of the
1991 : ! cycle. We cannot just skip the whole cycle.
1992 56 : ABI_MALLOC(gf,(mpw,nproj(isppol)))
1993 42 : ABI_MALLOC(gft_lm,(lmax2(isppol)))
1994 42 : ABI_MALLOC(gsum2,(nproj(isppol)))
1995 42 : ABI_MALLOC(kpg2,(mpw))
1996 42 : ABI_MALLOC(radial,(lmax2(isppol)))
1997 28 : ABI_MALLOC(ylmcp,(lmax2(isppol)))
1998 14 : ikg=0
1999 670 : do ikpt=1, nkpt
2000 : ! MPI: cycle over kpts not treated by this node
2001 656 : if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
2002 656 : write(msg, '(a,i6,a,2(i0,1x))' )' processor',rank,' will compute k-point,spin=',ikpt,isppol
2003 656 : call wrtout(std_out, msg)
2004 :
2005 : ! Initialize variables
2006 656 : npw_k=npwarr(ikpt)
2007 4496 : gsum2(:)=0.d0
2008 929680 : gf(:,:) = (0.d0,0.d0)
2009 2624 : kpt(:)=dtset%kptns(:,ikpt)
2010 663088 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2011 :
2012 166264 : do ipw=1, npw_k
2013 165608 : kpg(1)= (kpt(1) + real(kg_k(1,ipw),dp)) !k+G
2014 165608 : kpg(2)= (kpt(2) + real(kg_k(2,ipw),dp))
2015 165608 : kpg(3)= (kpt(3) + real(kg_k(3,ipw),dp))
2016 : ! Calculate modulus of k+G
2017 165608 : xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3)
2018 165608 : yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3)
2019 165608 : zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3)
2020 165608 : kpg2(ipw)= two_pi*sqrt(xx**2+yy**2+zz**2)
2021 : ! Complex Y_lm for k+G
2022 165608 : if(lmax(isppol)==0) then
2023 30814 : ylmcp(1)=c1/sqrt(four_pi)
2024 : else
2025 134794 : call ylm_cmplx(lmax(isppol),ylmcp,xx,yy,zz)
2026 : end if
2027 :
2028 1042449 : do iproj=1,nproj(isppol)
2029 : ! In PAW, we can use proj_radial > 4 to indicate that we just want the in-sphere contribution
2030 876185 : if( psps%usepaw==1) then
2031 406749 : if( just_augmentation(iproj,isppol)) cycle
2032 : end if
2033 : ! obtain radial part
2034 : call mlwfovlp_radial(proj_zona(iproj,isppol),lmax(isppol),lmax2(isppol), &
2035 876185 : radial,proj_radial(iproj,isppol),kpg2(ipw))
2036 : ! scale complex representation of projector orbital with radial functions of appropriate l
2037 5601682 : gft_lm(:)=radial(:)*ylmc_fac(1:lmax2(isppol),iproj,isppol)
2038 : ! complex structure factor for projector orbital position
2039 : arg = ( kpg(1)*proj_site(1,iproj,isppol) + &
2040 : kpg(2)*proj_site(2,iproj,isppol) + &
2041 876185 : kpg(3)*proj_site(3,iproj,isppol) ) * 2*pi
2042 876185 : cstr_fact = cmplx(cos(arg), -sin(arg) )
2043 :
2044 : ! obtain guiding functions
2045 5601682 : gf(ipw,iproj)=cstr_fact*dot_product(ylmcp,gft_lm)
2046 1041793 : gsum2(iproj)=gsum2(iproj)+real(gf(ipw,iproj))**2+aimag(gf(ipw,iproj))**2
2047 : end do !iproj
2048 : end do !ipw
2049 :
2050 4496 : do iproj=1,nproj(isppol)
2051 : ! In PAW, we can use proj_radial > 4 to indicate that we just want the in-sphere contribution
2052 3840 : if(psps%usepaw==1 ) then
2053 2272 : if (just_augmentation(iproj,isppol)) cycle
2054 : end if
2055 3840 : gsum2(iproj)=16._dp*pi**2*gsum2(iproj)/ucvol
2056 929024 : gf(:,iproj)=gf(:,iproj)/sqrt(gsum2(iproj))
2057 3840 : norm_error=max(abs(gsum2(iproj)-one),norm_error)
2058 4496 : norm_error_bar=norm_error_bar+(gsum2(iproj)-one)**2
2059 : end do !iproj
2060 : !
2061 : ! Guiding functions are computed.
2062 : ! compute overlaps of gaussian projectors and wave functions
2063 4496 : do iproj=1,nproj(isppol)
2064 : !
2065 : ! In PAW, we can use proj_radial > 4 to indicate that we just
2066 : ! want the in-sphere contribution
2067 : !
2068 3840 : if(psps%usepaw==1 ) then
2069 2272 : if ( just_augmentation(iproj,isppol)) cycle
2070 : end if
2071 : !
2072 3840 : jband=0
2073 56528 : do iband=1,mband
2074 55872 : if(band_in(iband,isppol)) then
2075 44800 : icg_shift=npw_k*nspinor*(iband-1)+icg
2076 44800 : jband=jband+1
2077 103680 : amn_tmp(:)=cmplx(0.d0,0.d0)
2078 103680 : do ispinor=1,nspinor
2079 10012872 : do ipw=1,npw_k
2080 : !
2081 : ! The case of spinors is tricky, we have nproj = nwan/2
2082 : ! so we project to spin up and spin down separately, to have at
2083 : ! the end an amn matrix with nwan projections.
2084 9909192 : idx=ipw*nspinor - (nspinor-ispinor)
2085 58880 : select type(mywfc)
2086 : type is (cg_cprj)
2087 8152856 : amn_tmp(ispinor)=amn_tmp(ispinor)+gf(ipw,iproj)*cmplx(mywfc%cg(1,idx+icg_shift),-mywfc%cg(2,idx+icg_shift))
2088 : type is (wfd_wf)
2089 1756336 : amn_tmp(ispinor)=amn_tmp(ispinor)+gf(ipw,iproj)*conjg(mywfc%cg_elem_complex(ipw, ispinor, iband, ikpt, isppol))
2090 : end select
2091 : end do !ipw
2092 : end do !ispinor
2093 103680 : do ispinor=1,nspinor
2094 58880 : iwan=(iproj*nspinor)- (nspinor-ispinor)
2095 103680 : A_matrix(jband,iwan,ikpt,isppol)=amn_tmp(ispinor)
2096 : end do
2097 : end if !band_in
2098 : end do !iband
2099 : end do !iproj
2100 656 : icg=icg+npw_k*nspinor*mband
2101 670 : ikg=ikg+npw_k
2102 : end do !ikpt
2103 : ! Deallocations
2104 14 : ABI_FREE(gf)
2105 14 : ABI_FREE(gft_lm)
2106 14 : ABI_FREE(gsum2)
2107 14 : ABI_FREE(kpg2)
2108 14 : ABI_FREE(radial)
2109 27 : ABI_FREE(ylmcp)
2110 : end do !isppol
2111 : !
2112 : ! if(isppol==1) then
2113 : ! norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)),dp))
2114 : ! else
2115 : ! norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)+nwan(2)),dp))
2116 : ! end if
2117 : ! if(norm_error>0.05_dp) then
2118 : ! write(msg, '(6a,f6.3,a,f6.3,12a)' )ch10,&
2119 : ! & ' mlwfovlp_proj : WARNING',ch10,&
2120 : ! & ' normalization error for wannier projectors',ch10,&
2121 : ! & ' is',norm_error_bar,' (average) and',norm_error,' (max).',ch10,&
2122 : ! & ' this may indicate more cell-to-cell overlap of the radial functions',ch10,&
2123 : ! & ' than you want.',ch10,&
2124 : ! & ' Action : modify zona (inverse range of radial functions)',ch10,&
2125 : ! ' under "begin projectors" in ',trim(filew90_win),' file',ch10
2126 : ! call wrtout(std_out,msg)
2127 : ! end if
2128 : !
2129 13 : ABI_FREE(ylmc_fac)
2130 : end if !lproj==2
2131 :
2132 :
2133 : !*************** computes projection from PROJECTORS ********************
2134 19 : if(lproj==3) then !! if LPROJPRJ
2135 : ! ----- set values for projections --------------------- ! INPUT
2136 : ! nbprjn:number of different l-values for projectors
2137 : ! lprjn: value of l for each projectors par ordre croissant
2138 : ! npprjn: number of projectors for each lprjn
2139 0 : natprjn=1 ! atoms with wannier functions are first
2140 : if(natprjn/=1) then ! in this case lprjn should depend on iatprjn
2141 : ABI_ERROR("natprjn/=1")
2142 : end if
2143 0 : nbprjn=2
2144 0 : ABI_MALLOC(lprjn,(nbprjn))
2145 0 : lprjn(1)=0
2146 0 : lprjn(2)=1
2147 0 : ABI_MALLOC(npprjn,(0:lprjn(nbprjn)))
2148 0 : npprjn(0)=1
2149 0 : npprjn(1)=1
2150 : ! --- test coherence of nbprjn and nwan
2151 0 : sumtmp=0
2152 0 : do iatprjn=1,natprjn
2153 0 : do libprjn=0,lprjn(nbprjn)
2154 0 : sumtmp=sumtmp+(2*libprjn+1)*npprjn(libprjn)
2155 : end do
2156 : end do
2157 0 : if(sumtmp/=nwan(1)) then
2158 0 : write(std_out,*) "Number of Wannier orbitals is not equal to number of projections"
2159 0 : write(std_out,*) "Action: check values of lprjn,npprjn % nwan"
2160 0 : write(std_out,*) "nwan, sumtmp=",nwan,sumtmp
2161 0 : ABI_ERROR("Aborting now")
2162 : end if
2163 : ! --- end test of coherence
2164 0 : ABI_MALLOC(amn2,(2,natom,nsppol,nkpt,mband,nspinor,nwan(1)))
2165 0 : if(psps%usepaw==1) then
2166 0 : amn2=zero
2167 0 : ibg=0
2168 0 : do isppol=1,nsppol
2169 0 : do ikpt=1,nkpt !TODO : hexu: check if it should be mkmem, or should skip if the kpt is not in this node.
2170 0 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
2171 0 : do iband=1,nband_k
2172 : ! write(std_out,*)"amn2",iband,ibg,ikpt
2173 0 : do ispinor=1,nspinor
2174 0 : icat=1
2175 0 : do itypat=1,dtset%ntypat
2176 0 : lmn_size=pawtab(itypat)%lmn_size
2177 0 : do iatom=icat,icat+nattyp(itypat)-1
2178 0 : jj1=0
2179 0 : do ilmn=1,lmn_size
2180 0 : if(iatom.le.natprjn) then
2181 : ! do iwan=1,nwan
2182 0 : do libprjn=0,lprjn(nbprjn)
2183 : ! if (psps%indlmn(1,ilmn,itypat)==proj_l(iwan)) then
2184 : ! if (psps%indlmn(2,ilmn,itypat)==mtransfo(proj_l(iwan),proj_m(iwan))) then
2185 0 : if (psps%indlmn(1,ilmn,itypat)==libprjn) then
2186 0 : if (psps%indlmn(3,ilmn,itypat)<=npprjn(libprjn)) then
2187 0 : if(band_in(iband,isppol)) then
2188 0 : jj1=jj1+1
2189 0 : if(jj1>nwan(isppol)) then
2190 0 : write(std_out,*) "number of wannier orbitals is lower than lmn_size"
2191 0 : write(std_out,*) jj1,nwan(isppol)
2192 0 : ABI_ERROR("Aborting now")
2193 : end if
2194 : !amn2(1,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(1,ilmn)
2195 : !amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(2,ilmn)
2196 : amn2(1,iatom,isppol,ikpt,iband,ispinor,jj1)= &
2197 0 : &mywfc%cprj_elem(1, ispinor, iband, ikpt, isppol, iatom, ilmn)
2198 : amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)= &
2199 0 : &mywfc%cprj_elem(2, ispinor, iband, ikpt, isppol, iatom, ilmn)
2200 :
2201 : !amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(2,ilmn)
2202 : end if
2203 : end if
2204 : end if
2205 : end do ! libprjn
2206 : ! endif
2207 : ! endif
2208 : ! enddo ! iwan
2209 : end if ! natprjn
2210 : end do !ilmn
2211 : end do ! iatom
2212 0 : icat=icat+nattyp(itypat)
2213 : end do ! itypat
2214 : end do ! ispinor
2215 : end do !iband
2216 0 : ibg=ibg+nband_k*nspinor
2217 : ! write(std_out,*)'amn2b',iband,ibg,ikpt
2218 : end do !ikpt
2219 : end do ! isppol
2220 :
2221 : ! ----------------------- Save Amn --------------------
2222 0 : do isppol=1,nsppol
2223 0 : do ikpt=1,nkpt
2224 0 : do iband2=1,nwan(isppol)
2225 : jband=0
2226 0 : do iband1=1,mband
2227 0 : if(band_in(iband1,isppol)) then
2228 0 : jband=jband+1
2229 : A_matrix(jband,iband2,ikpt,isppol)=&
2230 0 : cmplx(amn2(1,1,1,ikpt,iband1,1,iband2),amn2(2,1,1,ikpt,iband1,1,iband2))
2231 : end if
2232 : end do
2233 : end do
2234 : end do
2235 : end do
2236 : end if !usepaw
2237 0 : ABI_FREE(amn2)
2238 0 : ABI_FREE(npprjn)
2239 0 : ABI_FREE(lprjn)
2240 :
2241 : end if ! lproj==3
2242 :
2243 19 : end subroutine mlwfovlp_proj
2244 : !!***
2245 :
2246 : !!****f* m_mlwfovlp/mlwfovlp_projpaw
2247 : !! NAME
2248 : !! mlwfovlp_projpaw
2249 : !!
2250 : !! FUNCTION
2251 : !! Calculates the functions that are given to Wannier90 as an starting guess.
2252 : !! Here we project them inside the PAW spheres
2253 : !!
2254 : !! INPUTS
2255 : !! band_in(mband)= logical array which indicates the bands to be excluded from the calculation
2256 : !! cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
2257 : !! and each |p_lmn> non-local projector
2258 : !! just_augmentation= flag used to indicate that we are just going
2259 : !! to compute augmentation part of the matrix
2260 : !! and we are excluding the plane wave part.
2261 : !! mband= maximum number of bands
2262 : !! mkmem= number of k points which can fit in memory; set to 0 if use disk
2263 : !! natom= number of atoms in cell.
2264 : !! nband(nkpt*nsppol)= array cointaining number of bands at each k-point and isppol
2265 : !! nkpt=number of k points.
2266 : !! num_bands=number of bands actually used to construct the wannier function (NOT USED IN 6.7.1 SO WAS TEMPORARILY REMOVED)
2267 : !! nspinor=number of spinorial components of the wavefunctions
2268 : !! nsppol=1 for unpolarized, 2 for spin-polarized
2269 : !! ntypat=number of types of atoms in unit cell.
2270 : !! nwan= number of wannier fonctions (read in wannier90.win).
2271 : !! pawrad(ntypat)= type(pawrad_type) radial information of paw objects
2272 : !! pawtab(ntypat)= For PAW, TABulated data initialized at start
2273 : !! proj_l(mband)= angular part of the projection function (quantum number l)
2274 : !! proj_m(mband)= angular part of the projection function (quantum number m)
2275 : !! proj_radial(mband)= radial part of the projection.
2276 : !! proj_site(3,mband)= site of the projection.
2277 : !! proj_x(3,mband)= x axis for the projection.
2278 : !! proj_z(3,mband)= z axis for the projection.
2279 : !! proj_zona(mband)= extension of the radial part.
2280 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
2281 : !! rprimd(3,3)= Direct lattice vectors, Bohr units.
2282 : !! typat(natom)= atom type
2283 : !! xred(3,natom)=reduced dimensionless atomic coordinates
2284 : !!
2285 : !! OUTPUT
2286 : !! A_paw(max_num_bands,nwan,nkpt) = A matrix containing initial guess for MLWFs
2287 : !! (augmentation part of the matrix)
2288 : !!
2289 : !! NOTES
2290 : !! This routine is still under developement
2291 : !!
2292 : !! SOURCE
2293 :
2294 0 : subroutine mlwfovlp_projpaw(A_paw,band_in,mywfc,just_augmentation,max_num_bands,mband,mkmem,&
2295 0 : &mwan,natom,nband,nkpt,&
2296 0 : &nspinor,nsppol,ntypat,nwan,pawrad,pawtab,&
2297 0 : &proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,&
2298 0 : &rprimd,typat,xred)
2299 :
2300 : !Arguments ------------------------------------
2301 : integer,intent(in) :: max_num_bands,mband,mkmem,mwan,natom,nkpt
2302 : integer,intent(in) :: nspinor,nsppol,ntypat
2303 : !arrays
2304 : integer,intent(in) :: nband(nsppol*nkpt),nwan(nsppol)
2305 : integer,intent(in) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol)
2306 : integer,intent(in) :: typat(natom)
2307 : real(dp),intent(in):: proj_site(3,mband,nsppol)
2308 : real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol)
2309 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2310 : complex(dp),intent(out) :: A_paw(max_num_bands,mwan,nkpt,nsppol)
2311 : logical,intent(in) :: band_in(mband,nsppol)
2312 : logical,intent(in)::just_augmentation(mwan,nsppol)
2313 : !type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
2314 : type(abstract_wf), intent(inout) :: mywfc
2315 : type(pawrad_type),intent(in) :: pawrad(ntypat)
2316 : type(pawtab_type),intent(in) :: pawtab(ntypat)
2317 : type(pseudopotential_type),intent(in) :: psps
2318 :
2319 : !Local variables-------------------------------
2320 : !local variables
2321 : integer :: basis_size,iatom,iband,ii
2322 : integer :: ikpt,ir,isppol,itypat,iwan,jband
2323 : integer :: ll,lm,ln,mm,ilmn
2324 : integer :: lmn_size,max_lmax2, mesh_size,nn
2325 0 : integer :: lmax(nsppol),lmax2(nsppol)
2326 : real(dp):: aa,int_rad2,prod_real,prod_imag
2327 : real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0
2328 : real(dp):: sum,wan_lm_fac,x
2329 : complex(dp)::prod
2330 : character(len=500) :: msg
2331 : !arrays
2332 0 : integer :: index(mband,nkpt,nsppol)
2333 0 : real(dp) :: dist,norm(mwan,nsppol)
2334 0 : real(dp) :: proj_cart(3,mwan,nsppol),proj_site_unit(3,mwan,nsppol)
2335 0 : real(dp) :: xcart_unit(3,natom),xred_unit(3,natom)
2336 0 : real(dp),allocatable :: aux(:),ff(:),r(:),int_rad(:),rad_int(:)
2337 0 : real(dp),allocatable :: ylmr_fac(:,:,:)
2338 : integer,parameter :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/) ! Tables 3.1 & 3.2, User guide
2339 :
2340 : ! *************************************************************************
2341 :
2342 : ABI_UNUSED(mkmem)
2343 : ABI_UNUSED(nspinor)
2344 :
2345 0 : write(msg, '(a,a)' )ch10,'** mlwfovlp_proj: compute in-sphere part of A_matrix'
2346 0 : call wrtout(std_out,msg)
2347 :
2348 : !Check input variables
2349 0 : do isppol=1,nsppol
2350 0 : do iwan=1,nwan(nsppol)
2351 0 : if(proj_radial(iwan,isppol)<1 .or. proj_radial(iwan,isppol)>4)then
2352 : write(msg,'(a,a,a,i0)')&
2353 0 : ' proj_radial should be between 1 and 4,',ch10,&
2354 0 : ' however, proj_radial=',proj_radial(iwan,isppol)
2355 0 : ABI_BUG(msg)
2356 : end if
2357 : end do
2358 : end do
2359 :
2360 : !Initialize
2361 0 : A_paw(:,:,:,:)=cmplx(0.d0,0.d0)
2362 :
2363 : !Get index for cprj
2364 0 : ii=0
2365 0 : do isppol=1,nsppol
2366 0 : do ikpt=1,nkpt
2367 0 : do iband=1,nband(ikpt)
2368 0 : ii=ii+1
2369 0 : index(iband,ikpt,isppol)=ii
2370 : end do
2371 : end do
2372 : end do
2373 :
2374 : !obtain lmax and lmax2
2375 0 : lmax(:)=0
2376 0 : lmax2(:)=0
2377 0 : do isppol=1,nsppol
2378 0 : do iwan=1,nwan(isppol)
2379 0 : lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iwan,isppol)))
2380 : end do !iwan
2381 0 : lmax2(isppol)=(lmax(isppol)+1)**2
2382 : end do
2383 0 : max_lmax2=maxval(lmax2(:))
2384 : !
2385 : !get ylmfac, factor used for rotations and hybrid orbitals
2386 : !
2387 0 : ABI_MALLOC(ylmr_fac,(max_lmax2,mwan,nsppol))
2388 :
2389 :
2390 0 : do isppol=1,nsppol
2391 : call mlwfovlp_ylmfar(ylmr_fac(1:lmax2(isppol),1:nwan(isppol),isppol),&
2392 : & lmax(isppol),lmax2(isppol),mband,nwan(isppol),proj_l(:,isppol),proj_m(:,isppol),&
2393 0 : & proj_x(:,:,isppol),proj_z(:,:,isppol))
2394 : !
2395 : ! Shift projection centers and atom centers to the primitive cell
2396 : ! This will be useful after, when we check if the Wannier function
2397 : ! lies on one specific atom
2398 : !
2399 0 : proj_site_unit(:,:,:)=0.d0
2400 0 : do iwan=1,nwan(isppol)
2401 0 : do ii=1,3
2402 0 : proj_site_unit(ii,iwan,isppol)=ABS(proj_site(ii,iwan,isppol)-AINT(proj_site(ii,iwan,isppol)) )
2403 : end do
2404 : end do
2405 0 : do iatom=1,natom
2406 0 : do ii=1,3
2407 0 : xred_unit(ii,iatom)=ABS(xred(ii,iatom)-AINT(xred(ii,iatom)) )
2408 : end do
2409 : end do
2410 0 : call xred2xcart(natom,rprimd,xcart_unit,xred_unit)
2411 0 : call xred2xcart(mwan,rprimd,proj_cart(:,:,isppol),proj_site_unit(:,:,isppol))
2412 : !
2413 : ! Normalize the Wannier functions
2414 : !
2415 : ! Radial part
2416 0 : mesh_size= nint((rmax - xmin ) / dx + 1)
2417 0 : ABI_MALLOC( ff,(mesh_size))
2418 0 : ABI_MALLOC(r,(mesh_size))
2419 0 : ABI_MALLOC(rad_int,(mesh_size))
2420 0 : ABI_MALLOC(aux,(mesh_size))
2421 0 : do ir=1, mesh_size
2422 0 : x=xmin+DBLE(ir-1)*dx
2423 0 : r(ir)=x
2424 : end do !ir
2425 0 : do iwan=1,nwan(isppol)
2426 : ! write(std_out,*)'iwan',iwan
2427 : ! radial functions shown in table 3.3 of wannier90 manual
2428 0 : if(proj_radial(iwan,isppol)==1) ff(:) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) * exp(-proj_zona(iwan,isppol)*r(:))
2429 0 : if(proj_radial(iwan,isppol)==2) ff(:) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *&
2430 0 : & (2.d0 - proj_zona(iwan,isppol)*r(:))*exp(-proj_zona(iwan,isppol)*r(:)/2.d0)
2431 0 : if(proj_radial(iwan,isppol)==3) ff(:) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)&
2432 : & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*r(:)/3.d0 + 2.d0*proj_zona(iwan,isppol)**2*r(:)**2/27.d0)&
2433 0 : & * exp(-proj_zona(iwan,isppol) * r(:)/3.d0)
2434 :
2435 0 : if(proj_radial(iwan,isppol)/=4) then
2436 0 : aux(:)=ff(:)**2*r(:)**2
2437 0 : call simpson_int(mesh_size,dx,aux,rad_int)
2438 0 : sum=0.d0
2439 0 : do ir=1,mesh_size
2440 0 : sum=sum+rad_int(ir)
2441 : end do
2442 0 : int_rad2=sum/real(mesh_size,dp)
2443 : !
2444 : ! do ir=1,mesh_size
2445 : ! if(iwan==1) write(400,*)r(ir),aux(ir),rad_int(ir)
2446 : ! end do
2447 : else
2448 : !
2449 : ! ==4: gaussian function
2450 : ! f(x)=\exp(-1/4(x/aa)**2)
2451 : ! \int f(x)f(x) dx = \int \exp(-1/2(x/aa)**2) = aa*sqrt(2pi)
2452 : !
2453 0 : int_rad2=sqrt(2.d0*pi)*proj_zona(iwan,isppol)
2454 : end if
2455 :
2456 : !
2457 : ! Now angular part
2458 : !
2459 0 : prod_real=0.d0
2460 0 : do lm=1,lmax2(isppol)
2461 0 : wan_lm_fac=ylmr_fac(lm,iwan,isppol)
2462 : ! write(std_out,*)'wan_lm_fac',wan_lm_fac
2463 : ! write(std_out,*)'int_rad2',int_rad2
2464 0 : prod_real= prod_real + wan_lm_fac**2 * int_rad2
2465 : end do
2466 0 : norm(iwan,isppol)=sqrt(prod_real)
2467 : end do !iwan
2468 0 : ABI_FREE(ff)
2469 0 : ABI_FREE(r)
2470 0 : ABI_FREE(rad_int)
2471 0 : ABI_FREE(aux)
2472 : !
2473 : ! Now that we found our guiding functions
2474 : ! We proceed with the internal product of
2475 : ! our guiding functions and the wave function
2476 : ! Amn=<G_m|\Psi_n> inside the sphere.
2477 : ! The term <G_m|\Psi_n> inside the sphere is:
2478 : ! = \sum_i <G_n | \phi_i - \tphi_i> <p_im|\Psi_m>
2479 : !
2480 : !
2481 : ! G_n \phi and \tphi can be decomposed in
2482 : ! a radial function times an angular function.
2483 : !
2484 : !
2485 : ! Big loop on iwan and iatom
2486 : !
2487 0 : do iwan=1,nwan(isppol)
2488 0 : do iatom=1,natom
2489 : !
2490 : ! check if center of wannier function coincides
2491 : ! with the center of the atom
2492 : !
2493 : dist=((proj_cart(1,iwan,isppol)-xcart_unit(1,iatom))**2 + &
2494 : (proj_cart(2,iwan,isppol)-xcart_unit(2,iatom))**2 + &
2495 0 : (proj_cart(3,iwan,isppol)-xcart_unit(3,iatom))**2)**0.5
2496 :
2497 : ! if the distance between the centers is major than 0.1 angstroms skip
2498 0 : if( dist > 0.188972613) cycle
2499 0 : write(msg, '(2a,i4,a,i4,2a)')ch10, ' Wannier function center',iwan,' is on top of atom',&
2500 0 : iatom,ch10,' Calculating in-sphere contribution'
2501 0 : call wrtout(ab_out,msg)
2502 0 : call wrtout(std_out,msg)
2503 : ! Get useful quantities
2504 0 : itypat=typat(iatom)
2505 0 : lmn_size=pawtab(itypat)%lmn_size
2506 0 : basis_size=pawtab(itypat)%basis_size
2507 0 : mesh_size=pawtab(itypat)%mesh_size
2508 0 : ABI_MALLOC(int_rad,(basis_size))
2509 0 : ABI_MALLOC(ff,(mesh_size))
2510 0 : ABI_MALLOC(aux,(mesh_size))
2511 :
2512 : ! Integrate first the radial part and save it into an array
2513 : ! radial functions shown in table 3.3 of wannier90 manual
2514 0 : if(proj_radial(iwan,isppol)==1) aux(1:mesh_size) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) *&
2515 0 : & exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size))
2516 0 : if(proj_radial(iwan,isppol)==2) aux(1:mesh_size) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *&
2517 : & (2.d0 - proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) &
2518 0 : & * exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/2.d0)
2519 0 : if(proj_radial(iwan,isppol)==3) aux(1:mesh_size) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)&
2520 : & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/3.d0 &
2521 : & + 2.d0*proj_zona(iwan,isppol)**2 *pawrad(itypat)%rad(1:mesh_size)**2/27.d0)&
2522 0 : & * exp(-proj_zona(iwan,isppol) * pawrad(itypat)%rad(1:mesh_size)/3.d0)
2523 : !
2524 : ! ==4: gaussian function
2525 : ! f(x)=\exp(-1/4(x/aa)**2)
2526 : !
2527 0 : if(proj_radial(iwan,isppol)==4) then
2528 0 : aa=1.d0/proj_zona(iwan,isppol)
2529 0 : aux(1:mesh_size)= exp(-0.25d0*(pawrad(itypat)%rad(1:mesh_size)*aa)**2)
2530 : end if
2531 : !
2532 : ! Normalize aux
2533 0 : aux(:)=aux(:)/norm(iwan,isppol)
2534 : !
2535 0 : do ln=1,basis_size
2536 0 : if(just_augmentation(iwan,isppol)) then
2537 : !
2538 : ! just augmentation region contribution
2539 : ! In this case there is no need to use \tphi
2540 : ! ff= \int R_wan(r) (R_phi(ln;r)/r ) r^2 dr
2541 : !
2542 : ff(1:mesh_size)= aux(1:mesh_size) * pawtab(itypat)%phi(1:mesh_size,ln) &
2543 0 : & * pawrad(itypat)%rad(1:mesh_size)
2544 : else
2545 : ! Inside sphere contribution = \phi - \tphi
2546 : ! ff= \int R_wan(r) (R_phi(ln;r)/r - R_tphi(ln;r)/r) r^2 dr
2547 : ff(1:mesh_size)= aux(1:mesh_size) * (pawtab(itypat)%phi(1:mesh_size,ln)-pawtab(itypat)%tphi(1:mesh_size,ln)) &
2548 0 : & * pawrad(itypat)%rad(1:mesh_size)
2549 : end if
2550 : !
2551 : ! Integration with simpson routine
2552 : !
2553 0 : call simp_gen(int_rad(ln),ff,pawrad(itypat))
2554 : ! do ii=1,mesh_size
2555 : ! unit_ln=400+ln
2556 : ! if( iwan==1 ) write(unit_ln,*)pawrad(itypat)%rad(ii),ff(ii),int_rad(ln)
2557 : ! end do
2558 : end do !ln
2559 0 : ABI_FREE(ff)
2560 0 : ABI_FREE(aux)
2561 : !
2562 : ! Now integrate the angular part
2563 : ! Cycle on i indices
2564 : !
2565 : ! prod_real=0.d0
2566 0 : do ilmn=1, lmn_size
2567 0 : ll=Psps%indlmn(1,ilmn,itypat)
2568 0 : mm=Psps%indlmn(2,ilmn,itypat)
2569 0 : nn=Psps%indlmn(3,ilmn,itypat)
2570 0 : lm=Psps%indlmn(4,ilmn,itypat)
2571 0 : ln=Psps%indlmn(5,ilmn,itypat)
2572 : ! write(std_out,*)'ll ',ll,' mm ',mm,'nn',nn,"lm",lm,"ln",ln
2573 : !
2574 : ! Get wannier factor for that lm component
2575 0 : if(lm <=lmax2(isppol)) then
2576 0 : wan_lm_fac=ylmr_fac(lm,iwan,isppol)
2577 : ! Make delta product
2578 : ! Here we integrate the angular part
2579 : ! Since the integral of the product of two spherical harmonics
2580 : ! is a delta function
2581 0 : if( abs(wan_lm_fac) > 0.0d0) then
2582 : ! write(std_out,*) 'll',ll,'mm',mm,'lm',lm,'ln',ln,'factor',wan_lm_fac !lm index for wannier function
2583 : !
2584 : ! Calculate Amn_paw, now that the radial and angular integrations are done
2585 : !
2586 0 : prod=cmplx(0.d0,0.d0)
2587 0 : do ikpt=1,nkpt
2588 0 : jband=0
2589 : ! NOTE: hexu: this doesn't seem right for nspinor=2
2590 : ! NOTE: also nband size is (nsppol*nkpt)
2591 0 : do iband=1,nband(ikpt)
2592 0 : if(band_in(iband,isppol)) then
2593 0 : jband=jband+1
2594 :
2595 : !prod_real= cprj(iatom,index(iband,ikpt,isppol))%cp(1,ilmn) * int_rad(ln) * wan_lm_fac
2596 : !prod_imag= cprj(iatom,index(iband,ikpt,isppol))%cp(2,ilmn) * int_rad(ln) * wan_lm_fac
2597 : ! FIXME: here ispinor is set to 1
2598 : ! There should be a loop over ispinor
2599 0 : prod_real= mywfc%cprj_elem(1, 1, iband, ikpt, isppol, iatom, ilmn ) * int_rad(ln) * wan_lm_fac
2600 0 : prod_imag= mywfc%cprj_elem(2, 1, iband, ikpt, isppol, iatom, ilmn ) * int_rad(ln) * wan_lm_fac
2601 0 : prod=cmplx(prod_real,prod_imag)
2602 :
2603 0 : A_paw(jband,iwan,ikpt,isppol)=A_paw(jband,iwan,ikpt,isppol)+prod
2604 : end if !band_in
2605 : end do !iband
2606 : end do !ikpt
2607 : !
2608 : end if !lm<=lmax2
2609 : end if ! abs(wan_lm_fac) > 0.0d0
2610 : end do !ilmn=1, lmn_size
2611 0 : ABI_FREE(int_rad)
2612 : end do !iatom
2613 : end do !iwan
2614 : end do !isppol
2615 :
2616 : !Deallocate quantities
2617 0 : ABI_FREE(ylmr_fac)
2618 :
2619 0 : end subroutine mlwfovlp_projpaw
2620 : !!***
2621 :
2622 : !!****f* m_mlwfovlp/mlwfovlp_radial
2623 : !! NAME
2624 : !! mlwfovlp_radial
2625 : !!
2626 : !! FUNCTION
2627 : !! Calculates the radial part of the initial functions given to Wannier90
2628 : !! as an starting point for the minimization.
2629 : !! The trial functions are a set of solutions to the radial part of the hydrogenic
2630 : !! Schrodinger equation as it is explained in Table 3.3 of the Wannier90 user guide.
2631 : !!
2632 : !! INPUTS
2633 : !! alpha= Z/a = zona
2634 : !! lmax= maximum value of l
2635 : !! rvalue= integer defining the choice for radial functions R(r).
2636 : !! It can take values from 1-3.
2637 : !! It is associted to the radial part of the hydrogenic Schrodinger equation for l=0,
2638 : !! See the manual of Wannier90 for more information. (www.wannier.org)
2639 : !! xx= scalar number used to calculate the spherical bessel function. J_il(xx)
2640 : !!
2641 : !! OUTPUT
2642 : !! mlwfovlp_radial= radial part for initial projections used to construct MLWF
2643 : !!
2644 : !! SIDE EFFECTS
2645 : !! None
2646 : !!
2647 : !! NOTES
2648 : !! Calculates the radial part of the initial functions given as an initial
2649 : !! guess by the user to construct the MLWF.
2650 : !!
2651 : !! SOURCE
2652 :
2653 876185 : subroutine mlwfovlp_radial(alpha,lmax,lmax2,radial,rvalue,xx)
2654 :
2655 : !Arguments ------------------------------------
2656 : !scalars
2657 : integer,intent(in) :: lmax,lmax2,rvalue
2658 : real(dp),intent(in) :: alpha,xx
2659 : !arrays
2660 : real(dp),intent(out) :: radial(lmax2)
2661 :
2662 : !Local variables
2663 : !scalars
2664 : integer :: ir,ll,lm,mesh,mm
2665 : real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0
2666 : real(dp) :: aa,ftmp,gauss,rtmp,x
2667 : character(len=500) :: msg
2668 : !arrays
2669 : real(dp),parameter :: dblefact(4)=(/1_dp,3_dp,15_dp,105_dp/)
2670 876185 : real(dp),allocatable :: aux(:),bes(:),cosr(:),func_r(:),r(:),rad_int(:),sinr(:)
2671 :
2672 : ! *************************************************************************
2673 :
2674 : !Radial functions in the form of hydrogenic orbitals as defined in the
2675 : !wannier90 manual.
2676 876185 : if(( rvalue > 0 ).and.(rvalue < 4)) then
2677 :
2678 : ! mesh
2679 710529 : mesh= nint((rmax - xmin ) / dx + 1)
2680 710529 : ABI_MALLOC( bes,(mesh))
2681 710529 : ABI_MALLOC(func_r,(mesh))
2682 710529 : ABI_MALLOC(r,(mesh))
2683 710529 : ABI_MALLOC(rad_int,(mesh))
2684 710529 : ABI_MALLOC( aux,(mesh))
2685 710529 : ABI_MALLOC(cosr,(mesh))
2686 710529 : ABI_MALLOC(sinr,(mesh))
2687 475343901 : do ir=1, mesh
2688 474633372 : x=xmin+DBLE(ir-1)*dx
2689 475343901 : r(ir)=x
2690 : end do !ir
2691 :
2692 : ! radial functions shown in table 3.3 of wannier90 manual
2693 475343901 : if (rvalue==1) func_r(:) = 2.d0 * alpha**(3.d0/2.d0) * exp(-alpha*r(:))
2694 710529 : if (rvalue==2) func_r(:) = 1.d0/(2.d0*sqrt(2.d0))*alpha**(3.d0/2.d0) *&
2695 0 : & (2.d0 - alpha*r(:))*exp(-alpha*r(:)/2.d0)
2696 710529 : if (rvalue==3) func_r(:) = sqrt(4.d0/27.d0)*alpha**(3.d0/2.d0)&
2697 : & * (1.d0 - 2.d0*alpha*r(:)/3.d0 + 2.d0*alpha**2*r(:)**2/27.d0)&
2698 0 : & * exp(-alpha * r(:)/3.d0)
2699 :
2700 : ! compute spherical bessel functions
2701 475343901 : cosr(:)=cos(xx*r(:))
2702 475343901 : sinr(:)=sin(xx*r(:))
2703 710529 : lm=0
2704 2374640 : do ll=0,lmax
2705 1664111 : call besjm(xx,bes,cosr,ll,mesh,sinr,r)
2706 1113290259 : aux(:)=bes(:)*func_r(:)*r(:)
2707 : ! do ir=1,mesh
2708 : ! write(310,*) r(ir),bes(ir)
2709 : ! end do
2710 1664111 : call simpson_int(mesh,dx,aux,rad_int)
2711 1664111 : rtmp=rad_int(mesh)/mesh
2712 6582125 : do mm=-ll,ll
2713 4207485 : lm=lm+1
2714 5871596 : radial(lm)=rtmp
2715 : end do !mm
2716 : end do !ll
2717 710529 : ABI_FREE(bes)
2718 710529 : ABI_FREE(func_r)
2719 710529 : ABI_FREE(r)
2720 710529 : ABI_FREE(aux)
2721 710529 : ABI_FREE(rad_int)
2722 710529 : ABI_FREE(cosr)
2723 710529 : ABI_FREE(sinr)
2724 :
2725 : ! Radial part in the form of Gaussian functions of a given width
2726 : ! Taken by code of made by drh.
2727 165656 : elseif ( rvalue == 4) then
2728 165656 : aa=1._dp/alpha
2729 165656 : gauss=exp(-0.25_dp*(aa*xx)**2)
2730 165656 : lm=0
2731 448764 : do ll=0,lmax
2732 283108 : ftmp=(0.5_dp*pi)**(0.25_dp)*aa*sqrt(aa/dblefact(ll+1))*(aa*xx)**ll*gauss
2733 966776 : do mm=-ll,ll
2734 518012 : lm=lm+1
2735 801120 : radial(lm)=ftmp
2736 : end do
2737 : end do
2738 : else ! rvalue < 0 of rvalue > 4
2739 : write(msg,'(a,i6,5a)')&
2740 0 : ' Radial function r=',rvalue,ch10,&
2741 0 : ' is not defined',ch10,&
2742 0 : ' Modify .win file',ch10
2743 0 : ABI_BUG(msg)
2744 : end if !rvalue
2745 :
2746 876185 : end subroutine mlwfovlp_radial
2747 : !!***
2748 :
2749 : !!****f* m_mlwfovlp/mlwfovlp_ylmfac
2750 : !! NAME
2751 : !! mlwfovlp_ylmfac
2752 : !!
2753 : !! FUNCTION
2754 : !! Routine that produces a factor by which the initial
2755 : !! guess of functions will be multiplied for the Wannier90 interface.
2756 : !! It is just used if there are rotations, or if the functions required
2757 : !! are linear combinations of the ylm real functions.
2758 : !!
2759 : !! Example,
2760 : !! For a function G(r)= 1/2 s + 1/3 px - 1/2 pz
2761 : !! it would produce a matrix of the following form:
2762 : !! [1/2,-1/2,1/3,0,0...0]
2763 : !!
2764 : !! The real spherical harmonics are given as factors of complex spherical harmonics
2765 : !! The real spherical harmonics are given in table 3.1 of Wannier90 user guide.
2766 : !!
2767 : !! INPUTS
2768 : !! lmax= maximum l value for spherical harmonics
2769 : !! lmax2=number of ylm functions
2770 : !! mband=maximum number of bands
2771 : !! nwan = number of wannier functions
2772 : !! proj_l(mband)= angular part of the projection function (quantum number l)
2773 : !! proj_m(mband)= angular part of the projection function (quantum number m)
2774 : !! proj_x(3,mband)= x axis for the projection.
2775 : !! proj_z(3,mband)= z axis for the projection.
2776 : !!
2777 : !! OUTPUT
2778 : !! ylmc_fac(lmax2,nwan)=matrix containig a factor for ylm hybrid orbitals
2779 : !!
2780 : !! SIDE EFFECTS
2781 : !! (only writing, printing)
2782 : !!
2783 : !! SOURCE
2784 :
2785 :
2786 14 : subroutine mlwfovlp_ylmfac(ylmc_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z)
2787 :
2788 : !Arguments ------------------------------------
2789 : integer, intent(in):: lmax,lmax2,nwan,mband
2790 : ! arrays
2791 : integer,intent(in) :: proj_l(mband),proj_m(mband)
2792 : real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband)
2793 : complex(dp),intent(out)::ylmc_fac(lmax2,nwan)
2794 : !
2795 : !Local variables-------------------------------
2796 : !
2797 : integer :: orb_idx(16)=(/1,3,4,2,7,8,6,9,5,13,14,12,15,11,16,10/) !Tab3.1 Wannier90 user guide
2798 : integer :: idum,ii,info,inversion_flag
2799 : integer :: ir,iwan,jj,ll,lm,lmc,mm,mr
2800 : real(dp):: onem,test
2801 : ! arrays
2802 28 : integer:: ipiv(lmax2)
2803 28 : real(dp)::r(3,lmax2),rp(3,lmax2)
2804 : real(dp)::rs2,rs3,rs6,rs12,umat(3,3)
2805 28 : complex(dp)::crot(lmax2,lmax2),ctor(lmax2,lmax2),orb_lm(lmax2,-5:3,7)
2806 28 : complex(dp):: ylmcp(lmax2)
2807 28 : complex(dp):: ylmc_rr(lmax2,lmax2),ylmc_rr_save(lmax2,lmax2)
2808 28 : complex(dp):: ylmc_rrinv(lmax2,lmax2),ylmc_rp(lmax2,lmax2)
2809 : complex(dp),parameter :: c0=(0._dp,0._dp),c1=(1._dp,0._dp),ci=(0._dp,1._dp)
2810 : character(len=500) :: msg
2811 :
2812 : ! *************************************************************************
2813 :
2814 :
2815 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2816 : !DEBUG
2817 : !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2
2818 : !write(std_out,*)'mband ',mband,'nwan ',nwan
2819 : !
2820 : !do iwan=1,nwan
2821 : !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan)
2822 : !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan)
2823 : !end do
2824 : !!END DEBUG
2825 :
2826 : !constants for linear combinations of ylm's
2827 14 : rs2=1._dp/sqrt(2._dp)
2828 14 : rs3=1._dp/sqrt(3._dp)
2829 14 : rs6=1._dp/sqrt(6._dp)
2830 14 : rs12=1._dp/sqrt(12._dp)
2831 :
2832 : !complex lm coefficients for real spherical harmonics in conventional order
2833 : !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2,
2834 : !fz(x2-y2),fx(x2-3y2)
2835 468 : ctor(:,:)=c0
2836 43 : do ll=0,lmax
2837 29 : mm=0
2838 29 : lm= ll**2+ll+mm+1
2839 29 : ctor(lm,lm)=c1
2840 43 : if(ll>0) then
2841 : onem=one
2842 33 : do mm=1,ll
2843 18 : onem=-onem !(-1^mm)
2844 18 : lm= ll**2+ll+mm+1
2845 18 : lmc=ll**2+ll-mm+1
2846 18 : ctor(lm ,lm )=rs2*c1
2847 18 : ctor(lmc,lm )=onem*rs2*c1
2848 18 : ctor(lm ,lmc)=rs2*ci
2849 33 : ctor(lmc,lmc)=-onem*rs2*ci
2850 : end do
2851 : end if
2852 : end do
2853 :
2854 : lm=0
2855 43 : do ll=0,lmax
2856 108 : do mm=-ll,ll
2857 65 : lm=lm+1
2858 483 : ctor(:,lm)=ctor(:,lm)*conjg(ci)**ll
2859 : end do !mm
2860 : end do !ll
2861 :
2862 :
2863 : !coefficients for basic wannier orbitals in Table 3.1 order
2864 5089 : orb_lm(:,:,:)=c0
2865 : ii=0
2866 43 : do ll=0,lmax
2867 108 : do mr=1,2*ll+1
2868 65 : ii=ii+1
2869 483 : orb_lm(:,ll,mr)=ctor(:,orb_idx(ii))
2870 : end do
2871 : end do
2872 :
2873 :
2874 :
2875 : !coefficients for linear combinations in table 3.2 order
2876 14 : if(lmax>=1) then
2877 : ! s px
2878 75 : orb_lm(:,-1,1)=rs2*ctor(:,1)+rs2*ctor(:,4)
2879 75 : orb_lm(:,-1,2)=rs2*ctor(:,1)-rs2*ctor(:,4)
2880 : ! s px py
2881 75 : orb_lm(:,-2,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2)
2882 75 : orb_lm(:,-2,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2)
2883 75 : orb_lm(:,-2,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4)
2884 : ! s px py pz
2885 75 : orb_lm(:,-3,1)=half*(ctor(:,1)+ctor(:,4)+ctor(:,2)+ctor(:,3))
2886 75 : orb_lm(:,-3,2)=half*(ctor(:,1)+ctor(:,4)-ctor(:,2)-ctor(:,3))
2887 75 : orb_lm(:,-3,3)=half*(ctor(:,1)-ctor(:,4)+ctor(:,2)-ctor(:,3))
2888 75 : orb_lm(:,-3,4)=half*(ctor(:,1)-ctor(:,4)-ctor(:,2)+ctor(:,3))
2889 : end if
2890 14 : if(lmax>=2) then
2891 : ! s px py
2892 30 : orb_lm(:,-4,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2)
2893 30 : orb_lm(:,-4,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2)
2894 30 : orb_lm(:,-4,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4)
2895 : ! pz dz2
2896 30 : orb_lm(:,-4,4)= rs2*ctor(:,3)+rs2*ctor(:,7)
2897 30 : orb_lm(:,-4,5)=-rs2*ctor(:,3)+rs2*ctor(:,7)
2898 : ! s px dz2 dx2-y2
2899 30 : orb_lm(:,-5,1)=rs6*ctor(:,1)-rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9)
2900 30 : orb_lm(:,-5,2)=rs6*ctor(:,1)+rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9)
2901 : ! s py dz2 dx2-y2
2902 30 : orb_lm(:,-5,3)=rs6*ctor(:,1)-rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9)
2903 30 : orb_lm(:,-5,4)=rs6*ctor(:,1)+rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9)
2904 : ! s pz dz2
2905 30 : orb_lm(:,-5,5)=rs6*ctor(:,1)-rs2*ctor(:,3)+rs3*ctor(:,7)
2906 30 : orb_lm(:,-5,6)=rs6*ctor(:,1)+rs2*ctor(:,3)+rs3*ctor(:,7)
2907 : end if
2908 :
2909 : !stuff complex wannier orbital coefficient array
2910 89 : do iwan=1,nwan
2911 500 : ylmc_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan))
2912 : end do
2913 :
2914 :
2915 : !setup to rotate ylmc_fac to new axes if called for
2916 : !skip if only s projectors are used
2917 14 : if ( lmax>0 ) then
2918 : ! generate a set of nr=lmax2 random vectors
2919 : ! idum=123456
2920 75 : do ir=1,lmax2
2921 252 : do ii=1,3
2922 252 : r(ii,ir) = uniformrandom(idum)-0.5d0
2923 : end do !ii
2924 63 : call ylm_cmplx(lmax,ylmcp,r(1,ir),r(2,ir),r(3,ir))
2925 450 : ylmc_rr(ir,:)=conjg(ylmcp(:))
2926 462 : ylmc_rr_save(ir,:)=conjg(ylmcp(:))
2927 : end do !ir
2928 :
2929 462 : ylmc_rrinv(:,:)=c0
2930 75 : do ii=1,lmax2
2931 75 : ylmc_rrinv(ii,ii)=c1
2932 : end do !ii
2933 : ! calculate inverse of ylmc(ir,lm) matrix
2934 12 : call ZGESV(lmax2,lmax2,ylmc_rr,lmax2,ipiv,ylmc_rrinv,lmax2,info)
2935 :
2936 : ! check that r points are independent (ie., that matrix inversion wasn't
2937 : ! too close to singular)
2938 3675 : ylmc_rr=matmul(ylmc_rrinv,ylmc_rr_save)
2939 12 : test=zero
2940 75 : do ii=1,lmax2
2941 63 : ylmc_rr(ii,ii)=ylmc_rr(ii,ii)-c1
2942 462 : do jj=1,lmax2
2943 450 : test=max(abs(ylmc_rr(ii,jj)),test)
2944 : end do !ii
2945 : end do !jj
2946 12 : if(test>tol8) then
2947 : write(msg, '(5a)' )&
2948 0 : & ' matrix inversion error for wannier rotations',ch10,&
2949 0 : & ' random vectors r(j,1:nr) are not all independent !! ',ch10,&
2950 0 : & ' Action : re-seed uniformrandom or maybe just try again'
2951 0 : ABI_ERROR(msg)
2952 : end if !test>tol8
2953 :
2954 : ! end of the preliminaries, now to the rotations of the wannier orbitals
2955 79 : do iwan=1,nwan
2956 : ! don't bother for s orbitals
2957 67 : if(proj_l(iwan)==0) cycle
2958 : ! check for default axes and cycle if found
2959 : if(proj_z(1,iwan)==zero .and. proj_z(2,iwan)==zero .and.&
2960 : & proj_z(3,iwan)== one .and. proj_x(1,iwan)==one .and.&
2961 63 : & proj_x(2,iwan)==zero .and. proj_x(3,iwan)==zero) cycle
2962 :
2963 : ! get the u matrix that rotates the reference frame
2964 4 : call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat)
2965 :
2966 : ! find rotated r-vectors. Optional inversion
2967 : ! operation is an extension of the wannier90 axis-setting options
2968 : ! which only allow for proper axis rotations
2969 4 : if(inversion_flag==1) then
2970 0 : rp(:,:)= -matmul ( umat(:,:), r(:,:) )
2971 : else
2972 276 : rp(:,:) = matmul ( umat(:,:) , r(:,:) )
2973 : end if !inversion_flag
2974 :
2975 20 : do ir=1,lmax2
2976 : ! get the ylm representation of the rotated vectors
2977 16 : call ylm_cmplx(lmax,ylmcp,rp(1,ir),rp(2,ir),rp(3,ir))
2978 84 : ylmc_rp(ir,:)=conjg(ylmcp(:))
2979 : end do !ir
2980 : ! the matrix product sum(ir) ylmc_rrinv(lm,ir)*ylmc_rp(ir,lm') gives the
2981 : ! the complex lmXlm matrix representation of the coordinate rotation
2982 420 : crot(:,:)=matmul(ylmc_rrinv(:,:),ylmc_rp(:,:))
2983 :
2984 : ! now rotate the current wannier orbital
2985 100 : ylmcp(:)=matmul(crot(:,:),ylmc_fac(:,iwan))
2986 91 : ylmc_fac(:,iwan)=ylmcp(:)
2987 :
2988 : ! write(std_out,*)'ylmc_fac',ylmc_fac(:,iwan)
2989 : end do !iwan
2990 : end if !lmax>0
2991 :
2992 14 : end subroutine mlwfovlp_ylmfac
2993 : !!***
2994 :
2995 : !!****f* m_mlwfovlp/mlwfovlp_ylmfar
2996 : !! NAME
2997 : !! mlwfovlp_ylmfar
2998 : !!
2999 : !! FUNCTION
3000 : !! Routine that produces a fator by which the initial
3001 : !! guess of functions will be multiplied for the Wannier90 interface.
3002 : !! It is just used if there are rotations, or if the functions required
3003 : !! are linear combinations of the ylm real functions.
3004 : !!
3005 : !! Example,
3006 : !! For a function G(r)= 1/2 s + 1/3 px - 1/2 pz
3007 : !! it would produce a matrix of the following form:
3008 : !! [1/2,-1/2,1/3,0,0...0]
3009 : !!
3010 : !! This function is similar to mlwfovlp_ylmfac, but the factors it uses
3011 : !! real spherical harmonics instead of complex
3012 : !! spherical harmonics. Remember that real spherical harmonics
3013 : !! are linear combinations of complex
3014 : !! spherical harmonics
3015 : !!
3016 : !! INPUTS
3017 : !! lmax= maximum l value for spherical harmonics
3018 : !! lmax2=number of ylm functions
3019 : !! mband=maximum number of bands
3020 : !! nwan = number of wannier functions
3021 : !! proj_l(mband)= angular part of the projection function (quantum number l)
3022 : !! proj_m(mband)= angular part of the projection function (quantum number m)
3023 : !! proj_x(3,mband)= x axis for the projection.
3024 : !! proj_z(3,mband)= z axis for the projection.
3025 : !!
3026 : !! OUTPUT
3027 : !! ylmc_fac(lmax2,nwan)=matrix containig a factor for ylm hybrid orbitals
3028 : !!
3029 : !! SOURCE
3030 :
3031 0 : subroutine mlwfovlp_ylmfar(ylmr_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z)
3032 :
3033 : !Arguments ------------------------------------
3034 : integer, intent(in):: lmax,lmax2,nwan,mband
3035 : ! arrays
3036 : integer,intent(in) :: proj_l(mband),proj_m(mband)
3037 : real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband)
3038 : real(dp),intent(out)::ylmr_fac(lmax2,nwan)
3039 : !
3040 : !Local variables-------------------------------
3041 : !
3042 : integer :: idum,ii,inversion_flag
3043 : integer :: ir,iwan,jj,ll,lm,mm,mr
3044 : real(dp) :: onem,test
3045 : ! arrays
3046 0 : real(dp),allocatable:: dummy(:,:),nrm(:)
3047 0 : real(dp) :: r(3,lmax2),rp(3,lmax2)
3048 : real(dp) :: rs2,rs3,rs6,rs12,umat(3,3)
3049 0 : real(dp) :: rot(lmax2,lmax2),tor(lmax2,lmax2),orb_lm(lmax2,-5:3,7)
3050 0 : real(dp) :: ylmrp(lmax2)
3051 0 : real(dp) :: ylmr_rr(lmax2,lmax2),ylmr_rr_save(lmax2,lmax2)
3052 0 : real(dp) :: ylmr_rrinv(lmax2,lmax2),ylmr_rp(lmax2,lmax2)
3053 : character(len=500) :: msg ! to be uncommented, if needed
3054 : !integer :: orb_idx(16)=(/1,3,4,2,7,8,6,9,5,13,14,12,15,11,16,10/) !Tab3.1 Wannier90 user guide
3055 :
3056 : ! *************************************************************************
3057 :
3058 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3059 : !DEBUG
3060 : !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2
3061 : !write(std_out,*)'mband ',mband,'nwan ',nwan
3062 : !
3063 : !do iwan=1,nwan
3064 : !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan)
3065 : !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan)
3066 : !end do
3067 : !!END DEBUG
3068 :
3069 : !constants for linear combinations of ylm's
3070 0 : rs2=1._dp/sqrt(2._dp)
3071 0 : rs3=1._dp/sqrt(3._dp)
3072 0 : rs6=1._dp/sqrt(6._dp)
3073 0 : rs12=1._dp/sqrt(12._dp)
3074 :
3075 : !
3076 : !mapping lm coefficients for real spherical harmonics
3077 : !table 3.1 of Wannier90 user guide with real spherical harmonics in routine initylmr
3078 : !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2,
3079 : !fz(x2-y2),fx(x2-3y2)
3080 : !note: check ordering of f orbitals, it might be wrong
3081 :
3082 0 : tor(:,:)=0.d0
3083 0 : lm=0
3084 0 : do ll=0,lmax
3085 0 : do mm=-ll,ll
3086 0 : onem=(-1.d0)**mm
3087 0 : lm=lm+1
3088 0 : if(ll == 0) then
3089 0 : tor(lm,lm)=1.d0
3090 : else
3091 0 : tor(lm,lm)=onem*1.d0
3092 : end if
3093 : end do !mm
3094 : end do !ll
3095 : !do lm=1,16
3096 : !write(std_out,*)'tor lm=',lm,tor(:,lm)
3097 : !end do
3098 :
3099 : !coefficients for basic wannier orbitals in Table 3.1 order
3100 0 : orb_lm(:,:,:)=0.d0
3101 : ii=0
3102 0 : do ll=0,lmax
3103 0 : do mr=1,2*ll+1
3104 0 : ii=ii+1
3105 0 : orb_lm(:,ll,mr)= tor(:,ii)
3106 : ! write(std_out,*)'ii',ii,'orb_lm',orb_lm(:,ll,mr)
3107 : end do
3108 : end do
3109 :
3110 :
3111 :
3112 : !coefficients for linear combinations in table 3.2 order
3113 0 : if(lmax>=1) then
3114 : ! s px
3115 0 : orb_lm(:,-1,1)=rs2*tor(:,1)+rs2*tor(:,4)
3116 0 : orb_lm(:,-1,2)=rs2*tor(:,1)-rs2*tor(:,4)
3117 : ! s px py
3118 0 : orb_lm(:,-2,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2)
3119 0 : orb_lm(:,-2,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2)
3120 0 : orb_lm(:,-2,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4)
3121 : ! s px py pz
3122 0 : orb_lm(:,-3,1)=half*(tor(:,1)+tor(:,4)+tor(:,2)+tor(:,3))
3123 0 : orb_lm(:,-3,2)=half*(tor(:,1)+tor(:,4)-tor(:,2)-tor(:,3))
3124 0 : orb_lm(:,-3,3)=half*(tor(:,1)-tor(:,4)+tor(:,2)-tor(:,3))
3125 0 : orb_lm(:,-3,4)=half*(tor(:,1)-tor(:,4)-tor(:,2)+tor(:,3))
3126 : end if
3127 0 : if(lmax>=2) then
3128 : ! s px py
3129 0 : orb_lm(:,-4,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2)
3130 0 : orb_lm(:,-4,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2)
3131 0 : orb_lm(:,-4,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4)
3132 : ! pz dz2
3133 0 : orb_lm(:,-4,4)= rs2*tor(:,3)+rs2*tor(:,7)
3134 0 : orb_lm(:,-4,5)=-rs2*tor(:,3)+rs2*tor(:,7)
3135 : ! s px dz2 dx2-y2
3136 0 : orb_lm(:,-5,1)=rs6*tor(:,1)-rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9)
3137 0 : orb_lm(:,-5,2)=rs6*tor(:,1)+rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9)
3138 : ! s py dz2 dx2-y2
3139 0 : orb_lm(:,-5,3)=rs6*tor(:,1)-rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9)
3140 0 : orb_lm(:,-5,4)=rs6*tor(:,1)+rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9)
3141 : ! s pz dz2
3142 0 : orb_lm(:,-5,5)=rs6*tor(:,1)-rs2*tor(:,3)+rs3*tor(:,7)
3143 0 : orb_lm(:,-5,6)=rs6*tor(:,1)+rs2*tor(:,3)+rs3*tor(:,7)
3144 : end if
3145 :
3146 : !real wannier orbital coefficient array
3147 0 : do iwan=1,nwan
3148 0 : ylmr_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan))
3149 : end do
3150 :
3151 :
3152 : !setup to rotate ylmr_fac to new axes if called for
3153 : !skip if only s projetors are used
3154 0 : if ( lmax>0 ) then
3155 : ! generate a set of nr=lmax2 random vectors
3156 0 : idum=123456
3157 0 : do ir=1,lmax2
3158 0 : do ii=1,3
3159 0 : r(ii,ir) = uniformrandom(idum)-0.5d0
3160 : end do !ii
3161 : end do !ir
3162 0 : ABI_MALLOC(nrm,(lmax2))
3163 0 : nrm(:)=sqrt(r(1,:)**2+r(2,:)**2+r(3,:)**2)**0.5
3164 0 : call initylmr(lmax+1,1,lmax2,nrm,1,r(:,:),ylmr_rr_save(:,:),dummy)
3165 0 : ylmr_rr(:,:)=ylmr_rr_save(:,:)
3166 0 : do ir=1,lmax2
3167 0 : ylmr_rr_save(ir,:)=ylmr_rr(:,ir)
3168 : end do
3169 0 : ABI_FREE(nrm)
3170 :
3171 0 : ylmr_rrinv(:,:)=0.d0
3172 0 : do ii=1,lmax2
3173 0 : ylmr_rrinv(ii,ii)=1.d0
3174 : end do !ii
3175 : ! calculate inverse of ylmr(ir,lm) matrix
3176 0 : ylmr_rrinv(:,:)=ylmr_rr_save(:,:)
3177 0 : call matrginv(ylmr_rrinv,lmax2,lmax2)
3178 :
3179 : ! check that r points are independent (ie., that matrix inversion wasn't too close to singular)
3180 0 : ylmr_rr=matmul(ylmr_rrinv,ylmr_rr_save)
3181 0 : test=0.d0
3182 0 : do ii=1,lmax2
3183 0 : ylmr_rr(ii,ii)=ylmr_rr(ii,ii)-1.d0
3184 0 : do jj=1,lmax2
3185 0 : test=max(abs(ylmr_rr(ii,jj)),test)
3186 : end do !ii
3187 : end do !jj
3188 0 : if(test>tol8) then
3189 : write(msg, '(5a)' )&
3190 0 : ' matrix inversion error for wannier rotations',ch10,&
3191 0 : ' random vectors r(j,1:nr) are not all independent !! ',ch10,&
3192 0 : ' Action : re-seed uniformrandom or maybe just try again'
3193 0 : ABI_ERROR(msg)
3194 : end if !test>tol8
3195 :
3196 : ! end of the preliminaries, now to the rotations of the wannier orbitals
3197 0 : do iwan=1,nwan
3198 : ! don't bother for s orbitals
3199 0 : if(proj_l(iwan)==0) cycle
3200 : ! check for default axes and cycle if found
3201 : if(proj_z(1,iwan)==0.d0 .and. proj_z(2,iwan)==0.d0 .and.&
3202 : & proj_z(3,iwan)== 1.d0 .and. proj_x(1,iwan)==1.d0 .and.&
3203 0 : & proj_x(2,iwan)==0.d0 .and. proj_x(3,iwan)==0.d0) cycle
3204 :
3205 : ! get the u matrix that rotates the reference frame
3206 0 : call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat)
3207 : !
3208 : ! find rotated r-vectors. Optional inversion
3209 : ! operation is an extension of the wannier90 axis-setting options
3210 : ! which only allow for proper axis rotations
3211 0 : if(inversion_flag==1) then
3212 0 : rp(:,:)= -matmul ( umat(:,:), r(:,:) )
3213 : else
3214 0 : rp(:,:) = matmul ( umat(:,:) , r(:,:) )
3215 : end if !inversion_flag
3216 :
3217 : ! get the ylm representation of the rotated vectors
3218 0 : ABI_MALLOC(nrm,(lmax2))
3219 0 : nrm(:)=sqrt(rp(1,:)**2+rp(2,:)**2+rp(3,:)**2)**0.5
3220 0 : call initylmr(lmax+1,1,lmax2,nrm,1,rp(:,:),ylmr_rp(:,:),dummy)
3221 0 : ylmr_rr(:,:)=ylmr_rp(:,:)
3222 0 : do ir=1,lmax2
3223 0 : ylmr_rp(ir,:)=ylmr_rr(:,ir)
3224 : end do
3225 0 : ABI_FREE(nrm)
3226 : ! the matrix product sum(ir) ylmr_rrinv(lm,ir)*ylmr_rp(ir,lm') gives the
3227 : ! the lmXlm matrix representation of the coordinate rotation
3228 :
3229 0 : rot(:,:)=matmul(ylmr_rrinv(:,:),ylmr_rp(:,:))
3230 : !
3231 : ! now rotate the current wannier orbital
3232 0 : ylmrp(:)=matmul(rot(:,:),ylmr_fac(:,iwan))
3233 0 : ylmr_fac(:,iwan)=ylmrp(:)
3234 : end do !iwan
3235 : end if !lmax>0
3236 :
3237 0 : end subroutine mlwfovlp_ylmfar
3238 : !!***
3239 :
3240 : !!****f* m_mlwfovlp/wan_from_abiwan
3241 : !! NAME
3242 : !! wan_from_abiwan
3243 : !!
3244 : !! FUNCTION
3245 : !! Initialize a wan_t instance from the ABIWAN.nc netcf file
3246 : !!
3247 : !! INPUTS
3248 : !!
3249 : !! OUTPUT
3250 : !!
3251 : !! SOURCE
3252 :
3253 0 : subroutine wan_from_abiwan(wan, abiwan_filepath, spin, nsppol, keep_umats, out_prefix, comm)
3254 :
3255 : !Arguments ------------------------------------
3256 : class(wan_t),intent(out) :: wan
3257 : character(len=*),intent(in) :: abiwan_filepath, out_prefix
3258 : logical,intent(in) :: keep_umats
3259 : integer,intent(in) :: spin, nsppol, comm
3260 :
3261 : !Local variables-------------------------------
3262 : !scalars
3263 : integer,parameter :: master = 0
3264 : integer :: ncid, my_rank, ii, ir, ik, ib, jb, mb, nsppol_, mband, nwan, nkbz, ount, num_bands, nr_h, nextbands
3265 : character(len=500) :: msg
3266 : character(len=fnlen) :: out_path
3267 0 : type(crystal_t) :: cryst
3268 : complex(dp) :: ctmp
3269 : !arrays
3270 : integer :: kptrlatt(3,3)
3271 0 : integer,allocatable :: int_1d(:), int_2d(:,:) !, shiftk(:,:)
3272 0 : real(dp),allocatable :: u_mat(:,:,:,:), u_mat_opt(:,:,:,:), et_opt(:,:)
3273 0 : complex(dp),allocatable :: chs(:,:,:), chw(:,:,:)
3274 : !************************************************************************
3275 :
3276 0 : my_rank = xmpi_comm_rank(comm)
3277 :
3278 : !if (my_rank == master) then
3279 0 : wan%spin = spin
3280 0 : NCF_CHECK(nctk_open_read(ncid, abiwan_filepath, comm))
3281 :
3282 : ! Get dimensions.
3283 0 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spins", nsppol_))
3284 0 : ABI_CHECK_IEQ(nsppol, nsppol_, "Inconsistent number of spins")
3285 :
3286 0 : call cryst%ncread(ncid)
3287 :
3288 : ! NB: mband is the value of nband and not the number of bands for Wannier that is called num_bands!
3289 0 : NCF_CHECK(nctk_get_dim(ncid, "max_number_of_states", mband))
3290 0 : NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", nkbz))
3291 0 : NCF_CHECK(nctk_get_dim(ncid, "nrpts", nr_h))
3292 0 : NCF_CHECK(nctk_get_dim(ncid, "mwan", wan%max_nwan))
3293 0 : NCF_CHECK(nf90_get_var(ncid, vid("nwan"), nwan, start=[spin]))
3294 0 : wan%nkbz = nkbz; wan%nwan = nwan; wan%nr_h = nr_h
3295 :
3296 : ! Read variables for this spin.
3297 0 : NCF_CHECK(nf90_get_var(ncid, vid("num_bands"), num_bands, start=[spin]))
3298 0 : wan%num_bands = num_bands
3299 0 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "spread"), wan%spread, start=[1,spin]))
3300 :
3301 0 : ABI_MALLOC(wan%kbz, (3, nkbz))
3302 0 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_coordinates_of_kpoints"), wan%kbz))
3303 :
3304 0 : ABI_MALLOC(wan%r_h, (3, nr_h))
3305 0 : ABI_MALLOC(wan%ndegen_h, (nr_h))
3306 0 : NCF_CHECK(nf90_get_var(ncid, vid("irvec"), wan%r_h))
3307 0 : NCF_CHECK(nf90_get_var(ncid, vid("ndegen"), wan%ndegen_h))
3308 0 : NCF_CHECK(nf90_get_var(ncid, vid("have_disentangled_spin"), ii, start=[spin]))
3309 :
3310 0 : ABI_MALLOC(wan%exclude_bands, (mband))
3311 0 : NCF_CHECK(nf90_get_var(ncid, vid("exclude_bands"), wan%exclude_bands, start=[1,spin]))
3312 :
3313 0 : wan%have_disentangled = (ii /= 0)
3314 :
3315 : ! Read U matrices using real arrays.
3316 : ! TODO: This should be tested more carefully, especially when we are excluding bands.
3317 0 : ABI_MALLOC(u_mat, (2, nwan, nwan, nkbz))
3318 0 : NCF_CHECK(nf90_get_var(ncid, vid("U_matrix"), u_mat, start=[1,1,1,1,spin], count=[2, nwan, nwan, nkbz, 1]))
3319 0 : ABI_MALLOC(u_mat_opt, (2, num_bands, nwan, nkbz))
3320 0 : NCF_CHECK(nf90_get_var(ncid, vid("U_matrix_opt"), u_mat_opt, start=[1,1,1,1,spin], count=[2, num_bands, nwan, nkbz, 1]))
3321 :
3322 : ! Copy data: real --> complex
3323 0 : ABI_MALLOC(wan%u_mat, (num_bands, nwan, nkbz))
3324 0 : ABI_MALLOC(wan%u_mat_opt, (nwan, nwan, nkbz))
3325 0 : wan%u_mat = u_mat(1,:,:,:) + j_dpc * u_mat(2,:,:,:)
3326 0 : wan%u_mat_opt = u_mat_opt(1,:,:,:) + j_dpc * u_mat_opt(2,:,:,:)
3327 :
3328 0 : ABI_FREE(u_mat)
3329 0 : ABI_FREE(u_mat_opt)
3330 :
3331 0 : ABI_MALLOC(wan%centres, (3, nwan))
3332 0 : ABI_MALLOC(wan%spreads, (nwan))
3333 0 : NCF_CHECK(nf90_get_var(ncid, vid("wann_centres"), wan%centres, start=[1,1,spin], count=[3, nwan, 1]))
3334 0 : NCF_CHECK(nf90_get_var(ncid, vid("wann_spreads"), wan%spreads, start=[1,spin], count=[nwan]))
3335 :
3336 0 : ABI_MALLOC(int_1d, (num_bands))
3337 0 : NCF_CHECK(nf90_get_var(ncid, vid("band_in_int"), int_1d, start=[1,spin], count=[num_bands, 1]))
3338 0 : ABI_MALLOC(wan%band_in, (num_bands))
3339 0 : wan%band_in = (int_1d /= 0)
3340 0 : ABI_FREE(int_1d)
3341 :
3342 0 : ABI_MALLOC(int_2d, (num_bands, wan%nkbz))
3343 0 : NCF_CHECK(nf90_get_var(ncid, vid("lwindow_int"), int_2d, start=[1,1,spin], count=[num_bands, nkbz, 1]))
3344 0 : ABI_MALLOC(wan%lwindow, (wan%num_bands, wan%nkbz))
3345 0 : wan%lwindow = (int_2d /= 0)
3346 0 : ABI_FREE(int_2d)
3347 :
3348 0 : NCF_CHECK(nf90_get_var(ncid, vid("kptrlatt"), kptrlatt))
3349 0 : ABI_CHECK(isdiagmat(kptrlatt), "kptrlatt should be diagonal! Please use nkgpt with one shift")
3350 0 : wan%ngkpt = get_diag(kptrlatt)
3351 : !NCF_CHECK(nctk_get_dim(ncid, "nshiftk", wan%nshiftk))
3352 : !if nshiftk
3353 :
3354 : ! Read all KS eigenvalues and trasfer data to %eigs_w (note mband here)
3355 0 : ABI_MALLOC(wan%all_eigens, (mband, nkbz))
3356 0 : NCF_CHECK(nf90_get_var(ncid, vid("eigenvalues"), wan%all_eigens, start=[1,1,spin]))
3357 :
3358 0 : NCF_CHECK(nf90_close(ncid))
3359 : !end if ! master
3360 :
3361 : ! Compute dimwin, winstart, bmin and bmax from lwindow.
3362 0 : ABI_ICALLOC(wan%dimwin, (nkbz))
3363 0 : ABI_ICALLOC(wan%winstart, (nkbz))
3364 :
3365 0 : wan%bmin = huge(1); wan%bmax = -1
3366 0 : do ik=1,nkbz
3367 0 : do ib=1,wan%num_bands
3368 0 : if (wan%lwindow(ib, ik)) then
3369 0 : wan%dimwin(ik) = wan%dimwin(ik) + 1
3370 0 : if (wan%winstart(ik) == 0) wan%winstart(ik) = ib
3371 0 : wan%bmin = min(wan%bmin, ib)
3372 0 : wan%bmax = max(wan%bmax, ib)
3373 : end if
3374 : end do
3375 : end do
3376 :
3377 0 : call wan%krank%from_kptrlatt(nkbz, wan%kbz, kptrlatt, compute_invrank=.True.)
3378 :
3379 0 : ABI_MALLOC(wan%rmod_h, (wan%nr_h))
3380 0 : do ir=1,wan%nr_h
3381 0 : wan%rmod_h(ir) = sqrt(dot_product(wan%r_h(:,ir), matmul(cryst%rmet, wan%r_h(:,ir))))
3382 : end do
3383 :
3384 : ! Get total rotation matrix: the product of the optimal subspace x the rotation among the nwan Wannier functions.
3385 0 : ii = maxval(wan%dimwin)
3386 0 : ABI_CALLOC(wan%u_k, (ii, nwan, nkbz))
3387 0 : do ik=1,nkbz
3388 0 : wan%u_k(1:wan%dimwin(ik), 1:nwan, ik) = matmul(wan%u_mat_opt(1:wan%dimwin(ik), :, ik), wan%u_mat(:, 1:nwan, ik))
3389 : end do
3390 :
3391 0 : wan%keep_umats = keep_umats
3392 :
3393 0 : if (.not. keep_umats) then
3394 0 : ABI_FREE(wan%u_mat)
3395 0 : ABI_FREE(wan%u_mat_opt)
3396 : end if
3397 :
3398 : ! ====================================================
3399 : ! Build the Hamiltonian in the Wannier representation
3400 : ! ====================================================
3401 : !call wan%get_window_eig(et_opt)
3402 0 : nextbands = count(wan%exclude_bands /= 0)
3403 : !REAL(KIND = DP) :: et_opt(nbndep, nks)
3404 : ! KS eigenvalues within the outer window in the first dimwin(ik) entries
3405 0 : ii = wan%num_bands ! TODO: Check
3406 0 : ABI_MALLOC(et_opt, (ii ,nkbz))
3407 :
3408 : !print *, "nextbands", nextbands; print *, "exclude_bands:", wan%exclude_bands
3409 : !print *, "band_in:", wan%band_in; print *, "lwindow:", wan%lwindow
3410 :
3411 0 : if (nextbands /= 0) then
3412 0 : do ik=1,nkbz
3413 : jb = 0; mb = 0
3414 0 : do ib=1,wan%num_bands
3415 : !if (wan%exclude_bands(ib) /= 0) cycle
3416 0 : if (.not. wan%band_in(ib)) cycle
3417 0 : jb = jb + 1
3418 0 : if (wan%lwindow(jb, ik)) then
3419 0 : mb = mb + 1; et_opt(mb, ik) = wan%all_eigens(ib, ik)
3420 : end if
3421 : end do
3422 : end do
3423 :
3424 : else
3425 0 : do ik=1,nkbz
3426 0 : mb = 0
3427 0 : do ib=1,wan%dimwin(ik)
3428 0 : if (wan%lwindow(ib, ik)) then
3429 0 : mb = mb + 1; et_opt(mb, ik) = wan%all_eigens(ib, ik)
3430 : end if
3431 : end do
3432 : end do
3433 : end if
3434 :
3435 0 : ABI_CALLOC(chs, (nwan, nwan, nkbz))
3436 :
3437 0 : do ik=1,nkbz
3438 0 : do jb=1,nwan
3439 0 : do ib=1,jb
3440 0 : ctmp = czero
3441 0 : do mb=1,wan%dimwin(ik)
3442 0 : ctmp = ctmp + conjg(wan%u_k(mb, ib, ik)) * et_opt(mb, ik) * wan%u_k(mb, jb, ik)
3443 : end do
3444 0 : chs(ib, jb, ik) = ctmp
3445 0 : chs(jb, ib, ik) = conjg(ctmp)
3446 : end do
3447 : end do
3448 : end do ! ik
3449 0 : ABI_FREE(et_opt)
3450 :
3451 0 : ABI_CALLOC(chw, (nwan, nwan, nr_h))
3452 0 : do ir=1,nr_h
3453 0 : do ik=1,nkbz
3454 0 : chw(:,:,ir) = chw(:,:,ir) + chs(:,:,ik) * exp(-j_dpc * two_pi * dot_product(wan%kbz(:, ik), wan%r_h(:, ir))) / dble(nkbz)
3455 : end do
3456 : end do
3457 0 : ABI_FREE(chs)
3458 :
3459 : ! Now rearrange the data to have R_e in the first dimension.
3460 0 : ABI_CALLOC(wan%hwan_r, (nr_h, nwan, nwan))
3461 0 : do ir=1,nr_h
3462 0 : wan%hwan_r(ir,:,:) = chw(:,:,ir)
3463 : end do
3464 0 : ABI_FREE(chw)
3465 :
3466 : ! Write spatial decay to file.
3467 0 : if (my_rank == master .and. len_trim(out_prefix) > 0) then
3468 0 : out_path = strcat(out_prefix, "_spin", itoa(spin), "_HRWAN.txt")
3469 0 : if (open_file(out_path, msg, newunit=ount, form="formatted", action="write", status='unknown') /= 0) then
3470 0 : ABI_ERROR(msg)
3471 : end if
3472 0 : write(ount, "(a)")"# Decay of Hamiltonian in the Wannier representation"
3473 0 : write(ount, "(a)")"# |R| [Bohr] Max_{m,n} |H(R,m,n)| [Ha]"
3474 0 : do ir=1,nr_h
3475 0 : write(ount, *) wan%rmod_h(ir), maxval(abs(wan%hwan_r(ir,:,:)))
3476 : end do
3477 0 : close(ount)
3478 : end if
3479 :
3480 0 : call cryst%free()
3481 :
3482 : contains
3483 0 : integer function vid(var_name)
3484 : character(len=*),intent(in) :: var_name
3485 0 : vid = nctk_idname(ncid, var_name)
3486 : end function vid
3487 :
3488 : end subroutine wan_from_abiwan
3489 : !!***
3490 :
3491 : !!****f* m_mlwfovlp/wan_print
3492 : !! NAME
3493 : !! wan_print
3494 : !!
3495 : !! FUNCTION
3496 : !!
3497 : !! SOURCE
3498 :
3499 0 : subroutine wan_print(wan, units)
3500 :
3501 : !Arguments ------------------------------------
3502 : class(wan_t),intent(in) :: wan
3503 : integer,intent(in) :: units(:)
3504 :
3505 : !Local variables-------------------------------
3506 : type(yamldoc_t) :: ydoc
3507 : !************************************************************************
3508 :
3509 0 : ydoc = yamldoc_open("WANNIER_PARAMS")
3510 :
3511 : !call ydoc%add_string("method", "Gaussian")
3512 : call ydoc%add_ints("spin, nwan, num_bands, bmin, bmax, nkbz, nr_h, nr_e, nr_p", &
3513 0 : [wan%spin, wan%nwan, wan%num_bands, wan%bmin, wan%bmax, wan%nkbz, wan%nr_h, wan%nr_e, wan%nr_p])
3514 : !call ydoc%add_reals("nelect, wan_mesh_step_eV", &
3515 : ! [wan%nelect, wan%step * Ha_eV])
3516 : !call ydoc%add_real("", efermi * Ha_eV)
3517 :
3518 : ! Write header in Yaml format but prepend # so that one can still use tools such as gnuplot or xmgrace.
3519 0 : call ydoc%write_units_and_free(units)
3520 :
3521 0 : end subroutine wan_print
3522 : !!***
3523 :
3524 : !!****f* m_mlwfovlp/wan_interp_ham
3525 : !! NAME
3526 : !! wan_interp_ham
3527 : !!
3528 : !! FUNCTION
3529 : !! Interpolate the Hamiltonian at an arbitrary k-point and return the rotation
3530 : !! matrix. If vcart is present, also return diagonal Cartesian group velocities
3531 : !! obtained from the analytic derivative of the real-space Hamiltonian.
3532 : !! Velocities are in atomic units.
3533 : !!
3534 : !! TODO: This routine presently implements the legacy Wannier90 interpolation
3535 : !! obtained with use_ws_distance=.false. Modern Wannier90 versions enable
3536 : !! use_ws_distance by default and use pair-dependent translated lattice vectors
3537 : !! and degeneracies from w90_ws_distance:ws_translate_dist (irdist_ws and
3538 : !! wdist_ndeg). For an Al test, the legacy and default Wannier90 interpolations
3539 : !! differed by about 0.05 eV on average and up to 0.26 eV along a k-path.
3540 : !!
3541 : !! The current Wannier90 library interface does not return these arrays and
3542 : !! deallocates its internal parameters before wannier_run returns. A possible
3543 : !! implementation is to request write_hr=.true., read the resulting *_hr.dat
3544 : !! and *_wsvec.dat files after wannier_run, and store the Hamiltonian and
3545 : !! pair-dependent WS data in ABIWAN.nc for use here and by AbiPy.
3546 : !!
3547 : !! INPUTS
3548 : !! cryst: Crystal structure providing the dimensional primitive vectors.
3549 : !! kpt: Reduced coordinates of the interpolation point.
3550 : !!
3551 : !! OUTPUTS
3552 : !! uk_wan: Eigenvectors of the interpolated Hamiltonian, stored by columns.
3553 : !! eigens: Interpolated eigenvalues in Hartree.
3554 : !! vcart: Optional diagonal Cartesian group velocities in atomic units.
3555 : !! Individual values inside an exactly degenerate subspace depend on the
3556 : !! eigenvectors selected by the diagonalization.
3557 : !!
3558 : !! SOURCE
3559 :
3560 0 : subroutine wan_interp_ham(wan, cryst, kpt, uk_wan, eigens, vcart)
3561 :
3562 : !Arguments ------------------------------------
3563 : class(wan_t),intent(in) :: wan
3564 : class(crystal_t),intent(in) :: cryst
3565 : real(dp),intent(in) :: kpt(3)
3566 : real(dp),intent(out) :: eigens(wan%nwan)
3567 : complex(dp),intent(out) :: uk_wan(wan%nwan, wan%nwan)
3568 : real(dp),optional,intent(out) :: vcart(3, wan%nwan)
3569 :
3570 : !Local variables-------------------------------
3571 : integer :: ir, idir, ib
3572 0 : real(dp) :: rcart(3,wan%nr_h)
3573 0 : complex(dp) :: eikr(wan%nr_h),deikr(wan%nr_h), dham(wan%nwan, wan%nwan), vmat(wan%nwan, wan%nwan)
3574 : !************************************************************************
3575 :
3576 0 : do ir=1,wan%nr_h
3577 0 : eikr(ir) = exp(j_dpc * two_pi * dot_product(kpt, wan%r_h(:, ir))) / wan%ndegen_h(ir)
3578 : end do
3579 :
3580 : ! H_ij(k) = sum_R e^{+ik.R} * H_ij(R)
3581 0 : call ZGEMV("T", wan%nr_h, wan%nwan**2, cone, wan%hwan_r, wan%nr_h, eikr, 1, czero, uk_wan, 1)
3582 :
3583 : ! Hermitianize and diagonalize. xheev returns eigenvectors as columns.
3584 0 : uk_wan = half * (uk_wan + transpose(conjg(uk_wan)))
3585 0 : call xheev("V", "U", wan%nwan, uk_wan, eigens)
3586 :
3587 0 : if (present(vcart)) then
3588 0 : do ir=1,wan%nr_h
3589 0 : rcart(:,ir) = matmul(cryst%rprimd, real(wan%r_h(:,ir), kind=dp))
3590 : end do
3591 0 : do idir=1,3
3592 0 : do ir=1,wan%nr_h
3593 0 : deikr(ir) = j_dpc * rcart(idir,ir) * eikr(ir)
3594 : end do
3595 0 : call ZGEMV("T", wan%nr_h, wan%nwan**2, cone, wan%hwan_r, wan%nr_h, deikr, 1, czero, dham, 1)
3596 0 : dham = half * (dham + transpose(conjg(dham)))
3597 0 : vmat = matmul(transpose(conjg(uk_wan)), matmul(dham, uk_wan))
3598 0 : do ib=1,wan%nwan
3599 0 : vcart(idir,ib) = real(vmat(ib,ib), kind=dp)
3600 : end do
3601 : end do
3602 : end if
3603 :
3604 0 : end subroutine wan_interp_ham
3605 : !!***
3606 :
3607 : !!****f* m_mlwfovlp/wan_free
3608 : !! NAME
3609 : !! wan_free
3610 : !!
3611 : !! FUNCTION
3612 : !! Free dynamic memory.
3613 : !!
3614 : !! SOURCE
3615 :
3616 32 : subroutine wan_free(wan)
3617 :
3618 : !Arguments ------------------------------------
3619 : class(wan_t),intent(inout) :: wan
3620 : !************************************************************************
3621 :
3622 : ! integer
3623 32 : ABI_SFREE(wan%dimwin)
3624 32 : ABI_SFREE(wan%winstart)
3625 32 : ABI_SFREE(wan%exclude_bands)
3626 32 : ABI_SFREE(wan%r_h)
3627 32 : ABI_SFREE(wan%r_e)
3628 32 : ABI_SFREE(wan%r_p)
3629 :
3630 : ! real
3631 32 : ABI_SFREE(wan%ndegen_h)
3632 32 : ABI_SFREE(wan%ndegen_e)
3633 32 : ABI_SFREE(wan%ndegen_p)
3634 32 : ABI_SFREE(wan%rmod_h)
3635 32 : ABI_SFREE(wan%rmod_e)
3636 32 : ABI_SFREE(wan%rmod_p)
3637 32 : ABI_SFREE(wan%all_eigens)
3638 32 : ABI_SFREE(wan%centres)
3639 32 : ABI_SFREE(wan%spreads)
3640 32 : ABI_SFREE(wan%kbz)
3641 32 : ABI_SFREE(wan%band_in)
3642 32 : ABI_SFREE(wan%lwindow)
3643 :
3644 : ! Complex
3645 32 : ABI_SFREE(wan%u_mat)
3646 32 : ABI_SFREE(wan%u_mat_opt)
3647 32 : ABI_SFREE(wan%u_k)
3648 32 : ABI_SFREE(wan%hwan_r)
3649 32 : ABI_SFREE(wan%grpe_wwp)
3650 :
3651 32 : call wan%krank%free()
3652 :
3653 32 : end subroutine wan_free
3654 : !!***
3655 :
3656 : !!****f* m_mlwfovlp/wan_setup_eph_ws_kq
3657 : !! NAME
3658 : !! wan_setup_eph_ws_kq
3659 : !!
3660 : !! FUNCTION
3661 : !! Prepare the interpolation of the e-ph matrix elements.
3662 : !!
3663 : !! SOURCE
3664 :
3665 0 : subroutine wan_setup_eph_ws_kq(wan, cryst, shiftk, kptrlatt, qptrlatt, my_pert_start, my_npert, pert_comm)
3666 :
3667 : !Arguments ------------------------------------
3668 : class(wan_t),intent(inout) :: wan
3669 : type(crystal_t),intent(in) :: cryst
3670 : real(dp),intent(in) :: shiftk(3)
3671 : integer,intent(in) :: kptrlatt(3,3), qptrlatt(3,3), my_pert_start, my_npert
3672 : type(xcomm_t),target,intent(in) :: pert_comm
3673 :
3674 : !Local variables-------------------------------
3675 : integer,parameter :: lmax(3) = [2,2,2]
3676 : real(dp),parameter :: center(3) = zero
3677 : !************************************************************************
3678 :
3679 : ABI_UNUSED(shiftk)
3680 :
3681 : ! See Appendix A of Phys. Rev. Research 3, 043022 for possible improvements.
3682 0 : if (.not. allocated(wan%r_h)) then
3683 0 : call wigner_seitz(center, lmax, kptrlatt, cryst%rmet, wan%nr_h, wan%r_h, wan%ndegen_h, wan%rmod_h)
3684 : end if
3685 :
3686 0 : if (.not. allocated(wan%r_e)) then
3687 0 : call wigner_seitz(center, lmax, kptrlatt, cryst%rmet, wan%nr_e, wan%r_e, wan%ndegen_e, wan%rmod_e)
3688 : end if
3689 :
3690 0 : if (.not. allocated(wan%r_p)) then
3691 0 : call wigner_seitz(center, lmax, qptrlatt, cryst%rmet, wan%nr_p, wan%r_p, wan%ndegen_p, wan%rmod_p)
3692 : end if
3693 :
3694 0 : if (wan%my_pert_start == -1) then
3695 0 : wan%my_pert_start = my_pert_start; wan%my_npert = my_npert; wan%pert_comm => pert_comm
3696 : else
3697 0 : ABI_CHECK_IEQ(wan%my_pert_start, my_pert_start, "different values for my_pert_start")
3698 0 : ABI_CHECK_IEQ(wan%my_npert, my_npert, "different values for my_npert")
3699 0 : ABI_CHECK(associated(wan%pert_comm), "wan%per_comm is not associated!")
3700 0 : ABI_CHECK_IEQ(wan%pert_comm%value, pert_comm%value, "different values for pert_comm%value")
3701 : end if
3702 :
3703 : ! Allocate g in the Wannier representation.
3704 0 : ABI_CALLOC(wan%grpe_wwp, (wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, my_npert))
3705 :
3706 0 : end subroutine wan_setup_eph_ws_kq
3707 : !!***
3708 :
3709 : !!****f* m_mlwfovlp/wan_interp_eph_manyq
3710 : !! NAME
3711 : !! wan_interp_eph_manyq
3712 : !!
3713 : !! FUNCTION
3714 : !! Interpolate the e-ph matrix elements for one k-point and nq q-points.
3715 : !! Returns matrix elements in the atomic-representation and, optionally, the
3716 : !! eigenvalues obtained while diagonalizing the interpolated Hamiltonians.
3717 : !!
3718 : !! SOURCE
3719 :
3720 0 : subroutine wan_interp_eph_manyq(wan, cryst, nq, qpts, kpt, g_atm, out_eigens_k, out_eigens_kq)
3721 :
3722 : !Arguments ------------------------------------
3723 : class(wan_t),intent(in) :: wan
3724 : class(crystal_t),intent(in) :: cryst
3725 : integer,intent(in) :: nq
3726 : real(dp),intent(in) :: qpts(3,nq), kpt(3)
3727 : complex(dp),intent(out) :: g_atm(wan%nwan, wan%nwan, wan%my_npert, nq)
3728 : real(dp),optional,intent(out) :: out_eigens_k(wan%nwan), out_eigens_kq(wan%nwan,nq)
3729 :
3730 : !Local variables-------------------------------
3731 : integer :: ir, nr_e, nr_p, nwan, iq, my_npert, ipc, ncols_e, ncols_w
3732 : !arrays
3733 0 : real(dp) :: kq(3), eigens_k(wan%nwan), eigens_kq(wan%nwan)
3734 0 : complex(dp),allocatable :: eikr(:), eiqr(:), u_k(:,:), u_kq(:,:), cbuf_e(:,:,:,:), cbuf_w(:,:,:), cmat_w(:,:)
3735 : !************************************************************************
3736 :
3737 : ! TODO: Handle long-range part.
3738 0 : nr_p = wan%nr_p; nr_e = wan%nr_e; nwan = wan%nwan; my_npert = wan%my_npert
3739 :
3740 0 : ABI_MALLOC(eikr, (nr_e))
3741 0 : ABI_MALLOC(eiqr, (nr_p))
3742 0 : ABI_MALLOC(u_k, (nwan, nwan))
3743 0 : ABI_MALLOC(u_kq, (nwan, nwan))
3744 0 : ABI_MALLOC(cmat_w, (nwan, nwan))
3745 :
3746 0 : do ir=1,nr_e
3747 0 : eikr(ir) = exp(+j_dpc * two_pi * dot_product(kpt, wan%r_e(:, ir))) / wan%ndegen_e(ir)
3748 : end do
3749 0 : call wan%interp_ham(cryst, kpt, u_k, eigens_k)
3750 0 : if (present(out_eigens_k)) out_eigens_k = eigens_k
3751 :
3752 : ! grpe_wwp has shape: (nr_p, nr_e, nwan, nwan, my_npert))
3753 0 : ncols_e = nr_e * nwan **2 * my_npert
3754 0 : ABI_MALLOC(cbuf_e, (nr_e, nwan, nwan, my_npert))
3755 :
3756 0 : ncols_w = nwan ** 2 * my_npert
3757 0 : ABI_MALLOC(cbuf_w, (nwan, nwan, my_npert))
3758 :
3759 0 : do iq=1,nq
3760 0 : kq = kpt + qpts(:,iq)
3761 0 : call wan%interp_ham(cryst, kq, u_kq, eigens_kq)
3762 0 : if (present(out_eigens_kq)) out_eigens_kq(:,iq) = eigens_kq
3763 0 : do ir=1,nr_p
3764 0 : eiqr(ir) = exp(+j_dpc * two_pi * dot_product(qpts(:,iq), wan%r_p(:, ir))) / wan%ndegen_p(ir)
3765 : end do
3766 :
3767 : ! Transform along r_p
3768 : ! [Eqn. 22 of PRB 76, 165108 (2007)]
3769 : ! g~(R_e,q') = 1/ndegen(R_p) sum_R_p e^{iq'R_p} g(R_e,R_p)
3770 0 : call ZGEMV("T", nr_p, ncols_e, cone, wan%grpe_wwp, nr_p, eiqr, 1, czero, cbuf_e, 1)
3771 :
3772 : ! Transfor along r_e
3773 : ! [Eqn. 22 of PRB 76, 165108 (2007)]
3774 : ! g~(k',q') = 1/ndegen(R_e) sum_R_e e^{ik'R_e} g(R_e,q')
3775 0 : call ZGEMV("T", nr_e, ncols_w, cone, cbuf_e, nr_e, eikr, 1, czero, cbuf_w, 1)
3776 :
3777 : ! Go to Bloch space.
3778 : ! [Eqn. 22 of PRB 76, 165108 (2007)]
3779 : ! g(k',q') = U(k'+q') * g~(k',q') * U(k')^\dagger
3780 : !
3781 : ! RM - this is what is calculated
3782 : ! g(k',q') = U(k'+q')^\dagger * g~(k',q') * U(k')
3783 :
3784 : ! the two zgemm calls perform the following operations:
3785 : ! epmatf = [ cufkq * epmatf ] * cufkk^\dagger
3786 :
3787 0 : do ipc=1,my_npert
3788 0 : call ZGEMM('N', 'N', nwan, nwan, nwan, cone, u_kq, nwan, cbuf_w(:,:,ipc), nwan, czero, cmat_w, nwan)
3789 0 : call ZGEMM('N', 'C', nwan, nwan, nwan, cone, cmat_w, nwan, u_k, nwan, czero, g_atm(:,:,ipc,iq), nwan)
3790 : end do
3791 : end do ! iq
3792 :
3793 0 : ABI_FREE(cbuf_e)
3794 0 : ABI_FREE(cbuf_w)
3795 0 : ABI_FREE(eikr)
3796 0 : ABI_FREE(eiqr)
3797 0 : ABI_FREE(u_k)
3798 0 : ABI_FREE(u_kq)
3799 0 : ABI_FREE(cmat_w)
3800 :
3801 0 : end subroutine wan_interp_eph_manyq
3802 : !!***
3803 :
3804 : !!****f* m_mlwfovlp/wan_interp_eph_manyk
3805 : !! NAME
3806 : !! wan_interp_eph_manyk
3807 : !!
3808 : !! FUNCTION
3809 : !! Interpolate the e-ph matrix elements for one q-point and nk k-points.
3810 : !! The q-dependent R_p -> q Fourier transform is performed only once, then
3811 : !! reused for all k-points. Work arrays are allocated once per k batch.
3812 : !! Results are returned in the atomic-perturbation representation.
3813 : !!
3814 : !! SOURCE
3815 :
3816 0 : subroutine wan_interp_eph_manyk(wan, cryst, nk, kpts, qpt, g_atm, out_eigens_k, out_eigens_kq)
3817 :
3818 : !Arguments ------------------------------------
3819 : class(wan_t),intent(in) :: wan
3820 : class(crystal_t),intent(in) :: cryst
3821 : integer,intent(in) :: nk
3822 : real(dp),intent(in) :: kpts(3,nk), qpt(3)
3823 : complex(dp),intent(out) :: g_atm(wan%nwan, wan%nwan, wan%my_npert, nk)
3824 : real(dp),optional,intent(out) :: out_eigens_k(wan%nwan,nk), out_eigens_kq(wan%nwan,nk)
3825 :
3826 : !Local variables-------------------------------
3827 : integer :: ir, ik, ipc, nr_e, nr_p, nwan, my_npert, ncols_e, ncols_w
3828 0 : real(dp) :: kq(3), eigens_k(wan%nwan), eigens_kq(wan%nwan)
3829 0 : complex(dp),allocatable :: eikr(:), eiqr(:), u_k(:,:), u_kq(:,:), cbuf_e(:,:,:,:), cbuf_w(:,:,:), cmat_w(:,:)
3830 : !************************************************************************
3831 :
3832 0 : nr_p = wan%nr_p; nr_e = wan%nr_e; nwan = wan%nwan; my_npert = wan%my_npert
3833 :
3834 0 : ABI_MALLOC(eikr, (nr_e))
3835 0 : ABI_MALLOC(eiqr, (nr_p))
3836 0 : ABI_MALLOC(u_k, (nwan, nwan))
3837 0 : ABI_MALLOC(u_kq, (nwan, nwan))
3838 0 : ABI_MALLOC(cmat_w, (nwan, nwan))
3839 0 : ABI_MALLOC(cbuf_e, (nr_e, nwan, nwan, my_npert))
3840 0 : ABI_MALLOC(cbuf_w, (nwan, nwan, my_npert))
3841 :
3842 : ! Prepare the q-dependent vertex once for the entire k batch:
3843 : ! g(R_e,q) = sum_R_p exp(i q.R_p) g(R_e,R_p) / ndegen(R_p).
3844 0 : do ir=1,nr_p
3845 0 : eiqr(ir) = exp(+j_dpc * two_pi * dot_product(qpt, wan%r_p(:,ir))) / wan%ndegen_p(ir)
3846 : end do
3847 0 : ncols_e = nr_e * nwan ** 2 * my_npert
3848 0 : call ZGEMV("T", nr_p, ncols_e, cone, wan%grpe_wwp, nr_p, eiqr, 1, czero, cbuf_e, 1)
3849 :
3850 0 : ncols_w = nwan ** 2 * my_npert
3851 0 : do ik=1,nk
3852 0 : do ir=1,nr_e
3853 0 : eikr(ir) = exp(+j_dpc * two_pi * dot_product(kpts(:,ik), wan%r_e(:,ir))) / wan%ndegen_e(ir)
3854 : end do
3855 :
3856 0 : call wan%interp_ham(cryst, kpts(:,ik), u_k, eigens_k)
3857 0 : kq = kpts(:,ik) + qpt
3858 0 : call wan%interp_ham(cryst, kq, u_kq, eigens_kq)
3859 0 : if (present(out_eigens_k)) out_eigens_k(:,ik) = eigens_k
3860 0 : if (present(out_eigens_kq)) out_eigens_kq(:,ik) = eigens_kq
3861 :
3862 : ! g(k,q) in the Wannier gauge.
3863 0 : call ZGEMV("T", nr_e, ncols_w, cone, cbuf_e, nr_e, eikr, 1, czero, cbuf_w, 1)
3864 :
3865 : ! Rotate from the Wannier gauge to the interpolated electronic eigenstates.
3866 0 : do ipc=1,my_npert
3867 0 : call ZGEMM('N', 'N', nwan, nwan, nwan, cone, u_kq, nwan, cbuf_w(:,:,ipc), nwan, czero, cmat_w, nwan)
3868 0 : call ZGEMM('N', 'C', nwan, nwan, nwan, cone, cmat_w, nwan, u_k, nwan, czero, g_atm(:,:,ipc,ik), nwan)
3869 : end do
3870 : end do
3871 :
3872 0 : ABI_FREE(eikr)
3873 0 : ABI_FREE(eiqr)
3874 0 : ABI_FREE(u_k)
3875 0 : ABI_FREE(u_kq)
3876 0 : ABI_FREE(cmat_w)
3877 0 : ABI_FREE(cbuf_e)
3878 0 : ABI_FREE(cbuf_w)
3879 :
3880 0 : end subroutine wan_interp_eph_manyk
3881 : !!***
3882 :
3883 : !!****f* m_mlwfovlp/wan_prepare_eph_q
3884 : !! NAME
3885 : !! wan_prepare_eph_q
3886 : !!
3887 : !! FUNCTION
3888 : !! Fourier transform g(R_e,R_p) along R_p for one q-point. The result can be
3889 : !! reused by multiple calls to wan_interp_eph_manyk_from_q with bounded k
3890 : !! batches.
3891 : !!
3892 : !! SOURCE
3893 :
3894 0 : subroutine wan_prepare_eph_q(wan, qpt, g_req)
3895 :
3896 : class(wan_t),intent(in) :: wan
3897 : real(dp),intent(in) :: qpt(3)
3898 : complex(dp),intent(out) :: g_req(wan%nr_e, wan%nwan, wan%nwan, wan%my_npert)
3899 :
3900 : integer :: ir, ncols_e
3901 0 : complex(dp) :: eiqr(wan%nr_p)
3902 : !************************************************************************
3903 :
3904 0 : do ir=1,wan%nr_p
3905 0 : eiqr(ir) = exp(+j_dpc * two_pi * dot_product(qpt, wan%r_p(:,ir))) / wan%ndegen_p(ir)
3906 : end do
3907 0 : ncols_e = wan%nr_e * wan%nwan ** 2 * wan%my_npert
3908 0 : call ZGEMV("T", wan%nr_p, ncols_e, cone, wan%grpe_wwp, wan%nr_p, eiqr, 1, czero, g_req, 1)
3909 :
3910 0 : end subroutine wan_prepare_eph_q
3911 : !!***
3912 :
3913 : !!****f* m_mlwfovlp/wan_interp_eph_manyk_from_q
3914 : !! NAME
3915 : !! wan_interp_eph_manyk_from_q
3916 : !!
3917 : !! FUNCTION
3918 : !! Interpolate the e-ph matrix elements for a k-point batch from a vertex
3919 : !! previously transformed to q by wan_prepare_eph_q.
3920 : !!
3921 : !! SOURCE
3922 :
3923 0 : subroutine wan_interp_eph_manyk_from_q(wan, cryst, nk, kpts, qpt, g_req, g_atm, out_eigens_k, out_eigens_kq)
3924 :
3925 : class(wan_t),intent(in) :: wan
3926 : class(crystal_t),intent(in) :: cryst
3927 : integer,intent(in) :: nk
3928 : real(dp),intent(in) :: kpts(3,nk), qpt(3)
3929 : complex(dp),intent(in) :: g_req(wan%nr_e, wan%nwan, wan%nwan, wan%my_npert)
3930 : complex(dp),intent(out) :: g_atm(wan%nwan, wan%nwan, wan%my_npert, nk)
3931 : real(dp),optional,intent(out) :: out_eigens_k(wan%nwan,nk), out_eigens_kq(wan%nwan,nk)
3932 :
3933 : integer :: ir, ik, ipc, ncols_w
3934 0 : real(dp) :: kq(3), eigens_k(wan%nwan), eigens_kq(wan%nwan)
3935 0 : complex(dp),allocatable :: eikr(:), u_k(:,:), u_kq(:,:), g_wan(:,:,:), cmat_w(:,:)
3936 : !************************************************************************
3937 :
3938 0 : ABI_MALLOC(eikr, (wan%nr_e))
3939 0 : ABI_MALLOC(u_k, (wan%nwan, wan%nwan))
3940 0 : ABI_MALLOC(u_kq, (wan%nwan, wan%nwan))
3941 0 : ABI_MALLOC(g_wan, (wan%nwan, wan%nwan, wan%my_npert))
3942 0 : ABI_MALLOC(cmat_w, (wan%nwan, wan%nwan))
3943 :
3944 0 : ncols_w = wan%nwan ** 2 * wan%my_npert
3945 0 : do ik=1,nk
3946 0 : do ir=1,wan%nr_e
3947 0 : eikr(ir) = exp(+j_dpc * two_pi * dot_product(kpts(:,ik), wan%r_e(:,ir))) / wan%ndegen_e(ir)
3948 : end do
3949 :
3950 0 : call wan%interp_ham(cryst, kpts(:,ik), u_k, eigens_k)
3951 0 : kq = kpts(:,ik) + qpt
3952 0 : call wan%interp_ham(cryst, kq, u_kq, eigens_kq)
3953 0 : if (present(out_eigens_k)) out_eigens_k(:,ik) = eigens_k
3954 0 : if (present(out_eigens_kq)) out_eigens_kq(:,ik) = eigens_kq
3955 :
3956 0 : call ZGEMV("T", wan%nr_e, ncols_w, cone, g_req, wan%nr_e, eikr, 1, czero, g_wan, 1)
3957 0 : do ipc=1,wan%my_npert
3958 : call ZGEMM('N', 'N', wan%nwan, wan%nwan, wan%nwan, cone, u_kq, wan%nwan, &
3959 0 : g_wan(:,:,ipc), wan%nwan, czero, cmat_w, wan%nwan)
3960 : call ZGEMM('N', 'C', wan%nwan, wan%nwan, wan%nwan, cone, cmat_w, wan%nwan, &
3961 0 : u_k, wan%nwan, czero, g_atm(:,:,ipc,ik), wan%nwan)
3962 : end do
3963 : end do
3964 :
3965 0 : ABI_FREE(eikr)
3966 0 : ABI_FREE(u_k)
3967 0 : ABI_FREE(u_kq)
3968 0 : ABI_FREE(g_wan)
3969 0 : ABI_FREE(cmat_w)
3970 :
3971 0 : end subroutine wan_interp_eph_manyk_from_q
3972 : !!***
3973 :
3974 : !!****f* m_mlwfovlp/wan_eph_kbatch_size
3975 : !! NAME
3976 : !! wan_eph_kbatch_size
3977 : !!
3978 : !! FUNCTION
3979 : !! Return a k-batch size that bounds the combined complex local-perturbation
3980 : !! and complete atomic-perturbation buffers. memory_mb defaults to 64 MiB.
3981 : !!
3982 : !! SOURCE
3983 :
3984 0 : integer function wan_eph_kbatch_size(wan, nk, natom3, memory_mb) result(nk_batch)
3985 :
3986 : class(wan_t),intent(in) :: wan
3987 : integer,intent(in) :: nk, natom3
3988 : real(dp),optional,intent(in) :: memory_mb
3989 :
3990 : real(dp) :: limit_mb
3991 : !************************************************************************
3992 :
3993 0 : limit_mb = 64.0_dp
3994 0 : if (present(memory_mb)) limit_mb = memory_mb
3995 0 : ABI_CHECK(limit_mb > zero, "The Wannier e-ph k-batch memory limit must be positive")
3996 : nk_batch = int(limit_mb * 1024.0_dp ** 2 / &
3997 0 : (16.0_dp * wan%nwan ** 2 * (wan%my_npert + natom3)))
3998 0 : nk_batch = max(1, min(nk, nk_batch))
3999 :
4000 0 : end function wan_eph_kbatch_size
4001 : !!***
4002 :
4003 : !!****f* m_mlwfovlp/wan_ncwrite_gwan
4004 : !! NAME
4005 : !! wan_ncwrite_gwan
4006 : !!
4007 : !! FUNCTION
4008 : !! Write the e-ph matrix elements in the Wannier representation g(R_e,R_ph) to the GWAN.nc netcdf file.
4009 : !!
4010 : !! SOURCE
4011 :
4012 0 : subroutine wan_ncwrite_gwan(wan, dtfil, cryst, ebands, pert_comm)
4013 :
4014 : !Arguments ------------------------------------
4015 : class(wan_t),target,intent(in) :: wan
4016 : type(datafiles_type),intent(in) :: dtfil
4017 : type(crystal_t),intent(in) :: cryst
4018 : type(ebands_t),intent(in) :: ebands
4019 : type(xcomm_t),intent(in) :: pert_comm
4020 :
4021 : !Local variables-------------------------------
4022 : !scalars
4023 : integer :: spin, root_ncid, spin_ncid, ncerr, natom3, ount, ir, var_id, units(2), batch_size, idat, ndat
4024 0 : real(dp), contiguous, pointer :: rpt_d4(:,:,:,:), rpt_d6(:,:,:,:,:,:)
4025 : character(len=fnlen) :: gwan_filepath, txt_path
4026 0 : complex(dp),target,allocatable :: cbuf5(:,:,:,:,:)
4027 : character(len=500) :: msg
4028 : !************************************************************************
4029 :
4030 0 : units = [std_out, ab_out]
4031 0 : spin = wan%spin; natom3 = 3 * cryst%natom
4032 :
4033 0 : gwan_filepath = strcat(dtfil%filnam_ds(4), "_GWAN.nc")
4034 0 : call wrtout(units, sjoin("- Writing e-ph matrix elements in the Wannier representation to file:", gwan_filepath))
4035 :
4036 0 : if (spin == 1) then
4037 0 : NCF_CHECK(nctk_open_create(root_ncid, gwan_filepath, pert_comm%value))
4038 0 : NCF_CHECK(cryst%ncwrite(root_ncid))
4039 : !NCF_CHECK(hdr%ncwrite(root_ncid))
4040 0 : NCF_CHECK(ebands%ncwrite(root_ncid))
4041 : else
4042 0 : NCF_CHECK(nctk_open_modify(root_ncid, gwan_filepath, pert_comm%value))
4043 : end if
4044 :
4045 : ! Create group for this spin and define dimensions.
4046 0 : NCF_CHECK(nf90_def_grp(root_ncid, strcat("gwan", "_spin", itoa(spin)), spin_ncid))
4047 :
4048 : ncerr = nctk_def_dims(spin_ncid, [ &
4049 : nctkdim_t("nwan", wan%nwan), &
4050 : nctkdim_t("natom3", natom3), &
4051 : nctkdim_t("nr_h", wan%nr_h), &
4052 : nctkdim_t("nr_e", wan%nr_e), &
4053 : nctkdim_t("nr_p", wan%nr_p) &
4054 0 : ], defmode=.True.)
4055 0 : NCF_CHECK(ncerr)
4056 :
4057 : ncerr = nctk_def_arrays(spin_ncid, [ &
4058 : nctkarr_t("r_h", "dp", "three, nr_h"), &
4059 : nctkarr_t("r_e", "dp", "three, nr_e"), &
4060 : nctkarr_t("r_p", "dp", "three, nr_p"), &
4061 : nctkarr_t("ndegen_h", "int", "nr_h"), &
4062 : nctkarr_t("ndegen_e", "int", "nr_e"), &
4063 : nctkarr_t("ndegen_p", "int", "nr_p"), &
4064 : nctkarr_t("hwan_r", "dp", "two, nr_h, nwan, nwan"), &
4065 : nctkarr_t("grpe_wwp", "dp", "two, nr_p, nr_e, nwan, nwan, natom3") &
4066 0 : ])
4067 0 : NCF_CHECK(ncerr)
4068 :
4069 : ! Write data.
4070 0 : NCF_CHECK(nctk_set_datamode(spin_ncid))
4071 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("r_h"), wan%r_h))
4072 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("r_e"), wan%r_e))
4073 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("r_p"), wan%r_p))
4074 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ndegen_h"), wan%ndegen_h))
4075 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ndegen_e"), wan%ndegen_e))
4076 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("ndegen_p"), wan%ndegen_p))
4077 :
4078 0 : call c_f_pointer(c_loc(wan%hwan_r), rpt_d4, [2, wan%nr_h, wan%nwan, wan%nwan])
4079 0 : NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("hwan_r"), rpt_d4))
4080 :
4081 : ! Take into account that the array might be distributed over perturbations.
4082 : !NCF_CHECK(nctk_set_collective(spin_ncid, vid_spin("foo")))
4083 0 : call c_f_pointer(c_loc(wan%grpe_wwp), rpt_d6, [2, wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert])
4084 : ncerr = nf90_put_var(spin_ncid, vid_spin("grpe_wwp"), rpt_d6, &
4085 : start=[1,1,1,1,1,wan%my_pert_start], &
4086 0 : count=[2, wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert])
4087 0 : NCF_CHECK(ncerr)
4088 0 : NCF_CHECK(nf90_close(root_ncid))
4089 :
4090 0 : call xmpi_barrier(pert_comm%value)
4091 :
4092 : ! Check spatial decay of the EP matrix elements in the wannier basis
4093 : ! We plot: R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)|
4094 0 : if (pert_comm%me == 0) then
4095 0 : NCF_CHECK(nctk_open_read(root_ncid, gwan_filepath, xmpi_comm_self))
4096 : ! Get group for this spin.
4097 0 : NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("gwan", "_spin", itoa(spin)), spin_ncid))
4098 :
4099 0 : txt_path = strcat(dtfil%filnam_ds(4), "_spin", itoa(spin), "_GWAN.txt")
4100 0 : if (open_file(txt_path, msg, newunit=ount, form="formatted", action="write", status='unknown') /= 0) then
4101 0 : ABI_ERROR(msg)
4102 : end if
4103 0 : write(ount, '(a)') '# R_e [Bohr] max_{m,n,nu} |g(m,n,nu R_e,:)| [Ha/Bohr] '
4104 0 : var_id = vid_spin("grpe_wwp")
4105 : ! Perform IO in batches to keep memory at bay.
4106 : ! FIXME
4107 0 : batch_size = 1
4108 0 : ABI_MALLOC(cbuf5, (wan%nr_p, batch_size, wan%nwan, wan%nwan, natom3))
4109 0 : call c_f_pointer(c_loc(cbuf5), rpt_d6, [2, wan%nr_p, batch_size, wan%nwan, wan%nwan, natom3])
4110 0 : do ir=1,wan%nr_e, batch_size
4111 0 : ndat = blocked_loop(ir, wan%nr_e, batch_size)
4112 : !nctkarr_t("grpe_wwp", "dp", "two, nr_p, nr_e, nwan, nwan, natom3") &
4113 : ncerr = nf90_get_var(spin_ncid, var_id, rpt_d6, &
4114 0 : start=[1,1,ir,1,1,1], count=[2, wan%nr_p, batch_size, wan%nwan, wan%nwan, natom3])
4115 0 : NCF_CHECK(ncerr)
4116 0 : do idat=1,ndat
4117 0 : write(ount, *) wan%rmod_e(ir+idat-1), maxval(abs(cbuf5(:,idat,:,:,:)))
4118 : end do
4119 : end do
4120 0 : ABI_FREE(cbuf5)
4121 :
4122 0 : close(ount)
4123 0 : NCF_CHECK(nf90_close(root_ncid))
4124 : end if
4125 :
4126 : contains
4127 0 : integer function vid_spin(var_name)
4128 : character(len=*),intent(in) :: var_name
4129 0 : vid_spin = nctk_idname(spin_ncid, var_name)
4130 : end function vid_spin
4131 :
4132 : end subroutine wan_ncwrite_gwan
4133 : !!***
4134 :
4135 : !!****f* m_mlwfovlp/wan_load_gwan
4136 : !! NAME
4137 : !! wan_load_gwan
4138 : !!
4139 : !! FUNCTION
4140 : !! Read g(R_e, R_p) in the Wannier representation from an ABIWAN.nc file
4141 : !!
4142 : !! SOURCE
4143 :
4144 0 : subroutine wan_load_gwan(wan, gwan_filepath, cryst, spin, nsppol, all_comm)
4145 :
4146 : !Arguments ------------------------------------
4147 : class(wan_t),target,intent(inout) :: wan
4148 : character(len=*),intent(in) :: gwan_filepath
4149 : integer,intent(in) :: spin, nsppol
4150 : type(crystal_t),intent(in) :: cryst
4151 : type(xcomm_t),intent(in) :: all_comm ! , pert_comm
4152 :
4153 : !Local variables-------------------------------
4154 : !scalars
4155 : integer :: root_ncid, spin_ncid, ncerr, units(2)
4156 : real(dp) :: cpu, wall, gflops
4157 : logical,parameter :: keep_umats = .False.
4158 0 : type(crystal_t) :: gwan_cryst
4159 0 : real(dp), contiguous, pointer :: rpt_d6(:,:,:,:,:,:) !, rpt_d4(:,:,:,:)
4160 : !************************************************************************
4161 :
4162 0 : units = [std_out, ab_out]
4163 0 : call cwtime(cpu, wall, gflops, "start")
4164 0 : if (nsppol == 2) then
4165 0 : call wrtout(units, sjoin(" Reading g(R_e, R_p) for spin:", itoa(spin), " from GWAN file:", gwan_filepath))
4166 : else
4167 0 : call wrtout(units, sjoin(" Reading g(R_e, R_p) from GWAN file:", gwan_filepath))
4168 : end if
4169 :
4170 0 : NCF_CHECK(nctk_open_read(root_ncid, gwan_filepath, all_comm%value))
4171 :
4172 0 : call gwan_cryst%ncread(root_ncid)
4173 0 : if (cryst%compare(gwan_cryst, header=" Comparing input crystal with GWAN crystal.") /= 0) then
4174 0 : ABI_ERROR("Crystal structure from input and GWAN file do not agree! Check messages above!")
4175 : end if
4176 0 : call gwan_cryst%free()
4177 :
4178 : ! Get netcdf group for this spin.
4179 0 : NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("gwan", "_spin", itoa(spin)), spin_ncid))
4180 :
4181 : ! Read supercell lattice vectors.
4182 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "nr_e", wan%nr_e))
4183 0 : NCF_CHECK(nctk_get_dim(spin_ncid, "nr_p", wan%nr_p))
4184 :
4185 0 : ABI_MALLOC(wan%r_e, (3, wan%nr_e))
4186 0 : ABI_MALLOC(wan%r_p, (3, wan%nr_p))
4187 0 : ABI_MALLOC(wan%ndegen_e, (wan%nr_e))
4188 0 : ABI_MALLOC(wan%ndegen_p, (wan%nr_p))
4189 :
4190 0 : NCF_CHECK(nf90_get_var(spin_ncid, vid_spin("r_e"), wan%r_e))
4191 0 : NCF_CHECK(nf90_get_var(spin_ncid, vid_spin("r_p"), wan%r_p))
4192 0 : NCF_CHECK(nf90_get_var(spin_ncid, vid_spin("ndegen_e"), wan%ndegen_e))
4193 0 : NCF_CHECK(nf90_get_var(spin_ncid, vid_spin("ndegen_p"), wan%ndegen_p))
4194 :
4195 : !call c_f_pointer(c_loc(wan%hwan_r), rpt_d4, [2, wan%nr_h, wan%nwan, wan%nwan])
4196 : !NCF_CHECK(nf90_get_var(spin_ncid, vid_spin("hwan_r"), rpt_d4))
4197 :
4198 : ! TODO:
4199 : ! Take into account that the array might be distributed over perturbations.
4200 : !print *, "wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert", wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert
4201 0 : ABI_MALLOC(wan%grpe_wwp, (wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert))
4202 0 : call c_f_pointer(c_loc(wan%grpe_wwp), rpt_d6, [2, wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert])
4203 :
4204 : ! nctkarr_t("grpe_wwp", "dp", "two, nr_p, nr_e, nwan, nwan, natom3") &
4205 0 : if (all_comm%nproc > 1) then
4206 0 : NCF_CHECK(nctk_set_collective(spin_ncid, vid_spin("grpe_wwp")))
4207 : end if
4208 : ncerr = nf90_get_var(spin_ncid, vid_spin("grpe_wwp"), rpt_d6, &
4209 : start=[1, 1, 1, 1, 1, wan%my_pert_start], &
4210 0 : count=[2, wan%nr_p, wan%nr_e, wan%nwan, wan%nwan, wan%my_npert])
4211 0 : NCF_CHECK(ncerr)
4212 :
4213 0 : NCF_CHECK(nf90_close(root_ncid))
4214 0 : call wrtout(units, " Reading of GWAN.nc file completed.")
4215 0 : call cwtime_report(" Reading Wannier e-ph matrix elements from GWAN.nc", cpu, wall, gflops)
4216 :
4217 : contains
4218 0 : integer function vid_spin(var_name)
4219 : character(len=*),intent(in) :: var_name
4220 0 : vid_spin = nctk_idname(spin_ncid, var_name)
4221 : end function vid_spin
4222 : end subroutine wan_load_gwan
4223 : !!***
4224 :
4225 : !!****f* m_mlwfovlp/wan_interp_ebands
4226 : !! NAME
4227 : !! wan_interp_ebands
4228 : !!
4229 : !! FUNCTION
4230 : !! Build new ebands_t object on a k-mesh via Wannier interpolation.
4231 : !!
4232 : !! INPUT
4233 : !! cryst<crystal_t> = Crystalline structure.
4234 : !! intp_kptrlatt(3,3) = New k-mesh
4235 : !! intp_nshiftk= Number of shifts in new k-mesh.
4236 : !! intp_shiftk(3,intp_nshiftk) = Shifts in new k-mesh.
4237 : !! band_block(2)=Initial and final band index. If [0,0], all bands are used
4238 : !!
4239 : !! OUTPUT
4240 : !! out_ebands: object with interpolated energies.
4241 : !!
4242 : !! NOTES
4243 : !! Fermi level of the interpolated bands is not recomputed by this routine
4244 : !! but copied from in_ebands.
4245 : !!
4246 : !! SOURCE
4247 :
4248 0 : subroutine wan_interp_ebands(wan_spin, cryst, in_ebands, intp_kptrlatt, intp_nshiftk, intp_shiftk, out_ebands, comm)
4249 :
4250 : !Arguments ------------------------------------
4251 : type(ebands_t),intent(in) :: in_ebands
4252 : type(wan_t),intent(in) :: wan_spin(in_ebands%nsppol)
4253 : type(crystal_t),intent(in) :: cryst
4254 : integer,intent(in) :: intp_kptrlatt(3,3), intp_nshiftk, comm
4255 : real(dp),intent(in) :: intp_shiftk(3,intp_nshiftk)
4256 : type(ebands_t),intent(out) :: out_ebands
4257 :
4258 : !Local variables-------------------------------
4259 : !scalars
4260 : integer :: spin, ik, nwan, ierr, cnt, my_rank, nproc
4261 : real(dp) :: cpu, wall, gflops
4262 : !arrays
4263 : integer :: band_block(2)
4264 : real(dp) :: params(4)
4265 0 : real(dp),allocatable :: eigens_k(:)
4266 0 : complex(dp),allocatable :: u_k(:,:)
4267 : !************************************************************************
4268 :
4269 0 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
4270 :
4271 0 : call cwtime(cpu, wall, gflops, "start")
4272 :
4273 : ! Build new ebands object with memory to be filled.
4274 0 : band_block(:) = [1, wan_spin(1)%max_nwan]
4275 : out_ebands = in_ebands%interp_kmesh(cryst, params, intp_kptrlatt, intp_nshiftk, intp_shiftk, &
4276 0 : band_block, comm, malloc_only=.True.)
4277 0 : out_ebands%eig = zero
4278 :
4279 : cnt = 0
4280 0 : do spin=1,in_ebands%nsppol
4281 0 : associate (wan => wan_spin(spin))
4282 0 : nwan = wan%nwan
4283 0 : ABI_MALLOC(u_k, (nwan, nwan))
4284 0 : ABI_MALLOC(eigens_k, (nwan))
4285 0 : do ik=1,out_ebands%nkpt
4286 0 : cnt = cnt + 1; if (mod(cnt - 1, nproc) /= my_rank) cycle ! MPI parallelism inside comm.
4287 0 : call wan%interp_ham(cryst, out_ebands%kptns(:,ik), u_k, eigens_k)
4288 0 : out_ebands%eig(1:nwan, ik, spin) = eigens_k
4289 : end do ! ik
4290 0 : ABI_FREE(u_k)
4291 0 : ABI_FREE(eigens_k)
4292 : end associate
4293 : end do ! spin
4294 :
4295 0 : call xmpi_sum(out_ebands%eig, comm, ierr)
4296 :
4297 : ! Copy Fermi energies.
4298 0 : out_ebands%fermie = in_ebands%fermie
4299 0 : out_ebands%fermih = in_ebands%fermih
4300 :
4301 0 : call cwtime_report(" Wannier interpolation of electronic bands", cpu, wall, gflops)
4302 :
4303 0 : end subroutine wan_interp_ebands
4304 : !!***
4305 :
4306 9909192 : end module m_mlwfovlp
4307 : !!***
|
