Line data Source code
1 : !!****m* ABINIT/m_gwls_hamiltonian
2 : !! NAME
3 : !! m_gwls_hamiltonian
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 :
23 : module m_gwls_hamiltonian
24 :
25 : use m_gwls_utility
26 : use m_gwls_wf
27 : use m_dtset
28 : use m_bandfft_kpt
29 : use m_cgtools
30 : use defs_basis
31 : use m_abicore
32 : use m_xmpi
33 : use m_pawang
34 : use m_errors
35 : use m_abi_mixing
36 : use m_mpinfo
37 : use m_crystal
38 :
39 : use defs_abitypes, only : MPI_type
40 : use m_io_tools, only : get_unit
41 : use m_hamiltonian, only : gs_hamiltonian_type
42 : use m_pawcprj, only : pawcprj_type
43 : use m_vcoul, only : vcoul_t
44 : use m_gsphere, only : gsphere_t
45 : use m_bz_mesh, only : kmesh_t
46 : use m_fft, only : fftpac, fourwf
47 : use m_getghc, only : getghc
48 : use m_io_kss, only : make_gvec_kss
49 :
50 : implicit none
51 : save
52 : private
53 : !!***
54 :
55 : !Data passed in argument to build_H and build_vxc that we need to copy for the GW calculation.
56 : !build_H
57 : type(dataset_type) :: dtset !Public
58 : type(MPI_type) :: mpi_enreg !Public
59 : integer :: cpopt
60 : real(dp), allocatable :: cg(:,:) !Public
61 : integer :: dimffnl
62 : type(gs_hamiltonian_type) :: gs_hamk
63 : real(dp), allocatable, target :: ffnl(:,:,:,:)
64 : integer, allocatable, target :: kg_k(:,:)
65 : real(dp), allocatable, target :: kinpw(:)
66 : real(dp), allocatable, target :: ph3d(:,:,:)
67 : real(dp), allocatable, target :: vlocal(:,:,:,:)
68 : !build_vxc
69 : real(dp), allocatable :: vxc_dg(:,:)
70 :
71 : !Other variables initialized in build_H
72 : integer :: ndat
73 : integer :: sij_opt
74 : integer :: tim_getghc
75 : integer :: type_calc
76 : integer :: nfft !Public
77 : integer :: nline !Public
78 : integer :: ngfft(18)
79 : integer :: e !Public
80 : integer :: nspinor !Public
81 : integer :: n1, n2, n3
82 : integer :: n4, n5, n6 !Public
83 : integer :: npw_k !Public
84 : integer :: ckpt !Public
85 : integer :: mgfft
86 : integer :: tim_fourwf
87 : integer :: i
88 : integer :: nband !Public
89 : integer :: ispden !Public
90 : integer :: v
91 : integer :: nbandv !Public
92 : integer :: mcg
93 : integer :: ktot !Public
94 : !integer :: tmp2i(2)
95 : integer, parameter :: iovar=6 !137 !Public
96 : integer, allocatable :: gbound(:,:)
97 : integer, allocatable :: istwfk(:) !Public
98 : real(dp) :: eshift
99 : real(dp) :: tolwfr !Public
100 : real(dp) :: ucvol !Public
101 : real(dp) :: weight
102 : real(dp) :: tmpc(2)
103 : real(dp), allocatable :: vxc(:,:,:,:) !Public (Just the transcript of vxc_dg on the "small" (wfk) real grid)
104 : real(dp), allocatable :: psik1(:,:) !Working wf in k space for module subtourines.
105 : real(dp), allocatable :: psik2(:,:) !Working wf in k space for module subtourines.
106 : real(dp), allocatable :: psik3(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
107 : real(dp), allocatable :: psik4(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
108 : real(dp), allocatable :: psikb1(:,:) !Working wf in k space for module subtourines.
109 : real(dp), allocatable :: psikb2(:,:) !Working wf in k space for module subtourines.
110 : real(dp), allocatable :: psikb3(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
111 : real(dp), allocatable :: psikb4(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
112 : real(dp), allocatable :: psig1(:,:) !Working wf in k space for module subtourines.
113 : real(dp), allocatable :: psig2(:,:) !Working wf in k space for module subtourines.
114 : real(dp), allocatable :: psig3(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
115 : real(dp), allocatable :: psig4(:,:) !Working wf in k space for module subtourines. RESERVED FOR Hpsik.
116 : real(dp), allocatable :: psir1(:,:,:,:) !Working wf in real space (sg) for module subtourines.
117 : real(dp), allocatable :: psir2(:,:,:,:) !Working wf in real space (sg) for module subtourines.
118 : real(dp), allocatable :: psir3(:,:,:,:) !Working wf in real space (sg) for module subtourines.
119 : real(dp), allocatable :: psidg(:,:) !Working wf in real space (dg) for module subtourines.
120 :
121 : real(dp), allocatable :: denpot(:,:,:) !Working array for real space product of wf with fourwf.
122 : real(dp), allocatable :: pcon(:)
123 : real(dp), allocatable :: eig(:) !Public
124 : real(dp), allocatable :: dummy2(:,:), dummy3(:,:,:)
125 : real(dp), allocatable :: scprod2(:,:)
126 : type(pawcprj_type), allocatable :: conjgrprj(:,:) !Not allocated (currently).
127 :
128 : !To enable the use of ABINIT GW tools for constructing the square root of the coulomb operator.
129 : integer :: timrev
130 : real(dp) :: ecut_eff
131 : integer, pointer :: gvec(:,:)
132 : type(vcoul_t) :: Vcp
133 : type(crystal_t) :: Cryst
134 : type(gsphere_t) :: Gsphere
135 : type(kmesh_t) :: Kmesh, Qmesh
136 : character(len=132),pointer :: title(:) !SET2NULL
137 : complex(dp), allocatable :: vc_sqrt(:)
138 :
139 : !MPI over bands requires :
140 : integer :: blocksize !Public
141 : integer :: nbdblock !Public
142 : integer :: npw_g !Public
143 : integer :: ikpt_this_proc !Public
144 : integer :: npw_kb
145 : integer :: iblock
146 : integer :: iband
147 : integer :: n
148 : integer :: npw_serial
149 : integer, pointer :: kg_k_gather(:,:)
150 : real(dp), pointer :: kinpw_gather(:)
151 : real(dp), pointer :: ph3d_gather(:,:,:)
152 : real(dp), pointer :: ffnl_gather(:,:,:,:)
153 : !!***
154 :
155 :
156 : public :: dtset, mpi_enreg, cg, kinpw
157 : public :: nfft, n4, n5, n6, nline, tolwfr, npw_k, nbandv, eig, nspinor, ckpt
158 : public :: istwfk, nband, iovar, e, ucvol, vxc, ispden, ktot, kg_k
159 : public :: Hpsik, Hpsikc, g_to_r, gr_to_g, kbkb_to_kb
160 : public :: build_H, destroy_H, set_precondition, build_vxc
161 : public :: unset_precondition, precondition, exchange, dft_xc_energy, sqrt_vc_k, precondition_cplx
162 : public :: pcon, wf_block_distribute
163 : public :: blocksize, nbdblock, ikpt_this_proc, npw_g, npw_kb, kg_k_gather, kinpw_gather
164 :
165 : public :: dummy2, dummy3, mgfft, ngfft, weight, gbound, psir1, psir3, tim_fourwf
166 :
167 : real(dp), public, allocatable :: kernel_wavefunctions_FFT(:,:,:)
168 : real(dp), public, allocatable :: valence_wavefunctions_FFT(:,:,:)
169 :
170 : public :: pc_k_valence_kernel
171 : ! public :: CleanupSQMRKernel
172 : !!***
173 :
174 : contains
175 :
176 :
177 : !!****f* m_gwls_hamiltonian/DistributeValenceWavefunctions
178 : !! NAME
179 : !! DistributeValenceWavefunctions
180 : !!
181 : !! FUNCTION
182 : !! .
183 : !!
184 : !! INPUTS
185 : !!
186 : !! OUTPUT
187 : !!
188 : !! SOURCE
189 :
190 8 : subroutine DistributeValenceWavefunctions()
191 : !--------------------------------------------------------------------------------
192 : !
193 : ! This subroutine distributes, once and for all, the valence wavefunctions to
194 : ! the FFT configuration. They will thus be ready to be used within the
195 : ! susceptibility operator.
196 : !
197 : !--------------------------------------------------------------------------------
198 : integer :: iblk, mb, v
199 :
200 8 : real(dp), allocatable :: psik_v(:,:) !wavefunctions in LA format
201 8 : real(dp), allocatable :: psik_v_alltoall(:,:) !wavefunctions in FFT format
202 :
203 : ! *************************************************************************
204 :
205 :
206 :
207 : !===================================================================
208 : ! Allocate the global array which will contain the valence states
209 : !===================================================================
210 32 : ABI_MALLOC( valence_wavefunctions_FFT, (2,npw_g,nbdblock))
211 :
212 38618 : valence_wavefunctions_FFT = zero
213 :
214 :
215 : !====================================================
216 : ! Allocate working arrays
217 : !====================================================
218 :
219 24 : ABI_MALLOC(psik_v, (2,npw_kb))
220 24 : ABI_MALLOC(psik_v_alltoall,(2,npw_g))
221 :
222 :
223 :
224 : ! loop on all blocks of states,
225 68 : do iblk = 1, nbdblock
226 :
227 : ! loop on valence states for this block; if the state is conduction, fill with zeros
228 140 : do mb = 1, blocksize
229 :
230 80 : v = (iblk-1)*blocksize+mb
231 :
232 140 : if (v <= nbandv) then
233 15452 : psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(v-1)*npw_k+1:v*npw_k)
234 : else
235 23178 : psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = zero
236 : end if
237 :
238 : end do
239 :
240 : ! change configuration of the data
241 60 : call wf_block_distribute(psik_v, psik_v_alltoall,1) ! LA -> FFT
242 :
243 : ! copy data in global array
244 :
245 38618 : valence_wavefunctions_FFT(:,:,iblk) = psik_v_alltoall(:,:)
246 :
247 : end do
248 :
249 : !====================================================
250 : ! cleanup
251 : !====================================================
252 :
253 8 : ABI_FREE(psik_v)
254 8 : ABI_FREE(psik_v_alltoall)
255 :
256 :
257 8 : end subroutine DistributeValenceWavefunctions
258 : !!***
259 :
260 : !!****f* m_gwls_hamiltonian/DistributeValenceKernel
261 : !! NAME
262 : !! DistributeValenceKernel
263 : !!
264 : !! FUNCTION
265 : !! .
266 : !!
267 : !! INPUTS
268 : !!
269 : !! OUTPUT
270 : !!
271 : !! SOURCE
272 :
273 8 : subroutine DistributeValenceKernel()
274 : !--------------------------------------------------------------------------------
275 : !
276 : ! This subroutine distributes, once and for all, the kernel of the static
277 : ! SQMR operator.
278 : !
279 : ! In this first, quick and dirty implementation, we simply distribute
280 : ! ALL the valence bands on ALL the FFT groups. This is not efficient, but
281 : ! kernel projections are not a bottleneck of the computation.
282 : !
283 : ! A better (forthcoming) algorithm would only distribute the actual kernel,
284 : ! not all valence bands.
285 : !--------------------------------------------------------------------------------
286 :
287 : integer :: mb, n
288 :
289 8 : real(dp), allocatable :: psik_n(:,:) !wavefunctions in LA format
290 8 : real(dp), allocatable :: psik_n_alltoall(:,:) !wavefunctions in FFT format
291 :
292 : ! *************************************************************************
293 :
294 :
295 : !===================================================================
296 : ! Allocate the global array which will contain the valence states
297 : !===================================================================
298 32 : ABI_MALLOC( kernel_wavefunctions_FFT, (2,npw_g, nband))
299 :
300 46348 : kernel_wavefunctions_FFT = zero
301 :
302 : !====================================================
303 : ! Allocate working arrays
304 : !====================================================
305 :
306 24 : ABI_MALLOC(psik_n, (2,npw_kb))
307 24 : ABI_MALLOC(psik_n_alltoall,(2,npw_g))
308 :
309 :
310 : ! loop on all valence states
311 88 : do n = 1, nband
312 :
313 : ! Copy multiple instances of this valence state in the array;
314 : ! this way, each FFT group will have ALL the valence states!
315 200 : do mb = 1, blocksize
316 46460 : psik_n(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
317 : end do
318 :
319 : ! change configuration of the data
320 80 : call wf_block_distribute(psik_n, psik_n_alltoall,1) ! LA -> FFT
321 :
322 : ! copy data in global array
323 46348 : kernel_wavefunctions_FFT(:,:,n) = psik_n_alltoall(:,:)
324 :
325 : end do
326 :
327 : !====================================================
328 : ! cleanup
329 : !====================================================
330 :
331 8 : ABI_FREE(psik_n)
332 8 : ABI_FREE(psik_n_alltoall)
333 :
334 8 : end subroutine DistributeValenceKernel
335 : !!***
336 :
337 : !!****f* m_gwls_hamiltonian/pc_k_valence_kernel
338 : !! NAME
339 : !! pc_k_valence_kernel
340 : !!
341 : !! FUNCTION
342 : !! .
343 : !!
344 : !! INPUTS
345 : !!
346 : !! OUTPUT
347 : !!
348 : !! SOURCE
349 :
350 38570 : subroutine pc_k_valence_kernel(psi_inout,n)
351 :
352 : !================================================================================
353 : ! This routine projects out of the valence kernel. It assumes the input/output
354 : ! array is distributed in the FFT configuration, and that the global
355 : ! array containing the kernel (defined in this module) is already prepared
356 : ! and ready to be used.
357 : !================================================================================
358 :
359 : real(dp), intent(inout) :: psi_inout(2,npw_g)
360 : integer , intent(in), optional :: n
361 :
362 :
363 38570 : real(dp), allocatable :: psi_projected(:,:)
364 :
365 : real(dp) :: tmpc(2)
366 :
367 : integer :: v, n_max
368 :
369 : integer :: mpi_communicator, ierr
370 :
371 : ! *************************************************************************
372 :
373 :
374 115710 : ABI_MALLOC( psi_projected, (2,npw_g))
375 :
376 38570 : mpi_communicator = mpi_enreg%comm_fft
377 :
378 38570 : if (present(n)) then
379 0 : n_max = n
380 : else
381 38570 : n_max = nbandv
382 : end if
383 :
384 : !====================================================
385 : ! Compute projection
386 : !====================================================
387 :
388 23370572 : psi_projected(:,:) = zero
389 :
390 192850 : do v = 1, n_max
391 :
392 : ! compute overlap of kernel member with function
393 154280 : tmpc = cg_zdotc(npw_g ,kernel_wavefunctions_FFT(:,:,v),psi_inout(:,:))
394 :
395 : ! Communicate results
396 154280 : call xmpi_sum(tmpc, mpi_communicator, ierr) ! sum on all processors working on FFT!
397 :
398 : ! add overlap with array to project
399 31263616 : psi_projected(1,:) = psi_projected(1,:) + (tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v))
400 31302186 : psi_projected(2,:) = psi_projected(2,:) + (tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v))
401 :
402 : !psi_inout(1,:) = psi_inout(1,:) -( tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v) )
403 : !psi_inout(2,:) = psi_inout(2,:) -( tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v) )
404 :
405 :
406 : end do
407 :
408 38570 : if (present(n)) then
409 0 : psi_inout(:,:) = psi_projected(:,:)
410 : else
411 23370572 : psi_inout(:,:) = psi_inout(:,:) - psi_projected(:,:)
412 : end if
413 :
414 :
415 :
416 38570 : ABI_FREE( psi_projected)
417 :
418 38570 : end subroutine pc_k_valence_kernel
419 : !!***
420 :
421 : !!****f* m_gwls_hamiltonian/wf_block_distribute
422 : !! NAME
423 : !! wf_block_distribute
424 : !!
425 : !! FUNCTION
426 : !! .
427 : !!
428 : !! INPUTS
429 : !!
430 : !! OUTPUT
431 : !!
432 : !! SOURCE
433 :
434 6722 : subroutine wf_block_distribute(psik, psik_alltoall, direction)
435 :
436 : !================================================================================
437 : !
438 : ! This subroutine distributes a block of wavefunctions in the "linear algebra"
439 : ! configuration, to a configuration appropriate to perform FFT (or apply Hamiltonian).
440 : !
441 : ! The code below is inspired by the subroutine prep_getghc, which rearranges
442 : ! data over MPI processors prior to applying the Hamiltonian.
443 : !
444 : ! input:
445 : ! npw_kb : dimension of wavefunctions in the "linear algebra" configuration,
446 : ! which means the G-vectors are distributed over all
447 : ! processors, and every processor has information about all the bands.
448 : !
449 : !
450 : ! npw_g : dimension of wavefunctions in the "FFT" configuration,
451 : ! which means the G-vectors are distributed along rows of the
452 : ! MPI topology, and a given row only has one band.
453 : !
454 : ! direction : Are we going from LA to FFT, or from FFT to LA?
455 : ! 1 : LA -> FFT
456 : ! 2 : LA <- FFT
457 : !
458 : ! input/output:
459 : ! psik : block of wavefunctions, in LA configuration
460 : ! psik_alltoall : wavefunction, in FFT configuration
461 : !================================================================================
462 :
463 : integer , intent(in) :: direction ! flag which determines the direction of the transfer
464 : real(dp), intent(inout) :: psik(2,npw_kb) ! block of wavefunctions, in "linear algebra" configuration
465 : real(dp), intent(inout) :: psik_alltoall(2,npw_g) ! wavefunction in "FFT" configuration; a single band, but more G-vectors
466 :
467 :
468 : integer :: ier, spaceComm
469 :
470 6722 : integer,allocatable :: rdisplsloc(:)
471 6722 : integer,allocatable :: recvcountsloc(:)
472 6722 : integer,allocatable :: sdisplsloc(:)
473 6722 : integer,allocatable :: sendcountsloc(:)
474 :
475 : integer :: nproc_band, bandpp
476 :
477 6722 : integer,pointer :: rdispls(:)
478 6722 : integer,pointer :: recvcounts(:)
479 6722 : integer,pointer :: sdispls(:)
480 6722 : integer,pointer :: sendcounts(:)
481 :
482 : ! *************************************************************************
483 :
484 :
485 : ! extract information in order to perform MPI communication.
486 : ! This code comes from prep_getghc
487 6722 : nproc_band = mpi_enreg%nproc_band
488 6722 : bandpp = mpi_enreg%bandpp
489 :
490 6722 : if(mpi_enreg%nproc_band*mpi_enreg%bandpp > 1) then
491 :
492 2664 : spaceComm=mpi_enreg%comm_fft
493 2664 : if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_band
494 :
495 7992 : ABI_MALLOC(sendcountsloc,(nproc_band))
496 5328 : ABI_MALLOC(sdisplsloc ,(nproc_band))
497 5328 : ABI_MALLOC(recvcountsloc,(nproc_band))
498 5328 : ABI_MALLOC(rdisplsloc ,(nproc_band))
499 :
500 2664 : recvcounts =>bandfft_kpt(ikpt_this_proc)%recvcounts(:)
501 2664 : sendcounts =>bandfft_kpt(ikpt_this_proc)%sendcounts(:)
502 2664 : rdispls =>bandfft_kpt(ikpt_this_proc)%rdispls (:)
503 2664 : sdispls =>bandfft_kpt(ikpt_this_proc)%sdispls (:)
504 :
505 7992 : recvcountsloc(:)= recvcounts(:)*2*nspinor*bandpp
506 7992 : rdisplsloc(:) = rdispls(:)*2*nspinor*bandpp
507 7992 : sendcountsloc(:)= sendcounts(:)*2*nspinor
508 7992 : sdisplsloc(:) = sdispls(:)*2*nspinor
509 :
510 : ! use MPI to communicate information!
511 2664 : if (direction == 1) then
512 : ! LA -> FFT
513 652 : call xmpi_alltoallv(psik, sendcountsloc, sdisplsloc, psik_alltoall, recvcountsloc,rdisplsloc, spaceComm, ier)
514 :
515 2012 : else if (direction == 2) then
516 : ! FFT -> LA
517 :
518 2012 : call xmpi_alltoallv(psik_alltoall,recvcountsloc,rdisplsloc, psik, sendcountsloc,sdisplsloc,spaceComm,ier)
519 :
520 : end if
521 :
522 2664 : ABI_FREE(sendcountsloc)
523 2664 : ABI_FREE(sdisplsloc )
524 2664 : ABI_FREE(recvcountsloc)
525 2664 : ABI_FREE(rdisplsloc )
526 :
527 : else
528 :
529 674926 : if(direction == 1) psik_alltoall = psik
530 :
531 2465966 : if(direction == 2) psik = psik_alltoall
532 :
533 : end if
534 :
535 13444 : end subroutine wf_block_distribute
536 : !!***
537 :
538 :
539 :
540 : !!****f* m_gwls_hamiltonian/exchange
541 : !! NAME
542 : !! exchange
543 : !!
544 : !! FUNCTION
545 : !! .
546 : !!
547 : !! INPUTS
548 : !!
549 : !! OUTPUT
550 : !!
551 : !!
552 : !!
553 : !! SOURCE
554 :
555 9 : function exchange(e, Lbasis_lanczos)
556 :
557 : !use m_bandfft_kpt
558 : use m_cgtools
559 : !================================================================================
560 : ! This subroutine computes the exchange energy in band+FFT parallel
561 : !
562 : !================================================================================
563 : real(dp) :: exchange
564 :
565 : integer, intent(in) :: e
566 :
567 : ! If these arguments are provided, the exchange energy is to be projected on this subspace
568 : complex(dp), optional, intent(in) :: Lbasis_lanczos(:,:) ! complex array which contains the Lanczos basis
569 :
570 9 : real(dp), allocatable :: psik_e(:,:) !Working array to store the wavefunction
571 :
572 9 : real(dp), allocatable :: psik_v(:,:) !Working array to store the wavefunction
573 :
574 9 : real(dp), allocatable :: psik_out(:,:) !Working array to store the wavefunction
575 :
576 :
577 : integer :: iblk, mb
578 :
579 : integer :: l, lmax
580 :
581 : real(dp) :: tmpc(2)
582 :
583 : logical :: project
584 :
585 : ! *************************************************************************
586 :
587 : !--------------------------------------------------------------------------------
588 : ! Determine if the exhcange energy must be projected on the Lanczos basis
589 : ! a truncated Coulomb potential
590 : !--------------------------------------------------------------------------------
591 :
592 9 : project = .false.
593 9 : if (present(Lbasis_lanczos)) then
594 1 : project = .true.
595 1 : lmax = size(Lbasis_lanczos, 2)
596 : end if
597 :
598 : !--------------------------------------------------------------------------------
599 : ! The goal of this routine is to compute the exchange energy using
600 : ! a truncated Coulomb potential
601 : !--------------------------------------------------------------------------------
602 :
603 9 : exchange = 0.0_dp
604 :
605 : !====================================================
606 : ! build the block of wavefunctions which will
607 : ! contain copies the e-state
608 : !====================================================
609 :
610 27 : ABI_MALLOC(psik_e,(2,npw_kb))
611 :
612 : ! fill psik_e with as many copies of the e-state as there are band processors;
613 : ! that way, upon LA -> FFT, each row of fft processors will have the e-state
614 22 : do v = 1, blocksize
615 5419 : psik_e(:,(v-1)*npw_k+1:v*npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k)
616 : end do
617 :
618 : !====================================================
619 : ! Allocate the block of wavefunctions which will
620 : ! contain the valence states
621 : !====================================================
622 :
623 18 : ABI_MALLOC(psik_v, (2,npw_kb))
624 18 : ABI_MALLOC(psik_out, (2,npw_kb))
625 :
626 : ! loop on all blocks of states,
627 79 : do iblk = 1, nbdblock
628 :
629 : ! loop on valence states for this block; if the state is conduction, fill with zeros
630 160 : do mb = 1, blocksize
631 :
632 90 : v = (iblk-1)*blocksize+mb
633 :
634 160 : if (v <= nbandv) then
635 18540 : psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(v-1)*npw_k+1:v*npw_k)
636 : else
637 27810 : psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = zero
638 : end if
639 :
640 : end do
641 :
642 70 : call kbkb_to_kb(psik_out,psik_v,psik_e)
643 :
644 : ! apply Coulomb potential, and take norm: cumulate the exchange energy
645 169 : do mb = 1, blocksize
646 :
647 46440 : psik1 = psik_out(:,(mb-1)*npw_k+1:mb*npw_k)
648 90 : call sqrt_vc_k(psik1)
649 :
650 160 : if (project) then
651 : ! project on the Lanczos basis
652 90 : do l = 1, lmax
653 20640 : psik2(1,:) = dble(Lbasis_lanczos(:,l))
654 20640 : psik2(2,:) = dimag(Lbasis_lanczos(:,l))
655 80 : tmpc = scprod_k(psik2,psik1)
656 90 : exchange = exchange - (tmpc(1)**2+tmpc(2)**2)
657 : end do
658 : else
659 : ! compute the "exact" exchange energy
660 80 : exchange = exchange - norm_k(psik1)**2
661 : end if
662 :
663 : end do
664 :
665 : end do
666 :
667 9 : ABI_FREE(psik_e)
668 :
669 9 : ABI_FREE(psik_v)
670 :
671 9 : ABI_FREE(psik_out)
672 :
673 9 : end function exchange
674 : !!***
675 :
676 :
677 : !!****f* m_gwls_hamiltonian/dft_xc_energy
678 : !! NAME
679 : !! dft_xc_energy
680 : !!
681 : !! FUNCTION
682 : !! .
683 : !!
684 : !! INPUTS
685 : !!
686 : !! OUTPUT
687 : !!
688 : !!
689 : !!
690 : !! SOURCE
691 :
692 8 : function dft_xc_energy(e)
693 :
694 : real(dp) :: dft_xc_energy
695 : integer, intent(in) :: e
696 :
697 :
698 : integer :: cplex, option, ierr
699 8 : real(dp), allocatable :: psik_e(:,:) !Working array to store the wavefunction
700 8 : real(dp), allocatable :: psik_e_alltoall(:,:) !Working array to store the wavefunction
701 :
702 8 : real(dp), allocatable :: psik_out(:,:) !Working array to store the wavefunction
703 :
704 : real(dp) :: tmpc(2)
705 : integer :: mpi_communicator
706 : real(dp) :: dft_xc_energy_tmp
707 :
708 : ! *************************************************************************
709 :
710 : !--------------------------------------------------------------------------------
711 : ! The goal of this routine is to compute the xc energy using
712 : ! the DFT Vxc array.
713 : !--------------------------------------------------------------------------------
714 :
715 : ! Parallel-MPI code. This is inspired from code within the getghc routine, which
716 : ! applies a complex potential to a wavefunction.
717 :
718 :
719 8 : cplex = 1 ! real potential
720 8 : option = 2 ! multiply wavefunction by potential
721 8 : nspinor= 1
722 : ! Allocate wavefunction arrays which contains the coefficients of the e-state
723 24 : ABI_MALLOC(psik_e, (2,npw_kb))
724 24 : ABI_MALLOC(psik_e_alltoall,(2,npw_g))
725 24 : ABI_MALLOC(psik_out, (3,npw_g))
726 :
727 :
728 : ! Only fetch the e-state, setting other states in block to zero!
729 4634 : psik_e(:,:) = zero
730 3863 : psik_e(:,1:npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k)
731 6176 : psik_out(:,:) = zero
732 :
733 :
734 : ! change the configuration of the data
735 8 : call wf_block_distribute(psik_e, psik_e_alltoall,1) ! LA -> FFT
736 :
737 :
738 : ! Call fourwf to generate the product, in k-space
739 : ! Computation:
740 : ! psik_e_alltoall(k) -> psi(r)
741 : ! res(r) = (vxc(r) x psi(r))
742 : ! psik_out(k) <- res(r)
743 : ! psir3 is a dummy, not used here.
744 :
745 : call fourwf(cplex,vxc(:,:,:,ispden),psik_e_alltoall,psik_out,psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,&
746 8 : & mpi_enreg,1,ngfft,npw_g,npw_g,n4,n5,n6,option,tim_fourwf,weight,weight)
747 :
748 :
749 8 : tmpc = cg_zdotc(npw_g, psik_e_alltoall,psik_out)
750 :
751 8 : mpi_communicator = mpi_enreg%comm_fft
752 8 : call xmpi_sum(tmpc,mpi_communicator, ierr) ! sum on all processors working on FFT!
753 :
754 8 : dft_xc_energy_tmp = tmpc(1)
755 :
756 8 : mpi_communicator = mpi_enreg%comm_band
757 8 : call xmpi_sum(dft_xc_energy_tmp,mpi_communicator, ierr) ! sum on all processors working on FFT!
758 :
759 8 : dft_xc_energy = dft_xc_energy_tmp
760 :
761 8 : ABI_FREE(psik_e)
762 8 : ABI_FREE(psik_e_alltoall)
763 8 : ABI_FREE(psik_out)
764 :
765 8 : end function dft_xc_energy
766 : !!***
767 :
768 : !!****f* m_gwls_hamiltonian/set_precondition
769 : !! NAME
770 : !! set_precondition
771 : !!
772 : !! FUNCTION
773 : !! .
774 : !!
775 : !! INPUTS
776 : !!
777 : !! OUTPUT
778 : !!
779 : !! SOURCE
780 :
781 736 : subroutine set_precondition(lambda,omega)
782 : !--------------------------------------------------------------------------------
783 : ! This subroutine preconditions the problem
784 : !
785 : ! A x = b,
786 : !
787 : ! with:
788 : !
789 : ! omega lambda Operator
790 : ! ------------------------------------------------------
791 : ! absent absent A = (H - lambda_0) (value of lambda_0 not important)
792 : ! present present A = (H - lambda)^2 + omega^2
793 : ! other cases not implemented
794 : !
795 : !
796 : ! In the above, b = psik.
797 : !
798 : !--------------------------------------------------------------------------------
799 :
800 : ! TODO :
801 : ! - eliminate the 2 "if(kinpw(i) < huge(zero)*1.0d-11)"
802 : ! since ecutsm = 0.0 always (check if that's true in this gw_sternheimer subroutine).
803 :
804 : real(dp), intent(in), optional :: lambda, omega
805 :
806 : real(dp) :: poly, x
807 : logical :: omega_imaginary
808 :
809 :
810 : !integer,save :: counter = 0
811 : !integer :: io_unit
812 : !character(18):: filename
813 : !logical :: file_exists
814 :
815 : ! *************************************************************************
816 :
817 :
818 736 : if (present(lambda) .and. present(omega)) then
819 : omega_imaginary =.true.
820 : else
821 736 : omega_imaginary =.false.
822 : end if
823 :
824 : !io_unit = get_unit()
825 : !filename = "PRECONDITIONER.log"
826 : !inquire(file=filename,exist=file_exists)
827 :
828 : !if (file_exists) then
829 : ! open(io_unit,file=filename,position='append',status=files_status_old)
830 : !else
831 : ! open(io_unit,file=filename,status=files_status_new)
832 : ! write(io_unit,10) "#======================================================================================="
833 : ! write(io_unit,10) "# "
834 : ! write(io_unit,10) "# This file contains information regarding the preconditioning scheme for SQMR. "
835 : ! write(io_unit,10) "# "
836 : ! write(io_unit,10) "#======================================================================================="
837 : !end if
838 :
839 : !counter = counter + 1
840 : !write(io_unit,10) "# "
841 : !write(io_unit,15) "# Call #:", counter
842 : !if (present(lambda) .and. present(omega)) then
843 : ! write(io_unit,10) "# lambda and omega are present: imaginary frequency case"
844 : !else
845 : ! write(io_unit,10) "# lambda and omega are absent: omega = 0 case"
846 : !end if
847 :
848 : !do i=1,npw_k
849 :
850 148768 : do i = 1, npw_g
851 :
852 148032 : if(omega_imaginary) then
853 0 : x = (kinpw_gather(i)-lambda)**2 + omega**2
854 : else
855 148032 : x = kinpw_gather(i)
856 : end if
857 :
858 148768 : if(x < huge(zero)*1.0d-11) then
859 148032 : poly = 27.0 + x*(18.0 + x*(12.0 + 8.0*x))
860 148032 : pcon(i) = poly/(poly + 16.0*(x**4))
861 : !pcon(i) = 1.0/(1.0+x) !I don't know why, it gives better results for Silane than the above polynomial.
862 : else
863 0 : pcon(i) = zero
864 : end if
865 : end do
866 :
867 : !write(io_unit,30) " prec( 1: 10) = ",pcon(1:10)
868 : !write(io_unit,30) " prec(npw_k-10:npw_k) = ",pcon(npw_k-10:npw_k)
869 :
870 : !close(io_unit)
871 :
872 : !10 format(A)
873 : !15 format(A,I8)
874 : !30 format(A,1000(ES24.12))
875 :
876 736 : end subroutine set_precondition
877 : !!***
878 :
879 : !!****f* m_gwls_hamiltonian/unset_precondition
880 : !! NAME
881 : !! unset_precondition
882 : !!
883 : !! FUNCTION
884 : !! .
885 : !!
886 : !! INPUTS
887 : !!
888 : !! OUTPUT
889 : !!
890 : !! SOURCE
891 :
892 128 : subroutine unset_precondition()
893 :
894 : ! *************************************************************************
895 :
896 33024 : pcon = one
897 128 : end subroutine unset_precondition
898 : !!***
899 :
900 : !!****f* m_gwls_hamiltonian/precondition
901 : !! NAME
902 : !! precondition
903 : !!
904 : !! FUNCTION
905 : !! .
906 : !!
907 : !! INPUTS
908 : !!
909 : !! OUTPUT
910 : !!
911 : !! SOURCE
912 :
913 16285 : subroutine precondition(psi_out,psi_in)
914 :
915 : real(dp), intent(out) :: psi_out(2,npw_g)
916 : real(dp), intent(in) :: psi_in(2,npw_g)
917 :
918 : ! *************************************************************************
919 :
920 3274274 : do i=1,npw_g
921 9790252 : psi_out(:,i) = psi_in(:,i)*pcon(i)
922 : end do
923 :
924 16285 : end subroutine precondition
925 : !!***
926 :
927 : !!****f* m_gwls_hamiltonian/precondition_cplx
928 : !! NAME
929 : !! precondition_cplx
930 : !!
931 : !! FUNCTION
932 : !! .
933 : !!
934 : !! INPUTS
935 : !!
936 : !! OUTPUT
937 : !!
938 : !! SOURCE
939 :
940 7352 : subroutine precondition_cplx(psi_out,psi_in)
941 :
942 : complex(dp), intent(out) :: psi_out(npw_g)
943 : complex(dp), intent(in) :: psi_in(npw_g)
944 :
945 : ! *************************************************************************
946 :
947 1896816 : do i=1,npw_g
948 1896816 : psi_out(i) = psi_in(i)*pcon(i)
949 : end do
950 7352 : end subroutine precondition_cplx
951 : !!***
952 :
953 : !!****f* m_gwls_hamiltonian/sqrt_vc_k
954 : !! NAME
955 : !! sqrt_vc_k
956 : !!
957 : !! FUNCTION
958 : !! .
959 : !!
960 : !! INPUTS
961 : !!
962 : !! OUTPUT
963 : !!
964 : !! SOURCE
965 :
966 1165 : subroutine sqrt_vc_k(psi_inout)
967 :
968 : !External variables
969 : real(dp), intent(inout) :: psi_inout(2,npw_k)
970 :
971 : !Internal variable
972 : complex(dp) :: c
973 :
974 : ! *************************************************************************
975 :
976 177210 : do i=1,npw_k
977 176045 : c = vc_sqrt(i) * cmplx(psi_inout(1,i),psi_inout(2,i),dp)
978 176045 : psi_inout(1,i) = dble (c)
979 177210 : psi_inout(2,i) = dimag(c)
980 : end do
981 1165 : end subroutine sqrt_vc_k
982 : !!***
983 :
984 : !!****f* m_gwls_hamiltonian/Hpsik
985 : !! NAME
986 : !! Hpsik
987 : !!
988 : !! FUNCTION
989 : !! .
990 : !!
991 : !! INPUTS
992 : !!
993 : !! OUTPUT
994 : !!
995 : !! SOURCE
996 :
997 17569 : subroutine Hpsik(psi_out,psi_in,cte)
998 :
999 : !External variables
1000 : real(dp), intent(inout) :: psi_out(2,npw_g)
1001 : real(dp), intent(inout), optional :: psi_in(2,npw_g)
1002 : real(dp), intent(in), optional :: cte
1003 :
1004 : ! *************************************************************************
1005 :
1006 : !psikb1 is allocated in init_hamiltonian and destroyed in destroy_hamiltonian, to avoid doing it at each application of
1007 : !the hamiltonian...
1008 :
1009 : !We are keeping track of how the module treat each argument of the hamiltonian (subroutine getghc) here.
1010 : !Legend : T = transmitted from 79_seqpar_mpi/vtowfk.F90 through the call gwls_sternheimer
1011 : ! I = input argument (allocated and filled)
1012 : ! O = output argument (allocated)
1013 : ! D = dummy argument (declared but not allocated. Read/Write attempts should trigger a segfault.)
1014 : !call getghc(cpopt, T
1015 : !psi, I
1016 : !conjgrprj, D
1017 : !Hpsik, O
1018 : !gsc, D (output for PAW : <G|S|C>)
1019 : !gs_hamk, T
1020 : !gvnlxc, D (<G|Vnonlocal+VFockACE|C>)
1021 : !eshift, I (<G|H-eshift.S|C>)
1022 : !mpi_enreg, T
1023 : !ndat, Fixed to 1 (# of FFTs to do in //)
1024 : !dtset%prtvol, T
1025 : !sij_opt, Fixed to 0 (PAW dependent : 0-><G|H|C> ; 1-><G|H|C> & <G|S|C> ; -1-><G|H-cte.S|C>)
1026 : !tim_getghc, Fixed to 0 (identity of the timer of the calling subroutine. 1:cgwf 5:lobpcg)
1027 : !type_calc) Fixed to 0 (0:whole H N:some part only)
1028 :
1029 : ! It is necessary to pass this optional argument to getghc if paral_kgb=1, or else the code exits cleanly with a "BUG" message...
1030 : ! It will be important to properly understand what this means when we start using kpoints...
1031 :
1032 26977 : if(present(psi_in)) then
1033 :
1034 : call getghc(cpopt,psi_in,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,&
1035 17569 : & mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1036 :
1037 : else
1038 6925434 : psig3 = psi_out
1039 :
1040 : call getghc(cpopt,psig3,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,&
1041 9408 : & mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1042 :
1043 : end if
1044 26977 : if(present(cte)) then
1045 17821 : if(present(psi_in)) then
1046 10531628 : psi_out = psi_out - cte*psi_in
1047 : else
1048 197696 : psi_out = psi_out - cte*psig3
1049 : end if
1050 : end if
1051 :
1052 26977 : end subroutine Hpsik
1053 : !!***
1054 :
1055 : !!****f* m_gwls_hamiltonian/Hpsikc
1056 : !! NAME
1057 : !! Hpsikc
1058 : !!
1059 : !! FUNCTION
1060 : !! .
1061 : !!
1062 : !! INPUTS
1063 : !!
1064 : !! OUTPUT
1065 : !!
1066 : !! SOURCE
1067 :
1068 7352 : subroutine Hpsikc(psi_out,psi_in,cte)
1069 :
1070 : !External variables
1071 : complex(dp), intent(out) :: psi_out(npw_g)
1072 : complex(dp), intent(in) :: psi_in(npw_g)
1073 : complex(dp), intent(in), optional :: cte
1074 :
1075 : ! *************************************************************************
1076 :
1077 :
1078 1896816 : psig1(1,:) = dble(psi_in)
1079 1896816 : psig1(2,:) = dimag(psi_in)
1080 :
1081 7352 : call Hpsik(psig1)
1082 :
1083 1896816 : psi_out = dcmplx(psig1(1,:),psig1(2,:))
1084 :
1085 1896816 : if(present(cte)) psi_out = psi_out - cte*psi_in
1086 7352 : end subroutine Hpsikc
1087 : !!***
1088 :
1089 : !!!****f* m_gwls_hamiltonian/pc_k
1090 : !!! NAME
1091 : !!! pc_k
1092 : !!!
1093 : !!! FUNCTION
1094 : !!! .
1095 : !!!
1096 : !!! INPUTS
1097 : !!!
1098 : !!! OUTPUT
1099 : !!!
1100 : !!! SOURCE
1101 : !
1102 : !subroutine pc_k(psi_inout,n,eig_e,above)
1103 : !
1104 : !real(dp), intent(inout) :: psi_inout(2,npw_kb)
1105 : !integer , intent(in), optional :: n
1106 : !real(dp), intent(in), optional :: eig_e
1107 : !logical, intent(in), optional :: above !Has an effect only if n is also given in argument
1108 : !
1109 : !!Local variables
1110 : !real(dp),parameter :: degeneracy_tolerance = 2.0e-8
1111 : !real(dp),parameter :: projection_tolerance = 1.0e-16
1112 : !real(dp) :: z(2)
1113 : !
1114 : !integer :: mpi_communicator
1115 : !
1116 : !! *************************************************************************
1117 : !
1118 : !
1119 : !mpi_communicator = mpi_enreg%comm_bandfft
1120 : !
1121 : !if( present(eig_e) ) then
1122 : ! do i = 1, nband
1123 : ! if (abs(eig_e-eig(i)) < degeneracy_tolerance) then
1124 : ! z(:) = scprod_k(cg(:,(i-1)*npw_k+1:i*npw_k),psi_inout)
1125 : ! ! project it out!
1126 : ! if ( sqrt(z(1)**2+z(2)**2) > projection_tolerance) then
1127 : ! ! only project if z is large enough; if we project when z is
1128 : ! ! very small, we introduce roundoff error!
1129 : ! psi_inout(1,1:npw_k) = psi_inout(1,:) - ( z(1)*cg(1,(i-1)*npw_k+1:i*npw_k)-z(2)*cg(2,(i-1)*npw_k+1:i*npw_k) )
1130 : ! psi_inout(2,1:npw_k) = psi_inout(2,:) - ( z(1)*cg(2,(i-1)*npw_k+1:i*npw_k)+z(2)*cg(1,(i-1)*npw_k+1:i*npw_k) )
1131 : ! end if
1132 : ! end if
1133 : ! end do
1134 : ! elseif( present(n) ) then
1135 : ! !If there is a state "n" given in argument, then project on subspace <=n.
1136 : !
1137 : ! if ( n == 0) then
1138 : ! ! if n = 0, there are no states with a lower energy! Thus, the projection vanishes.
1139 : ! if(.not. (present(above) .and. above)) psi_inout = zero
1140 : ! else
1141 : ! if(present(above) .and. above) then
1142 : ! !call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,mpi_enreg,0,n,npw_k,nspinor,dtset%ortalg,1,dummy2,scprod2,0,0,0)
1143 : ! call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,0,n,npw_k,nspinor,dummy2,scprod2,0,0,0,mpi_enreg%me_g0,mpi_communicator)
1144 : ! else
1145 : ! psikb4 = psi_inout
1146 : ! !call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,mpi_enreg,0,n,npw_k,nspinor,dtset%ortalg,1,dummy2,scprod2,0,0,0)
1147 : ! call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,0,n,npw_k,nspinor,dummy2,scprod2,0,0,0,mpi_enreg%me_g0,mpi_communicator)
1148 : ! psi_inout = psikb4 - psi_inout
1149 : ! end if
1150 : ! end if
1151 : !else
1152 : ! !If there is no "n" given in argument, then project on conduction states.
1153 : ! !call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,mpi_enreg,0,nbandv,npw_k,nspinor,dtset%ortalg,1,dummy2,scprod2,0,0,0)
1154 : ! call projbd(cg,psi_inout,-1,0,0,istwfk(ckpt),mcg,0,nbandv,npw_k,nspinor,dummy2,scprod2,0,0,0,mpi_enreg%me_g0,mpi_communicator)
1155 : !end if
1156 : !end subroutine pc_k
1157 : !!!***
1158 :
1159 : !!****f* m_gwls_hamiltonian/g_to_r
1160 : !! NAME
1161 : !! g_to_r
1162 : !!
1163 : !! FUNCTION
1164 : !! .
1165 : !!
1166 : !! INPUTS
1167 : !!
1168 : !! OUTPUT
1169 : !!
1170 : !! SOURCE
1171 :
1172 3385 : subroutine g_to_r(psi_out,psi_in)
1173 :
1174 : real(dp), intent(out) :: psi_out(2,n4,n5,n6)
1175 : real(dp), intent(in) :: psi_in(2,npw_g)
1176 : integer :: option, cplex
1177 : integer :: nproc_fft, me_fft, nd3 !, ierr
1178 :
1179 : ! *************************************************************************
1180 :
1181 3385 : option = 0 ! fft wavefunction to real space
1182 3385 : cplex = 2 ! complex potential
1183 :
1184 2132417 : psig4 = psi_in
1185 67209175 : psi_out = zero
1186 : call fourwf(cplex,dummy3,psig4,dummy2,psi_out,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg, &
1187 3385 : 1,ngfft,npw_g,npw_g,n4,n5,n6,option,tim_fourwf,weight,weight)
1188 67209175 : psi_out = psi_out/sqrt(ucvol)
1189 :
1190 : !! This comes from prep_fourwf
1191 : !nproc_fft=mpi_enreg%nproc_fft
1192 : !if (nproc_fft>1) then
1193 : ! me_fft=mpi_enreg%me_fft
1194 : ! if (me_fft>0) then
1195 : ! nd3=(ngfft(3)-1)/nproc_fft+1
1196 : ! psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3)
1197 : ! psi_out(:,:,:,1:nd3)=zero
1198 : ! end if
1199 : ! call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr)
1200 : !end if
1201 :
1202 : !Instead of communications the whole real space vector on all comm_fft CPUs,
1203 : !it is possible to let each CPU keep it's part only and communicate only the sums.
1204 : !Then, however, we need to clean the trash in the part of the real space
1205 : !vector that's not attributed to the given comm_fft CPU.
1206 3385 : nproc_fft=mpi_enreg%nproc_fft
1207 3385 : if (nproc_fft>1) then
1208 1256 : me_fft=mpi_enreg%me_fft
1209 1256 : if (me_fft>0) then
1210 628 : nd3=(ngfft(3)-1)/nproc_fft+1
1211 12468940 : psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3)
1212 6234784 : psi_out(:,:,:,1:nd3)=zero
1213 : end if
1214 : !!! call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr)
1215 : end if
1216 :
1217 3385 : end subroutine g_to_r
1218 : !!***
1219 :
1220 : !!****f* m_gwls_hamiltonian/gr_to_g
1221 : !! NAME
1222 : !! gr_to_g
1223 : !!
1224 : !! FUNCTION
1225 : !! .
1226 : !!
1227 : !! INPUTS
1228 : !!
1229 : !! OUTPUT
1230 : !!
1231 : !! SOURCE
1232 :
1233 5199 : subroutine gr_to_g(psig_out,psir_in,psig_in)
1234 :
1235 : real(dp), intent(in) :: psir_in(2,n4,n5,n6)
1236 : real(dp), intent(in), optional :: psig_in(2,npw_g)
1237 : real(dp), intent(out) :: psig_out(2,npw_g)
1238 :
1239 : integer :: i1, i2, i3
1240 : integer :: cplex, option
1241 :
1242 : ! *************************************************************************
1243 :
1244 5269 : cplex = 2 ! complex potential
1245 5269 : option= 2 ! multiply wavefunction by potential
1246 :
1247 5269 : if(.not. present(psig_in)) then
1248 46330 : psig4(:,:) = psig_out(:,:)
1249 : else
1250 3302766 : psig4(:,:) = psig_in(:,:)
1251 : end if
1252 :
1253 100111 : do i3=1,n6
1254 1902109 : do i2=1,n5
1255 36134802 : do i1=1,n4
1256 34237962 : denpot(2*i1-1,i2,i3)= psir_in(1,i1,i2,i3)
1257 36039960 : denpot(2*i1 ,i2,i3)= psir_in(2,i1,i2,i3)
1258 : end do
1259 : end do
1260 : end do
1261 :
1262 : call fourwf( cplex, & ! complex potential
1263 : denpot, & ! real space wavefunction, in denpot format
1264 : psig4, & ! fourier space wavefunction
1265 : psig_out, & ! result, in FFT configuration
1266 : psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg,1, & ! Various other arguments
1267 5269 : ngfft,npw_g,npw_g,n4,n5,n6,option,tim_fourwf,weight,weight)
1268 :
1269 5269 : end subroutine gr_to_g
1270 : !!***
1271 :
1272 : !!****f* m_gwls_hamiltonian/kbkb_to_kb
1273 : !! NAME
1274 : !! kbkb_to_kb
1275 : !!
1276 : !! FUNCTION
1277 : !! .
1278 : !!
1279 : !! INPUTS
1280 : !!
1281 : !! OUTPUT
1282 : !!
1283 : !! SOURCE
1284 :
1285 70 : subroutine kbkb_to_kb(psik_out,psik_in_1,psik_in_2)
1286 : !----------------------------------------------------------------------------------------------------
1287 : ! This function computes the direct product of two wavefunctions in real space,
1288 : ! psi_out(r) = psi_in_1^*(r)*psi_in_2(r) (without complex conjugating any of the 2 wavefunctions)
1289 : ! and returns the result in fourier space, psi_out(k).
1290 : !
1291 : ! The two input wavefunctions are in k space.
1292 : !
1293 : !
1294 : !----------------------------------------------------------------------------------------------------
1295 : real(dp), intent(out) :: psik_out(2,npw_kb)
1296 : real(dp), intent(inout) :: psik_in_1(2,npw_kb), psik_in_2(2,npw_kb)
1297 :
1298 : ! *************************************************************************
1299 :
1300 : ! change configuration of the data : LA -> FFT
1301 70 : call wf_block_distribute(psik_in_1, psig1,1)
1302 70 : call wf_block_distribute(psik_in_2, psig2,1)
1303 :
1304 : ! Fourier transform the first input
1305 70 : call g_to_r(psir1, psig1)
1306 :
1307 : ! Complex conjugate in real space the first input
1308 480130 : psir1(2,:,:,:) = -psir1(2,:,:,:)
1309 :
1310 : ! Fourrier transform the second input, multiply component-wise
1311 70 : call gr_to_g(psig2,psir1)
1312 :
1313 : ! change configuration of the output : FFT -> LA
1314 70 : call wf_block_distribute(psik_out, psig2,2)
1315 :
1316 70 : end subroutine kbkb_to_kb
1317 : !!***
1318 :
1319 : !!****f* m_gwls_hamiltonian/build_vxc
1320 : !! NAME
1321 : !! build_vxc
1322 : !!
1323 : !! FUNCTION
1324 : !! .
1325 : !!
1326 : !! INPUTS
1327 : !!
1328 : !! OUTPUT
1329 : !!
1330 : !! SOURCE
1331 :
1332 8 : subroutine build_vxc(vxc2,nfft2,nspden2)
1333 : !Only transcribe the argument vxc2 in the module; the change from dg to sg is done in build_H (and stored in vxc), since the
1334 : !arguments of fftpac are built in build_H.
1335 :
1336 : !We need the dimensions of vxc since they don't exist yet in the module; build_vxc being called before build_H.
1337 : integer, intent(in) :: nfft2, nspden2
1338 : real(dp), intent(in) :: vxc2(nfft2,nspden2)
1339 :
1340 : ! *************************************************************************
1341 :
1342 32 : ABI_MALLOC(vxc_dg,(nfft2,nspden2))
1343 35016 : vxc_dg = vxc2
1344 :
1345 8 : end subroutine build_vxc
1346 : !!***
1347 :
1348 : !!****f* m_gwls_hamiltonian/destroy_H
1349 : !! NAME
1350 : !! destroy_H
1351 : !!
1352 : !! FUNCTION
1353 : !! .
1354 : !!
1355 : !! INPUTS
1356 : !!
1357 : !! OUTPUT
1358 : !!
1359 : !! SOURCE
1360 :
1361 8 : subroutine destroy_H
1362 :
1363 8 : call dtset%free()
1364 8 : call gs_hamk%free()
1365 :
1366 8 : call cryst%free()
1367 8 : call Kmesh%free()
1368 8 : call Qmesh%free()
1369 8 : call Gsphere%free()
1370 8 : call Vcp%free()
1371 :
1372 8 : call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
1373 8 : call destroy_mpi_enreg(mpi_enreg)
1374 :
1375 8 : ABI_SFREE(cg)
1376 8 : ABI_SFREE(gbound)
1377 8 : ABI_SFREE(kg_k)
1378 8 : ABI_SFREE(ffnl)
1379 8 : ABI_SFREE(ph3d)
1380 8 : ABI_SFREE(kinpw)
1381 8 : ABI_SFREE(vxc)
1382 8 : ABI_SFREE(vlocal)
1383 8 : ABI_SFREE(conjgrprj)
1384 8 : ABI_SFREE(istwfk)
1385 8 : ABI_SFREE(dummy2)
1386 8 : ABI_SFREE(dummy3)
1387 8 : ABI_SFREE(eig)
1388 8 : ABI_SFREE(scprod2)
1389 8 : ABI_SFREE(pcon)
1390 8 : ABI_SFREE(psik1)
1391 8 : ABI_SFREE(psik2)
1392 8 : ABI_SFREE(psik3)
1393 8 : ABI_SFREE(psik4)
1394 8 : ABI_SFREE(psikb1)
1395 8 : ABI_SFREE(psikb2)
1396 8 : ABI_SFREE(psikb3)
1397 8 : ABI_SFREE(psikb4)
1398 8 : ABI_SFREE(psig1)
1399 8 : ABI_SFREE(psig2)
1400 8 : ABI_SFREE(psig3)
1401 8 : ABI_SFREE(psig4)
1402 8 : ABI_SFREE(psir1)
1403 8 : ABI_SFREE(psir2)
1404 8 : ABI_SFREE(psir3)
1405 8 : ABI_SFREE(psidg)
1406 8 : ABI_SFREE(vxc_dg)
1407 8 : ABI_SFREE(denpot)
1408 8 : ABI_SFREE(kernel_wavefunctions_FFT)
1409 8 : ABI_SFREE(valence_wavefunctions_FFT)
1410 8 : if(associated(gvec)) then
1411 8 : ABI_FREE(gvec)
1412 : end if
1413 8 : ABI_SFREE(vc_sqrt)
1414 :
1415 8 : end subroutine destroy_H
1416 : !!***
1417 :
1418 : !!****f* m_gwls_hamiltonian/build_H
1419 : !! NAME
1420 : !! build_H
1421 : !!
1422 : !! FUNCTION
1423 : !! Arguments of gw_sternheimer, received as argument by build_H
1424 : !!
1425 : !! INPUTS
1426 : !!
1427 : !! OUTPUT
1428 : !!
1429 : !! SOURCE
1430 :
1431 8 : subroutine build_H(dtset2,mpi_enreg2,cpopt2,cg2,gs_hamk2,kg_k2,kinpw2)
1432 :
1433 : use m_wfutils
1434 :
1435 : type(dataset_type), intent(in) :: dtset2
1436 : type(MPI_type), intent(in) :: mpi_enreg2
1437 : type(gs_hamiltonian_type), intent(inout) :: gs_hamk2
1438 :
1439 : integer, intent(in) :: cpopt2
1440 : integer, intent(in) :: kg_k2(3,gs_hamk2%npw_k)
1441 : !Since mcg = npw_k*nband when there is only gamma, the next line works. If there is not only Gamma, mcg is the proper size of cg.
1442 : real(dp), intent(in) :: cg2(2,dtset2%mpw*dtset2%nspinor*dtset2%mband*dtset2%mkmem*dtset2%nsppol)
1443 : real(dp), intent(in) :: kinpw2(gs_hamk2%npw_k)
1444 : !Local variables : most of them are in the module for now.
1445 : real(dp) :: commutation_error
1446 : integer :: io_unit
1447 : integer :: io_unit_debug
1448 : integer :: ierr
1449 : integer :: cplx
1450 : integer :: j, k
1451 : integer :: dimph3d
1452 : integer :: mb
1453 :
1454 : integer :: mpi_communicator
1455 :
1456 : character(128):: filename_debug
1457 :
1458 :
1459 8 : real(dp), allocatable :: wfk_tmp1(:,:) ,wfk_tmp2(:,:)
1460 :
1461 : ! *************************************************************************
1462 :
1463 : ! Hartree-Fock cannot be used with GWLS.
1464 8 : if(dtset2%usefock==1 .and. associated(gs_hamk2%fockcommon)) then
1465 0 : ABI_ERROR(' build_H : Hartree-Fock option can not be used with optdriver==66 (GWLS calculations).')
1466 : end if
1467 :
1468 : !First we copy the data structure types
1469 8 : dtset = dtset2%copy()
1470 :
1471 8 : call copy_mpi_enreg(mpi_enreg2,mpi_enreg)
1472 :
1473 8 : call gs_hamk2%copy(gs_hamk)
1474 :
1475 : !Then we copy the standard types
1476 8 : cpopt = cpopt2
1477 8 : npw_k = gs_hamk2%npw_k
1478 :
1479 24 : ABI_MALLOC(cg,(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol))
1480 38566 : cg = cg2
1481 :
1482 48 : ABI_MALLOC(vlocal,(gs_hamk2%n4,gs_hamk2%n5,gs_hamk2%n6,gs_hamk2%nvloc))
1483 54888 : vlocal = gs_hamk2%vlocal
1484 8 : gs_hamk%vlocal => vlocal
1485 :
1486 24 : ABI_MALLOC(kg_k,(3,npw_k))
1487 5156 : kg_k = kg_k2
1488 24 : ABI_MALLOC(kinpw,(npw_k))
1489 1301 : kinpw = kinpw2
1490 :
1491 8 : dimffnl=0; if (blocksize<=1) dimffnl=size(gs_hamk2%ffnl_k,2)
1492 48 : ABI_MALLOC(ffnl,(npw_k,dimffnl,gs_hamk%lmnmax,gs_hamk%ntypat))
1493 8 : dimph3d=0; if (blocksize<=1) dimph3d=gs_hamk2%matblk
1494 32 : ABI_MALLOC(ph3d,(2,npw_k,dimph3d))
1495 :
1496 : !Initializing variables from dataset
1497 8 : nfft = dtset%nfft
1498 8 : nline = dtset%nline
1499 8 : tolwfr = dtset%tolwfr
1500 152 : ngfft = dtset%ngfft
1501 8 : mcg = dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
1502 8 : nspinor=dtset%nspinor
1503 8 : n1=ngfft(1)
1504 8 : n2=ngfft(2)
1505 8 : n3=ngfft(3)
1506 8 : n4=ngfft(4)
1507 8 : n5=ngfft(5)
1508 8 : n6=ngfft(6)
1509 32 : mgfft=maxval(ngfft(1:3))
1510 8 : nbandv = int(dtset%nelect)/2
1511 24 : ABI_MALLOC(istwfk,(dtset%nkpt))
1512 16 : istwfk(:)=dtset%istwfk
1513 :
1514 : !Initializing variables from gs_hamk
1515 8 : ucvol = gs_hamk%ucvol
1516 32 : ABI_MALLOC(gbound,(2*mgfft+8,2))
1517 : !gbound = gs_hamk%gbound_k !Must be done later for bandft parallelism
1518 :
1519 : !Parameters which need to be set by hand for now...
1520 8 : weight = 1 ! The weight of the k-pts, which sum to 1.
1521 8 : tim_fourwf = 0
1522 8 : ckpt = 1 ! The kpt of the wf to FFT. We only consider gamma point for now.
1523 :
1524 : ! ndat = 1 used to be hard coded, which works fine for FFT only parallelism.
1525 8 : ndat = 1 ! # of FFT to do in parallel in fourwf
1526 :
1527 8 : eshift = 0.0 ! For PAW, opt_sij=-1 gives <G|H-lambda.S|C> in ghc instead of <G|H|C>
1528 8 : sij_opt = 0 ! Option in PAW to tell getghc what to compute...
1529 8 : tim_getghc = 0 ! timing code of the calling subroutine(can be set to 0 if not attributed)
1530 8 : type_calc = 0 ! option governing which part of Hamitonian is to be applied: 0: whole Hamiltonian
1531 8 : nband = dtset%mband
1532 8 : ispden = 1 !When required, the spin index to be used. We don't support spin polarized calculations yet...
1533 :
1534 :
1535 :
1536 : !========================================================================================================================
1537 : ! Trying to implement band parallelism (Bruno, 14/10/2013)
1538 : !------------------------------------------------------------------------------------------------------------------------
1539 : !
1540 : ! Band parallelism is described in the article on the LOBPCG method by F. Bottin et al,
1541 : ! "Large scale ab initio calculations based on three levels of parallelisation", Computational Materials Science 42 (2008) 329-336
1542 : !
1543 : ! The article describes the ABINIT implementation specifically, suggesting that the information there is directly related
1544 : ! to the code.
1545 : !
1546 : ! The main points in order to understand band + FFT parallelism are:
1547 : !
1548 : ! 1) The processors are arranged in a 2D Cartesian MPI topology, or dimensions M x P, with
1549 : ! M = nproc_bands
1550 : ! P = nproc_fft
1551 : ! There are of course M x P processors
1552 : !
1553 : ! 2) the wavefunction coefficients (namely the array cg) are distributed on these processors. Unfortunately,
1554 : ! two different distribution schemes are necessary: one for linear algebra (matrix products, etc), and one
1555 : ! consistent with the Goedecker FFTs. The code must thus change the data distribution back and forth.
1556 : !
1557 : ! 3) the "linear algebra" distribution describes the array cg directly. The G vectors are distributed among all
1558 : ! M x P processors, such that npw = npw_tot / (nproc_bands x nproc_fft). Every processor has information about
1559 : ! all bands. Schematically:
1560 : !
1561 : !
1562 : ! <- P = nproc_fft ->
1563 : ! ----------------------------------
1564 : ! | n | n | n | n | ( for all band indices n)
1565 : ! ^ | G_11 | ... | ... | G_P1 |
1566 : ! | ----------------------------------
1567 : ! M = nproc_bands | n | n | n | n |
1568 : ! | | G_1. | ... | ... | G_P. |
1569 : ! v ----------------------------------
1570 : ! | n | n | n | n |
1571 : ! | G_1M | ... | ... | G_PM |
1572 : ! ----------------------------------
1573 : !
1574 : !
1575 : ! the LOBPCG algorithm acts on blocks of bands. Thus, although each processor has information about all bands,
1576 : ! we must conceptually view the bands grouped in blocks of size M, as {c_{1}(G),...,c_{M}(G)},
1577 : ! {c_{M+1}(G)...,c_{2M}(G)}, ...
1578 : !
1579 : ! With this conceptual grouping, for a given block each processor has M x NG/(M x P) coefficients, where NG is the
1580 : ! total number of G-vectors.
1581 : !
1582 : ! 4) The distribution of information above is not appropriate for parallel FFTs. The data for one block is thus
1583 : ! transposed and takes the form
1584 : !
1585 : ! <- P = nproc_fft ->
1586 : ! ----------------------------------
1587 : ! | 1 | 1 | ... | 1 |
1588 : ! ^ | G_1 | G_2 | ... | G_P |
1589 : ! | ----------------------------------
1590 : ! M = nproc_bands | .. | ... | ... | ... |
1591 : ! | | G_1 | G_2 | ... | G_P |
1592 : ! v ----------------------------------
1593 : ! | M | M | ... | M |
1594 : ! | G_1 | G_2 | ... | G_P |
1595 : ! ----------------------------------
1596 : !
1597 : ! where the set of G vectors G_i = (G_{i1}, G_{i2}, ..., G_{iM}). Thus, a given processor has information about a single
1598 : ! band, but more G vectors. Each processor has
1599 : ! NG/(M x P) x M = NG/P coefficients, just like before.
1600 : ! Also, it is clear that the information is only communicated in the columns of the diagram above, avoiding
1601 : ! "all processors to all processors" communication.
1602 : !
1603 : ! The data is now distributed properly to do parallel FFTs! Each row in the diagram above corresponds to FFTs done
1604 : ! in parallel over nproc_fft processors, on a given band. There are M rows running thus in parallel!
1605 : !
1606 : ! The underlying ABINIT routines are equipped to handle FFT parallelism, not band distributed parallelism.
1607 : ! prep_getghc.F90 and lobpgcwf.F90 show how to rearange information to be able to apply basic ABINIT routines.
1608 : !
1609 : ! The code below is inspired / guessed from lobpcgwf and prep_getghc. WE ASSUME THERE IS NO SPINORS! CODE SHOULD
1610 : ! BE HEAVILY REVIEWED for k-points and spinors, if and when they come to be useful.
1611 : !
1612 : !========================================================================================================================
1613 :
1614 : ! variables for band parallelism, following the language in lobpcgwf
1615 8 : nspinor = 1 ! how many spinors are present. Hardcoded to 1 for now.
1616 8 : blocksize = mpi_enreg%nproc_band*mpi_enreg%bandpp ! how many bands treated in a block; # of FFT done in parallel in getghc
1617 8 : nbdblock = nband/blocksize ! how many blocks of bands are there
1618 8 : ikpt_this_proc = 1 ! Assuming only 1 kpoint, for molecules.
1619 : ! This will have to be reviewed for crystals!
1620 :
1621 8 : npw_kb = npw_k*blocksize
1622 8 : npw_g = gs_hamk2%npw_fft_k
1623 :
1624 8 : if (blocksize>1) then
1625 4 : kg_k_gather => bandfft_kpt(ikpt_this_proc)%kg_k_gather
1626 4 : ffnl_gather => bandfft_kpt(ikpt_this_proc)%ffnl_gather
1627 4 : ph3d_gather => bandfft_kpt(ikpt_this_proc)%ph3d_gather
1628 4 : kinpw_gather => bandfft_kpt(ikpt_this_proc)%kinpw_gather
1629 : else
1630 6232 : ffnl = gs_hamk2%ffnl_k
1631 15448 : ph3d = gs_hamk2%ph3d_k
1632 4 : kg_k_gather => kg_k
1633 4 : ffnl_gather => ffnl
1634 4 : ph3d_gather => ph3d
1635 4 : kinpw_gather => kinpw
1636 : endif
1637 8 : call gs_hamk%load_k(kinpw_k=kinpw_gather,kg_k=kg_k_gather,ffnl_k=ffnl_gather,ph3d_k=ph3d_gather)
1638 :
1639 736 : gbound = gs_hamk%gbound_k
1640 :
1641 : !Set up wfk routines
1642 8 : call set_wf(ucvol/nfft,n4,n5,n6,npw_k,blocksize,npw_g,mpi_enreg%comm_bandfft,mpi_enreg%comm_fft,mpi_enreg%comm_band)
1643 :
1644 : ! prepare the valence wavefunctions and projection operator to work in band+FFT parallel
1645 8 : call DistributeValenceWavefunctions()
1646 8 : call DistributeValenceKernel()
1647 :
1648 8 : cplx = 2 ! wavefunctions have complex coefficients
1649 :
1650 : ! Initialize the dummy variables
1651 8 : ABI_MALLOC(conjgrprj,(0,0))
1652 8 : ABI_MALLOC(dummy2,(0,0))
1653 8 : ABI_MALLOC(dummy3,(0,0,0))
1654 :
1655 : !Initialisation of the total counter for iterations of SQMR :
1656 8 : ktot = 0
1657 :
1658 : !Allocations of working arrays for the module
1659 : !- for pc_k function
1660 24 : ABI_MALLOC(scprod2,(2,nband))
1661 :
1662 : !- for precondition function (and set_precondition subroutine)
1663 24 : ABI_MALLOC(pcon,(npw_g))
1664 1550 : pcon = one
1665 :
1666 : !- for (private) working wf
1667 24 : ABI_MALLOC(psik1,(2,npw_k))
1668 16 : ABI_MALLOC(psik2,(2,npw_k))
1669 16 : ABI_MALLOC(psik3,(2,npw_k))
1670 16 : ABI_MALLOC(psik4,(2,npw_k))
1671 24 : ABI_MALLOC(psikb1,(2,npw_kb))
1672 16 : ABI_MALLOC(psikb2,(2,npw_kb))
1673 16 : ABI_MALLOC(psikb3,(2,npw_kb))
1674 16 : ABI_MALLOC(psikb4,(2,npw_kb))
1675 24 : ABI_MALLOC(psig1,(2,npw_g))
1676 16 : ABI_MALLOC(psig2,(2,npw_g))
1677 16 : ABI_MALLOC(psig3,(2,npw_g))
1678 16 : ABI_MALLOC(psig4,(2,npw_g))
1679 40 : ABI_MALLOC(psir1,(2,n4,n5,n6))
1680 32 : ABI_MALLOC(psir2,(2,n4,n5,n6))
1681 32 : ABI_MALLOC(psir3,(2,n4,n5,n6))
1682 24 : ABI_MALLOC(psidg,(2,nfft))
1683 40 : ABI_MALLOC(denpot,(2*n4,n5,n6))
1684 :
1685 158840 : psir1 = zero
1686 158840 : psir2 = zero
1687 158840 : psir3 = zero
1688 106856 : denpot = zero
1689 :
1690 : !Construct the vector of eigenvalues and write then to std output
1691 24 : ABI_MALLOC(eig,(nband))
1692 88 : eig=zero
1693 :
1694 8 : write(std_out,*) ch10,"Eigenvalues computation check, routine build_H:",ch10
1695 8 : io_unit_debug = get_unit()
1696 8 : write(filename_debug,'(A,I0.4,A)') "DEBUG_PROC=",mpi_enreg%me,".dat"
1697 :
1698 :
1699 :
1700 8 : mpi_communicator = mpi_enreg%comm_bandfft
1701 :
1702 8 : open(file=filename_debug,status=files_status_old,unit=io_unit_debug)
1703 :
1704 8 : write(io_unit_debug,'(A)') " Parameters:"
1705 8 : write(io_unit_debug,'(A,I5)') " nband = ", nband
1706 8 : write(io_unit_debug,'(A,I5)') " blocksize = ", blocksize
1707 8 : write(io_unit_debug,'(A,I5)') " npw_k = ", npw_k
1708 8 : write(io_unit_debug,'(A,I5)') " nbdblock = ", nbdblock
1709 :
1710 : ! temporary wavefunction array, for data in the "linear algebra" distribution
1711 24 : ABI_MALLOC( wfk_tmp1, (2,npw_k))
1712 16 : ABI_MALLOC( wfk_tmp2, (2,npw_k))
1713 :
1714 88 : do n=1, nband
1715 : ! Extract i^t/h wavefunction
1716 38630 : wfk_tmp1(:,1:npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
1717 :
1718 : ! Are the wavefunctions normalized?
1719 : !tmpc = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp1)
1720 : !call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors
1721 :
1722 : ! DEBUGGING CODE
1723 80 : write(std_out,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1)
1724 80 : write(io_unit_debug,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1)
1725 88 : flush(io_unit_debug)
1726 : end do
1727 :
1728 :
1729 : ! loop on blocks of bands. All bands in a given block will be done in parallel!
1730 :
1731 : ! loop on block of bands
1732 68 : do iblock = 1, nbdblock
1733 :
1734 : ! loop on states for this block
1735 140 : do iband = 1, blocksize
1736 :
1737 80 : n = (iblock-1)*blocksize+iband
1738 :
1739 38690 : psikb1(:,(iband-1)*npw_k+1:iband*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
1740 :
1741 : end do
1742 :
1743 : ! change configuration of the data
1744 60 : call wf_block_distribute(psikb1, psig1,1) ! LA -> FFT
1745 :
1746 : ! Apply hamiltonian on wavefunction. In bandFFT parallel, Hpsik calls prep_getghc, which knows how to transform the
1747 : ! distribution of the data from "linear algebra"-like to "FFT"-like. The output is *probably* in the "linear algebra"-like
1748 : ! distribution.
1749 60 : call Hpsik(psig2,psig1)
1750 :
1751 : ! change configuration of the data
1752 60 : call wf_block_distribute(psikb2, psig2,2) ! LA -> FFT
1753 :
1754 : ! extract the wavefunctions band by band
1755 148 : do iband=1, blocksize
1756 :
1757 80 : n = blocksize*(iblock-1) + iband ! band index
1758 :
1759 38630 : wfk_tmp1(:,1:npw_k) = psikb1(:,(iband-1)*npw_k+1:iband*npw_k)
1760 38630 : wfk_tmp2(:,1:npw_k) = psikb2(:,(iband-1)*npw_k+1:iband*npw_k)
1761 :
1762 80 : tmpc = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp2)
1763 :
1764 80 : call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors
1765 80 : eig(n) = tmpc(1)
1766 :
1767 80 : write(std_out,'(A,I5,A,F24.16,A)') "(build_H) band ", n, ", eig =", eig(n), " Ha."
1768 :
1769 : ! DEBUGGING CODE
1770 80 : write(io_unit_debug,'(A,I5,A,F24.16,A)') "band ", n, ", eig =", eig(n), " Ha."
1771 80 : flush(io_unit_debug)
1772 140 : flush(io_unit_debug)
1773 : end do
1774 :
1775 : end do
1776 :
1777 8 : ABI_FREE( wfk_tmp1 )
1778 8 : ABI_FREE( wfk_tmp2 )
1779 8 : close(io_unit_debug)
1780 :
1781 :
1782 :
1783 : ! Finishing the construction of vxc (transcribing it from the double real grid (for the density)
1784 : ! to the single real grid (for the wfs).
1785 : ! Assumes we only need one spin component; one transcription per spin being needed.
1786 8 : if(allocated(vxc_dg)) then
1787 48 : ABI_MALLOC(vxc,(n4,n5,n6,dtset%nspden))
1788 54880 : vxc = zero
1789 8 : call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vxc_dg(:,ispden),vxc(:,:,:,ispden),2)
1790 : end if
1791 :
1792 8 : timrev = 1 !Assumes istwfk *1. 1:time reversal symmetry NOT present | 2:" " " present
1793 24 : ABI_MALLOC(title,(dtset%ntypat))
1794 24 : do i=1,dtset%ntypat
1795 24 : title(i) = "Bloup" ! The clean way would be to get the psps structure in this module
1796 : ! (build_vxc is called from a place in GS calculations where it is available;
1797 : ! should be the easiest way). For now, this allows the code to run.
1798 : end do
1799 : call cryst%init(dtset%amu_orig(:,1),dtset%spgroup,dtset%natom,dtset%npsp,&
1800 : dtset%ntypat,dtset%nsym,dtset%rprimd_orig(:,:,1),dtset%typat,&
1801 : dtset%xred_orig(:,:,1),dtset%ziontypat,dtset%znucl,timrev,.false.,.false.,title,&
1802 8 : dtset%symrel,dtset%tnons,dtset%symafm)
1803 8 : ABI_FREE(title)
1804 8 : call Cryst%print()
1805 :
1806 : !TODO : Should be put in a separate build_vc constructor, and should be called right after build_H in the context of optdriver 66.
1807 8 : if(dtset%optdriver==66) then
1808 :
1809 : !Set up of the k-points and tables in the whole BZ
1810 8 : call Kmesh%init(Cryst,dtset%nkpt,dtset%kptns,dtset%kptopt,wrap_1zone=.false.)
1811 16 : call Kmesh%print([std_out], header="K-mesh for the wavefunctions")
1812 8 : call Qmesh%find_qmesh(Cryst,Kmesh)
1813 16 : call Qmesh%print([std_out], header="Q-mesh for the screening function")
1814 :
1815 : !------------------------------
1816 : !Building the vc_sqrt structure
1817 : !------------------------------
1818 : ! The sqrt_vc code relies on a different parallelism scheme than Vanilla ABINIT.
1819 : ! The Gsphere object must be built on a single processor having all the G-vectors
1820 : ! available.
1821 :
1822 : !gsph_init need gvec built according to the KSS convention, so that kg_k is not suitable (not properly sorted).
1823 8 : npw_serial=npw_k
1824 8 : call xmpi_sum(npw_serial,mpi_enreg%comm_bandfft,ierr)
1825 :
1826 8 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
1827 : call make_gvec_kss(dtset%nkpt,dtset%kptns,ecut_eff,dtset%symmorphi,dtset%nsym,dtset%symrel,dtset%tnons,Cryst%gprimd,&
1828 8 : dtset%prtvol,npw_serial,gvec,ierr)
1829 :
1830 8 : call Gsphere%init(Cryst,npw_serial,gvec=gvec)
1831 16 : call Gsphere%print([std_out], 0)
1832 :
1833 : call Vcp%init(Gsphere,Cryst,Qmesh,Kmesh,dtset%gw_rcut,dtset%gw_icutcoul,dtset%vcutgeo,dtset%ecutsigx,npw_serial,&
1834 8 : dtset%nkpt,dtset%kptns,mpi_enreg%comm_world)
1835 :
1836 : ! Since Vcp%vc_sqrt is sorted according to the KSS convention for G vectors
1837 : ! BUT will be applied to GS wavefunctions (where G vectors are sorted otherwise)
1838 : ! It is necessary to construct a vc_sqrt vector with the GS sorting.
1839 : ! Moreover, if we are in parallel (over band / FFTs), Vcp%vc_sqrt is the serial version
1840 : ! and need to be distributed according to the GS scheme (Linear Algebra configuration).
1841 24 : ABI_MALLOC(vc_sqrt,(npw_k))
1842 1293 : vc_sqrt=zero
1843 8 : k=0
1844 1301 : do i=1,npw_k
1845 331530 : do j=1,npw_serial
1846 390865 : if(all(kg_k(:,i)==gvec(:,j))) k=j
1847 : end do
1848 1293 : vc_sqrt(i)=Vcp%vc_sqrt(k,1)
1849 : end do
1850 : end if
1851 :
1852 : !--------------------------------------------------------------------------------
1853 : ! Security check : The eigenstates of the projector and the hamiltonian need to
1854 : ! agree down to the precision requested (tolwfr). Otherwise, SQMR is doomed.
1855 : ! Now that the density is read and the eigenstates are calculated in the GW run,
1856 : ! it is less useful. A test on tolwfr could be sufficient.
1857 : !
1858 : ! This remains a good tool to debug band+fft parallelism
1859 : !--------------------------------------------------------------------------------
1860 :
1861 : ! only write on the head node!
1862 8 : if (mpi_enreg%me == 0) then
1863 5 : io_unit = get_unit()
1864 5 : open(file='build_H.log',status=files_status_old,unit=io_unit)
1865 5 : write(io_unit,10) '#----------------------------------------------------------------------------'
1866 5 : write(io_unit,10) '#'
1867 5 : write(io_unit,10) '# This file presents the results of a small test to check '
1868 5 : write(io_unit,10) '# how well the Hamiltonian commutes with the projection '
1869 5 : write(io_unit,10) '# operator. '
1870 5 : write(io_unit,10) '#'
1871 5 : write(io_unit,10) '# Definitions:'
1872 5 : write(io_unit,10) '#'
1873 5 : write(io_unit,10) '# P : projections on conduction states'
1874 5 : write(io_unit,10) '# H : Hamiltonian operator'
1875 5 : write(io_unit,10) '# | psi > : eigenstate'
1876 5 : write(io_unit,10) '#'
1877 5 : flush(io_unit)
1878 : end if
1879 :
1880 4634 : psikb1 = zero
1881 :
1882 : ! sum all valence states, on copy in every band block.
1883 40 : do i=1,nbandv
1884 :
1885 88 : do mb = 1, blocksize
1886 :
1887 18584 : psikb1(:,(mb-1)*npw_k+1:mb*npw_k) = psikb1(:,(mb-1)*npw_k+1:mb*npw_k) + cg(:,(i-1)*npw_k+1:i*npw_k)
1888 :
1889 : end do
1890 :
1891 : end do
1892 :
1893 : ! normalize to one! Only sum the fist band block
1894 8 : tmpc = cg_zdotc(npw_k,psikb1,psikb1)
1895 8 : call xmpi_sum(tmpc,mpi_communicator ,ierr) ! sum on all processors
1896 4634 : psikb1 = psikb1/sqrt(tmpc(1))
1897 :
1898 :
1899 : ! change data distribution
1900 8 : call wf_block_distribute(psikb1, psig1,1) ! LA -> FFT
1901 :
1902 : ! Apply P.H operator
1903 8 : call Hpsik(psig2 ,psig1)
1904 8 : call pc_k_valence_kernel(psig2)
1905 :
1906 : ! Apply H.P operator
1907 8 : call pc_k_valence_kernel(psig1)
1908 8 : call Hpsik(psig1)
1909 :
1910 : ! compute error
1911 4642 : psig3 = psig1 - psig2 ! explicitly write the difference in an array
1912 8 : tmpc = cg_zdotc(npw_g,psig3,psig3)
1913 8 : commutation_error = tmpc(1)
1914 :
1915 8 : call xmpi_sum(commutation_error , mpi_enreg%comm_fft, ierr) ! sum on all processors working on FFT!
1916 :
1917 : ! only write on the head node!
1918 8 : if (mpi_enreg%me == 0) then
1919 5 : write(io_unit,20) ' || (PH -HP) |b> ||^2 = ',commutation_error
1920 5 : write(io_unit,20) ' tolwfr = ',tolwfr
1921 : end if
1922 :
1923 8 : if(commutation_error > tolwfr) then
1924 : !write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> Decision taken exit!'
1925 :
1926 0 : if (mpi_enreg%me == 0) write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> This must be fixed!'
1927 :
1928 0 : write(std_out,20) "WARNING-build_H: The tolerance tolwfr=",tolwfr
1929 0 : write(std_out,20) " is smaller than the commutation error ||(PH-HP)|b>||^2=",commutation_error
1930 0 : write(std_out,10) " Either set tolwfr to a less stringent value in this calculation"
1931 0 : write(std_out,10) " or to a more stringent value in the wavefunction calculation."
1932 : end if
1933 :
1934 8 : if (mpi_enreg%me == 0) close(io_unit)
1935 :
1936 : 10 format(A)
1937 : 20 format(A,ES12.3)
1938 8 : end subroutine build_H
1939 : !!***
1940 :
1941 : end module m_gwls_hamiltonian
1942 : !!***
|
