Branch data 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 [ - + ]: 8601579 : subroutine getgh1c(berryopt,cwave,cwaveprj,gh1c,grad_berry,gs1c,gs_hamkq,&
133 [ - + ]: 8601579 : 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 : 17203158 : real(dp) :: enlout(ndat),tsec(2),svectout_dum(1,1),vectout_dum(1,1)
165 : 8601579 : real(dp),allocatable :: cwave_sp(:,:),cwavef1(:,:),cwavef2(:,:)
166 : 8601579 : real(dp),allocatable :: gh1c_sp(:,:),gh1c1(:,:),gh1c2(:,:),gh1c3(:,:),gh1c4(:,:)
167 : 8601579 : real(dp),allocatable :: gh1c_mGGA(:,:),gh1ndc(:,:),gvnl2(:,:)
168 : 8601579 : real(dp),target,allocatable :: nonlop_out(:,:),vlocal1_tmp(:,:,:)
169 : : !#define _DEV_USE_WORK
170 : : !#ifdef _DEV_USE_WORK
171 : 8601579 : real(dp),target,allocatable :: work(:,:,:,:)
172 : : !#endif
173 : 8601579 : real(dp),contiguous, pointer :: gvnlx1_(:,:), dkinpw(:),kinpw1(:)
174 : 8601579 : type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:)
175 : 8601579 : 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 : 8601579 : call timab(196+tim_getgh1c,1,tsec)
185 : : ABI_NVTX_START_RANGE(NVTX_GETGH1C)
186 : :
187 : : !======================================================================
188 : : !== Initialisations and compatibility tests
189 : : !======================================================================
190 : :
191 : 8601579 : npw = gs_hamkq%npw_k
192 : 8601579 : npw1 = gs_hamkq%npw_kp
193 : 8601579 : natom = gs_hamkq%natom
194 : :
195 : : ! Compatibility tests
196 [ + + + - : 8601579 : 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 [ + + + + ]: 8601579 : 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 [ + + - + ]: 8601579 : 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 [ - + ]: 8601579 : 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 : 8601579 : my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
220 [ + + - + : 25804737 : ABI_CHECK_IGEQ(size(cwave), 2*npw*my_nspinor*ndat, 'wrong size for cwave!')
- - ]
221 [ + + - + : 25804737 : ABI_CHECK_IGEQ(size(gh1c), 2*npw1*my_nspinor*ndat, 'wrong size for gh1c!')
- - ]
222 [ + + ]: 8601579 : if (usevnl/=0) then
223 [ + + - + : 13997439 : ABI_CHECK_IGEQ(size(gvnlx1), 2*npw1*my_nspinor*ndat, 'wrong size for gvnlx1!')
- - ]
224 : : end if
225 [ + + ]: 8601579 : if (sij_opt==1) then
226 [ + + - + : 3019632 : ABI_CHECK_IGEQ(size(gs1c), 2*npw1*my_nspinor*ndat, 'wrong size for gs1c!')
- - ]
227 : : end if
228 [ + + ]: 8601579 : 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 : 8601579 : usecprj=gs_hamkq%usecprj ; ncpgr=0
234 [ + + ]: 8601579 : if(gs_hamkq%usepaw==1) then
235 [ + + + + ]: 3110652 : if (size(cwaveprj)==0) usecprj=0
236 [ - + ]: 1008084 : if (usecprj/=0) then
237 : 1008084 : ncpgr=cwaveprj(1,1)%ncpgr
238 [ + + - + : 3024252 : ABI_CHECK_IGEQ(size(cwaveprj), gs_hamkq%natom*my_nspinor*ndat, 'wrong size for cwaveprj!')
- - ]
239 [ + - - + : 1008084 : 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 : 8601579 : tim_nonlop=8
252 [ + + + + ]: 8601579 : if (tim_getgh1c==1.and.ipert<=natom) tim_nonlop=7
253 [ + + + + ]: 8601579 : if (tim_getgh1c==2.and.ipert<=natom) tim_nonlop=5
254 [ + + + + ]: 8601579 : if (tim_getgh1c==1.and.ipert> natom) tim_nonlop=8
255 [ + + + + ]: 8601579 : if (tim_getgh1c==2.and.ipert> natom) tim_nonlop=5
256 [ + + ]: 8601579 : if (tim_getgh1c==3 ) tim_nonlop=0
257 : :
258 : 8601579 : compute_conjugate = .false.
259 [ + + ]: 8601579 : 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 : 8601579 : 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 [ + + + - : 8601579 : 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 [ + + + - ]: 4342052 : 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 [ + - ]: 4342052 : if(gs_hamkq%gpu_option/=ABI_GPU_OPENMP) then
454 : : !$OMP PARALLEL DO
455 [ + + ]: 982625138 : do ipw=1,npw1*my_nspinor*ndat
456 [ + + ]: 2939191310 : 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 [ + + ]: 8601579 : if (usevnl==1) then
472 : 4665813 : gvnlx1_ => gvnlx1
473 : : else
474 [ - + + - : 11807298 : 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 [ + + - + : 8601579 : 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 [ + + + + : 2368186 : else if (ipert==natom+1.and.(optnl>0.or.sij_opt/=0)) then
+ - ]
572 : :
573 : : ABI_NVTX_START_RANGE(NVTX_GETGH1_DDK)
574 : :
575 : 1437654 : tim_nonlop=8 ; signs=2 ; choice=5
576 [ + + ]: 1437654 : if (gs_hamkq%usepaw==1) then
577 [ + + ]: 258940 : if (usecprj==1) then
578 : 255484 : 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 [ + - ]: 258940 : 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 [ + + ]: 258940 : if(gs_hamkq%nspinor==2) cpopt=-1
588 [ + + ]: 258940 : if(associated(gs_hamkq%vectornd)) cpopt=-1
589 : : call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,lambda,mpi_enreg,ndat,nnlout,&
590 : 258940 : paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_)
591 [ + + ]: 258940 : if (usecprj==0) then
592 : 3456 : call pawcprj_free(cwaveprj_tmp)
593 [ + - + + : 10368 : ABI_FREE(cwaveprj_tmp)
- + - + -
+ ]
594 : : end if
595 : 258940 : 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 : 8601579 : usevnl2=allocated(gvnl2)
888 [ + + + + : 8601579 : has_kin=(ipert==natom+1.or.ipert==natom+3.or.ipert==natom+4)
+ + ]
889 [ + - ]: 8601579 : if (associated(gs_hamkq%kinpw_kp)) then
890 : 8601579 : 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 [ + + ]: 8601579 : if (associated(rf_hamkq%dkinpw_k)) then
895 : 5622660 : dkinpw => rf_hamkq%dkinpw_k
896 [ - + ]: 2978919 : else if (has_kin) then
897 : 0 : ABI_BUG('need dkinpw allocated!')
898 : : end if
899 : :
900 [ + + ]: 5622660 : 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 [ + + ]: 3677390 : do ispinor=1,my_nspinor*ndat
907 : : !$OMP PARALLEL DO PRIVATE(ipw,ipws)
908 [ + + ]: 273428640 : do ipw=1,npw
909 : 269751250 : ipws=ipw+npw*(ispinor-1)
910 [ + - ]: 271619033 : if(kinpw1(ipw)<threshold)then
911 : 269751250 : gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+dkinpw(ipw)*cwave(1,ipws)
912 : 269751250 : 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 [ + + + + ]: 8601579 : has_nd1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vectornd) )
928 : :
929 : : if (has_nd1) then
930 [ - + + - : 379188 : 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 [ + - - - : 126396 : & 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 [ + + ]: 285816 : do ispinor=1,my_nspinor*ndat
944 : : !$OMP PARALLEL DO PRIVATE(ipw,ipws)
945 [ + + ]: 29235012 : do ipw=1,npw
946 : 28949196 : ipws=ipw+npw*(ispinor-1)
947 : 28949196 : gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+gh1ndc(1,ipws)
948 : 29108616 : 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 [ - + ]: 126396 : 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 [ + + + + ]: 8601579 : 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 [ + + + + ]: 8601579 : 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 [ + + ]: 17813196 : do ispinor=1,my_nspinor*ndat
1002 : 9211801 : ipws=(ispinor-1)*npw1
1003 : : !$OMP PARALLEL DO PRIVATE(ipw)
1004 [ + + ]: 2538139373 : do ipw=1+ipws,npw1+ipws
1005 [ + - ]: 2529537794 : if(kinpw1(ipw-ipws)<threshold)then
1006 : 2520325993 : gh1c(1,ipw)=gh1c(1,ipw)+gvnlx1_(1,ipw)
1007 : 2520325993 : 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 [ + + ]: 8601579 : 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 [ + + ]: 8601579 : if (usevnl==1) then
1048 : 4665813 : nullify(gvnlx1_)
1049 : : else
1050 [ - + ]: 3935766 : ABI_FREE(gvnlx1_)
1051 : : end if
1052 : :
1053 : : ABI_NVTX_END_RANGE()
1054 : :
1055 : 8601579 : 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 [ - + - + : 17203158 : 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 : 386001 : subroutine rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,nfft,ngfft,nvloc,&
1118 : 386001 : & 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 : 386001 : real(dp), contiguous, pointer :: vtrial_ptr(:,:),vtrial1_ptr(:,:)
1138 : 386001 : real(dp),allocatable :: cgrvtrial(:,:),cgrvtrial1(:,:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:)
1139 : : ! *************************************************************************
1140 : :
1141 : 386001 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
1142 : 386001 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
1143 : 386001 : paral_kgb = mpi_enreg%paral_kgb
1144 : :
1145 : 386001 : if (nspden/=4) then
1146 : 385706 : vtrial_ptr => vtrial
1147 [ + + + + ]: 385706 : 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 [ - + + - : 12140 : ABI_MALLOC(cgrvtrial,(nfft,nspden))
- + + - +
- + - - +
- + ]
1152 : 3035 : call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr)
1153 : 3035 : call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal(:,:,:,1),2)
1154 [ + - - + : 12140 : ABI_REMALLOC(cgrvtrial, (cplex*nfft, nspden))
+ - - + +
- + - + -
- + - + -
+ ]
1155 : 3035 : call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial1)
1156 : 3035 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,cgrvtrial,vlocal1(:,:,:,1),2)
1157 [ - + ]: 3035 : ABI_FREE(cgrvtrial)
1158 : : end if
1159 : 386001 : 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 [ + + - + ]: 772002 : 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 : 916900 : 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 : 916900 : 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 : 916900 : call timab(195, 1, tsec)
1251 : : !call cwtime(cpu, wall, gflops, "start")
1252 : :
1253 [ + + ]: 916900 : reuse_ffnlk_ = 0; if (present(reuse_ffnlk)) reuse_ffnlk_ = reuse_ffnlk
1254 [ + + ]: 916900 : reuse_ffnl1_ = 0; if (present(reuse_ffnl1)) reuse_ffnl1_ = reuse_ffnl1
1255 [ - + ]: 916900 : reuse_kpg_k_ = 0; if (present(reuse_kpg_k)) reuse_kpg_k_ = reuse_kpg_k
1256 [ - + ]: 916900 : reuse_kpg1_k_ = 0; if (present(reuse_kpg1_k)) reuse_kpg1_k_ = reuse_kpg1_k
1257 : :
1258 [ + + - + ]: 916900 : 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 [ + + - + : 916900 : 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 [ + + - + : 916900 : 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 : 916900 : ntypat = psps%ntypat
1269 : 916900 : qne0 = ((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14)
1270 : :
1271 : : ! Compute k+G vectors
1272 [ + - + + ]: 916900 : nkpg = 0; if (ipert >= 1 .and. ipert <= natom) nkpg = 3*dtset%nloalg(3)
1273 [ + - ]: 916900 : if (reuse_kpg_k_ == 0) then
1274 [ - + + - : 3667600 : ABI_MALLOC(kpg_k, (npw_k, nkpg))
+ + + - +
- + + - +
- + - + ]
1275 [ + + ]: 916900 : 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 [ + - + + ]: 916900 : nkpg1 = 0; if (ipert >= 1 .and. ipert <= natom) nkpg1 = 3*dtset%nloalg(3)
1282 [ + - ]: 916900 : if (reuse_kpg1_k_ == 0) then
1283 [ - + + - : 3667600 : ABI_MALLOC(kpg1_k, (npw1_k, nkpg1))
+ + + - +
- + + - +
- + - + ]
1284 [ + + ]: 916900 : 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 [ + + ]: 916900 : dimffnlk =0; if (ipert<=natom) dimffnlk=1
1291 : :
1292 : : ! Compute nonlocal form factors ffnlk at (k+G)
1293 : : ! (only for atomic displacement perturbation)
1294 [ + + ]: 916900 : if (reuse_ffnlk_ == 0) then
1295 [ - + + - : 5210640 : ABI_MALLOC(ffnlk, (npw_k, dimffnlk, psps%lmnmax, ntypat))
+ + + - -
+ + - - +
+ - + - +
+ - + - +
- + ]
1296 [ + + ]: 868440 : 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 [ + + ]: 916900 : if (ipert<=natom) then
1310 : 510820 : ider=0;idir0=0
1311 : : !-- k-point perturbation (1st-derivative)
1312 [ + + ]: 406080 : else if (ipert==natom+1) then
1313 : 223094 : 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 : 916900 : dimffnl1=1+ider
1336 [ + + - + ]: 916900 : if (ider==1.and.idir0==0) dimffnl1=2+2*psps%useylm
1337 [ + + + - ]: 916900 : if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm
1338 : :
1339 [ + + ]: 916900 : if (reuse_ffnl1_ == 0) then
1340 [ - + + - : 5210640 : 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 : 868440 : 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 : 916900 : print_info = 0
1351 [ + - + - : 916900 : 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 [ + + + + : 916900 : 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 [ - + + - : 2750700 : ABI_MALLOC(kinpw1, (npw1_k))
+ - + - -
+ - + -
+ ]
1367 [ + + ]: 127975165 : kinpw1(:)=zero
1368 : 916900 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
1369 : :
1370 [ - + + - : 2750700 : ABI_MALLOC(dkinpw,(npw_k)) ! 1st derivative (1st direction)
+ - + - -
+ - + -
+ ]
1371 [ + + ]: 127967009 : dkinpw(:)=zero
1372 [ + + + + ]: 916900 : 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 [ + + ]: 916900 : 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 [ + + ]: 916900 : if (ipert==natom+1) then
1383 : : ! Compute the derivative of the kinetic operator vs k
1384 : 223094 : 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 [ + + + + ]: 916900 : 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 [ + + + + ]: 916900 : 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 [ - + + - : 3667600 : 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 : 916900 : ph3d_k=ph3d, compute_ph3d=.true., compute_gbound=.true.)
1447 : :
1448 [ + + + + ]: 4584500 : if (size(ffnlk)>0) then
1449 : 510820 : call gs_hamkq%load_k(ffnl_k=ffnlk)
1450 : : else
1451 : 406080 : 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 : 916900 : kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,compute_gbound=.true.)
1458 : :
1459 [ + + ]: 916900 : 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 : 916900 : call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw)
1466 : :
1467 [ + + ]: 916900 : 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 : 916900 : call timab(195, 2, tsec)
1473 : :
1474 : 916900 : 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 : 816160 : subroutine getdc1(band,band_procs,bands_treated_now,cgq,cprjq,dcwavef,dcwaveprj,&
1515 : : & ibgq,icgq,istwfk,mcgq,mcprjq,&
1516 : 816160 : & 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 : 816160 : integer, allocatable :: nlmn(:)
1541 : 816160 : real(dp), allocatable :: dummy(:,:),scprod(:,:,:)
1542 : : real(dp), allocatable :: dcwavef_tmp(:,:)
1543 : 816160 : type(pawcprj_type),allocatable :: dcwaveprj_tmp(:,:)
1544 : : ! *********************************************************************
1545 : :
1546 : : DBG_ENTER("COLL")
1547 : :
1548 [ + - ]: 816160 : l_gpu_option = ABI_GPU_DISABLED; if (present(gpu_option)) l_gpu_option = gpu_option
1549 : 816160 : nproc_band=mpi_enreg%nproc_band
1550 : :
1551 [ - + ]: 816160 : ABI_MALLOC(dummy,(0,0))
1552 [ - + + - : 3264640 : ABI_MALLOC(scprod,(2,nband_me,ndat))
- + + - +
- + - - +
- + ]
1553 [ - + + - : 2448480 : 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 [ + + ]: 816160 : 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 [ - + - - : 6301456 : do band_ = 1, nband, ndat
- + + + ]
1572 [ + + ]: 5485296 : if (bands_treated_now(band_, 1) == 0) cycle
1573 [ + - ]: 825376 : if(l_gpu_option==ABI_GPU_DISABLED) then
1574 [ + + + + ]: 430797214 : 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 [ + + ]: 825376 : if (band_ == band) then
1582 [ + - ]: 816160 : if(l_gpu_option==ABI_GPU_DISABLED) then
1583 : : !$OMP PARALLEL DO
1584 [ + + ]: 125889402 : do ipw=1,npw1*nspinor*ndat
1585 [ + + ]: 376035886 : 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 : 825376 : 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 [ + + ]: 1650752 : 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 : 1650752 : & 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 : 825376 : 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 [ + + ]: 825376 : 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 [ + - ]: 816160 : if(l_gpu_option==ABI_GPU_DISABLED) then
1624 : : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(dcwavef,s1cwave0,dcwavef_tmp,npw1,nspinor)
1625 [ + + ]: 125889402 : do ipw=1,npw1*nspinor*ndat
1626 [ + + ]: 376035886 : 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 [ + + + - ]: 2466912 : 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 [ - + ]: 816160 : ABI_FREE(dummy)
1668 [ - + ]: 816160 : ABI_FREE(scprod)
1669 [ - + ]: 816160 : ABI_FREE(dcwavef_tmp)
1670 [ + + ]: 816160 : 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 [ - + - + : 816160 : 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 : 200736 : 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 : 200736 : & 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 : 200736 : real(dp),allocatable :: cwavein1(:,:),cwavein2(:,:)
2180 : 200736 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:)
2181 : 200736 : real(dp),allocatable :: vectornd_dir(:,:,:,:),work(:,:,:,:),zk(:,:,:)
2182 : : ! *********************************************************************
2183 : :
2184 [ + + + + ]: 141565464 : gh1ndc(:,:)=zero
2185 [ + - ]: 200736 : if (nvloc/=1) return
2186 : :
2187 : 200736 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
2188 [ + - ]: 200736 : if (mpi_enreg%paral_spinor==0) then
2189 : 200736 : shift=npw_k
2190 : 200736 : nspinor1TreatedByThisProc=.true.
2191 : 200736 : 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 : 200736 : usezora=((zora.EQ.1).OR.(zora.EQ.3))
2199 [ - + ]: 200736 : 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 [ - + + - : 1003680 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
2205 : :
2206 [ + + ]: 200736 : if (nspinortot==1) then
2207 : :
2208 [ - + + - : 431712 : 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 [ - + + - : 863424 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
2214 [ - + ]: 143904 : 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 [ + + + + : 879178464 : 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 : 143904 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
2223 [ - + ]: 143904 : ABI_FREE(vectornd_dir)
2224 : :
2225 : : ! scale by 2\pi
2226 [ + - ]: 143904 : if(gpu_option==ABI_GPU_DISABLED) then
2227 [ + + + + ]: 91638552 : 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 [ - + ]: 143904 : 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 [ - + ]: 200736 : ABI_FREE(work)
2376 [ - + ]: 200736 : if (usezora) then
2377 [ # # ]: 0 : ABI_FREE(zk)
2378 : : end if
2379 : :
2380 [ - + ]: 200736 : 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 : : !!***
|
