Branch data Line data Source code
1 : : !!****m* ABINIT/m_opernlb_ylm_mv
2 : : !! NAME
3 : : !! m_opernlb_ylm_mv
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2021-2026 ABINIT group (LB,MT)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_opernlb_ylm_mv
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_errors
26 : : #if defined HAVE_OPENMP
27 : : use OMP_LIB
28 : : #endif
29 : :
30 : : implicit none
31 : :
32 : : private
33 : : !!***
34 : :
35 : : public :: opernlb_ylm_mv
36 : : integer,public,save :: opernlb_mv_counter = -1
37 : : integer,public,save :: opernlb_mv_dgemv_counter = -1
38 : : !!***
39 : :
40 : : contains
41 : : !!***
42 : :
43 : : !!****f* ABINIT/opernlb_ylm_mv
44 : : !! NAME
45 : : !! opernlb_ylm_mv
46 : : !!
47 : : !! FUNCTION
48 : : !! "matrix-vector" alternative implementation of "opernlrb_ylm".
49 : : !!
50 : : !! * Operate with the non-local part of the hamiltonian,
51 : : !! from projected scalars to reciprocal space.
52 : : !! * Operate with the non-local projectors and the overlap matrix,
53 : : !! from projected scalars to reciprocal space.
54 : : !!
55 : : !! The input is gxfac (gxfac_sij):
56 : : !! gxfac(lmn) = Sum_l'm'n' D_l'm'n'.<p_l'm'n|c> (or S_l'm'n' for gxfac_sij)
57 : : !! and here we compute :
58 : : !! Sum_lmn <g|p_lmn> gxfac(lmn) = 4pi/sqrt(vol) exp(-2pi.i.g.R) Sum_lmn (-i)^l f_nl(g).Y_lm(g) gxfac(lmn)
59 : : !! Here this is done in 3 steps:
60 : : !! (1) compute for every lmn : gxfac_(lmn) = 4pi/sqrt(vol).(-i)^l.gxfac(lmn)
61 : : !! (2) compute for every g : scal(g) = Sum_lmn f_nl(g).Y_lm(g).gxfac_(lmn)
62 : : !! (3) compute for every g : vect(g) = exp(-2pi.i.g.R).scal(g)
63 : : !!
64 : : !! Step (2) is a real-matrix/complex-vector multiplication, here two options are possible:
65 : : !! - case nloalg(1)=2 : compute the real and imaginary parts separately using two calls of DGMEV
66 : : !! - case nloalg(1)=3 : in order to read the matrix only once, we compute both real and imaginary parts at the same time "by hand"
67 : : !!
68 : : !! Depending on the achitecture and the available blas library, one option could be more interesting than an other...
69 : : !!
70 : : !! INPUTS
71 : : !! choice=chooses possible output (see below)
72 : : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
73 : : !! 2 if <p_lmn|c> scalars are complex
74 : : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
75 : : !! dimffnl=second dimension of ffnl
76 : : !! ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
77 : : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
78 : : !! gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))= reduced projected scalars related to Sij (overlap)
79 : : !! ia3=gives the number of the first atom in the subset presently treated
80 : : !! indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
81 : : !! matblk=dimension of the array ph3d
82 : : !! nincat=number of atoms in the subset here treated
83 : : !! nlmn=number of (l,m,n) numbers for current type of atom
84 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
85 : : !! npw=number of plane waves in reciprocal space
86 : : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
87 : : !! paw_opt= define the nonlocal operator concerned with:
88 : : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
89 : : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
90 : : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
91 : : !! paw_opt=3 : PAW overlap matrix (Sij)
92 : : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
93 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
94 : : !! ucvol=unit cell volume (bohr^3)
95 : : !!
96 : : !! OUTPUT
97 : : !! (see side effects)
98 : : !!
99 : : !! SIDE EFFECTS
100 : : !! --if (paw_opt=0)
101 : : !! vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
102 : : !! or one of its derivatives to the input vect.
103 : : !! --if (paw_opt=0, 1 or 4)
104 : : !! vect(2,npwout*nspinor)=result of the aplication of the concerned operator
105 : : !! or one of its derivatives to the input vect.:
106 : : !! if (choice=1) <G|V_nonlocal|vect_in>
107 : : !! if (paw_opt=2)
108 : : !! vect(2,npwout*nspinor)=final vector in reciprocal space:
109 : : !! if (choice=1) <G|V_nonlocal-lamdba.(I+S)|vect_in> (note: not including <G|I|c>)
110 : : !! --if (paw_opt=3 or 4)
111 : : !! svect(2,npwout*nspinor)=result of the aplication of Sij (overlap matrix)
112 : : !! or one of its derivatives to the input vect.:
113 : : !! if (choice=1) <G|I+S|vect_in> (note: not including <G|I|c>)
114 : : !!
115 : : !! NOTES
116 : : !! 1-No openMP available for now
117 : : !! 2-Operate for one type of atom, and within this given type of atom,
118 : : !! for a subset of at most nincat atoms.
119 : : !! 3-projector derivatives (abs(choice)>1) are not implemented yet
120 : : !!
121 : : !! SOURCE
122 : :
123 : 0 : subroutine opernlb_ylm_mv(choice,cplex,cplex_fac,&
124 : 0 : & dimffnl,ffnl,gxfac,gxfac_sij,&
125 : 0 : & ia3,indlmn,matblk,nincat,nlmn,nloalg,npw,&
126 [ # # # # ]: 0 : & nspinor,paw_opt,ph3d,svect,ucvol,vect)
127 : :
128 : : !Arguments ------------------------------------
129 : : !scalars
130 : : integer,intent(in) :: choice,cplex,cplex_fac,dimffnl,ia3,matblk,nincat
131 : : integer,intent(in) :: nlmn,npw,nspinor,paw_opt
132 : : real(dp),intent(in) :: ucvol
133 : : !arrays
134 : : integer,intent(in) :: indlmn(6,nlmn),nloalg(3)
135 : : real(dp),intent(in),target :: ffnl(npw,dimffnl,nlmn)
136 : : real(dp),intent(in) :: gxfac(cplex_fac,nlmn,nincat,nspinor)
137 : : real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
138 : : real(dp),intent(in) :: ph3d(2,npw,matblk)
139 : : real(dp),intent(inout) :: svect(:,:),vect(:,:)
140 : : !Local variables-------------------------------
141 : : !Arrays
142 : : !scalars
143 : : logical :: use_dgemv
144 : : integer :: ia,iaph3d
145 : : integer :: il,ilmn,ipw,jpw,ipwshft,ispinor
146 : : real(dp) :: buffer_r,buffer_i,wt
147 : : !arrays
148 : : ! real(dp) :: tsec(2)
149 : 0 : real(dp) :: gxfac_(nlmn,2),gxfacs_(nlmn,2)
150 : 0 : real(dp),allocatable :: scalr(:),scali(:)
151 : 0 : real(dp),pointer :: ffnl_loc(:,:)
152 : : complex(dp) :: ctmp, cil(4)
153 : :
154 : : ! *************************************************************************
155 : :
156 : : DBG_ENTER("COLL")
157 : :
158 : : !Some checks
159 : : ! nthreads=1
160 : : !#if defined HAVE_OPENMP
161 : : ! nthreads=OMP_GET_NUM_THREADS()
162 : : !#endif
163 : : ! if (nthreads>1) then
164 : : ! ABI_ERROR('Only nthreads=1 is available for now.')
165 : : ! end if
166 : :
167 : 0 : if (abs(choice)>1) then
168 : 0 : ABI_ERROR('Only abs(choice)<=1 is available for now.')
169 : : end if
170 [ # # ]: 0 : if (nloalg(1)<2.or.nloalg(1)>10) then
171 : 0 : ABI_ERROR('nloalg(1) should be between 2 or 10.')
172 : : end if
173 : :
174 [ # # ]: 0 : use_dgemv = nloalg(1)==2.or.nloalg(1)==6.or.nloalg(1)==10
175 : : if (use_dgemv) then
176 [ # # ]: 0 : if(opernlb_mv_dgemv_counter>=0) then
177 : 0 : opernlb_mv_dgemv_counter = opernlb_mv_dgemv_counter + 1
178 [ # # ]: 0 : if (paw_opt==4) opernlb_mv_dgemv_counter = opernlb_mv_dgemv_counter + 1
179 : : end if
180 : : else
181 [ # # ]: 0 : if(opernlb_mv_counter>=0) then
182 : 0 : opernlb_mv_counter = opernlb_mv_counter + 1
183 [ # # ]: 0 : if (paw_opt==4) opernlb_mv_counter = opernlb_mv_counter + 1
184 : : end if
185 : : end if
186 : :
187 : : !Inits
188 : 0 : wt=four_pi/sqrt(ucvol)
189 : :
190 : 0 : ffnl_loc => ffnl(:,1,:)
191 : :
192 : :
193 [ # # # # : 0 : ABI_MALLOC(scalr,(npw))
# # # # #
# # # ]
194 [ # # # # : 0 : ABI_MALLOC(scali,(npw))
# # # # ]
195 : :
196 : : !$OMP PARALLEL PRIVATE(buffer_r,buffer_i,cil,il,ilmn,ipw,jpw,ctmp), &
197 : : !$OMP PRIVATE(ispinor,ipwshft,ia,iaph3d,gxfac_,gxfacs_)
198 : :
199 : : ! (-i)^l
200 : 0 : cil(1) = ( 1.0_DP, 0.0_DP) * wt
201 : 0 : cil(2) = ( 0.0_DP,-1.0_DP) * wt
202 : 0 : cil(3) = (-1.0_DP, 0.0_DP) * wt
203 : 0 : cil(4) = ( 0.0_DP, 1.0_DP) * wt
204 : :
205 : : ! if (paw_opt/=3) then
206 : : ! ABI_MALLOC(gxfac_,(nlmn,2))
207 : : ! end if
208 : : ! if (paw_opt>=3) then
209 : : ! ABI_MALLOC(gxfacs_,(nlmn,2))
210 : : ! end if
211 : :
212 : : !Loop on spinorial components
213 [ # # ]: 0 : do ispinor=1,nspinor
214 : 0 : ipwshft=(ispinor-1)*npw
215 : :
216 : : ! Loop on atoms (blocking)
217 [ # # ]: 0 : do ia=1,nincat
218 [ # # ]: 0 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
219 : : ! Step (1) : scale gxfac with 4pi/sqr(omega).(-i)^l
220 [ # # ]: 0 : if (paw_opt/=3) then
221 [ # # ]: 0 : if (cplex_fac==2) then
222 [ # # ]: 0 : do ilmn=1,nlmn
223 : 0 : il=mod(indlmn(1,ilmn),4)+1
224 : 0 : ctmp = cil(il) * cmplx( gxfac(1,ilmn,ia,ispinor), gxfac(2,ilmn,ia,ispinor), kind=DP )
225 : 0 : gxfac_(ilmn,1) = real(ctmp)
226 : 0 : gxfac_(ilmn,2) = aimag(ctmp)
227 : : end do
228 [ # # ]: 0 : else if (cplex_fac==1) then
229 [ # # ]: 0 : do ilmn=1,nlmn
230 : 0 : il=mod(indlmn(1,ilmn),4)+1
231 : 0 : ctmp = cil(il) * gxfac(1,ilmn,ia,ispinor)
232 : 0 : gxfac_(ilmn,1) = real(ctmp)
233 : 0 : gxfac_(ilmn,2) = aimag(ctmp)
234 : : end do
235 : : else
236 : 0 : ABI_BUG('Error : should not be possible to be here')
237 : : end if
238 : : end if
239 : :
240 : : ! Step (1) bis: Scale gxfac_sij with 4pi/sqr(omega).(-i)^l
241 [ # # ]: 0 : if (paw_opt>=3) then
242 [ # # ]: 0 : if (cplex==2) then
243 [ # # ]: 0 : do ilmn=1,nlmn
244 : 0 : il=mod(indlmn(1,ilmn),4)+1
245 : 0 : ctmp = cil(il) * cmplx( gxfac_sij(1,ilmn,ia,ispinor), gxfac_sij(2,ilmn,ia,ispinor), kind=DP )
246 : 0 : gxfacs_(ilmn,1) = real(ctmp)
247 : 0 : gxfacs_(ilmn,2) = aimag(ctmp)
248 : : end do
249 [ # # ]: 0 : else if (cplex==1) then
250 [ # # ]: 0 : do ilmn=1,nlmn
251 : 0 : il=mod(indlmn(1,ilmn),4)+1
252 : 0 : ctmp = cil(il) * gxfac_sij(1,ilmn,ia,ispinor)
253 : 0 : gxfacs_(ilmn,1) = real(ctmp)
254 : 0 : gxfacs_(ilmn,2) = aimag(ctmp)
255 : : end do
256 : : else
257 : 0 : ABI_BUG('Error : should not be possible to be here')
258 : : end if
259 : : end if
260 : :
261 : : ! Compute <g|Vnl|c> (or derivatives) for each plane wave:
262 [ # # ]: 0 : if (paw_opt/=3) then
263 : :
264 : : ! Step (2) scal(g) = Sum_lmn f_nl(g).Y_lm(g).gxfac_(lmn)
265 [ # # ]: 0 : if (use_dgemv) then
266 [ # # ]: 0 : call DGEMV('N',npw,nlmn,1.0_DP,ffnl_loc,npw,gxfac_(:,1),1,0.0_DP,scalr,1)
267 [ # # ]: 0 : call DGEMV('N',npw,nlmn,1.0_DP,ffnl_loc,npw,gxfac_(:,2),1,0.0_DP,scali,1)
268 : : else
269 : : ! scalr(:) = zero
270 : : ! scali(:) = zero
271 : : ! do ilmn=1,nlmn
272 : : ! do ipw=1,npw
273 : : ! scalr(ipw) = scalr(ipw) + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,1)
274 : : ! scali(ipw) = scali(ipw) + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,2)
275 : : ! end do
276 : : ! end do
277 : : !$OMP DO
278 [ # # ]: 0 : do ipw=1,npw
279 : : buffer_r = zero
280 : : buffer_i = zero
281 [ # # ]: 0 : do ilmn=1,nlmn
282 : 0 : buffer_r = buffer_r + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,1)
283 : 0 : buffer_i = buffer_i + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,2)
284 : : ! scalr(ipw) = scalr(ipw) + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,1)
285 : : ! scali(ipw) = scali(ipw) + ffnl_loc(ipw,ilmn) * gxfac_(ilmn,2)
286 : : end do
287 : 0 : scalr(ipw) = buffer_r
288 : 0 : scali(ipw) = buffer_i
289 : : end do
290 : : !$OMP END DO
291 : : end if
292 : :
293 : : ! Step (3) : vect(g) = exp(-2pi.i.g.R).scal(g)
294 : : !$OMP DO
295 [ # # ]: 0 : do ipw=1,npw
296 : 0 : jpw=ipw+ipwshft
297 : 0 : vect(1,jpw)=vect(1,jpw)+scalr(ipw)*ph3d(1,ipw,iaph3d)+scali(ipw)*ph3d(2,ipw,iaph3d)
298 : 0 : vect(2,jpw)=vect(2,jpw)-scalr(ipw)*ph3d(2,ipw,iaph3d)+scali(ipw)*ph3d(1,ipw,iaph3d)
299 : : end do
300 : : !$OMP END DO
301 : :
302 : : end if
303 : :
304 : : ! Compute <g|S|c> (or derivatives) for each plane wave:
305 [ # # ]: 0 : if (paw_opt>=3) then
306 : :
307 : : ! Step (2) (bis) scal(g) = Sum_lmn f_nl(g).Y_lm(g).gxfacs_(lmn)
308 [ # # ]: 0 : if (nloalg(1)==3) then
309 : : ! scalr(:) = zero
310 : : ! scali(:) = zero
311 : : ! do ilmn=1,nlmn
312 : : ! do ipw=1,npw
313 : : ! scalr(ipw) = scalr(ipw) + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,1)
314 : : ! scali(ipw) = scali(ipw) + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,2)
315 : : ! end do
316 : : ! end do
317 : : !$OMP DO
318 [ # # ]: 0 : do ipw=1,npw
319 : : buffer_r = zero
320 : : buffer_i = zero
321 [ # # ]: 0 : do ilmn=1,nlmn
322 : 0 : buffer_r = buffer_r + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,1)
323 : 0 : buffer_i = buffer_i + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,2)
324 : : ! scalr(ipw) = scalr(ipw) + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,1)
325 : : ! scali(ipw) = scali(ipw) + ffnl_loc(ipw,ilmn) * gxfacs_(ilmn,2)
326 : : end do
327 : 0 : scalr(ipw) = buffer_r
328 : 0 : scali(ipw) = buffer_i
329 : : end do
330 : : !$OMP END DO
331 [ # # ]: 0 : else if (nloalg(1)==2) then
332 [ # # ]: 0 : call DGEMV('N',npw,nlmn,1.0_DP,ffnl_loc,npw,gxfacs_(:,1),1,0.0_DP,scalr,1)
333 [ # # ]: 0 : call DGEMV('N',npw,nlmn,1.0_DP,ffnl_loc,npw,gxfacs_(:,2),1,0.0_DP,scali,1)
334 : : end if
335 : :
336 : : ! Step (3) (bis) : svect(g) = exp(-2pi.i.g.R).scal(g)
337 : : !$OMP DO
338 [ # # ]: 0 : do ipw=1,npw
339 : 0 : jpw=ipw+ipwshft
340 : 0 : svect(1,jpw)=svect(1,jpw)+scalr(ipw)*ph3d(1,ipw,iaph3d)+scali(ipw)*ph3d(2,ipw,iaph3d)
341 : 0 : svect(2,jpw)=svect(2,jpw)-scalr(ipw)*ph3d(2,ipw,iaph3d)+scali(ipw)*ph3d(1,ipw,iaph3d)
342 : : end do
343 : : !$OMP END DO
344 : :
345 : : end if
346 : :
347 : : ! End loop on atoms
348 : : end do
349 : : end do ! End loop on spinors
350 : : ! if (paw_opt/=3) then
351 : : ! ABI_FREE(gxfac_)
352 : : ! end if
353 : : ! if (paw_opt>=3) then
354 : : ! ABI_FREE(gxfacs_)
355 : : ! end if
356 : : !$OMP END PARALLEL
357 : :
358 [ # # ]: 0 : ABI_FREE(scalr)
359 [ # # ]: 0 : ABI_FREE(scali)
360 : :
361 : :
362 : : DBG_EXIT("COLL")
363 : :
364 [ # # ]: 0 : end subroutine opernlb_ylm_mv
365 : : !!***
366 : :
367 : : end module m_opernlb_ylm_mv
368 : : !!***
|
