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