Branch data 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
+ + - + +
+ - + - +
+ + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - +
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
|
