Line data Source code
1 : !!****m* ABINIT/m_getgh1c
2 : !! NAME
3 : !! m_getgh1c
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (XG, DRH, MT, SPr)
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_getgh1c
26 :
27 : use defs_basis
28 : use m_abicore
29 : use m_errors
30 : use m_dtset
31 : use m_xmpi
32 : use m_xomp
33 : use m_gputk
34 : use m_abi_linalg
35 :
36 : use, intrinsic :: iso_c_binding, only : c_ptr,c_loc,c_size_t
37 : use defs_abitypes, only : MPI_type
38 : use defs_datatypes, only : pseudopotential_type
39 : use m_time, only : timab
40 : use m_fstrings, only : sjoin, ltoa, itoa
41 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, &
42 : pawcprj_copy, pawcprj_lincom, pawcprj_axpby, pawcprj_mpi_sum
43 : use m_kg, only : kpgstr, mkkin, mkkpg, mkkin_metdqdq
44 : use m_mkffnl, only : mkffnl
45 : use m_pawfgr, only : pawfgr_type
46 : use m_fft, only : fftpac, fourwf
47 : use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type
48 : use m_cgtools, only : projbd
49 : use m_nonlop, only : nonlop
50 : use m_fourier_interpol, only : transgrid
51 :
52 : #if defined(HAVE_GPU_MARKERS)
53 : use m_nvtx_data
54 : #endif
55 :
56 : implicit none
57 :
58 : private
59 : !!***
60 :
61 : public :: getgh1c
62 : public :: rf_transgrid_and_pack
63 : public :: getgh1c_setup
64 : public :: getdc1
65 : public :: getgh1dqc
66 : public :: getgh1dqc_setup
67 : public :: getgh1ndc
68 : public :: getgh1c_mGGA
69 : !!***
70 :
71 : contains
72 : !!***
73 :
74 : !!****f* ABINIT/getgh1c
75 : !!
76 : !! NAME
77 : !! getgh1c
78 : !!
79 : !! FUNCTION
80 : !! Compute <G|H^(1)|C> (or <G|H^(1)-lambda.S^(1)|C>) for input vector |C> expressed in reciprocal space.
81 : !! (H^(1) is the 1st-order pertubed Hamiltonian, S^(1) is the 1st-order perturbed overlap operator).
82 : !! Result is put in array gh1c.
83 : !! If required, part of <G|K(1)+Vnonlocal^(1)|C> not depending on VHxc^(1) is also returned in gvnlx1c.
84 : !! If required, <G|S^(1)|C> is returned in gs1c (S=overlap - PAW only)
85 : !!
86 : !! INPUTS
87 : !! berryopt=option for Berry phase
88 : !! cwave(2,npw*nspinor)=input wavefunction, in reciprocal space
89 : !! cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives)
90 : !! if not allocated or size=0, they are locally computed (and not sorted)
91 : !! dkinpw(npw)=derivative of the (modified) kinetic energy for each plane wave at k (Hartree)
92 : !! grad_berry(2,npw1*nspinor*(berryopt/4))= the gradient of the Berry phase term
93 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
94 : !! idir=direction of the perturbation
95 : !! ipert=type of the perturbation
96 : !! lambda=real use to apply H^(1)-lambda.S^(1)
97 : !! mpi_enreg=information about MPI parallelization
98 : !! ndat=number of bands to compute in parallel
99 : !! npw=number of planewaves in basis sphere at given k.
100 : !! npw1=number of planewaves in basis sphere at k+q
101 : !! optlocal=0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
102 : !! 1: local part of H^(1) is computed in gh1c=<G|H^(1)|C>
103 : !! optnl=0: non-local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
104 : !! 1: non-local part of H^(1) depending on VHxc^(1) is not computed in gh1c=<G|H^(1)|C>
105 : !! 2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
106 : !! opt_gvnlx1=option controlling the use of gvnlx1 array:
107 : !! 0: used as an output
108 : !! 1: used as an input: (only for ipert=natom+2)
109 : !! NCPP: contains the ddk 1-st order WF
110 : !! PAW: contains frozen part of 1st-order hamiltonian
111 : !! 2: used as input/output: - used only for PAW and ipert=natom+2
112 : !! At input: contains the ddk 1-st order WF (times i)
113 : !! At output: contains frozen part of 1st-order hamiltonian
114 : !! rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,k+q
115 : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H^(1)|C> have to be computed
116 : !! (S=overlap) if 1, matrix elements <G|S^(1)|C> have to be computed in gs1c in addition to gh1c
117 : !! if -1, matrix elements <G|H^(1)-lambda.S^(1)|C> have to be computed in gh1c (gs1c not used)
118 : !! tim_getgh1c=timing code of the calling subroutine (can be set to 0 if not attributed)
119 : !! usevnl=1 if gvnlx1=(part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1)) has to be input/output
120 : !!
121 : !! OUTPUT
122 : !! gh1c(2,npw1*nspinor)= <G|H^(1)|C> or <G|H^(1)-lambda.S^(1)|C> on the k+q sphere
123 : !! (only kinetic+non-local parts if optlocal=0)
124 : !! if (usevnl==1)
125 : !! gvnlx1(2,npw1*nspinor*ndat)= part of <G|K^(1)+Vnl^(1)|C> not depending on VHxc^(1) (sij_opt/=-1)
126 : !! or part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1) (sij_opt==-1)
127 : !! if (sij_opt=1)
128 : !! gs1c(2,npw1*nspinor)=<G|S^(1)|C> (S=overlap) on the k+q sphere.
129 : !!
130 : !! SOURCE
131 :
132 8601471 : subroutine getgh1c(berryopt,cwave,cwaveprj,gh1c,grad_berry,gs1c,gs_hamkq,&
133 8601471 : gvnlx1,idir,ipert,lambda,mpi_enreg,ndat,optlocal,optnl,opt_gvnlx1,&
134 : rf_hamkq,sij_opt,tim_getgh1c,usevnl, &
135 : conj) ! optional
136 :
137 : !Arguments ------------------------------------
138 : !scalars
139 : logical,intent(in),optional :: conj
140 : integer,intent(in) :: berryopt,idir,ipert,optlocal,optnl,opt_gvnlx1,sij_opt,tim_getgh1c,usevnl,ndat
141 : type(MPI_type),intent(in) :: mpi_enreg
142 : type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
143 : type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq
144 : !arrays
145 : real(dp),intent(in) :: lambda(ndat)
146 : real(dp),intent(in) :: grad_berry(:,:)
147 : real(dp),intent(inout),target :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor*ndat)
148 : real(dp),intent(out),target :: gh1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor*ndat)
149 : real(dp),intent(out),target :: gs1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor*ndat)
150 : real(dp),intent(inout),target :: gvnlx1(2,gs_hamkq%npw_kp*gs_hamkq%nspinor*ndat)
151 : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
152 :
153 : !Local variables-------------------------------
154 : !scalars
155 : integer,parameter :: level=16
156 : integer :: choice,cplex1,cpopt,ipw,ipws,ispinor,istr,i1,i2,i3,idat, ndat__
157 : integer :: my_nspinor,natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs
158 : integer :: tim_fourwf,tim_nonlop,usecprj
159 : logical :: compute_conjugate,has_kin,has_mGGA1,has_nd1,usevnl2
160 : real(dp) :: weight !, cpu, wall, gflops
161 : real(dp),parameter :: threshold=huge(zero)*1.d-11
162 : !character(len=500) :: msg
163 : !arrays
164 17202942 : real(dp) :: enlout(ndat),tsec(2),svectout_dum(1,1),vectout_dum(1,1)
165 8601471 : real(dp),allocatable :: cwave_sp(:,:),cwavef1(:,:),cwavef2(:,:)
166 8601471 : real(dp),allocatable :: gh1c_sp(:,:),gh1c1(:,:),gh1c2(:,:),gh1c3(:,:),gh1c4(:,:)
167 8601471 : real(dp),allocatable :: gh1c_mGGA(:,:),gh1ndc(:,:),gvnl2(:,:)
168 8601471 : real(dp),target,allocatable :: nonlop_out(:,:),vlocal1_tmp(:,:,:)
169 : !#define _DEV_USE_WORK
170 : !#ifdef _DEV_USE_WORK
171 8601471 : real(dp),target,allocatable :: work(:,:,:,:)
172 : !#endif
173 8601471 : real(dp),contiguous, pointer :: gvnlx1_(:,:), dkinpw(:),kinpw1(:)
174 8601471 : type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:)
175 8601471 : type(pawcprj_type),pointer :: cwaveprj_ptr(:,:)
176 : #ifdef HAVE_OPENMP_OFFLOAD
177 : logical :: map_gh1c,map_gs1c,map_cwave,map_gvnlx1_
178 : #endif
179 : ! *********************************************************************
180 :
181 : DBG_ENTER("COLL")
182 :
183 : ! Keep track of total time spent in getgh1c
184 8601471 : call timab(196+tim_getgh1c,1,tsec)
185 : ABI_NVTX_START_RANGE(NVTX_GETGH1C)
186 :
187 : !======================================================================
188 : !== Initialisations and compatibility tests
189 : !======================================================================
190 :
191 8601471 : npw = gs_hamkq%npw_k
192 8601471 : npw1 = gs_hamkq%npw_kp
193 8601471 : natom = gs_hamkq%natom
194 :
195 : ! Compatibility tests
196 8601471 : if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then
197 671080 : if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))) .or. &
198 : (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc)))) then
199 0 : ABI_BUG('ekb derivatives must be allocated for ipert<=natom or natom+3/4 !')
200 : end if
201 : end if
202 8601471 : if(gs_hamkq%usepaw==1.and.(ipert==natom+2)) then
203 106864 : if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))) .or. &
204 : (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc)))) then
205 0 : ABI_BUG('ekb derivatives must be allocated for ipert=natom+2 !')
206 : end if
207 106864 : if (usevnl==0) then
208 0 : ABI_BUG('gvnlx1 must be allocated for ipert=natom+2 !')
209 : end if
210 : end if
211 8601471 : if(ipert==natom+2.and.opt_gvnlx1==0) then
212 0 : ABI_BUG('opt_gvnlx1=0 not compatible with ipert=natom+2 !')
213 : end if
214 8601471 : if (mpi_enreg%paral_spinor==1) then
215 0 : ABI_BUG('Not compatible with parallelization over spinorial components !')
216 : end if
217 :
218 : ! Check sizes
219 8601471 : my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
220 25804413 : ABI_CHECK_IGEQ(size(cwave), 2*npw*my_nspinor*ndat, 'wrong size for cwave!')
221 25804413 : ABI_CHECK_IGEQ(size(gh1c), 2*npw1*my_nspinor*ndat, 'wrong size for gh1c!')
222 8601471 : if (usevnl/=0) then
223 13997223 : ABI_CHECK_IGEQ(size(gvnlx1), 2*npw1*my_nspinor*ndat, 'wrong size for gvnlx1!')
224 : end if
225 8601471 : if (sij_opt==1) then
226 3019308 : ABI_CHECK_IGEQ(size(gs1c), 2*npw1*my_nspinor*ndat, 'wrong size for gs1c!')
227 : end if
228 8601471 : if (berryopt>=4) then
229 9312 : ABI_CHECK_IGEQ(size(grad_berry), 2*npw1*my_nspinor*ndat, 'wrong size for grad_berry!')
230 : end if
231 :
232 : ! PAW: specific treatment for usecprj input arg. force it to zero if cwaveprj is not allocated
233 8601471 : usecprj=gs_hamkq%usecprj ; ncpgr=0
234 8601471 : if(gs_hamkq%usepaw==1) then
235 3110328 : if (size(cwaveprj)==0) usecprj=0
236 1007976 : if (usecprj/=0) then
237 1007976 : ncpgr=cwaveprj(1,1)%ncpgr
238 3023928 : ABI_CHECK_IGEQ(size(cwaveprj), gs_hamkq%natom*my_nspinor*ndat, 'wrong size for cwaveprj!')
239 1007976 : if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then
240 663400 : if (ncpgr/=1)then
241 0 : ABI_BUG('Projected WFs (cprj) derivatives are not correctly stored !')
242 : end if
243 : end if
244 : end if
245 : else
246 7564695 : if(usecprj==1)then
247 0 : ABI_BUG('usecprj==1 not allowed for NC psps !')
248 : end if
249 : end if
250 :
251 8601471 : tim_nonlop=8
252 8601471 : if (tim_getgh1c==1.and.ipert<=natom) tim_nonlop=7
253 8601471 : if (tim_getgh1c==2.and.ipert<=natom) tim_nonlop=5
254 8601471 : if (tim_getgh1c==1.and.ipert> natom) tim_nonlop=8
255 8601471 : if (tim_getgh1c==2.and.ipert> natom) tim_nonlop=5
256 8601471 : if (tim_getgh1c==3 ) tim_nonlop=0
257 :
258 8601471 : compute_conjugate = .false.
259 8601471 : if(present(conj)) compute_conjugate = conj
260 :
261 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) then
262 : #ifdef HAVE_OPENMP_OFFLOAD
263 : map_gh1c = .not. ( xomp_target_is_present(c_loc(gh1c)))
264 : map_gs1c = .not. ( xomp_target_is_present(c_loc(gs1c))) .and. sij_opt==1
265 : map_cwave = .not. ( xomp_target_is_present(c_loc(cwave)))
266 : !print *, "getgh1c: map_gh1c=", map_gh1c, ", map_gs1c=", map_gs1c, ", map_cwave=", map_cwave
267 :
268 : !$OMP TARGET ENTER DATA MAP(alloc:gh1c) IF(map_gh1c)
269 : !$OMP TARGET ENTER DATA MAP(alloc:gs1c) IF(map_gs1c)
270 : !$OMP TARGET ENTER DATA MAP(to:cwave) IF(map_cwave)
271 : #endif
272 : end if
273 8601471 : has_mGGA1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vxctaulocal) )
274 :
275 : !======================================================================
276 : !== Apply the 1st-order local potential to the wavefunction
277 : !======================================================================
278 : ! Phonon perturbation or
279 : ! Electric field perturbation or
280 : ! Strain perturbation
281 : !-------------------------------------------
282 8601471 : if ((ipert<=natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)).and.ipert/=natom+1.and.optlocal>0) then !SPr deb
283 :
284 : ABI_NVTX_START_RANGE(NVTX_GETGH1_LOCPOT)
285 :
286 : ! Important: work (aka fofr) must be allocated with "ndat", but fourwf with option 2 never accesses it, so this allocation wastes memory.
287 : ! but it is neeeded to run on GPUs.
288 : ! Note that we use ndat__, since when nspinor 2 with nvloc 1, we can compute <g|vlocal1|u> for all ndat bands and the two spinor components
289 : ! with a single call to fourwf.
290 4259527 : ndat__ = ndat
291 4259527 : if (gs_hamkq%nvloc==1) ndat__ = ndat * gs_hamkq%nspinor
292 :
293 : !#ifdef _DEV_USE_WORK
294 21297635 : ABI_MALLOC(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6*ndat__))
295 : #ifdef HAVE_OPENMP_OFFLOAD
296 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(work,2*gs_hamkq%n4*gs_hamkq%n5*gs_hamkq%n6*ndat__)
297 : #endif
298 : !#endif
299 :
300 : !#else
301 : ! call gs_hamkq%alloc_fofr_work(ndat__)
302 : ! associate (work => gs_hamkq%fofr_work)
303 : !#endif
304 :
305 4259527 : if (gs_hamkq%nvloc==1) then
306 :
307 4170687 : weight=one; tim_fourwf=4
308 : call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gh1c,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
309 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat__,gs_hamkq%ngfft,&
310 4170687 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight, gpu_option=gs_hamkq%gpu_option)
311 :
312 4170687 : if(gs_hamkq%nspinor==2 .and. ndat__ == ndat)then
313 : ! Note: when ndat__ = ndat*nspinor (nvloc==1), fourwf above already handles both spinors
314 : ! This block is skipped.
315 0 : ABI_CHECK_IEQ(ndat, 1, "ndat > 1 with nspinor 2 and nspden 1 is buggy")
316 0 : ABI_MALLOC(cwave_sp,(2,npw))
317 0 : ABI_MALLOC(gh1c_sp,(2,npw1))
318 : !$OMP PARALLEL DO
319 0 : do ipw=1,npw
320 0 : cwave_sp(1,ipw)=cwave(1,ipw+npw)
321 0 : cwave_sp(2,ipw)=cwave(2,ipw+npw)
322 : end do
323 : call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave_sp,gh1c_sp,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
324 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
325 0 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight, gpu_option=gs_hamkq%gpu_option)
326 : !$OMP PARALLEL DO
327 0 : do ipw=1,npw1
328 0 : gh1c(1,ipw+npw1)=gh1c_sp(1,ipw)
329 0 : gh1c(2,ipw+npw1)=gh1c_sp(2,ipw)
330 : end do
331 0 : ABI_FREE(cwave_sp)
332 0 : ABI_FREE(gh1c_sp)
333 : end if ! nspinor == 2
334 :
335 : else
336 : ! Non-Collinear magnetism for nvloc=4
337 88840 : if (gs_hamkq%nspinor==2) then
338 88840 : weight=one ; tim_fourwf=4
339 266520 : ABI_MALLOC(gh1c1,(2,npw1*ndat))
340 177680 : ABI_MALLOC(gh1c2,(2,npw1*ndat))
341 177680 : ABI_MALLOC(gh1c3,(2,npw1*ndat))
342 177680 : ABI_MALLOC(gh1c4,(2,npw1*ndat))
343 205912456 : gh1c1(:,:)=zero; gh1c2(:,:)=zero; gh1c3(:,:)=zero ; gh1c4(:,:)=zero
344 444200 : ABI_MALLOC(vlocal1_tmp,(rf_hamkq%cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
345 : !SPr: notation/dimension corrected vlocal_tmp -> vlocal1_tmp
346 266520 : ABI_MALLOC(cwavef1,(2,npw*ndat))
347 177680 : ABI_MALLOC(cwavef2,(2,npw*ndat))
348 177680 : do idat=1,ndat
349 17329648 : do ipw=1,npw
350 51455904 : cwavef1(1:2,ipw+(idat-1)*npw)=cwave(1:2,ipw+(idat-1)*my_nspinor*npw)
351 51544744 : cwavef2(1:2,ipw+(idat-1)*npw)=cwave(1:2,ipw+(idat-1)*my_nspinor*npw+npw)
352 : end do
353 : end do
354 : ! gh1c1=v11*phi1
355 748753432 : vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,1)
356 : call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef1,gh1c1,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
357 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
358 88840 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight, gpu_option=gs_hamkq%gpu_option)
359 : ! gh1c2=v22*phi2
360 748753432 : vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,2)
361 : call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef2,gh1c2,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
362 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
363 88840 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight, gpu_option=gs_hamkq%gpu_option)
364 88840 : ABI_FREE(vlocal1_tmp)
365 88840 : cplex1=2
366 444200 : ABI_MALLOC(vlocal1_tmp,(cplex1*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
367 : ! gh1c3=(re(v12)-im(v12))*phi1 => v^21*phi1
368 88840 : if(rf_hamkq%cplex==1) then
369 816984 : do i3=1,gs_hamkq%n6
370 12716424 : do i2=1,gs_hamkq%n5
371 193744656 : do i1=1,gs_hamkq%n4
372 181062000 : vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(i1,i2,i3,3)
373 192961440 : vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(i1,i2,i3,4)
374 : end do
375 : end do
376 : end do
377 : else
378 : !SPr: modified definition of local potential components for cplex=2 (see dotprod_vn)
379 : !also, v21==v12* not always holds (e.g. magnetic field perturbation)
380 747808 : do i3=1,gs_hamkq%n6
381 14195488 : do i2=1,gs_hamkq%n5
382 284530176 : do i1=1,gs_hamkq%n4
383 270389760 : vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,4)
384 283837440 : vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(2*i1-1,i2,i3,4)
385 : end do
386 : end do
387 : end do
388 : end if
389 : call fourwf(cplex1,vlocal1_tmp,cwavef1,gh1c3,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
390 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
391 88840 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,gpu_option=gs_hamkq%gpu_option)
392 : ! gh1c4=(re(v12)+im(v12))*phi2 => v^12*phi2
393 88840 : if(rf_hamkq%cplex==1) then
394 816984 : do i3=1,gs_hamkq%n6
395 12716424 : do i2=1,gs_hamkq%n5
396 193744656 : do i1=1,gs_hamkq%n4
397 192961440 : vlocal1_tmp(2*i1,i2,i3)=-vlocal1_tmp(2*i1,i2,i3)
398 : end do
399 : end do
400 : end do
401 : else
402 : !for cplex=2 and time-reversal breaking perturbations,v21/=v12*
403 747808 : do i3=1,gs_hamkq%n6
404 14195488 : do i2=1,gs_hamkq%n5
405 284530176 : do i1=1,gs_hamkq%n4
406 270389760 : vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1-1,i2,i3,3)
407 283837440 : vlocal1_tmp(2*i1 ,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,3)
408 : end do
409 : end do
410 : end do
411 : end if
412 : call fourwf(cplex1,vlocal1_tmp,cwavef2,gh1c4,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
413 : gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
414 88840 : npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight, gpu_option=gs_hamkq%gpu_option)
415 88840 : ABI_FREE(vlocal1_tmp)
416 : ! Build gh1c from pieces
417 : ! gh1c_1 = (v11, v12) (psi1) matrix vector product
418 : ! gh1c_2 = (v12*,v22) (psi2)
419 177680 : do idat=1,ndat
420 17329648 : do ipw=1,npw1
421 51455904 : gh1c(1:2,ipw+(idat-1)*my_nspinor*npw1) =gh1c1(1:2,ipw+(idat-1)*my_nspinor*npw1)+gh1c4(1:2,ipw+(idat-1)*my_nspinor*npw1)
422 51544744 : gh1c(1:2,ipw+(idat-1)*my_nspinor*npw1+npw1)=gh1c3(1:2,ipw+(idat-1)*my_nspinor*npw1)+gh1c2(1:2,ipw+(idat-1)*my_nspinor*npw1)
423 : end do
424 : end do
425 88840 : ABI_FREE(gh1c1)
426 88840 : ABI_FREE(gh1c2)
427 88840 : ABI_FREE(gh1c3)
428 88840 : ABI_FREE(gh1c4)
429 88840 : ABI_FREE(cwavef1)
430 88840 : ABI_FREE(cwavef2)
431 : else
432 0 : ABI_BUG('nspinor/=1 for Non-collinear calculations!')
433 : end if
434 : end if ! nvloc
435 :
436 : !#ifndef _DEV_USE_WORK
437 : ! end associate
438 : !#endif
439 :
440 : !#ifdef _DEV_USE_WORK
441 : #ifdef HAVE_OPENMP_OFFLOAD
442 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(work,2*gs_hamkq%n4*gs_hamkq%n5*gs_hamkq%n6*ndat__)
443 : #endif
444 4259527 : ABI_FREE(work)
445 : !#endif
446 : ABI_NVTX_END_RANGE()
447 :
448 : ! k-point perturbation (or no local part, i.e. optlocal=0)
449 : ! -------------------------------------------
450 4341944 : else if (ipert==natom+1.or.optlocal==0) then
451 :
452 : ! In the case of ddk operator, no local contribution (also because no self-consistency)
453 4341944 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
454 : !$OMP PARALLEL DO
455 980236178 : do ipw=1,npw1*my_nspinor*ndat
456 2932024646 : gh1c(:,ipw)=zero
457 : end do
458 : else
459 : #ifdef HAVE_OPENMP_OFFLOAD
460 : call gpu_set_to_zero(gh1c,int(2,c_size_t)*npw1*my_nspinor*ndat)
461 : #endif
462 : end if
463 :
464 : end if
465 :
466 : !======================================================================
467 : !== Apply the 1st-order non-local potential to the wavefunction
468 : !======================================================================
469 :
470 : !Use of gvnlx1 depends on usevnl
471 8601471 : if (usevnl==1) then
472 4665741 : gvnlx1_ => gvnlx1
473 : else
474 11807190 : ABI_MALLOC(gvnlx1_,(2,npw1*my_nspinor*ndat))
475 : end if
476 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) then
477 : #ifdef HAVE_OPENMP_OFFLOAD
478 : map_gvnlx1_ = .not. ( xomp_target_is_present(c_loc(gvnlx1_)))
479 : if(map_gvnlx1_) call ompgpu_enter_map_alloc(gvnlx1_,2*npw1*my_nspinor*ndat)
480 : !$OMP TARGET UPDATE TO(gvnlx1_) IF(usevnl==1 .and. gs_hamkq%gpu_option==ABI_GPU_OPENMP)
481 : #endif
482 : end if
483 :
484 : !Phonon perturbation
485 : !-------------------------------------------
486 8601471 : if (ipert<=natom.and.(optnl>0.or.sij_opt/=0)) then
487 :
488 : ABI_NVTX_START_RANGE(NVTX_GETGH1_NLOCPOT)
489 :
490 : ! PAW:
491 6233393 : if (gs_hamkq%usepaw==1) then
492 :
493 608116 : if (usecprj==1) then
494 600436 : cwaveprj_ptr => cwaveprj
495 : else
496 53760 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor*ndat))
497 7680 : call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
498 7680 : cwaveprj_ptr => cwaveprj_tmp
499 : end if
500 :
501 : ! 1- Compute derivatives due to projectors |p_i>^(1)
502 : ! Only displaced atom contributes
503 608116 : cpopt=-1+5*usecprj ; choice=2 ; signs=2
504 608116 : paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3
505 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
506 608116 : paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_,iatom_only=ipert)
507 :
508 : ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1))
509 : ! All atoms contribute
510 608116 : if (optnl>=1) then
511 1824348 : ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor*ndat))
512 : #ifdef HAVE_OPENMP_OFFLOAD
513 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(nonlop_out,2*npw1*my_nspinor*ndat)
514 : #endif
515 608116 : cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1
516 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
517 608116 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
518 608116 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
519 : !$OMP PARALLEL DO
520 77194050 : do ipw=1,npw1*my_nspinor*ndat
521 383537786 : gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw)
522 : end do
523 : else
524 0 : call abi_xaxpy(2*npw1*my_nspinor*ndat, cone, nonlop_out, 1, gvnlx1_, 1, x_cplx=1, gpu_option=gs_hamkq%gpu_option)
525 : end if
526 : #ifdef HAVE_OPENMP_OFFLOAD
527 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(nonlop_out,2*npw1*my_nspinor*ndat)
528 : #endif
529 608116 : ABI_FREE(nonlop_out)
530 : end if
531 :
532 : ! 3- Compute derivatives due to self-consistent part of D_ij^(1) (depending on VHxc^(1))
533 : ! All atoms contribute
534 608116 : if (optnl==2) then
535 1095480 : ABI_MALLOC(gvnl2,(2,npw1*my_nspinor*ndat))
536 : #ifdef HAVE_OPENMP_OFFLOAD
537 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(gvnl2,2*npw1*my_nspinor*ndat)
538 : #endif
539 365160 : cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1
540 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
541 365160 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
542 : end if
543 :
544 608116 : if (usecprj==0) then
545 7680 : call pawcprj_free(cwaveprj_tmp)
546 23040 : ABI_FREE(cwaveprj_tmp)
547 : end if
548 608116 : nullify(cwaveprj_ptr)
549 :
550 : else
551 : ! Norm-conserving psps:
552 : ! Compute only derivatives due to projectors |p_i>^(1)
553 5625277 : cpopt=-1 ; choice=2 ; signs=2 ; paw_opt=0
554 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
555 5625277 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_,iatom_only=ipert)
556 5625277 : if (sij_opt==1) then
557 0 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
558 : !$OMP PARALLEL DO
559 0 : do ipw=1,npw1*my_nspinor*ndat
560 0 : gs1c(:,ipw)=zero
561 : end do
562 : else
563 0 : call gpu_set_to_zero(gs1c,int(2,c_size_t)*npw1*my_nspinor*ndat)
564 : end if
565 : end if
566 : end if
567 : ABI_NVTX_END_RANGE()
568 :
569 : ! k-point perturbation
570 : ! -------------------------------------------
571 2368078 : else if (ipert==natom+1.and.(optnl>0.or.sij_opt/=0)) then
572 :
573 : ABI_NVTX_START_RANGE(NVTX_GETGH1_DDK)
574 :
575 1437546 : tim_nonlop=8 ; signs=2 ; choice=5
576 1437546 : if (gs_hamkq%usepaw==1) then
577 258832 : if (usecprj==1) then
578 255376 : cwaveprj_ptr => cwaveprj
579 : else
580 24192 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor*ndat))
581 3456 : call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
582 3456 : cwaveprj_ptr => cwaveprj_tmp
583 : end if
584 258832 : cpopt=-1+5*usecprj; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3
585 : ! JLJ: BUG (wrong result) of H^(1) if stored cprj are used in PAW DDKs with nspinor==2 (==1 works fine).
586 : ! To be debugged, if someone has time...
587 258832 : if(gs_hamkq%nspinor==2) cpopt=-1
588 258832 : if(associated(gs_hamkq%vectornd)) cpopt=-1
589 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
590 258832 : paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_)
591 258832 : if (usecprj==0) then
592 3456 : call pawcprj_free(cwaveprj_tmp)
593 10368 : ABI_FREE(cwaveprj_tmp)
594 : end if
595 258832 : nullify(cwaveprj_ptr)
596 : else
597 1178714 : cpopt=-1 ; paw_opt=0
598 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
599 1178714 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_)
600 : end if
601 :
602 : ABI_NVTX_END_RANGE()
603 :
604 : !DEBUG
605 : ! gvnlx1_=zero
606 : !ENDDEBUG
607 :
608 :
609 : ! Electric field perturbation without Berry phase
610 : ! -------------------------------------------
611 : else if (ipert==natom+2 .and. &
612 : (berryopt/=4 .and. berryopt/=6 .and. berryopt/=7 .and. &
613 1441287 : berryopt/=14 .and. berryopt/=16 .and. berryopt/=17) .and.(optnl>0.or.sij_opt/=0))then
614 : ! gvnlx1 was already initialized in the calling routine, by reading a ddk file. It contains |i du^(0)/dk_band>
615 :
616 : ABI_NVTX_START_RANGE(NVTX_GETGH1_DDE)
617 :
618 510755 : if (gs_hamkq%usepaw==1) then
619 106864 : if (usecprj==1) then
620 89200 : cwaveprj_ptr => cwaveprj
621 : else
622 123648 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor*ndat))
623 17664 : call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
624 17664 : cwaveprj_ptr => cwaveprj_tmp
625 : end if
626 106864 : if (opt_gvnlx1==2.and.optnl>=1) then
627 :
628 : ! PAW: Compute application of S^(0) to ddk WF
629 106864 : cpopt=-1 ; choice=1 ; paw_opt=3 ; signs=2
630 320592 : ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor*ndat))
631 : #ifdef HAVE_OPENMP_OFFLOAD
632 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(nonlop_out,2*npw1*my_nspinor*ndat)
633 : #endif
634 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,lambda,mpi_enreg,ndat,nnlout,&
635 106864 : paw_opt,signs,nonlop_out,tim_nonlop,gvnlx1_,vectout_dum)
636 106864 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
637 : !$OMP PARALLEL DO
638 15681444 : do ipw=1,npw1*my_nspinor*ndat
639 77979764 : gvnlx1_(:,ipw)=nonlop_out(:,ipw)
640 : end do
641 : else
642 : #ifdef HAVE_OPENMP_OFFLOAD
643 : call gpu_copy(gvnlx1_, nonlop_out, int(2,c_size_t)*npw1*my_nspinor*ndat)
644 : #endif
645 : end if
646 :
647 : ! PAW: Compute part of H^(1) due to derivative of S
648 106864 : cpopt=4*usecprj ; choice=51 ; paw_opt=3 ; signs=2
649 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
650 106864 : paw_opt,signs,nonlop_out,tim_nonlop,cwave,vectout_dum)
651 :
652 106864 : if (gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
653 106864 : if(compute_conjugate) then
654 : !$OMP PARALLEL DO
655 0 : do ipw=1,npw1*my_nspinor*ndat ! Note the multiplication by -i
656 0 : gvnlx1_(1,ipw)=gvnlx1_(1,ipw)+nonlop_out(2,ipw)
657 0 : gvnlx1_(2,ipw)=gvnlx1_(2,ipw)-nonlop_out(1,ipw)
658 : end do
659 : else
660 : !$OMP PARALLEL DO
661 15681444 : do ipw=1,npw1*my_nspinor*ndat ! Note the multiplication by i
662 15574580 : gvnlx1_(1,ipw)=gvnlx1_(1,ipw)-nonlop_out(2,ipw)
663 15681444 : gvnlx1_(2,ipw)=gvnlx1_(2,ipw)+nonlop_out(1,ipw)
664 : end do
665 : end if
666 : else
667 : #ifdef HAVE_OPENMP_OFFLOAD
668 : if(compute_conjugate) then
669 : #ifdef HAVE_OPENMP_OFFLOAD
670 : !$OMP TARGET PARALLEL DO MAP(to:gvnlx1_,nonlop_out)
671 : do ipw=1,npw1*my_nspinor*ndat ! Note the multiplication by -i
672 : gvnlx1_(1,ipw)=gvnlx1_(1,ipw)+nonlop_out(2,ipw)
673 : gvnlx1_(2,ipw)=gvnlx1_(2,ipw)-nonlop_out(1,ipw)
674 : end do
675 : #endif
676 : else
677 : #ifdef HAVE_OPENMP_OFFLOAD
678 : !$OMP TARGET PARALLEL DO MAP(to:gvnlx1_,nonlop_out)
679 : do ipw=1,npw1*my_nspinor*ndat ! Note the multiplication by i
680 : gvnlx1_(1,ipw)=gvnlx1_(1,ipw)-nonlop_out(2,ipw)
681 : gvnlx1_(2,ipw)=gvnlx1_(2,ipw)+nonlop_out(1,ipw)
682 : end do
683 : #endif
684 : end if
685 : #endif
686 : end if
687 :
688 : ! PAW: Compute part of H^(1) due to derivative of electric field part of Dij
689 106864 : cpopt=2 ; choice=1 ; paw_opt=1 ; signs=2
690 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,lambda,mpi_enreg,ndat,nnlout,&
691 106864 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
692 :
693 106864 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
694 : !$OMP PARALLEL DO
695 15681444 : do ipw=1,npw1*my_nspinor*ndat
696 77979764 : gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw)
697 : end do
698 : else
699 : call abi_xaxpy(2*npw1*my_nspinor*ndat, cone, &
700 0 : & nonlop_out, 1, gvnlx1_, 1, x_cplx=1, gpu_option=gs_hamkq%gpu_option)
701 : end if
702 : #ifdef HAVE_OPENMP_OFFLOAD
703 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(nonlop_out,2*npw1*my_nspinor*ndat)
704 : #endif
705 106864 : ABI_FREE(nonlop_out)
706 :
707 : end if ! opt_gvnlx1==2
708 :
709 : ! PAW: Compute derivatives due to part of D_ij^(1) depending on VHxc^(1)
710 106864 : if (optnl>=2) then
711 229728 : ABI_MALLOC(gvnl2,(2,npw1*my_nspinor*ndat))
712 : #ifdef HAVE_OPENMP_OFFLOAD
713 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(gvnl2,2*npw1*my_nspinor*ndat)
714 : #endif
715 76576 : cpopt=-1+3*usecprj;if (opt_gvnlx1==2) cpopt=2
716 76576 : choice=1 ; paw_opt=1 ; signs=2
717 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,lambda,mpi_enreg,ndat,nnlout,&
718 76576 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
719 : end if
720 :
721 106864 : if (sij_opt==1) then
722 76576 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
723 : !$OMP PARALLEL DO
724 9716880 : do ipw=1,npw1*my_nspinor*ndat
725 28997488 : gs1c(:,ipw)=zero
726 : end do
727 : else
728 0 : call gpu_set_to_zero(gs1c,int(2,c_size_t)*npw1*my_nspinor*ndat)
729 : end if
730 : end if
731 106864 : if (usecprj==0) then
732 17664 : call pawcprj_free(cwaveprj_tmp)
733 52992 : ABI_FREE(cwaveprj_tmp)
734 : end if
735 106864 : nullify(cwaveprj_ptr)
736 : end if ! PAW
737 :
738 : ABI_NVTX_END_RANGE()
739 :
740 : ! Electric field perturbation with Berry phase
741 : ! -------------------------------------------
742 : else if (ipert==natom+2 .and. &
743 : (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
744 419777 : berryopt==14 .or. berryopt==16 .or. berryopt==17 ) .and.(optnl>0.or.sij_opt/=0))then
745 :
746 1152 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) then
747 0 : ABI_BUG("Not implemented for OpenMP GPU (gs_hamkq%gpu_option==2)")
748 : end if
749 1152 : if (optnl>=1) then
750 48960 : do ipw=1,npw1*my_nspinor*ndat
751 47808 : gvnlx1_(1,ipw)=-grad_berry(2,ipw)
752 48960 : gvnlx1_(2,ipw)= grad_berry(1,ipw)
753 : end do
754 : end if
755 1152 : if (sij_opt==1) then
756 0 : ABI_BUG("toto")
757 : !$OMP PARALLEL DO
758 0 : do ipw=1,npw1*my_nspinor*ndat
759 0 : gs1c(:,ipw)=zero
760 : end do
761 : end if
762 :
763 : ! Strain perturbation
764 : ! -------------------------------------------
765 418625 : else if ((ipert==natom+3.or.ipert==natom+4).and.(optnl>0.or.sij_opt/=0)) then
766 :
767 : ABI_NVTX_START_RANGE(NVTX_GETGH1_STRAIN)
768 :
769 371953 : istr=idir;if(ipert==natom+4) istr=istr+3
770 :
771 : ! PAW:
772 371953 : if (gs_hamkq%usepaw==1) then
773 :
774 62964 : if (usecprj==1) then
775 62964 : cwaveprj_ptr => cwaveprj
776 : else
777 0 : ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor*ndat))
778 0 : call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
779 0 : cwaveprj_ptr => cwaveprj_tmp
780 : end if
781 :
782 : ! 1- Compute derivatives due to projectors |p_i>^(1) All atoms contribute
783 62964 : cpopt=-1+5*usecprj ; choice=3 ; signs=2
784 62964 : paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3
785 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,lambda,mpi_enreg,ndat,nnlout,&
786 62964 : paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_)
787 :
788 : ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1)). All atoms contribute
789 62964 : if (optnl>=1) then
790 188892 : ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor*ndat))
791 : #ifdef HAVE_OPENMP_OFFLOAD
792 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(nonlop_out,2*npw1*my_nspinor*ndat)
793 : #endif
794 62964 : cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1
795 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,lambda,mpi_enreg,ndat,nnlout,&
796 62964 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
797 62964 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
798 : !$OMP PARALLEL DO
799 19198804 : do ipw=1,npw1*my_nspinor*ndat
800 95742164 : gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw)
801 : end do
802 : else
803 0 : call abi_xaxpy(2*npw1*my_nspinor*ndat, cone, nonlop_out, 1, gvnlx1_, 1, x_cplx=1, gpu_option=gs_hamkq%gpu_option)
804 : end if
805 : #ifdef HAVE_OPENMP_OFFLOAD
806 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(nonlop_out,2*npw1*my_nspinor*ndat)
807 : #endif
808 62964 : ABI_FREE(nonlop_out)
809 : end if
810 :
811 : ! 3- Compute derivatives due to part of D_ij^(1) depending on VHxc^(1). All atoms contribute
812 62964 : if (optnl>=2) then
813 73512 : ABI_MALLOC(gvnl2,(2,npw1*my_nspinor*ndat))
814 : #ifdef HAVE_OPENMP_OFFLOAD
815 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(gvnl2,2*npw1*my_nspinor*ndat)
816 : #endif
817 24504 : cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1
818 : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,lambda,mpi_enreg,ndat,nnlout,&
819 24504 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
820 : end if
821 :
822 62964 : if (usecprj==0) then
823 0 : call pawcprj_free(cwaveprj_tmp)
824 0 : ABI_FREE(cwaveprj_tmp)
825 : end if
826 62964 : nullify(cwaveprj_ptr)
827 :
828 : else
829 : ! Norm-conserving psps:
830 : ! Compute only derivatives due to projectors |p_i>^(1)
831 308989 : choice=3 ; cpopt=-1 ; signs=2 ; paw_opt=0
832 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,ndat,nnlout,&
833 926967 : paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_)
834 308989 : if (sij_opt==1) then
835 0 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
836 : !$OMP PARALLEL DO
837 0 : do ipw=1,npw1*my_nspinor*ndat
838 0 : gs1c(:,ipw)=zero
839 : end do
840 : else
841 0 : call gpu_set_to_zero(gs1c,int(2,c_size_t)*npw1*my_nspinor*ndat)
842 : end if
843 : end if
844 : end if
845 :
846 : ABI_NVTX_END_RANGE()
847 :
848 : ! No non-local part
849 : ! -------------------------------------------
850 46672 : else if (usevnl>0.or.(sij_opt/=0)) then
851 :
852 46488 : if (optnl>=1) then
853 46488 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
854 : !$OMP PARALLEL DO
855 22046312 : do ipw=1,npw1*my_nspinor*ndat
856 66045960 : gvnlx1_(:,ipw)=zero
857 : end do
858 : else
859 0 : call gpu_set_to_zero(gvnlx1_,int(2,c_size_t)*npw1*my_nspinor*ndat)
860 : end if
861 : end if
862 46488 : if (sij_opt/=0) then
863 0 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
864 : !$OMP PARALLEL DO
865 0 : do ipw=1,npw1*my_nspinor*ndat
866 0 : gs1c(:,ipw)=zero
867 : end do
868 : else
869 0 : call gpu_set_to_zero(gs1c,int(2,c_size_t)*npw1*my_nspinor*ndat)
870 : end if
871 : end if
872 :
873 : end if
874 :
875 : !======================================================================
876 : !== Apply the 1st-order kinetic operator to the wavefunction
877 : !== (add it to nl contribution)
878 : !======================================================================
879 :
880 : !Phonon perturbation or Electric field perturbation
881 : !-------------------------------------------
882 : !No kinetic contribution
883 :
884 : !k-point perturbation or Strain perturbation
885 : !-------------------------------------------
886 :
887 8601471 : usevnl2=allocated(gvnl2)
888 8601471 : has_kin=(ipert==natom+1.or.ipert==natom+3.or.ipert==natom+4)
889 8601471 : if (associated(gs_hamkq%kinpw_kp)) then
890 8601471 : kinpw1 => gs_hamkq%kinpw_kp
891 0 : else if (optnl>=1.or.usevnl2.or.has_kin) then
892 0 : ABI_BUG('need kinpw1 allocated!')
893 : end if
894 8601471 : if (associated(rf_hamkq%dkinpw_k)) then
895 5622552 : dkinpw => rf_hamkq%dkinpw_k
896 2978919 : else if (has_kin) then
897 0 : ABI_BUG('need dkinpw allocated!')
898 : end if
899 :
900 5622552 : if (has_kin) then
901 : ! Remember that npw=npw1 for ddk perturbation
902 : #ifdef HAVE_OPENMP_OFFLOAD
903 : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(ipws) MAP(to:cwave,gvnlx1_,dkinpw,kinpw1) &
904 : !$OMP& IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
905 : #endif
906 3677174 : do ispinor=1,my_nspinor*ndat
907 : !$OMP PARALLEL DO PRIVATE(ipw,ipws)
908 271039572 : do ipw=1,npw
909 267362398 : ipws=ipw+npw*(ispinor-1)
910 269230073 : if(kinpw1(ipw)<threshold)then
911 267362398 : gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+dkinpw(ipw)*cwave(1,ipws)
912 267362398 : gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+dkinpw(ipw)*cwave(2,ipws)
913 : else
914 0 : gvnlx1_(1,ipws)=zero
915 0 : gvnlx1_(2,ipws)=zero
916 : end if
917 : end do
918 : end do
919 : end if
920 :
921 : !======================================================================
922 : !== Apply the 1st-order nuclear dipole operator to the wavefunction
923 : !== Only coded for DDK
924 : !== (add it to nl contribution)
925 : !======================================================================
926 :
927 8601471 : has_nd1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vectornd) )
928 :
929 : if (has_nd1) then
930 378864 : ABI_MALLOC(gh1ndc,(2,npw*my_nspinor*ndat))
931 : #ifdef HAVE_OPENMP_OFFLOAD
932 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_enter_map_alloc(gh1ndc,2*npw*my_nspinor*ndat)
933 : #endif
934 : call getgh1ndc(cwave,gh1ndc,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,&
935 : & gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,npw,gs_hamkq%nvloc,&
936 : & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,my_nspinor,rf_hamkq%vectornd,&
937 126288 : & gs_hamkq%vlocal,rf_hamkq%zora,gs_hamkq%gpu_option)
938 : #ifdef HAVE_OPENMP_OFFLOAD
939 : !$OMP TARGET TEAMS DISTRIBUTE &
940 : !$OMP& MAP(to:gvnlx1_,gh1ndc) PRIVATE(ispinor) &
941 : !$OMP& IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
942 : #endif
943 285600 : do ispinor=1,my_nspinor*ndat
944 : !$OMP PARALLEL DO PRIVATE(ipw,ipws)
945 26845944 : do ipw=1,npw
946 26560344 : ipws=ipw+npw*(ispinor-1)
947 26560344 : gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+gh1ndc(1,ipws)
948 26719656 : gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+gh1ndc(2,ipws)
949 : end do
950 : end do
951 : #ifdef HAVE_OPENMP_OFFLOAD
952 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(gh1ndc,2*npw*my_nspinor*ndat)
953 : #endif
954 126288 : ABI_FREE(gh1ndc)
955 : end if
956 :
957 : !======================================================================
958 : !== Apply the 1st-order mGGA operator to the wavefunction
959 : !== Only coded for DDK
960 : !== (add it to nl contribution)
961 : !======================================================================
962 :
963 8601471 : has_mGGA1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vxctaulocal) )
964 :
965 : if (has_mGGA1) then
966 34560 : ABI_MALLOC(gh1c_mGGA,(2,npw*my_nspinor*ndat))
967 : #ifdef HAVE_OPENMP_OFFLOAD
968 : !$OMP TARGET UPDATE FROM(cwave) IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
969 : #endif
970 : call getgh1c_mGGA(cwave,gs_hamkq%gbound_k,gh1c_mGGA,gs_hamkq%gmet,gs_hamkq%gprimd,idir,gs_hamkq%istwf_k,&
971 : gs_hamkq%kg_k,gs_hamkq%kpt_k,gs_hamkq%mgfft,mpi_enreg,my_nspinor,gs_hamkq%n4,gs_hamkq%n5,&
972 : gs_hamkq%n6,ndat,gs_hamkq%ngfft,npw,gs_hamkq%nvloc,rf_hamkq%vxctaulocal,&
973 11520 : gpu_option=gs_hamkq%gpu_option)
974 : #ifdef HAVE_OPENMP_OFFLOAD
975 : !$OMP TARGET TEAMS DISTRIBUTE &
976 : !$OMP& MAP(to:gvnlx1_,gh1c_mGGA) PRIVATE(ispinor) &
977 : !$OMP& IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
978 : #endif
979 30720 : do ispinor=1,my_nspinor*ndat
980 : !$OMP PARALLEL DO PRIVATE(ipw,ipws)
981 2838720 : do ipw=1,npw
982 2808000 : ipws=ipw+npw*(ispinor-1)
983 2808000 : gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+gh1c_mGGA(1,ipws)
984 2827200 : gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+gh1c_mGGA(2,ipws)
985 : end do
986 : end do
987 11520 : ABI_FREE(gh1c_mGGA)
988 : end if
989 :
990 : !======================================================================
991 : !== Sum contributions to get the application of H^(1) to the wf
992 : !======================================================================
993 : !Also filter the wavefunctions for large modified kinetic energy
994 :
995 : !Add non-local+kinetic to local part
996 8601471 : if (optnl>=1.or.has_kin) then
997 : #ifdef HAVE_OPENMP_OFFLOAD
998 : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(ipws) MAP(to:gh1c,gvnlx1_,kinpw1) &
999 : !$OMP& IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1000 : #endif
1001 17812980 : do ispinor=1,my_nspinor*ndat
1002 9211693 : ipws=(ispinor-1)*npw1
1003 : !$OMP PARALLEL DO PRIVATE(ipw)
1004 2535750305 : do ipw=1+ipws,npw1+ipws
1005 2527148834 : if(kinpw1(ipw-ipws)<threshold)then
1006 2517937141 : gh1c(1,ipw)=gh1c(1,ipw)+gvnlx1_(1,ipw)
1007 2517937141 : gh1c(2,ipw)=gh1c(2,ipw)+gvnlx1_(2,ipw)
1008 : else
1009 0 : gh1c(1,ipw)=zero
1010 0 : gh1c(2,ipw)=zero
1011 : end if
1012 : end do
1013 : end do
1014 : end if
1015 :
1016 : !PAW: add non-local part due to first order change of VHxc
1017 8601471 : if (usevnl2) then
1018 : #ifdef HAVE_OPENMP_OFFLOAD
1019 : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(ipws) MAP(to:gh1c,gvnl2,kinpw1) &
1020 : !$OMP& IF(gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1021 : #endif
1022 932480 : do ispinor=1,my_nspinor*ndat
1023 466240 : ipws=(ispinor-1)*npw1
1024 : !$OMP PARALLEL DO PRIVATE(ipw)
1025 57029528 : do ipw=1+ipws,npw1+ipws
1026 56563288 : if(kinpw1(ipw-ipws)<threshold)then
1027 56097048 : gh1c(1,ipw)=gh1c(1,ipw)+gvnl2(1,ipw)
1028 56097048 : gh1c(2,ipw)=gh1c(2,ipw)+gvnl2(2,ipw)
1029 : end if
1030 : end do
1031 : end do
1032 : #ifdef HAVE_OPENMP_OFFLOAD
1033 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) call ompgpu_exit_map_delete(gvnl2,2*npw1*my_nspinor*ndat)
1034 : #endif
1035 466240 : ABI_FREE(gvnl2)
1036 : end if
1037 :
1038 : #ifdef HAVE_OPENMP_OFFLOAD
1039 : if(gs_hamkq%gpu_option==ABI_GPU_OPENMP) then
1040 : !$OMP TARGET EXIT DATA MAP(from:gh1c) IF(map_gh1c)
1041 : !$OMP TARGET EXIT DATA MAP(from:gs1c) IF(map_gs1c)
1042 : !$OMP TARGET EXIT DATA MAP(from:cwave) IF(map_cwave)
1043 : !$OMP TARGET UPDATE FROM(gvnlx1_) IF(usevnl==1 .and. gs_hamkq%gpu_option==ABI_GPU_OPENMP)
1044 : if(map_gvnlx1_) call ompgpu_exit_map_delete(gvnlx1_,2*npw1*my_nspinor*ndat)
1045 : end if
1046 : #endif
1047 8601471 : if (usevnl==1) then
1048 4665741 : nullify(gvnlx1_)
1049 : else
1050 3935730 : ABI_FREE(gvnlx1_)
1051 : end if
1052 :
1053 : ABI_NVTX_END_RANGE()
1054 :
1055 8601471 : call timab(196+tim_getgh1c,2,tsec)
1056 :
1057 : DBG_EXIT("COLL")
1058 :
1059 : !NOTE: These subroutines were placed here because NVHPC would cause issues when
1060 : ! processing OpenMP directives below by creating artifacts :
1061 : !293, Generating target enter data map(to: kinpw1$sd50(:))
1062 : ! Generating target enter data map(create: work(:,:,:,:))
1063 : ! Generating target enter data map(to: cwave$sd(:),dkinpw$sd46(:),gvnlx1_$sd(:))
1064 : !
1065 : ! Tested with NVHPC 24.3
1066 : ! Issues occurred with arrays: work, gvnlx1, gvnl2, nonlop_out
1067 : ! The name 'getgh1c_array' is for debugging purposes.
1068 :
1069 : #ifdef HAVE_OPENMP_OFFLOAD
1070 : contains
1071 : subroutine ompgpu_enter_map_alloc(getgh1c_array,size)
1072 : integer,intent(in) :: size
1073 : real(dp),intent(inout) :: getgh1c_array(size)
1074 : !$OMP TARGET ENTER DATA MAP(alloc:getgh1c_array)
1075 : end subroutine ompgpu_enter_map_alloc
1076 :
1077 : subroutine ompgpu_exit_map_delete(getgh1c_array,size)
1078 : integer,intent(in) :: size
1079 : real(dp),intent(inout) :: getgh1c_array(size)
1080 : !$OMP TARGET EXIT DATA MAP(delete:getgh1c_array)
1081 : end subroutine ompgpu_exit_map_delete
1082 :
1083 : #endif
1084 :
1085 17202942 : end subroutine getgh1c
1086 : !!***
1087 :
1088 : !----------------------------------------------------------------------
1089 :
1090 : !!****f* m_hamiltonian/rf_transgrid_and_pack
1091 : !! NAME
1092 : !! rf_transgrid_and_pack
1093 : !!
1094 : !! FUNCTION
1095 : !! Set up local potential vlocal1 with proper dimensioning, from vtrial1
1096 : !! taking into account the spin. Same thing for vlocal from vtrial.
1097 : !!
1098 : !! INPUTS
1099 : !! isppol=Spin index.
1100 : !! nspden=Number of density components
1101 : !! usepaw=1 if PAW, 0 for NC.
1102 : !! cplex=1 if DFPT potential is real, 2 for complex
1103 : !! nfftf=Number of FFT points on the FINE grid treated by this processor
1104 : !! nfft=Number of FFT points on the COARSE grid treated by this processor
1105 : !! ngfft(18)=Info on the coarse grid.
1106 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
1107 : !! mpi_enreg=information about MPI parallelization
1108 : !! vtrial(nfftf,nspden)=GS Vtrial(r) on the DENSE mesh
1109 : !! vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r) on the DENSE mesh
1110 : !!
1111 : !! OUTPUT
1112 : !! vlocal(n4,n5,n6,nvloc)= GS local potential in real space, on the augmented coarse fft grid
1113 : !! vlocal1(cplex*n4,n5,n6,nvloc)= RF local potential in real space, on the augmented coarse fft grid
1114 : !!
1115 : !! SOURCE
1116 :
1117 385929 : subroutine rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,nfft,ngfft,nvloc,&
1118 385929 : & pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1)
1119 :
1120 : !Arguments ------------------------------------
1121 : !scalars
1122 : integer,intent(in) :: isppol,nspden,usepaw,cplex,nfftf,nfft,nvloc
1123 : type(pawfgr_type),intent(in) :: pawfgr
1124 : type(MPI_type),intent(in) :: mpi_enreg
1125 : !arrays
1126 : integer,intent(in) :: ngfft(18)
1127 : real(dp),intent(in),target :: vtrial(nfftf,nspden)
1128 : real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden)
1129 : real(dp),intent(out) :: vlocal(ngfft(4),ngfft(5),ngfft(6),nvloc)
1130 : real(dp),intent(out) :: vlocal1(cplex*ngfft(4),ngfft(5),ngfft(6),nvloc)
1131 :
1132 : !Local variables-------------------------------
1133 : !scalars
1134 : integer :: n1,n2,n3,n4,n5,n6,paral_kgb,ispden
1135 : !arrays
1136 : real(dp) :: rhodum(1) !, tsec(2)
1137 385929 : real(dp), contiguous, pointer :: vtrial_ptr(:,:),vtrial1_ptr(:,:)
1138 385929 : real(dp),allocatable :: cgrvtrial(:,:),cgrvtrial1(:,:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:)
1139 : ! *************************************************************************
1140 :
1141 385929 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
1142 385929 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
1143 385929 : paral_kgb = mpi_enreg%paral_kgb
1144 :
1145 385929 : if (nspden/=4) then
1146 385634 : vtrial_ptr => vtrial
1147 385634 : if (usepaw==0.or.pawfgr%usefinegrid==0) then
1148 382671 : call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal(:,:,:,1),2)
1149 382671 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1,vlocal1(:,:,:,1),2)
1150 : else
1151 11852 : ABI_MALLOC(cgrvtrial,(nfft,nspden))
1152 2963 : call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr)
1153 2963 : call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal(:,:,:,1),2)
1154 11852 : ABI_REMALLOC(cgrvtrial, (cplex*nfft, nspden))
1155 2963 : call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial1)
1156 2963 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,cgrvtrial,vlocal1(:,:,:,1),2)
1157 2963 : ABI_FREE(cgrvtrial)
1158 : end if
1159 385929 : nullify(vtrial_ptr)
1160 : else
1161 : ! nspden==4 non-collinear magnetism
1162 295 : vtrial_ptr => vtrial
1163 295 : vtrial1_ptr => vtrial1
1164 1475 : ABI_MALLOC(vlocal_tmp,(n4,n5,n6))
1165 1475 : ABI_MALLOC(vlocal1_tmp,(cplex*n4,n5,n6))
1166 295 : if (usepaw==0.or.pawfgr%usefinegrid==0) then
1167 1475 : do ispden=1,nspden
1168 1180 : call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2)
1169 8572372 : vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
1170 1180 : call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2)
1171 13762427 : vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:)
1172 : end do
1173 : else
1174 : ! TODO FR EB check the correctness of the following lines for PAW calculations
1175 0 : ABI_MALLOC(cgrvtrial,(nfft,nspden))
1176 0 : ABI_MALLOC(cgrvtrial1,(nfft,nspden))
1177 0 : call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr)
1178 0 : call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial1,vtrial1_ptr)
1179 0 : do ispden=1,nspden
1180 0 : call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2)
1181 0 : vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
1182 0 : call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2)
1183 0 : vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:)
1184 : end do
1185 0 : ABI_FREE(cgrvtrial)
1186 : end if
1187 295 : ABI_FREE(vlocal_tmp)
1188 295 : ABI_FREE(vlocal1_tmp)
1189 : end if ! nspden
1190 :
1191 771858 : end subroutine rf_transgrid_and_pack
1192 : !!***
1193 :
1194 : !----------------------------------------------------------------------
1195 :
1196 : !!****f* m_hamiltonian/getgh1c_setup
1197 : !! NAME
1198 : !! getgh1c_setup
1199 : !!
1200 : !! FUNCTION
1201 : !!
1202 : !! INPUTS
1203 : !!
1204 : !! OUTPUT
1205 : !!
1206 : !! SOURCE
1207 :
1208 916828 : subroutine getgh1c_setup(gs_hamkq, rf_hamkq, dtset, psps, kpoint, kpq, idir, ipert,& ! In
1209 : natom, rmet, gprimd, gmet, istwf_k, npw_k, npw1_k, & ! In
1210 916828 : useylmgr1, kg_k, ylm_k, kg1_k, ylm1_k, ylmgr1_k, & ! In
1211 : dkinpw, nkpg, nkpg1, kpg_k, kpg1_k, kinpw1, ffnlk, ffnl1, ph3d, ph3d1,& ! Out
1212 : ddkinpw, dkinpw2, rf_hamk_dir2, ffnl1_test, & ! Optional
1213 : reuse_kpg_k, reuse_kpg1_k, reuse_ffnlk, reuse_ffnl1) ! Optional
1214 :
1215 : !Arguments ------------------------------------
1216 : !scalars
1217 : integer,intent(in) :: idir,ipert,istwf_k,npw_k,npw1_k,natom,useylmgr1
1218 : integer,intent(out) :: nkpg,nkpg1
1219 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
1220 : type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
1221 : type(rf_hamiltonian_type),intent(inout),optional :: rf_hamk_dir2
1222 : type(dataset_type),intent(in) :: dtset
1223 : type(pseudopotential_type),intent(in) :: psps
1224 : integer,optional,intent(in) :: reuse_kpg_k, reuse_kpg1_k, reuse_ffnlk, reuse_ffnl1
1225 : !arrays
1226 : integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k)
1227 : real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3)
1228 : real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
1229 : real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
1230 : real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)
1231 : real(dp),allocatable,intent(out) :: dkinpw(:),kinpw1(:)
1232 : real(dp),allocatable,intent(inout) :: ffnlk(:,:,:,:),ffnl1(:,:,:,:), kpg_k(:,:), kpg1_k(:,:)
1233 : real(dp),allocatable,intent(out),optional :: dkinpw2(:),ddkinpw(:),ffnl1_test(:,:,:,:)
1234 : real(dp),allocatable,intent(out) :: ph3d(:,:,:),ph3d1(:,:,:)
1235 :
1236 : !Local variables-------------------------------
1237 : !scalars
1238 : integer :: dimffnl1,dimffnlk,ider,idir0,idir1,idir2,istr,ntypat,print_info
1239 : integer :: reuse_ffnlk_, reuse_ffnl1_, reuse_kpg_k_, reuse_kpg1_k_
1240 : logical :: qne0
1241 : !real(dp) :: cpu, wall, gflops
1242 : !arrays
1243 : real(dp) :: ylmgr_dum(1,1,1), tsec(2)
1244 : ! *************************************************************************
1245 :
1246 : ! MG: This routine is called **many times** in the EPH code for phonon and DDK perturbations
1247 : ! Please, be extremely careful when adding extra stuff that may affect performance.
1248 :
1249 : ! Keep track of total time spent in getgh1c_setup (use 195 slot)
1250 916828 : call timab(195, 1, tsec)
1251 : !call cwtime(cpu, wall, gflops, "start")
1252 :
1253 916828 : reuse_ffnlk_ = 0; if (present(reuse_ffnlk)) reuse_ffnlk_ = reuse_ffnlk
1254 916828 : reuse_ffnl1_ = 0; if (present(reuse_ffnl1)) reuse_ffnl1_ = reuse_ffnl1
1255 916828 : reuse_kpg_k_ = 0; if (present(reuse_kpg_k)) reuse_kpg_k_ = reuse_kpg_k
1256 916828 : reuse_kpg1_k_ = 0; if (present(reuse_kpg1_k)) reuse_kpg1_k_ = reuse_kpg1_k
1257 :
1258 916828 : if(.not.present(ddkinpw) .and. ipert==natom+10) then
1259 0 : ABI_BUG("ddkinpw is not optional for ipert=natom+10.")
1260 : end if
1261 916828 : if(.not.present(dkinpw2) .and. ipert==natom+10 .and. idir>3) then
1262 0 : ABI_BUG("dkinpw2 is not optional for ipert=natom+10 and idir>3.")
1263 : end if
1264 916828 : if(.not.present(rf_hamk_dir2) .and. ((ipert==natom+10 .and. idir>3) .or. ipert==natom+11)) then
1265 0 : ABI_BUG("rf_hamk_dir2 is not optional for ipert=natom+10 (with idir>3) or ipert=natom+11.")
1266 : end if
1267 :
1268 916828 : ntypat = psps%ntypat
1269 916828 : qne0 = ((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14)
1270 :
1271 : ! Compute k+G vectors
1272 916828 : nkpg = 0; if (ipert >= 1 .and. ipert <= natom) nkpg = 3*dtset%nloalg(3)
1273 916828 : if (reuse_kpg_k_ == 0) then
1274 3667312 : ABI_MALLOC(kpg_k, (npw_k, nkpg))
1275 916828 : if (nkpg > 0) call mkkpg(kg_k, kpg_k, kpoint, nkpg, npw_k)
1276 : else
1277 0 : ABI_CHECK(all(shape(kpg_k) == [npw_k, nkpg]), sjoin("Wrong shape in input kpg_k", ltoa(shape(kpg_k))))
1278 : endif
1279 :
1280 : ! Compute k+q+G vectors
1281 916828 : nkpg1 = 0; if (ipert >= 1 .and. ipert <= natom) nkpg1 = 3*dtset%nloalg(3)
1282 916828 : if (reuse_kpg1_k_ == 0) then
1283 3667312 : ABI_MALLOC(kpg1_k, (npw1_k, nkpg1))
1284 916828 : if (nkpg1 > 0) call mkkpg(kg1_k, kpg1_k, kpq(:), nkpg1, npw1_k)
1285 : else
1286 0 : ABI_CHECK(all(shape(kpg1_k) == [npw1_k, nkpg1]), sjoin("Wrong shape in input kpg1_k:", ltoa(shape(kpg1_k))))
1287 : endif
1288 :
1289 : ! ===== Preparation of the non-local contributions
1290 916828 : dimffnlk =0; if (ipert<=natom) dimffnlk=1
1291 :
1292 : ! Compute nonlocal form factors ffnlk at (k+G)
1293 : ! (only for atomic displacement perturbation)
1294 916828 : if (reuse_ffnlk_ == 0) then
1295 5210208 : ABI_MALLOC(ffnlk, (npw_k, dimffnlk, psps%lmnmax, ntypat))
1296 868368 : if (ipert<=natom) then
1297 510820 : ider=0;idir0=0
1298 : call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,&
1299 : gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
1300 : psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,&
1301 510820 : psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
1302 : end if
1303 : else
1304 242300 : ABI_CHECK(all(shape(ffnlk) == [npw_k, dimffnlk, psps%lmnmax, ntypat]), sjoin("Wrong shape in input ffnlk:", ltoa(shape(ffnlk))))
1305 : end if
1306 :
1307 : ! Compute nonlocal form factors ffnl1 at (k+q+G)
1308 : !-- Atomic displacement perturbation
1309 916828 : if (ipert<=natom) then
1310 510820 : ider=0;idir0=0
1311 : !-- k-point perturbation (1st-derivative)
1312 406008 : else if (ipert==natom+1) then
1313 223022 : ider=1;idir0=idir
1314 : !-- k-point perturbation (2nd-derivative)
1315 182986 : else if (ipert==natom+10.or.ipert==natom+11) then
1316 14616 : ider=2;idir0=4
1317 : !-- Electric field perturbation
1318 168370 : else if (ipert==natom+2) then
1319 89874 : if (psps%usepaw==1) then
1320 8464 : ider=1;idir0=idir
1321 : else
1322 81410 : ider=0;idir0=0
1323 : end if
1324 : !-- Strain perturbation
1325 78496 : else if (ipert==natom+3.or.ipert==natom+4) then
1326 76900 : if (ipert==natom+3) istr=idir
1327 76900 : if (ipert==natom+4) istr=idir+3
1328 76900 : ider=1;idir0=-istr
1329 : !-- Magnetic field perturbation ( SPr, Zeeman ) or scalar potential
1330 1596 : else if(ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11))then
1331 1596 : ider=0;idir0=0
1332 : end if
1333 :
1334 : ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
1335 916828 : dimffnl1=1+ider
1336 916828 : if (ider==1.and.idir0==0) dimffnl1=2+2*psps%useylm
1337 916828 : if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm
1338 :
1339 916828 : if (reuse_ffnl1_ == 0) then
1340 5210208 : ABI_MALLOC(ffnl1, (npw1_k, dimffnl1, psps%lmnmax, ntypat))
1341 :
1342 : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,&
1343 : psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
1344 868368 : npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
1345 : else
1346 242300 : ABI_CHECK(all(shape(ffnl1) == [npw1_k, dimffnl1, psps%lmnmax, ntypat]), sjoin("Wrong shape in input ffnl1", ltoa(shape(ffnl1))))
1347 : end if
1348 :
1349 : ! Compute ffnl for nonlop with signs = 1
1350 916828 : print_info = 0
1351 916828 : if (dtset%prtvol==-19.or.dtset%prtvol==-20.or.dtset%prtvol==-21.or.dtset%nonlinear_info>=3) then
1352 240 : print_info = 1
1353 : end if
1354 916828 : if (present(ffnl1_test).and.print_info/=0.and.(ipert==natom+10.or.ipert==natom+11)) then
1355 1440 : ABI_MALLOC(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
1356 240 : idir0 = 0 ! for nonlop with signs = 1
1357 : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
1358 : psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
1359 240 : npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
1360 : end if
1361 :
1362 : !===== Preparation of the kinetic contributions
1363 : ! Note that not all these arrays should be allocated in the general case when wtk_k vanishes
1364 :
1365 : ! Compute (1/2) (2 Pi)**2 (k+q+G)**2:
1366 2750484 : ABI_MALLOC(kinpw1, (npw1_k))
1367 126382525 : kinpw1(:)=zero
1368 916828 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
1369 :
1370 2750484 : ABI_MALLOC(dkinpw,(npw_k)) ! 1st derivative (1st direction)
1371 126374369 : dkinpw(:)=zero
1372 916828 : if(ipert==natom+10 .and. idir>3) then
1373 14112 : ABI_MALLOC(dkinpw2,(npw_k)) ! 1st derivative (2nd directions)
1374 990680 : dkinpw2(:)=zero
1375 : end if
1376 916828 : if(ipert==natom+10) then
1377 22320 : ABI_MALLOC(ddkinpw,(npw_k)) ! 2nd derivative
1378 1524172 : ddkinpw(:)=zero
1379 : end if
1380 :
1381 : ! -- k-point perturbation (1st-derivative)
1382 916828 : if (ipert==natom+1) then
1383 : ! Compute the derivative of the kinetic operator vs k
1384 223022 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative
1385 : end if
1386 :
1387 : !-- k-point perturbation (2nd-derivative)
1388 916828 : if (ipert==natom+10.or.ipert==natom+11) then
1389 : ! Compute the derivative of the kinetic operator vs k in kinpw, second and first orders
1390 14616 : if(ipert==natom+10 .and. idir<=3) then
1391 4104 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative
1392 4104 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir,idir) ! 2nd derivative
1393 : else
1394 11072 : select case(idir)
1395 : ! Diagonal terms:
1396 : case (1)
1397 560 : idir1 = 1
1398 560 : idir2 = 1
1399 : case (2)
1400 560 : idir1 = 2
1401 560 : idir2 = 2
1402 : case (3)
1403 32 : idir1 = 3
1404 32 : idir2 = 3
1405 : ! Upper triangular terms:
1406 : case (4)
1407 2104 : idir1 = 2
1408 2104 : idir2 = 3
1409 : case (5)
1410 1048 : idir1 = 1
1411 1048 : idir2 = 3
1412 : case (6)
1413 2104 : idir1 = 1
1414 2104 : idir2 = 2
1415 : ! Lower triangular terms:
1416 : case (7)
1417 1544 : idir1 = 3
1418 1544 : idir2 = 2
1419 : case (8)
1420 1016 : idir1 = 3
1421 1016 : idir2 = 1
1422 : case (9)
1423 1544 : idir1 = 2
1424 10512 : idir2 = 1
1425 : end select
1426 10512 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) ! 1st derivative, idir1
1427 10512 : if(ipert==natom+10) then
1428 7056 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw2,kpoint,npw_k,idir2,0) ! 1st derivative, idir2
1429 7056 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir1,idir2) ! 2nd derivative
1430 : end if
1431 : end if
1432 : end if
1433 :
1434 : !-- Strain perturbation
1435 916828 : if (ipert==natom+3.or.ipert==natom+4) then
1436 76900 : if (ipert==natom+3) istr=idir
1437 76900 : if (ipert==natom+4) istr=idir+3
1438 : ! Compute the derivative of the kinetic operator vs strain
1439 76900 : call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,kg_k,kpoint,npw_k)
1440 : end if
1441 :
1442 : !===== Load the k/k+q dependent parts of the Hamiltonian
1443 : ! Load k-dependent part in the Hamiltonian datastructure
1444 3667312 : ABI_MALLOC(ph3d, (2,npw_k,gs_hamkq%matblk))
1445 : call gs_hamkq%load_k(kpt_k=kpoint, npw_k=npw_k, istwf_k=istwf_k, kg_k=kg_k, kpg_k=kpg_k,&
1446 916828 : ph3d_k=ph3d, compute_ph3d=.true., compute_gbound=.true.)
1447 :
1448 4584140 : if (size(ffnlk)>0) then
1449 510820 : call gs_hamkq%load_k(ffnl_k=ffnlk)
1450 : else
1451 406008 : call gs_hamkq%load_k(ffnl_k=ffnl1)
1452 : end if
1453 :
1454 : ! Load k+q-dependent part in the Hamiltonian datastructure
1455 : ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
1456 : call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
1457 916828 : kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,compute_gbound=.true.)
1458 :
1459 916828 : if (qne0) then
1460 1639484 : ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk))
1461 409871 : call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.)
1462 : end if
1463 :
1464 : ! Load k-dependent part in the 1st-order Hamiltonian datastructure
1465 916828 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
1466 :
1467 916828 : if (ipert==natom+10) then
1468 11160 : call rf_hamkq%load_k(ddkinpw_k=ddkinpw)
1469 11160 : if (idir>3) call rf_hamk_dir2%load_k(dkinpw_k=dkinpw2,ddkinpw_k=ddkinpw)
1470 : end if
1471 :
1472 916828 : call timab(195, 2, tsec)
1473 :
1474 916828 : end subroutine getgh1c_setup
1475 : !!***
1476 :
1477 : !!****f* ABINIT/getdc1
1478 : !!
1479 : !! NAME
1480 : !! getdc1
1481 : !!
1482 : !! FUNCTION
1483 : !! Compute |delta_C^(1)> from one wave function C - PAW ONLY
1484 : !! Compute <G|delta_C^(1)> (dcwavef) and eventually <P_i| delta_C^(1)> (dcwaveprj) where P_i= non-local projector
1485 : !! delta_C^(1) is the variation of wavefunction only due to variation of overlap operator S.
1486 : !! delta_C^(1)=-1/2.Sum_j [ <C_j|S^(1)|C>.C_j
1487 : !! see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42 and 40, term 2)
1488 : !!
1489 : !! INPUTS
1490 : !! cgq(2,mcgq)=wavefunction coefficients for all bands j on present processor, at k+Q: cgq=< G |Cnk+q>
1491 : !! cprjq(natom,mcprjq)= wave functions j at k+q projected with non-local projectors: cprjq=<P_i|Cnk+q>
1492 : !! ibgq=shift to be applied on the location of data in the array cprjq
1493 : !! icgq=shift to be applied on the location of data in the array cgq
1494 : !! istwfk=option parameter that describes the storage of wfs
1495 : !! mcgq=second dimension of the cgq array
1496 : !! mcprjq=second dimension of the cprjq array
1497 : !! mpi_enreg=information about MPI parallelization
1498 : !! ndat=number of bands to compute in parallel
1499 : !! natom= number of atoms in cell
1500 : !! nband=number of bands
1501 : !! npw1=number of planewaves in basis sphere at k+Q
1502 : !! nspinor=number of spinorial components of the wavefunctions
1503 : !! opt_cprj=flag governing the computation of <P_i|delta_C^(1)> (P_i= non-local projector)
1504 : !! s1cwave0(2,npw1*nspinor)=<G|S^(1)|C> where S^(1) is the first-order overlap operator
1505 : !!
1506 : !! OUTPUT
1507 : !! dcwavef(2,npw1*nspinor)=change of wavefunction due to change of overlap PROJECTED ON PLANE-WAVES:
1508 : !! dcwavef is delta_C(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
1509 : !! === if optcprj=1 ===
1510 : !! dcwaveprj(natom,nspinor*optcprj)=change of wavefunction due to change of overlap PROJECTED ON NL-PROJECTORS:
1511 : !!
1512 : !! SOURCE
1513 :
1514 816088 : subroutine getdc1(band,band_procs,bands_treated_now,cgq,cprjq,dcwavef,dcwaveprj,&
1515 : & ibgq,icgq,istwfk,mcgq,mcprjq,&
1516 816088 : & mpi_enreg,ndat,natom,nband,nband_me,npw1,nspinor,optcprj,s1cwave0,&
1517 : & gpu_option)
1518 :
1519 : !Arguments ------------------------------------
1520 : !scalars
1521 : integer,intent(in) :: ibgq,icgq,istwfk,mcgq,mcprjq,natom,nband,npw1,nspinor,optcprj
1522 : integer,intent(in) :: band, nband_me, ndat
1523 : type(MPI_type),intent(in) :: mpi_enreg
1524 : integer,optional,intent(in) :: gpu_option
1525 : !arrays
1526 : integer,intent(in) :: band_procs(nband),bands_treated_now(nband,ndat)
1527 : real(dp),intent(in) :: cgq(2,mcgq),s1cwave0(2,npw1*nspinor*ndat)
1528 : real(dp),intent(out) :: dcwavef(2,npw1*nspinor*ndat)
1529 : type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq)
1530 : type(pawcprj_type),intent(inout) :: dcwaveprj(natom,nspinor*ndat*optcprj)
1531 :
1532 : !Local variables-------------------------------
1533 : !scalars
1534 : integer, parameter :: tim_projbd=0
1535 : integer :: ipw, idat
1536 : integer :: band_, ierr, nproc_band
1537 : integer :: l_gpu_option
1538 : real(dp),parameter :: scal=-half
1539 : !arrays
1540 816088 : integer, allocatable :: nlmn(:)
1541 816088 : real(dp), allocatable :: dummy(:,:),scprod(:,:,:)
1542 : real(dp), allocatable :: dcwavef_tmp(:,:)
1543 816088 : type(pawcprj_type),allocatable :: dcwaveprj_tmp(:,:)
1544 : ! *********************************************************************
1545 :
1546 : DBG_ENTER("COLL")
1547 :
1548 816088 : l_gpu_option = ABI_GPU_DISABLED; if (present(gpu_option)) l_gpu_option = gpu_option
1549 816088 : nproc_band=mpi_enreg%nproc_band
1550 :
1551 816088 : ABI_MALLOC(dummy,(0,0))
1552 3264352 : ABI_MALLOC(scprod,(2,nband_me,ndat))
1553 2448264 : ABI_MALLOC(dcwavef_tmp,(2,npw1*nspinor*ndat))
1554 : #ifdef HAVE_OPENMP_OFFLOAD
1555 : !$OMP TARGET ENTER DATA MAP(alloc:dcwavef_tmp,scprod) IF(l_gpu_option==ABI_GPU_OPENMP)
1556 : #endif
1557 816088 : if (optcprj == 1) then
1558 844092 : ABI_MALLOC(nlmn,(natom))
1559 2100264 : ABI_MALLOC(dcwaveprj_tmp,(natom,nspinor*ndat*optcprj))
1560 974808 : nlmn(:)=dcwaveprj(:,1)%nlmn
1561 281364 : call pawcprj_alloc(dcwaveprj_tmp, 0, nlmn)
1562 281364 : ABI_FREE(nlmn)
1563 : end if
1564 :
1565 : !=== 1- COMPUTE: <G|S^(1)|C_k> - Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>]
1566 : !! using the projb routine
1567 : !Note the subtlety: projbd is called with useoverlap=0 and s1cwave0
1568 : !in order to get Sum[<cgq|s1|c>|cgq>]=Sum[<cgq|gs1>|cgq>]
1569 :
1570 : ! run over procs in my pool which have a dcwavef to projbd
1571 6301096 : do band_ = 1, nband, ndat
1572 5485008 : if (bands_treated_now(band_, 1) == 0) cycle
1573 825088 : if(l_gpu_option==ABI_GPU_DISABLED) then
1574 411686110 : dcwavef_tmp = zero
1575 0 : else if(l_gpu_option==ABI_GPU_OPENMP) then
1576 0 : call gpu_set_to_zero(dcwavef_tmp,int(2,c_size_t)*npw1*nspinor*ndat)
1577 : end if
1578 :
1579 : ! distribute dcwavef_tmp to my band pool
1580 : ! everyone works on a single band s1cwave0 = <G|S^(1)|C_k>
1581 825088 : if (band_ == band) then
1582 816088 : if(l_gpu_option==ABI_GPU_DISABLED) then
1583 : !$OMP PARALLEL DO
1584 124296762 : do ipw=1,npw1*nspinor*ndat
1585 371258110 : dcwavef_tmp(1:2,ipw)=s1cwave0(1:2,ipw)
1586 : end do
1587 : else if(l_gpu_option==ABI_GPU_OPENMP) then
1588 : #ifdef HAVE_OPENMP_OFFLOAD
1589 : call gpu_copy(dcwavef_tmp, s1cwave0, int(2,c_size_t)*npw1*nspinor*ndat)
1590 : #endif
1591 : end if
1592 : end if
1593 : #ifdef HAVE_OPENMP_OFFLOAD
1594 : !$OMP TARGET UPDATE FROM (dcwavef_tmp) IF(l_gpu_option==ABI_GPU_OPENMP)
1595 : #endif
1596 825088 : call xmpi_bcast(dcwavef_tmp,band_procs(band_),mpi_enreg%comm_band,ierr)
1597 : #ifdef HAVE_OPENMP_OFFLOAD
1598 : !$OMP TARGET UPDATE TO (dcwavef_tmp) IF(l_gpu_option==ABI_GPU_OPENMP)
1599 : #endif
1600 :
1601 : ! get the projbd onto my processor's bands dcwavef = dcwavef - <cgq|dcwavef>|cgq>
1602 : ! dcwavef = <G|S^(1)|C_k> - Sum_{MYj} [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>]
1603 : ! scprod = <C_k+q,j|S^(1)|C_k> for {MYj}
1604 1650176 : do idat=1,ndat
1605 : call projbd(cgq,dcwavef_tmp(:,1+(idat-1)*npw1*nspinor:idat*npw1*nspinor),-1,&
1606 : & icgq,0,istwfk,mcgq,0,nband_me,npw1,nspinor,&
1607 : & dummy,scprod(:,:,idat),0,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft,&
1608 1650176 : & gpu_option=l_gpu_option)
1609 : end do
1610 : #ifdef HAVE_OPENMP_OFFLOAD
1611 : !$OMP TARGET UPDATE FROM (scprod) IF(l_gpu_option==ABI_GPU_OPENMP)
1612 : #endif
1613 :
1614 :
1615 : ! sum all of the corrections
1616 : ! dcwavef = Nprocband * <G|S^(1)|C_k> - Sum_{ALLj} [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>]
1617 825088 : call xmpi_sum(dcwavef_tmp,mpi_enreg%comm_band,ierr,use_omp_map=(l_gpu_option==ABI_GPU_OPENMP))
1618 :
1619 : ! save to my proc if it is my turn, and subtract Ntuple counted dcwavef
1620 825088 : if (band_ == band) then
1621 : !=== 2- COMPUTE: <G|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] by subtraction
1622 : ! tested this is equivalent to previous coding to within 1.e-18 accumulated error (probably in favor of this coding)
1623 816088 : if(l_gpu_option==ABI_GPU_DISABLED) then
1624 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(dcwavef,s1cwave0,dcwavef_tmp,npw1,nspinor)
1625 124296762 : do ipw=1,npw1*nspinor*ndat
1626 371258110 : dcwavef(1:2,ipw)= scal*(nproc_band*s1cwave0(1:2,ipw)-dcwavef_tmp(1:2,ipw))
1627 : end do
1628 : else if(l_gpu_option==ABI_GPU_OPENMP) then
1629 : #ifdef HAVE_OPENMP_OFFLOAD
1630 : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO PRIVATE(ipw) MAP(to:dcwavef,s1cwave0,dcwavef_tmp)
1631 : do ipw=1,npw1*nspinor*ndat
1632 : dcwavef(1,ipw)= scal*(nproc_band*s1cwave0(1,ipw)-dcwavef_tmp(1,ipw))
1633 : dcwavef(2,ipw)= scal*(nproc_band*s1cwave0(2,ipw)-dcwavef_tmp(2,ipw))
1634 : end do
1635 : #endif
1636 : end if
1637 : end if
1638 : !print*,"dcwavef ", dcwavef(1,1)
1639 :
1640 : !=== 3- COMPUTE: <P_i|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<P_i|C_k+q,j>]
1641 : ! as above everyone has to operate on each band band_
1642 2466264 : if (optcprj==1.and.mcprjq>0) then
1643 : ! cprjq = <P_i|C_k+q,j> for MYj
1644 : ! dcwaveprj_tmp = Sum_MYj [<C_k+q,j|S^(1)|C_k>.<P_i|C_k+q,j>]
1645 574392 : do idat=1,ndat
1646 574392 : call pawcprj_lincom(scprod(:,:,idat),cprjq(:,ibgq+1:ibgq+nspinor*nband_me),dcwaveprj_tmp(:,1+(idat-1)*nspinor:idat*nspinor),nband_me)
1647 : end do
1648 :
1649 : ! still need to mpisum the dcwaveprj to get linear combination of all bands, not just mine
1650 : ! dcwaveprj = Sum_ALLj [<C_k+q,j|S^(1)|C_k,i>.<P_i|C_k+q,j>]
1651 287196 : call pawcprj_mpi_sum(dcwaveprj_tmp,mpi_enreg%comm_band,ierr)
1652 :
1653 287196 : if (band_ == band) then
1654 : ! dcwaveprj = -1/2 dcwaveprj_tmp
1655 : !TODO: check the correct order of scal and zero (alpha / beta) coefficients.
1656 : ! Here dcwaveprj is squashed by the _tmp variable which is used in parallel
1657 281364 : call pawcprj_axpby(scal,zero,dcwaveprj_tmp,dcwaveprj)
1658 : end if
1659 : end if
1660 :
1661 : end do ! procs in my band pool
1662 :
1663 : #ifdef HAVE_OPENMP_OFFLOAD
1664 : !$OMP TARGET EXIT DATA MAP(delete:dcwavef_tmp,scprod) IF(l_gpu_option==ABI_GPU_OPENMP)
1665 : #endif
1666 :
1667 816088 : ABI_FREE(dummy)
1668 816088 : ABI_FREE(scprod)
1669 816088 : ABI_FREE(dcwavef_tmp)
1670 816088 : if (optcprj == 1) then
1671 281364 : call pawcprj_free(dcwaveprj_tmp)
1672 974808 : ABI_FREE(dcwaveprj_tmp)
1673 : end if
1674 :
1675 : DBG_EXIT("COLL")
1676 :
1677 816088 : end subroutine getdc1
1678 : !!***
1679 :
1680 : !!****f* ABINIT/getgh1dqc
1681 : !! NAME
1682 : !! getgh1dqc
1683 : !!
1684 : !! FUNCTION
1685 : !! Computes <G|dH^(1)/dq_{gamma}|C> or <G|d^2H^(1)/dq_{gamma}dq_{delta}|C>
1686 : !! for input vector |C> expressed in reciprocal space.
1687 : !! dH^(1)/dq_{gamma} and d^2H^(1)/dq_{gamma}dq_{delta} are the first
1688 : !! and second q-gradient (at q=0) of the 1st-order perturbed Hamiltonian.
1689 : !! The first (second) derivative direction is inferred from idir (qdir1).
1690 : !!
1691 : !! INPUTS
1692 : !! cwave(2,npw*nspinor)=input wavefunction, in reciprocal space
1693 : !! cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives)
1694 : !! if not allocated or size=0, they are locally computed (and not sorted)
1695 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
1696 : !! idir=first index of the perturbation
1697 : !! ipert=type of the perturbation
1698 : !! mpi_enreg=information about MPI parallelization
1699 : !! npw=number of planewaves in basis sphere at given k.
1700 : !! npw1=number of planewaves in basis sphere at k+q
1701 : !! optlocal=0: local part of H^(1) is not computed
1702 : !! 1: local part of H^(1) is computed in gvloc1dqc
1703 : !! optnl=0: non-local part of H^(1) is not computed
1704 : !! 1: non-local part of H^(1) depending on VHxc^(1) is not computed in gvloc1dqc
1705 : !! 2: non-local part of H^(1) is totally computed in gvloc1dqc
1706 : !! qdir1= direction of the 1st q-gradient
1707 : !! rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,k+q
1708 : !! qdir2= (optional) direction of the 2nd q-gradient
1709 : !!
1710 : !! OUTPUT
1711 : !! gh1dqc(2,npw1*nspinor)= <G|dH^(1)/dq_{\gamma}|C> on the k+q sphere
1712 : !! gvloc1dqc(2,npw1*nspinor)= local potential part of gh1dqc
1713 : !! gvnl1dqc(2,npw1*nspinor)= non local potential part of gh1dqc
1714 : !!
1715 : !! NOTES
1716 : !! Currently two Hamiltonian gradients at (q=0) are implemented:
1717 : !! ipert<=natom -> first q-derivative along reduced coordinates directions
1718 : !! of the atomic displacement perturbation hamiltonian
1719 : !! ipert==natom+3 or natom+4 -> second q-derivative along cartesian coordinates
1720 : !! of the metric perturbation hamiltonian.
1721 : !! Which is equivalent (except for an i factor) to the first
1722 : !! q-derivative along cartesian coordinates of the strain
1723 : !! perturbation hamiltonian.
1724 : !!
1725 : !! SOURCE
1726 :
1727 1002208 : subroutine getgh1dqc(cwave,cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq,&
1728 : & idir,ipert,mpi_enreg,optlocal,optnl,qdir1,rf_hamkq,&
1729 : & qdir2) !optional
1730 :
1731 : !Arguments ------------------------------------
1732 : !scalars
1733 : integer,intent(in) :: idir,ipert,optlocal,optnl,qdir1
1734 : integer,intent(in),optional :: qdir2
1735 : type(MPI_type),intent(in) :: mpi_enreg
1736 : type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
1737 : type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq
1738 :
1739 : !arrays
1740 : real(dp),intent(inout) :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor)
1741 : real(dp),intent(out) :: gh1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
1742 : real(dp),intent(out) :: gvloc1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
1743 : real(dp),intent(out) :: gvnl1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
1744 1002208 : real(dp),contiguous, pointer :: dqdqkinpw(:),kinpw1(:)
1745 : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
1746 :
1747 : !Local variables-------------------------------
1748 : !scalars
1749 : integer :: choice,cpopt,iidir,ipw,ipws,ispinor,my_nspinor,natom,nnlout
1750 : integer :: npw,npw1,paw_opt,signs,tim_fourwf,tim_nonlop
1751 : logical :: has_kin
1752 : !character(len=500) :: msg
1753 : real(dp) :: lambda,weight
1754 :
1755 : !arrays
1756 : integer,parameter :: ngamma(3,3)=reshape((/1,6,5,9,2,4,8,7,3/),(/3,3/))
1757 : real(dp) :: enlout(1),svectout_dum(1,1)
1758 1002208 : real(dp),contiguous, pointer :: gvnl1dqc_(:,:)
1759 1002208 : real(dp), allocatable :: work(:,:,:,:)
1760 : ! *************************************************************************
1761 :
1762 : DBG_ENTER("COLL")
1763 :
1764 : !======================================================================
1765 : !== Initialisations and compatibility tests
1766 : !======================================================================
1767 :
1768 1002208 : npw =gs_hamkq%npw_k
1769 1002208 : npw1 =gs_hamkq%npw_kp
1770 1002208 : natom=gs_hamkq%natom
1771 :
1772 : !Compatibility tests
1773 1002208 : if (mpi_enreg%paral_spinor==1) then
1774 0 : ABI_BUG('Not compatible with parallelization over spinor components !')
1775 : end if
1776 :
1777 : !Check sizes
1778 1002208 : my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
1779 3006624 : if (size(cwave)<2*npw*my_nspinor) then
1780 0 : ABI_BUG('wrong size for cwave!')
1781 : end if
1782 3006624 : if (size(gh1dqc)<2*npw1*my_nspinor) then
1783 0 : ABI_BUG('wrong size for gh1dqc!')
1784 : end if
1785 :
1786 : !=============================================================================
1787 : !== Apply the q-gradients of the 1st-order local potential to the wavefunction
1788 : !=============================================================================
1789 :
1790 : !Phonon and metric (strain) perturbation
1791 1002208 : if (ipert<=natom+6.and.ipert/=natom+1.and.ipert/=natom+2.and.optlocal>0) then
1792 :
1793 5011040 : ABI_MALLOC(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
1794 :
1795 1002208 : weight=one ; tim_fourwf=4
1796 : call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gvloc1dqc,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
1797 : & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
1798 : & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,&
1799 1002208 : & gpu_option=gs_hamkq%gpu_option)
1800 :
1801 1002208 : ABI_FREE(work)
1802 :
1803 : else
1804 :
1805 : !$OMP PARALLEL DO
1806 0 : do ipw=1,npw1*my_nspinor
1807 0 : gvloc1dqc(:,ipw)=zero
1808 : end do
1809 :
1810 : end if
1811 :
1812 : !================================================================================
1813 : !== Apply the q-gradients of the 1st-order non-local potential to the wavefunction
1814 : !================================================================================
1815 :
1816 : !Initializations
1817 1002208 : lambda=zero
1818 1002208 : nnlout=1
1819 1002208 : tim_nonlop=0
1820 :
1821 : !Allocations
1822 3006624 : ABI_MALLOC(gvnl1dqc_,(2,npw1*my_nspinor))
1823 :
1824 : !Phonon perturbation
1825 : !-------------------------------------------
1826 : !1st q-gradient
1827 1002208 : if (ipert<=natom.and..not.present(qdir2).and.optnl>0) then
1828 380128 : cpopt=-1 ; choice=22 ; signs=2 ; paw_opt=0
1829 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
1830 760256 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,iatom_only=ipert,qdir=qdir1)
1831 :
1832 : !$OMP PARALLEL DO
1833 39529196 : do ipw=1,npw1*my_nspinor
1834 39149068 : gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw)
1835 39529196 : gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw)
1836 : end do
1837 :
1838 : !2nd q-gradient
1839 622080 : else if (ipert<=natom.and.present(qdir2).and.optnl>0) then
1840 248832 : iidir=ngamma(idir,qdir2)
1841 248832 : cpopt=-1 ; choice=25 ; signs=2 ; paw_opt=0
1842 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,iidir,(/lambda/),mpi_enreg,1,nnlout,&
1843 497664 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,iatom_only=ipert,qdir=qdir1)
1844 :
1845 : !$OMP PARALLEL DO
1846 24828768 : do ipw=1,npw1*my_nspinor
1847 24579936 : gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw)
1848 24828768 : gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw)
1849 : end do
1850 :
1851 : !Metric (strain) perturbation
1852 : !-------------------------------------------
1853 373248 : else if ((ipert==natom+3.or.ipert==natom+4).and.optnl>0) then
1854 373248 : cpopt=-1 ; choice=33 ; signs=2 ; paw_opt=0
1855 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
1856 746496 : & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,qdir=qdir1)
1857 :
1858 : !$OMP PARALLEL DO
1859 37243152 : do ipw=1,npw1*my_nspinor
1860 36869904 : gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw)
1861 37243152 : gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw)
1862 : end do
1863 :
1864 : else
1865 :
1866 : !$OMP PARALLEL DO
1867 0 : do ipw=1,npw1*my_nspinor
1868 0 : gvnl1dqc(:,ipw)=zero
1869 : end do
1870 :
1871 : end if
1872 :
1873 : !==============================================================================
1874 : !== Apply the q-gradients of the 1st-order kinetic operator to the wavefunction
1875 : !== (add it to nl contribution)
1876 : !==============================================================================
1877 :
1878 : !Strain (metric) perturbation
1879 : !-------------------------------------------
1880 1002208 : has_kin=(ipert==natom+3.or.ipert==natom+4)
1881 1002208 : if (associated(gs_hamkq%kinpw_kp)) then
1882 1002208 : kinpw1 => gs_hamkq%kinpw_kp
1883 0 : else if (has_kin) then
1884 0 : ABI_BUG('need kinpw1 allocated!')
1885 : end if
1886 1002208 : if (associated(rf_hamkq%dkinpw_k)) then
1887 1002208 : dqdqkinpw => rf_hamkq%dkinpw_k
1888 0 : else if (has_kin) then
1889 0 : ABI_BUG('need dqdqkinpw allocated!')
1890 : end if
1891 :
1892 1002208 : if (has_kin) then
1893 : ! Remember that npw=npw1
1894 746496 : do ispinor=1,my_nspinor
1895 : !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwave,ispinor,gvnl1dqc,dqdqkinpw,kinpw1,npw,my_nspinor)
1896 37616400 : do ipw=1,npw
1897 36869904 : ipws=ipw+npw*(ispinor-1)
1898 37243152 : if(kinpw1(ipw)<huge(zero)*1.d-11)then
1899 36869904 : gvnl1dqc(1,ipws)=gvnl1dqc(1,ipws)+dqdqkinpw(ipw)*cwave(1,ipws)
1900 36869904 : gvnl1dqc(2,ipws)=gvnl1dqc(2,ipws)+dqdqkinpw(ipw)*cwave(2,ipws)
1901 : else
1902 0 : gvnl1dqc(1,ipws)=zero
1903 0 : gvnl1dqc(2,ipws)=zero
1904 : end if
1905 : end do
1906 : end do
1907 : end if
1908 :
1909 : !===================================================================================
1910 : !== Sum contributions to get the application of dH^(1)/dq or d^2H^(1)/dqdq to the wf
1911 : !===================================================================================
1912 :
1913 2004416 : do ispinor=1,my_nspinor
1914 1002208 : ipws=(ispinor-1)*npw1
1915 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1dqc,gvnl1dqc,kinpw1,ipws,npw1)
1916 102603324 : do ipw=1+ipws,npw1+ipws
1917 101601116 : if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then
1918 100598908 : gh1dqc(1,ipw)=gvloc1dqc(1,ipw)+gvnl1dqc(1,ipw)
1919 100598908 : gh1dqc(2,ipw)=gvloc1dqc(2,ipw)+gvnl1dqc(2,ipw)
1920 : else
1921 0 : gh1dqc(1,ipw)=zero
1922 0 : gh1dqc(2,ipw)=zero
1923 : end if
1924 : end do
1925 : end do
1926 :
1927 1002208 : ABI_FREE(gvnl1dqc_)
1928 : DBG_EXIT("COLL")
1929 :
1930 1002208 : end subroutine getgh1dqc
1931 : !!***
1932 :
1933 : !!****f* m_hamiltonian/getgh1dqc_setup
1934 : !! NAME
1935 : !! getgh1dqc_setup
1936 : !!
1937 : !! FUNCTION
1938 : !!
1939 : !! INPUTS
1940 : !!
1941 : !!
1942 : !! OUTPUT
1943 : !!
1944 : !! SOURCE
1945 :
1946 250168 : subroutine getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpoint,kpq,idir,ipert,qdir1,& ! In
1947 : & natom,rmet,rprimd,gprimd,gmet,istwf_k,npw_k,npw1_k,nylmgr,& ! In
1948 250168 : & useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In
1949 : & nkpg,nkpg1,kpg_k,kpg1_k,dqdqkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out
1950 : & reuse_ffnlk,reuse_ffnl1,qdir2) ! Optional
1951 :
1952 : !Arguments ------------------------------------
1953 : !scalars
1954 : integer,intent(in) :: idir,ipert,istwf_k,natom,npw_k,npw1_k,nylmgr,qdir1,useylmgr1
1955 : integer,intent(in),optional :: reuse_ffnlk,reuse_ffnl1,qdir2
1956 : integer,intent(out) :: nkpg,nkpg1
1957 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
1958 : type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
1959 : type(dataset_type),intent(in) :: dtset
1960 : type(pseudopotential_type),intent(in) :: psps
1961 : !arrays
1962 : integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k)
1963 : real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1964 : real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
1965 : ! real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
1966 : real(dp),intent(in) :: ylmgr1_k(npw1_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
1967 : real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)
1968 : real(dp),allocatable,intent(out) :: dqdqkinpw(:),kinpw1(:)
1969 : real(dp),allocatable,intent(inout) :: ffnlk(:,:,:,:),ffnl1(:,:,:,:)
1970 : real(dp),allocatable,intent(out) :: kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:)
1971 :
1972 : !Local variables-------------------------------
1973 : !scalars
1974 : integer :: dimffnl1,dimffnlk,ider,idir0,ig,mu,mua,mub,ntypat
1975 : integer :: nu,nua,nub
1976 : integer :: reuse_ffnlk_,reuse_ffnl1_
1977 : logical :: qne0
1978 : !arrays
1979 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
1980 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
1981 : real(dp) :: ylmgr_dum(1,1,1)
1982 250168 : real(dp),allocatable :: ffnl1_tmp(:,:,:,:)
1983 : ! *************************************************************************
1984 :
1985 250168 : reuse_ffnlk_ = 0; if (present(reuse_ffnlk)) reuse_ffnlk_ = reuse_ffnlk
1986 250168 : reuse_ffnl1_ = 0; if (present(reuse_ffnl1)) reuse_ffnl1_ = reuse_ffnl1
1987 :
1988 250168 : ntypat = psps%ntypat
1989 250168 : qne0=((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14)
1990 :
1991 : !Compute (k+G) vectors
1992 250168 : nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3)
1993 1000672 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
1994 250168 : if (nkpg>0) then
1995 0 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1996 : end if
1997 :
1998 : !Compute (k+q+G) vectors
1999 250168 : nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3)
2000 1000672 : ABI_MALLOC(kpg1_k,(npw1_k,nkpg1))
2001 250168 : if (nkpg1>0) then
2002 0 : call mkkpg(kg1_k,kpg1_k,kpq(:),nkpg1,npw1_k)
2003 : end if
2004 :
2005 : !===== Preparation of the non-local contributions
2006 :
2007 250168 : dimffnlk=0;if (ipert<=natom) dimffnlk=1
2008 :
2009 : !Compute nonlocal form factors ffnlk at (k+G)
2010 250168 : if (reuse_ffnlk_ == 0) then
2011 493056 : ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,ntypat))
2012 82176 : if (ipert<=natom) then
2013 51072 : ider=0;idir0=0
2014 : call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,&
2015 : gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
2016 : psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,&
2017 51072 : psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
2018 : end if
2019 : else
2020 839960 : ABI_CHECK(all(shape(ffnlk) == [npw_k, dimffnlk, psps%lmnmax, ntypat]), "Wrong shape in input ffnlk")
2021 : end if
2022 :
2023 : !Compute nonlocal form factors ffnl1 at (k+q+G)
2024 : !TODO: For the second order gradients, this routine is called for each 3 directions of the
2025 : !derivative and every time it calculates all the form factors derivatives. This could be
2026 : !done just once.
2027 : !-- 1st q-grad of atomic displacement perturbation
2028 250168 : if (ipert<=natom.and..not.present(qdir2)) then
2029 94648 : ider=1;idir0=qdir1
2030 : !-- 2nd q-grad of atomic displacement perturbation
2031 155520 : else if (ipert<=natom.and.present(qdir2)) then
2032 62208 : ider=2;idir0=4
2033 : !-- 2nd q-grad of metric (1st q-grad of strain) perturbation
2034 93312 : else if (ipert==natom+3.or.ipert==natom+4) then
2035 93312 : ider=2;idir0=4
2036 : end if
2037 :
2038 : !Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
2039 250168 : dimffnl1=1+ider
2040 250168 : if (ider==2.and.(idir0==0.or.idir0==4)) dimffnl1=3+7*psps%useylm
2041 :
2042 250168 : if (reuse_ffnl1_ == 0) then
2043 493056 : ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,ntypat))
2044 : call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,&
2045 : psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
2046 82176 : npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
2047 : else
2048 839960 : ABI_CHECK(all(shape(ffnl1) == [npw1_k, dimffnl1, psps%lmnmax, ntypat]), "Wrong shape in input ffnl1")
2049 : end if
2050 :
2051 :
2052 : !Convert nonlocal form factors to cartesian coordinates.
2053 : !For metric (strain) perturbation only.
2054 250168 : if (ipert==natom+3.or.ipert==natom+4) then
2055 559872 : ABI_MALLOC(ffnl1_tmp,(npw1_k,dimffnl1,psps%lmnmax,ntypat))
2056 1118694240 : ffnl1_tmp=ffnl1
2057 :
2058 : !First q-derivative
2059 336494736 : ffnl1(:,2:4,:,:)=zero
2060 373248 : do mu=1,3
2061 28025676 : do ig=1,npw1_k
2062 110889648 : do nu=1,3
2063 1189054404 : ffnl1(ig,1+mu,:,:)=ffnl1(ig,1+mu,:,:)+ffnl1_tmp(ig,1+nu,:,:)*rprimd(mu,nu)
2064 : end do
2065 : end do
2066 : end do
2067 :
2068 : !Second q-derivative
2069 671683104 : ffnl1(:,5:10,:,:)=zero
2070 653184 : do mu=1,6
2071 559872 : mua=alpha(mu);mub=beta(mu)
2072 55958040 : do ig=1,npw1_k
2073 221779296 : do nua=1,3
2074 718963128 : do nub=1,3
2075 497743704 : nu=gamma(nua,nub)
2076 : ffnl1(ig,4+mu,:,:)=ffnl1(ig,4+mu,:,:)+ &
2077 7134326424 : & ffnl1_tmp(ig,4+nu,:,:)*rprimd(mua,nua)*rprimd(mub,nub)
2078 : end do
2079 : end do
2080 : end do
2081 : end do
2082 :
2083 93312 : ABI_FREE(ffnl1_tmp)
2084 : end if
2085 :
2086 : !===== Preparation of the kinetic contributions
2087 : ! Compute (1/2) (2 Pi)**2 (k+q+G)**2:
2088 750504 : ABI_MALLOC(kinpw1,(npw1_k))
2089 25316567 : kinpw1(:)=zero
2090 250168 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
2091 :
2092 750504 : ABI_MALLOC(dqdqkinpw,(npw_k))
2093 : !-- Metric (strain) perturbation
2094 250168 : if (ipert==natom+3.or.ipert==natom+4) then
2095 93312 : call mkkin_metdqdq(dqdqkinpw,dtset%effmass_free,gprimd,idir,kg_k,kpoint,npw_k,qdir1)
2096 : else
2097 16005779 : dqdqkinpw(:)=zero
2098 : end if
2099 :
2100 : !===== Load the k/k+q dependent parts of the Hamiltonian
2101 :
2102 : !Load k-dependent part in the Hamiltonian datastructure
2103 1000672 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk))
2104 : call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
2105 250168 : & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
2106 :
2107 1250840 : if (size(ffnlk)>0) then
2108 156856 : call gs_hamkq%load_k(ffnl_k=ffnlk)
2109 : else
2110 93312 : call gs_hamkq%load_k(ffnl_k=ffnl1)
2111 : end if
2112 :
2113 : !Load k+q-dependent part in the Hamiltonian datastructure
2114 : ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
2115 : call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
2116 : & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,&
2117 250168 : & compute_gbound=.true.)
2118 :
2119 250168 : if (qne0) then
2120 0 : ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk))
2121 0 : call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.)
2122 : end if
2123 :
2124 : !Load k-dependent part in the 1st-order Hamiltonian datastructure
2125 250168 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dqdqkinpw)
2126 :
2127 250168 : end subroutine getgh1dqc_setup
2128 : !!***
2129 :
2130 : !----------------------------------------------------------------------
2131 :
2132 : !!****f* ABINIT/getgh1ndc
2133 : !!
2134 : !! NAME
2135 : !! getgh1ndc
2136 : !!
2137 : !! FUNCTION
2138 : !! Compute 1st order magnetic nuclear dipole moment contribution to <G|H|C>
2139 : !! for input vector |C> expressed in reciprocal space.
2140 : !! Only for DDK perturbation
2141 : !!
2142 : !! INPUTS
2143 : !!
2144 : !! OUTPUT
2145 : !! gh1ndc(2,npw_k*my_nspinor*ndat)=1st order A.p contribution to <G|H|C> for array of nuclear dipoles
2146 : !!
2147 : !! NOTES
2148 : !! This codes only the DDK response for A.p, so effectively A_ipert|C>. The nuclear dipole Hamiltonian
2149 : !! (to first order in the nuclear dipole strength) is A.p where in atomic units
2150 : !! A.p=\alpha^2 m x (r-R)/(r-R)^3 . p. Here the components of A have been precomputed in real space
2151 : !! by make_vectornd. The first-order DDK contribution is i[A.p,r] = A_idir where idir is the
2152 : !! direction of the DDK perturbation, or 2\pi A_idir when A, p, and r are in reduced coords
2153 : !!
2154 : !! SOURCE
2155 :
2156 200556 : subroutine getgh1ndc(cwavein,gh1ndc,gbound_k,istwf_k,kg_k,mgfft,mpi_enreg,&
2157 : & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,&
2158 200556 : & vectornd,vlocal,zora,gpu_option)
2159 :
2160 : !Arguments ------------------------------------
2161 : !scalars
2162 : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,zora
2163 : integer,intent(in) :: gpu_option
2164 : type(MPI_type),intent(in) :: mpi_enreg
2165 : !arrays
2166 : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
2167 : real(dp),intent(inout) :: cwavein(2,npw_k*my_nspinor*ndat)
2168 : real(dp),intent(inout) :: gh1ndc(2,npw_k*my_nspinor*ndat)
2169 : real(dp),intent(inout) :: vectornd(n4,n5,n6,nvloc),vlocal(n4,n5,n6,nvloc)
2170 :
2171 : !Local variables-------------------------------
2172 : !scalars
2173 : integer,parameter :: tim_fourwf=1
2174 : integer :: idat,ipw,iv1,iv2,nspinortot,shift
2175 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc,usezora
2176 : real(dp) :: weight=one
2177 : real(dp), parameter :: HalfFineStruct2=half/InvFineStruct**2
2178 : !arrays
2179 200556 : real(dp),allocatable :: cwavein1(:,:),cwavein2(:,:)
2180 200556 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:)
2181 200556 : real(dp),allocatable :: vectornd_dir(:,:,:,:),work(:,:,:,:),zk(:,:,:)
2182 : ! *********************************************************************
2183 :
2184 129621024 : gh1ndc(:,:)=zero
2185 200556 : if (nvloc/=1) return
2186 :
2187 200556 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
2188 200556 : if (mpi_enreg%paral_spinor==0) then
2189 200556 : shift=npw_k
2190 200556 : nspinor1TreatedByThisProc=.true.
2191 200556 : nspinor2TreatedByThisProc=(nspinortot==2)
2192 : else
2193 0 : shift=0
2194 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
2195 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
2196 : end if
2197 :
2198 200556 : usezora=((zora.EQ.1).OR.(zora.EQ.3))
2199 200556 : if(usezora) then
2200 0 : ABI_MALLOC(zk,(n4,n5,n6))
2201 0 : zk(1:n4,1:1:n5,1:n6)=1.0/(1.0-HalfFineStruct2*vlocal(1:n4,1:n5,1:n6,nvloc))
2202 : end if
2203 :
2204 1002780 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
2205 :
2206 200556 : if (nspinortot==1) then
2207 :
2208 431172 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
2209 : #ifdef HAVE_OPENMP_OFFLOAD
2210 : !$OMP TARGET ENTER DATA MAP(alloc:ghc1) IF(gpu_option==ABI_GPU_OPENMP)
2211 : #endif
2212 :
2213 862344 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
2214 143724 : if (usezora) then
2215 0 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=zk(1:n4,1:n5,1:n6)*vectornd(1:n4,1:n5,1:n6,nvloc)
2216 : else
2217 809155404 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=vectornd(1:n4,1:n5,1:n6,nvloc)
2218 : end if
2219 : ! apply vector potential in direction ipert to input wavefunction
2220 : call fourwf(1,vectornd_dir,cwavein,ghc1,work,gbound_k,gbound_k,&
2221 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2222 143724 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2223 143724 : ABI_FREE(vectornd_dir)
2224 :
2225 : ! scale by 2\pi
2226 143724 : if(gpu_option==ABI_GPU_DISABLED) then
2227 79694112 : gh1ndc=two_pi*ghc1
2228 : else if(gpu_option==ABI_GPU_OPENMP) then
2229 : #ifdef HAVE_OPENMP_OFFLOAD
2230 : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO MAP(to:gh1ndc,ghc1) PRIVATE(ipw)
2231 : do ipw=1,npw_k*ndat
2232 : gh1ndc(1,ipw)=two_pi*ghc1(1,ipw)
2233 : gh1ndc(2,ipw)=two_pi*ghc1(2,ipw)
2234 : end do
2235 : #endif
2236 : end if
2237 :
2238 : #ifdef HAVE_OPENMP_OFFLOAD
2239 : !$OMP TARGET EXIT DATA MAP(delete:ghc1) IF(gpu_option==ABI_GPU_OPENMP)
2240 : #endif
2241 143724 : ABI_FREE(ghc1)
2242 :
2243 : else ! nspinortot==2
2244 :
2245 56832 : if (nspinor1TreatedByThisProc) then
2246 :
2247 170496 : ABI_MALLOC(cwavein1,(2,npw_k*ndat))
2248 113664 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
2249 : #ifdef HAVE_OPENMP_OFFLOAD
2250 : !$OMP TARGET ENTER DATA MAP(alloc:ghc1,cwavein1) IF(gpu_option==ABI_GPU_OPENMP)
2251 : #endif
2252 :
2253 56832 : if(gpu_option==ABI_GPU_DISABLED) then
2254 113664 : do idat=1,ndat
2255 8425344 : do ipw=1,npw_k
2256 24991872 : cwavein1(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k)
2257 : end do
2258 : end do
2259 : else if(gpu_option==ABI_GPU_OPENMP) then
2260 : #ifdef HAVE_OPENMP_OFFLOAD
2261 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavein1,cwavein) PRIVATE(idat)
2262 : do idat=1,ndat
2263 : !$OMP PARALLEL DO PRIVATE(ipw)
2264 : do ipw=1,npw_k
2265 : cwavein1(1,ipw+(idat-1)*npw_k)=cwavein(1,ipw+(idat-1)*my_nspinor*npw_k)
2266 : cwavein1(2,ipw+(idat-1)*npw_k)=cwavein(2,ipw+(idat-1)*my_nspinor*npw_k)
2267 : end do
2268 : end do
2269 : #endif
2270 : end if
2271 :
2272 340992 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
2273 56832 : if (usezora) then
2274 0 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=zk(1:n4,1:n5,1:n6)*vectornd(1:n4,1:n5,1:n6,nvloc)
2275 : else
2276 279215616 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=vectornd(1:n4,1:n5,1:n6,nvloc)
2277 : end if
2278 : call fourwf(1,vectornd_dir,cwavein1,ghc1,work,gbound_k,gbound_k,&
2279 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2280 56832 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2281 56832 : ABI_FREE(vectornd_dir)
2282 :
2283 56832 : if(gpu_option==ABI_GPU_DISABLED) then
2284 113664 : do idat=1,ndat
2285 56832 : iv1=1+(idat-1)*npw_k; iv2=npw_k+(idat-1)*npw_k
2286 25048704 : gh1ndc(1:2,iv1:iv2)=two_pi*ghc1(1:2,iv1:iv2)
2287 : end do
2288 : else if(gpu_option==ABI_GPU_OPENMP) then
2289 : #ifdef HAVE_OPENMP_OFFLOAD
2290 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:gh1ndc,ghc1) PRIVATE(idat)
2291 : do idat=1,ndat
2292 : !$OMP PARALLEL DO PRIVATE(ipw)
2293 : do ipw=1,npw_k
2294 : gh1ndc(1,ipw+(idat-1)*npw_k)=two_pi*ghc1(1,ipw+(idat-1)*npw_k)
2295 : gh1ndc(2,ipw+(idat-1)*npw_k)=two_pi*ghc1(2,ipw+(idat-1)*npw_k)
2296 : end do
2297 : end do
2298 : #endif
2299 : end if
2300 :
2301 : #ifdef HAVE_OPENMP_OFFLOAD
2302 : !$OMP TARGET EXIT DATA MAP(delete:ghc1,cwavein1) IF(gpu_option==ABI_GPU_OPENMP)
2303 : #endif
2304 56832 : ABI_FREE(ghc1)
2305 56832 : ABI_FREE(cwavein1)
2306 :
2307 : end if ! end spinor 1
2308 :
2309 56832 : if (nspinor2TreatedByThisProc) then
2310 :
2311 170496 : ABI_MALLOC(cwavein2,(2,npw_k*ndat))
2312 113664 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
2313 : #ifdef HAVE_OPENMP_OFFLOAD
2314 : !$OMP TARGET ENTER DATA MAP(alloc:ghc2,cwavein2) IF(gpu_option==ABI_GPU_OPENMP)
2315 : #endif
2316 :
2317 56832 : if(gpu_option==ABI_GPU_DISABLED) then
2318 113664 : do idat=1,ndat
2319 8425344 : do ipw=1,npw_k
2320 24991872 : cwavein2(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift)
2321 : end do
2322 : end do
2323 : else if(gpu_option==ABI_GPU_OPENMP) then
2324 : #ifdef HAVE_OPENMP_OFFLOAD
2325 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavein2,cwavein) PRIVATE(idat)
2326 : do idat=1,ndat
2327 : !$OMP PARALLEL DO PRIVATE(ipw)
2328 : do ipw=1,npw_k
2329 : cwavein2(1,ipw+(idat-1)*npw_k)=cwavein(1,ipw+(idat-1)*my_nspinor*npw_k+shift)
2330 : cwavein2(2,ipw+(idat-1)*npw_k)=cwavein(2,ipw+(idat-1)*my_nspinor*npw_k+shift)
2331 : end do
2332 : end do
2333 : #endif
2334 : end if
2335 :
2336 340992 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
2337 56832 : if (usezora) then
2338 0 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=zk(1:n4,1:n5,1:n6)*vectornd(1:n4,1:n5,1:n6,nvloc)
2339 : else
2340 279215616 : vectornd_dir(1:n4,1:n5,1:n6,nvloc)=vectornd(1:n4,1:n5,1:n6,nvloc)
2341 : end if
2342 : call fourwf(1,vectornd_dir,cwavein2,ghc2,work,gbound_k,gbound_k,&
2343 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2344 56832 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2345 56832 : ABI_FREE(vectornd_dir)
2346 :
2347 56832 : if(gpu_option==ABI_GPU_DISABLED) then
2348 113664 : do idat=1,ndat
2349 56832 : iv1=1+(idat-1)*npw_k; iv2=npw_k+(idat-1)*npw_k
2350 25048704 : gh1ndc(1:2,iv1+shift:iv2+shift)=two_pi*ghc2(1:2,iv1:iv2)
2351 : end do
2352 : else if(gpu_option==ABI_GPU_OPENMP) then
2353 : #ifdef HAVE_OPENMP_OFFLOAD
2354 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:gh1ndc,ghc2) PRIVATE(idat)
2355 : do idat=1,ndat
2356 : !$OMP PARALLEL DO PRIVATE(ipw)
2357 : do ipw=1,npw_k
2358 : gh1ndc(1,ipw+(idat-1)*npw_k+shift)=two_pi*ghc2(1,ipw+(idat-1)*npw_k)
2359 : gh1ndc(2,ipw+(idat-1)*npw_k+shift)=two_pi*ghc2(2,ipw+(idat-1)*npw_k)
2360 : end do
2361 : end do
2362 : #endif
2363 : end if
2364 :
2365 : #ifdef HAVE_OPENMP_OFFLOAD
2366 : !$OMP TARGET EXIT DATA MAP(delete:ghc2,cwavein2) IF(gpu_option==ABI_GPU_OPENMP)
2367 : #endif
2368 56832 : ABI_FREE(ghc2)
2369 56832 : ABI_FREE(cwavein2)
2370 :
2371 : end if ! end spinor 2
2372 :
2373 : end if ! nspinortot
2374 :
2375 200556 : ABI_FREE(work)
2376 200556 : if (usezora) then
2377 0 : ABI_FREE(zk)
2378 : end if
2379 :
2380 200556 : end subroutine getgh1ndc
2381 : !!***
2382 :
2383 : !!****f* ABINIT/getgh1c_mGGA
2384 : !!
2385 : !! NAME
2386 : !! getgh1c_mGGA
2387 : !!
2388 : !! FUNCTION
2389 : !! Compute first order metaGGA contribution to <G|H|C> for input vector |C> expressed in reciprocal space.
2390 : !! ONLY FOR DDK PERTURBATION
2391 : !!
2392 : !! INPUTS
2393 : !!
2394 : !! OUTPUT
2395 : !! gh1c_mGGA(2,npw_k*my_nspinor*ndat)=metaGGA contribution to <G|H1|C>
2396 : !!
2397 : !! SIDE EFFECTS
2398 : !!
2399 : !! SOURCE
2400 :
2401 19584 : subroutine getgh1c_mGGA(cwavein,gbound_k,gh1c_mGGA,gmet,gprimd,idir,istwf_k,kg_k,&
2402 19584 : & kpt,mgfft,mpi_enreg,my_nspinor,n4,n5,n6,ndat,ngfft,npw_k,nvloc,vxctaulocal,gpu_option)
2403 :
2404 : !Arguments ------------------------------------
2405 : !scalars
2406 : integer,intent(in) :: idir,istwf_k,mgfft,my_nspinor,n4,n5,n6,ndat,npw_k,nvloc
2407 : integer,intent(in),optional :: gpu_option
2408 : type(MPI_type),intent(in) :: mpi_enreg
2409 : !arrays
2410 : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
2411 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3)
2412 : real(dp),intent(inout) :: cwavein(2,npw_k*my_nspinor*ndat)
2413 : real(dp),intent(inout) :: gh1c_mGGA(2,npw_k*my_nspinor*ndat)
2414 : real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvloc,4)
2415 :
2416 : !Local variables-------------------------------
2417 : !scalars
2418 : integer :: idat,ii,ipw,nspinortot,shift,gpu_option_
2419 : integer,parameter :: tim_fourwf=1
2420 : real(dp) :: weight=one
2421 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
2422 : !arrays
2423 19584 : real(dp),allocatable :: cwavein1(:,:),cwavein2(:,:),dgcwavef(:,:,:),dlcwavef(:,:)
2424 19584 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:),work(:,:,:,:)
2425 :
2426 : if(present(gpu_option)) then
2427 : gpu_option_=gpu_option
2428 : else
2429 : gpu_option_=0
2430 : end if
2431 :
2432 14340384 : gh1c_mGGA(:,:)=zero
2433 :
2434 19584 : if (nvloc/=1) return
2435 :
2436 19584 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
2437 19584 : if (mpi_enreg%paral_spinor==0) then
2438 19584 : shift=npw_k
2439 19584 : nspinor1TreatedByThisProc=.true.
2440 19584 : nspinor2TreatedByThisProc=(nspinortot==2)
2441 : else
2442 0 : shift=0
2443 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
2444 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
2445 : end if
2446 :
2447 97920 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
2448 :
2449 19584 : if (nspinortot==1) then
2450 :
2451 19584 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
2452 26112 : ABI_MALLOC(dgcwavef,(2,npw_k*ndat,3))
2453 13056 : ABI_MALLOC(dlcwavef,(2,npw_k*ndat))
2454 :
2455 : !!$OMP PARALLEL DO
2456 11482752 : dgcwavef = zero; dlcwavef = zero
2457 13056 : do idat=1,ndat
2458 967776 : do ipw=1,npw_k
2459 3818880 : dgcwavef(1,ipw+(idat-1)*npw_k,1:3)=+half*gprimd(idir,1:3)*cwavein(2,ipw+(idat-1)*npw_k)
2460 3818880 : dgcwavef(2,ipw+(idat-1)*npw_k,1:3)=-half*gprimd(idir,1:3)*cwavein(1,ipw+(idat-1)*npw_k)
2461 : dlcwavef(1:2,ipw+(idat-1)*npw_k)=two_pi*two_pi*DOT_PRODUCT(gmet(idir,:),(kpt(:)+kg_k(:,ipw)))*&
2462 5734848 : &cwavein(1:2,ipw+(idat-1)*npw_k)
2463 : end do
2464 : end do
2465 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
2466 2870688 : ghc1=zero
2467 : call fourwf(1,vxctaulocal(:,:,:,:,1),dlcwavef,ghc1,work,gbound_k,gbound_k,&
2468 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2469 6528 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2470 : !!$OMP PARALLEL DO
2471 13056 : do idat=1,ndat
2472 967776 : do ipw=1,npw_k
2473 2870688 : gh1c_mGGA(:,ipw+(idat-1)*npw_k)=gh1c_mGGA(:,ipw+(idat-1)*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
2474 : end do
2475 : end do
2476 : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
2477 26112 : do ii=1,3
2478 8612064 : ghc1=zero
2479 : call fourwf(1,vxctaulocal(:,:,:,:,1+ii),dgcwavef(:,:,ii),ghc1,work,gbound_k,gbound_k,&
2480 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2481 19584 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2482 : !!$OMP PARALLEL DO
2483 45696 : do idat=1,ndat
2484 2903328 : do ipw=1,npw_k
2485 8612064 : gh1c_mGGA(:,ipw+(idat-1)*npw_k)=gh1c_mGGA(:,ipw+(idat-1)*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
2486 : end do
2487 : end do
2488 : end do ! ii
2489 6528 : ABI_FREE(dgcwavef)
2490 6528 : ABI_FREE(dlcwavef)
2491 6528 : ABI_FREE(ghc1)
2492 :
2493 : else ! nspinortot==2
2494 :
2495 39168 : ABI_MALLOC(cwavein1,(2,npw_k*ndat))
2496 26112 : ABI_MALLOC(cwavein2,(2,npw_k*ndat))
2497 26112 : do idat=1,ndat
2498 1935552 : do ipw=1,npw_k
2499 5728320 : cwavein1(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k)
2500 5741376 : cwavein2(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift)
2501 : end do
2502 : end do
2503 :
2504 13056 : if (nspinor1TreatedByThisProc) then
2505 26112 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
2506 52224 : ABI_MALLOC(dgcwavef,(2,npw_k*ndat,3))
2507 26112 : ABI_MALLOC(dlcwavef,(2,npw_k*ndat))
2508 :
2509 22965504 : dgcwavef = zero; dlcwavef = zero
2510 26112 : do idat=1,ndat
2511 1935552 : do ipw=1,npw_k
2512 7637760 : dgcwavef(1,ipw+(idat-1)*npw_k,1:3)=+half*gprimd(idir,1:3)*cwavein1(2,ipw+(idat-1)*npw_k)
2513 7637760 : dgcwavef(2,ipw+(idat-1)*npw_k,1:3)=-half*gprimd(idir,1:3)*cwavein1(1,ipw+(idat-1)*npw_k)
2514 : dlcwavef(1:2,ipw+(idat-1)*npw_k)=two_pi*two_pi*DOT_PRODUCT(gmet(idir,:),(kpt(:)+kg_k(:,ipw)))*&
2515 11469696 : &cwavein1(1:2,ipw+(idat-1)*npw_k)
2516 : end do
2517 : end do
2518 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
2519 5741376 : ghc1=zero
2520 : call fourwf(1,vxctaulocal(:,:,:,:,1),dlcwavef,ghc1,work,gbound_k,gbound_k,&
2521 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2522 13056 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2523 26112 : do idat=1,ndat
2524 1935552 : do ipw=1,npw_k
2525 5741376 : gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
2526 : end do
2527 : end do
2528 : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
2529 52224 : do ii=1,3
2530 17224128 : ghc1=zero
2531 : call fourwf(1,vxctaulocal(:,:,:,:,1+ii),dgcwavef(:,:,ii),ghc1,work,gbound_k,gbound_k,&
2532 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2533 39168 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2534 91392 : do idat=1,ndat
2535 5806656 : do ipw=1,npw_k
2536 17224128 : gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
2537 : end do
2538 : end do
2539 : end do ! ii
2540 13056 : ABI_FREE(dgcwavef)
2541 13056 : ABI_FREE(dlcwavef)
2542 :
2543 13056 : ABI_FREE(ghc1)
2544 : end if ! end spinor 1
2545 :
2546 13056 : if (nspinor2TreatedByThisProc) then
2547 26112 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
2548 52224 : ABI_MALLOC(dgcwavef,(2,npw_k*ndat,3))
2549 26112 : ABI_MALLOC(dlcwavef,(2,npw_k*ndat))
2550 :
2551 22965504 : dgcwavef = zero; dlcwavef = zero
2552 26112 : do idat=1,ndat
2553 1935552 : do ipw=1,npw_k
2554 7637760 : dgcwavef(1,ipw+(idat-1)*npw_k,1:3)=+half*gprimd(idir,1:3)*cwavein2(2,ipw+(idat-1)*npw_k)
2555 7637760 : dgcwavef(2,ipw+(idat-1)*npw_k,1:3)=-half*gprimd(idir,1:3)*cwavein2(1,ipw+(idat-1)*npw_k)
2556 : dlcwavef(1:2,ipw+(idat-1)*npw_k)=two_pi*two_pi*DOT_PRODUCT(gmet(idir,:),(kpt(:)+kg_k(:,ipw)))*&
2557 11469696 : &cwavein2(1:2,ipw+(idat-1)*npw_k)
2558 : end do
2559 : end do
2560 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
2561 5741376 : ghc2=zero
2562 : call fourwf(1,vxctaulocal(:,:,:,:,1),dlcwavef,ghc2,work,gbound_k,gbound_k,&
2563 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2564 13056 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2565 26112 : do idat=1,ndat
2566 1935552 : do ipw=1,npw_k
2567 : gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)+&
2568 5741376 : & ghc2(:,ipw+(idat-1)*npw_k)
2569 : end do
2570 : end do
2571 : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
2572 52224 : do ii=1,3
2573 17224128 : ghc2=zero
2574 : call fourwf(1,vxctaulocal(:,:,:,:,1+ii),dgcwavef(:,:,ii),ghc2,work,gbound_k,gbound_k,&
2575 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
2576 39168 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2577 91392 : do idat=1,ndat
2578 5806656 : do ipw=1,npw_k
2579 : gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=gh1c_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)+&
2580 17224128 : & ghc2(:,ipw+(idat-1)*npw_k)
2581 : end do
2582 : end do
2583 : end do ! ii
2584 13056 : ABI_FREE(dgcwavef)
2585 13056 : ABI_FREE(dlcwavef)
2586 :
2587 13056 : ABI_FREE(ghc2)
2588 : end if ! end spinor 2
2589 :
2590 13056 : ABI_FREE(cwavein1)
2591 13056 : ABI_FREE(cwavein2)
2592 :
2593 : end if ! nspinortot
2594 :
2595 19584 : ABI_FREE(work)
2596 :
2597 : end subroutine getgh1c_mGGA
2598 : !!***
2599 :
2600 : end module m_getgh1c
2601 : !!***
|
