Branch data 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 [ - + - + : 286 : subroutine calc_vhxc_me(Wfd, Mflags, Mels, Cryst, Dtset, nfftf, ngfftf, &
- + - + -
+ - + - +
- + - + -
+ - + ]
114 [ - + ]: 286 : vtrial, vhartr, vxc, Psps, Pawtab, Paw_an, Pawang, Pawfgrtab, Paw_ij, dijexc_core, &
115 : 286 : rhor, usexcnhat, nhat, nhatgr, nhatgrdim, kstab, &
116 : 286 : 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 : 286 : 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 : 572 : integer :: got(Wfd%nproc)
158 : 286 : integer,allocatable :: kcalc2ibz(:),dimlmn(:),bbp_ks_distrb(:,:,:,:)
159 : 286 : integer,contiguous, pointer :: kg_k(:,:)
160 : 572 : real(dp) :: tmp_xc(2,Wfd%nspinor**2),tmp_xcval(2,Wfd%nspinor**2)
161 : 572 : real(dp) :: tmp_H(2,Wfd%nspinor**2),tmp_U(2,Wfd%nspinor**2)
162 : 572 : 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 [ - + ]: 286 : real(dp),allocatable :: kxc_(:,:),xccc3d_(:),vxc_val(:,:),vxc_val_hybrid(:,:), kinpw(:),veffh0(:,:)
165 : : complex(dp) :: tmp(3)
166 : 286 : complex(gwp),contiguous, pointer :: ur1_up(:),ur1_dwn(:),ur2_up(:),ur2_dwn(:),cg1(:),cg2(:)
167 : 286 : complex(gwp),target,allocatable :: ur1(:),ur2(:)
168 : 286 : complex(dp),allocatable :: vxcab(:),vxcab_val(:),vxcab_val_hybrid(:),u1cjg_u2dpc(:),kinwf2(:),veffh0_ab(:)
169 : 286 : logical,allocatable :: bbp_mask(:,:)
170 : 286 : type(pawcprj_type),allocatable :: Cprj_b1ks(:,:),Cprj_b2ks(:,:)
171 [ + + ]: 858 : type(libxc_functional_type) :: xc_funcs_hybrid(2)
172 : : ! *************************************************************************
173 : :
174 : : DBG_ENTER("COLL")
175 : :
176 [ - + + - : 1144 : 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 [ + + + + : 1135 : 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 : 286 : rank = Wfd%my_rank
185 : 286 : call initmpi_seq(MPI_enreg_seq)
186 : 286 : call MPI_enreg_seq%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
187 : :
188 : 286 : nspinor=Wfd%nspinor; nsppol =Wfd%nsppol; nspden =Wfd%nspden
189 [ + + ]: 286 : 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 [ - + + - : 858 : ABI_MALLOC(kcalc2ibz,(Wfd%nkibz))
+ - + - -
+ - + ]
195 [ + + ]: 1998 : kcalc2ibz=0
196 : :
197 : : ! Index in the IBZ of the GW k-points.
198 : : ! Only these points will be considered.
199 : 286 : nk_calc=0
200 [ + + ]: 1998 : do ik_ibz=1,Wfd%nkibz
201 [ + + + + : 4148 : if ( ALL(kstab(1,ik_ibz,:)/=0) .and. ALL(kstab(2,ik_ibz,:)/=0) ) then
+ + + + +
+ ]
202 : 1065 : nk_calc=nk_calc+1; kcalc2ibz(nk_calc) = ik_ibz
203 : : end if
204 : : end do
205 : :
206 : 286 : call Mels%init(Mflags, nsppol, nspden, Wfd%nspinor, Wfd%nkibz, Wfd%kibz, kstab)
207 : :
208 [ - + ]: 286 : 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 : 286 : call wrtout(std_out," calc_vhxc_braket: calculating v_xc[n_val] (excluding non-linear core corrections)")
214 : :
215 [ + + ]: 578 : do isppol=1,nsppol
216 [ + - - + : 5235865 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density rhor = ',MINVAL(rhor(:,isppol))
- - + + ]
217 : 292 : call wrtout(std_out, msg)
218 [ + + ]: 578 : if (Wfd%usepaw==1) then
219 [ + - - + : 223486 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density nhat = ',MINVAL(nhat(:,isppol))
- - + + ]
220 : 7 : call wrtout(std_out, msg)
221 [ + - - + : 223486 : write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density trho-nhat = ',MINVAL(rhor(:,isppol)-nhat(:,isppol))
- - + + ]
222 : 7 : call wrtout(std_out, msg)
223 : 7 : write(msg,'(a,i2)')' using usexcnhat = ',usexcnhat
224 : 7 : call wrtout(std_out, msg)
225 : : end if
226 : : end do
227 : :
228 : 286 : option = 0 ! Only exc, vxc
229 : 286 : nkxc = 0 ! No computation of XC kernel
230 : 286 : n3xccc_= 0 ! No core
231 : 286 : nk3xc = 0 ! k3xc not needed
232 [ + + + - ]: 286 : nmxc = Dtset%usepaw==1 .and. mod(abs(Dtset%usepawu),10) == 4
233 : :
234 [ - + ]: 286 : ABI_MALLOC(xccc3d_,(n3xccc_))
235 [ - + + - : 572 : ABI_MALLOC(kxc_,(nfftf,nkxc))
- + - + ]
236 [ - + + - : 1144 : ABI_MALLOC(vxc_val,(nfftf,nspden))
- + + - +
- + - - +
- + ]
237 : :
238 : 286 : 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 [ + - - + ]: 286 : usexcnhat,vxc_val,vxcval_avg,xccc3d_,xcdata,taur=taur)
243 : :
244 : : ! FABIEN's development
245 : : ! Hybrid functional treatment
246 [ + + ]: 286 : 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 [ - + ]: 286 : ABI_FREE(xccc3d_)
288 [ - + ]: 286 : ABI_FREE(kxc_)
289 : :
290 : 286 : write(msg,'(a,f8.4,2a,f8.4,a)')' E_xc[n_val] = ',bigexc_val, ' [Ha]. ','<V_xc[n_val]> = ',vxcval_avg,' [Ha]. '
291 : 286 : call wrtout(std_out, msg)
292 [ - + ]: 286 : 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 [ + + ]: 286 : 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 [ + + - + ]: 286 : 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 [ + + ]: 286 : 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 [ - + + - : 858 : ABI_MALLOC(ur1, (nfftf * nspinor))
+ - + - -
+ - + ]
340 [ - + + - : 572 : ABI_MALLOC(ur2, (nfftf * nspinor))
- + - + ]
341 [ - + + - : 572 : 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 [ + + ]: 726 : got=0; my_nmels=0
348 [ - + + - : 1716 : ABI_MALLOC(bbp_ks_distrb,(Wfd%mband,Wfd%mband,nk_calc,nsppol))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
349 [ + + ]: 578 : do is=1,nsppol
350 [ + + ]: 1653 : do ikc=1,nk_calc
351 : 1075 : ik_ibz=kcalc2ibz(ikc)
352 [ + + + + ]: 584723 : bbp_mask=.FALSE.
353 : 1075 : b_start=kstab(1,ik_ibz,is)
354 : 1075 : b_stop =kstab(2,ik_ibz,is)
355 [ + + ]: 1075 : if (Mflags%only_diago==1) then
356 : : !do jb=b1,b2
357 [ + + ]: 1724 : do jb=b_start,b_stop
358 : 1724 : 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 : 1075 : call wfd%distribute_bbp(ik_ibz,is,"Upper",my_nbbp,bbp_ks_distrb(:,:,ikc,is),got,bbp_mask)
365 : 1367 : my_nmels = my_nmels + my_nbbp
366 : : end do
367 : : end do
368 [ - + ]: 286 : ABI_FREE(bbp_mask)
369 : :
370 : 286 : write(msg,'(a,i0,a)')" Will calculate ",my_nmels," <b,k,s|O|b',k,s> matrix elements in calc_vhxc_me."
371 : 286 : call wrtout(std_out, msg)
372 : :
373 : : ! =====================================
374 : : ! ==== Loop over required k-points ====
375 : : ! =====================================
376 : 286 : nfftfm1=one/nfftf
377 : :
378 [ + + ]: 578 : do is=1,nsppol
379 [ + - + - : 15432 : if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE
+ + + + +
- ]
380 [ + + ]: 1653 : do ikc=1,nk_calc
381 [ + - + + : 20994 : if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE
+ + + - ]
382 : :
383 : 1075 : ik_ibz = kcalc2ibz(ikc)
384 : 1075 : b_start = kstab(1,ik_ibz,is)
385 : 1075 : b_stop = kstab(2,ik_ibz,is)
386 : 1075 : npw_k = Wfd%Kdata(ik_ibz)%npw
387 [ + + ]: 4300 : kpt = Wfd%kibz(:,ik_ibz)
388 : 1075 : kg_k => Wfd%kdata(ik_ibz)%kg_k
389 : 1075 : istwf_k = wfd%istwfk(ik_ibz)
390 : :
391 : : ! Calculate |k+G|^2 needed by hbareme and kineticme
392 : : ! MRM: Solved ecut problem
393 [ + + ]: 1075 : 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 [ + + ]: 10726 : do jb=b_start,b_stop
404 [ + + + + : 24049 : if (ALL(bbp_ks_distrb(:,jb,ikc,is)/=rank)) CYCLE
+ + ]
405 : :
406 [ - + ]: 9095 : ABI_CHECK(wfd%get_wave_ptr(jb, ik_ibz, is, wave_jb, msg) == 0, msg)
407 : :
408 [ + + ]: 9095 : 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 : 9095 : call wfd%get_ur(jb,ik_ibz,is,ur2)
415 : :
416 : : !do ib=b1,jb ! Upper triangle
417 [ + + ]: 63417 : do ib=b_start,jb
418 [ + + ]: 53247 : if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE
419 : : ! Off-diagonal elements only for QPSCGW.
420 [ + + + - ]: 44255 : if (Mflags%only_diago==1.and.ib/=jb) CYCLE
421 : :
422 : 44255 : call wfd%get_ur(ib,ik_ibz,is,ur1)
423 [ + + ]: 261478919 : u1cjg_u2dpc(:) = CONJG(ur1) *ur2
424 : :
425 [ + + ]: 44255 : if (Mflags%has_vxc == 1) then
426 [ + + ]: 137582690 : Mels%vxc(ib, jb, ik_ibz, is) = sum(u1cjg_u2dpc(1:nfftf) * vxc(1:nfftf, is)) * nfftfm1
427 [ + + ]: 23143 : 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 [ + + ]: 44255 : if (Mflags%has_vxcval == 1) then
431 [ + + ]: 137582690 : Mels%vxcval(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vxc_val(1:nfftf, is)) * nfftfm1
432 [ + + ]: 23143 : 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 [ + + ]: 44255 : 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 [ + - ]: 44255 : if (Mflags%has_vhartree==1) then
441 [ + + ]: 260093959 : Mels%vhartree(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vhartr(1:nfftf)) * nfftfm1
442 [ + + ]: 44255 : 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 [ + + ]: 44255 : 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 [ + + ]: 44255 : 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 [ + + + + ]: 53906 : 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 [ + + ]: 1367 : 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 [ - + ]: 286 : ABI_FREE(ur1)
525 [ - + ]: 286 : ABI_FREE(ur2)
526 [ - + ]: 286 : ABI_FREE(vxc_val)
527 [ - + ]: 286 : ABI_FREE(u1cjg_u2dpc)
528 [ + + ]: 286 : if(Mflags%has_vxcval_hybrid==1) then
529 [ - + ]: 23 : ABI_FREE(vxc_val_hybrid)
530 : : end if
531 [ + + ]: 286 : 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 [ + + ]: 286 : 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 [ + + ]: 286 : if (Wfd%usepaw==1) then
550 : : ! Tests if needed pointers in Paw_ij are allocated.
551 [ + - + - : 7 : 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 [ - + ]: 7 : ABI_CHECK(nspinor == 1, "PAW with nspinor not tested")
554 : :
555 : : ! For DFT+U
556 [ + + ]: 42 : do iat=1,Cryst%natom
557 : 35 : itypat=Cryst%typat(iat)
558 [ - + ]: 42 : 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 [ - + ]: 7 : 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 [ + + ]: 19 : lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
570 : :
571 [ + + ]: 7 : 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 [ + - + - : 2 : .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 : 7 : nsploop=nspinor**2
580 : :
581 : : ! ====================================
582 : : ! === Assemble PAW matrix elements ===
583 : : ! ====================================
584 [ - + + - : 21 : ABI_MALLOC(dimlmn, (Cryst%natom))
+ - + - -
+ - + ]
585 [ + + ]: 42 : do iat=1,Cryst%natom
586 : 42 : dimlmn(iat)=Pawtab(Cryst%typat(iat))%lmn_size
587 : : end do
588 : :
589 [ - + + - : 70 : ABI_MALLOC(Cprj_b1ks, (Cryst%natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
590 [ - + + - : 63 : ABI_MALLOC(Cprj_b2ks, (Cryst%natom,nspinor))
- + + - -
+ - + + +
+ + ]
591 : 7 : call pawcprj_alloc(Cprj_b1ks, 0, dimlmn)
592 : 7 : call pawcprj_alloc(Cprj_b2ks, 0, dimlmn)
593 : :
594 [ + + ]: 14 : do is=1,nsppol
595 [ + - + - : 2269 : if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE
+ + + + +
- ]
596 : :
597 : : ! Loop over required k-points
598 [ + + ]: 25 : do ikc=1,nk_calc
599 [ + - + + : 2273 : if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE
+ + + - ]
600 : 11 : ik_ibz=kcalc2ibz(ikc)
601 : 11 : b_start=kstab(1,ik_ibz,is)
602 : 11 : b_stop =kstab(2,ik_ibz,is)
603 : :
604 : : !do jb=b1,b2
605 [ + + ]: 72 : do jb=b_start,b_stop
606 [ + - + + : 369 : 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 : 54 : call wfd%get_cprj(jb,ik_ibz,is,Cryst,Cprj_b2ks,sorted=.FALSE.)
611 : :
612 : : !do ib=b1,jb ! Upper triangle
613 [ + + ]: 248 : do ib=b_start,jb
614 [ + + ]: 183 : if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE
615 : :
616 : : ! Off-diagonal elements only for QPSCGW.
617 [ + - + - ]: 54 : if (Mflags%only_diago==1.and.ib/=jb) CYCLE
618 : :
619 : 54 : 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 [ + + + + : 1296 : tmp_xc = zero; tmp_xcval = zero; tmp_H = zero; tmp_U = zero; tmp_h0ij = zero; tmp_sigcx = zero
+ + + + +
+ + + + +
+ + + + +
+ + + +
+ ]
625 : :
626 [ + + ]: 204 : do iat=1,Cryst%natom
627 : 150 : itypat =Cryst%typat(iat)
628 : 150 : lmn_size =Pawtab(itypat)%lmn_size
629 : 150 : cplex_dij=Paw_ij(iat)%cplex_dij
630 : 150 : klmn1=1
631 : :
632 [ + + ]: 1574 : do jlmn=1,lmn_size
633 : 1370 : j0lmn=jlmn*(jlmn-1)/2
634 [ + + ]: 8990 : do ilmn=1,jlmn
635 : 7470 : 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 [ + + ]: 7470 : fact=one; if (ilmn==jlmn) fact=half
639 : :
640 : : ! Loop over four components if nspinor==2
641 : : ! If collinear nsploop==1
642 [ + + ]: 14940 : do iab=1,nsploop
643 : 7470 : 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 : 7470 : +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 : 7470 : -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 [ + - ]: 14940 : if (nspinor==1) then
659 : :
660 [ - + ]: 7470 : 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 [ + + ]: 7470 : if (Mflags%has_sxcore==1) then ! * Fock operator generated by core electrons.
667 : 2268 : dijsigcx = dijexc_core(klmn,1,itypat)
668 : 2268 : tmp_sigcx(1,iab)=tmp_sigcx(1,iab) + dijsigcx*re_p*fact
669 : 2268 : tmp_sigcx(2,iab)=tmp_sigcx(2,iab) + dijsigcx*im_p*fact
670 : : end if
671 : :
672 [ + - ]: 7470 : if (Mflags%has_vxc==1) then ! * Accumulate vxc[n1+nc] + vxc[n1+tn+nc].
673 : 7470 : vxc1 = Paw_ij(iat)%dijxc(klmn,is)+Paw_ij(iat)%dijxc_hat(klmn,is)
674 : 7470 : tmp_xc(1,iab)=tmp_xc(1,iab) + vxc1*re_p*fact
675 : 7470 : tmp_xc(2,iab)=tmp_xc(2,iab) + vxc1*im_p*fact
676 : : end if
677 : :
678 [ + - ]: 7470 : if (Mflags%has_vxcval==1) then ! * Accumulate valence-only XC.
679 : 7470 : vxc1_val=Paw_ij(iat)%dijxc_val(klmn,is)
680 : 7470 : tmp_xcval(1,1)=tmp_xcval(1,1) + vxc1_val*re_p*fact
681 : 7470 : tmp_xcval(2,1)=tmp_xcval(2,1) + vxc1_val*im_p*fact
682 : : end if
683 : :
684 [ + - ]: 7470 : if (Mflags%has_vhartree==1) then ! * Accumulate Hartree term of the PAW Hamiltonian.
685 : 7470 : DijH=Paw_ij(iat)%dijhartree(klmn)
686 : 7470 : tmp_H(1,1)=tmp_H(1,1) + DijH*re_p*fact
687 : 7470 : 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 [ - + ]: 7470 : 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 : 8840 : 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 [ + - ]: 108 : if (nspinor==1) then
755 : :
756 [ - + ]: 54 : 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 [ + - ]: 54 : if (Mflags%has_vxc==1) &
760 : 54 : & 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 [ + - ]: 54 : if (Mflags%has_vxcval==1) &
763 : 54 : & 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 [ - + ]: 54 : 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 [ + - ]: 54 : if (Mflags%has_vhartree==1) &
769 : 54 : & 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 [ - + ]: 54 : 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 [ + + ]: 54 : if (Mflags%has_sxcore==1) &
775 : 24 : & 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 [ - + ]: 7 : ABI_FREE(dimlmn)
805 : 7 : call pawcprj_free(Cprj_b1ks)
806 [ + - + + : 42 : ABI_FREE(Cprj_b1ks)
- + - + -
+ ]
807 : 7 : call pawcprj_free(Cprj_b2ks)
808 [ + - + + : 42 : ABI_FREE(Cprj_b2ks)
- + - + -
+ ]
809 : : end if !PAW
810 : :
811 [ - + ]: 286 : ABI_FREE(bbp_ks_distrb)
812 : :
813 : : ! Sum up contributions on each node
814 : : ! Set the corresponding has_* flags to 2.
815 : 286 : call Mels%mpisum(wfd%comm)
816 : :
817 : : ! Reconstruct lower triangle.
818 : 286 : call Mels%herm()
819 : :
820 [ - + ]: 286 : ABI_FREE(kcalc2ibz)
821 : 286 : call destroy_mpi_enreg(MPI_enreg_seq)
822 : :
823 : : DBG_EXIT("COLL")
824 : :
825 [ + - - + : 1144 : end subroutine calc_vhxc_me
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - + - -
- - - - ]
826 : : !!***
827 : :
828 : : end module m_vhxc_me
829 : : !!***
|
