Branch data Line data Source code
1 : : !!****m* ABINIT/m_getchc
2 : : !! NAME
3 : : !! m_getchc
4 : : !!
5 : : !! FUNCTION
6 : : !! Compute <G|H|C> for input vector |C> expressed in reciprocal space;
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, LSI, MT)
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_getchc
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_xmpi
28 : :
29 : : use defs_abitypes, only : mpi_type
30 : : use m_time, only : timab
31 : : use m_pawcprj, only : pawcprj_type!, pawcprj_alloc, pawcprj_free, pawcprj_getdim
32 : : ! use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_get_ikpt
33 : : use m_hamiltonian, only : gs_hamiltonian_type, KPRIME_H_K, K_H_KPRIME, K_H_K, KPRIME_H_KPRIME
34 : : use m_nonlop, only : nonlop
35 : : ! use m_fock, only : fock_common_type, fock_get_getghc_call
36 : : ! use m_fock_getghc, only : fock_getghc, fock_ACE_getghc
37 : : use m_cgtools, only : dotprod_g
38 : :
39 : : implicit none
40 : :
41 : : private
42 : : !!***
43 : :
44 : : public :: getchc
45 : : public :: getcsc
46 : : !!***
47 : :
48 : : contains
49 : : !!***
50 : :
51 : : !!****f* ABINIT/getchc
52 : : !!
53 : : !! NAME
54 : : !! getchc
55 : : !!
56 : : !! FUNCTION
57 : : !! Compute <C_left|H|C> for input vectors |C> and |C_left>.
58 : : !! Note that |C_left> can be an array of "ndat" wavefunctions if only the non-local part is computed
59 : : !!
60 : : !! INPUTS
61 : : !! cpopt=flag defining the status of cwaveprj%cp(:)=<Proj_i|Cnk> scalars (PAW only)
62 : : !! (same meaning as in nonlop.F90 routine)
63 : : !! if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
64 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
65 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
66 : : !! if cpopt= 2 <p_lmn|in> are already in memory;
67 : : !! if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved
68 : : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
69 : : !! cwavef(2,npw*my_nspinor)=planewave coefficients of wavefunction.
70 : : !! cwavef_left(2,npw*my_nspinor)=planewave coefficients of wavefunction left.
71 : : !! cwaveprj(natom,my_nspinor*(1+cpopt))= wave function projected on nl projectors
72 : : !! cwaveprj_left(natom,my_nspinor*(1+cpopt))= wave function projected on nl projectors (for left WF)
73 : : !! cwavef_r(2,n4,n5,n6,nspinor) = wave function in real space
74 : : !! cwavef_left_r(2,n4,n5,n6,nspinor) = wave function in real space (for left WF)
75 : : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
76 : : !! lambda=factor to be used when computing <G|H-lambda.S|C> - only for sij_opt=-1
77 : : !! Typically lambda is the eigenvalue (or its guess)
78 : : !! mpi_enreg=information about MPI parallelization
79 : : !! ndat=number of left wavefunctions
80 : : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H|C> have to be computed
81 : : !! (S=overlap) if 1, matrix elements <G|S|C> have to be computed in gsc in addition to ghc
82 : : !! if -1, matrix elements <G|H-lambda.S|C> have to be computed in ghc (gsc not used)
83 : : !! type_calc= option governing which part of Hamitonian is to be applied:
84 : : ! 0: whole Hamiltonian
85 : : !! 1: local part only
86 : : !! 2: non-local+Fock+kinetic only (added to the existing Hamiltonian)
87 : : !! 3: local + kinetic only (added to the existing Hamiltonian)
88 : : !! ===== Optional inputs =====
89 : : !! [kg_fft_k(3,:)]=optional, (k+G) vector coordinates to be used for the FFT tranformation
90 : : !! instead of the one contained in gs_ham datastructure.
91 : : !! Typically used for real WF (in parallel) which are FFT-transformed 2 by 2.
92 : : !! [kg_fft_kp(3,:)]=optional, (k^prime+G) vector coordinates to be used for the FFT tranformation
93 : : !! [select_k]=optional, option governing the choice of k points to be used.
94 : : !! gs_ham datastructure contains quantities needed to apply Hamiltonian
95 : : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
96 : : !! if select_k=1, <k^prime|H|k> is applied [default]
97 : : !! if select_k=2, <k|H|k^prime> is applied
98 : : !! if select_k=3, <k|H|k> is applied
99 : : !! if select_k=4, <k^prime|H|k^prime> is applied
100 : : !!
101 : : !! OUTPUT
102 : : !! chc(2*ndat)=matrix elements <C_left|H|C> (if sij_opt>=0)
103 : : !! or <C_left|H-lambda.S|C> (if sij_opt=-1)
104 : : !! SIDE EFFECTS
105 : : !!
106 : : !! SOURCE
107 : :
108 [ # # # # : 0 : subroutine getchc(chc,cpopt,cwavef,cwavef_left,cwaveprj,cwaveprj_left,cwavef_r,cwavef_left_r,&
# # # # #
# # # ]
109 : : & gs_ham,lambda,mpi_enreg,ndat,&
110 : : & sij_opt,type_calc,&
111 [ # # # # : 0 : & kg_fft_k,kg_fft_kp,select_k) ! optional arguments
# # # # ]
112 : :
113 : : !Arguments ------------------------------------
114 : : !scalars
115 : : integer,intent(in) :: cpopt,ndat
116 : : integer,intent(in) :: sij_opt,type_calc
117 : : integer,intent(in),optional :: select_k
118 : : real(dp),intent(in) :: lambda
119 : : real(dp),intent(inout) :: chc(2*ndat)
120 : : type(MPI_type),intent(in) :: mpi_enreg
121 : : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
122 : : !arrays
123 : : integer,intent(in),optional,target :: kg_fft_k(:,:),kg_fft_kp(:,:)
124 : : real(dp),intent(inout) :: cwavef(:,:),cwavef_left(:,:),cwavef_r(:,:,:,:,:),cwavef_left_r(:,:,:,:,:)
125 : : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:),cwaveprj_left(:,:)
126 : :
127 : : !Local variables-------------------------------
128 : : !scalars
129 : : integer,parameter :: re=1,im=2
130 : : integer :: choice,cpopt_here,i1,i2,i3,idat,idir
131 : : integer :: ig,igspinor,ispinor,ispinor_left,my_nspinor
132 : : integer :: nnlout,nffttot,npw,npw_k1,npw_k2,nspinortot,n1,n2,n3
133 : : integer :: paw_opt,select_k_,shift1,shift2,signs,tim_nonlop
134 : : logical :: k1_eq_k2,has_fock
135 : : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
136 : : character(len=500) :: msg
137 : : !arrays
138 : 0 : integer, pointer :: gbound_k1(:,:),gbound_k2(:,:),kg_k1(:,:),kg_k2(:,:)
139 : 0 : real(dp) :: enlout(ndat),enlout_im(ndat),lambda_ndat(ndat),tsec(2),z_tmp(2)
140 : 0 : real(dp),allocatable :: gsc(:,:),gvnlxc(:,:)
141 : 0 : real(dp), pointer :: kinpw_k1(:),kinpw_k2(:),kpt_k1(:),kpt_k2(:)
142 : :
143 : : ! *********************************************************************
144 : :
145 : : DBG_ENTER("COLL")
146 : :
147 : : !Keep track of total time spent in getchc:
148 : 0 : call timab(1370,1,tsec)
149 : :
150 : : !Select k-dependent objects according to select_k input parameter
151 [ # # ]: 0 : select_k_=1;if (present(select_k)) select_k_=select_k
152 [ # # ]: 0 : if (select_k_==KPRIME_H_K) then
153 : : ! <k^prime|H|k>
154 : 0 : npw_k1 = gs_ham%npw_fft_k ; npw_k2 = gs_ham%npw_fft_kp
155 : 0 : kpt_k1 => gs_ham%kpt_k ; kpt_k2 => gs_ham%kpt_kp
156 : 0 : kg_k1 => gs_ham%kg_k ; kg_k2 => gs_ham%kg_kp
157 : 0 : gbound_k1 => gs_ham%gbound_k ; gbound_k2 => gs_ham%gbound_kp
158 : 0 : kinpw_k1 => gs_ham%kinpw_k ; kinpw_k2 => gs_ham%kinpw_kp
159 [ # # ]: 0 : else if (select_k_==K_H_KPRIME) then
160 : : ! <k|H|k^prime>
161 : 0 : npw_k1 = gs_ham%npw_fft_kp; npw_k2 = gs_ham%npw_fft_k
162 : 0 : kpt_k1 => gs_ham%kpt_kp ; kpt_k2 => gs_ham%kpt_k
163 : 0 : kg_k1 => gs_ham%kg_kp ; kg_k2 => gs_ham%kg_k
164 : 0 : gbound_k1 => gs_ham%gbound_kp ; gbound_k2 => gs_ham%gbound_k
165 : 0 : kinpw_k1 => gs_ham%kinpw_kp ; kinpw_k2 => gs_ham%kinpw_k
166 [ # # ]: 0 : else if (select_k_==K_H_K) then
167 : : ! <k|H|k>
168 : 0 : npw_k1 = gs_ham%npw_fft_k ; npw_k2 = gs_ham%npw_fft_k
169 : 0 : kpt_k1 => gs_ham%kpt_k ; kpt_k2 => gs_ham%kpt_k
170 : 0 : kg_k1 => gs_ham%kg_k ; kg_k2 => gs_ham%kg_k
171 : 0 : gbound_k1 => gs_ham%gbound_k ; gbound_k2 => gs_ham%gbound_k
172 : 0 : kinpw_k1 => gs_ham%kinpw_k ; kinpw_k2 => gs_ham%kinpw_k
173 [ # # ]: 0 : else if (select_k_==KPRIME_H_KPRIME) then
174 : : ! <k^prime|H|k^prime>
175 : 0 : npw_k1 = gs_ham%npw_fft_kp; npw_k2 = gs_ham%npw_fft_kp
176 : 0 : kpt_k1 => gs_ham%kpt_kp ; kpt_k2 => gs_ham%kpt_kp
177 : 0 : kg_k1 => gs_ham%kg_kp ; kg_k2 => gs_ham%kg_kp
178 : 0 : gbound_k1 => gs_ham%gbound_kp ; gbound_k2 => gs_ham%gbound_kp
179 : 0 : kinpw_k1 => gs_ham%kinpw_kp ; kinpw_k2 => gs_ham%kinpw_kp
180 : : end if
181 [ # # # # ]: 0 : k1_eq_k2=(all(abs(kpt_k1(:)-kpt_k2(:))<tol8))
182 [ # # ]: 0 : if (.not.k1_eq_k2) then
183 : 0 : ABI_ERROR('getchc is not implemented yet for k1/=k2')
184 : : end if
185 : :
186 : : !Check sizes
187 : 0 : my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor)
188 [ # # # # ]: 0 : if (size(cwavef)<2*npw_k1*my_nspinor) then
189 : 0 : msg='wrong size for cwavef!'
190 : 0 : ABI_BUG(msg)
191 : : end if
192 [ # # # # ]: 0 : if (gs_ham%usepaw==1.and.cpopt>=0) then
193 [ # # # # ]: 0 : if (size(cwaveprj)<gs_ham%natom*my_nspinor) then
194 : 0 : msg='wrong size for cwaveprj!'
195 : 0 : ABI_BUG(msg)
196 : : end if
197 : : end if
198 [ # # ]: 0 : if (gs_ham%usepaw==1) then
199 [ # # # # ]: 0 : if (size(cwaveprj_left)<gs_ham%natom*my_nspinor*ndat) then
200 : 0 : msg='wrong size for cwaveprj_left!'
201 : 0 : ABI_BUG(msg)
202 : : end if
203 : : end if
204 : :
205 : : !Eventually overwrite plane waves data for FFT
206 [ # # ]: 0 : if (present(kg_fft_k)) then
207 : 0 : kg_k1 => kg_fft_k ; kg_k2 => kg_fft_k
208 : 0 : npw_k1=size(kg_k1,2) ; npw_k2=size(kg_k2,2)
209 : : end if
210 [ # # ]: 0 : if (present(kg_fft_kp)) then
211 : 0 : kg_k2 => kg_fft_kp ; npw_k2=size(kg_k2,2)
212 : : end if
213 : :
214 : : !paral_kgb constraint
215 [ # # # # ]: 0 : if (mpi_enreg%paral_kgb==1.and.(.not.k1_eq_k2)) then
216 : 0 : msg='paral_kgb=1 not allowed for k/=k_^prime!'
217 : 0 : ABI_BUG(msg)
218 : : end if
219 : :
220 : : !Do we add Fock exchange term ?
221 : 0 : has_fock=(associated(gs_ham%fockcommon))
222 [ # # ]: 0 : if (has_fock) then
223 : 0 : ABI_BUG('Fock not implemented yet')
224 : : end if
225 : : ! if (has_fock) fock => gs_ham%fockcommon
226 : :
227 : : !Parallelization over spinors management
228 : : if (mpi_enreg%paral_spinor==0) then
229 : : shift1=npw_k1;shift2=npw_k2
230 : : nspinor1TreatedByThisProc=.true.
231 : : nspinor2TreatedByThisProc=(nspinortot==2)
232 : : else
233 : 0 : shift1=0;shift2=0
234 : 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
235 : 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
236 : : end if
237 : :
238 : 0 : npw=gs_ham%npw_k
239 : 0 : nspinortot=gs_ham%nspinor
240 : :
241 : : !============================================================
242 : : ! Application of the local potential
243 : : !============================================================
244 : :
245 [ # # ]: 0 : if ((type_calc==0).or.(type_calc==1).or.(type_calc==3)) then
246 : :
247 : 0 : call timab(1371,1,tsec)
248 : : ! Need a Vlocal
249 [ # # ]: 0 : if (.not.associated(gs_ham%vlocal)) then
250 : 0 : ABI_BUG("We need vlocal in gs_ham!")
251 : : end if
252 [ # # ]: 0 : if (ndat>1) then
253 : 0 : ABI_ERROR("ndat should be 1 for the local part")
254 : : end if
255 : :
256 : 0 : n1=gs_ham%ngfft(1)
257 : 0 : n2=gs_ham%ngfft(2)
258 : 0 : n3=gs_ham%ngfft(3)
259 : 0 : nffttot=n1*n2*n3
260 [ # # ]: 0 : chc = zero
261 : : ! Treat scalar local potentials
262 [ # # ]: 0 : if (gs_ham%nvloc==1) then
263 [ # # ]: 0 : if (gs_ham%istwf_k==2) then
264 [ # # ]: 0 : do i3=1,n3
265 [ # # ]: 0 : do i2=1,n2
266 [ # # ]: 0 : do i1=1,n1
267 : 0 : chc(1) = chc(1) + gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(1,i1,i2,i3,1)*cwavef_left_r(1,i1,i2,i3,1)
268 : : end do
269 : : end do
270 : : end do
271 : 0 : chc(2)=zero
272 : : else
273 [ # # ]: 0 : do ispinor=1,my_nspinor
274 [ # # ]: 0 : do i3=1,n3
275 [ # # ]: 0 : do i2=1,n2
276 [ # # ]: 0 : do i1=1,n1
277 : : z_tmp(1) = cwavef_r(1,i1,i2,i3,ispinor)*cwavef_left_r(1,i1,i2,i3,ispinor) &
278 : 0 : & +cwavef_r(2,i1,i2,i3,ispinor)*cwavef_left_r(2,i1,i2,i3,ispinor)
279 : : z_tmp(2) = cwavef_r(2,i1,i2,i3,ispinor)*cwavef_left_r(1,i1,i2,i3,ispinor) &
280 : 0 : & -cwavef_r(1,i1,i2,i3,ispinor)*cwavef_left_r(2,i1,i2,i3,ispinor)
281 : 0 : chc(1) = chc(1) + gs_ham%vlocal(i1,i2,i3,1)*z_tmp(1)
282 : 0 : chc(2) = chc(2) + gs_ham%vlocal(i1,i2,i3,1)*z_tmp(2)
283 : : end do
284 : : end do
285 : : end do
286 : : end do
287 : : end if
288 : : else ! nvloc = 4
289 [ # # ]: 0 : do ispinor=1,my_nspinor
290 [ # # ]: 0 : do ispinor_left=1,my_nspinor
291 [ # # ]: 0 : do i3=1,n3
292 [ # # ]: 0 : do i2=1,n2
293 [ # # ]: 0 : do i1=1,n1
294 : : z_tmp(1) = cwavef_r(1,i1,i2,i3,ispinor)*cwavef_left_r(1,i1,i2,i3,ispinor_left) &
295 : 0 : & +cwavef_r(2,i1,i2,i3,ispinor)*cwavef_left_r(2,i1,i2,i3,ispinor_left)
296 : : z_tmp(2) = cwavef_r(2,i1,i2,i3,ispinor)*cwavef_left_r(1,i1,i2,i3,ispinor_left) &
297 : 0 : & -cwavef_r(1,i1,i2,i3,ispinor)*cwavef_left_r(2,i1,i2,i3,ispinor_left)
298 [ # # ]: 0 : if (ispinor==ispinor_left) then
299 : : ! Then vloc is real : vloc_uu = vloc(1) and vloc_dd = vloc(2)
300 : 0 : chc(1) = chc(1) + gs_ham%vlocal(i1,i2,i3,ispinor)*z_tmp(1)
301 : 0 : chc(2) = chc(2) + gs_ham%vlocal(i1,i2,i3,ispinor)*z_tmp(2)
302 [ # # ]: 0 : else if (ispinor==1.and.ispinor_left==2) then ! Psi(left)_d^* Psi_u vloc_ud
303 : : ! Otherwise vloc is complex Re(vloc_ud) = vloc(3)
304 : : ! Im(vloc_ud) = vloc(4)
305 : : ! vloc_du = (vloc_ud)^*
306 : 0 : chc(1) = chc(1) + gs_ham%vlocal(i1,i2,i3,3)*z_tmp(1) + gs_ham%vlocal(i1,i2,i3,4)*z_tmp(2)
307 : 0 : chc(2) = chc(2) + gs_ham%vlocal(i1,i2,i3,3)*z_tmp(1) - gs_ham%vlocal(i1,i2,i3,4)*z_tmp(2)
308 [ # # ]: 0 : else if (ispinor==2.and.ispinor_left==1) then ! Psi(left)_u^* Psi_d vloc_ud
309 : 0 : chc(1) = chc(1) + gs_ham%vlocal(i1,i2,i3,3)*z_tmp(1) - gs_ham%vlocal(i1,i2,i3,4)*z_tmp(2)
310 : 0 : chc(2) = chc(2) + gs_ham%vlocal(i1,i2,i3,3)*z_tmp(1) + gs_ham%vlocal(i1,i2,i3,4)*z_tmp(2)
311 : : end if
312 : : end do
313 : : end do
314 : : end do
315 : : end do
316 : : end do
317 : : end if
318 [ # # ]: 0 : chc = chc / dble(nffttot)
319 : :
320 : 0 : call timab(1371,2,tsec)
321 : :
322 : : end if ! type_calc
323 : :
324 [ # # # # : 0 : if ((type_calc==0).or.(type_calc==2).or.(type_calc==3).or.(type_calc==4)) then
# # ]
325 : :
326 : : !============================================================
327 : : ! Application of the non-local potential and the Fock potential
328 : : !============================================================
329 : :
330 : : if ((type_calc==0).or.(type_calc==2).or.(type_calc==4)) then
331 : :
332 : 0 : signs=1 ; choice=1 ; nnlout=1 ; idir=0 ; tim_nonlop=15
333 [ # # ]: 0 : cpopt_here=-1;if (gs_ham%usepaw==1) cpopt_here=cpopt
334 [ # # # # ]: 0 : ABI_MALLOC(gvnlxc,(0,0))
335 [ # # # # ]: 0 : ABI_MALLOC(gsc,(0,0))
336 : : ! if (has_fock) then
337 : : ! if (gs_ham%usepaw==1) then
338 : : ! cpopt_here=max(cpopt,0)
339 : : ! if (cpopt<2) then
340 : : ! ABI_DATATYPE_ALLOCATE(cwaveprj_fock,(gs_ham%natom,my_nspinor*ndat))
341 : : ! ABI_ALLOCATE(dimcprj,(gs_ham%natom))
342 : : ! call pawcprj_getdim(dimcprj,gs_ham%natom,gs_ham%nattyp,gs_ham%ntypat,&
343 : : !& gs_ham%typat,fock%pawtab,'O')
344 : : ! call pawcprj_alloc(cwaveprj_fock,0,dimcprj)
345 : : ! ABI_DEALLOCATE(dimcprj)
346 : : ! else
347 : : ! cwaveprj_fock=>cwaveprj
348 : : ! end if
349 : : ! cwaveprj_nonlop=>cwaveprj_fock
350 : : ! else
351 : : ! cwaveprj_nonlop=>cwaveprj
352 : : ! cwaveprj_fock=>cwaveprj
353 : : ! end if
354 : : ! else
355 : : ! end if
356 [ # # ]: 0 : paw_opt=gs_ham%usepaw ; if (sij_opt/=0) paw_opt=sij_opt+3
357 [ # # ]: 0 : lambda_ndat = lambda
358 : :
359 [ # # ]: 0 : enlout=zero
360 [ # # ]: 0 : enlout_im=zero
361 : : call nonlop(choice,cpopt_here,cwaveprj,enlout,gs_ham,idir,lambda_ndat,mpi_enreg,1,&
362 : : & nnlout,paw_opt,signs,gsc,tim_nonlop,cwavef,gvnlxc,select_k=select_k_,&
363 : 0 : & cprjin_left=cwaveprj_left,enlout_im=enlout_im,ndat_left=ndat)
364 : :
365 [ # # ]: 0 : do idat=1,ndat
366 : 0 : chc(2*idat-1) = chc(2*idat-1) + enlout(idat)
367 : 0 : chc(2*idat ) = chc(2*idat ) + enlout_im(idat)
368 : : end do
369 : :
370 [ # # ]: 0 : ABI_FREE(gvnlxc)
371 [ # # ]: 0 : ABI_FREE(gsc)
372 : :
373 : : end if ! if(type_calc...
374 : :
375 : : !============================================================
376 : : ! Assemble kinetic, local, nonlocal and Fock contributions
377 : : !============================================================
378 : :
379 [ # # # # ]: 0 : if (type_calc==0.or.type_calc==2.or.type_calc==3) then
380 : :
381 [ # # ]: 0 : if (ndat>1) then
382 : 0 : ABI_ERROR("ndat should be 1 for the kinetic part")
383 : : end if
384 : :
385 : 0 : call timab(1372,1,tsec)
386 : : ! Add modified kinetic contributions
387 : : ! to <CP|H|C(n,k)>.
388 [ # # ]: 0 : if (gs_ham%istwf_k==1) then
389 : : ! !!$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2)
390 [ # # ]: 0 : do ispinor=1,my_nspinor
391 [ # # ]: 0 : do ig=1,npw_k2
392 : 0 : igspinor=ig+npw_k2*(ispinor-1)
393 [ # # ]: 0 : if (kinpw_k2(ig)<huge(zero)*1.d-11) then
394 : 0 : chc(1) = chc(1) + kinpw_k2(ig)*cwavef(re,igspinor)*cwavef_left(re,igspinor)
395 : 0 : chc(1) = chc(1) + kinpw_k2(ig)*cwavef(im,igspinor)*cwavef_left(im,igspinor)
396 : 0 : chc(2) = chc(2) + kinpw_k2(ig)*cwavef(im,igspinor)*cwavef_left(re,igspinor)
397 : 0 : chc(2) = chc(2) - kinpw_k2(ig)*cwavef(re,igspinor)*cwavef_left(im,igspinor)
398 : : end if
399 : : end do ! ig
400 : : end do ! ispinor
401 [ # # # # ]: 0 : else if (gs_ham%istwf_k==2.and.mpi_enreg%me_g0==1) then
402 [ # # ]: 0 : if (kinpw_k2(1)<huge(zero)*1.d-11) then
403 : 0 : chc(1) = chc(1) + kinpw_k2(1)*cwavef(re,1)*cwavef_left(re,1)
404 : : end if
405 [ # # ]: 0 : do ig=2,npw_k2
406 [ # # ]: 0 : if (kinpw_k2(ig)<huge(zero)*1.d-11) then
407 : 0 : chc(1) = chc(1) + 2*kinpw_k2(ig)*cwavef(re,ig)*cwavef_left(re,ig)
408 : 0 : chc(1) = chc(1) + 2*kinpw_k2(ig)*cwavef(im,ig)*cwavef_left(im,ig)
409 : : end if
410 : : end do
411 : : else
412 [ # # ]: 0 : do ig=1,npw_k2
413 [ # # ]: 0 : if (kinpw_k2(ig)<huge(zero)*1.d-11) then
414 : 0 : chc(1) = chc(1) + 2*kinpw_k2(ig)*cwavef(re,ig)*cwavef_left(re,ig)
415 : 0 : chc(1) = chc(1) + 2*kinpw_k2(ig)*cwavef(im,ig)*cwavef_left(im,ig)
416 : : end if
417 : : end do ! ig
418 : : end if
419 : : ! Special case of PAW + Fock : only return Fock operator contribution in gvnlxc
420 : : ! if (gs_ham%usepaw==1 .and. has_fock)then
421 : : ! gvnlxc=gvnlxc-gvnlc
422 : : ! ABI_DEALLOCATE(gvnlc)
423 : : ! endif
424 : : !
425 : : ! if ((type_calc==0).or.(type_calc==2)) then
426 : : ! if (has_fock.and.gs_ham%usepaw==1.and.cpopt<2) then
427 : : ! call pawcprj_free(cwaveprj_fock)
428 : : ! ABI_DATATYPE_DEALLOCATE(cwaveprj_fock)
429 : : ! end if
430 : : ! end if
431 : 0 : call timab(1372,2,tsec)
432 : : end if
433 : :
434 : : end if ! type_calc
435 : :
436 : 0 : call timab(1370,2,tsec)
437 : :
438 : : DBG_EXIT("COLL")
439 : :
440 [ # # # # ]: 0 : end subroutine getchc
441 : : !!***
442 : :
443 : : !----------------------------------------------------------------------
444 : :
445 : : !!****f* ABINIT/getcsc
446 : : !! NAME
447 : : !! getcsc
448 : : !!
449 : : !! FUNCTION
450 : : !! Compute <C_left|S|C> for input vectors |C> and |C_left>.
451 : : !! Note that |C_left> can be an array of "ndat" wavefunctions
452 : : !!
453 : : !! INPUTS
454 : : !! cpopt=flag defining the status of cwaveprj%cp(:)=<Proj_i|Cnk> scalars (PAW only)
455 : : !! (same meaning as in nonlop.F90 routine)
456 : : !! if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
457 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
458 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
459 : : !! if cpopt= 2 <p_lmn|in> are already in memory;
460 : : !! if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved
461 : : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
462 : : !! cwavef(2,npw*my_nspinor)=planewave coefficients of wavefunction.
463 : : !! cwavef_left(2,npw*my_nspinor)=planewave coefficients of wavefunction left.
464 : : !! cprj(natom,my_nspinor*(1+cpopt))= wave function projected on nl projectors
465 : : !! cprj_left(natom,my_nspinor*(1+cpopt))= wave function projected on nl projectors (for left WF)
466 : : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
467 : : !! lambda=factor to be used when computing <G|H-lambda.S|C> - only for sij_opt=-1
468 : : !! Typically lambda is the eigenvalue (or its guess)
469 : : !! mpi_enreg=information about MPI parallelization
470 : : !! ndat=number of left wavefunctions
471 : : !! mpi_enreg=information about MPI parallelization
472 : : !!
473 : : !! OUTPUT
474 : : !! csc(2*ndat)=matrix elements <C_left|S|C>
475 : : !!
476 : : !! SOURCE
477 : :
478 [ # # # # : 0 : subroutine getcsc(csc,cpopt,cwavef,cwavef_left,cprj,cprj_left,gs_ham,mpi_enreg,ndat,&
# # # # ]
479 : : & tim_getcsc,&
480 : : & select_k) ! optional arguments
481 : :
482 : : !Arguments ------------------------------------
483 : : !scalars
484 : : integer,intent(in) :: cpopt,ndat,tim_getcsc
485 : : integer,intent(in),optional :: select_k
486 : : real(dp),intent(out) :: csc(2*ndat)
487 : : type(MPI_type),intent(in) :: mpi_enreg
488 : : type(gs_hamiltonian_type),intent(inout) :: gs_ham
489 : : !arrays
490 : : real(dp),intent(inout) :: cwavef(:,:)
491 : : real(dp),intent(inout),target :: cwavef_left(:,:)
492 : : type(pawcprj_type),intent(inout) :: cprj(:,:)
493 : : type(pawcprj_type),intent(inout),target :: cprj_left(:,:)
494 : :
495 : : !Local variables-------------------------------
496 : : !scalars
497 : : integer :: choice,idat,idir,istwf_k
498 : : integer :: npw,nspinor,paw_opt,select_k_,signs,tim_nonlop,nnlout
499 : : !character(len=500) :: msg
500 : : !arrays
501 : 0 : real(dp) :: tsec(2),real_csc(ndat)
502 : 0 : real(dp),allocatable :: gsc(:,:),gvnlxc(:,:)
503 : 0 : real(dp),allocatable :: enlout(:),enlout_im(:)
504 : : ! real(dp) :: dum
505 : : ! real(dp),pointer :: cwavef_left_idat(:,:)
506 : : ! *********************************************************************
507 : :
508 : : DBG_ENTER("COLL")
509 : :
510 : 0 : call timab(1360+tim_getcsc,1,tsec)
511 : :
512 : 0 : istwf_k = gs_ham%istwf_k
513 : 0 : npw = gs_ham%npw_k
514 : 0 : nspinor = gs_ham%nspinor
515 : :
516 [ # # ]: 0 : if (size(cwavef,2)/=npw*nspinor) then
517 : 0 : ABI_BUG('Wrong size for cwavef')
518 : : end if
519 [ # # ]: 0 : if (size(cwavef_left,2)/=npw*nspinor*ndat) then
520 : 0 : ABI_BUG('Wrong size for cwavef_left')
521 : : end if
522 : :
523 : 0 : call timab(1361,1,tsec)
524 [ # # ]: 0 : if (istwf_k==1) then
525 [ # # # # : 0 : call zgemv('C',npw*nspinor,ndat,cone,cwavef_left,npw*nspinor,cwavef,1,czero,csc,1)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
526 : : else ! nspinor==1 in that case
527 [ # # # # : 0 : call dgemv('C',2*npw,ndat,one,cwavef_left,2*npw,cwavef,1,zero,real_csc,1)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
528 [ # # ]: 0 : do idat=1,ndat
529 : 0 : csc(2*idat-1) = two*real_csc(idat)
530 : 0 : csc(2*idat ) = zero
531 : : end do
532 [ # # # # ]: 0 : if (istwf_k==2 .and. mpi_enreg%me_g0==1) then ! Gamma k-point and I have G=0
533 [ # # ]: 0 : do idat=1,ndat
534 : 0 : csc(2*idat-1) = csc(2*idat-1) - cwavef_left(1,1+npw*(idat-1))*cwavef(1,1)
535 : : end do
536 : : end if
537 : : end if
538 : 0 : call timab(1361,2,tsec)
539 : :
540 : :
541 [ # # ]: 0 : if (gs_ham%usepaw==1) then
542 : :
543 [ # # ]: 0 : if (size(cprj,2)/=nspinor) then
544 : 0 : ABI_BUG('Wrong size for cprj')
545 : : end if
546 [ # # ]: 0 : if (size(cprj_left,2)/=nspinor*ndat) then
547 : 0 : ABI_BUG('Wrong size for cprj_left')
548 : : end if
549 : :
550 [ # # ]: 0 : select_k_=1;if (present(select_k)) select_k_=select_k
551 : 0 : choice=1 ; nnlout=1 ; idir=0 ; tim_nonlop=16 ; paw_opt=3
552 [ # # # # ]: 0 : ABI_MALLOC(gsc,(0,0))
553 [ # # # # ]: 0 : ABI_MALLOC(gvnlxc,(0,0))
554 [ # # # # : 0 : ABI_MALLOC(enlout ,(ndat))
# # # # #
# # # #
# ]
555 [ # # # # : 0 : ABI_MALLOC(enlout_im,(ndat))
# # # # ]
556 [ # # ]: 0 : enlout=zero
557 [ # # ]: 0 : enlout_im=zero
558 : 0 : signs=1
559 : : call nonlop(choice,cpopt,cprj,enlout,gs_ham,idir,(/zero/),mpi_enreg,1,&
560 : : & nnlout,paw_opt,signs,gsc,tim_nonlop,cwavef,gvnlxc,select_k=select_k_,&
561 : 0 : & cprjin_left=cprj_left,enlout_im=enlout_im,ndat_left=ndat)
562 [ # # ]: 0 : do idat=1,ndat
563 : 0 : csc(2*idat-1) = csc(2*idat-1) + enlout(idat)
564 : 0 : csc(2*idat ) = csc(2*idat ) + enlout_im(idat)
565 : : end do
566 [ # # ]: 0 : ABI_FREE(gsc)
567 [ # # ]: 0 : ABI_FREE(gvnlxc)
568 [ # # ]: 0 : ABI_FREE(enlout )
569 [ # # ]: 0 : ABI_FREE(enlout_im)
570 : :
571 : : end if
572 : :
573 : 0 : call timab(1360+tim_getcsc,2,tsec)
574 : :
575 : : DBG_EXIT("COLL")
576 : :
577 [ # # # # : 0 : end subroutine getcsc
# # # # ]
578 : : !!***
579 : :
580 : : end module m_getchc
581 : : !!***
|
