Branch data 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 : : !!***
|
