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