Line data Source code
1 : !!****m* ABINIT/m_getgh2c
2 : !! NAME
3 : !! m_getgh2c
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2015-2026 ABINIT group (MT,JLJ)
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 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_getgh2c
26 :
27 : use defs_basis
28 : use m_abicore
29 : use m_errors
30 :
31 : use defs_abitypes, only : mpi_type
32 : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_free
33 : use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type
34 : use m_nonlop, only : nonlop
35 :
36 : #if defined(HAVE_GPU_MARKERS)
37 : use m_nvtx_data
38 : #endif
39 :
40 : implicit none
41 :
42 : private
43 : !!***
44 :
45 : public :: getgh2c
46 : !!***
47 :
48 : contains
49 : !!***
50 :
51 : !!****f* ABINIT/getgh2c
52 : !! NAME
53 : !! getgh2c
54 : !!
55 : !! FUNCTION
56 : !! Compute <G|H^(2)|C> (or <G|H^(2)-Eps.S^(2)|C>) for input vector |C> expressed in reciprocal space.
57 : !! (H^(2) is the 2nd-order pertubed Hamiltonian, S^(2) is the 2nd-order perturbed overlap operator).
58 : !! Result is put in array gh2c.
59 : !! If required, part of <G|K(2)+Vnonlocal^(2)|C> not depending on VHxc^(2) is also returned in gvnl2.
60 : !! If required, <G|S^(2)|C> is returned in gs2c (S=overlap - PAW only)
61 : !! Available for the following cases :
62 : !! ipert = natom+10 (dkdk) : 2nd derivative w.r.t wavevector
63 : !! natom+11 (dkdE) : mixed 2nd derivative w.r.t wavector and eletric field
64 : !! also if natom+12<=ipert<=2*natom+11 :
65 : !! (dtaudE) : mixed 2nd derivative w.r.t atom. displ. and eletric field (nonlocal only)
66 : !!
67 : !! INPUTS
68 : !! cwavef(2,npw*nspinor)=input wavefunction, in reciprocal space
69 : !! cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C>
70 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian
71 : !! idir=direction of the perturbation
72 : !! ipert=type of the perturbation
73 : !! lambda=real use to apply H^(2)-lambda.S^(2)
74 : !! mpi_enreg=information about MPI parallelization
75 : !! optlocal=0: local part of H^(2) is not computed in gh2c=<G|H^(2)|C>
76 : !! 1: local part of H^(2) is computed in gh2c=<G|H^(2)|C>
77 : !! optnl=0: non-local part of H^(2) is not computed in gh2c=<G|H^(2)|C>
78 : !! 1: non-local part of H^(2) depending on VHxc^(2) is not computed in gh2c=<G|H^(2)|C>
79 : !! 2: non-local part of H^(2) is totally computed in gh2c=<G|H^(2)|C>
80 : !! opt_gvnl2=option controlling the use of gvnl2 array:
81 : !! 0: not used
82 : !! 1: used as input: - used only for PAW and ipert=natom+11/+12
83 : !! At input: contains the derivative w.r.t wavevector of cwavef (times i)
84 : !! rf_hamkq <type(rf_hamiltonian_type)>=all data for the 2nd-order Hamiltonian at k,k+q
85 : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H^(2)|C> have to be computed
86 : !! (S=overlap) if 1, matrix elements <G|S^(2)|C> have to be computed in gs2c in addition to gh2c
87 : !! if -1, matrix elements <G|H^(2)-lambda.S^(2)|C> have to be computed in gh2c (gs2c not used)
88 : !! tim_getgh2c=timing code of the calling subroutine (can be set to 0 if not attributed)
89 : !! usevnl=1 if gvnl2=(part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2)) has to be input/output
90 : !!
91 : !! OUTPUT
92 : !! gh2c(2,npw1*nspinor)= <G|H^(2)|C> or <G|H^(2)-lambda.S^(2)|C>
93 : !! (only kinetic+non-local parts if optlocal=0)
94 : !! if (usevnl==1)
95 : !! gvnl2(2,npw1*nspinor)= part of <G|K^(2)+Vnl^(2)|C> not depending on VHxc^(2) (sij_opt/=-1)
96 : !! or part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2) (sij_opt==-1)
97 : !! if (sij_opt=1)
98 : !! gs2c(2,npw1*nspinor)=<G|S^(2)|C> (S=overlap).
99 : !!
100 : !! SOURCE
101 :
102 72096 : subroutine getgh2c(cwavef,cwaveprj,gh2c,gs2c,gs_hamkq,gvnl2,idir,ipert,lambda,&
103 23040 : & mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamkq,sij_opt,tim_getgh2c,usevnl,conj,enl,optkin)
104 :
105 : !Arguments ------------------------------------
106 : !scalars
107 : logical,intent(in),optional :: conj
108 : integer,intent(in) :: idir,ipert,optlocal,optnl,opt_gvnl2,sij_opt,tim_getgh2c,usevnl
109 : integer,intent(in),optional :: optkin
110 : real(dp),intent(in) :: lambda
111 : type(MPI_type),intent(in) :: mpi_enreg
112 : type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
113 : type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq
114 : !arrays
115 : real(dp),intent(in),optional,target :: enl(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq)
116 : real(dp),intent(inout) :: cwavef(:,:)
117 : real(dp),intent(inout),target :: gvnl2(:,:)
118 : real(dp),intent(out) :: gh2c(:,:),gs2c(:,:)
119 : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
120 :
121 : !Local variables-------------------------------
122 : !scalars
123 : integer,parameter :: tim_nonlop=0
124 : integer :: choice,cpopt,iatm,idir1,idir2,idirc,ipw,ipws,ispinor,my_nspinor
125 : integer :: natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs,usecprj
126 : logical :: compute_conjugate,has_kin,has_vnl,pert_phon_elfd
127 : real(dp) :: enlout_dum(1)
128 : character(len=500) :: msg
129 : !arrays
130 : ! integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
131 : integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/)
132 : !real(dp) :: tsec(2)
133 : real(dp) :: svectout_dum(1,1),vectout_dum(1,1)
134 72096 : real(dp),allocatable :: nonlop_out(:,:)
135 72096 : real(dp), pointer :: gvnl2_(:,:)
136 72096 : real(dp), pointer :: ddkinpw(:),kinpw1(:),enl_ptr(:,:,:,:)
137 72096 : real(dp),allocatable,target :: enl_temp(:,:,:,:)
138 72096 : type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:)
139 72096 : type(pawcprj_type),pointer :: cwaveprj_ptr(:,:)
140 :
141 : ! *********************************************************************
142 :
143 : DBG_ENTER("COLL")
144 : ABI_UNUSED(tim_getgh2c)
145 : ABI_NVTX_START_RANGE(NVTX_GETGH2C)
146 :
147 : !Keep track of total time spent in getgh2c
148 : !call timab(196+tim_getgh2c,1,tsec)
149 :
150 : !======================================================================
151 : !== Initialisations and compatibility tests
152 : !======================================================================
153 :
154 72096 : npw =gs_hamkq%npw_k
155 72096 : npw1 =gs_hamkq%npw_kp
156 72096 : natom=gs_hamkq%natom
157 :
158 : !Compatibility tests
159 25920 : if(ipert/=natom+10.and.ipert/=natom+11.and.ipert>2*natom+11)then
160 0 : msg='only ipert==natom+10/+11 and natom+11<=ipert<=2*natom+11 implemented!'
161 0 : ABI_BUG(msg)
162 : end if
163 72096 : pert_phon_elfd = .false.
164 72096 : if (ipert>natom+11.and.ipert<=2*natom+11) pert_phon_elfd = .true.
165 72096 : if (mpi_enreg%paral_spinor==1) then
166 0 : msg='Not compatible with parallelization over spinorial components!'
167 0 : ABI_BUG(msg)
168 : end if
169 72096 : if (gs_hamkq%nvloc>1) then
170 0 : msg='Not compatible with nvloc=4 (non-coll. magnetism)!'
171 0 : ABI_BUG(msg)
172 : end if
173 72096 : if((ipert==natom+11.or.pert_phon_elfd).and.gs_hamkq%usepaw==1.and.optnl>=1) then
174 25920 : if (gs_hamkq%nvloc>1) then
175 0 : msg='Not compatible with nvloc=4 (non-coll. magnetism)!'
176 0 : ABI_BUG(msg)
177 : end if
178 25920 : if (present(enl)) then
179 23040 : enl_ptr => enl
180 2880 : else if (associated(rf_hamkq%e1kbfr).and.associated(rf_hamkq%e1kbsc).and.optnl==2) then
181 2880 : ABI_CHECK(size(rf_hamkq%e1kbfr,4)==1,'BUG in getgh2c: qphase>1!')
182 2880 : ABI_CHECK(size(rf_hamkq%e1kbsc,4)==1,'BUG in getgh2c: qphase>1!')
183 17280 : ABI_MALLOC(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq))
184 221760 : enl_temp(:,:,:,:) = rf_hamkq%e1kbfr(:,:,:,:) + rf_hamkq%e1kbsc(:,:,:,:)
185 2880 : enl_ptr => enl_temp
186 0 : else if (associated(rf_hamkq%e1kbfr)) then
187 0 : ABI_CHECK(size(rf_hamkq%e1kbfr,4)==1,'BUG in getgh2c: qphase>1!')
188 0 : enl_ptr => rf_hamkq%e1kbfr
189 : else
190 0 : msg='For ipert=natom+11/pert_phon_elfd : e1kbfr and/or e1kbsc must be associated or enl optional input must be present.'
191 0 : ABI_BUG(msg)
192 : end if
193 25920 : if (usevnl==0) then
194 0 : msg='gvnl2 must be allocated for ipert=natom+11/pert_phon_elfd !'
195 0 : ABI_BUG(msg)
196 : end if
197 25920 : if(opt_gvnl2==0) then
198 0 : msg='opt_gvnl2=0 not compatible with ipert=natom+11/pert_phon_elfd !'
199 0 : ABI_BUG(msg)
200 : end if
201 : end if
202 :
203 : !Check sizes
204 72096 : my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
205 216288 : if (size(cwavef)<2*npw*my_nspinor) then
206 0 : msg='wrong size for cwavef!'
207 0 : ABI_BUG(msg)
208 : end if
209 216288 : if (size(gh2c)<2*npw1*my_nspinor) then
210 0 : msg='wrong size for gh2c!'
211 0 : ABI_BUG(msg)
212 : end if
213 72096 : if (usevnl/=0) then
214 77760 : if (size(gvnl2)<2*npw1*my_nspinor) then
215 0 : msg='wrong size for gvnl2!'
216 0 : ABI_BUG(msg)
217 : end if
218 : end if
219 72096 : if (sij_opt==1) then
220 14400 : if (size(gs2c)<2*npw1*my_nspinor) then
221 0 : msg='wrong size for gs2c!'
222 0 : ABI_BUG(msg)
223 : end if
224 : end if
225 :
226 : !PAW: specific treatment for usecprj input arg
227 : ! force it to zero if cwaveprj is not allocated
228 72096 : usecprj=gs_hamkq%usecprj ; ncpgr=0
229 72096 : if(gs_hamkq%usepaw==1) then
230 83520 : if (size(cwaveprj)==0) usecprj=0
231 960 : if (usecprj/=0) then
232 960 : ncpgr=cwaveprj(1,1)%ncpgr
233 2880 : if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then
234 0 : msg='wrong size for cwaveprj!'
235 0 : ABI_BUG(msg)
236 : end if
237 : end if
238 : else
239 44256 : if(usecprj==1)then
240 0 : msg='usecprj==1 not allowed for NC psps !'
241 0 : ABI_BUG(msg)
242 : end if
243 : end if
244 :
245 : ! tim_nonlop=8
246 : ! if (tim_getgh2c==1.and.ipert<=natom) tim_nonlop=7
247 : ! if (tim_getgh2c==2.and.ipert<=natom) tim_nonlop=5
248 : ! if (tim_getgh2c==1.and.ipert> natom) tim_nonlop=8
249 : ! if (tim_getgh2c==2.and.ipert> natom) tim_nonlop=5
250 : ! if (tim_getgh2c==3 ) tim_nonlop=0
251 :
252 72096 : idir1=alpha(idir);idir2=beta(idir)
253 :
254 72096 : compute_conjugate = .false.
255 72096 : if(present(conj)) compute_conjugate = conj
256 :
257 : !======================================================================
258 : !== Apply the 2nd-order local potential to the wavefunction
259 : !======================================================================
260 :
261 72096 : if (ipert/=natom+10.and.ipert/=natom+11.and.optlocal>0) then
262 0 : msg='local part not implemented'
263 0 : ABI_BUG(msg)
264 : else
265 : ! In the case of ddk operator, no local contribution (also because no self-consistency)
266 : !$OMP PARALLEL DO
267 8760796 : do ipw=1,npw1*my_nspinor
268 26138196 : gh2c(:,ipw)=zero
269 : end do
270 :
271 : end if
272 :
273 : !======================================================================
274 : !== Apply the 2st-order non-local potential to the wavefunction
275 : !======================================================================
276 :
277 72096 : has_vnl=(ipert==natom+10.or.ipert==natom+11.or.pert_phon_elfd)
278 :
279 : !Use of gvnl2 depends on usevnl
280 72096 : if (usevnl==1) then
281 25920 : gvnl2_ => gvnl2
282 : else
283 138528 : ABI_MALLOC(gvnl2_,(2,npw1*my_nspinor))
284 : end if
285 :
286 72096 : if (has_vnl.and.(optnl>0.or.sij_opt/=0)) then
287 :
288 72096 : idirc=3*(idir1-1)+idir2 !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(idir1,idir2)
289 :
290 : ! d^2[H_nl]/dk1dk2
291 : ! -------------------------------------------
292 72096 : if (ipert==natom+10) then
293 46176 : if (gs_hamkq%usepaw==1) then
294 1920 : if (usecprj==1) then
295 384 : cwaveprj_ptr => cwaveprj
296 : else
297 10752 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor))
298 1536 : call pawcprj_alloc(cwaveprj_tmp,0,gs_hamkq%dimcprj)
299 1536 : cwaveprj_ptr => cwaveprj_tmp
300 : end if
301 1920 : cpopt=-1+5*usecprj
302 1920 : choice=8; signs=2; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3
303 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/lambda/),mpi_enreg,1,nnlout,&
304 3840 : & paw_opt,signs,gs2c,tim_nonlop,cwavef,gvnl2_)
305 1920 : if (usecprj==0) then
306 1536 : call pawcprj_free(cwaveprj_tmp)
307 4608 : ABI_FREE(cwaveprj_tmp)
308 : end if
309 1920 : nullify(cwaveprj_ptr)
310 : else
311 44256 : choice=8; signs=2; cpopt=-1 ; paw_opt=0
312 : call nonlop(choice,cpopt,cwaveprj,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
313 44256 : & paw_opt,signs,svectout_dum,tim_nonlop,cwavef,gvnl2_)
314 : end if
315 :
316 : ! d^2[H_nl]/dk1dE2 : Non-zero only in PAW
317 : ! -------------------------------------------
318 25920 : else if (ipert==natom+11.and.gs_hamkq%usepaw==1) then
319 :
320 8640 : ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor))
321 :
322 2880 : if (usecprj==1) then
323 576 : cwaveprj_ptr => cwaveprj
324 : else
325 16128 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor))
326 2304 : call pawcprj_alloc(cwaveprj_tmp,2,gs_hamkq%dimcprj)
327 2304 : cwaveprj_ptr => cwaveprj_tmp
328 : end if
329 :
330 2880 : if (opt_gvnl2==1.and.optnl>=1) then
331 :
332 : ! Compute application of dS/dk1 to i*d[cwavef]/dk2
333 : ! sum_{i,j} s_ij d(|p_i><p_j|)/dk(idir1) | i*psi^(k(idir2)) >
334 2880 : cpopt=-1 ; choice=5 ; paw_opt=3 ; signs=2
335 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
336 2880 : & paw_opt,signs,nonlop_out,tim_nonlop,gvnl2_,vectout_dum)
337 :
338 : !$OMP PARALLEL DO
339 274320 : do ipw=1,npw1*my_nspinor
340 817200 : gvnl2_(:,ipw)=nonlop_out(:,ipw)
341 : end do
342 :
343 : ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
344 : ! sum_{i,j} chi_ij(idir2) d(|p_i><p_j|)/dk(idir1) | psi^(0) >
345 2880 : cpopt=4*usecprj ; choice=5 ; paw_opt=1 ; signs=2
346 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
347 2880 : & paw_opt,signs,svectout_dum,tim_nonlop,cwavef,nonlop_out,enl=enl_ptr)
348 :
349 : !$OMP PARALLEL DO
350 274320 : do ipw=1,npw1*my_nspinor
351 817200 : gvnl2_(:,ipw)=gvnl2_(:,ipw)+nonlop_out(:,ipw)
352 : end do
353 :
354 : else
355 :
356 : !$OMP PARALLEL DO
357 0 : do ipw=1,npw1*my_nspinor
358 0 : gvnl2_(:,ipw)=zero
359 : end do
360 :
361 : end if ! opt_gvnl2==1
362 :
363 : ! Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2>
364 : ! i * sum_{i,j} (d(|p_i><dp_j/dk(idir2)|)/dk(idir1) | psi^(0) >
365 2880 : cpopt=-1+5*usecprj ; choice=81 ; paw_opt=3 ; signs=2
366 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
367 2880 : & paw_opt,signs,nonlop_out,tim_nonlop,cwavef,vectout_dum)
368 :
369 2880 : if(compute_conjugate) then
370 : !$OMP PARALLEL DO
371 0 : do ipw=1,npw1*my_nspinor ! Note the multiplication by -i
372 0 : gvnl2_(1,ipw)=gvnl2_(1,ipw)+nonlop_out(2,ipw)
373 0 : gvnl2_(2,ipw)=gvnl2_(2,ipw)-nonlop_out(1,ipw)
374 : end do
375 : else
376 : !$OMP PARALLEL DO
377 274320 : do ipw=1,npw1*my_nspinor ! Note the multiplication by i
378 271440 : gvnl2_(1,ipw)=gvnl2_(1,ipw)-nonlop_out(2,ipw)
379 274320 : gvnl2_(2,ipw)=gvnl2_(2,ipw)+nonlop_out(1,ipw)
380 : end do
381 : end if
382 :
383 2880 : ABI_FREE(nonlop_out)
384 820080 : if (sij_opt==1) gs2c=zero
385 2880 : if (usecprj==0) then
386 2304 : call pawcprj_free(cwaveprj_tmp)
387 6912 : ABI_FREE(cwaveprj_tmp)
388 : end if
389 2880 : nullify(cwaveprj_ptr)
390 :
391 : ! d^2[H_nl]/dtau1dE2 : Non-zero only in PAW
392 : ! -------------------------------------------
393 23040 : else if (pert_phon_elfd.and.gs_hamkq%usepaw==1) then
394 :
395 23040 : iatm = ipert-(natom+11)
396 23040 : if (iatm<1.or.iatm>natom) then
397 0 : ABI_BUG(" iatm must be between 1 and natom")
398 : end if
399 :
400 69120 : ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor))
401 :
402 23040 : if (usecprj==1) then
403 0 : cwaveprj_ptr => cwaveprj
404 : else
405 161280 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor))
406 23040 : call pawcprj_alloc(cwaveprj_tmp,2,gs_hamkq%dimcprj)
407 23040 : cwaveprj_ptr => cwaveprj_tmp
408 : end if
409 :
410 23040 : if (opt_gvnl2==1) then
411 :
412 : ! Compute application of dS/dtau1 to i*d[cwavef]/dk2
413 : ! sum_{i,j} s_ij d(|p_i><p_j|)/dtau(idir1) | i*psi^(k(idir2)) >
414 23040 : cpopt=-1 ; choice=2 ; paw_opt=3 ; signs=2
415 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
416 23040 : & paw_opt,signs,nonlop_out,tim_nonlop,gvnl2_,vectout_dum,iatom_only=iatm)
417 :
418 : !$OMP PARALLEL DO
419 2194560 : do ipw=1,npw1*my_nspinor
420 6537600 : gvnl2_(:,ipw)=nonlop_out(:,ipw)
421 : end do
422 :
423 : ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
424 : ! sum_{i,j} chi_ij(idir2) d(|p_i><p_j|)/dtau(idir1) | psi^(0) >
425 23040 : cpopt=4*usecprj ; choice=2 ; paw_opt=1 ; signs=2
426 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
427 23040 : & paw_opt,signs,svectout_dum,tim_nonlop,cwavef,nonlop_out,enl=enl_ptr,iatom_only=iatm)
428 :
429 : !$OMP PARALLEL DO
430 2194560 : do ipw=1,npw1*my_nspinor
431 6537600 : gvnl2_(:,ipw)=gvnl2_(:,ipw)+nonlop_out(:,ipw)
432 : end do
433 :
434 : else
435 :
436 : !$OMP PARALLEL DO
437 0 : do ipw=1,npw1*my_nspinor
438 0 : gvnl2_(:,ipw)=zero
439 : end do
440 :
441 : end if ! opt_gvnl2==1
442 :
443 : ! Compute derivatives due to projectors |d^2[p_i]/dtau1dk2>,|d[p_i]/dtau1>,|d[p_i]/dk2>
444 : ! i * sum_{i,j} (d(|p_i><dp_j/dk(idir2)|)/dtau(idir1) | psi^(0) >
445 23040 : cpopt=-1+5*usecprj ; choice=54 ; paw_opt=3 ; signs=2
446 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
447 23040 : & paw_opt,signs,nonlop_out,tim_nonlop,cwavef,vectout_dum,iatom_only=iatm)
448 :
449 23040 : if(compute_conjugate) then
450 : !$OMP PARALLEL DO
451 0 : do ipw=1,npw1*my_nspinor ! Note the multiplication by -i
452 0 : gvnl2_(1,ipw)=gvnl2_(1,ipw)+nonlop_out(2,ipw)
453 0 : gvnl2_(2,ipw)=gvnl2_(2,ipw)-nonlop_out(1,ipw)
454 : end do
455 : else
456 : !$OMP PARALLEL DO
457 2194560 : do ipw=1,npw1*my_nspinor ! Note the multiplication by i
458 2171520 : gvnl2_(1,ipw)=gvnl2_(1,ipw)-nonlop_out(2,ipw)
459 2194560 : gvnl2_(2,ipw)=gvnl2_(2,ipw)+nonlop_out(1,ipw)
460 : end do
461 : end if
462 :
463 23040 : ABI_FREE(nonlop_out)
464 23040 : if (sij_opt==1) gs2c=zero
465 23040 : if (usecprj==0) then
466 23040 : call pawcprj_free(cwaveprj_tmp)
467 69120 : ABI_FREE(cwaveprj_tmp)
468 : end if
469 23040 : nullify(cwaveprj_ptr)
470 :
471 : end if
472 :
473 : !No non-local part
474 : !-------------------------------------------
475 : else
476 :
477 0 : if (optnl>=1) then
478 : !$OMP PARALLEL DO
479 0 : do ipw=1,npw1*my_nspinor
480 0 : gvnl2_(:,ipw)=zero
481 : end do
482 : end if
483 0 : if (sij_opt/=0) then
484 : !$OMP PARALLEL DO
485 0 : do ipw=1,npw1*my_nspinor
486 0 : gs2c(:,ipw)=zero
487 : end do
488 : end if
489 :
490 : end if
491 :
492 72096 : if (associated(enl_ptr)) then
493 37716 : nullify(enl_ptr)
494 : end if
495 72096 : if (allocated(enl_temp)) then
496 2880 : ABI_FREE(enl_temp)
497 : end if
498 :
499 : !======================================================================
500 : !== Apply the 2nd-order kinetic operator to the wavefunction
501 : !======================================================================
502 :
503 72096 : if (present(optkin)) then
504 7872 : has_kin=(optkin/=0.and.ipert==natom+10)
505 : else
506 64224 : has_kin=(ipert==natom+10)
507 : end if
508 :
509 : !k-point perturbation
510 : !-------------------------------------------
511 72096 : if (associated(gs_hamkq%kinpw_kp)) then
512 72096 : kinpw1 => gs_hamkq%kinpw_kp
513 0 : else if (optnl>=1.or.has_kin) then
514 0 : msg='need kinpw1 allocated!'
515 0 : ABI_BUG(msg)
516 : end if
517 72096 : if (associated(rf_hamkq%ddkinpw_k)) then
518 46176 : ddkinpw => rf_hamkq%ddkinpw_k
519 25920 : else if (has_kin) then
520 0 : msg='need ddkinpw allocated!'
521 0 : ABI_BUG(msg)
522 : end if
523 :
524 46176 : if (has_kin) then
525 91584 : do ispinor=1,my_nspinor
526 : !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwavef,ispinor,gvnl2_,ddkinpw,kinpw1,npw,my_nspinor)
527 6301132 : do ipw=1,npw
528 6209548 : ipws=ipw+npw*(ispinor-1)
529 6255340 : if(kinpw1(ipw)<huge(zero)*1.d-11)then
530 6209548 : gvnl2_(1,ipws)=gvnl2_(1,ipws)+ddkinpw(ipw)*cwavef(1,ipws)
531 6209548 : gvnl2_(2,ipws)=gvnl2_(2,ipws)+ddkinpw(ipw)*cwavef(2,ipws)
532 : else
533 0 : gvnl2_(1,ipws)=zero
534 0 : gvnl2_(2,ipws)=zero
535 : end if
536 : end do
537 : end do
538 : end if
539 :
540 : !======================================================================
541 : !== Sum contributions to get the application of H^(2) to the wf
542 : !======================================================================
543 : !Also filter the wavefunctions for large modified kinetic energy
544 :
545 : !Add non-local+kinetic to local part
546 72096 : if (optnl>=1.or.has_kin) then
547 144192 : do ispinor=1,my_nspinor
548 72096 : ipws=(ispinor-1)*npw1
549 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh2c,gvnl2_,kinpw1,ipws,npw1)
550 8832892 : do ipw=1+ipws,npw1+ipws
551 8760796 : if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then
552 8688700 : gh2c(1,ipw)=gh2c(1,ipw)+gvnl2_(1,ipw)
553 8688700 : gh2c(2,ipw)=gh2c(2,ipw)+gvnl2_(2,ipw)
554 : else
555 0 : gh2c(1,ipw)=zero
556 0 : gh2c(2,ipw)=zero
557 : end if
558 : end do
559 : end do
560 : end if
561 :
562 72096 : if (usevnl==1) then
563 25920 : nullify(gvnl2_)
564 : else
565 46176 : ABI_FREE(gvnl2_)
566 : end if
567 :
568 : ABI_NVTX_END_RANGE()
569 :
570 : !call timab(196+tim_getgh2c,2,tsec)
571 :
572 : DBG_EXIT("COLL")
573 :
574 216288 : end subroutine getgh2c
575 : !!***
576 :
577 : end module m_getgh2c
578 : !!***
|
