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