Line data Source code
1 : !!****m* ABINIT/m_abstract_wf
2 : !! NAME
3 : !! m_abstract_wf
4 : !!
5 : !! FUNCTION
6 : !! Interface with Wannier90.
7 : !! This module contains the abstract type abstract_wf and its children.
8 : !! The abstract_wf type is used to store the wavefunctions either in wfd or wfk format.
9 : !! And it provide the common interface for both format.
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2005-2026 ABINIT group (BAmadon, CEspejo, FJollet, TRangel, DRH)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public License, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !!
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_abstract_wf
25 :
26 : use defs_basis
27 : use defs_wannier90
28 : use m_abicore
29 : use m_errors
30 : use m_atomdata
31 : use m_xmpi
32 : use m_sort
33 : #ifdef FC_NAG
34 : use f90_unix_dir
35 : #endif
36 : use netcdf
37 : use m_nctk
38 : use m_hdr
39 : use m_dtset
40 : use m_dtfil
41 :
42 : use m_build_info, only : abinit_version
43 : use defs_wvltypes, only : wvl_internal_type
44 : use defs_datatypes, only : pseudopotential_type
45 : use defs_abitypes, only : MPI_type
46 : use m_io_tools, only : delete_file
47 : use m_fstrings, only : strcat, sjoin, itoa
48 : use m_fftcore, only : ngfft_seq, get_kg
49 : use m_crystal, only : crystal_t
50 : use m_ebands, only : ebands_t
51 : use m_pawtab, only : pawtab_type
52 : use m_pawrhoij, only : pawrhoij_copy
53 : use m_pawcprj, only : pawcprj_type
54 : use m_wfd, only : wfd_t, wave_t, WFD_STORED
55 :
56 : implicit none
57 :
58 : private
59 :
60 :
61 : !!***
62 :
63 :
64 : !type wann_ksetting_t:
65 : ! This type is used to store the information about the kpoints
66 : ! and the MPI distribution of the kpoints.
67 : ! It is used in the wannier90 interface.
68 : type, public:: wann_ksetting_t
69 : logical :: has_ovikp = .False. !.True. if ovikp is allocated
70 : type(crystal_t), pointer :: cryst => null() ! crystal structure
71 : !nkpt: number of kpoints
72 : !mband: number of bands
73 : !num_nnmax: maximum number of nearest neighbors
74 : !nntot: total number of nearest neighbors
75 : !nsppol: number of spin polarizations
76 : integer :: nkpt=0, mband=0, num_nnmax=0, nntot=0, nsppol=0
77 : !rank: rank of the current process
78 : !comm: communicator
79 : !nprocs: number of processes
80 : integer :: rank=-999, comm=-999, nprocs=-999
81 : !ovikp: array of size (nkpt, num_nnmax) containing the indices of the nearest neighbors
82 : integer, allocatable :: ovikp(:, :)
83 : !my_nspin: number of spin polarizations for the current process
84 : !my_nkpt: number of kpoints for the current process
85 : !my_nkpt_pnn: number of kpoints for the current process including the nearest neighbors
86 : integer :: my_nspin, my_nkpt, my_nkpt_pnn
87 : !my_spins: array of size (my_nspin) containing the spin indices for the current process
88 : !my_ikpts: array of size (my_nkpt) containing the indices of the kpoints for the current process
89 : !my_ikpts_pnn: array of size (my_nkpt_pnn) containing the indices of the kpoints for the current process including the nearest neighbors
90 : integer, allocatable :: my_spins(:), my_ikpts(:), my_ikpts_pnn(:)
91 : !kkpts: array of size (nkpt, 3) containing the kpoints
92 : !my_kkpts: array of size (my_nkpt, 3) containing the kpoints for the current process
93 : !my_kkpts_pnn: array of size (my_nkpt_pnn, 3) containing the kpoints for the current process including the nearest neighbors
94 : real(dp), pointer :: kkpts(:, :)=> null()
95 : real(dp), allocatable :: my_kkpts(:, :), my_kkpts_pnn(:, :)
96 : contains
97 : procedure :: init => wann_ksetting_init
98 : procedure :: set_ovikp => wann_ksetting_set_ovikp
99 : procedure :: free => wann_ksetting_free
100 : procedure :: get_mpw_gmax => wann_ksetting_get_mpw_gmax
101 : procedure :: distribute_mpi => wann_ksetting_distribute_mpi
102 : procedure :: get_bks_mask => wann_ksetting_get_bks_mask
103 : end type wann_ksetting_t
104 :
105 :
106 :
107 : !type abstract_wf:
108 : ! This type is used to store the information about the wavefunctions
109 : ! and the MPI distribution of the wavefunctions.
110 : ! It provide the common interface for the wavefunctions.
111 : type, public :: abstract_wf
112 : !has_paw: .True. if the wavefunctions are PAW wavefunctions
113 : logical :: has_paw = .True.
114 : !cryst: crystal structure
115 : type(crystal_t), pointer :: cryst => null()
116 : type(datafiles_type),pointer :: dtfil => null()
117 : type(dataset_type),pointer :: dtset => null()
118 : !hdr: header, it contains the information about the wavefunctions
119 : type(hdr_type), pointer :: hdr => null()
120 : !MPI_enreg: MPI type for the wavefunctions
121 : type(mpi_type), pointer :: MPI_enreg => null()
122 : !psps: pseudopotential information
123 : type(pseudopotential_type), pointer :: psps => null()
124 : !pawtab: PAW tabulated data
125 : type(pawtab_type),pointer :: pawtab(:) => null()
126 : !ebands: eigenvalues
127 : type(ebands_t), pointer :: ebands => null()
128 : !kset: kpoints information
129 : type(wann_ksetting_t) :: kset
130 : !natom: number of atoms
131 : !nspinor: number of spinors
132 : !nsppol: number of spin polarizations
133 : !mband: number of bands
134 : !mkmem: number of kpoints in this process
135 : !nkpt: number of kpoints
136 :
137 : integer :: natom=0, nspinor=0, nsppol=0, mband=0, &
138 : & mkmem=0, nkpt=0, rank=-999, nprocs=-999, comm=-999
139 : contains
140 : !procedure :: init => abstract_wf_init
141 : procedure :: abstract_init
142 : procedure :: free => abstract_wf_free
143 : procedure :: cg_elem => abstract_wf_cg_elem
144 : procedure :: cg_elem_complex => abstract_wf_cg_elem_complex
145 : procedure :: cprj_elem =>abstract_wf_cprj_elem
146 : procedure :: get_cg_ptr => abstract_wf_get_cprj_ptr
147 : procedure :: get_cprj_ptr => abstract_wf_get_cprj_ptr
148 : procedure :: load_cg => abstract_wf_load_cg
149 : !procedure :: show_info
150 : procedure :: get_kgs=> abstract_wf_get_kgs
151 :
152 : end type abstract_wf
153 :
154 :
155 : !type cg_cprj:
156 : ! This type is used to store the information about the wavefunctions
157 : ! and the MPI distribution of the wavefunctions.
158 : ! It is a child of the abstract_wf type and it provides the interface
159 : ! for the wavefunctions in the CG and CPRJ format.
160 : type, public,extends(abstract_wf) :: cg_cprj
161 : real(dp), pointer :: cg(:, :)=>null()
162 : type(pawcprj_type), pointer :: cprj(:,:)=>null()
163 : integer, pointer :: iwav(:,:,:, :)=>null()
164 : integer, allocatable :: icprj(:, :, :)
165 : contains
166 : procedure :: init => cg_cprj_init !initialize the cg_cprj type
167 : procedure :: free => cg_cprj_free !free the cg_cprj type
168 : procedure :: compute_index_cprj !compute the index of the cprj
169 : procedure :: cg_elem
170 : procedure :: cg_elem_complex
171 : procedure :: cprj_elem
172 : procedure :: get_cg_ptr => cg_cprj_get_cprj_ptr
173 : procedure :: get_cprj_ptr => cg_cprj_get_cprj_ptr
174 : procedure :: write_cg_and_cprj_tmpfile
175 : procedure :: remove_tmpfile
176 : procedure :: load_cg
177 : end type cg_cprj
178 :
179 :
180 : !type wfd_t:
181 : ! wavefunctions in the WFD format.
182 : type, public, extends(abstract_wf) :: wfd_wf
183 : ! The working wfd, ebands, etc
184 : type(wfd_t), pointer :: wfd => null()
185 :
186 : type(wave_t), pointer :: waveprt => null()
187 :
188 : ! if wfd is in IBZ, expand to fullBZ
189 : logical :: expanded = .False.
190 : integer, allocatable :: bz2ibz(:, :)
191 : type(wfd_t), pointer :: wfd_ibz => null()
192 : type(ebands_t), pointer :: ebands_ibz => null()
193 : type(hdr_type), pointer :: hdr_ibz => null()
194 :
195 : logical, allocatable :: bks_mask(:, :, :)
196 : logical, allocatable :: keep_ur(:, :, :)
197 :
198 : type(wfd_t) :: wfd_bz
199 : type(ebands_t) :: ebands_bz
200 : type(hdr_type) :: hdr_bz
201 : type(dataset_type) :: dtset_bz
202 : type(MPI_type) :: mpi_enreg_bz
203 : contains
204 : procedure :: init => wfd_wf_init
205 : procedure :: free => wfd_wf_free
206 : !procedure :: get_ug => wfd_wf_ug
207 : procedure :: cg_elem => wfd_cg_elem
208 : procedure :: cg_elem_complex => wfd_cg_elem_complex
209 : procedure :: cprj_elem =>wfd_cprj_elem
210 : procedure :: load_cg => wfd_load_cg
211 : end type wfd_wf
212 :
213 : public :: init_mywfc, compute_iwav, write_cg_and_cprj
214 :
215 : contains
216 :
217 : !-----------------------------------------------------------------------------
218 : !> Initialize the wann_ksetting_t type
219 : !> @param self: the wann_ksetting_t type
220 : !> @param cryst: the crystal structure
221 : !> @param nkpt: number of kpoints
222 : !> @param mband: number of bands
223 : !> @param nsppol: number of spin polarizations
224 : !> @param kkpts: the kpoints
225 : !> @param comm: MPI communicator
226 : !> @param nprocs: number of processes
227 : !> @param rank: rank of the process
228 : !-----------------------------------------------------------------------------
229 5 : subroutine wann_ksetting_init(self, cryst, nkpt, mband, &
230 5 : & nsppol, kkpts, comm, nprocs, rank)
231 : class(wann_ksetting_t), intent(inout) :: self
232 : type(crystal_t), target, intent(in):: cryst
233 : integer, intent(in) :: nkpt, mband, nsppol, comm, nprocs, rank
234 : real(dp), target, intent(in) :: kkpts(:, :)
235 5 : self%comm=comm
236 5 : self%nprocs=nprocs
237 5 : self%rank=rank
238 5 : self%cryst=>cryst
239 5 : self%nkpt=nkpt
240 5 : self%mband=mband
241 5 : self%nsppol=nsppol
242 5 : self%my_nspin = nsppol
243 5 : self%kkpts=> kkpts
244 5 : end subroutine wann_ksetting_init
245 :
246 :
247 :
248 : !-----------------------------------------------------------------------------
249 : !> compute the ovikp array (the overlap information of the kpoints)
250 : !> @param self: the wann_ksetting_t type
251 : !> @param ovikp: the overlap information of the kpoints
252 : !> @param nntot: number of nearest neighbors
253 : !> @param num_nnmax: maximum number of nearest neighbors
254 : !> @param mpi_enreg: MPI information
255 : !-----------------------------------------------------------------------------
256 5 : subroutine wann_ksetting_set_ovikp(self, ovikp, nntot, num_nnmax, mpi_enreg)
257 : class(wann_ksetting_t), intent(inout) :: self
258 : integer, intent(in) :: nntot, num_nnmax
259 : integer, intent(in) :: ovikp(:, :)
260 : type(mpi_type), intent(inout) :: mpi_enreg
261 5 : if (self%has_ovikp) then
262 0 : ABI_ERROR("ovikp already set!")
263 : end if
264 5 : self%has_ovikp = .True.
265 5 : self%nntot = nntot
266 5 : self%num_nnmax=num_nnmax
267 20 : ABI_MALLOC(self%ovikp, (self%nkpt, num_nnmax))
268 314 : self%ovikp(:,:) = ovikp(:, :)
269 5 : call self%distribute_mpi(mpi_enreg)
270 5 : end subroutine wann_ksetting_set_ovikp
271 :
272 :
273 : !-----------------------------------------------------------------------------
274 : !> free the wann_ksetting_t type
275 : !> @param self: the wann_ksetting_t type
276 : !-----------------------------------------------------------------------------
277 33 : subroutine wann_ksetting_free(self)
278 : class(wann_ksetting_t), intent(inout) :: self
279 33 : nullify(self%cryst)
280 33 : nullify(self%kkpts)
281 33 : if (self%has_ovikp) then
282 5 : ABI_SFREE(self%ovikp)
283 5 : ABI_FREE(self%my_spins)
284 5 : ABI_FREE(self%my_ikpts)
285 5 : ABI_FREE(self%my_kkpts)
286 5 : ABI_FREE(self%my_ikpts_pnn)
287 5 : ABI_FREE(self%my_kkpts_pnn)
288 : end if
289 33 : end subroutine wann_ksetting_free
290 :
291 :
292 : !-----------------------------------------------------------------------------
293 : !> set the mpi distribution using the information in the wann_ksetting_t type
294 : !> @param self: the wann_ksetting_t type
295 : !> @param mpi_enreg: MPI information
296 : !-----------------------------------------------------------------------------
297 5 : subroutine wann_ksetting_distribute_mpi(self, mpi_enreg)
298 : class(wann_ksetting_t), intent(inout) :: self
299 : type(MPI_type), intent(inout) :: mpi_enreg
300 : integer :: ikpt, inn, ik_me, ik_nn, ispin
301 10 : logical :: belongs(self%nkpt)
302 : integer :: counter
303 5 : MPI_enreg%comm_cell=self%comm
304 5 : mpi_enreg%me=self%rank
305 5 : mpi_enreg%me_kpt=self%rank
306 5 : mpi_enreg%nproc = self%nprocs
307 5 : MPI_enreg%paral_spinor=0
308 : !write(std_out,*) "Distributed mpi:", "rank:", self%rank, "comm", self%comm
309 : !write(std_out,*) "nprocs:", self%nprocs, "me:", mpi_enreg%me, "me_kpt:", mpi_enreg%me_kpt
310 5 : if (.not. allocated(mpi_enreg%proc_distrb))then
311 25 : ABI_MALLOC(mpi_enreg%proc_distrb, (self%nkpt, self%mband, self%nsppol) )
312 : end if
313 3386 : mpi_enreg%proc_distrb(:, :, :) =-999
314 :
315 :
316 5 : self%my_nkpt=0
317 5 : self%my_nkpt_pnn=0
318 :
319 15 : ABI_MALLOC(self%my_spins, (self%nsppol))
320 5 : self%my_nspin=self%nsppol
321 10 : do ispin =1, self%nsppol
322 10 : self%my_spins(ispin)=ispin
323 : end do
324 : ! split kpoints
325 5 : call xmpi_split_block(self%nkpt, self%comm, self%my_nkpt, self%my_ikpts)
326 :
327 15 : ABI_MALLOC(self%my_kkpts, (3, self%my_nkpt))
328 :
329 309 : belongs(:)=.False.
330 309 : do ikpt=1, self%my_nkpt
331 304 : ik_me=self%my_ikpts(ikpt)
332 1216 : self%my_kkpts(:, ikpt) = self%kkpts(:, ik_me)
333 304 : belongs(ik_me) = .True.
334 3936 : MPI_enreg%proc_distrb(ik_me,:,:)= self%rank
335 613 : do inn=1, self%nntot
336 304 : ik_nn=self%ovikp(ikpt, inn)
337 608 : belongs(ik_nn) = .True.
338 : end do
339 : end do
340 :
341 5 : self%my_nkpt_pnn=0
342 309 : do ikpt =1, self%nkpt
343 309 : if (belongs(ikpt)) self%my_nkpt_pnn=self%my_nkpt_pnn+1
344 : end do
345 :
346 15 : ABI_MALLOC(self%my_ikpts_pnn, (self%my_nkpt_pnn))
347 15 : ABI_MALLOC(self%my_kkpts_pnn, (3, self%my_nkpt_pnn))
348 :
349 5 : counter=0
350 309 : do ikpt =1, self%nkpt
351 309 : if (belongs(ikpt)) then
352 304 : counter =counter+1
353 304 : self%my_ikpts_pnn(counter) = ikpt
354 1216 : self%my_kkpts_pnn(:, counter) = self%kkpts(:, ikpt)
355 : end if
356 : end do
357 :
358 :
359 5 : end subroutine wann_ksetting_distribute_mpi
360 :
361 : !-----------------------------------------------------------------------------
362 : !> using the kpoint information in the wann_ksetting_t type, set the
363 : !> bks_mask and keep_ur arrays.
364 : !> @param self: the wann_ksetting_t type
365 : !> @param bks_mask: the bks_mask array
366 : !> @param keep_ur: the keep_ur array
367 : !> @param nband: the number of bands
368 : !> @param nsppol: the number of spins
369 : !> @param keep_ur_value: the value to set the keep_ur array to
370 : !-----------------------------------------------------------------------------
371 5 : subroutine wann_ksetting_get_bks_mask(self, bks_mask, keep_ur, nband, nsppol, keep_ur_value)
372 : class(wann_ksetting_t), intent(inout) :: self
373 : integer, intent(in) :: nband, nsppol
374 : logical, allocatable, intent(inout) :: bks_mask(:, :, :) ! bank, kpt, spin
375 : logical, allocatable, intent(inout) :: keep_ur(:, :, :)
376 : logical, intent(in) :: keep_ur_value
377 : integer :: ikpt, ispin, inn, ik_nn, ik_me
378 :
379 25 : ABI_MALLOC(bks_mask, (nband, self%nkpt, nsppol))
380 20 : ABI_MALLOC(keep_ur, (nband, self%nkpt, nsppol))
381 3642 : bks_mask(:, :, :) = .False.
382 3642 : keep_ur(:, :, :) = .False.
383 10 : do ispin =1, nsppol
384 314 : do ikpt =1 , self%my_nkpt
385 304 : ik_me = self%my_ikpts(ikpt)
386 3936 : bks_mask(:, ik_me, :) =.True.
387 3936 : keep_ur(:, ik_me, :) = keep_ur_value
388 613 : do inn =1, self%nntot
389 304 : ik_nn=self%ovikp(ikpt, inn)
390 3632 : bks_mask(:, ik_nn, ispin)=.True.
391 3936 : keep_ur(:, ik_nn, ispin) = keep_ur_value
392 : end do
393 : end do
394 : end do
395 5 : end subroutine wann_ksetting_get_bks_mask
396 :
397 :
398 : !-----------------------------------------------------------------------------
399 : !> using the kpoint information in the wann_ksetting_t type, set the
400 : !> ecut, mpw, and gmax values.
401 : !> @param self: the wann_ksetting_t type
402 : !> @param ecut: the ecut value
403 : !> @param mpw: the mpw value
404 : !> @param gmax: the gmax array
405 : !-----------------------------------------------------------------------------
406 5 : subroutine wann_ksetting_get_mpw_gmax(self, ecut, mpw, gmax)
407 : class(wann_ksetting_t), intent(inout) :: self
408 : real(dp),intent(in) :: ecut
409 : integer,intent(out) :: mpw, gmax(3)
410 : integer,parameter :: istwfk1 = 1
411 : !real(dp) :: cpu, wall, gflops !weight_k,
412 : !type(gqk_t),pointer :: gqk
413 : !arrays
414 : integer :: my_gmax(3), onpw, ipw, ii, my_mpw, ierr
415 5 : integer,allocatable :: gtmp(:,:)
416 : real(dp) :: kk(3)
417 :
418 : integer :: spin, my_ik
419 :
420 : !----------------------------------------------------------------------
421 :
422 5 : mpw = 0; gmax = 0
423 :
424 : ! TODO: This is an hotspot due to the double loop over k and q.
425 : ! Should use a geometrical approach to compute mpw and gmax.
426 :
427 5 : call wrtout(std_out, " Computing mpw. This may take some time for dense k/q meshes...")
428 : !call cwtime(cpu, wall, gflops, "start")
429 :
430 : !do my_is=1,gstore%my_nspins
431 : ! gqk => gstore%gqk(my_is)
432 : ! spin = gstore%my_spins(my_is)
433 10 : do spin =1, self%nsppol
434 :
435 314 : do my_ik=1,self%my_nkpt_pnn
436 1216 : kk = self%my_kkpts_pnn(:, my_ik)
437 :
438 : ! Compute G sphere, returning npw. Note istwfk == 1.
439 304 : call get_kg(kk, istwfk1, ecut, self%cryst%gmet, onpw, gtmp)
440 304 : mpw = max(mpw, onpw)
441 69743 : do ipw=1,onpw
442 278060 : do ii=1,3
443 277756 : gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw)))
444 : end do
445 : end do
446 309 : ABI_FREE(gtmp)
447 : end do ! my_ik
448 : end do ! my_is
449 :
450 5 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, self%comm, ierr)
451 5 : my_gmax = gmax; call xmpi_max(my_gmax, gmax, self%comm, ierr)
452 :
453 5 : call wrtout(std_out, sjoin(' Optimal value of mpw: ', itoa(mpw)))
454 : !call cwtime_report(" gmax and mpw", cpu, wall, gflops)
455 :
456 5 : end subroutine wann_ksetting_get_mpw_gmax
457 :
458 :
459 : !-----------------------------------------------------------------------------
460 : !> Initialize the abstract_wf type, depending on the type of wavefunction
461 : !> If the optional wfd, cg/cprj input are present, then the abstract_wf
462 : !> is initialized as the wfd_wf type or the cg_cprj_wf type.
463 : !> @param mywfc: the abstract_wf type
464 : !> @param ebands: the ebands type
465 : !> @param wfd: the wfd type (optional)
466 : !> @param cg: the cg array (optional)
467 : !> @param cprj: the cprj array (optional)
468 : !> @param cryst: the crystal_t type
469 : !> @param dtset: the dataset_type type
470 : !> @param dtfil: the datafiles_type type
471 : !> @param hdr: the header_type type
472 : !> @param MPI_enreg: the mpi_type type
473 : !> @param nprocs: the number of processors
474 : !> @param psps: the psps_type type
475 : !> @param pawtab: the pawtab_type type
476 : !> @param rank: the rank of the processor
477 : !> @param comm: the communicator
478 : !-----------------------------------------------------------------------------
479 19 : subroutine init_mywfc(mywfc, ebands, wfd , cg, cprj, cryst, &
480 19 : & dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
481 : class(abstract_wf), pointer, intent(inout) :: mywfc
482 : type(crystal_t), target, intent(in) :: cryst
483 : type(ebands_t), target, optional, intent(in) :: ebands
484 : type(wfd_t), target, optional, intent(inout) :: wfd
485 : real(dp), target, optional, intent(in):: cg(:, :)
486 : type(pawcprj_type), target, optional, intent(in):: cprj(:,:)
487 : type(dataset_type),target, intent(in) :: dtset
488 : type(datafiles_type),target, intent(in) :: dtfil
489 : type(mpi_type), target, intent(inout) :: MPI_enreg
490 : type(pseudopotential_type), target, intent(in) :: psps
491 : type(pawtab_type), target, optional, intent(in) :: pawtab(:)
492 : type(hdr_type), target, intent(in) :: hdr
493 : integer, intent(in) :: nprocs, rank, comm
494 :
495 19 : if(present(cg)) then
496 14 : ABI_MALLOC_TYPE_SCALAR(cg_cprj, mywfc)
497 5 : else if (present(wfd)) then
498 815 : ABI_MALLOC_TYPE_SCALAR(wfd_wf, mywfc)
499 : end if
500 : select type(mywfc)
501 : type is(cg_cprj)
502 : call mywfc%init( ebands=ebands, cg=cg, cprj=cprj, cryst=cryst, dtset=dtset, &
503 : & dtfil=dtfil, hdr=hdr, MPI_enreg=mpi_enreg, nprocs=nprocs, &
504 14 : & psps=psps, pawtab=pawtab, rank=rank, comm=comm)
505 : type is(wfd_wf)
506 5 : call mywfc%init( ebands, wfd, cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
507 : !call wfd_print_norm(mywfc%wfd, mywfc%hdr)
508 : end select
509 19 : end subroutine init_mywfc
510 :
511 : !-----------------------------------------------------------------------------
512 : !> Initialize the abstract_wf type, common part for both wfd_wf and cg_cprj
513 : !> @param self: the abstract_wf type
514 : !> @param ebands: the ebands type
515 : !> @param cryst: the crystal_t type
516 : !> @param dtset: the dataset_type type
517 : !> @param dtfil: the datafiles_type type
518 : !> @param hdr: the header_type type
519 : !> @param MPI_enreg: the mpi_type type
520 : !> @param nprocs: the number of processors
521 : !> @param psps: the psps_type type
522 : !> @param pawtab: the pawtab_type type
523 : !> @param rank: the rank of the processor
524 : !> @param comm: the communicator
525 : !-----------------------------------------------------------------------------
526 19 : subroutine abstract_init(self, ebands, cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
527 : class(abstract_wf), intent(inout) :: self
528 : type(crystal_t), target, intent(in) :: cryst
529 : type(ebands_t), target, intent(in) :: ebands
530 : type(dataset_type),target, intent(in) :: dtset
531 : type(datafiles_type),target, intent(in) :: dtfil
532 : type(mpi_type), target, intent(inout) :: MPI_enreg
533 : type(pseudopotential_type), target, intent(in) :: psps
534 : type(pawtab_type), target, optional, intent(in) :: pawtab(:)
535 : type(hdr_type), target, intent(in) :: hdr
536 : integer, intent(in) :: nprocs, rank, comm
537 19 : if(present(pawtab)) self%pawtab => pawtab
538 19 : self%ebands => ebands
539 19 : self%cryst => cryst
540 19 : self%dtset => dtset
541 19 : self%dtfil => dtfil
542 19 : self%hdr => hdr
543 19 : self%MPI_enreg => MPI_enreg
544 19 : self%psps => psps
545 19 : self%natom = hdr%natom
546 19 : self%nspinor = hdr%nspinor
547 19 : self%nsppol = hdr%nsppol
548 19 : self%mband= hdr%mband
549 19 : self%mkmem = dtset%mkmem
550 19 : self%nkpt = hdr%nkpt
551 19 : self%rank = rank
552 19 : self%nprocs = nprocs
553 19 : self%comm=comm
554 19 : end subroutine abstract_init
555 :
556 : !-----------------------------------------------------------------------------
557 : !> get one element of the wavefunction for any abstract_wf type
558 : !> @param self: the abstract_wf type
559 : !> @param icplx: the complex flag
560 : !> @param ig: the index of the G vector
561 : !> @param ispinor: the index of the spinor
562 : !> @param iband: the index of the band
563 : !> @param ikpt: the index of the k point
564 : !> @param isppol: the index of the spin polarization
565 : !> @return: the element of the wavefunction
566 : !-----------------------------------------------------------------------------
567 0 : function abstract_wf_cg_elem(self, icplx, ig, ispinor, iband, ikpt, isppol ) result(res)
568 : class(abstract_wf), intent(inout) :: self
569 : integer, intent(in) :: icplx, ig, ispinor, iband, ikpt, isppol
570 : real(dp) :: res
571 0 : ABI_UNUSED_A(self)
572 : ABI_UNUSED(icplx)
573 : ABI_UNUSED(ig)
574 : ABI_UNUSED(ispinor)
575 : ABI_UNUSED(iband)
576 : ABI_UNUSED(ikpt)
577 : ABI_UNUSED(isppol)
578 : ABI_UNUSED(res)
579 0 : ABI_ERROR("Function should be overrided:")
580 0 : end function abstract_wf_cg_elem
581 :
582 : !-----------------------------------------------------------------------------
583 : !> get one element of the wavefunction for any abstract_wf type, complex version
584 : !-----------------------------------------------------------------------------
585 0 : function abstract_wf_cg_elem_complex(self, ig, ispinor, iband, ikpt, isppol ) result(res)
586 : class(abstract_wf), intent(inout) :: self
587 : integer, intent(in) :: ig, ispinor, iband, ikpt, isppol
588 : complex(dp) :: res
589 0 : ABI_UNUSED_A(self)
590 : ABI_UNUSED(ig)
591 : ABI_UNUSED(ispinor)
592 : ABI_UNUSED(iband)
593 : ABI_UNUSED(ikpt)
594 : ABI_UNUSED(isppol)
595 : ABI_UNUSED(res)
596 0 : ABI_ERROR("Function should be overrided:")
597 0 : end function abstract_wf_cg_elem_complex
598 :
599 : !-----------------------------------------------------------------------------
600 : !> get a block of the wavefunction for any abstract_wf type.
601 : !> It should be overrided for each type.
602 : !> @param self: the abstract_wf type
603 : !> @param ikpt2: the index of the k point
604 : !> @param isppol: the index of the spin polarization
605 : !> @param cg_read: the wavefunction block
606 : !-----------------------------------------------------------------------------
607 0 : subroutine abstract_wf_load_cg(self, ikpt2, isppol, cg_read)
608 : class(abstract_wf), intent(inout) :: self
609 : integer, intent(in) :: ikpt2, isppol
610 : real(dp), intent(inout) :: cg_read(:, :) ! (2, nspinor*mpw*mband )
611 0 : ABI_UNUSED_A(self)
612 : ABI_UNUSED(ikpt2)
613 : ABI_UNUSED(isppol)
614 : ABI_UNUSED(cg_read)
615 0 : ABI_ERROR("This function abstarct_wf_load_cg should be overrided.")
616 0 : end subroutine abstract_wf_load_cg
617 :
618 :
619 :
620 : !-----------------------------------------------------------------------------
621 : !> get one cprj element of the wavefunction for any abstract_wf type, complex version,
622 : !> @param self: the abstract_wf type
623 : !> @param icplx: the complex flag
624 : !> @param ispinor: the index of the spinor
625 : !> @param iband: the index of the band
626 : !> @param ikpt: the index of the k point
627 : !> @param isppol: the index of the spin polarization
628 : !> @param iatom: the index of the atom
629 : !> @param ilmn: the index of the lmn
630 : !> @return: the element of the wavefunction
631 : !-----------------------------------------------------------------------------
632 0 : function abstract_wf_cprj_elem(self,icplx,ispinor, iband, ikpt, isppol, iatom, ilmn) result(res)
633 : class(abstract_wf), intent(inout) :: self
634 : integer, intent(in) :: icplx, ispinor, iband, ikpt, isppol, ilmn, iatom
635 : real(dp) :: res
636 0 : ABI_UNUSED_A(self)
637 : ABI_UNUSED(icplx)
638 : ABI_UNUSED(ispinor)
639 : ABI_UNUSED(iband)
640 : ABI_UNUSED(ikpt)
641 : ABI_UNUSED(isppol)
642 : ABI_UNUSED(iatom)
643 : ABI_UNUSED(ilmn)
644 : ABI_UNUSED(res)
645 0 : ABI_ERROR("Function should be overrided:")
646 0 : end function abstract_wf_cprj_elem
647 :
648 : !-----------------------------------------------------------------------------
649 : !> get a pointer to the cg part of wavefunction for any abstract_wf type.
650 : !> It should be overrided for each type.
651 : !> @param self: the abstract_wf type
652 : !> @return: the pointer to the wavefunction block, in the form of a 2D array.
653 : !-----------------------------------------------------------------------------
654 : function abstract_wf_get_cg_ptr(self) result(cg)
655 : class(abstract_wf), target, intent(inout) :: self
656 : real(dp), pointer :: cg(:,:)
657 : ABI_UNUSED_A(self)
658 : ABI_UNUSED(cg)
659 : ABI_ERROR("The function abstract_wf%get_cg_ptr is not implemented")
660 : end function abstract_wf_get_cg_ptr
661 :
662 : !-----------------------------------------------------------------------------
663 : !> get a pointer to the cprj part of wavefunction for any abstract_wf type.
664 : !> It should be overrided for each type.
665 : !> @param self: the abstract_wf type
666 : !> @return: the pointer to the wavefunction block, in the form of a 2D array.
667 : !-----------------------------------------------------------------------------
668 0 : function abstract_wf_get_cprj_ptr(self) result(cprj)
669 : class(abstract_wf), target, intent(inout) :: self
670 : type(pawcprj_type), pointer :: cprj(:, :)
671 0 : ABI_UNUSED_A(self)
672 0 : ABI_UNUSED_A(cprj)
673 0 : ABI_ERROR("The function abstract_wf%get_cprj_ptr is not implemented")
674 0 : end function abstract_wf_get_cprj_ptr
675 :
676 : !-----------------------------------------------------------------------------
677 : !> get a pointer to the kg of wavefunction for any abstract_wf type.
678 : !> @param self: the abstract_wf type
679 : !> @return: the pointer to the kg, in the form of a 2D array.
680 : !-----------------------------------------------------------------------------
681 5 : subroutine abstract_wf_get_kgs(self, ptr_kg)
682 : class(abstract_wf), intent(inout) :: self
683 : integer, intent(inout) :: ptr_kg(:, :)
684 : integer :: npw_k, ik, ikg, ik_me
685 : real(dp) :: ecut_eff
686 5 : integer, allocatable :: kg_k(:,:)
687 : integer, parameter :: istwfk_1=1
688 288165 : ptr_kg(:,:)=zero
689 : !ecut_eff = dtset%ecut * dtset%dilatmx **2
690 5 : ecut_eff=self%hdr%ecut_eff
691 5 : ikg=0
692 309 : do ik=1, self%kset%my_nkpt
693 304 : ik_me = self%kset%my_ikpts(ik)
694 304 : npw_k = self%hdr%npwarr(ik_me)
695 : call get_kg(self%ebands%kptns(:,ik_me),istwfk_1,ecut_eff, &
696 304 : & self%cryst%gmet,npw_k,kg_k)
697 278060 : ptr_kg(:,1+ikg:npw_k+ikg)=kg_k(:, :) !wfd%Kdata(ik)%kg_k(:,:)
698 304 : ikg =ikg+npw_k
699 309 : ABI_FREE(kg_k)
700 : end do
701 5 : end subroutine abstract_wf_get_kgs
702 :
703 :
704 : !-----------------------------------------------------------------------------
705 : !> free the memory of the abstract_wf type
706 : !-----------------------------------------------------------------------------
707 33 : subroutine abstract_wf_free(self)
708 : class(abstract_wf), intent(inout) :: self
709 33 : call self%kset%free()
710 33 : nullify(self%cryst)
711 33 : nullify(self%dtset)
712 33 : nullify(self%dtfil)
713 33 : nullify(self%hdr)
714 33 : nullify(self%MPI_enreg)
715 33 : nullify(self%psps)
716 33 : nullify(self%pawtab)
717 33 : end subroutine abstract_wf_free
718 :
719 : ! subroutine show_info(self)
720 : ! class(abstract_wf), intent(inout) :: self
721 : ! ! cg_elem
722 : ! integer :: ipw=3, ispinor=1, ikpt=2, iband=4, isppol=1
723 : ! print *, "========showing wf info=============="
724 : ! print *, "ipw:", ipw, " ispinor:", ispinor, " ikpt:", ikpt, " iband:", iband, " isppol:", isppol
725 : ! print *, "cg_elem:", self%cg_elem_complex(ipw, ispinor, ikpt, iband, isppol)
726 : ! print *, "====end showing wf info=============="
727 : ! end subroutine show_info
728 :
729 : !-----------------------------------------------------------------------------
730 : !> free the memory of the wfd_wf type
731 : !-----------------------------------------------------------------------------
732 5 : subroutine wfd_wf_free(self)
733 : class(wfd_wf), intent(inout) :: self
734 : ! TODO reenable this
735 5 : if (self%expanded) then
736 5 : call self%wfd_bz%free()
737 5 : call self%hdr_bz%free()
738 5 : call self%ebands_bz%free()
739 5 : call self%dtset_bz%free()
740 : !call self%mpi_enreg_bz%free()
741 5 : ABI_FREE(self%bz2ibz)
742 5 : ABI_FREE(self%bks_mask)
743 5 : ABI_FREE(self%keep_ur)
744 : end if
745 5 : call self%abstract_wf%free()
746 :
747 5 : end subroutine wfd_wf_free
748 :
749 :
750 : !-----------------------------------------------------------------------------
751 : !> Initialize the wfd_wf type
752 : !-----------------------------------------------------------------------------
753 5 : subroutine wfd_wf_init(self, ebands, wfd,cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
754 : class(wfd_wf), target, intent(inout) :: self
755 : type(ebands_t), target, intent(in) :: ebands
756 : type(crystal_t), target, intent(in) :: cryst
757 : type(wfd_t), target, intent(inout) :: wfd
758 : type(dataset_type),target, intent(in) :: dtset
759 : type(datafiles_type),target, intent(in) :: dtfil
760 : type(mpi_type), target, intent(inout) :: MPI_enreg
761 : type(pseudopotential_type), target, intent(in) :: psps
762 : type(pawtab_type), target, optional, intent(in) :: pawtab(:)
763 : type(hdr_type), target, intent(in) :: hdr
764 : integer, intent(in) :: nprocs, rank, comm
765 :
766 :
767 :
768 5 : self%comm=comm
769 5 : self%rank=rank
770 5 : self%nprocs=nprocs
771 : !print *, "set mpi info to wfd_wf", "self%comm", "self%rank", "self%nprocs"
772 :
773 5 : self%expanded=(dtset%kptopt==1 .or. dtset%kptopt==2)
774 5 : if (self%expanded) then
775 : ! expand the kpoints to the full BZ
776 5 : self%expanded=.True.
777 5 : self%ebands_ibz => ebands
778 5 : self%wfd_ibz => wfd
779 5 : self%hdr_ibz => hdr
780 5 : call ebands_and_hdr_expandk()
781 : call self%kset%init(cryst=cryst,nkpt=self%hdr_bz%nkpt, mband=self%hdr_bz%mband, &
782 : & nsppol=self%hdr_bz%nsppol,kkpts= self%hdr_bz%kptns, &
783 5 : & comm=self%comm, nprocs=self%nprocs, rank=self%rank)
784 :
785 5 : self%hdr=> self%hdr_bz
786 5 : call set_fake_ovikp()
787 :
788 5 : call wfd_expandk()
789 5 : self%wfd=> self%wfd_bz
790 5 : call dtset_expandk()
791 : call self%abstract_wf%abstract_init(self%ebands_bz, cryst, self%dtset_bz, dtfil, &
792 5 : & self%hdr_bz, MPI_enreg, nprocs, psps, pawtab, rank, comm)
793 : else ! the kpoints are already in the full BZ and does not need to be expanded
794 0 : self%expanded=.False.
795 0 : self%wfd => wfd
796 0 : self%ebands => ebands
797 0 : self%hdr=> hdr
798 : call self%kset%init(cryst=cryst,nkpt=hdr%nkpt, mband=hdr%mband, &
799 : & nsppol=hdr%nsppol,kkpts= hdr%kptns, &
800 0 : & comm=self%comm, nprocs=self%nprocs, rank=self%rank)
801 0 : call set_fake_ovikp()
802 0 : call self%abstract_wf%abstract_init(ebands,cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
803 : end if
804 : contains
805 :
806 : !-----------------------------------------------------------------------------
807 : !> subsubroutine of wfd_wf_init
808 : !> set the ovikp of the kset. fake means that all the kpoints are assumed to be neighbors
809 : !-----------------------------------------------------------------------------
810 5 : subroutine set_fake_ovikp()
811 10 : integer :: ovikp(self%hdr%nkpt, 1), i
812 309 : do i=1, self%hdr%nkpt
813 309 : ovikp(i, 1)=i
814 : end do
815 5 : call self%kset%set_ovikp(ovikp, 1 , 1, mpi_enreg)
816 5 : end subroutine set_fake_ovikp
817 :
818 :
819 :
820 : !-----------------------------------------------------------------------------
821 : !> subsubroutine of wfd_wf_init
822 : !> expand the kpoints to the full BZ
823 : !-----------------------------------------------------------------------------
824 5 : subroutine dtset_expandk()
825 : integer :: nkpt
826 5 : self%dtset_bz = dtset%copy()
827 5 : ABI_FREE(self%dtset_bz%kpt)
828 5 : ABI_FREE(self%dtset_bz%kptns)
829 5 : ABI_FREE(self%dtset_bz%istwfk)
830 5 : ABI_FREE(self%dtset_bz%nband)
831 5 : ABI_FREE(self%dtset_bz%wtk)
832 5 : nkpt=self%ebands_bz%nkpt
833 5 : self%dtset_bz%kptopt=3
834 5 : self%dtset_bz%nkpt = nkpt
835 15 : ABI_MALLOC(self%dtset_bz%kpt, (3,nkpt))
836 10 : ABI_MALLOC(self%dtset_bz%kptns, (3,nkpt ))
837 15 : ABI_MALLOC(self%dtset_bz%istwfk, (nkpt ))
838 10 : ABI_MALLOC(self%dtset_bz%nband, (nkpt ))
839 15 : ABI_MALLOC(self%dtset_bz%wtk, (nkpt))
840 1221 : self%dtset_bz%kpt(:,:) = self%ebands_bz%kptns(:,:)
841 1221 : self%dtset_bz%kptns(:,:) = self%ebands_bz%kptns(:,:)
842 309 : self%dtset_bz%istwfk(:) = 1
843 309 : self%dtset_bz%nband(:) = self%hdr_bz%nband(:)
844 309 : self%dtset_bz%wtk(:) = self%hdr_bz%wtk(:)
845 5 : self%dtset_bz%mkmem = self%kset%my_nkpt
846 5 : end subroutine dtset_expandk
847 :
848 : !-----------------------------------------------------------------------------
849 : !> subsubroutine of wfd_wf_init
850 : !> expand the ebands and hdr to the full BZ
851 : !-----------------------------------------------------------------------------
852 5 : subroutine ebands_and_hdr_expandk()
853 : real(dp) :: ecut_eff, dksqmax
854 5 : type(wvl_internal_type) :: dummy_wvl
855 : integer:: kptopt3=3
856 : character(len=200) :: msg
857 : ! NOTE: is this OK to assume so?
858 5 : ecut_eff = dtset%ecut * dtset%dilatmx **2
859 : call ebands%expandk(cryst=cryst, ecut_eff=ecut_eff, &
860 : & force_istwfk1=.True., dksqmax=dksqmax, &
861 5 : & bz2ibz=self%bz2ibz, outb=self%ebands_bz)
862 : ! Note: test if force_istwfk1 is not set to True, force rotate
863 5 : if (dksqmax > tol12) then
864 : write(msg, '(3a,es16.6,4a)' )&
865 0 : 'At least one of the k points could not be generated from a symmetrical one.',ch10,&
866 0 : 'dksqmax=',dksqmax,ch10,&
867 0 : 'Action: check your WFK file and k-point input variables',ch10,&
868 0 : ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
869 0 : ABI_ERROR(msg)
870 : end if
871 :
872 : call self%hdr_bz%init_lowlvl(self%ebands_bz,psps,pawtab,dummy_wvl,abinit_version,&
873 : hdr%pertcase,hdr%natom,hdr%nsym,hdr%nspden,hdr%ecut,dtset%pawecutdg,hdr%ecutsm,dtset%dilatmx,&
874 : hdr%intxc,hdr%ixc,hdr%stmbias,hdr%usewvl,dtset%pawcpxocc,dtset%pawspnorb,dtset%ngfft,dtset%ngfftdg,hdr%so_psp,&
875 : hdr%qptn,cryst%rprimd,cryst%xred,hdr%symrel,hdr%tnons,hdr%symafm,hdr%typat,hdr%amu,hdr%icoulomb,&
876 : kptopt3,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,dtset%cellcharge(1),&
877 : dtset%kptrlatt_orig,dtset%kptrlatt,&
878 5 : dtset%nshiftk_orig,dtset%nshiftk,dtset%shiftk_orig,dtset%shiftk)
879 : ! End CP modified
880 :
881 5 : if (psps%usepaw == 1) call pawrhoij_copy(hdr%pawrhoij, self%hdr_bz%pawrhoij)
882 5 : end subroutine ebands_and_hdr_expandk
883 :
884 : !-----------------------------------------------------------------------------
885 : !> subsubroutine of wfd_wf_init
886 : !> whether the kpoint is irreducible
887 : !@param[in] ik the index of the kpoint in the full BZ
888 : !@return true if the kpoint is irreducible
889 : !-----------------------------------------------------------------------------
890 : logical function isirr(ik)
891 : integer, intent(in) :: ik
892 : integer :: isym, itimrev, g0(3)
893 : !ik_ibz = bz2ibz(ikf,1)
894 : isym = self%bz2ibz(ik,2)
895 : itimrev = self%bz2ibz(ik,6)
896 : g0 = self%bz2ibz(ik,3:5) ! IS(k_ibz) + g0 = k_bz
897 : isirr = (isym == 1 .and. itimrev == 0 .and. all(g0 == 0))
898 : end function isirr
899 :
900 :
901 : !-----------------------------------------------------------------------------
902 : !> subsubroutine of wfd_wf_init
903 : !> expand the wavefunctions to the full BZ
904 : !-----------------------------------------------------------------------------
905 5 : subroutine wfd_expandk()
906 5 : integer, allocatable :: istwfk(:)
907 : integer :: ik, spin, band
908 5 : complex(gwp), allocatable :: ug(:)
909 : integer :: work_ngfft(18),gmax(3),indkk(6,1)
910 5 : real(dp),allocatable :: work(:,:,:,:), cg_kbz(:, :, :)
911 : integer ::mpw, mband, npw_kbz, npw_check, size, ik_ibz
912 5 : integer,allocatable :: kg_kbz(:,:), kg_check(:,:)
913 : real(dp):: kk_bz(3), kk_ibz(3)
914 : character(len=500) :: msg
915 :
916 5 : mband= dtset%mband
917 15 : ABI_MALLOC(istwfk, (self%ebands_bz%nkpt))
918 309 : istwfk(:) = 1
919 :
920 5 : call self%kset%get_mpw_gmax(dtset%ecut, mpw, gmax)
921 :
922 : !mpw = maxval(self%hdr_bz%npwarr)
923 5 : mband = self%hdr_bz%mband
924 : !call gstore%get_mpw_gmax(ecut, mpw, gmax)
925 20 : gmax = gmax + 4 ! FIXME: this is to account for umklapp
926 20 : gmax = 2*gmax + 1
927 5 : call ngfft_seq(work_ngfft, gmax)
928 :
929 :
930 25 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
931 15 : ABI_MALLOC(kg_kbz, (3, mpw))
932 20 : ABI_MALLOC(cg_kbz, (2, mpw*self%hdr_bz%nspinor, self%hdr_bz%mband))
933 :
934 :
935 : call self%kset%get_bks_mask( bks_mask=self%bks_mask, keep_ur=self%keep_ur, &
936 5 : & nband=self%hdr_bz%mband, nsppol=hdr%nsppol, keep_ur_value=.False.)
937 3642 : self%bks_mask(:,:,:)=.True.
938 : call self%wfd_bz%init(Cryst=cryst,Pawtab=pawtab,Psps=psps, &
939 : & keep_ur=self%keep_ur,mband=self%hdr_bz%mband,nband=self%ebands_bz%nband, &
940 : &nkibz=self%ebands_bz%nkpt,nsppol=dtset%nsppol,bks_mask=self%bks_mask,&
941 : &nspden=dtset%nspden,nspinor=hdr%nspinor,ecut=dtset%ecut, &
942 : &ecutsm=dtset%ecutsm,dilatmx=dtset%dilatmx, &
943 : &istwfk=istwfk,kibz=self%ebands_bz%kptns,ngfft=wfd%ngfft, &
944 : &nloalg=wfd%nloalg,prtvol=dtset%prtvol,pawprtvol=dtset%pawprtvol,comm=comm,&
945 5 : & use_fnl_dir0der0 = .False.) ! optional
946 :
947 10 : do spin =1, dtset%nsppol
948 314 : do ik=1, self%ebands_bz%nkpt
949 4566752 : work(:, :, :, :) =0.0_dp
950 288464 : kg_kbz(:, :)=0
951 2299856 : cg_kbz(:, :, :)=0.0_dp
952 304 : npw_kbz=self%hdr_bz%npwarr(ik)
953 304 : size =self%hdr_bz%nspinor*npw_kbz
954 : !print *, "size", npw_kbz
955 : !print *, "size", size
956 2128 : indkk(:, 1) = self%bz2ibz(ik,: )
957 : !TODO: ecut or ecut_eff?
958 : !print *, "indkk: ", indkk(:, 1)
959 :
960 304 : ik_ibz =self%bz2ibz(ik, 1)
961 1216 : kk_bz=self%ebands_bz%kptns(:, ik)
962 1216 : kk_ibz=self%ebands_ibz%kptns(:,ik_ibz )
963 :
964 :
965 : !Note that we use force_rotate here.
966 : !Otherwise the sym_ug_kg gives different npw_kbz as in ebands_bz or hdr_bz.
967 : !if the kpoint is in the IBZ.
968 : call wfd%sym_ug_kg(ecut=dtset%ecut, &
969 : & kk_bz=kk_bz, kk_ibz=kk_ibz, bstart=1, nband=mband, &
970 : & spin=spin, mpw=mpw, indkk=indkk, cryst=cryst, &
971 : & work_ngfft=work_ngfft, work=work, istwf_kbz=istwfk(ik), &
972 : & npw_kbz=npw_kbz, kg_kbz=kg_kbz, cgs_kbz=cg_kbz, &
973 304 : & force_rotate=.True.)
974 :
975 : ! The rotated coefficients are stored in the order defined by kg_kbz.
976 : ! get_kgs regenerates this list later for mlwfovlp, so both the set and
977 : ! its ordering must agree exactly.
978 304 : call get_kg(kk_bz, 1, self%hdr_bz%ecut_eff, cryst%gmet, npw_check, kg_check)
979 304 : write(msg, '(a,i0,a,i0,a,i0)') "Inconsistent plane-wave count at full-BZ k-point ", ik, &
980 608 : ": sym_ug_kg gives ", npw_kbz, ", get_kg gives ", npw_check
981 304 : ABI_CHECK(npw_check == npw_kbz, msg)
982 304 : write(msg, '(a,i0)') "Inconsistent G-vector ordering at full-BZ k-point ", ik
983 278060 : ABI_CHECK(all(kg_check(:,1:npw_check) == kg_kbz(:,1:npw_kbz)), msg)
984 304 : ABI_FREE(kg_check)
985 :
986 304 : self%hdr_bz%npwarr(ik)=npw_kbz
987 304 : self%ebands_bz%npwarr(ik)=npw_kbz
988 304 : self%wfd_bz%npwarr(ik)=npw_kbz
989 304 : size =self%hdr_bz%nspinor*npw_kbz
990 912 : ABI_MALLOC(ug, (size))
991 3632 : do band = 1, self%ebands_bz%nband(ik)
992 739804 : ug(:)= cmplx(cg_kbz(1, 1:size,band), cg_kbz(2, 1:size, band), kind=gwp)
993 : !ug(:) = ug(:) / sqrt(sum(cg_kbz(:, 1:size, band)**2))
994 : call self%wfd_bz%push_ug(band, ik, spin, Cryst,ug, &
995 3632 : & update_ur=.True., update_cprj=.False.)
996 : end do
997 309 : ABI_FREE(ug)
998 : end do
999 : end do
1000 :
1001 5 : ABI_FREE(cg_kbz)
1002 5 : ABI_FREE(kg_kbz)
1003 5 : ABI_FREE(work)
1004 5 : ABI_FREE(istwfk)
1005 5 : end subroutine wfd_expandk
1006 :
1007 :
1008 : end subroutine wfd_wf_init
1009 :
1010 :
1011 :
1012 : !
1013 : ! subroutine wfd_build_cache(self, iband, ik, isppol)
1014 : ! class(wfd_wf), intent(inout) :: self
1015 : ! integer, intent(in) :: iband, ik, isppol
1016 : ! integer :: ik_ibz, size
1017 : ! ik_ibz=ik
1018 : ! size = self%wfd%npwarr(ik_ibz) * self%wfd%nspinor
1019 : ! if(allocated(self%cg_cache)) then
1020 : ! ABI_FREE(self%cg_cache)
1021 : ! end if
1022 : ! ABI_MALLOC(self%cg_cache, (2, size))
1023 : ! call self%wfd%copy_cg(iband, ik_ibz, isppol, self%cg_cache)
1024 : ! self%iband_c = iband
1025 : ! self%ikpt_c = ik
1026 : ! self%isppol_c = isppol
1027 : ! end subroutine wfd_build_cache
1028 :
1029 :
1030 : ! subroutine wfd_wf_ug(self, iband, ikpt, isppol, ug )
1031 : ! class(wfd_Wf), intent(inout) :: self
1032 : ! integer, intent(in) :: iband, ikpt, isppol
1033 : ! real(dp), intent(inout) :: ug(:, :)
1034 : ! integer :: ik_ibz
1035 : ! integer :: npw_k
1036 : ! type(wave_t),pointer :: wave
1037 : ! character(len=500) :: msg
1038 : ! ik_ibz=ikpt
1039 : ! !if(.not. (iband==self%iband_c .and. ikpt==self%ikpt_c .and. isppol==self%isppol_c)) then
1040 : ! !print *, "Building cache for : ", iband, ikpt, isppol
1041 : ! !call self%build_cache(iband, ikpt, isppol)
1042 : ! !end if
1043 : ! ABI_CHECK(self%wfd%get_wave_ptr(iband, ik_ibz, isppol, wave, msg) == 0, msg)
1044 : ! if (.not. wave%has_ug == WFD_STORED) then
1045 : ! write(msg,'(a,i0,a,3i0)')" Node ",self%wfd%my_rank," doesn't have (band,ik_ibz,spin): ",iband,ik_ibz,isppol
1046 : ! ABI_BUG(msg)
1047 : ! end if
1048 : ! npw_k = self%Wfd%npwarr(ik_ibz)
1049 : ! call xcopy(npw_k*self%Wfd%nspinor, wave%ug, 1, ug, 1)
1050 : ! end subroutine wfd_wf_ug
1051 :
1052 :
1053 : !---------------------------------------------------------------------
1054 : !> get one element of the wavefunction in the wfd_Wf object
1055 : !> override the abstract_wf method
1056 : !---------------------------------------------------------------------
1057 459503616 : function wfd_cg_elem(self, icplx, ig, ispinor, iband, ikpt, isppol ) result(res)
1058 : class(wfd_Wf), intent(inout) :: self
1059 : integer, intent(in) :: icplx, ig, ispinor, iband, ikpt, isppol
1060 : integer :: ik_ibz
1061 : real(dp) :: res
1062 : complex(dp) :: t
1063 : integer :: npw_k
1064 : type(wave_t),pointer :: wave
1065 : character(len=500) :: msg
1066 459503616 : ik_ibz=ikpt
1067 : !if(.not. (iband==self%iband_c .and. ikpt==self%ikpt_c .and. isppol==self%isppol_c)) then
1068 : !print *, "Building cache for : ", iband, ikpt, isppol
1069 : !call self%build_cache(iband, ikpt, isppol)
1070 : !end if
1071 459503616 : ABI_CHECK(self%wfd%get_wave_ptr(iband, ik_ibz, isppol, wave, msg) == 0, msg)
1072 459503616 : if (.not. wave%has_ug == WFD_STORED) then
1073 0 : write(msg,'(a,i0,a,3i0)')" Node ",self%wfd%my_rank," doesn't have (band,ik_ibz,spin): ",iband,ik_ibz,isppol
1074 0 : ABI_BUG(msg)
1075 : end if
1076 459503616 : npw_k = self%Wfd%npwarr(ik_ibz)
1077 : !call xcopy(npw_k*Wfd%nspinor, wave%ug, 1, ug, 1)
1078 : !res = self%cg_cache(icplx, ig+self%hdr%npwarr(ik_ibz)*(ispinor-1))
1079 459503616 : t = wave%ug(ig+self%hdr%npwarr(ik_ibz)*(ispinor-1))
1080 689255424 : select case(icplx)
1081 : case(1)
1082 229751808 : res = real(t)
1083 : case(2)
1084 229751808 : res = aimag(t)
1085 : case default
1086 459503616 : res=-999999.99_dp
1087 : end select
1088 459503616 : end function wfd_cg_elem
1089 :
1090 : !---------------------------------------------------------------------
1091 : !> get one element of the wavefunction in the wfd_Wf object, complex version
1092 : !> override the abstract_wf method
1093 : !---------------------------------------------------------------------
1094 1756336 : function wfd_cg_elem_complex(self, ig, ispinor, iband, ikpt, isppol ) result(res)
1095 : class(wfd_wf), intent(inout) :: self
1096 : integer, intent(in) :: ig, ispinor, iband, ikpt, isppol
1097 : complex(dp) :: res
1098 : integer :: ik_ibz
1099 : integer :: npw_k
1100 : type(wave_t),pointer :: wave
1101 : character(len=500) :: msg
1102 1756336 : ik_ibz=ikpt
1103 1756336 : ABI_CHECK(self%wfd%get_wave_ptr(iband, ik_ibz, isppol, wave, msg) == 0, msg)
1104 1756336 : if (.not. wave%has_ug == WFD_STORED) then
1105 0 : write(msg,'(a,i0,a,3i0)')" Node ",self%wfd%my_rank," doesn't have (band,ik_ibz,spin): ",iband,ik_ibz,isppol
1106 0 : ABI_BUG(msg)
1107 : end if
1108 1756336 : npw_k = self%Wfd%npwarr(ik_ibz)
1109 1756336 : res = wave%ug(ig+self%hdr%npwarr(ik_ibz)*(ispinor-1))
1110 1756336 : end function wfd_cg_elem_complex
1111 :
1112 : !---------------------------------------------------------------------
1113 : !> get a block of the wavefunction in the wfd_Wf object
1114 : !> override the abstract_wf method
1115 : !---------------------------------------------------------------------
1116 0 : subroutine wfd_load_cg(self, ikpt2, isppol, cg_read)
1117 : class(wfd_wf), intent(inout) :: self
1118 : integer, intent(in) :: ikpt2, isppol
1119 : real(dp), intent(inout) :: cg_read(:, :)
1120 : integer :: iband, iblk, size
1121 0 : iblk=0
1122 0 : size=self%hdr%npwarr(ikpt2) * self%nspinor
1123 0 : do iband =1, self%mband
1124 0 : call self%wfd%copy_cg(iband, ikpt2, isppol, cg_read(:, iblk+1:iblk+size))
1125 0 : iblk = iblk + size
1126 : end do
1127 0 : end subroutine wfd_load_cg
1128 :
1129 :
1130 :
1131 : !---------------------------------------------------------------------
1132 : !> get one cprj element of the wavefunction in the wfd_Wf object
1133 : !> override the abstract_wf method
1134 : !---------------------------------------------------------------------
1135 0 : function wfd_cprj_elem(self,icplx,ispinor, iband, ikpt, isppol, iatom, ilmn) result(res)
1136 : class(wfd_wf), intent(inout) :: self
1137 : integer, intent(in) :: icplx, ispinor, iband, ikpt, isppol, ilmn, iatom
1138 : real(dp) :: res
1139 0 : type(pawcprj_type) :: cprj_out(self%natom,self%nspinor)
1140 :
1141 : integer :: ik_ibz
1142 : !TODO:
1143 0 : ik_ibz = ikpt
1144 : !call self%wfd%ug2cprj(band=iband,ik_ibz=ik_ibz,spin=ispin,choice=1,idir=0,natom=self%natom,Cryst=self%Cryst ,cwaveprj,sorted=.False.)
1145 : call self%wfd%get_cprj( band=iband, ik_ibz=ik_ibz, spin=isppol, &
1146 0 : & Cryst=self%cryst, Cprj_out=cprj_out, sorted=.False.)
1147 0 : res=cprj_out(iatom, ispinor)%cp(icplx, ilmn)
1148 0 : end function wfd_cprj_elem
1149 :
1150 :
1151 : !---------------------------------------------------------------------
1152 : !> get one cpaw element of the wavefunction in the wfd_Wf object, complex version
1153 : !> override the abstract_wf method
1154 : !---------------------------------------------------------------------
1155 14 : subroutine cg_cprj_init(self, ebands, cg, cprj, cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
1156 : class(cg_cprj), intent(inout) :: self
1157 : type(crystal_t), target, intent(in) :: cryst
1158 : type(ebands_t), target, intent(in) :: ebands
1159 : real(dp), target, optional, intent(in):: cg(:, :)
1160 : type(pawcprj_type), target, optional, intent(in):: cprj(:,:)
1161 : type(dataset_type),target, intent(in) :: dtset
1162 : type(datafiles_type),target, intent(in) :: dtfil
1163 : type(mpi_type), target, intent(inout) :: MPI_enreg
1164 : type(pseudopotential_type), target, intent(in) :: psps
1165 : type(pawtab_type), target, optional, intent(in) :: pawtab(:)
1166 : type(hdr_type), target, intent(in) :: hdr
1167 : integer, intent(in) :: nprocs, rank, comm
1168 14 : call self%abstract_wf%abstract_init(ebands, cryst, dtset, dtfil, hdr, MPI_enreg, nprocs, psps, pawtab, rank, comm)
1169 14 : if(present(cg)) self%cg=> cg
1170 14 : if(present(cprj) .and. present(pawtab)) then
1171 14 : self%has_paw=.True.
1172 : end if
1173 14 : if(present(cprj)) then
1174 : !print *, "cprj linked"
1175 : !print *, "before linked", cprj(1,1)%cp(:,:)
1176 14 : self%cprj => cprj
1177 : !print *, "after linked", self%cprj(1,1)%cp(:,:)
1178 : end if
1179 84 : ABI_MALLOC(self%iwav,(self%nspinor, self%mband,self%nkpt,self%nsppol))
1180 14 : call compute_iwav(MPI_enreg, dtset, hdr, self%iwav, nprocs, rank)
1181 14 : call self%compute_index_cprj()
1182 14 : if(nprocs>1) then
1183 0 : call self%write_cg_and_cprj_tmpfile()
1184 : end if
1185 : !call self%show_info()
1186 14 : end subroutine cg_cprj_init
1187 :
1188 : !---------------------------------------------------------------------
1189 : !> free the cg_cprj object ant the temporary files
1190 : !---------------------------------------------------------------------
1191 14 : subroutine cg_cprj_free(self)
1192 : class(cg_cprj), intent(inout) :: self
1193 14 : if(self%nprocs>1) then
1194 0 : call self%remove_tmpfile()
1195 : end if
1196 14 : nullify(self%cg)
1197 14 : nullify(self%cprj)
1198 14 : call self%abstract_wf%free()
1199 14 : ABI_FREE(self%iwav)
1200 14 : nullify(self%iwav)
1201 14 : ABI_FREE(self%icprj)
1202 14 : call self%abstract_wf%free()
1203 14 : end subroutine cg_cprj_free
1204 :
1205 : ! return one entry of cg.
1206 : ! parameters:
1207 : ! icplx: 1 for real part, 2 for imaginary part
1208 : ! ig: index of G vector
1209 : ! ispinor: index of spinor
1210 : ! iband: index of band
1211 : ! ikpt: index of k point
1212 : ! isppol: index of spin
1213 1570800548 : function cg_elem(self, icplx, ig, ispinor, iband, ikpt, isppol ) result(res)
1214 : class(cg_cprj), intent(inout) :: self
1215 : integer, intent(in) :: icplx, ig, ispinor, iband, ikpt, isppol
1216 : integer :: ind
1217 : real(dp) :: res
1218 1570800548 : ind=ig+self%iwav(ispinor, iband,ikpt,isppol)
1219 1570800548 : res=self%cg(icplx,ind)
1220 1570800548 : end function cg_elem
1221 :
1222 :
1223 :
1224 : ! return one entry of cg in complex form.
1225 : ! Parameters:
1226 : ! same as cg_elem, except that icplx is not needed.
1227 0 : function cg_elem_complex(self, ig,ispinor, iband, ikpt, isppol) result(res)
1228 : class(cg_cprj), intent(inout) :: self
1229 : integer, intent(in) :: ig, ispinor, iband, ikpt, isppol
1230 : integer :: ind
1231 : complex(dp) :: res
1232 0 : ind=ig+self%iwav(ispinor, iband,ikpt,isppol)
1233 0 : res=CMPLX(self%cg(1,ind), self%cg(2,ind), kind=dp)
1234 0 : end function cg_elem_complex
1235 :
1236 : !return a pointer to the cg array in the cg_cprj object.
1237 : function cg_cprj_get_cg_ptr(self) result(cg)
1238 : class(cg_cprj), target, intent(inout) :: self
1239 : real(dp), pointer :: cg(:,:)
1240 : cg=> self%cg
1241 : end function cg_cprj_get_cg_ptr
1242 :
1243 : !return a pointer to the cprj array in the cg_cprj object.
1244 2944 : function cg_cprj_get_cprj_ptr(self) result(cprj)
1245 : class(cg_cprj), target, intent(inout) :: self
1246 : type(pawcprj_type), pointer :: cprj(:, :)
1247 2944 : cprj=>self%cprj
1248 2944 : end function cg_cprj_get_cprj_ptr
1249 :
1250 : !---------------------------------------------------------------------
1251 : !> compute the index of the cprj array.
1252 : !---------------------------------------------------------------------
1253 14 : subroutine compute_index_cprj(self)
1254 : ! FIXME:hexu: this is modified from the m_mlwfovlp,
1255 : ! but I think it should be carefully checked.
1256 : ! mcprj=nspinor*mband*mkmem*nsppol
1257 : ! 1. nspinor=2 case seems to be wrong.
1258 : class(cg_cprj), intent(inout) :: self
1259 : integer :: ii, ikpt, isppol, iband
1260 70 : ABI_MALLOC(self%icprj, (self%mband,self%nkpt,self%nsppol))
1261 14 : ii=0
1262 29 : do isppol=1,self%nsppol
1263 : ! FIXME: check if it should be mkmem or nkpt.
1264 429 : do ikpt=1,self%nkpt
1265 : ! FIXME: nband has the shape of (nsppol*nkpt).
1266 : !nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1267 5615 : do iband=1,self%dtset%nband(ikpt)
1268 : ! FIXME: hexu: should cycle if the kpt is not in this node??
1269 5200 : ii=ii+1
1270 5600 : self%icprj(iband,ikpt,isppol)=ii
1271 : end do
1272 : end do
1273 : end do
1274 14 : end subroutine compute_index_cprj
1275 :
1276 :
1277 : ! get one element of cprj
1278 0 : function cprj_elem(self,icplx,ispinor, iband, ikpt, isppol, iatom, ilmn) result(res)
1279 : class(cg_cprj), intent(inout) :: self
1280 : integer, intent(in) :: icplx, ispinor, iband, ikpt, isppol, ilmn, iatom
1281 : real(dp) :: res
1282 : integer :: ig
1283 : ! TODO: this seems to be better than compute_index_cprj,
1284 : ! But should it be mband or nband(ikpt)
1285 :
1286 : ! mcprj=nspinor*mband*mkmem*nsppol
1287 : ! this is the original version in m_mlwfovlp
1288 : !ig=iband+(ikpt-1)*self%mband*self%nspinor + &
1289 : ! &(isppol-1)*self%mkmem*self%mband*self%nspinor
1290 :
1291 : ig=ispinor+(iband-1)*self%nspinor+(ikpt-1)*self%mband*self%nspinor + &
1292 0 : &(isppol-1)*self%mkmem*self%mband*self%nspinor
1293 0 : res= self%cprj(iatom, ig)%cp(icplx, ilmn)
1294 0 : end function cprj_elem
1295 :
1296 : !---------------------------------------------------------------------
1297 : !> write cg and cprj to tmpfile
1298 : !---------------------------------------------------------------------
1299 0 : subroutine write_cg_and_cprj_tmpfile(self)
1300 : class(cg_cprj), intent(inout) :: self
1301 : call write_cg_and_cprj(self%dtset, self%cg, self%cprj, self%dtfil, self%iwav, &
1302 : & self%hdr%npwarr, self%mband, self%natom, &
1303 : & self%nsppol, self%nkpt, self%MPI_enreg, &
1304 0 : & self%rank, self%psps, self%pawtab)
1305 0 : end subroutine write_cg_and_cprj_tmpfile
1306 :
1307 : !---------------------------------------------------------------------
1308 : !> remove tmpfile of cg and cprj
1309 : !---------------------------------------------------------------------
1310 0 : subroutine remove_tmpfile(self)
1311 : class(cg_cprj), intent(inout) :: self
1312 : integer :: isppol, ikpt, ierr
1313 : integer :: master=1
1314 : character(len=fnlen) :: wfnname
1315 : character(len=500) :: message
1316 0 : if(self%dtset%prtvol>0) then
1317 0 : write(message, '(3a)' ) ch10,&
1318 0 : & ' mlwfovlp : Removing temporary files with cg and cprj (PAW)',ch10
1319 0 : call wrtout(ab_out,message,'COLL')
1320 0 : call wrtout(std_out, message,'COLL')
1321 : end if
1322 : !
1323 : ! Just master node will remove the files
1324 : !
1325 0 : if(self%rank==master) then
1326 0 : do isppol=1,self%nsppol
1327 0 : do ikpt=1,self%nkpt
1328 0 : write(wfnname,'(a,I5.5,".",I1)') trim(self%dtfil%fnametmp_cg),ikpt,isppol
1329 0 : call delete_file(wfnname,ierr)
1330 0 : if(self%psps%usepaw==1) then
1331 0 : write(wfnname,'(a,I5.5,".",I1)') trim(self%dtfil%fnametmp_cprj),ikpt,isppol
1332 0 : call delete_file(wfnname,ierr)
1333 : end if
1334 : end do !ikpt
1335 : end do !isppol
1336 : end if
1337 0 : end subroutine remove_tmpfile
1338 :
1339 :
1340 : !---------------------------------------------------------------------
1341 : !> read cg from tmpfile
1342 : !@param[in] ikpt2 index of kpt
1343 : !@param[in] isppol index of spin
1344 : !@param[out] cg_read cg read from tmpfile
1345 : !---------------------------------------------------------------------
1346 0 : subroutine load_cg(self, ikpt2, isppol, cg_read)
1347 : class(cg_cprj), intent(inout) :: self
1348 : integer, intent(in) :: ikpt2, isppol
1349 : real(dp), intent(inout) :: cg_read(:, :)
1350 : character(len=fnlen) :: cg_file
1351 : integer :: npw_k2, ios, ii
1352 : character(len=500) :: message
1353 :
1354 : integer :: iunit , iband2, ipw, index
1355 0 : write(cg_file,'(a,I5.5,".",I1)') trim(self%dtfil%fnametmp_cg),ikpt2,isppol
1356 0 : iunit=1000+ikpt2+ikpt2*(isppol-1)
1357 0 : npw_k2=self%hdr%npwarr(ikpt2)
1358 :
1359 0 : open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios)
1360 0 : if(ios /= 0) then
1361 0 : write(message,*) " mlwfovlp_pw: file",trim(cg_file), "not found"
1362 0 : ABI_ERROR(message)
1363 : end if
1364 : !
1365 0 : do iband2=1,self%mband
1366 0 : do ipw=1,npw_k2*self%nspinor
1367 0 : index=ipw+(iband2-1)*npw_k2*self%nspinor
1368 0 : read(iunit) (cg_read(ii,index),ii=1,2)
1369 : ! if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
1370 : ! if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
1371 : end do
1372 : end do
1373 0 : close(iunit)
1374 0 : end subroutine load_cg
1375 :
1376 :
1377 : !---------------------------------------------------------------------
1378 : !> compute the index of the wavefunctions in the cg array
1379 : !---------------------------------------------------------------------
1380 14 : subroutine compute_iwav(MPI_enreg, dtset, hdr, iwav, nprocs, rank)
1381 : type(mpi_type), intent(in) :: MPI_enreg
1382 : type(dataset_type), intent(in) :: dtset
1383 : type(hdr_type), intent(in) :: hdr
1384 : integer, intent(inout) :: iwav(:, :, :, :)
1385 : ! dimension: (nspinor, mband,nkpt,nsppol))
1386 : integer, intent(in) :: nprocs, rank
1387 : character(len=500) :: message
1388 28 : integer :: icg(hdr%nsppol, hdr%nkpt)
1389 : integer :: icgtemp, isppol, ikpt, nband_k, npw_k, iband, ispinor
1390 :
1391 14 : write(message, '(a,a)' ) ch10,&
1392 28 : & ' mlwfovlp : compute shifts for g-points '
1393 14 : call wrtout(std_out, message,'COLL')
1394 : !----------------------------------------------------------------------
1395 : !Compute shifts for g points (icg,iwav)
1396 : !(here mband is not used, because shifts are internal variables of abinit)
1397 : !----------------------------------------------------------------------
1398 : !write(std_out,*) mpw*dtset%nspinor*mband*mkmem*nsppol
1399 : !ABI_MALLOC(icg,(nsppol,nkpt))
1400 750 : icg=0
1401 29 : icgtemp=0
1402 : !ABI_MALLOC(iwav,(dtset%mband,nkpt,nsppol))
1403 12365 : iwav(:,:,:, :)=0
1404 29 : do isppol=1,hdr%nsppol
1405 429 : do ikpt=1,hdr%nkpt
1406 : ! MPI:cycle over k-points not treated by this node
1407 400 : if (nprocs>1 ) then !sometimes we can have just one processor
1408 0 : if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank) /=0) CYCLE
1409 : end if
1410 :
1411 : ! write(std_out,*)'rank',rank,'ikpt',ikpt,'isppol',isppol
1412 400 : nband_k=dtset%nband(ikpt+(isppol-1)*hdr%nkpt)
1413 : ! write(std_out,*) ikpt+(isppol-1)*nkpt,nkpt
1414 400 : npw_k=hdr%npwarr(ikpt)
1415 5615 : do iband=1,nband_k
1416 5200 : if(iband.gt. dtset%mband) then
1417 0 : write(message,'(a,3i0)')" mband",iband,dtset%mband,nband_k
1418 0 : ABI_ERROR(message)
1419 : end if
1420 11936 : do ispinor =1, dtset%nspinor
1421 : iwav(ispinor, iband,ikpt,isppol)= &
1422 : !& (iband-1)*npw_k*dtset%nspinor+icgtemp &
1423 11936 : & (ispinor-1)*npw_k +icgtemp
1424 : end do
1425 : !icgtemp=icgtemp+ npw_k*dtset%nspinor*nband_k
1426 5600 : icgtemp = icgtemp + npw_k
1427 : end do ! iband
1428 : ! icg(isppol,ikpt)=icgtemp
1429 : ! write(std_out,*) "icg", isppol,ikpt,icg(isppol,ikpt)
1430 : end do ! ikpt
1431 : end do ! isppol
1432 : !write(std_out,*) "shift for cg computed"
1433 : !
1434 : !Shifts computed.
1435 14 : end subroutine compute_iwav
1436 :
1437 : !---------------------------------------------------------------------
1438 : !> write the cg and cprj arrays into a list of files
1439 : !> The name of the files are given by dtset%fnametmp_cg and dtset%fnametmp_cprj
1440 : !> The files are written in unformatted format
1441 : !---------------------------------------------------------------------
1442 0 : subroutine write_cg_and_cprj(dtset, cg, cprj, dtfil, iwav, npwarr, mband, natom, &
1443 0 : &nsppol, nkpt, MPI_enreg, rank, psps, pawtab)
1444 :
1445 : type(dataset_type),intent(in) :: dtset
1446 : type(MPI_type),intent(in) :: mpi_enreg
1447 : type(datafiles_type),intent(in) :: dtfil
1448 : integer,intent(in):: iwav(:,:,:,:)
1449 : real(dp),intent(in) :: cg(:, :)
1450 : type(pawcprj_type), intent(in) :: cprj(:, :)
1451 : integer, intent(in) :: rank, nsppol, nkpt, mband, natom
1452 : integer, intent(in) :: npwarr(nkpt)
1453 : type(pseudopotential_type),intent(in) :: psps
1454 : !type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1455 : type(pawtab_type),intent(in) :: pawtab(:)
1456 :
1457 : integer :: ikpt, ikpt2, isppol, iun_plot, npw_k, iband, i, ipw, ispinor, ig
1458 : integer :: iatom, itypat, lmn_size, ilmn
1459 : character(len=fnlen) :: wfnname
1460 : character(len=1000) :: message
1461 :
1462 0 : if(dtset%prtvol>0) then
1463 0 : write(message, '(3a)' ) ch10,&
1464 0 : & ' mlwfovlp : Creating temporary files with cg and cprj (PAW)',ch10
1465 0 : call wrtout(ab_out,message,'COLL')
1466 0 : call wrtout(std_out, message,'COLL')
1467 : end if
1468 : !
1469 0 : do isppol=1,nsppol
1470 0 : do ikpt=1,nkpt
1471 : !
1472 : ! MPI:cycle over k-points not treated by this node
1473 : !
1474 0 : if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank) /=0) CYCLE
1475 :
1476 : ! write(std_out,*)'writing kpt ',ikpt,'isppol',isppol,' by node ', rank
1477 0 : write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cg),ikpt,isppol
1478 0 : iun_plot=1000+ikpt+ikpt*(isppol-1)
1479 :
1480 0 : open (unit=iun_plot, file=wfnname,form='unformatted')
1481 0 : npw_k=npwarr(ikpt)
1482 0 : do iband=1,mband
1483 0 : do ispinor = 1, dtset%nspinor
1484 0 : do ipw=1,npw_k
1485 0 : write(iun_plot) (cg(i,ipw+iwav(ispinor, iband,ikpt,isppol)),i=1,2)
1486 : end do
1487 : end do
1488 : end do
1489 0 : close(iun_plot)
1490 : end do !ikpt
1491 : end do !isppol
1492 : !
1493 : ! In the PAW case we also need to write out cprj into files
1494 : !
1495 0 : if(psps%usepaw==1) then
1496 : !
1497 : ! big loop on atoms, kpts, bands and lmn
1498 : !print *, "nsppol", nsppol
1499 : !print *, "nkpt", nkpt
1500 : !print *, "mkmem", dtset%mkmem
1501 : !
1502 0 : ikpt2=0
1503 0 : do isppol=1,nsppol
1504 0 : do ikpt=1,nkpt
1505 : !
1506 : ! MPI:cycle over k-points not treated by this node
1507 : !
1508 0 : if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-MPI_enreg%me) /=0) CYCLE
1509 :
1510 0 : ikpt2=ikpt2+1 !sums just on the k-points treated by this node
1511 : !
1512 0 : write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cprj),ikpt,isppol
1513 0 : iun_plot=1000+ikpt
1514 0 : open (unit=iun_plot, file=wfnname,form='unformatted')
1515 : !
1516 0 : do iband=1,mband*dtset%nspinor
1517 : !ig=iband+(ikpt2-1)*mband*dtset%nspinor +(isppol-1)*nkpt*mband*dtset%nspinor !index for cprj(:,ig)
1518 : ! cprj: only mkmem k-points are stored in this node.
1519 0 : ig=iband+(ikpt2-1)*mband*dtset%nspinor +(isppol-1)*dtset%mkmem*mband*dtset%nspinor !index for cprj(:,ig)
1520 : !
1521 0 : do iatom=1,natom
1522 0 : itypat=dtset%typat(iatom)
1523 0 : lmn_size=pawtab(itypat)%lmn_size
1524 : !
1525 0 : do ilmn=1,lmn_size
1526 0 : write(iun_plot) (( cprj(iatom,ig)%cp(i,ilmn)),i=1,2)
1527 : end do !ilmn
1528 : end do !iatom
1529 : end do !iband
1530 :
1531 0 : close(iun_plot)
1532 : end do !ikpt
1533 : end do !isppol
1534 : end if !usepaw==1
1535 0 : end subroutine write_cg_and_cprj
1536 :
1537 :
1538 : ! subroutine wfd_print_norm(wfd, hdr)
1539 : ! type(wfd_t), intent(in) :: wfd
1540 : ! type(hdr_type), intent(in) :: hdr
1541 : ! integer :: spin, band, ikpt, size
1542 : ! real(dp), allocatable :: cgtemp(:, :)
1543 : ! do spin=1, wfd%nsppol
1544 : ! do ikpt=1, wfd%nkibz
1545 : ! size=wfd%nspinor * hdr%npwarr(ikpt)
1546 : ! ABI_MALLOC(cgtemp, (2, size))
1547 : ! do band=1, wfd%mband
1548 : ! cgtemp(:, :)=0
1549 : ! !if(isirr(ik)) then
1550 : ! call wfd%copy_cg(band,ikpt, spin, cgtemp)
1551 : ! !print *, "spin:", spin, "band:", band, "ik:", ikpt, "delta:", sum(cgtemp**2)
1552 : ! end do
1553 : ! ABI_FREE(cgtemp)
1554 : ! end do
1555 : ! end do
1556 : ! end subroutine wfd_print_norm
1557 :
1558 :
1559 119 : end module m_abstract_wf
|
