Line data Source code
1 : !!****m* ABINIT/m_vhxc_me
2 : !! NAME
3 : !! m_vhxc_me
4 : !!
5 : !! FUNCTION
6 : !! Evaluate the matrix elements of $v_H$ and $v_{xc}$ and $v_U$
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_vhxc_me
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xcdata
28 : use libxc_functionals
29 : use m_dtset
30 : use m_distribfft
31 :
32 : use defs_datatypes,only : pseudopotential_type
33 : use defs_abitypes, only : MPI_type
34 : use m_pawang, only : pawang_type
35 : use m_pawtab, only : pawtab_type
36 : use m_paw_an, only : paw_an_type
37 : use m_paw_ij, only : paw_ij_type
38 : use m_pawfgrtab, only : pawfgrtab_type
39 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
40 : use m_paw_denpot, only : paw_mknewh0
41 : use m_hide_blas, only : xdotc
42 : use m_wfd, only : wfdgw_t, wave_t
43 : use m_crystal, only : crystal_t
44 : use m_melemts, only : melflags_t, melements_t
45 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
46 : use m_kg, only : mkkin
47 : use m_rhotoxc, only : rhotoxc
48 :
49 : implicit none
50 :
51 : private
52 : !!***
53 :
54 : public :: calc_vhxc_me
55 : !!***
56 :
57 : contains
58 : !!***
59 :
60 : !!****f* ABINIT/calc_vhxc_me
61 : !! NAME
62 : !! calc_vhxc_me
63 : !!
64 : !! FUNCTION
65 : !! Evaluate the matrix elements of $v_H$ and $v_{xc}$ and $v_U$
66 : !! both in case of NC pseudopotentials and PAW (DFT+U, presently, is only available in PAW)
67 : !! The matrix elements of $v_{xc}$ are calculated with and without the core contribution.
68 : !! The later quantity is required in case of GW calculations.
69 : !!
70 : !! INPUTS
71 : !! Wfd <type (wfdgw_t)>=Structure gathering information on the wavefunctions.
72 : !! Mflags: Flags specifying the quantities to be computed.
73 : !! Dtset <type(dataset_type)>=all input variables in this dataset
74 : !! ngfftf(18)contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
75 : !! nfftf=number of points in the fine FFT mesh (for this processor)
76 : !! Pawtab(Dtset%ntypat*Dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
77 : !! Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
78 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
79 : !! Paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
80 : !! Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
81 : !! Cryst<crystal_t>=unit cell and symmetries
82 : !! vhartr(nfftf)= Hartree potential in real space on the fine FFT mesh
83 : !! vxc(nfftf,nspden)= xc potential in real space on the fine FFT grid
84 : !! rhor(nfftf,nspden)=density in real space (smooth part if PAW).
85 : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
86 : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
87 : !! kstab(2,Wfd%nkibz,Wfd%nsppol)=Table temporary used to be compatible with the old implementation.
88 : !!
89 : !! OUTPUT
90 : !! Mels
91 : !! %kinetic=matrix elements of $t$.
92 : !! %vxc =matrix elements of $v_{xc}[nv+nc]$.
93 : !! %vxcval =matrix elements of $v_{xc}[nv]$.
94 : !! %vxcval_hybrid=matrix elements of $v_{xc}[nv]^{hybrid functional}$.
95 : !! %vhartr =matrix elements of $v_H$.
96 : !! %vu =matrix elements of $v_U$.
97 : !!
98 : !! SIDE EFFECTS
99 : !! Paw_ij= In case of self-Consistency it is changed. It will contain the new H0
100 : !! Hamiltonian calculated using the QP densities. The valence contribution to XC
101 : !! is removed.
102 : !!
103 : !! NOTES
104 : !! All the quantities ($v_H$, $v_{xc}$ and $\psi$ are evaluated on the "fine" FFT mesh.
105 : !! In case of calculations with NC pseudopotentials the usual mesh is defined by ecut.
106 : !! For PAW calculations the dense FFT grid defined by pawecutdg is used
107 : !! Besides, in case of PAW, the matrix elements of V_hartree do not contain the onsite
108 : !! contributions due to the Coulomb potential generated by ncore and tncore.
109 : !! These quantities, as well as the onsite kinetic terms, are stored in Paw_ij%dij0.
110 : !!
111 : !! SOURCE
112 :
113 281 : subroutine calc_vhxc_me(Wfd, Mflags, Mels, Cryst, Dtset, nfftf, ngfftf, &
114 281 : vtrial, vhartr, vxc, Psps, Pawtab, Paw_an, Pawang, Pawfgrtab, Paw_ij, dijexc_core, &
115 281 : rhor, usexcnhat, nhat, nhatgr, nhatgrdim, kstab, &
116 281 : taur) ! optional arguments
117 :
118 : !Arguments ------------------------------------
119 : !scalars
120 : integer,intent(in) :: nhatgrdim,usexcnhat,nfftf
121 : type(Dataset_type),intent(in) :: Dtset
122 : type(Pseudopotential_type),intent(in) :: Psps
123 : type(wfdgw_t),target,intent(inout) :: Wfd
124 : type(Pawang_type),intent(in) :: Pawang
125 : type(crystal_t),intent(in) :: Cryst
126 : type(melflags_t),intent(in) :: Mflags
127 : type(melements_t),intent(out) :: Mels
128 : !arrays
129 : integer,intent(in) :: ngfftf(18)
130 : integer,intent(in) :: kstab(2, Wfd%nkibz, Wfd%nsppol)
131 : real(dp),intent(in) :: vhartr(nfftf), vxc(nfftf, Wfd%nspden), vtrial(nfftf, Wfd%nspden)
132 : real(dp),intent(in) :: rhor(nfftf, Wfd%nspden)
133 : real(dp),intent(in) :: nhat(nfftf, Wfd%nspden * Wfd%usepaw)
134 : real(dp),intent(in) :: nhatgr(nfftf, Wfd%nspden, 3 * nhatgrdim)
135 : real(dp),intent(in),optional :: taur(nfftf, Wfd%nspden * Dtset%usekden)
136 : real(dp),intent(in) :: dijexc_core(:,:,:) ! (cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)
137 : type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat * Wfd%usepaw)
138 : type(Paw_an_type),intent(in) :: Paw_an(Cryst%natom)
139 : type(Paw_ij_type),intent(inout) :: Paw_ij(Cryst%natom)
140 : type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom)
141 :
142 : !Local variables-------------------------------
143 : !scalars
144 : integer :: auxc_ixc,iat,ikc,ik_ibz,ib,jb,is,b_start,b_stop,istwf_k
145 : integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,lmn2_size_max
146 : integer :: isppol,cplex_dij,npw_k,nspinor,nsppol,nspden,nk_calc,rank
147 : integer :: iab,isp1,isp2,ixc_sigma,nsploop,nkxc,option,n3xccc_,nk3xc,my_nbbp,my_nmels
148 : real(dp) :: nfftfm1,fact,DijH,bigexc_val,bigsxc_val,bigexc_hybrid_val,vxcval_avg,vxcval_hybrid_avg,h0dij,vxc1,vxc1_val,re_p,im_p,dijsigcx,dum
149 : complex(dp) :: cdot
150 : logical :: ltest,nmxc
151 : character(len=500) :: msg
152 281 : type(MPI_type) :: MPI_enreg_seq
153 : type(xcdata_type) :: xcdata,xcdata_hybrid
154 : type(wave_t),pointer :: wave_jb, wave_ib
155 : !arrays
156 : integer,parameter :: spinor_idxs(2,4)=RESHAPE([1,1,2,2,1,2,2,1], [2,4])
157 562 : integer :: got(Wfd%nproc)
158 281 : integer,allocatable :: kcalc2ibz(:),dimlmn(:),bbp_ks_distrb(:,:,:,:)
159 281 : integer,contiguous, pointer :: kg_k(:,:)
160 562 : real(dp) :: tmp_xc(2,Wfd%nspinor**2),tmp_xcval(2,Wfd%nspinor**2)
161 562 : real(dp) :: tmp_H(2,Wfd%nspinor**2),tmp_U(2,Wfd%nspinor**2)
162 562 : real(dp) :: tmp_h0ij(2,Wfd%nspinor**2),tmp_sigcx(2,Wfd%nspinor**2)
163 : real(dp) :: dijU(2),kpt(3),vxc1ab(2),vxc1ab_val(2)
164 281 : real(dp),allocatable :: kxc_(:,:),xccc3d_(:),vxc_val(:,:),vxc_val_hybrid(:,:), kinpw(:),veffh0(:,:)
165 : complex(dp) :: tmp(3)
166 281 : complex(gwp),contiguous, pointer :: ur1_up(:),ur1_dwn(:),ur2_up(:),ur2_dwn(:),cg1(:),cg2(:)
167 281 : complex(gwp),target,allocatable :: ur1(:),ur2(:)
168 281 : complex(dp),allocatable :: vxcab(:),vxcab_val(:),vxcab_val_hybrid(:),u1cjg_u2dpc(:),kinwf2(:),veffh0_ab(:)
169 281 : logical,allocatable :: bbp_mask(:,:)
170 281 : type(pawcprj_type),allocatable :: Cprj_b1ks(:,:),Cprj_b2ks(:,:)
171 843 : type(libxc_functional_type) :: xc_funcs_hybrid(2)
172 : ! *************************************************************************
173 :
174 : DBG_ENTER("COLL")
175 :
176 1124 : ABI_MALLOC(bbp_mask,(Wfd%mband, Wfd%mband))
177 :
178 : ! Usually FFT meshes for wavefunctions and potentials are not equal. Two approaches are possible:
179 : ! Either we Fourier interpolate potentials on the coarse WF mesh or we FFT the wfs on the dense mesh.
180 : ! The later approach is used, more CPU demanding but more accurate.
181 1115 : if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst, Psps, ngfftf)
182 :
183 : ! Fake MPI_type for sequential part
184 281 : rank = Wfd%my_rank
185 281 : call initmpi_seq(MPI_enreg_seq)
186 281 : call MPI_enreg_seq%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
187 :
188 281 : nspinor=Wfd%nspinor; nsppol =Wfd%nsppol; nspden =Wfd%nspden
189 281 : if (nspinor == 2) ABI_WARNING("Remember to ADD SO")
190 :
191 : ! TODO not used for the time being but it should be a standard input of the routine.
192 : ! bbks_mask(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)=Logical mask used to select
193 : ! the matrix elements to be calculated.
194 843 : ABI_MALLOC(kcalc2ibz,(Wfd%nkibz))
195 1978 : kcalc2ibz=0
196 :
197 : ! Index in the IBZ of the GW k-points.
198 : ! Only these points will be considered.
199 281 : nk_calc=0
200 1978 : do ik_ibz=1,Wfd%nkibz
201 4118 : if ( ALL(kstab(1,ik_ibz,:)/=0) .and. ALL(kstab(2,ik_ibz,:)/=0) ) then
202 1060 : nk_calc=nk_calc+1; kcalc2ibz(nk_calc) = ik_ibz
203 : end if
204 : end do
205 :
206 281 : call Mels%init(Mflags, nsppol, nspden, Wfd%nspinor, Wfd%nkibz, Wfd%kibz, kstab)
207 :
208 281 : if (Mflags%has_lexexch==1) then
209 0 : ABI_ERROR("Local EXX not coded!")
210 : end if
211 :
212 : ! Evaluate $v_\xc$ using only the valence charge.
213 281 : call wrtout(std_out," calc_vhxc_braket: calculating v_xc[n_val] (excluding non-linear core corrections)")
214 :
215 568 : do isppol=1,nsppol
216 5171730 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density rhor = ',MINVAL(rhor(:,isppol))
217 287 : call wrtout(std_out, msg)
218 568 : if (Wfd%usepaw==1) then
219 169482 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density nhat = ',MINVAL(nhat(:,isppol))
220 5 : call wrtout(std_out, msg)
221 169482 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density trho-nhat = ',MINVAL(rhor(:,isppol)-nhat(:,isppol))
222 5 : call wrtout(std_out, msg)
223 5 : write(msg,'(a,i2)')' using usexcnhat = ',usexcnhat
224 5 : call wrtout(std_out, msg)
225 : end if
226 : end do
227 :
228 281 : option = 0 ! Only exc, vxc
229 281 : nkxc = 0 ! No computation of XC kernel
230 281 : n3xccc_= 0 ! No core
231 281 : nk3xc = 0 ! k3xc not needed
232 281 : nmxc = Dtset%usepaw==1 .and. mod(abs(Dtset%usepawu),10) == 4
233 :
234 281 : ABI_MALLOC(xccc3d_,(n3xccc_))
235 562 : ABI_MALLOC(kxc_,(nfftf,nkxc))
236 1124 : ABI_MALLOC(vxc_val,(nfftf,nspden))
237 :
238 281 : call xcdata_init(xcdata,dtset=Dtset)
239 :
240 : call rhotoxc(bigexc_val,bigsxc_val,kxc_,MPI_enreg_seq,nfftf,ngfftf,&
241 : nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,&
242 281 : usexcnhat,vxc_val,vxcval_avg,xccc3d_,xcdata,taur=taur)
243 :
244 : ! FABIEN's development
245 : ! Hybrid functional treatment
246 281 : if (Mflags%has_vxcval_hybrid ==1 ) then
247 :
248 23 : call wrtout(std_out,' Hybrid functional xc potential is being set')
249 23 : ixc_sigma=Dtset%ixc_sigma
250 23 : call get_auxc_ixc(auxc_ixc,ixc_sigma)
251 23 : call xcdata_init(xcdata_hybrid,dtset=Dtset,auxc_ixc=auxc_ixc,ixc=ixc_sigma)
252 :
253 23 : if(ixc_sigma<0)then
254 18 : if(libxc_functionals_check()) then
255 18 : call libxc_functionals_init(ixc_sigma,Dtset%nspden,xc_functionals=xc_funcs_hybrid)
256 : ! Do not forget, negative values of hyb_mixing(_sr),hyb_range_* means that they have been user-defined.
257 18 : if (dtset%ixc==-402.or.dtset%ixc==-406.or.dtset%ixc==-427.or.dtset%ixc==-428 .or. dtset%ixc==-456 .or. &
258 : & min(Dtset%hyb_mixing,Dtset%hyb_mixing_sr,Dtset%hyb_range_dft,Dtset%hyb_range_fock)<-tol8)then
259 : call libxc_functionals_set_hybridparams(hyb_range=abs(Dtset%hyb_range_dft),&
260 1 : & hyb_mixing=abs(Dtset%hyb_mixing),hyb_mixing_sr=abs(Dtset%hyb_mixing_sr),xc_functionals=xc_funcs_hybrid)
261 : endif
262 : else
263 0 : call wrtout(std_out, 'LIBXC is not present: hybrid functionals are not available')
264 : end if
265 : end if
266 :
267 23 : write(msg, '(a, f4.2)') ' Fock fraction = ', max(abs(Dtset%hyb_mixing),abs(Dtset%hyb_mixing_sr))
268 23 : call wrtout(std_out, msg)
269 23 : write(msg, '(a, f5.2, a)') ' Fock inverse screening length = ',abs(Dtset%hyb_range_dft), ' (bohr^-1)'
270 23 : call wrtout(std_out, msg)
271 :
272 69 : ABI_MALLOC(vxc_val_hybrid,(nfftf,nspden))
273 :
274 46 : if(ixc_sigma<0)then
275 : call rhotoxc(bigexc_hybrid_val,dum,kxc_,MPI_enreg_seq,nfftf,ngfftf,&
276 : nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,&
277 18 : usexcnhat,vxc_val_hybrid,vxcval_hybrid_avg,xccc3d_,xcdata_hybrid,xc_funcs=xc_funcs_hybrid)
278 18 : call libxc_functionals_end(xc_functionals=xc_funcs_hybrid)
279 : else
280 : call rhotoxc(bigexc_hybrid_val,dum,kxc_,MPI_enreg_seq,nfftf,ngfftf,&
281 : nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,&
282 5 : usexcnhat,vxc_val_hybrid,vxcval_hybrid_avg,xccc3d_,xcdata_hybrid)
283 : end if
284 :
285 : endif
286 :
287 281 : ABI_FREE(xccc3d_)
288 281 : ABI_FREE(kxc_)
289 :
290 281 : write(msg,'(a,f8.4,2a,f8.4,a)')' E_xc[n_val] = ',bigexc_val, ' [Ha]. ','<V_xc[n_val]> = ',vxcval_avg,' [Ha]. '
291 281 : call wrtout(std_out, msg)
292 281 : if(bigsxc_val>zero) then
293 0 : write(msg,'(a,f8.4)')' S_xc[n_val] = ',bigsxc_val
294 0 : call wrtout(std_out, msg)
295 : end if
296 :
297 : ! has_hbare uses veffh0. Why only use it with usepaw=1? Let it be available always
298 281 : if (Mflags%has_hbare == 1) then
299 9 : if (Mflags%has_kinetic/=1) then
300 0 : ABI_ERROR("Kinetic energy mels are required for the construction of Hbare mels!")
301 : end if
302 : ! Effective potential of the bare Hamiltonian: valence term is subtracted.
303 27 : ABI_MALLOC(veffh0,(nfftf,nspden))
304 30027 : veffh0=vtrial-vxc_val
305 : !veffh0=vtrial !this is to retrieve the KS Hamiltonian
306 : endif
307 :
308 : ! If PAW and qp-SCGW then update Paw_ij and calculate the matrix elements ===
309 : ! We cannot simply rely on gwcalctyp because I need KS vxc in sigma.
310 281 : if (Wfd%usepaw==1.and.Mflags%has_hbare==1) then
311 0 : ABI_CHECK(Mflags%only_diago==0,"Wrong only_diago")
312 :
313 : call paw_mknewh0(Cryst%natom,nsppol,nspden,nfftf,Dtset%pawspnorb,Dtset%pawprtvol,Cryst,&
314 0 : Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial)
315 :
316 : ! Effective potential of the bare Hamiltonian: valence term is subtracted.
317 0 : veffh0=vtrial-vxc_val
318 : !veffh0=vtrial !this is to retrieve the KS Hamiltonian
319 : end if
320 :
321 : ! Setup of the hermitian operator vxcab ===
322 : ! if nspden==4 vxc contains (v^11, v^22, Re[V^12], Im[V^12].
323 : ! Cannot use directly Re and Im since we also need off-diagonal elements.
324 281 : if (wfd%nspden == 4) then
325 15 : ABI_MALLOC(vxcab, (nfftf))
326 10 : ABI_MALLOC(vxcab_val, (nfftf))
327 40005 : vxcab (:) = DCMPLX(vxc (:,3), vxc (:,4))
328 40005 : vxcab_val(:) = DCMPLX(vxc_val(:,3), vxc_val(:,4))
329 5 : if (Mflags%has_vxcval_hybrid==1) then
330 0 : ABI_MALLOC(vxcab_val_hybrid,(nfftf))
331 0 : vxcab_val_hybrid(:)=DCMPLX(vxc_val_hybrid(:,3),vxc_val_hybrid(:,4))
332 : end if
333 5 : if (Mflags%has_hbare==1) then
334 0 : ABI_MALLOC(veffh0_ab,(nfftf))
335 0 : veffh0_ab(:)=DCMPLX(veffh0(:,3),veffh0(:,4))
336 : end if
337 : end if
338 :
339 843 : ABI_MALLOC(ur1, (nfftf * nspinor))
340 562 : ABI_MALLOC(ur2, (nfftf * nspinor))
341 562 : ABI_MALLOC(u1cjg_u2dpc, (nfftf * nspinor))
342 :
343 : ! Create distribution table for tasks.
344 : ! This section is parallelized inside wfd%comm
345 : ! as all processors are calling the routine with all GW wavefunctions
346 : ! TODO the table can be calculated at each (k,s) to save some memory.
347 716 : got=0; my_nmels=0
348 1686 : ABI_MALLOC(bbp_ks_distrb,(Wfd%mband,Wfd%mband,nk_calc,nsppol))
349 568 : do is=1,nsppol
350 1638 : do ikc=1,nk_calc
351 1070 : ik_ibz=kcalc2ibz(ikc)
352 583968 : bbp_mask=.FALSE.
353 1070 : b_start=kstab(1,ik_ibz,is)
354 1070 : b_stop =kstab(2,ik_ibz,is)
355 1070 : if (Mflags%only_diago==1) then
356 : !do jb=b1,b2
357 1705 : do jb=b_start,b_stop
358 1705 : bbp_mask(jb,jb)=.TRUE.
359 : end do
360 : else
361 97974 : bbp_mask(b_start:b_stop,b_start:b_stop)=.TRUE.
362 : end if
363 :
364 1070 : call wfd%distribute_bbp(ik_ibz,is,"Upper",my_nbbp,bbp_ks_distrb(:,:,ikc,is),got,bbp_mask)
365 1357 : my_nmels = my_nmels + my_nbbp
366 : end do
367 : end do
368 281 : ABI_FREE(bbp_mask)
369 :
370 281 : write(msg,'(a,i0,a)')" Will calculate ",my_nmels," <b,k,s|O|b',k,s> matrix elements in calc_vhxc_me."
371 281 : call wrtout(std_out, msg)
372 :
373 : ! =====================================
374 : ! ==== Loop over required k-points ====
375 : ! =====================================
376 281 : nfftfm1=one/nfftf
377 :
378 568 : do is=1,nsppol
379 15127 : if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE
380 1638 : do ikc=1,nk_calc
381 20689 : if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE
382 :
383 1070 : ik_ibz = kcalc2ibz(ikc)
384 1070 : b_start = kstab(1,ik_ibz,is)
385 1070 : b_stop = kstab(2,ik_ibz,is)
386 1070 : npw_k = Wfd%Kdata(ik_ibz)%npw
387 4280 : kpt = Wfd%kibz(:,ik_ibz)
388 1070 : kg_k => Wfd%kdata(ik_ibz)%kg_k
389 1070 : istwf_k = wfd%istwfk(ik_ibz)
390 :
391 : ! Calculate |k+G|^2 needed by hbareme and kineticme
392 : ! MRM: Solved ecut problem
393 1070 : if (Mflags%has_kinetic == 1) then
394 162 : ABI_MALLOC(kinpw, (npw_k))
395 162 : ABI_MALLOC(kinwf2, (npw_k*nspinor))
396 54 : call mkkin(Dtset%ecutwfn,Dtset%ecutsm,Dtset%effmass_free,Cryst%gmet,kg_k,kinpw,kpt,npw_k,0,0)
397 6423 : where (kinpw>HUGE(zero)*1.d-11)
398 : kinpw=zero
399 : end where
400 : end if
401 :
402 : !do jb=b1,b2
403 10707 : do jb=b_start,b_stop
404 23954 : if (ALL(bbp_ks_distrb(:,jb,ikc,is)/=rank)) CYCLE
405 :
406 9081 : ABI_CHECK(wfd%get_wave_ptr(jb, ik_ibz, is, wave_jb, msg) == 0, msg)
407 :
408 9081 : if (Mflags%has_kinetic == 1) then
409 412 : cg2 => wave_jb%ug
410 45522 : kinwf2(1:npw_k) = cg2(1:npw_k)*kinpw(:)
411 412 : if (nspinor==2) kinwf2(npw_k+1:)=cg2(npw_k+1:)*kinpw(:)
412 : end if
413 :
414 9081 : call wfd%get_ur(jb,ik_ibz,is,ur2)
415 :
416 : !do ib=b1,jb ! Upper triangle
417 63369 : do ib=b_start,jb
418 53218 : if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE
419 : ! Off-diagonal elements only for QPSCGW.
420 44241 : if (Mflags%only_diago==1.and.ib/=jb) CYCLE
421 :
422 44241 : call wfd%get_ur(ib,ik_ibz,is,ur1)
423 261242655 : u1cjg_u2dpc(:) = CONJG(ur1) *ur2
424 :
425 44241 : if (Mflags%has_vxc == 1) then
426 137346426 : Mels%vxc(ib, jb, ik_ibz, is) = sum(u1cjg_u2dpc(1:nfftf) * vxc(1:nfftf, is)) * nfftfm1
427 23129 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) &
428 745024 : Mels%vxc(ib, jb, ik_ibz, 2) = sum(u1cjg_u2dpc(nfftf+1:) * vxc(1:nfftf, is)) * nfftfm1
429 : end if
430 44241 : if (Mflags%has_vxcval == 1) then
431 137346426 : Mels%vxcval(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vxc_val(1:nfftf, is)) * nfftfm1
432 23129 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) &
433 745024 : Mels%vxcval(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1:) * vxc_val(1:nfftf, is)) * nfftfm1
434 : end if
435 44241 : if (Mflags%has_vxcval_hybrid == 1) then
436 37260828 : Mels%vxcval_hybrid(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vxc_val_hybrid(1:nfftf, is)) * nfftfm1
437 12024 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) &
438 0 : Mels%vxcval_hybrid(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1) * vxc_val_hybrid(1:nfftf, is)) * nfftfm1
439 : end if
440 44241 : if (Mflags%has_vhartree==1) then
441 259857695 : Mels%vhartree(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vhartr(1:nfftf)) * nfftfm1
442 44241 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) &
443 745024 : Mels%vhartree(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1:) * vhartr(1:nfftf)) * nfftfm1
444 : end if
445 44241 : if (Mflags%has_kinetic==1) then
446 1813 : ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, is, wave_ib, msg) == 0, msg)
447 1813 : cg1 => wave_ib%ug(1:npw_k)
448 192893 : cdot = DOT_PRODUCT(cg1, kinwf2(1:npw_k))
449 : !if (istwf_k /= 1) then
450 : ! cdot = two * cdot; if (istwf_k == 2) cdot = cdot - GWPC_CONJG(cg1(1)) * kinwf2(1)
451 : !end if
452 1813 : Mels%kinetic(ib, jb, ik_ibz, is) = cdot
453 1813 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) then
454 0 : cg1 => wave_ib%ug(npw_k+1:)
455 0 : Mels%kinetic(ib, jb, ik_ibz, 2) = DOT_PRODUCT(cg1, kinwf2(npw_k+1:))
456 : end if
457 : end if
458 44241 : if (Mflags%has_hbare==1) then
459 : Mels%hbare(ib, jb, ik_ibz, is) = Mels%kinetic(ib, jb, ik_ibz, is) &
460 5433813 : + SUM(u1cjg_u2dpc(1:nfftf) * veffh0(1:nfftf, is)) * nfftfm1
461 1813 : if (wfd%nspinor == 2 .and. wfd%nspden == 1) then
462 0 : cg1 => wave_ib%ug(npw_k+1:)
463 : Mels%hbare(ib, jb, ik_ibz, 2) = &
464 0 : Mels%kinetic(ib, jb, ik_ibz, 2) + SUM(u1cjg_u2dpc(nfftf+1:) * veffh0(1:nfftf, is)) * nfftfm1
465 : end if
466 : end if
467 :
468 53878 : if (nspinor == 2 .and. wfd%nspden == 4) then
469 : ! Here I can skip 21 if ib==jb
470 80 : ur1_up => ur1(1:nfftf)
471 80 : ur1_dwn => ur1(nfftf+1:2*nfftf)
472 80 : ur2_up => ur2(1:nfftf)
473 80 : ur2_dwn => ur2(nfftf+1:2*nfftf)
474 :
475 80 : if (Mflags%has_kinetic==1) then
476 0 : cg1 => wave_ib%ug(npw_k+1:)
477 0 : tmp(1)=DOT_PRODUCT(cg1,kinwf2(npw_k+1:))
478 0 : Mels%kinetic(ib,jb,ik_ibz,2 )=tmp(1)
479 0 : Mels%kinetic(ib,jb,ik_ibz,3:4)=czero
480 : end if
481 80 : if (Mflags%has_hbare==1) then
482 0 : cg1 => wave_ib%ug(npw_k+1:)
483 0 : tmp(1)=SUM(CONJG(ur1_dwn)*veffh0(:,2)*ur2_dwn)*nfftfm1 + Mels%kinetic(ib,jb,ik_ibz,2)
484 0 : tmp(2)=SUM(CONJG(ur1_dwn)* veffh0_ab(:) *ur2_dwn)*nfftfm1
485 0 : tmp(3)=SUM(CONJG(ur1_dwn)*CONJG(veffh0_ab(:))*ur2_dwn)*nfftfm1
486 0 : Mels%hbare(ib,jb,ik_ibz,2:4)=tmp(:)
487 : end if
488 80 : if (Mflags%has_vxc==1) then
489 640080 : tmp(1) = SUM(CONJG(ur1_dwn)* vxc(:,2) *ur2_dwn)*nfftfm1
490 640080 : tmp(2) = SUM(CONJG(ur1_up )* vxcab(:) *ur2_dwn)*nfftfm1
491 640080 : tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab(:))*ur2_up )*nfftfm1
492 320 : Mels%vxc(ib,jb,ik_ibz,2:4)=tmp(:)
493 : end if
494 80 : if (Mflags%has_vxcval==1) then
495 640080 : tmp(1) = SUM(CONJG(ur1_dwn)* vxc_val(:,2) *ur2_dwn)*nfftfm1
496 640080 : tmp(2) = SUM(CONJG(ur1_up )* vxcab_val(:) *ur2_dwn)*nfftfm1
497 640080 : tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab_val(:))*ur2_up )*nfftfm1
498 320 : Mels%vxcval(ib,jb,ik_ibz,2:4)=tmp(:)
499 : end if
500 80 : if (Mflags%has_vxcval_hybrid==1) then
501 0 : tmp(1) = SUM(CONJG(ur1_dwn)* vxc_val_hybrid(:,2) *ur2_dwn)*nfftfm1
502 0 : tmp(2) = SUM(CONJG(ur1_up )* vxcab_val_hybrid(:) *ur2_dwn)*nfftfm1
503 0 : tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab_val_hybrid(:))*ur2_up )*nfftfm1
504 0 : Mels%vxcval_hybrid(ib,jb,ik_ibz,2:4)=tmp(:)
505 : end if
506 80 : if (Mflags%has_vhartree==1) then
507 640080 : tmp(1) = SUM(CONJG(ur1_dwn)*vhartr(:)*ur2_dwn)*nfftfm1
508 80 : Mels%vhartree(ib,jb,ik_ibz,2 )=tmp(1)
509 240 : Mels%vhartree(ib,jb,ik_ibz,3:4)=czero
510 : end if
511 : end if !nspinor==2
512 :
513 : end do !ib
514 : end do !jb
515 :
516 1357 : if (Mflags%has_kinetic==1) then
517 54 : ABI_FREE(kinpw)
518 54 : ABI_FREE(kinwf2)
519 : end if
520 :
521 : end do !ikc
522 : end do !is
523 :
524 281 : ABI_FREE(ur1)
525 281 : ABI_FREE(ur2)
526 281 : ABI_FREE(vxc_val)
527 281 : ABI_FREE(u1cjg_u2dpc)
528 281 : if(Mflags%has_vxcval_hybrid==1) then
529 23 : ABI_FREE(vxc_val_hybrid)
530 : end if
531 281 : if (wfd%nspden == 4) then
532 5 : ABI_FREE(vxcab)
533 5 : ABI_FREE(vxcab_val)
534 5 : if(Mflags%has_vxcval_hybrid==1) then
535 0 : ABI_FREE(vxcab_val_hybrid)
536 : end if
537 : end if
538 :
539 281 : if (Mflags%has_hbare==1) then
540 9 : ABI_FREE(veffh0)
541 9 : if (nspinor==2) then
542 0 : ABI_FREE(veffh0_ab)
543 : end if
544 : end if
545 :
546 : ! ====================================
547 : ! ===== Additional terms for PAW =====
548 : ! ====================================
549 281 : if (Wfd%usepaw==1) then
550 : ! Tests if needed pointers in Paw_ij are allocated.
551 5 : ltest=(allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_hat).and.allocated(Paw_ij(1)%dijxc_val))
552 0 : ABI_CHECK(ltest,"dijxc, dijxc_hat or dijxc_val not allocated")
553 5 : ABI_CHECK(nspinor == 1, "PAW with nspinor not tested")
554 :
555 : ! For DFT+U
556 36 : do iat=1,Cryst%natom
557 31 : itypat=Cryst%typat(iat)
558 36 : if (Pawtab(itypat)%usepawu/=0) then
559 0 : ltest=(allocated(Paw_ij(iat)%dijU))
560 0 : ABI_CHECK(ltest,"DFT+U but dijU not allocated")
561 : end if
562 : end do
563 :
564 5 : if (Dtset%pawspnorb>0) then
565 0 : ltest=(allocated(Paw_ij(1)%dijso))
566 0 : ABI_CHECK(ltest,"dijso not allocated")
567 : end if
568 :
569 13 : lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
570 :
571 5 : if (Mflags%has_sxcore==1) then
572 : if ( SIZE(dijexc_core,DIM=1) /= lmn2_size_max &
573 : .or. SIZE(dijexc_core,DIM=2) /= 1 &
574 1 : .or. SIZE(dijexc_core,DIM=3) /= Cryst%ntypat ) then
575 0 : ABI_BUG("Wrong sizes in dijexc_core")
576 : end if
577 : end if
578 :
579 5 : nsploop=nspinor**2
580 :
581 : ! ====================================
582 : ! === Assemble PAW matrix elements ===
583 : ! ====================================
584 15 : ABI_MALLOC(dimlmn, (Cryst%natom))
585 36 : do iat=1,Cryst%natom
586 36 : dimlmn(iat)=Pawtab(Cryst%typat(iat))%lmn_size
587 : end do
588 :
589 56 : ABI_MALLOC(Cprj_b1ks, (Cryst%natom,nspinor))
590 51 : ABI_MALLOC(Cprj_b2ks, (Cryst%natom,nspinor))
591 5 : call pawcprj_alloc(Cprj_b1ks, 0, dimlmn)
592 5 : call pawcprj_alloc(Cprj_b2ks, 0, dimlmn)
593 :
594 10 : do is=1,nsppol
595 2075 : if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE
596 :
597 : ! Loop over required k-points
598 19 : do ikc=1,nk_calc
599 2079 : if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE
600 9 : ik_ibz=kcalc2ibz(ikc)
601 9 : b_start=kstab(1,ik_ibz,is)
602 9 : b_stop =kstab(2,ik_ibz,is)
603 :
604 : !do jb=b1,b2
605 60 : do jb=b_start,b_stop
606 301 : if (ALL(bbp_ks_distrb(:,jb,ikc,is)/=rank)) CYCLE
607 :
608 : ! Load projected wavefunctions for this k-point, spin and band ===
609 : ! Cprj are unsorted, full correspondence with xred. See ctocprj.F90!!
610 46 : call wfd%get_cprj(jb,ik_ibz,is,Cryst,Cprj_b2ks,sorted=.FALSE.)
611 :
612 : !do ib=b1,jb ! Upper triangle
613 218 : do ib=b_start,jb
614 163 : if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE
615 :
616 : ! Off-diagonal elements only for QPSCGW.
617 46 : if (Mflags%only_diago==1.and.ib/=jb) CYCLE
618 :
619 46 : call wfd%get_cprj(ib,ik_ibz,is,Cryst,Cprj_b1ks,sorted=.FALSE.)
620 : !
621 : ! === Get onsite matrix elements summing over atoms and channels ===
622 : ! * Spin is external and fixed (1,2) if collinear.
623 : ! * if noncollinear loop internally over the four components ab.
624 1104 : tmp_xc = zero; tmp_xcval = zero; tmp_H = zero; tmp_U = zero; tmp_h0ij = zero; tmp_sigcx = zero
625 :
626 180 : do iat=1,Cryst%natom
627 134 : itypat =Cryst%typat(iat)
628 134 : lmn_size =Pawtab(itypat)%lmn_size
629 134 : cplex_dij=Paw_ij(iat)%cplex_dij
630 134 : klmn1=1
631 :
632 1342 : do jlmn=1,lmn_size
633 1162 : j0lmn=jlmn*(jlmn-1)/2
634 7110 : do ilmn=1,jlmn
635 5814 : klmn=j0lmn+ilmn
636 : ! TODO Be careful, here I assume that the onsite terms ij are symmetric
637 : ! should check the spin-orbit case!
638 5814 : fact=one; if (ilmn==jlmn) fact=half
639 :
640 : ! Loop over four components if nspinor==2
641 : ! If collinear nsploop==1
642 11628 : do iab=1,nsploop
643 5814 : isp1=spinor_idxs(1,iab); isp2=spinor_idxs(2,iab)
644 :
645 : re_p= Cprj_b1ks(iat,isp1)%cp(1,ilmn) * Cprj_b2ks(iat,isp2)%cp(1,jlmn) &
646 : +Cprj_b1ks(iat,isp1)%cp(2,ilmn) * Cprj_b2ks(iat,isp2)%cp(2,jlmn) &
647 : +Cprj_b1ks(iat,isp1)%cp(1,jlmn) * Cprj_b2ks(iat,isp2)%cp(1,ilmn) &
648 5814 : +Cprj_b1ks(iat,isp1)%cp(2,jlmn) * Cprj_b2ks(iat,isp2)%cp(2,ilmn)
649 :
650 : im_p= Cprj_b1ks(iat,isp1)%cp(1,ilmn) * Cprj_b2ks(iat,isp2)%cp(2,jlmn) &
651 : -Cprj_b1ks(iat,isp1)%cp(2,ilmn) * Cprj_b2ks(iat,isp2)%cp(1,jlmn) &
652 : +Cprj_b1ks(iat,isp1)%cp(1,jlmn) * Cprj_b2ks(iat,isp2)%cp(2,ilmn) &
653 5814 : -Cprj_b1ks(iat,isp1)%cp(2,jlmn) * Cprj_b2ks(iat,isp2)%cp(1,ilmn)
654 :
655 : ! ==================================================
656 : ! === Load onsite matrix elements and accumulate ===
657 : ! ==================================================
658 11628 : if (nspinor==1) then
659 :
660 5814 : if (Mflags%has_hbare==1) then ! * Get new dij of h0 and accumulate.
661 0 : h0dij=Paw_ij(iat)%dij(klmn,is)
662 0 : tmp_h0ij(1,iab)=tmp_h0ij(1,iab) + h0dij*re_p*fact
663 0 : tmp_h0ij(2,iab)=tmp_h0ij(2,iab) + h0dij*im_p*fact
664 : end if
665 :
666 5814 : if (Mflags%has_sxcore==1) then ! * Fock operator generated by core electrons.
667 1440 : dijsigcx = dijexc_core(klmn,1,itypat)
668 1440 : tmp_sigcx(1,iab)=tmp_sigcx(1,iab) + dijsigcx*re_p*fact
669 1440 : tmp_sigcx(2,iab)=tmp_sigcx(2,iab) + dijsigcx*im_p*fact
670 : end if
671 :
672 5814 : if (Mflags%has_vxc==1) then ! * Accumulate vxc[n1+nc] + vxc[n1+tn+nc].
673 5814 : vxc1 = Paw_ij(iat)%dijxc(klmn,is)+Paw_ij(iat)%dijxc_hat(klmn,is)
674 5814 : tmp_xc(1,iab)=tmp_xc(1,iab) + vxc1*re_p*fact
675 5814 : tmp_xc(2,iab)=tmp_xc(2,iab) + vxc1*im_p*fact
676 : end if
677 :
678 5814 : if (Mflags%has_vxcval==1) then ! * Accumulate valence-only XC.
679 5814 : vxc1_val=Paw_ij(iat)%dijxc_val(klmn,is)
680 5814 : tmp_xcval(1,1)=tmp_xcval(1,1) + vxc1_val*re_p*fact
681 5814 : tmp_xcval(2,1)=tmp_xcval(2,1) + vxc1_val*im_p*fact
682 : end if
683 :
684 5814 : if (Mflags%has_vhartree==1) then ! * Accumulate Hartree term of the PAW Hamiltonian.
685 5814 : DijH=Paw_ij(iat)%dijhartree(klmn)
686 5814 : tmp_H(1,1)=tmp_H(1,1) + DijH*re_p*fact
687 5814 : tmp_H(2,1)=tmp_H(2,1) + DijH*im_p*fact
688 : end if
689 :
690 : ! Accumulate U term of the PAW Hamiltonian (only onsite AE contribution)
691 5814 : if (Mflags%has_vu==1) then
692 0 : if (Pawtab(itypat)%usepawu/=0) then
693 0 : dijU(1)=Paw_ij(iat)%dijU(klmn,is)
694 0 : tmp_U(1,1)=tmp_U(1,1) + dijU(1)*re_p*fact
695 0 : tmp_U(2,1)=tmp_U(2,1) + dijU(1)*im_p*fact
696 : end if
697 : end if
698 :
699 : else
700 : ! Spinorial case
701 :
702 : ! FIXME H0 + spinor not implemented
703 0 : if (Mflags%has_hbare==1.or.Mflags%has_sxcore==1) then
704 0 : ABI_ERROR("not implemented")
705 : end if
706 :
707 0 : if (Mflags%has_vxc==1) then ! * Accumulate vxc[n1+nc] + vxc[n1+tn+nc].
708 0 : vxc1ab(1) = Paw_ij(iat)%dijxc(klmn1, iab)+Paw_ij(iat)%dijxc_hat(klmn1, iab)
709 0 : vxc1ab(2) = Paw_ij(iat)%dijxc(klmn1+1,iab)+Paw_ij(iat)%dijxc_hat(klmn1+1,iab)
710 0 : tmp_xc(1,iab) = tmp_xc(1,iab) + (vxc1ab(1)*re_p - vxc1ab(2)*im_p)*fact
711 0 : tmp_xc(2,iab) = tmp_xc(2,iab) + (vxc1ab(2)*re_p + vxc1ab(1)*im_p)*fact
712 : end if
713 :
714 0 : if (Mflags%has_vxcval==1) then ! * Accumulate valence-only XC.
715 0 : vxc1ab_val(1) = Paw_ij(iat)%dijxc_val(klmn1, iab)
716 0 : vxc1ab_val(2) = Paw_ij(iat)%dijxc_val(klmn1+1,iab)
717 0 : tmp_xcval(1,iab) = tmp_xcval(1,iab) + (vxc1ab_val(1)*re_p - vxc1ab_val(2)*im_p)*fact
718 0 : tmp_xcval(2,iab) = tmp_xcval(2,iab) + (vxc1ab_val(2)*re_p + vxc1ab_val(1)*im_p)*fact
719 : end if
720 :
721 : ! * In GW, dijhartree is always real.
722 0 : if (Mflags%has_vhartree==1) then ! * Accumulate Hartree term of the PAW Hamiltonian.
723 0 : if (iab==1.or.iab==2) then
724 0 : DijH = Paw_ij(iat)%dijhartree(klmn)
725 0 : tmp_H(1,iab) = tmp_H(1,iab) + DijH*re_p*fact
726 0 : tmp_H(2,iab) = tmp_H(2,iab) + DijH*im_p*fact
727 : end if
728 : end if
729 :
730 : ! TODO "ADD DFT+U and SO"
731 : ! check this part
732 0 : if (Mflags%has_vu==1) then
733 0 : if (Pawtab(itypat)%usepawu/=0) then
734 : ! Accumulate the U term of the PAW Hamiltonian (only onsite AE contribution)
735 0 : dijU(1)=Paw_ij(iat)%dijU(klmn1 ,iab)
736 0 : dijU(2)=Paw_ij(iat)%dijU(klmn1+1,iab)
737 0 : tmp_U(1,iab) = tmp_U(1,iab) + (dijU(1)*re_p - dijU(2)*im_p)*fact
738 0 : tmp_U(2,iab) = tmp_U(2,iab) + (dijU(2)*re_p + dijU(1)*im_p)*fact
739 : end if
740 : end if
741 :
742 : end if
743 : end do !iab
744 :
745 6976 : klmn1=klmn1+cplex_dij
746 :
747 : end do !ilmn
748 : end do !jlmn
749 : end do !iat
750 :
751 : ! ========================================
752 : ! ==== Add to plane wave contribution ====
753 : ! ========================================
754 92 : if (nspinor==1) then
755 :
756 46 : if (Mflags%has_hbare==1) &
757 0 : & Mels%hbare(ib,jb,ik_ibz,is) = Mels%hbare(ib,jb,ik_ibz,is) + DCMPLX(tmp_h0ij(1,1),tmp_h0ij(2,1))
758 :
759 46 : if (Mflags%has_vxc==1) &
760 46 : & Mels%vxc(ib,jb,ik_ibz,is) = Mels%vxc(ib,jb,ik_ibz,is) + DCMPLX(tmp_xc(1,1),tmp_xc(2,1))
761 :
762 46 : if (Mflags%has_vxcval==1) &
763 46 : & Mels%vxcval(ib,jb,ik_ibz,is) = Mels%vxcval(ib,jb,ik_ibz,is) + DCMPLX(tmp_xcval(1,1),tmp_xcval(2,1))
764 :
765 46 : if (Mflags%has_vxcval_hybrid==1) &
766 0 : & Mels%vxcval_hybrid(ib,jb,ik_ibz,is) = Mels%vxcval_hybrid(ib,jb,ik_ibz,is) + DCMPLX(tmp_xcval(1,1),tmp_xcval(2,1))
767 :
768 46 : if (Mflags%has_vhartree==1) &
769 46 : & Mels%vhartree(ib,jb,ik_ibz,is) = Mels%vhartree(ib,jb,ik_ibz,is) + DCMPLX(tmp_H (1,1),tmp_H (2,1))
770 :
771 46 : if (Mflags%has_vu==1) &
772 0 : & Mels%vu(ib,jb,ik_ibz,is) = DCMPLX(tmp_U(1,1),tmp_U(2,1))
773 :
774 46 : if (Mflags%has_sxcore==1) &
775 20 : & Mels%sxcore(ib,jb,ik_ibz,is) = DCMPLX(tmp_sigcx(1,1),tmp_sigcx(2,1))
776 :
777 : else
778 :
779 0 : if (Mflags%has_hbare==1) &
780 0 : & Mels%hbare(ib,jb,ik_ibz,:) = Mels%hbare(ib,jb,ik_ibz,:) + DCMPLX(tmp_h0ij(1,:),tmp_h0ij(2,:))
781 :
782 0 : if (Mflags%has_vxc==1) &
783 0 : & Mels%vxc(ib,jb,ik_ibz,:) = Mels%vxc(ib,jb,ik_ibz,:) + DCMPLX(tmp_xc(1,:),tmp_xc(2,:))
784 :
785 0 : if (Mflags%has_vxcval==1) &
786 0 : & Mels%vxcval(ib,jb,ik_ibz,:) = Mels%vxcval(ib,jb,ik_ibz,:) + DCMPLX(tmp_xcval(1,:),tmp_xcval(2,:))
787 :
788 0 : if (Mflags%has_vxcval_hybrid==1) &
789 0 : & Mels%vxcval_hybrid(ib,jb,ik_ibz,:) = Mels%vxcval_hybrid(ib,jb,ik_ibz,:) + DCMPLX(tmp_xcval(1,:),tmp_xcval(2,:))
790 :
791 0 : if (Mflags%has_vhartree==1) &
792 0 : & Mels%vhartree(ib,jb,ik_ibz,:) = Mels%vhartree(ib,jb,ik_ibz,:) + DCMPLX(tmp_H (1,:),tmp_H (2,:))
793 :
794 0 : if (Mflags%has_vu==1) &
795 0 : & Mels%vu(ib,jb,ik_ibz,:) = DCMPLX(tmp_U(1,:),tmp_U(2,:))
796 : end if
797 :
798 : end do !ib
799 : end do !jb
800 :
801 : end do !is
802 : end do !ikc
803 :
804 5 : ABI_FREE(dimlmn)
805 5 : call pawcprj_free(Cprj_b1ks)
806 36 : ABI_FREE(Cprj_b1ks)
807 5 : call pawcprj_free(Cprj_b2ks)
808 36 : ABI_FREE(Cprj_b2ks)
809 : end if !PAW
810 :
811 281 : ABI_FREE(bbp_ks_distrb)
812 :
813 : ! Sum up contributions on each node
814 : ! Set the corresponding has_* flags to 2.
815 281 : call Mels%mpisum(wfd%comm)
816 :
817 : ! Reconstruct lower triangle.
818 281 : call Mels%herm()
819 :
820 281 : ABI_FREE(kcalc2ibz)
821 281 : call destroy_mpi_enreg(MPI_enreg_seq)
822 :
823 : DBG_EXIT("COLL")
824 :
825 1124 : end subroutine calc_vhxc_me
826 : !!***
827 :
828 : end module m_vhxc_me
829 : !!***
|
