Line data Source code
1 : !!****m* ABINIT/m_exc_build
2 : !! NAME
3 : !! m_exc_build
4 : !!
5 : !! FUNCTION
6 : !! Build the BSE Hamiltonian in the e-h representation with MPI
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
10 : !! Copyright (C) 2009-2026 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_exc_build
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_bs_defs
28 : use m_bse_io
29 : USE_MPI
30 : use m_xmpi
31 : use m_errors
32 : use m_screen
33 : use m_hdr
34 :
35 : use m_wfd, only : wfdgw_t, wave_t, WFD_STORED
36 : use defs_datatypes, only : pseudopotential_type
37 : use m_gwdefs, only : czero_gw, cone_gw, GW_TOLQ0
38 : use m_time, only : cwtime, timab
39 : use m_io_tools, only : get_unit, open_file
40 : use m_hide_blas, only : xdotc, xgemv
41 : use m_geometry, only : normv
42 : use m_crystal, only : crystal_t
43 : use m_gsphere, only : gsphere_t
44 : use m_vcoul, only : vcoul_t
45 : use m_bz_mesh, only : kmesh_t, findqg0
46 : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g
47 : use m_pawang, only : pawang_type
48 : use m_pawtab, only : pawtab_type
49 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
50 : use m_paw_sym, only : paw_symcprj_op
51 : use m_oscillators, only : rho_tw_g, sym_rhotwgq0
52 :
53 : implicit none
54 :
55 : #if defined HAVE_MPI1
56 : include 'mpif.h'
57 : #endif
58 :
59 : private
60 : !!***
61 :
62 : public :: exc_build_ham ! Calculate and write the excitonic Hamiltonian to file.
63 : !!***
64 :
65 : contains
66 : !!***
67 :
68 : !!****f* m_exc_build/exc_build_block
69 : !! NAME
70 : !! exc_build_block
71 : !!
72 : !! FUNCTION
73 : !! Calculate and write the excitonic Hamiltonian on an external binary file (Fortran file open
74 : !! in random mode) for subsequent treatment in the Bethe-Salpeter code.
75 : !!
76 : !! INPUTS
77 : !! BSp<excparam>=The parameters for the Bethe-Salpeter calculation.
78 : !! Cryst<crystal_t>=Info on the crystalline structure.
79 : !! Kmesh<kmesh_t>=The list of k-points in the BZ, IBZ and symmetry tables.
80 : !! Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
81 : !! ktabr(nfftot_osc,BSp%nkbz)=The FFT index of $(R^{-1}(r-\tau))$ where R is symmetry needed to obtains
82 : !! the k-points from the irreducible image. Used to symmetrize u_Sk where S = \transpose R^{-1}
83 : !! Gsph_x<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored).
84 : !! Gsph_c<gsphere_t>=Info on the G-sphere used to describe the correlation part.
85 : !! Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
86 : !! screen<screen_t>=Data type gathering info and data for W.
87 : !! nfftot_osc=Total Number of FFT points used for the oscillator matrix elements.
88 : !! ngfft_osc(18)=Info on the FFT algorithm used to calculate the oscillator matrix elements.
89 : !! Psps<Pseudopotential_type>=Variables related to pseudopotentials
90 : !! Pawtab(Psps%ntypat)<pawtab_type>=PAW tabulated starting data.
91 : !! Pawang<pawang_type>=PAW angular mesh and related data.
92 : !! Paw_pwff(Cryst%ntypat*Wfd%usepaw)<pawpwff_t>=Form factor used to calculate the onsite matrix
93 : !! elements of a plane wave.
94 : !! Wfd<wfdgw_t>=Handler for the wavefunctions.
95 : !! prtvol=Verbosity level.
96 : !! rhxtwg_q0
97 : !! is_resonant
98 : !! fname
99 : !! comm=MPI communicator.
100 : !!
101 : !! OUTPUT
102 : !! The excitonic Hamiltonian is saved on an external binary file (see below).
103 : !!
104 : !! NOTES
105 : !! *) Version for K_V = K_C (q=0), thus KP_V = KP_C
106 : !! *) No exchange limit: use DFT energies in case.
107 : !! *) Symmetry of H(-k-k') = H*(k k') not used.
108 : !! *) Coulomb term can be approximated as diagonal in G.
109 : !! *) Valence bands treated from lomo on.
110 : !! *) Symmetries of the sub-blocks are used to reduce the number of elements to calculate.
111 : !!
112 : !! ____________
113 : !! |_(cv)__(vc)_|
114 : !! H_exc = | R C |
115 : !! | -C* -R* |
116 : !!
117 : !! where C is symmetric and R is Hermitian provided that the QP energies are real.
118 : !!
119 : !! For nsppol=1 ==> R = diag-W+2v; C = -W+2v
120 : !! since the Hamiltonian can be diagonalized in the spin-singlet basis set thanks to
121 : !! the fact that spin triplet does not contribute to the optical limit of epsilon.
122 : !!
123 : !! For nsppol=2 ==> R = diag-W+v; C = -W+v
124 : !! Now the matrix elements depend on the spin of the transitions but only those
125 : !! transitions in which the spin of the electron and of the hole are equal contribute
126 : !! to the macroscopic dielectric function. Moreover only the exchange term can connect transitions of different spin.
127 : !! When nsppol==2 the transitions are ordered using | (cv up) | (cv dwn) | (vc up) | (vc down) |
128 : !!
129 : !! The resonant block is given by:
130 : !!
131 : !! | (v'c' up) | (v'c' dwn) |
132 : !! ----------------------------------- where v_{-+} = v_{+-}^H when the momentum of the photon is neglected.
133 : !! | [diag-W+v]++ | v+- | (vc up) Note that v_{+-} is not Hermitian due to the presence of different spins.
134 : !! R = ----------------------------------- Actually it reduces to a Hermitian matrix when the system is not spin polarized.
135 : !! | v-+ | [diag-W+v]-- | (vc dwn) but in this case one should use nsppol=1.
136 : !! ----------------------------------- As a consequence the entire matrix is calculated and stored on file.
137 : !!
138 : !! The coupling block is given by:
139 : !!
140 : !! | (c'v' up) | (c'v dwn) |
141 : !! ----------------------------------- where v_{-+} = v_{+-}^t when the momentum of the photon is neglected.
142 : !! | [-W+v]++ | v+- | (vc up) Also in this case the entire matrix v_{+-} has to be calculated
143 : !! C = ----------------------------------- and stored on file.
144 : !! | v-+ | [-W+v]-- | (vc dwn)
145 : !! -----------------------------------
146 : !!
147 : !! SOURCE
148 :
149 23 : subroutine exc_build_block(BSp,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,Wfd,screen,Hdr_bse,&
150 23 : nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff,rhxtwg_q0,is_resonant,fname)
151 :
152 : !Arguments ------------------------------------
153 : !scalars
154 : integer,intent(in) :: nfftot_osc
155 : character(len=*),intent(in) :: fname
156 : logical,intent(in) :: is_resonant
157 : type(excparam),intent(in) :: BSp
158 : type(screen_t),intent(inout) :: screen
159 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
160 : type(crystal_t),intent(in) :: Cryst
161 : type(vcoul_t),intent(in) :: Vcp
162 : type(gsphere_t),intent(in) :: Gsph_x,Gsph_c
163 : type(Pseudopotential_type),intent(in) :: Psps
164 : type(Hdr_type),intent(inout) :: Hdr_bse
165 : type(pawang_type),intent(in) :: Pawang
166 : type(wfdgw_t),target,intent(inout) :: Wfd
167 : !arrays
168 : integer,intent(in) :: ngfft_osc(18)
169 : integer,intent(in) :: ktabr(nfftot_osc,Kmesh%nbz)
170 : complex(gwp),intent(in) :: rhxtwg_q0(BSp%npweps,BSp%lomo_min:BSp%humo_max,BSp%lomo_min:BSp%humo_max,Wfd%nkibz,Wfd%nsppol)
171 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Wfd%usepaw)
172 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Wfd%usepaw)
173 :
174 : !Local variables ------------------------------
175 : !scalars
176 : integer,parameter :: map2sphere=1,ndat1=1,master=0
177 : integer(i8b) :: bsize_my_block
178 : integer :: nspinor,nsppol,ISg,mpi_err,tmp_size,ngx
179 : integer :: ik_bz,ikp_bz,col_glob,itpk_min,itpk_max
180 : integer :: dim_rtwg,bsh_unt,ncol,dump_unt,npweps
181 : #ifdef HAVE_MPI_IO
182 : integer :: amode,mpi_fh,hmat_type,offset_err,old_type
183 : integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset
184 : logical,parameter :: is_fortran_file=.TRUE.
185 : #endif
186 : integer :: neh1,neh2,ig,nblocks
187 : integer :: ik_ibz,itim_k,ikp_ibz,itim_kp,isym_k,isym_kp
188 : integer :: iq_bz,iq_ibz,isym_q,itim_q,iqbz0,rank
189 : integer :: iv,ivp,ic,icp,jj,nrows,sender,my_ncols
190 : integer :: use_padfft,prev_nrows,spin1,spin2,block
191 : integer :: ierr,nproc,my_rank,mgfft_osc,fftalga_osc,comm
192 : integer(i8b) :: tot_nels,prev_nels,prev_ncols,nels,ir,it,itp,ist,iend,my_hsize
193 : real(dp) :: faq,kx_fact,cputime,walltime,gflops
194 : complex(sp) :: http,ctemp
195 : complex(dp) :: ph_mkpt,ph_mkt,ene_t,ene_tp
196 : logical,parameter :: with_umklp=.FALSE.
197 : logical :: use_mpiio,do_coulomb_term,do_exchange_term,w_is_diagonal,isirred
198 : logical :: is_qeq0
199 : character(len=500) :: msg
200 : type(wave_t),pointer :: wave_ck, wave_ckp, wave_vk, wave_vkp
201 : !arrays
202 : integer :: bidx(2,4),g0(3),spin_ids(2,3)
203 : integer(i8b) :: nels_block(3)
204 : integer :: my_cols(2),my_rows(2),proc_end(2),proc_start(2)
205 : integer :: my_extrema(2,2),sender_extrema(2,2),my_starts(2),my_ends(2)
206 23 : integer,allocatable :: igfftg0(:),ktabr_k(:),ktabr_kp(:),id_tab(:)
207 23 : integer,allocatable :: ncols_of(:)
208 23 : integer(i8b),allocatable :: t_start(:),t_stop(:),hsize_of(:)
209 23 : integer,allocatable :: col_start(:),col_stop(:)
210 23 : integer,allocatable :: gbound(:,:)
211 : real(dp) :: kbz(3),kpbz(3),qbz(3),spinrot_k(4),spinrot_kp(4),kmkp(3),tsec(2)
212 23 : complex(dp),allocatable :: my_bsham(:),buffer(:),buffer_2d(:,:),my_kxssp(:,:),prev_col(:)
213 : !DBYG
214 23 : complex(dp),allocatable :: acoeffs(:),bcoeffs(:),ccoeffs(:) ! Coeff of W = a/q^2 + b/q + c
215 : integer :: a_unt, b_unt, c_unt
216 : complex(dp) :: aatmp, bbtmp, cctmp
217 23 : complex(gwp),allocatable :: aa_vpv(:),aa_cpc(:),aa_ctccp(:)
218 23 : complex(gwp),allocatable :: bb_vpv1(:),bb_cpc1(:),bb_ctccp1(:)
219 23 : complex(gwp),allocatable :: bb_vpv2(:),bb_cpc2(:),bb_ctccp2(:)
220 23 : complex(gwp),allocatable :: cc_vpv(:),cc_cpc(:),cc_ctccp(:)
221 23 : complex(dp),allocatable :: abuffer(:),aprev_col(:)
222 23 : complex(dp),allocatable :: bbuffer(:),bprev_col(:)
223 23 : complex(dp),allocatable :: cbuffer(:),cprev_col(:)
224 : character(len=fnlen) :: tmpfname
225 : integer :: ii
226 : !END DBYG
227 23 : complex(gwp),allocatable :: vc_sqrt_qbz(:), rhotwg1(:),rhotwg2(:),rhxtwg_vpv(:),rhxtwg_cpc(:),ctccp(:)
228 23 : complex(gwp),target,allocatable :: ur_ckp(:),ur_vkp(:),ur_vk(:),ur_ck(:)
229 23 : complex(gwp),contiguous, pointer :: ptur_ckp(:),ptur_vkp(:),ptur_vk(:),ptur_ck(:)
230 23 : type(pawcprj_type),target,allocatable :: Cp_tmp1(:,:),Cp_tmp2(:,:)
231 23 : type(pawcprj_type),target,allocatable :: Cp_tmp3(:,:),Cp_tmp4(:,:)
232 23 : type(pawcprj_type),allocatable :: Cp_ckp(:,:),Cp_vkp(:,:)
233 23 : type(pawcprj_type),allocatable :: Cp_vk(:,:),Cp_ck(:,:)
234 23 : type(pawcprj_type),pointer :: ptcp_ckp(:,:),ptcp_vkp(:,:),ptcp_vk(:,:),ptcp_ck(:,:)
235 23 : type(pawpwij_t),allocatable :: Pwij_q(:)
236 : #ifdef HAVE_MPI_IO
237 : integer(XMPI_OFFSET_KIND) :: tmp_off,my_offpad
238 23 : integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:),offset_of_block(:)
239 : #endif
240 : #ifdef DEV_MG_DEBUG_MODE
241 : integer,allocatable :: ttp_check(:,:)
242 : #endif
243 : !************************************************************************
244 :
245 23 : call timab(680,1,tsec)
246 23 : call timab(681,1,tsec)
247 :
248 : DBG_ENTER("COLL")
249 :
250 23 : ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
251 92 : ABI_CHECK(nfftot_osc==PRODUCT(ngfft_osc(1:3)),"mismatch in FFT size")
252 :
253 23 : if (Wfd%nsppol==2) then
254 1 : ABI_WARNING("nsppol==2 is still under testing")
255 : end if
256 : ! MPI variables.
257 23 : comm = Wfd%comm
258 23 : nproc = Wfd%nproc
259 23 : my_rank = Wfd%my_rank
260 :
261 : ! Basic constants.
262 23 : nspinor = Wfd%nspinor
263 23 : nsppol = Wfd%nsppol
264 23 : dim_rtwg=1; faq = one/(Cryst%ucvol*Kmesh%nbz)
265 23 : npweps = Bsp%npweps
266 :
267 : ! Prepare the FFT tables to have u(r) on the ngfft_osc mesh.
268 92 : mgfft_osc = MAXVAL(ngfft_osc(1:3))
269 23 : fftalga_osc = ngfft_osc(7)/100
270 92 : if ( ANY(ngfft_osc(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_osc)
271 :
272 69 : ABI_MALLOC(igfftg0,(npweps))
273 69 : ABI_MALLOC(ktabr_k,(nfftot_osc))
274 46 : ABI_MALLOC(ktabr_kp,(nfftot_osc))
275 46 : ABI_MALLOC(id_tab,(nfftot_osc))
276 482671 : id_tab = (/(ic, ic=1,nfftot_osc)/)
277 :
278 : ! Workspace arrays for wavefunctions and oscillator matrix elements.
279 69 : ABI_MALLOC(rhxtwg_vpv,(npweps))
280 46 : ABI_MALLOC(rhxtwg_cpc,(npweps))
281 :
282 23 : if (BSp%prep_interp) then
283 1 : call wrtout(std_out,"Preparing BSE interpolation")
284 2 : ABI_MALLOC(aa_vpv,(npweps))
285 2 : ABI_MALLOC(bb_vpv1,(npweps))
286 2 : ABI_MALLOC(bb_vpv2,(npweps))
287 2 : ABI_MALLOC(cc_vpv,(npweps))
288 2 : ABI_MALLOC(aa_cpc,(npweps))
289 2 : ABI_MALLOC(bb_cpc1,(npweps))
290 2 : ABI_MALLOC(bb_cpc2,(npweps))
291 2 : ABI_MALLOC(cc_cpc,(npweps))
292 : end if
293 :
294 69 : ABI_MALLOC(ur_ckp,(nspinor*nfftot_osc))
295 46 : ABI_MALLOC(ur_vkp,(nspinor*nfftot_osc))
296 46 : ABI_MALLOC(ur_ck ,(nspinor*nfftot_osc))
297 46 : ABI_MALLOC(ur_vk ,(nspinor*nfftot_osc))
298 :
299 23 : if (Wfd%usepaw==1) then
300 14 : ABI_MALLOC(Cp_vk,(Wfd%natom,nspinor))
301 2 : call pawcprj_alloc(Cp_vk,0,Wfd%nlmn_atm)
302 14 : ABI_MALLOC(Cp_ck,(Wfd%natom,nspinor))
303 2 : call pawcprj_alloc(Cp_ck,0,Wfd%nlmn_atm)
304 14 : ABI_MALLOC(Cp_ckp,(Wfd%natom,nspinor))
305 2 : call pawcprj_alloc(Cp_ckp,0,Wfd%nlmn_atm)
306 14 : ABI_MALLOC(Cp_vkp,(Wfd%natom,nspinor))
307 2 : call pawcprj_alloc(Cp_vkp,0,Wfd%nlmn_atm)
308 :
309 14 : ABI_MALLOC(Cp_tmp1,(Wfd%natom,nspinor))
310 2 : call pawcprj_alloc(Cp_tmp1,0,Wfd%nlmn_atm)
311 14 : ABI_MALLOC(Cp_tmp2,(Wfd%natom,nspinor))
312 2 : call pawcprj_alloc(Cp_tmp2,0,Wfd%nlmn_atm)
313 14 : ABI_MALLOC(Cp_tmp3,(Wfd%natom,nspinor))
314 2 : call pawcprj_alloc(Cp_tmp3,0,Wfd%nlmn_atm)
315 14 : ABI_MALLOC(Cp_tmp4,(Wfd%natom,nspinor))
316 2 : call pawcprj_alloc(Cp_tmp4,0,Wfd%nlmn_atm)
317 : end if
318 : !
319 : ! Identify the index of q==0
320 23 : iqbz0=0
321 1079 : do iq_bz=1,Qmesh%nbz
322 1446 : if (ALL(ABS(Qmesh%bz(:,iq_bz))<tol3)) iqbz0 = iq_bz
323 : end do
324 23 : ABI_CHECK(iqbz0/=0,"q=0 not found")
325 : !
326 : ! Treat the spin polarization.
327 69 : spin_ids(:,1) = (/1,1/)
328 69 : spin_ids(:,2) = (/2,2/)
329 69 : spin_ids(:,3) = (/1,2/)
330 :
331 23 : nblocks=1
332 23 : kx_fact=two
333 23 : nels_block(:)=0
334 23 : nels_block(1)=BSp%nreh(1)*(BSp%nreh(1)+1_i8b)/2
335 23 : tot_nels=nels_block(1)
336 :
337 23 : if (nsppol==2) then
338 1 : nblocks=3
339 1 : kx_fact=one
340 : nels_block(1) = BSp%nreh(1)*(BSp%nreh(1)+1_i8b)/2 ! Only the upper triangle for block 1 and 2
341 1 : nels_block(2) = BSp%nreh(2)*(BSp%nreh(2)+1_i8b)/2
342 1 : nels_block(3) = BSp%nreh(1)*BSp%nreh(2)*1_i8b ! Note: Block 3 does not have symmetries.
343 4 : tot_nels= SUM(nels_block)
344 : end if
345 : !
346 : ! Distribute the calculation of the matrix elements among the nodes.
347 : ! * tstart and t_stop give the initial and final transition index treated by each node.
348 : ! * my_hsize is the number of transitions treated by this processor
349 : ! * my_cols(1:2) gives the initial and final column treated by this node.
350 : !
351 : use_mpiio=.FALSE.
352 : #ifdef HAVE_MPI_IO
353 : use_mpiio = (nproc>1)
354 : #endif
355 23 : use_mpiio=.FALSE.
356 : !use_mpiio=.TRUE.
357 :
358 23 : if (is_resonant) then
359 : if (use_mpiio) then
360 : write(msg,'(2a,f6.2,a)')&
361 : ". Writing resonant excitonic Hamiltonian on file "//TRIM(fname)," via MPI-IO; file size= ",two*tot_nels*dp*b2Gb," [Gb]."
362 : else
363 : write(msg,'(2a,f6.2,a)')&
364 22 : ". Writing resonant excitonic Hamiltonian on file "//TRIM(fname),"; file size= ",two*dp*tot_nels*b2Gb," [Gb]."
365 : end if
366 : else
367 : if (use_mpiio) then
368 : write(msg,'(2a,f6.2,a)')&
369 : ". Writing coupling excitonic Hamiltonian on file "//TRIM(fname)," via MPI-IO; file size= ",tot_nels*2*dp*b2Gb," [Gb]."
370 : else
371 : write(msg,'(2a,f6.2,a)')&
372 1 : ". Writing coupling excitonic Hamiltonian on file "//TRIM(fname),"; file size= ",two*dp*tot_nels*b2Gb," [Gb]."
373 : end if
374 : end if
375 69 : call wrtout([std_out, ab_out], msg, do_flush=.True.)
376 :
377 : ! Master writes the BSE header with Fortran IO.
378 23 : if (my_rank==master) then
379 19 : if (open_file(fname,msg,newunit=bsh_unt,form="unformatted",action="write") /= 0) then
380 0 : ABI_ERROR(msg)
381 : end if
382 19 : call exc_write_bshdr(bsh_unt,Bsp,Hdr_bse)
383 : ! To force the writing (needed for MPI-IO).
384 19 : close(bsh_unt)
385 :
386 19 : if (.not.use_mpiio) then ! Reopen the file and skip the header.
387 19 : if (open_file(fname,msg,newunit=bsh_unt,form="unformatted",action="readwrite") /= 0) then
388 0 : ABI_ERROR(msg)
389 : end if
390 19 : call exc_skip_bshdr(bsh_unt,ierr)
391 : end if
392 :
393 19 : if (BSp%prep_interp) then
394 1 : tmpfname = fname
395 1 : ii = LEN_TRIM(fname)
396 1 : tmpfname(ii-2:ii+1) = 'ABSR'
397 1 : if (open_file(tmpfname,msg,newunit=a_unt,form='unformatted',action="write") /= 0) then
398 0 : ABI_ERROR(msg)
399 : end if
400 1 : tmpfname(ii-2:ii+1) = 'BBSR'
401 1 : if (open_file(tmpfname,msg,newunit=b_unt,form='unformatted',action="write") /= 0) then
402 0 : ABI_ERROR(msg)
403 : end if
404 1 : tmpfname(ii-2:ii+1) = 'CBSR'
405 1 : if (open_file(tmpfname,msg,newunit=c_unt,form='unformatted',action="write") /= 0) then
406 0 : ABI_ERROR(msg)
407 : end if
408 1 : call exc_write_bshdr(a_unt,Bsp,Hdr_bse)
409 1 : call exc_write_bshdr(b_unt,Bsp,Hdr_bse)
410 1 : call exc_write_bshdr(c_unt,Bsp,Hdr_bse)
411 1 : close(a_unt)
412 1 : close(b_unt)
413 1 : close(c_unt)
414 3 : if (.not.use_mpiio) then ! Reopen the file and skip the header.
415 1 : tmpfname(ii-2:ii+1) = 'ABSR'
416 1 : if (open_file(tmpfname,msg,newunit=a_unt,form='unformatted',action="readwrite") /= 0) then
417 0 : ABI_ERROR(msg)
418 : end if
419 1 : call exc_skip_bshdr(a_unt,ierr)
420 1 : tmpfname(ii-2:ii+1) = 'BBSR'
421 1 : if (open_file(tmpfname,msg,newunit=b_unt,form='unformatted',action="readwrite") /= 0) then
422 0 : ABI_ERROR(msg)
423 : end if
424 1 : call exc_skip_bshdr(b_unt,ierr)
425 1 : tmpfname(ii-2:ii+1) = 'CBSR'
426 1 : if (open_file(tmpfname,msg,newunit=c_unt,form='unformatted',action="readwrite") /= 0) then
427 0 : ABI_ERROR(msg)
428 : end if
429 1 : call exc_skip_bshdr(c_unt,ierr)
430 : end if
431 : end if
432 : end if
433 :
434 23 : call xmpi_barrier(comm)
435 :
436 : if (use_mpiio) then
437 : #ifdef HAVE_MPI_IO
438 : ! Open the file with MPI-IO
439 : amode = MPI_MODE_RDWR
440 : !amode = MPI_MODE_CREATE + MPI_MODE_RDWR,
441 :
442 : call MPI_FILE_OPEN(comm, fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
443 : ABI_CHECK_MPI(mpi_err,"opening: "//TRIM(fname))
444 :
445 : ! Skip the header.
446 : call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
447 :
448 : ! Precompute the offset of the each block including the Fortran markers.
449 : ABI_MALLOC(offset_of_block,(nblocks))
450 : offset_of_block(1) = ehdr_offset
451 : do block=2,nblocks
452 : tmp_off = offset_of_block(block-1) + nels_block(block-1)*xmpi_bsize_dpc
453 : tmp_off = tmp_off + Bsp%nreh(block-1)*2*xmpio_bsize_frm ! markers.
454 : offset_of_block(block) = tmp_off
455 : end do
456 : #endif
457 : end if
458 :
459 23 : call timab(681,2,tsec)
460 :
461 47 : do block=1,nsppol
462 : !
463 : ! Indices used to loop over bands.
464 : ! bidx contains the starting and final indices used to loop over bands.
465 : !
466 : ! (b3,b4)
467 : ! |... ...|
468 : ! (b1,b2) |... ...|
469 : !
470 : ! Resonant matrix is given by
471 : ! (v',c')
472 : ! |... ...|
473 : ! (v,c) |... ...|
474 : !
475 : ! Coupling matrix is given by
476 : ! (c',v')
477 : ! |... ...|
478 : ! (v,c) |... ...|
479 :
480 24 : if (is_resonant) then
481 69 : bidx(:,1) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b1
482 69 : bidx(:,2) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b2
483 69 : bidx(:,3) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b3
484 69 : bidx(:,4) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b4
485 : else
486 3 : bidx(:,1) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b1
487 3 : bidx(:,2) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b2
488 3 : bidx(:,3) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b3
489 3 : bidx(:,4) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b4
490 : end if
491 :
492 24 : spin1 = spin_ids(1,block)
493 24 : spin2 = spin_ids(2,block)
494 :
495 24 : do_coulomb_term = (Bsp%use_coulomb_term .and. (spin1==spin2))
496 24 : do_exchange_term = (Bsp%exchange_term>0)
497 24 : w_is_diagonal = BSp%use_diagonal_Wgg
498 : !
499 : ! Distribution of the matrix elements among the nodes.
500 : ! Note that rank0 will get the first transitions.
501 24 : nels=nels_block(block)
502 72 : ABI_MALLOC(t_start,(0:nproc-1))
503 48 : ABI_MALLOC(t_stop,(0:nproc-1))
504 24 : call xmpi_split_work2_i8b(nels,nproc,t_start,t_stop)
505 :
506 48 : ABI_MALLOC(hsize_of,(0:nproc-1))
507 62 : hsize_of=0
508 62 : do rank=0,nproc-1
509 62 : if (t_stop(rank)>=t_start(rank)) hsize_of(rank) = t_stop(rank)-t_start(rank)+1
510 : !write(std_out,*)"nels",nels,hsize_of(rank)
511 : end do
512 :
513 24 : my_hsize = hsize_of(my_rank)
514 24 : if (my_hsize<=0) then
515 0 : write(msg,'(a,i0)')"Wrong number of transitions: my_hsize= ",my_hsize
516 0 : ABI_ERROR(msg)
517 : end if
518 24 : if (my_hsize /= INT(my_hsize,KIND=i4b)) then
519 0 : write(msg,'(a,i0)')"Size of local block too large for a default integer, Increase the number of CPUs: my_hsize= ",my_hsize
520 0 : ABI_ERROR(msg)
521 : end if
522 :
523 24 : my_cols=0
524 9272 : do itp=1,Bsp%nreh(block)
525 2182216 : do it=1,itp
526 2172944 : ir = it + itp*(itp-1_i8b)/2
527 2172944 : if (ir==t_start(my_rank)) then
528 24 : my_rows(1) = it
529 24 : my_cols(1) = itp
530 : end if
531 2182192 : if (ir==t_stop(my_rank)) then
532 24 : my_rows(2) = it
533 24 : my_cols(2) = itp
534 : end if
535 : end do
536 : end do
537 :
538 72 : my_starts = [my_rows(1),my_cols(1)]
539 24 : my_ends = [my_rows(2),my_cols(2)]
540 : !
541 : ! Announce the treatment of submatrix treated by each node.
542 24 : bsize_my_block = 2*dp*my_hsize
543 24 : write(msg,'(4(a,i0))')' Treating ',my_hsize,'/',nels,' matrix elements, from column ',my_cols(1),' up to column ',my_cols(2)
544 24 : call wrtout(std_out, msg)
545 :
546 24 : if (is_resonant) then
547 23 : write(msg,'(a,f8.1,a)')' Calculating resonant blocks. Memory required: ',bsize_my_block*b2Mb,' [Mb] <<< MEM'
548 : else
549 1 : write(msg,'(a,f8.1,a)')' Calculating coupling blocks. Memory required: ',bsize_my_block*b2Mb,' [Mb] <<< MEM'
550 : end if
551 24 : call wrtout(std_out, msg)
552 :
553 : ! Allocate big (scalable) buffer to store the BS matrix on this node.
554 72 : ABI_MALLOC_OR_DIE(my_bsham,(t_start(my_rank):t_stop(my_rank)), ierr)
555 :
556 24 : if (BSp%prep_interp) then
557 : ! Allocate big (scalable) buffers to store a,b,c coefficients
558 2 : ABI_MALLOC_OR_DIE(acoeffs,(t_start (my_rank):t_stop(my_rank)), ierr)
559 2 : ABI_MALLOC_OR_DIE(bcoeffs,(t_start(my_rank):t_stop(my_rank)), ierr)
560 2 : ABI_MALLOC_OR_DIE(ccoeffs,(t_start(my_rank):t_stop(my_rank)), ierr)
561 : end if
562 :
563 24 : if (do_coulomb_term) then ! Construct Coulomb term.
564 :
565 24 : call timab(682,1,tsec) ! exc_build_ham(Coulomb)
566 :
567 24 : write(msg,'(a,2i2,a)')" Calculating direct Coulomb term for (spin1,spin2) ",spin1,spin2," using full W_{GG'} ..."
568 24 : if (w_is_diagonal) then
569 : write(msg,'(a,2i2,a)')&
570 2 : " Calculating direct Coulomb term for (spin1, spin2) ",spin1,spin2," using diagonal approximation for W_{GG'} ..."
571 : end if
572 24 : call wrtout(std_out, msg)
573 :
574 72 : ABI_MALLOC(ctccp,(npweps))
575 :
576 24 : if (BSp%prep_interp) then
577 2 : ABI_MALLOC(aa_ctccp,(npweps))
578 2 : ABI_MALLOC(bb_ctccp1,(npweps))
579 2 : ABI_MALLOC(bb_ctccp2,(npweps))
580 2 : ABI_MALLOC(cc_ctccp,(npweps))
581 : end if
582 :
583 48 : ABI_MALLOC(vc_sqrt_qbz,(npweps))
584 :
585 : #ifdef DEV_MG_DEBUG_MODE
586 : ABI_MALLOC(ttp_check,(BSp%nreh(block),BSp%nreh(block)))
587 : ttp_check=0
588 : #endif
589 :
590 1112 : do ikp_bz=1,BSp%nkbz ! Loop over kp
591 : ! NOTE: this way of looping is good for bulk but it's not optimal in the
592 : ! case of systems sampled only at Gamma e.g. isolated systems in which
593 : ! one should take advantage of Hermiticity by looping over c-v !!!!
594 :
595 : ! Check whether (vp,cp,ikp_bz,spin2) belongs to the set of columns treated by me for some vp,cp
596 : ! Be careful since vcks2t contains zeros corresponding to transitions that should be skipped.
597 13456 : itpk_min = MINVAL(Bsp%vcks2t(:,:,ikp_bz,spin2), MASK=(Bsp%vcks2t(:,:,ikp_bz,spin2)>0) )
598 13456 : itpk_max = MAXVAL(Bsp%vcks2t(:,:,ikp_bz,spin2))
599 1088 : if (my_cols(2)<itpk_min .or. my_cols(1)>itpk_max) CYCLE
600 :
601 836 : write(msg,'(3(a,i0))')" status: ",ikp_bz,"/",BSp%nkbz," done by node ",my_rank
602 836 : call wrtout(std_out, msg, do_flush=.True.)
603 :
604 : ! * Get ikp_ibz, non-symmorphic phase, ph_mkpt, and symmetries from ikp_bz.
605 836 : call kmesh%get_BZ_item(ikp_bz,kpbz,ikp_ibz,isym_kp,itim_kp,ph_mkpt,isirred=isirred)
606 : !ABI_CHECK(isirred,"not irred!")
607 : !ABI_CHECK(ph_mkpt == cone, "Wrong phase!")
608 :
609 6059836 : ktabr_kp(:) = ktabr(:,ikp_bz)
610 4180 : spinrot_kp(:)=Cryst%spinrot(:,isym_kp)
611 : !ABI_CHECK(ALL(ktabr_kp == id_tab), "wrong tab")
612 :
613 23600 : do ik_bz=1,ikp_bz ! Loop over k
614 : !
615 : ! * Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz
616 22740 : call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt,isirred=isirred)
617 : !ABI_CHECK(isirred,"not irred!")
618 : !ABI_CHECK(ph_mkt == cone, "Wrong phase!")
619 :
620 172008240 : ktabr_k(:) = ktabr(:,ik_bz)
621 113700 : spinrot_k(:)=Cryst%spinrot(:,isym_k)
622 : !ABI_CHECK(ALL(ktabr_k == id_tab), "wrong tab")
623 : !if(itim_k==2) CYCLE ! time-reversal or not
624 : !
625 : ! * Find q = K-KP-G0 in the full BZ.
626 90960 : kmkp = Kmesh%bz(:,ik_bz) - Kmesh%bz(:,ikp_bz)
627 22740 : call findqg0(iq_bz,g0,kmkp,Qmesh%nbz,Qmesh%bz,BSp%mG0)
628 :
629 : ! Evaluate the tables needed for the padded FFT performed in rhotwg. Note that we have
630 : ! to pass G-G0 to sphereboundary instead of G as we need FFT results on the shifted G-sphere,
631 : ! If Gamma is not inside G-G0 one has to disable FFT padding as sphereboundary will give wrong tables.
632 : ! * Get the G-G0 shift for the FFT of the oscillators.
633 : !
634 90960 : ABI_MALLOC(gbound,(2*mgfft_osc+8,2))
635 22740 : call Gsph_c%fft_tabs(g0,mgfft_osc,ngfft_osc,use_padfft,gbound,igfftg0)
636 22740 : if ( ANY(fftalga_osc == (/2,4/)) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
637 22740 : if (use_padfft==0) then
638 28 : ABI_FREE(gbound)
639 56 : ABI_MALLOC(gbound,(2*mgfft_osc+8,2*use_padfft))
640 : end if
641 : !
642 : ! Get iq_ibz, and symmetries from iq_bz
643 22740 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
644 22740 : is_qeq0 = (normv(qbz,Cryst%gmet,'G')<GW_TOLQ0)
645 :
646 : ! Symmetrize em1(omega=0)
647 22740 : call screen%rotate_iqbz(iq_bz, Cryst, Gsph_c, Qmesh, Vcp)
648 :
649 : ! Set up table of |q_BZ+G|
650 22740 : if (iq_ibz==1) then
651 40208 : do ig=1,npweps
652 39372 : isg = Gsph_c%rottb(ig,itim_q,isym_q)
653 40208 : vc_sqrt_qbz(isg)=Vcp%vcqlwl_sqrt(ig,1)
654 : end do
655 : else
656 1014192 : do ig=1,npweps
657 992288 : isg = Gsph_c%rottb(ig,itim_q,isym_q)
658 1014192 : vc_sqrt_qbz(isg) = Vcp%vc_sqrt(ig,iq_ibz)
659 : end do
660 : end if
661 :
662 : ! === Evaluate oscillator matrix elements ===
663 : ! * $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
664 22740 : if (Wfd%usepaw==1.and.ik_bz/=ikp_bz) then
665 4108 : ABI_MALLOC(Pwij_q,(Cryst%ntypat))
666 1020 : call pawpwij_init(Pwij_q,npweps,Qmesh%bz(:,iq_bz),Gsph_c%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
667 : end if
668 :
669 : ! =======================================
670 : ! === Loop over the four band indices ===
671 : ! =======================================
672 90388 : do ic=bidx(1,2),bidx(2,2) !do ic=BSp%lumo,BSp%nbnds
673 :
674 67648 : ABI_CHECK(wfd%get_wave_ptr(ic, ik_ibz, spin1, wave_ck, msg) == 0, msg)
675 67648 : if (wave_ck%has_ur == WFD_STORED) then
676 67420 : ptur_ck => wave_ck%ur
677 : else
678 228 : call wfd%get_ur(ic,ik_ibz,spin1,ur_ck)
679 228 : ptur_ck => ur_ck
680 : end if
681 : !
682 : ! Get cprj for this (c,kbz,s1) in the BZ.
683 : ! phase due to the umklapp G0 in k-q is already included.
684 67648 : if (Wfd%usepaw==1) then
685 4296 : if (wave_ck%has_cprj == WFD_STORED) then
686 4296 : ptcp_ck => wave_ck%cprj
687 : else
688 0 : call wfd%get_cprj(ic,ik_ibz,spin1,Cryst,Cp_tmp1,sorted=.FALSE.)
689 0 : ptcp_ck => Cp_tmp1
690 : end if
691 4296 : call paw_symcprj_op(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,ptcp_ck,Cp_ck)
692 : end if
693 :
694 306004 : do icp=bidx(1,4),bidx(2,4) !do icp=BSp%lumo,BSp%nbnds
695 : ! Calculate matrix-elements rhxtwg_cpc
696 : !
697 215616 : if (ik_bz==ikp_bz) then ! Already in memory.
698 8688 : rhxtwg_cpc(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,icp,ic,ik_ibz,spin1),Gsph_c)
699 :
700 : else
701 : ! Calculate matrix element from wfr.
702 : ! TODO: change the order of the loops.
703 :
704 206928 : ABI_CHECK(wfd%get_wave_ptr(icp, ikp_ibz, spin2, wave_ckp, msg) == 0, msg)
705 206928 : if (wave_ckp%has_ur == WFD_STORED) then
706 206724 : ptur_ckp => wave_ckp%ur
707 : else
708 204 : call wfd%get_ur(icp,ikp_ibz,spin2,ur_ckp)
709 204 : ptur_ckp => ur_ckp
710 : end if
711 :
712 : ! Load cprj for this (c,k,s2) in the BZ.
713 : ! Do not care about umklapp G0 in k-q as the phase is already included.
714 206928 : if (Wfd%usepaw==1) then
715 15984 : if (wave_ckp%has_cprj == WFD_STORED) then
716 15984 : ptcp_ckp => wave_ckp%cprj
717 : else
718 0 : call wfd%get_cprj(icp,ikp_ibz,spin2,Cryst,Cp_tmp2,sorted=.FALSE.)
719 0 : ptcp_ckp => Cp_tmp2
720 : end if
721 15984 : call paw_symcprj_op(ikp_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,ptcp_ckp,Cp_ckp)
722 : end if
723 :
724 : call rho_tw_g(nspinor,npweps,nfftot_osc,ndat1,ngfft_osc,map2sphere,use_padfft,igfftg0,gbound,&
725 : ptur_ckp,itim_kp,ktabr_kp,ph_mkpt,spinrot_kp,&
726 : ptur_ck ,itim_k ,ktabr_k ,ph_mkt ,spinrot_k ,&
727 206928 : dim_rtwg,rhxtwg_cpc)
728 :
729 206928 : if (Wfd%usepaw==1) then ! Add PAW onsite contribution.
730 15984 : call paw_rho_tw_g(cryst,Pwij_q,npweps,dim_rtwg,nspinor,Gsph_c%gvec, Cp_ckp,Cp_ck,rhxtwg_cpc)
731 : end if
732 : end if
733 :
734 215616 : if (BSp%prep_interp) then
735 30528 : aa_cpc = rhxtwg_cpc
736 29376 : aa_cpc(2:) = czero
737 30528 : bb_cpc1 = vc_sqrt_qbz*rhxtwg_cpc
738 576 : bb_cpc1(1) = czero
739 30528 : bb_cpc2 = rhxtwg_cpc
740 29376 : bb_cpc2(2:) = czero
741 :
742 576 : if(ik_bz == ikp_bz) then
743 : ! Enforce orthogonality of the wavefunctions.
744 128 : if(icp == ic) then
745 32 : aa_cpc(1) = cone
746 32 : bb_cpc2(1) = cone
747 : else
748 96 : aa_cpc(1) = czero
749 96 : bb_cpc2(1) = czero
750 : end if
751 : end if
752 :
753 : ! MG TODO: a does not require a call to w0gemv
754 576 : call screen%w0gemv("C",npweps,nspinor,w_is_diagonal,cone_gw,czero_gw,aa_cpc,aa_ctccp)
755 576 : call screen%w0gemv("C",npweps,nspinor,w_is_diagonal,cone_gw,czero_gw,bb_cpc1,bb_ctccp1)
756 576 : call screen%w0gemv("C",npweps,nspinor,w_is_diagonal,cone_gw,czero_gw,bb_cpc2,bb_ctccp2)
757 :
758 30528 : cc_cpc = vc_sqrt_qbz*rhxtwg_cpc
759 576 : cc_cpc(1) = czero
760 :
761 576 : call screen%w0gemv("C",npweps,nspinor,w_is_diagonal,cone_gw,czero_gw,cc_cpc,cc_ctccp)
762 : end if
763 :
764 : ! Prepare sum_GG' rho_c'c*(G) W_qbz(G,G') rho_v'v(G')
765 : ! First sum on G: sum_G rho_c'c(G) W_qbz*(G,G') (W_qbz conjugated)
766 9953856 : rhxtwg_cpc = rhxtwg_cpc * vc_sqrt_qbz
767 215616 : call screen%w0gemv("C",npweps,nspinor,w_is_diagonal,cone_gw,czero_gw,rhxtwg_cpc,ctccp)
768 :
769 921936 : do iv=bidx(1,1),bidx(2,1) !do iv=BSp%lomo,BSp%homo
770 638672 : it = BSp%vcks2t(iv,ic,ik_bz,spin1); if (it==0) CYCLE ! ir-uv-cutoff
771 638672 : ene_t = BSp%Trans(it,spin1)%en
772 :
773 : ! TODO: use this but change the order of the loops.
774 638672 : ABI_CHECK(wfd%get_wave_ptr(iv, ik_ibz, spin1, wave_vk, msg) == 0, msg)
775 :
776 638672 : if (wave_vk%has_ur == WFD_STORED) then
777 638453 : ptur_vk => wave_vk%ur
778 : else
779 219 : call wfd%get_ur(iv,ik_ibz,spin1,ur_vk)
780 219 : ptur_vk => ur_vk
781 : end if
782 : !
783 : ! Load cprj for this (v,k,s1) in the BZ.
784 : ! Do not care about umklapp G0 in k-q as the phase is already included.
785 638672 : if (Wfd%usepaw==1) then
786 51264 : if (wave_vk%has_cprj == WFD_STORED) then
787 51264 : ptcp_vk => wave_vk%cprj
788 : else
789 0 : call wfd%get_cprj(iv,ik_ibz,spin1,Cryst,Cp_tmp3,sorted=.FALSE.)
790 0 : ptcp_vk => Cp_tmp3
791 : end if
792 51264 : call paw_symcprj_op(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,ptcp_vk,Cp_vk)
793 : end if
794 :
795 2773248 : do ivp=bidx(1,3),bidx(2,3) !do ivp=BSp%lomo,BSp%homo
796 :
797 1918960 : if (is_resonant) then
798 1842928 : itp = BSp%vcks2t(ivp,icp,ikp_bz,spin2)
799 : else ! have to exchange band indices
800 76032 : itp = BSp%vcks2t(icp,ivp,ikp_bz,spin2)
801 : end if
802 :
803 1918960 : if (itp==0) CYCLE ! ir-uv-cutoff
804 :
805 : ! FIXME Temporary work around, when ikp_bz == ik it might happen that itp<it
806 : ! should rewrite the loops using contracted k-dependent indices for bands
807 1918960 : if (itp<it) CYCLE
808 :
809 1883684 : ir = it + itp*(itp-1)/2
810 1883684 : if (ir<t_start(my_rank).or.ir>t_stop(my_rank)) CYCLE
811 :
812 1877264 : ene_tp = BSp%Trans(itp,spin2)%en
813 :
814 : ! ============================================
815 : ! === Calculate matrix elements rhxtwg_vpv ===
816 : ! ============================================
817 1877264 : if (ik_bz==ikp_bz) then
818 : ! Already in memory.
819 42928 : rhxtwg_vpv(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,ivp,iv,ik_ibz,spin1),Gsph_c)
820 :
821 : else
822 :
823 1834336 : ABI_CHECK(wfd%get_wave_ptr(ivp, ikp_ibz, spin2, wave_vkp, msg) == 0, msg)
824 :
825 : ! Calculate matrix element from wfr.
826 1834336 : if (wave_vkp%has_ur == WFD_STORED) then
827 1834123 : ptur_vkp => wave_vkp%ur
828 : else
829 213 : call wfd%get_ur(ivp,ikp_ibz,spin2,ur_vkp)
830 213 : ptur_vkp => ur_vkp
831 : end if
832 : !
833 : ! Load cprj for this (vp,kp,s2) in the BZ.
834 : ! Do not care about umklapp G0 in k-q as the phase is already included.
835 1834336 : if (Wfd%usepaw==1) then
836 144640 : if (wave_vkp%has_cprj == WFD_STORED) then
837 144640 : ptcp_vkp => wave_vkp%cprj
838 : else
839 0 : call wfd%get_cprj(ivp,ikp_ibz,spin2,Cryst,Cp_tmp4,sorted=.FALSE.)
840 0 : ptcp_vkp => Cp_tmp4
841 : end if
842 144640 : call paw_symcprj_op(ikp_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,ptcp_vkp,Cp_vkp)
843 : end if
844 :
845 : call rho_tw_g(nspinor,npweps,nfftot_osc,ndat1,ngfft_osc,map2sphere,use_padfft,igfftg0,gbound,&
846 : ptur_vkp,itim_kp,ktabr_kp,ph_mkpt,spinrot_kp,&
847 : ptur_vk ,itim_k ,ktabr_k ,ph_mkt ,spinrot_k ,&
848 1834336 : dim_rtwg,rhxtwg_vpv)
849 :
850 1834336 : if (Wfd%usepaw==1) then ! Add PAW onsite contribution.
851 144640 : call paw_rho_tw_g(cryst,Pwij_q,npweps,dim_rtwg,nspinor,Gsph_c%gvec,Cp_vkp,Cp_vk,rhxtwg_vpv)
852 : end if
853 : end if
854 :
855 : ! Index in the global Hamiltonian matrix.
856 : ir = it + itp*(itp-1_i8b)/2
857 :
858 1877264 : if (ir<t_start(my_rank).or.ir>t_stop(my_rank)) then
859 0 : write(msg,'(a,3(1x,i0))')" Gonna SIGFAULT, ir, t_start, t_stop ",ir,t_start(my_rank),t_stop(my_rank)
860 0 : ABI_ERROR(msg)
861 : end if
862 : !ABI_CHECK(itp >= it,"itp < it")
863 :
864 1877264 : if (BSp%prep_interp) then
865 : ! Save a,b, c coefficients.
866 246768 : aa_vpv = rhxtwg_vpv
867 237456 : aa_vpv(2:) = czero
868 246768 : bb_vpv1 = rhxtwg_vpv
869 237456 : bb_vpv1(2:) = czero
870 246768 : bb_vpv2 = vc_sqrt_qbz*rhxtwg_vpv
871 4656 : bb_vpv2(1) = czero
872 :
873 4656 : if (ik_bz == ikp_bz) then
874 : ! Enforce orthogonality of the wavefunctions.
875 624 : if (ivp == iv) then
876 240 : aa_vpv(1) = cone
877 240 : bb_vpv1(1) = cone
878 : else
879 384 : aa_vpv(1) = czero
880 384 : bb_vpv1(1) = czero
881 : end if
882 : end if
883 :
884 246768 : cc_vpv = vc_sqrt_qbz*rhxtwg_vpv
885 4656 : cc_vpv(1) = czero
886 :
887 4656 : aatmp = -faq * xdotc(npweps,aa_ctccp,1,aa_vpv,1)
888 4656 : bbtmp = -faq * xdotc(npweps,bb_ctccp1,1,bb_vpv1,1)-faq*xdotc(npweps,bb_ctccp2,1,bb_vpv2,1)
889 4656 : cctmp = -faq * xdotc(npweps,cc_ctccp,1,cc_vpv,1)
890 :
891 4656 : acoeffs(ir) = aatmp
892 4656 : bcoeffs(ir) = bbtmp
893 4656 : ccoeffs(ir) = cctmp
894 : end if
895 :
896 : ! sum_G2 rho_c'c(G) W_qbz(G,G') rho_v'v(G')
897 87272848 : rhxtwg_vpv = vc_sqrt_qbz * rhxtwg_vpv
898 1877264 : http = - faq * xdotc(npweps,ctccp,1,rhxtwg_vpv,1)
899 :
900 : ! Save result taking into account the symmetry of the matrix.
901 : ! Note that the diagonal of the resonant block is not forced to be real
902 2557632 : my_bsham(ir) = http
903 :
904 : #ifdef DEV_MG_DEBUG_MODE
905 : ttp_check(it,itp) = ttp_check(it,itp)+1
906 : #endif
907 : end do !ivp
908 : end do !iv
909 : end do !icp
910 : end do !ic
911 :
912 22740 : ABI_FREE(gbound)
913 :
914 23828 : if (Wfd%usepaw==1.and.ik_bz/=ikp_bz) then ! Free the onsite contribution for this q.
915 1020 : call pawpwij_free(Pwij_q)
916 2068 : ABI_FREE(Pwij_q)
917 : end if
918 :
919 : end do ! ik_bz
920 : end do ! Fat loop over ikp_bz
921 :
922 : #ifdef DEV_MG_DEBUG_MODE
923 : do itp=1,BSp%nreh(block)
924 : do it=1,BSp%nreh(block)
925 : ir = it + itp*(itp-1_i8b)/2
926 : if (itp>=it .and. ttp_check(it,itp) /= 1) then
927 : if (ir>=t_start(my_rank).and.ir<=t_stop(my_rank)) then
928 : write(std_out,*)"WARN: upper triangle is not 1 ",it,itp,ttp_check(it,itp)
929 : write(std_out,*)TRIM(repr_trans(Bsp%Trans(it ,spin1)))
930 : write(std_out,*)TRIM(repr_trans(Bsp%Trans(itp,spin2)))
931 : end if
932 : end if
933 : if (itp< it .and. ttp_check(it,itp) /= 0) then
934 : write(std_out,*)"WARN: then lower triangle is not 0 ",it,itp,ttp_check(it,itp)
935 : write(std_out,*)TRIM(repr_trans(Bsp%Trans(it ,spin1)))
936 : write(std_out,*)TRIM(repr_trans(Bsp%Trans(itp,spin2)))
937 : end if
938 : end do
939 : end do
940 : ierr = SUM(SUM(ttp_check,DIM=2),DIM=1)
941 : if (ierr/=my_hsize) then
942 : write(msg,'(a,2i0)')"ierr/=my_hsize",ierr,my_hsize
943 : ABI_ERROR(msg)
944 : end if
945 : ABI_FREE(ttp_check)
946 : #endif
947 :
948 24 : ABI_FREE(ctccp)
949 24 : if(Bsp%prep_interp) then
950 1 : ABI_FREE(aa_ctccp)
951 1 : ABI_FREE(bb_ctccp1)
952 1 : ABI_FREE(bb_ctccp2)
953 1 : ABI_FREE(cc_ctccp)
954 : end if
955 :
956 24 : ABI_FREE(vc_sqrt_qbz)
957 24 : call wrtout(std_out,' Coulomb term completed')
958 :
959 24 : call timab(682,2,tsec) ! exc_build_ham(Coulomb)
960 : end if ! do_coulomb_term
961 : !
962 : ! =====================
963 : ! === Exchange term ===
964 : ! =====================
965 : ! TODO might add treatment of <psi|q+G|psi> for q+G -> 0
966 : ! TODO might used enlarged G-sphere for better convergence.
967 24 : if (do_exchange_term) then
968 :
969 : !call exc_build_v(spin1,spin2,nsppol,npweps,Bsp,Cryst,Kmesh,Qmesh,Gsph_x,Gsph_c,Vcp,&
970 : ! & is_resonant,rhxtwg_q0,nproc,my_rank,t_start,t_stop,my_bsham,comm)
971 :
972 24 : call timab(683,1,tsec) ! exc_build_ham(exchange)
973 :
974 24 : write(msg,'(a,2i2,a)')" Calculating exchange term for (spin1,spin2) ",spin1,spin2," ..."
975 24 : call wrtout(std_out, msg)
976 :
977 72 : ABI_MALLOC(rhotwg1,(npweps))
978 48 : ABI_MALLOC(rhotwg2,(npweps))
979 :
980 24 : ngx = Gsph_x%ng
981 72 : ABI_MALLOC(vc_sqrt_qbz,(ngx))
982 :
983 : ! * Get iq_ibz, and symmetries from iq_bz.
984 24 : iq_bz = iqbz0 ! q = 0 -> iqbz0
985 24 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
986 :
987 : ! * Set up table of |q(BZ)+G|
988 24 : if (iq_ibz==1) then
989 7078 : do ig=1,ngx
990 7054 : ISg = Gsph_x%rottb(ig,itim_q,isym_q)
991 7078 : vc_sqrt_qbz(ISg)=Vcp%vcqlwl_sqrt(ig,1)
992 : end do
993 : else
994 0 : ABI_ERROR("iq_ibz should be 1")
995 : end if
996 :
997 9272 : do itp=1,BSp%nreh(block) ! Loop over transition tp = (kp,vp,cp,spin2)
998 :
999 9248 : if (itp<my_cols(1) .or. itp>my_cols(2)) CYCLE ! I dont have this column.
1000 7716 : ene_tp = Bsp%Trans(itp,spin2)%en
1001 7716 : ikp_bz = Bsp%Trans(itp,spin2)%k
1002 7716 : ivp = Bsp%Trans(itp,spin2)%v
1003 7716 : icp = Bsp%Trans(itp,spin2)%c
1004 :
1005 7716 : ikp_ibz = Kmesh%tab (ikp_bz)
1006 7716 : isym_kp = Kmesh%tabo(ikp_bz)
1007 7716 : itim_kp = (3-Kmesh%tabi(ikp_bz))/2
1008 :
1009 7716 : if (is_resonant) then
1010 7332 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,ivp,icp,ikp_ibz,spin2),Gsph_c)
1011 : else ! Code for coupling block.
1012 384 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,icp,ivp,ikp_ibz,spin2),Gsph_c)
1013 : end if
1014 : !
1015 : ! Multiply by the Coulomb term.
1016 365516 : do ig=2,npweps
1017 365516 : rhotwg2(ig) = rhotwg2(ig) * vc_sqrt_qbz(ig) * vc_sqrt_qbz(ig)
1018 : end do
1019 :
1020 1886073 : do it=1,itp ! Loop over transition t = (k,v,c,spin1)
1021 1878333 : ir = it + itp*(itp-1_i8b)/2
1022 1878333 : if (ir<t_start(my_rank) .or. ir>t_stop(my_rank)) CYCLE
1023 :
1024 1877264 : ene_t = Bsp%Trans(it,spin1)%en
1025 1877264 : ik_bz = Bsp%Trans(it,spin1)%k
1026 1877264 : iv = Bsp%Trans(it,spin1)%v
1027 1877264 : ic = Bsp%Trans(it,spin1)%c
1028 :
1029 1877264 : ik_ibz = Kmesh%tab(ik_bz)
1030 1877264 : isym_k = Kmesh%tabo(ik_bz)
1031 1877264 : itim_k = (3-Kmesh%tabi(ik_bz))/2
1032 : !if (itim_k==2) CYCLE ! time-reversal or not
1033 :
1034 1877264 : rhotwg1(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,iv,ic,ik_ibz,spin1),Gsph_c)
1035 : !
1036 : ! sum over G/=0
1037 1877264 : ctemp = xdotc(npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
1038 1877264 : ctemp = faq * kx_fact * ctemp
1039 :
1040 : ! exchange term is non divergent !
1041 1877264 : if (BSp%prep_interp) then
1042 4656 : ccoeffs(ir) = ccoeffs(ir) + ctemp
1043 : end if
1044 :
1045 1887581 : my_bsham(ir) = my_bsham(ir) + ctemp
1046 : end do !it
1047 : end do !itp
1048 :
1049 24 : ABI_FREE(rhotwg1)
1050 24 : ABI_FREE(rhotwg2)
1051 24 : ABI_FREE(vc_sqrt_qbz)
1052 :
1053 24 : call timab(683,2,tsec) ! exc_build_ham(exchange)
1054 : end if ! do_exchange_term
1055 : !
1056 : ! =====================
1057 : ! === Diagonal term ===
1058 : ! =====================
1059 24 : if (is_resonant .and. spin1==spin2) then
1060 23 : write(msg,'(a,2i2,a)')" Adding diagonal term for (spin1,spin2) ",spin1,spin2," ..."
1061 23 : call wrtout(std_out, msg)
1062 8887 : do it=1,BSp%nreh(block)
1063 8864 : ir = it + it*(it-1_i8b)/2
1064 8887 : if (ir>=t_start(my_rank) .and. ir<=t_stop(my_rank)) my_bsham(ir) = my_bsham(ir) + Bsp%Trans(it,spin1)%en
1065 : end do
1066 : end if
1067 :
1068 : if (.FALSE.) then
1069 : dump_unt = get_unit()
1070 : msg=' Coupling Hamiltonian matrix elements: '
1071 : if (is_resonant) msg=' Reasonant Hamiltonian matrix elements: '
1072 : call wrtout(dump_unt, msg)
1073 : call wrtout(dump_unt,' k v c s k" v" c" s" H')
1074 : do itp=1,BSp%nreh(block)
1075 : ikp_bz = Bsp%Trans(itp,spin2)%k
1076 : ivp = Bsp%Trans(itp,spin2)%v
1077 : icp = Bsp%Trans(itp,spin2)%c
1078 : do it=1,itp
1079 : ik_bz = Bsp%Trans(it,spin1)%k
1080 : iv = Bsp%Trans(it,spin1)%v
1081 : ic = Bsp%Trans(it,spin1)%c
1082 : ir = it + itp*(itp-1_i8b)/2
1083 : if (ir>=t_start(my_rank).and.ir<=t_stop(my_rank)) then
1084 : http = my_bsham(ir)
1085 : !if (ABS(http) > tol3) then
1086 : write(msg,'(2(i0,1x),2(i5,3i3,3x),2f7.3)')it,itp, ik_bz,iv,ic,spin1, ikp_bz,ivp,icp,spin2, http
1087 : call wrtout(dump_unt, msg)
1088 : !end if
1089 : end if
1090 : end do
1091 : end do
1092 : end if
1093 :
1094 : !DBYG
1095 : if (.False.) then
1096 : dump_unt = get_unit()
1097 : dump_unt = 999
1098 : msg=' Coupling Hamiltonian matrix elements: '
1099 : if (is_resonant) msg=' Resonant Hamiltonian matrix elements: '
1100 : call wrtout(dump_unt, msg)
1101 : call wrtout(dump_unt,' k v c s k" v" c" s" H')
1102 : do itp=1,BSp%nreh(block)
1103 : ikp_bz = Bsp%Trans(itp,spin2)%k
1104 : ivp = Bsp%Trans(itp,spin2)%v
1105 : icp = Bsp%Trans(itp,spin2)%c
1106 : do it=1,BSp%nreh(block)
1107 : ik_bz = Bsp%Trans(it,spin1)%k
1108 : iv = Bsp%Trans(it,spin1)%v
1109 : ic = Bsp%Trans(it,spin1)%c
1110 : if(it > itp) then
1111 : ir = itp+it*(it-1_i8b)/2
1112 : else
1113 : ir = it + itp*(itp-1_i8b)/2
1114 : end if
1115 : if (ir>=t_start(my_rank).and.ir<=t_stop(my_rank)) then
1116 : if(it > itp) then
1117 : http = CONJG(my_bsham(ir))
1118 : if (BSp%prep_interp) then
1119 : aatmp = CONJG(acoeffs(ir))
1120 : bbtmp = CONJG(bcoeffs(ir))
1121 : cctmp = CONJG(ccoeffs(ir))
1122 : end if
1123 : else
1124 : http = my_bsham(ir)
1125 : if (BSp%prep_interp) then
1126 : aatmp = acoeffs(ir)
1127 : bbtmp = bcoeffs(ir)
1128 : cctmp = ccoeffs(ir)
1129 : end if
1130 : end if
1131 : if (it == itp) http = http - Bsp%Trans(it,spin1)%en
1132 : !if (ABS(http) > tol3) then
1133 : if (BSp%prep_interp) then
1134 : write(msg,'(2(i0,1x),2(i5,3i3,3x),2f24.20,2f24.20,2f24.20,2f24.20)')it,itp, ik_bz,iv,ic,spin1, ikp_bz,ivp,icp,&
1135 : & spin2, http, aatmp, bbtmp, cctmp
1136 : else
1137 : write(msg,'(2(i0,1x),2(i5,3i3,3x),2f24.20)')it,itp, ik_bz,iv,ic,spin1, ikp_bz,ivp,icp,spin2, http
1138 : end if
1139 : call wrtout(dump_unt, msg)
1140 : !end if
1141 : end if
1142 : end do
1143 : end do
1144 : end if
1145 :
1146 24 : call timab(684,1,tsec) ! exc_build_ham(synchro)
1147 24 : call xmpi_barrier(comm)
1148 24 : call timab(684,2,tsec) ! exc_build_ham(synchro)
1149 : !
1150 : ! =================================
1151 : ! === Write Hamiltonian on disk ===
1152 : ! =================================
1153 24 : call timab(685,1,tsec) ! exc_build_ham(write_ham)
1154 48 : if (use_mpiio) then
1155 : #ifdef HAVE_MPI_IO
1156 : ! Write the Hamiltonian with collective MPI-IO.
1157 : if (BSp%prep_interp) then
1158 : ABI_ERROR("Preparation of interpolation technique not yet coded with MPI-IO")
1159 : end if
1160 : ABI_CHECK(nsppol==1,"nsppol==2 not coded, offset is wrong")
1161 : !
1162 : old_type = MPI_DOUBLE_COMPLEX
1163 : call xmpio_create_fherm_packed(my_starts,my_ends,is_fortran_file,my_offset,old_type,hmat_type,offset_err)
1164 :
1165 : if (offset_err/=0) then
1166 : write(msg,"(3a)")&
1167 : "Global position index cannot be stored in a standard Fortran integer. ",ch10,&
1168 : "BSE matrix cannot be written with a single MPI-IO call. "
1169 : ABI_ERROR(msg)
1170 : end if
1171 : !
1172 : ! Each node uses a different offset to skip the header and the blocks written by the other CPUs.
1173 : my_offset = offset_of_block(block) + my_offset
1174 :
1175 : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, MPI_BYTE, hmat_type, 'native', MPI_INFO_NULL, mpi_err)
1176 : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1177 :
1178 : call MPI_TYPE_FREE(hmat_type,mpi_err)
1179 : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1180 :
1181 : if (hsize_of(my_rank) /= INT(hsize_of(my_rank),kind=i4b) ) then
1182 : ABI_ERROR("Wraparound error")
1183 : end if
1184 :
1185 : tmp_size = INT(hsize_of(my_rank))
1186 : call MPI_FILE_WRITE_ALL(mpi_fh, my_bsham, tmp_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1187 : ABI_CHECK_MPI(mpi_err,"FILE_WRITE")
1188 :
1189 : ! It seems that personal calls in make the code stuck
1190 : !if (is_fortran_file .and. my_rank==master) then ! Master writes the Fortran record markers.
1191 : ! Write the Fortran record markers.
1192 : neh2=BSp%nreh(block)
1193 : ABI_MALLOC(bsize_frecord,(neh2))
1194 : bsize_frecord = (/(col_glob * xmpi_bsize_dpc, col_glob=1,neh2)/)
1195 : ! ehdr_offset points to the end of the header.
1196 : !call xmpio_write_frmarkers(mpi_fh,ehdr_offset,xmpio_collective,neh2,bsize_frecord,mpi_err)
1197 : my_offset = offset_of_block(block)
1198 : call xmpio_write_frmarkers(mpi_fh,my_offset,xmpio_collective,neh2,bsize_frecord,ierr)
1199 : ABI_CHECK(ierr==0,"Error while writing Fortran markers")
1200 : ABI_FREE(bsize_frecord)
1201 : #else
1202 : ABI_BUG("You should not be here!")
1203 : #endif
1204 : else
1205 : ! Use FORTRAN IO with sequential access mode.
1206 : ! * Each node sends its data to master node.
1207 : ! * Blocks are distributed according to the rank of the node.
1208 : ! * Matrix is written by columns hence make sure that the last column is completely written.
1209 24 : call cwtime(cputime,walltime,gflops,"start")
1210 :
1211 24 : if (my_rank==master) then
1212 20 : prev_nrows=0; if (my_cols(2) /= my_rows(2)) prev_nrows = my_rows(2)
1213 20 : ncol = my_cols(2)-my_cols(1)+1
1214 20 : ist=1
1215 7428 : do jj=1,ncol
1216 7408 : col_glob = my_starts(2) + jj - 1
1217 7408 : nrows = col_glob; if (jj==ncol) nrows=my_rows(2)
1218 7408 : iend = ist + nrows -1
1219 7408 : write(bsh_unt) my_bsham(ist:iend)
1220 7408 : if (BSp%prep_interp) then
1221 96 : write(a_unt) acoeffs(ist:iend)
1222 96 : write(b_unt) bcoeffs(ist:iend)
1223 96 : write(c_unt) ccoeffs(ist:iend)
1224 : end if
1225 7428 : ist=iend+1
1226 : end do
1227 20 : write(msg,'(2(a,i0))')" Wraparound error: iend=",iend," my_hsize=",hsize_of(my_rank)
1228 20 : ABI_CHECK(iend == hsize_of(my_rank),msg)
1229 20 : ABI_FREE(my_bsham)
1230 20 : if (BSp%prep_interp) then
1231 1 : ABI_FREE(acoeffs)
1232 1 : ABI_FREE(bcoeffs)
1233 1 : ABI_FREE(ccoeffs)
1234 : end if
1235 : end if
1236 :
1237 24 : call xmpi_barrier(comm)
1238 : !
1239 : ! Collect data from the other nodes.
1240 38 : do sender=1,nproc-1
1241 : ! If I'm not involved, jump to the end of the loop and wait there (sequential IO? Of course!)
1242 14 : if (all(my_rank /= [sender, master])) goto 100
1243 :
1244 8 : if (my_rank==master) then
1245 12 : ABI_MALLOC(buffer,(hsize_of(sender)))
1246 4 : if (BSp%prep_interp) then
1247 0 : ABI_MALLOC(abuffer,(hsize_of(sender)))
1248 0 : ABI_MALLOC(bbuffer,(hsize_of(sender)))
1249 0 : ABI_MALLOC(cbuffer,(hsize_of(sender)))
1250 : end if
1251 : end if
1252 8 : tmp_size = INT(hsize_of(sender),kind=i4b)
1253 8 : call xmpi_exch(my_bsham,tmp_size,sender,buffer,master,comm,10*block+1,mpi_err)
1254 8 : if (BSp%prep_interp) then
1255 0 : call xmpi_exch(acoeffs,tmp_size,sender,abuffer,master,comm,10*block+2,mpi_err)
1256 0 : call xmpi_exch(bcoeffs,tmp_size,sender,bbuffer,master,comm,10*block+3,mpi_err)
1257 0 : call xmpi_exch(ccoeffs,tmp_size,sender,cbuffer,master,comm,10*block+4,mpi_err)
1258 : end if
1259 :
1260 : ! TODO Be careful with the MPI TAG here, add optional Arguments in xmpi_exch so that the TAG can be specified!
1261 24 : proc_start = (/my_rows(1),my_cols(1)/)
1262 24 : proc_end = (/my_rows(2),my_cols(2)/)
1263 24 : my_extrema(:,1) = proc_start
1264 24 : my_extrema(:,2) = proc_end
1265 :
1266 8 : sender_extrema = my_extrema ! just to avoid NAN on sender. xechh_mpi is not well designed
1267 8 : call xmpi_exch(my_extrema,4,sender,sender_extrema,master,comm,10*block+5,mpi_err)
1268 :
1269 8 : if (my_rank==master) then
1270 12 : proc_start = sender_extrema(:,1)
1271 12 : proc_end = sender_extrema(:,2)
1272 : !write(std_out,*)"proc_start, proc_end",proc_start,proc_end
1273 :
1274 4 : if (prev_nrows>0) then ! backspace the file if the last record written was not complete.
1275 : !write(std_out,*)" master node had to call backspace"
1276 4 : backspace(bsh_unt)
1277 12 : ABI_MALLOC(prev_col,(prev_nrows))
1278 4 : read(bsh_unt) prev_col
1279 4 : backspace(bsh_unt)
1280 :
1281 4 : if (BSp%prep_interp) then
1282 0 : backspace(a_unt)
1283 0 : ABI_MALLOC(aprev_col,(prev_nrows))
1284 0 : read(a_unt) aprev_col
1285 0 : backspace(a_unt)
1286 :
1287 0 : backspace(b_unt)
1288 0 : ABI_MALLOC(bprev_col,(prev_nrows))
1289 0 : read(b_unt) bprev_col
1290 0 : backspace(b_unt)
1291 :
1292 0 : backspace(c_unt)
1293 0 : ABI_MALLOC(cprev_col,(prev_nrows))
1294 0 : read(c_unt) cprev_col
1295 0 : backspace(c_unt)
1296 : end if
1297 : end if
1298 : !
1299 : ! Write the columns owned by sender.
1300 4 : ncol = proc_end(2)-proc_start(2)+1
1301 4 : ist=1
1302 312 : do jj=1,ncol
1303 308 : col_glob = proc_start(2) + jj-1
1304 308 : nrows = col_glob
1305 308 : if (jj==1 ) nrows=col_glob - proc_start(1) + 1
1306 308 : if (jj==ncol) then
1307 4 : nrows=proc_end(1)
1308 4 : if (ncol==1) nrows=proc_end(1) - proc_start(1) + 1
1309 : end if
1310 308 : iend = ist + nrows -1
1311 : !write(std_out,*)"Using nrows, ist, iend=",nrows,ist,iend
1312 308 : if (jj==1 .and. prev_nrows>0) then ! join prev_col and this subcolumn.
1313 1073 : write(bsh_unt) CMPLX(prev_col,kind=dp),CMPLX(buffer(ist:iend),kind=dp)
1314 4 : if (BSp%prep_interp) then
1315 0 : write(a_unt) CMPLX(aprev_col,kind=dp),CMPLX(abuffer(ist:iend),kind=dp)
1316 0 : write(b_unt) CMPLX(bprev_col,kind=dp),CMPLX(bbuffer(ist:iend),kind=dp)
1317 0 : write(c_unt) CMPLX(cprev_col,kind=dp),CMPLX(cbuffer(ist:iend),kind=dp)
1318 : end if
1319 4 : prev_nrows = prev_nrows + iend-ist+1
1320 : else
1321 92149 : write(bsh_unt) CMPLX(buffer(ist:iend),kind=dp)
1322 304 : if (BSp%prep_interp) then
1323 0 : write(a_unt) CMPLX(abuffer(ist:iend),kind=dp)
1324 0 : write(b_unt) CMPLX(bbuffer(ist:iend),kind=dp)
1325 0 : write(c_unt) CMPLX(cbuffer(ist:iend),kind=dp)
1326 : end if
1327 : prev_nrows=0
1328 : end if
1329 312 : ist=iend+1
1330 : end do
1331 4 : if (ncol>1) then ! Reset prev_nrows if a new column has begun.
1332 4 : prev_nrows = proc_end(1)
1333 4 : if (proc_end(1) == proc_end(2)) prev_nrows = 0
1334 : end if
1335 4 : if (iend/=hsize_of(sender)) then
1336 0 : write(msg,'(2(a,i0))')" Wraparound error: iend=",iend," my_hsize=",hsize_of(sender)
1337 0 : ABI_ERROR(msg)
1338 : end if
1339 4 : ABI_SFREE(prev_col)
1340 4 : if (BSp%prep_interp) then
1341 0 : ABI_SFREE(aprev_col)
1342 0 : ABI_SFREE(bprev_col)
1343 0 : ABI_SFREE(cprev_col)
1344 : end if
1345 4 : ABI_FREE(buffer)
1346 4 : if (BSp%prep_interp) then
1347 0 : ABI_FREE(abuffer)
1348 0 : ABI_FREE(bbuffer)
1349 0 : ABI_FREE(cbuffer)
1350 : end if
1351 : end if ! master
1352 : !
1353 38 : 100 call xmpi_barrier(comm)
1354 : end do ! sender
1355 :
1356 24 : call cwtime(cputime,walltime,gflops,"stop")
1357 24 : write(msg,'(2(a,f9.1),a)')" Fortran-IO completed. cpu_time: ",cputime,"[s], walltime: ",walltime," [s]"
1358 24 : call wrtout(std_out, msg, do_flush=.True.)
1359 : end if ! use_mpiio
1360 24 : call timab(685,2,tsec) ! exc_build_ham(write_ham)
1361 : !
1362 24 : ABI_SFREE(my_bsham)
1363 24 : if (BSp%prep_interp) then
1364 1 : ABI_SFREE(acoeffs)
1365 1 : ABI_SFREE(bcoeffs)
1366 1 : ABI_SFREE(ccoeffs)
1367 : end if
1368 24 : ABI_FREE(t_start)
1369 24 : ABI_FREE(t_stop)
1370 47 : ABI_FREE(hsize_of)
1371 : end do ! block
1372 : !
1373 : ! ===========================================
1374 : ! === Exchange term for spin_up spin_down ===
1375 : ! ===========================================
1376 :
1377 23 : if (nsppol==2) then
1378 1 : call timab(686,2,tsec) ! exc_build_ham(exch.spin)
1379 1 : block=3
1380 1 : neh1=BSp%nreh(1)
1381 1 : neh2=BSp%nreh(2)
1382 : !
1383 : ! The oscillators at q=0 are available on each node for both spin.
1384 : ! Here the calculation of the block is parallelized over columns.
1385 3 : ABI_MALLOC(col_start,(0:nproc-1))
1386 2 : ABI_MALLOC(col_stop,(0:nproc-1))
1387 1 : call xmpi_split_work2_i4b(neh2,nproc,col_start,col_stop)
1388 :
1389 1 : my_cols(1) = col_start(my_rank)
1390 1 : my_cols(2) = col_stop (my_rank)
1391 1 : if (my_cols(2)-my_cols(1)<=0) then
1392 0 : ABI_ERROR("One of the processors has zero columns!")
1393 : end if
1394 :
1395 2 : ABI_MALLOC(ncols_of,(0:nproc-1))
1396 2 : ncols_of=0
1397 2 : do rank=0,nproc-1
1398 2 : if (col_stop(rank)>=col_start(rank)) ncols_of(rank) = col_stop(rank)-col_start(rank)+1
1399 : end do
1400 :
1401 1 : ABI_FREE(col_start)
1402 1 : ABI_FREE(col_stop)
1403 : !
1404 : ! TODO might add treatment of <psi|q+G|psi> for q+G -> 0
1405 : ! TODO might used enlarged G-sphere for better convergence.
1406 : ! Note that my_kxssp is always written on file when nsppol=2, even when
1407 : ! non-local field effects are neglected.
1408 4 : ABI_MALLOC(my_kxssp,(neh1,my_cols(1):my_cols(2)))
1409 147841 : my_kxssp=czero
1410 :
1411 1 : if (do_exchange_term) then
1412 1 : spin1=1; spin2=2
1413 1 : write(msg,'(a,2i2,a)')" Calculating exchange term for (spin1,spin2) ",spin1,spin2," ..."
1414 1 : call wrtout(std_out, msg)
1415 :
1416 3 : ABI_MALLOC(rhotwg1,(npweps))
1417 2 : ABI_MALLOC(rhotwg2,(npweps))
1418 :
1419 1 : ngx = Gsph_x%ng
1420 3 : ABI_MALLOC(vc_sqrt_qbz,(ngx))
1421 : !
1422 : ! * Get iq_ibz, and symmetries from iq_bz.
1423 1 : iq_bz = iqbz0 ! q = 0 -> iqbz0
1424 1 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
1425 : !
1426 : ! * Set up table of |q(BZ)+G|
1427 1 : if (iq_ibz==1) then
1428 170 : do ig=1,ngx
1429 169 : ISg = Gsph_x%rottb(ig,itim_q,isym_q)
1430 170 : vc_sqrt_qbz(ISg)=Vcp%vcqlwl_sqrt(ig,1)
1431 : end do
1432 : else
1433 0 : ABI_ERROR("iq_ibz should be 1")
1434 : end if
1435 :
1436 385 : do itp=1,neh2 ! Loop over transition tp = (kp,vp,cp,spin2)
1437 :
1438 384 : if (itp<my_cols(1) .or. itp>my_cols(2)) CYCLE ! I dont have this column.
1439 384 : ene_tp = Bsp%Trans(itp,spin2)%en
1440 384 : ikp_bz = Bsp%Trans(itp,spin2)%k
1441 384 : ivp = Bsp%Trans(itp,spin2)%v
1442 384 : icp = Bsp%Trans(itp,spin2)%c
1443 :
1444 384 : ikp_ibz = Kmesh%tab (ikp_bz)
1445 384 : isym_kp = Kmesh%tabo(ikp_bz)
1446 384 : itim_kp = (3-Kmesh%tabi(ikp_bz))/2
1447 :
1448 384 : if (is_resonant) then
1449 384 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,ivp,icp,ikp_ibz,spin2),Gsph_c)
1450 : else ! Code for coupling block.
1451 0 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,icp,ivp,ikp_ibz,spin2),Gsph_c)
1452 : end if
1453 : !
1454 : ! Multiply by the Coulomb term.
1455 19584 : do ig=2,npweps
1456 19584 : rhotwg2(ig) = rhotwg2(ig) * vc_sqrt_qbz(ig) * vc_sqrt_qbz(ig)
1457 : end do
1458 :
1459 147841 : do it=1,neh1 ! Loop over transition t = (k,v,c,spin1) FULL matrix.
1460 :
1461 147456 : ene_t = Bsp%Trans(it,spin1)%en
1462 147456 : ik_bz = Bsp%Trans(it,spin1)%k
1463 147456 : iv = Bsp%Trans(it,spin1)%v
1464 147456 : ic = Bsp%Trans(it,spin1)%c
1465 :
1466 147456 : ik_ibz = Kmesh%tab(ik_bz)
1467 147456 : isym_k = Kmesh%tabo(ik_bz)
1468 147456 : itim_k = (3-Kmesh%tabi(ik_bz))/2
1469 : !if (itim_k==2) CYCLE ! time-reversal or not
1470 :
1471 147456 : rhotwg1(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,iv,ic,ik_ibz,spin1),Gsph_c)
1472 : !
1473 : ! sum over G/=0
1474 147456 : ctemp = XDOTC(npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
1475 147456 : ctemp = faq * kx_fact * ctemp
1476 :
1477 147840 : my_kxssp(it,itp) = ctemp
1478 : end do !it
1479 : end do !itp
1480 :
1481 1 : ABI_FREE(rhotwg1)
1482 1 : ABI_FREE(rhotwg2)
1483 1 : ABI_FREE(vc_sqrt_qbz)
1484 : end if ! do_exchange_term
1485 1 : call timab(686,2,tsec) ! exc_build_ham(exch.spin)
1486 : !
1487 : ! =====================================
1488 : ! === Write the Hamiltonian on disk ===
1489 : ! =====================================
1490 1 : call timab(685,1,tsec) ! exc_build_ham(write_ham)
1491 :
1492 : if (use_mpiio) then
1493 : #ifdef HAVE_MPI_IO
1494 : my_ncols=ncols_of(my_rank); old_type=MPI_DOUBLE_COMPLEX
1495 : call xmpio_create_fsubarray_2D((/neh1,my_ncols/),(/neh1,my_ncols/),(/1,1/),old_type,hmat_type,my_offpad,mpi_err)
1496 : ABI_CHECK_MPI(mpi_err,"fsubarray_2D")
1497 : !
1498 : ! Each node uses a different offset to skip the header and the blocks written by the other CPUs.
1499 : prev_nels=0
1500 : prev_ncols=0
1501 : if (my_rank>0) then
1502 : prev_ncols = SUM(ncols_of(0:my_rank-1))
1503 : prev_nels = neh1*prev_ncols
1504 : end if
1505 : tmp_off = prev_nels*xmpi_bsize_dpc + prev_ncols*2*xmpio_bsize_frm
1506 :
1507 : my_offset = offset_of_block(block) + tmp_off + my_offpad
1508 :
1509 : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, MPI_BYTE, hmat_type, 'native', MPI_INFO_NULL, mpi_err)
1510 : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1511 :
1512 : call MPI_TYPE_FREE(hmat_type,mpi_err)
1513 : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1514 :
1515 : tmp_size = INT(neh1*my_ncols)
1516 : call MPI_FILE_WRITE_ALL(mpi_fh, my_kxssp,tmp_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1517 : ABI_CHECK_MPI(mpi_err,"FILE_WRITE")
1518 :
1519 : ! It seems that personal calls in make the code stuck
1520 : ! Master writes the Fortran record markers.
1521 : ABI_MALLOC(bsize_frecord,(neh2))
1522 : bsize_frecord = neh1 * xmpi_bsize_dpc
1523 : ! ehdr_offset points to the end of the header.
1524 : !call xmpio_write_frmarkers(mpi_fh,ehdr_offset,xmpio_collective,neh2,bsize_frecord,mpi_err)
1525 : my_offset = offset_of_block(block)
1526 : call xmpio_write_frmarkers(mpi_fh,my_offset,xmpio_collective,neh2,bsize_frecord,ierr)
1527 : ABI_CHECK(ierr==0,"Error while writing Fortran markers")
1528 : ABI_FREE(bsize_frecord)
1529 : #else
1530 : ABI_BUG("You should not be here")
1531 : #endif
1532 : else
1533 : ! Use FORTRAN IO with sequential access mode.
1534 : ! * Each node sends its data to master node.
1535 : ! * Columns are distributed according to the rank of the node.
1536 1 : if (my_rank==master) then
1537 385 : do jj=my_cols(1),my_cols(2)
1538 385 : write(bsh_unt) my_kxssp(:,jj)
1539 : end do
1540 1 : ABI_FREE(my_kxssp)
1541 : end if
1542 :
1543 1 : call xmpi_barrier(comm)
1544 : !
1545 : ! Collect data from the other nodes.
1546 1 : do sender=1,nproc-1
1547 : ! If I'm not involved, jump to the end of the loop and wait there (sequential IO? Of course!)
1548 0 : if (all(my_rank /= [sender, master])) goto 200
1549 :
1550 0 : if (my_rank==master) then
1551 0 : ABI_MALLOC(buffer_2d,(neh1,ncols_of(sender)))
1552 : end if
1553 0 : call xmpi_exch(my_kxssp,neh1*ncols_of(sender),sender,buffer_2d,master,comm,5,mpi_err)
1554 : !
1555 0 : if (my_rank==master) then ! Write the columns owned by sender.
1556 0 : do jj=1,ncols_of(sender)
1557 0 : write(bsh_unt) buffer_2d(:,jj)
1558 : end do
1559 0 : ABI_FREE(buffer_2d)
1560 : end if ! master
1561 : !
1562 1 : 200 call xmpi_barrier(comm)
1563 : end do ! sender
1564 : end if
1565 1 : call timab(685,2,tsec) ! exc_build_ham(write_ham)
1566 :
1567 1 : ABI_FREE(ncols_of)
1568 1 : ABI_SFREE(my_kxssp)
1569 : end if
1570 :
1571 : ! Close the file.
1572 : if (use_mpiio) then
1573 : #ifdef HAVE_MPI_IO
1574 : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1575 : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1576 : ABI_FREE(offset_of_block)
1577 : #endif
1578 : end if
1579 :
1580 : ! master closes the Fortran files.
1581 23 : if (my_rank==master) then
1582 19 : close(bsh_unt)
1583 19 : if (BSp%prep_interp) then
1584 1 : close(a_unt)
1585 1 : close(b_unt)
1586 1 : close(c_unt)
1587 : end if
1588 : end if
1589 :
1590 : ! Free memory.
1591 23 : ABI_FREE(igfftg0)
1592 23 : ABI_FREE(ktabr_k)
1593 23 : ABI_FREE(id_tab)
1594 23 : ABI_FREE(ktabr_kp)
1595 23 : ABI_FREE(rhxtwg_vpv)
1596 23 : ABI_FREE(rhxtwg_cpc)
1597 23 : if (BSp%prep_interp) then
1598 1 : ABI_FREE(aa_vpv)
1599 1 : ABI_FREE(bb_vpv1)
1600 1 : ABI_FREE(bb_vpv2)
1601 1 : ABI_FREE(cc_vpv)
1602 1 : ABI_FREE(aa_cpc)
1603 1 : ABI_FREE(bb_cpc1)
1604 1 : ABI_FREE(bb_cpc2)
1605 1 : ABI_FREE(cc_cpc)
1606 : end if
1607 23 : ABI_FREE(ur_ckp)
1608 23 : ABI_FREE(ur_vkp)
1609 23 : ABI_FREE(ur_vk)
1610 23 : ABI_FREE(ur_ck)
1611 :
1612 : ! Deallocation for PAW.
1613 23 : if (Wfd%usepaw==1) then
1614 2 : call pawcprj_free(Cp_vk)
1615 6 : ABI_FREE(Cp_vk)
1616 2 : call pawcprj_free(Cp_ck)
1617 6 : ABI_FREE(Cp_ck)
1618 2 : call pawcprj_free(Cp_ckp)
1619 6 : ABI_FREE(Cp_ckp)
1620 2 : call pawcprj_free(Cp_vkp)
1621 6 : ABI_FREE(Cp_vkp)
1622 2 : call pawcprj_free(Cp_tmp1)
1623 6 : ABI_FREE(Cp_tmp1)
1624 2 : call pawcprj_free(Cp_tmp2)
1625 6 : ABI_FREE(Cp_tmp2)
1626 2 : call pawcprj_free(Cp_tmp3)
1627 6 : ABI_FREE(Cp_tmp3)
1628 2 : call pawcprj_free(Cp_tmp4)
1629 6 : ABI_FREE(Cp_tmp4)
1630 : end if
1631 :
1632 23 : call xmpi_barrier(comm)
1633 :
1634 : DBG_EXIT("COLL")
1635 :
1636 23 : call timab(680,2,tsec)
1637 :
1638 46 : end subroutine exc_build_block
1639 : !!***
1640 :
1641 : !!****f* m_exc_build/exc_build_v
1642 : !! NAME
1643 : !! exc_build_v
1644 : !!
1645 : !! FUNCTION
1646 : !! Calculate and write the excitonic Hamiltonian on an external binary file (Fortran file open
1647 : !! in random mode) for subsequent treatment in the Bethe-Salpeter code.
1648 : !!
1649 : !! INPUTS
1650 : !! BSp<excparam>=The parameters for the Bethe-Salpeter calculation.
1651 : !! Cryst<crystal_t>=Info on the crystalline structure.
1652 : !! Kmesh<kmesh_t>=The list of k-points in the BZ, IBZ and symmetry tables.
1653 : !! Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
1654 : !! Gsph_x<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored).
1655 : !! Gsph_c<gsphere_t>=Info on the G-sphere used to describe the correlation part.
1656 : !! Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
1657 : !! rhxtwg_q0
1658 : !! is_resonant
1659 : !! comm=MPI communicator.
1660 : !!
1661 : !! OUTPUT
1662 : !!
1663 : !! NOTES
1664 : !! *) Version for K_V = K_C (q=0), thus KP_V = KP_C
1665 : !! *) No exchange limit: use DFT energies in case.
1666 : !! *) Symmetry of H(-k-k') = H*(k k') not used.
1667 : !! *) Coulomb term can be approximated as diagonal in G.
1668 : !! *) Valence bands treated from lomo on.
1669 : !! *) Symmetries of the sub-blocks are used to reduce the number of elements to calculate.
1670 : !!
1671 : !! ____________
1672 : !! |_(cv)__(vc)_|
1673 : !! H_exc = | R C |
1674 : !! | -C* -R* |
1675 : !!
1676 : !! where C is symmetric and R is Hermitian provided that the QP energies are real.
1677 : !!
1678 : !! For nsppol=1 ==> R = diag-W+2v; C = -W+2v
1679 : !! since the Hamiltonian can be diagonalized in the spin-singlet basis set thanks to
1680 : !! the fact that spin triplet does not contribute to the optical limit of epsilon.
1681 : !!
1682 : !! For nsppol=2 ==> R = diag-W+v; C = -W+v
1683 : !! Now the matrix elements depend on the spin of the transitions but only those
1684 : !! transitions in which the spin of the electron and of the hole are equal contribute
1685 : !! to the macroscopic dielectric function. Moreover only the exchange term can connect transitions of different spin.
1686 : !! When nsppol==2 the transitions are ordered using | (cv up) | (cv dwn) | (vc up) | (vc down) |
1687 : !!
1688 : !! The resonant block is given by:
1689 : !!
1690 : !! | (v'c' up) | (v'c' dwn) |
1691 : !! ----------------------------------- where v_{-+} = v_{+-}^H when the momentum of the photon is neglected.
1692 : !! | [diag-W+v]++ | v+- | (vc up) Note that v_{+-} is not Hermitian due to the presence of different spins.
1693 : !! R = ----------------------------------- Actually it reduces to a Hermitian matrix when the system is not spin polarized.
1694 : !! | v-+ | [diag-W+v]-- | (vc dwn) but in this case one should use nsppol=1.
1695 : !! ----------------------------------- As a consequence the entire matrix is calculated and stored on file.
1696 : !!
1697 : !! The coupling block is given by:
1698 : !!
1699 : !! | (c'v' up) | (c'v dwn) |
1700 : !! ----------------------------------- where v_{-+} = v_{+-}^t when the momentum of the photon is neglected.
1701 : !! | [-W+v]++ | v+- | (vc up) Also in this case the entire matrix v_{+-} has to be calculated
1702 : !! C = ----------------------------------- and stored on file.
1703 : !! | v-+ | [-W+v]-- | (vc dwn)
1704 : !! -----------------------------------
1705 : !!
1706 : !! SOURCE
1707 :
1708 : subroutine exc_build_v(spin1,spin2,nsppol,npweps,Bsp,Cryst,Kmesh,Qmesh,Gsph_x,Gsph_c,Vcp,&
1709 : & is_resonant,rhxtwg_q0,nproc,my_rank,t_start,t_stop,my_bsham)
1710 :
1711 : !Arguments ------------------------------------
1712 : !scalars
1713 : integer,intent(in) :: spin1,spin2,nsppol,npweps,nproc,my_rank
1714 : logical,intent(in) :: is_resonant
1715 : type(excparam),intent(in) :: BSp
1716 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
1717 : type(crystal_t),intent(in) :: Cryst
1718 : type(vcoul_t),intent(in) :: Vcp
1719 : type(gsphere_t),intent(in) :: Gsph_x,Gsph_c
1720 : !arrays
1721 : integer(i8b),intent(in) :: t_start(0:nproc-1),t_stop(0:nproc-1)
1722 : complex(gwp),intent(in) :: rhxtwg_q0(npweps,BSp%lomo_min:BSp%humo_max,BSp%lomo_min:BSp%humo_max,Kmesh%nibz,nsppol)
1723 : complex(dp),intent(inout) :: my_bsham(t_start(my_rank):t_stop(my_rank))
1724 :
1725 : !Local variables ------------------------------
1726 : !scalars
1727 : integer :: ISg,ngx,ik_bz,ikp_bz,dim_rtwg
1728 : integer :: neh1,neh2,ig,nblocks
1729 : integer :: ik_ibz,itim_k,ikp_ibz,itim_kp,isym_k,isym_kp
1730 : integer :: iq_bz,iq_ibz,isym_q,itim_q,iqbz0,rank
1731 : integer :: iv,ivp,ic,icp
1732 : integer :: block
1733 : integer(i8b) :: tot_nels,ir,it,itp
1734 : real(dp) :: faq,kx_fact
1735 : complex(sp) :: ctemp
1736 : character(len=500) :: msg
1737 : !arrays
1738 : integer :: bidx(2,4),spin_ids(2,3), my_cols(2),my_rows(2) !,proc_end(2),proc_start(2)
1739 : integer(i8b) :: nels_block(3)
1740 : integer,allocatable :: ncols_of(:), col_start(:),col_stop(:)
1741 : real(dp) :: qbz(3),tsec(2) !kbz(3),kpbz(3),
1742 : complex(dp),allocatable :: my_kxssp(:,:)
1743 : complex(gwp),allocatable :: vc_sqrt_qbz(:),rhotwg1(:),rhotwg2(:)
1744 : !************************************************************************
1745 :
1746 : DBG_ENTER("COLL")
1747 :
1748 : write(msg,'(a,2i2,a)')" Calculating exchange term for (spin1,spin2) ",spin1,spin2," ..."
1749 : call wrtout(std_out, msg)
1750 :
1751 : ! Basic constants.
1752 : dim_rtwg=1; faq = one/(Cryst%ucvol*Kmesh%nbz)
1753 :
1754 : ! Identify the index of q==0
1755 : iqbz0=0
1756 : do iq_bz=1,Qmesh%nbz
1757 : if (ALL(ABS(Qmesh%bz(:,iq_bz))<tol3)) iqbz0 = iq_bz
1758 : end do
1759 : ABI_CHECK(iqbz0/=0,"q=0 not found")
1760 : !
1761 : ! Treat the spin polarization.
1762 : spin_ids(:,1) = (/1,1/)
1763 : spin_ids(:,2) = (/2,2/)
1764 : spin_ids(:,3) = (/1,2/)
1765 :
1766 : nblocks=1
1767 : kx_fact=two
1768 : nels_block(:)=0
1769 : nels_block(1)=BSp%nreh(1)*(BSp%nreh(1)+1_i8b)/2
1770 : tot_nels=nels_block(1)
1771 :
1772 : if (nsppol==2) then
1773 : nblocks=3
1774 : kx_fact=one
1775 : nels_block(1) = BSp%nreh(1)*(BSp%nreh(1)+1_i8b)/2 ! Only the upper triangle for block 1 and 2
1776 : nels_block(2) = BSp%nreh(2)*(BSp%nreh(2)+1_i8b)/2
1777 : nels_block(3) = BSp%nreh(1)*BSp%nreh(2)*1_i8b ! Note: Block 3 does not have symmetries.
1778 : tot_nels= SUM(nels_block)
1779 : end if
1780 : !
1781 : ! Distribute the calculation of the matrix elements among the nodes.
1782 : ! * tstart and t_stop give the initial and final transition index treated by each node.
1783 : ! * my_hsize is the number of transitions treated by this processor
1784 : ! * my_cols(1:2) gives the initial and final column treated by this node.
1785 : !
1786 : do block=1,nsppol
1787 : !
1788 : ! Indices used to loop over bands.
1789 : ! bidx contains the starting and final indices used to loop over bands.
1790 : !
1791 : ! (b3,b4)
1792 : ! |... ...|
1793 : ! (b1,b2) |... ...|
1794 : !
1795 : ! Resonant matrix is given by
1796 : ! (v',c')
1797 : ! |... ...|
1798 : ! (v,c) |... ...|
1799 : !
1800 : ! Coupling matrix is given by
1801 : ! (c',v')
1802 : ! |... ...|
1803 : ! (v,c) |... ...|
1804 :
1805 : if (is_resonant) then
1806 : bidx(:,1) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b1
1807 : bidx(:,2) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b2
1808 : bidx(:,3) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b3
1809 : bidx(:,4) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b4
1810 : else
1811 : bidx(:,1) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b1
1812 : bidx(:,2) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b2
1813 : bidx(:,3) = [BSp%lumo_spin(block),BSp%humo_spin(block)] ! range for b3
1814 : bidx(:,4) = [BSp%lomo_spin(block),BSp%homo_spin(block)] ! range for b4
1815 : end if
1816 :
1817 : !spin1 = spin_ids(1,block)
1818 : !spin2 = spin_ids(2,block)
1819 :
1820 : my_cols=0
1821 : do itp=1,Bsp%nreh(block)
1822 : do it=1,itp
1823 : ir = it + itp*(itp-1_i8b)/2
1824 : if (ir==t_start(my_rank)) then
1825 : my_rows(1) = it
1826 : my_cols(1) = itp
1827 : end if
1828 : if (ir==t_stop(my_rank)) then
1829 : my_rows(2) = it
1830 : my_cols(2) = itp
1831 : end if
1832 : end do
1833 : end do
1834 :
1835 : ! Allocate big (scalable) buffer to store the BS matrix on this node.
1836 : !ABI_MALLOC(my_bsham,(t_start(my_rank):t_stop(my_rank)))
1837 : !
1838 : ! =====================
1839 : ! === Exchange term ===
1840 : ! =====================
1841 : ! TODO might add treatment of <psi|q+G|psi> for q+G -> 0
1842 : ! TODO might used enlarged G-sphere for better convergence.
1843 : !if (do_exchange_term) then
1844 : call timab(683,1,tsec) ! exc_build_ham(exchange)
1845 :
1846 : ABI_MALLOC(rhotwg1,(npweps))
1847 : ABI_MALLOC(rhotwg2,(npweps))
1848 :
1849 : ngx = Gsph_x%ng
1850 : ABI_MALLOC(vc_sqrt_qbz,(ngx))
1851 :
1852 : ! * Get iq_ibz, and symmetries from iq_bz.
1853 : iq_bz = iqbz0 ! q = 0 -> iqbz0
1854 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
1855 :
1856 : ! * Set up table of |q(BZ)+G|
1857 : if (iq_ibz==1) then
1858 : do ig=1,ngx
1859 : ISg = Gsph_x%rottb(ig,itim_q,isym_q)
1860 : vc_sqrt_qbz(ISg)=Vcp%vcqlwl_sqrt(ig,1)
1861 : end do
1862 : else
1863 : ABI_ERROR("iq_ibz should be 1")
1864 : end if
1865 :
1866 : do itp=1,BSp%nreh(block) ! Loop over transition tp = (kp,vp,cp,spin2)
1867 :
1868 : if (itp<my_cols(1) .or. itp>my_cols(2)) CYCLE ! I dont have this column.
1869 : ikp_bz = Bsp%Trans(itp,spin2)%k
1870 : ivp = Bsp%Trans(itp,spin2)%v
1871 : icp = Bsp%Trans(itp,spin2)%c
1872 :
1873 : ikp_ibz = Kmesh%tab (ikp_bz)
1874 : isym_kp = Kmesh%tabo(ikp_bz)
1875 : itim_kp = (3-Kmesh%tabi(ikp_bz))/2
1876 :
1877 : if (is_resonant) then
1878 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,ivp,icp,ikp_ibz,spin2),Gsph_c)
1879 : else ! Code for coupling block.
1880 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,icp,ivp,ikp_ibz,spin2),Gsph_c)
1881 : end if
1882 : !
1883 : ! Multiply by the Coulomb term.
1884 : do ig=2,npweps
1885 : rhotwg2(ig) = rhotwg2(ig) * vc_sqrt_qbz(ig) * vc_sqrt_qbz(ig)
1886 : end do
1887 :
1888 : do it=1,itp ! Loop over transition t = (k,v,c,spin1)
1889 : ir = it + itp*(itp-1_i8b)/2
1890 : if (ir<t_start(my_rank) .or. ir>t_stop(my_rank)) CYCLE
1891 :
1892 : ik_bz = Bsp%Trans(it,spin1)%k
1893 : iv = Bsp%Trans(it,spin1)%v
1894 : ic = Bsp%Trans(it,spin1)%c
1895 :
1896 : ik_ibz = Kmesh%tab(ik_bz)
1897 : isym_k = Kmesh%tabo(ik_bz)
1898 : itim_k = (3-Kmesh%tabi(ik_bz))/2
1899 : !if (itim_k==2) CYCLE ! time-reversal or not
1900 :
1901 : rhotwg1(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,iv,ic,ik_ibz,spin1),Gsph_c)
1902 : !
1903 : ! sum over G/=0
1904 : ctemp = xdotc(npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
1905 : ctemp = faq * kx_fact * ctemp
1906 :
1907 : ! exchange term is non divergent !
1908 : !if (BSp%prep_interp) then
1909 : ! ccoeffs(ir) = ccoeffs(ir) + ctemp
1910 : !end if
1911 :
1912 : my_bsham(ir) = my_bsham(ir) + ctemp
1913 : end do !it
1914 : end do !itp
1915 :
1916 : ABI_FREE(rhotwg1)
1917 : ABI_FREE(rhotwg2)
1918 : ABI_FREE(vc_sqrt_qbz)
1919 :
1920 : call timab(683,2,tsec) ! exc_build_ham(exchange)
1921 : !end if ! do_exchange_term
1922 : end do ! block
1923 :
1924 : !
1925 : ! ===========================================
1926 : ! === Exchange term for spin_up spin_down ===
1927 : ! ===========================================
1928 :
1929 : if (nsppol==2) then
1930 : call timab(686,2,tsec) ! exc_build_ham(exch.spin)
1931 : block=3
1932 : neh1=BSp%nreh(1)
1933 : neh2=BSp%nreh(2)
1934 : !
1935 : ! The oscillators at q=0 are available on each node for both spin.
1936 : ! Here the calculation of the block is parallelized over columns.
1937 : ABI_MALLOC(col_start,(0:nproc-1))
1938 : ABI_MALLOC(col_stop,(0:nproc-1))
1939 : call xmpi_split_work2_i4b(neh2,nproc,col_start,col_stop) !check this but it should be OK.
1940 :
1941 : my_cols(1) = col_start(my_rank)
1942 : my_cols(2) = col_stop (my_rank)
1943 : if (my_cols(2)-my_cols(1)<=0) then
1944 : ABI_ERROR("One of the processors has zero columns!")
1945 : end if
1946 :
1947 : ABI_MALLOC(ncols_of,(0:nproc-1))
1948 : ncols_of=0
1949 : do rank=0,nproc-1
1950 : if (col_stop(rank)>=col_start(rank)) ncols_of(rank) = col_stop(rank)-col_start(rank)+1
1951 : end do
1952 :
1953 : ABI_FREE(col_start)
1954 : ABI_FREE(col_stop)
1955 : !
1956 : ! TODO might add treatment of <psi|q+G|psi> for q+G -> 0
1957 : ! TODO might used enlarged G-sphere for better convergence.
1958 : ! Note that my_kxssp is always written on file when nsppol=2, even when
1959 : ! non-local field effects are neglected.
1960 : ABI_MALLOC(my_kxssp,(neh1,my_cols(1):my_cols(2)))
1961 : my_kxssp=czero
1962 :
1963 : !if (do_exchange_term) then
1964 : !spin1=1; spin2=2
1965 : ABI_MALLOC(rhotwg1,(npweps))
1966 : ABI_MALLOC(rhotwg2,(npweps))
1967 :
1968 : ngx = Gsph_x%ng
1969 : ABI_MALLOC(vc_sqrt_qbz,(ngx))
1970 : !
1971 : ! * Get iq_ibz, and symmetries from iq_bz.
1972 : iq_bz = iqbz0 ! q = 0 -> iqbz0
1973 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
1974 : !
1975 : ! * Set up table of |q(BZ)+G|
1976 : if (iq_ibz==1) then
1977 : do ig=1,ngx
1978 : ISg = Gsph_x%rottb(ig,itim_q,isym_q)
1979 : vc_sqrt_qbz(ISg)=Vcp%vcqlwl_sqrt(ig,1)
1980 : end do
1981 : else
1982 : ABI_ERROR("iq_ibz should be 1")
1983 : end if
1984 :
1985 : do itp=1,neh2 ! Loop over transition tp = (kp,vp,cp,spin2)
1986 :
1987 : if (itp<my_cols(1) .or. itp>my_cols(2)) CYCLE ! I dont have this column.
1988 : ikp_bz = Bsp%Trans(itp,spin2)%k
1989 : ivp = Bsp%Trans(itp,spin2)%v
1990 : icp = Bsp%Trans(itp,spin2)%c
1991 :
1992 : ikp_ibz = Kmesh%tab (ikp_bz)
1993 : isym_kp = Kmesh%tabo(ikp_bz)
1994 : itim_kp = (3-Kmesh%tabi(ikp_bz))/2
1995 :
1996 : if (is_resonant) then
1997 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,ivp,icp,ikp_ibz,spin2),Gsph_c)
1998 : else ! Code for coupling block.
1999 : rhotwg2(:) = sym_rhotwgq0(itim_kp,isym_kp,dim_rtwg,npweps,rhxtwg_q0(:,icp,ivp,ikp_ibz,spin2),Gsph_c)
2000 : end if
2001 : !
2002 : ! Multiply by the Coulomb term.
2003 : do ig=2,npweps
2004 : rhotwg2(ig) = rhotwg2(ig) * vc_sqrt_qbz(ig) * vc_sqrt_qbz(ig)
2005 : end do
2006 :
2007 : do it=1,neh1 ! Loop over transition t = (k,v,c,spin1) FULL matrix.
2008 : ik_bz = Bsp%Trans(it,spin1)%k
2009 : iv = Bsp%Trans(it,spin1)%v
2010 : ic = Bsp%Trans(it,spin1)%c
2011 :
2012 : ik_ibz = Kmesh%tab(ik_bz)
2013 : isym_k = Kmesh%tabo(ik_bz)
2014 : itim_k = (3-Kmesh%tabi(ik_bz))/2
2015 : !if (itim_k==2) CYCLE ! time-reversal or not
2016 :
2017 : rhotwg1(:) = sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npweps,rhxtwg_q0(:,iv,ic,ik_ibz,spin1),Gsph_c)
2018 : !
2019 : ! sum over G/=0
2020 : ctemp = XDOTC(npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
2021 : ctemp = faq * kx_fact * ctemp
2022 :
2023 : my_kxssp(it,itp) = ctemp
2024 : end do !it
2025 : end do !itp
2026 :
2027 : ABI_FREE(rhotwg1)
2028 : ABI_FREE(rhotwg2)
2029 : ABI_FREE(vc_sqrt_qbz)
2030 : !end if ! do_exchange_term
2031 : call timab(686,2,tsec) ! exc_build_ham(exch.spin)
2032 :
2033 : ABI_FREE(ncols_of)
2034 : ABI_SFREE(my_kxssp)
2035 : end if
2036 :
2037 : DBG_EXIT("COLL")
2038 :
2039 : end subroutine exc_build_v
2040 : !!***
2041 :
2042 : !!****f* m_exc_build/exc_build_ham
2043 : !! NAME
2044 : !! exc_build_ham
2045 : !!
2046 : !! FUNCTION
2047 : !! Calculate and write the excitonic Hamiltonian on an external binary file (Fortran binary file)
2048 : !! for subsequent treatment in the Bethe-Salpeter code.
2049 : !!
2050 : !! INPUTS
2051 : !! BSp<excparam>=The parameters for the Bethe-Salpeter calculation.
2052 : !! BS_files<excfiles>=File names internally used in the BS code.
2053 : !! Cryst<crystal_t>=Info on the crystalline structure.
2054 : !! Kmesh<kmesh_t>=The list of k-points in the BZ, IBZ and symmetry tables.
2055 : !! Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
2056 : !! ktabr(nfftot_osc,BSp%nkbz)=The FFT index of $(R^{-1}(r-\tau))$ where R is symmetry needed to obtains
2057 : !! the k-points from the irreducible image. Used to symmetrize u_Sk where S = \transpose R^{-1}
2058 : !! Gsph_x<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored).
2059 : !! Gsph_c<gsphere_t>=Info on the G-sphere used to describe the correlation part.
2060 : !! Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
2061 : !! screen<screen_t>=Data type gathering info and data for W.
2062 : !! nfftot_osc=Total Number of FFT points used for the oscillator matrix elements.
2063 : !! ngfft_osc(18)=Info on the FFT algorithm used to calculate the oscillator matrix elements.
2064 : !! Psps<Pseudopotential_type>=Variables related to pseudopotentials
2065 : !! Pawtab(Psps%ntypat)<pawtab_type>=PAW tabulated starting data.
2066 : !! Pawang<pawang_type>=PAW angular mesh and related data.
2067 : !! Paw_pwff(Cryst%ntypat*Wfd%usepaw)<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
2068 : !! Wfd<wfdgw_t>=Handler for the wavefunctions.
2069 : !!
2070 : !! OUTPUT
2071 : !! The excitonic Hamiltonian is saved on an external binary file (see below).
2072 : !!
2073 : !! SOURCE
2074 :
2075 29 : subroutine exc_build_ham(BSp,BS_files,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
2076 29 : Wfd,screen,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff)
2077 :
2078 : !Arguments ------------------------------------
2079 : !scalars
2080 : integer,intent(in) :: nfftot_osc
2081 : type(excparam),intent(in) :: BSp
2082 : type(excfiles),intent(in) :: BS_files
2083 : type(screen_t),intent(inout) :: screen
2084 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
2085 : type(crystal_t),intent(in) :: Cryst
2086 : type(vcoul_t),intent(in) :: Vcp
2087 : type(gsphere_t),intent(in) :: Gsph_x,Gsph_c
2088 : type(Pseudopotential_type),intent(in) :: Psps
2089 : type(Hdr_type),intent(inout) :: Hdr_bse
2090 : type(pawang_type),intent(in) :: Pawang
2091 : type(wfdgw_t),target,intent(inout) :: Wfd
2092 : !arrays
2093 : integer,intent(in) :: ngfft_osc(18)
2094 : integer,intent(in) :: ktabr(nfftot_osc,Kmesh%nbz)
2095 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Wfd%usepaw)
2096 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Wfd%usepaw)
2097 :
2098 : !Local variables ------------------------------
2099 : !scalars
2100 : logical :: do_resonant,do_coupling
2101 : !character(len=500) :: msg
2102 : !arrays
2103 : real(dp) :: tsec(2)
2104 29 : complex(gwp),allocatable :: all_mgq0(:,:,:,:,:)
2105 : !************************************************************************
2106 :
2107 29 : call timab(670,1,tsec)
2108 :
2109 29 : ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
2110 116 : ABI_CHECK(nfftot_osc==PRODUCT(ngfft_osc(1:3)),"mismatch in FFT size")
2111 :
2112 29 : if (BSp%have_complex_ene) then
2113 0 : ABI_ERROR("Complex energies are not supported yet")
2114 : end if
2115 :
2116 : ! Do we have to compute some block?
2117 29 : do_resonant = (BS_files%in_hreso == BSE_NOFILE)
2118 29 : do_coupling = (BS_files%in_hcoup == BSE_NOFILE)
2119 :
2120 29 : if (BSp%use_coupling == 0) then
2121 28 : if (.not.do_resonant) then
2122 6 : call wrtout(std_out,"Will skip the calculation of resonant block (will use BSR file)")
2123 6 : goto 100
2124 : end if
2125 : else
2126 1 : if (.not. do_resonant .and. .not. do_coupling) then
2127 0 : call wrtout(std_out,"Will skip the calculation of both resonant and coupling block (will use BSR and BSC files)")
2128 0 : goto 100
2129 : end if
2130 : end if
2131 :
2132 : ! Compute M_{k,q=0}^{b,b}(G) for all k-points in the IBZ and each pair b, b'
2133 : ! used for the exchange part and part of the Coulomb term.
2134 23 : call wrtout(std_out," Calculating all matrix elements for q=0 to save CPU time")
2135 :
2136 : call wfd_all_mgq0(Wfd,Cryst,Qmesh,Gsph_x,Vcp,Psps,Pawtab,Paw_pwff,&
2137 23 : Bsp%lomo_spin,Bsp%homo_spin,Bsp%humo_spin,nfftot_osc,ngfft_osc,Bsp%npweps,all_mgq0)
2138 :
2139 : ! ========================
2140 : ! ==== Resonant Block ====
2141 : ! ========================
2142 23 : if (do_resonant) then
2143 22 : call timab(672,1,tsec)
2144 : call exc_build_block(BSp,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
2145 22 : Wfd,screen,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff,all_mgq0,.TRUE.,BS_files%out_hreso)
2146 22 : call timab(672,2,tsec)
2147 : end if
2148 :
2149 : ! ========================
2150 : ! ==== Coupling Block ====
2151 : ! ========================
2152 23 : if (do_coupling .and. BSp%use_coupling > 0) then
2153 1 : call timab(673,1,tsec)
2154 : call exc_build_block(BSp,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
2155 1 : Wfd,screen,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff,all_mgq0,.FALSE.,BS_files%out_hcoup)
2156 1 : call timab(673,2,tsec)
2157 : end if
2158 :
2159 : ! Free memory.
2160 23 : ABI_FREE(all_mgq0)
2161 :
2162 29 : 100 call timab(670,2,tsec)
2163 :
2164 29 : end subroutine exc_build_ham
2165 : !!***
2166 :
2167 : !!****f* m_exc_build/wfd_all_mgq0
2168 : !! NAME
2169 : !! wfd_all_mgq0
2170 : !!
2171 : !! FUNCTION
2172 : !!
2173 : !! INPUTS
2174 : !! Wfd<wfdgw_t>=Handler for the wavefunctions.
2175 : !! Cryst<crystal_t>=Info on the crystalline structure.
2176 : !! Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
2177 : !! Gsph_x<gsphere_t>=G-sphere with the G-vectors in mgq0.
2178 : !! Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
2179 : !! Psps<Pseudopotential_type>=Variables related to pseudopotentials
2180 : !! Pawtab(Psps%ntypat)<pawtab_type>=PAW tabulated starting data.
2181 : !! Paw_pwff(Cryst%ntypat*Wfd%usepaw)<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
2182 : !! lomo_spin(Wfd%nsppol)=Lowest occupied band for each spin
2183 : !! homo_spin(Wfd%nsppol)=Highest occupied band for each spin
2184 : !! humo_spin(Wfd%nsppol)=Highest unoccupied band for each spin
2185 : !! nfftot_osc=Total Number of FFT points used for the oscillator matrix elements.
2186 : !! ngfft_osc(18)=Info on the FFT algorithm used to calculate the oscillator matrix elements.
2187 : !! npweps=Number of G-vectors in mgq0.
2188 : !!
2189 : !! OUTPUT
2190 : !! mgq0(npweps,lomo_min:humo_max,lomo_min:humo_max,Wfd%nkibz,Wfd%nsppol)
2191 : !! Allocated here and filled with the matrix elements on each node.
2192 : !!
2193 : !! SOURCE
2194 :
2195 23 : subroutine wfd_all_mgq0(Wfd,Cryst,Qmesh,Gsph_x,Vcp,&
2196 23 : & Psps,Pawtab,Paw_pwff,lomo_spin,homo_spin,humo_spin,nfftot_osc,ngfft_osc,npweps,mgq0)
2197 :
2198 : !Arguments ------------------------------------
2199 : !scalars
2200 : integer,intent(in) :: nfftot_osc,npweps
2201 : type(kmesh_t),intent(in) :: Qmesh
2202 : type(crystal_t),intent(in) :: Cryst
2203 : type(vcoul_t),intent(in) :: Vcp
2204 : type(gsphere_t),intent(in) :: Gsph_x
2205 : type(Pseudopotential_type),intent(in) :: Psps
2206 : type(wfdgw_t),target,intent(inout) :: Wfd
2207 : !arrays
2208 : integer,intent(in) :: lomo_spin(Wfd%nsppol),homo_spin(Wfd%nsppol),humo_spin(Wfd%nsppol)
2209 : integer,intent(in) :: ngfft_osc(18)
2210 : complex(gwp),allocatable,intent(out) :: mgq0(:,:,:,:,:)
2211 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
2212 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Wfd%usepaw)
2213 :
2214 : !Local variables ------------------------------
2215 : !scalars
2216 : integer,parameter :: map2sphere1=1,dim_rtwg1=1,ndat1=1
2217 : integer :: use_padfft,mgfft_osc,fftalga_osc,ii
2218 : integer :: ik_ibz,itim_k,isym_k,iq_bz,iq_ibz,isym_q,itim_q,iqbz0
2219 : integer :: ierr,iv,ic,spin,lomo_min,humo_max !,inv_ipw,ipw
2220 : real(dp) :: cpu,wall,gflops !q0vol,fcc_const
2221 : complex(dp) :: ph_mkt
2222 : character(len=500) :: msg
2223 : type(wave_t),pointer :: wave_v, wave_c
2224 : !arrays
2225 23 : integer,allocatable :: igfftg0(:),task_distrib(:,:,:,:)
2226 23 : integer,allocatable :: gbound(:,:),id_tab(:)
2227 : real(dp) :: qbz(3),spinrot_k(4),tsec(2)
2228 23 : complex(gwp),allocatable :: rhotwg1(:)
2229 23 : complex(gwp),target,allocatable :: ur1(:),ur2(:)
2230 23 : complex(gwp),contiguous, pointer :: ptr_ur1(:),ptr_ur2(:)
2231 23 : type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:)
2232 23 : type(pawpwij_t),allocatable :: Pwij_q0(:)
2233 : !************************************************************************
2234 :
2235 23 : call timab(671,1,tsec)
2236 :
2237 23 : ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
2238 92 : ABI_CHECK(nfftot_osc==PRODUCT(ngfft_osc(1:3)),"mismatch in FFT size")
2239 :
2240 71 : lomo_min = MINVAL(lomo_spin); humo_max = MAXVAL(humo_spin)
2241 :
2242 50 : if ( ANY(ngfft_osc(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_osc)
2243 :
2244 92 : mgfft_osc = MAXVAL(ngfft_osc(1:3))
2245 23 : fftalga_osc = ngfft_osc(7)/100 !; fftalgc_osc=MOD(ngfft_osc(7),10)
2246 :
2247 : ! (temporary) Table used for the wavefunction in the IBZ.
2248 69 : ABI_MALLOC(id_tab, (Wfd%nfft))
2249 482671 : id_tab = (/(ii, ii=1,Wfd%nfft)/)
2250 :
2251 : ! Analytic integration of 4pi/q^2 over the volume element:
2252 : ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
2253 : ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
2254 : ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255
2255 : ! (see gwa.pdf, appendix A.4)
2256 :
2257 : ! If q=0 and C=V then set up rho-twiddle(G=0) to reflect an
2258 : ! analytic integration of q**-2 over the volume element:
2259 : ! <q**-2> = 7.44 V**(-2/3) (for fcc cell)
2260 :
2261 : ! q0vol = (8.0*pi**3) / (Cryst%ucvol*Kmesh%nbz)
2262 : ! fcc_const = SQRT(7.44*q0vol**(-2.0/3.0))
2263 : ! rtw = (6.0*pi**2/(Cryst%ucvol*Kmesh%nkbz))**(1./3.)
2264 : ! Average of (q+q')**-2 integration for head of Coulomb matrix
2265 : ! INTRTW(QL) = (2*pi*rtw + pi*(rtw**2/QL-QL)*LOG((QL+rtw)/(QL-rtw)))
2266 : ! & * (Cryst%ucvol*Kmesh%nbz)/(2*pi)**3. * QL*QL
2267 :
2268 23 : if (Wfd%usepaw==1) then
2269 14 : ABI_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor))
2270 2 : call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm)
2271 14 : ABI_MALLOC(Cp2,(Wfd%natom,Wfd%nspinor))
2272 2 : call pawcprj_alloc(Cp2,0,Wfd%nlmn_atm)
2273 : end if
2274 :
2275 69 : ABI_MALLOC(ur1,(nfftot_osc*Wfd%nspinor))
2276 46 : ABI_MALLOC(ur2,(nfftot_osc*Wfd%nspinor))
2277 :
2278 : ! Identify q==0
2279 23 : iqbz0=0
2280 1079 : do iq_bz=1,Qmesh%nbz
2281 1446 : if (ALL(ABS(Qmesh%bz(:,iq_bz))<tol3)) iqbz0 = iq_bz
2282 : end do
2283 23 : ABI_CHECK(iqbz0/=0,"q=0 not found in q-point list!")
2284 :
2285 : ! * Get iq_ibz, and symmetries from iqbz0.
2286 23 : call qmesh%get_BZ_item(iqbz0,qbz,iq_ibz,isym_q,itim_q)
2287 :
2288 23 : if (Wfd%usepaw==1) then ! Prepare onsite contributions at q==0
2289 9 : ABI_MALLOC(Pwij_q0,(Cryst%ntypat))
2290 2 : call pawpwij_init(Pwij_q0,npweps,Qmesh%bz(:,iqbz0),Gsph_x%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
2291 : end if
2292 : !
2293 : ! Tables for the FFT of the oscillators.
2294 : ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
2295 : ! b) gbound table for the zero-padded FFT performed in rhotwg.
2296 69 : ABI_MALLOC(igfftg0,(Gsph_x%ng))
2297 92 : ABI_MALLOC(gbound,(2*mgfft_osc+8,2))
2298 23 : call Gsph_x%fft_tabs((/0,0,0/),mgfft_osc,ngfft_osc,use_padfft,gbound,igfftg0)
2299 23 : if ( ANY(fftalga_osc == (/2,4/)) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
2300 23 : if (use_padfft==0) then
2301 0 : ABI_FREE(gbound)
2302 0 : ABI_MALLOC(gbound,(2*mgfft_osc+8,2*use_padfft))
2303 : end if
2304 :
2305 69 : ABI_MALLOC(rhotwg1,(npweps))
2306 :
2307 161 : ABI_MALLOC_OR_DIE(mgq0, (npweps,lomo_min:humo_max,lomo_min:humo_max,Wfd%nkibz,Wfd%nsppol), ierr)
2308 1317080 : mgq0 = czero
2309 :
2310 23 : call cwtime(cpu,wall,gflops,"start")
2311 :
2312 47 : do spin=1,Wfd%nsppol
2313 : ! Distribute the calculation of the matrix elements.
2314 : ! processors have the entire set of wavefunctions hence we divide the workload
2315 : ! without checking if the pair of states is available. Last dimension is fake.
2316 120 : ABI_MALLOC(task_distrib,(lomo_spin(spin):humo_spin(spin),lomo_spin(spin):humo_spin(spin),Wfd%nkibz,1))
2317 24 : call xmpi_distab(Wfd%nproc,task_distrib)
2318 :
2319 : ! loop over the k-points in IBZ
2320 887 : do ik_ibz=1,Wfd%nkibz
2321 8543 : if ( ALL(task_distrib(:,:,ik_ibz,1)/= Wfd%my_rank) ) CYCLE
2322 :
2323 : ! Don't need to symmetrize the wavefunctions.
2324 3035 : itim_k=1; isym_k=1; ph_mkt=cone; spinrot_k=Cryst%spinrot(:,isym_k)
2325 :
2326 4109 : do iv=lomo_spin(spin),humo_spin(spin) ! Loop over band V
2327 3478 : if ( ALL(task_distrib(:,iv,ik_ibz,1)/=Wfd%my_rank) ) CYCLE
2328 :
2329 3478 : ABI_CHECK(wfd%get_wave_ptr(iv, ik_ibz, spin, wave_v, msg) == 0, msg)
2330 :
2331 3478 : if (wave_v%has_ur == WFD_STORED) then
2332 2989 : ptr_ur1 => wave_v%ur
2333 : else
2334 489 : call wfd%get_ur(iv,ik_ibz,spin,ur1)
2335 489 : ptr_ur1 => ur1
2336 : end if
2337 :
2338 3478 : if (Wfd%usepaw==1) call wfd%get_cprj(iv,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.)
2339 :
2340 : ! Loop over band C
2341 24829 : do ic=lomo_spin(spin),humo_spin(spin)
2342 20488 : if ( task_distrib(ic,iv,ik_ibz,1)/=Wfd%my_rank ) CYCLE
2343 :
2344 20488 : ABI_CHECK(wfd%get_wave_ptr(ic, ik_ibz, spin, wave_c, msg) == 0, msg)
2345 :
2346 20488 : if (wave_c%has_ur == WFD_STORED) then
2347 17853 : ptr_ur2 => wave_c%ur
2348 : else
2349 2635 : call wfd%get_ur(ic,ik_ibz,spin,ur2)
2350 2635 : ptr_ur2 => ur2
2351 : end if
2352 :
2353 20488 : if (Wfd%usepaw==1) call wfd%get_cprj(ic,ik_ibz,spin,Cryst,Cp2,sorted=.FALSE.)
2354 :
2355 : call rho_tw_g(Wfd%nspinor,npweps,nfftot_osc,ndat1,ngfft_osc,map2sphere1,use_padfft,igfftg0,gbound,&
2356 : ptr_ur1,1,id_tab,ph_mkt,spinrot_k,&
2357 : ptr_ur2,1,id_tab,ph_mkt,spinrot_k,&
2358 20488 : dim_rtwg1,rhotwg1)
2359 :
2360 20488 : if (Wfd%usepaw==1) then
2361 : ! Add PAW onsite contribution.
2362 304 : call paw_rho_tw_g(cryst,Pwij_q0,npweps,dim_rtwg1,Wfd%nspinor,Gsph_x%gvec,Cp1,Cp2,rhotwg1)
2363 : end if
2364 :
2365 : ! If q=0 treat Exchange and Coulomb-term independently
2366 20488 : if (iv <= homo_spin(spin) .and. ic <= homo_spin(spin) .or. &
2367 : iv > homo_spin(spin) .and. ic > homo_spin(spin)) then
2368 :
2369 10424 : if (iv/=ic) then !COULOMB term: C/=V: ignore them
2370 6946 : rhotwg1(1) = czero_gw
2371 : else
2372 : ! If q=0 and C=V then set up rho-twiddle(G=0) to reflect an
2373 : ! analytic integration of q**-2 over the volume element:
2374 : ! <q**-2> = 7.44 V**(-2/3) (for fcc cell)
2375 : !rhotwg1(1) = fcc_const * qpg(1,iqbz0)
2376 3478 : rhotwg1(1) = SQRT(GWPC_CMPLX(Vcp%i_sz,zero)) / Vcp%vcqlwl_sqrt(1,1)
2377 : !if (vcut) rhotwg1(1) = 1.0
2378 : end if
2379 :
2380 : else
2381 : ! At present this term is set to zero
2382 : ! EXCHANGE term: limit value.
2383 : ! Set up rho-twiddle(G=0) using small vector q instead of zero and k.p perturbation theory (see notes)
2384 10064 : rhotwg1(1) = czero_gw
2385 : end if
2386 :
2387 930890 : mgq0(:,iv,ic,ik_ibz,spin) = rhotwg1(:)
2388 : end do !ic
2389 : end do !iv
2390 : end do !ik_ibz
2391 :
2392 47 : ABI_FREE(task_distrib)
2393 : end do !spin
2394 :
2395 : ! TODO: One can speedup the calculation by computing the upper triangle of the
2396 : ! matrix in (b,b') space and then take advantage of the symmetry property:
2397 : !
2398 : ! M_{k,0}{{bb'}(G)^* = M{k,0}{b'b'}(-G)
2399 :
2400 : #if 0
2401 : !!!! $OMP PARALLEL DO COLLAPSE(3) PRIVATE(inv_ipw)
2402 : do spin=1,Wfd%nsppol
2403 : do ik_ibz=1,Wfd%nkibz
2404 : do iv=lomo_spin(spin),humo_spin(spin)
2405 : do ic=1,iv-1
2406 : do ipw=1,npweps
2407 : inv_ipw = gsph_x%g2mg(ipw)
2408 : mgq0(inv_ipw,ic,iv,ik_ibz,spin) = mgq0(ipw,iv,ic,ik_ibz,spin)
2409 : end do
2410 : end do
2411 : end do
2412 : end do
2413 : end do
2414 : #endif
2415 : !
2416 : ! Gather matrix elements on each node.
2417 23 : call xmpi_sum(mgq0,Wfd%comm,ierr)
2418 :
2419 23 : call cwtime(cpu,wall,gflops,"stop")
2420 23 : write(msg,'(2(a,f9.6))')"cpu_time = ",cpu,", wall_time = ",wall
2421 23 : call wrtout(std_out, msg)
2422 :
2423 23 : ABI_FREE(rhotwg1)
2424 23 : ABI_FREE(igfftg0)
2425 23 : ABI_FREE(gbound)
2426 23 : ABI_FREE(ur1)
2427 23 : ABI_FREE(ur2)
2428 23 : ABI_FREE(id_tab)
2429 :
2430 23 : if (Wfd%usepaw==1) then
2431 : ! Deallocation for PAW.
2432 2 : call pawpwij_free(Pwij_q0)
2433 5 : ABI_FREE(Pwij_q0)
2434 2 : call pawcprj_free(Cp1)
2435 6 : ABI_FREE(Cp1)
2436 2 : call pawcprj_free(Cp2)
2437 6 : ABI_FREE(Cp2)
2438 : end if
2439 :
2440 23 : call timab(671,2,tsec)
2441 :
2442 92 : end subroutine wfd_all_mgq0
2443 : !!***
2444 :
2445 : end module m_exc_build
2446 : !!***
|
