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