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