Line data Source code
1 : !!****m* ABINIT/m_mkffnl
2 : !! NAME
3 : !! m_mkffnl
4 : !!
5 : !! FUNCTION
6 : !! Make FFNL, nonlocal form factors, for each type of atom up to ntypat
7 : !! and for each angular momentum.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MT, DRH)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_mkffnl
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_splines
29 : use m_xmpi
30 :
31 : use m_time, only : timab
32 : use m_kg, only : mkkin
33 : use m_sort, only : sort_dp
34 : use defs_datatypes, only : pseudopotential_type
35 : use m_crystal, only : crystal_t
36 :
37 : implicit none
38 :
39 : private
40 : !!***
41 :
42 : public :: mkffnl_objs
43 : public :: mkffnl
44 : !!***
45 :
46 : contains
47 : !!***
48 :
49 : !!****f* ABINIT/mkffnl_objs
50 : !! NAME
51 : !! mkffnl_objs
52 : !!
53 : !! FUNCTION
54 : !! Simplified wrapper around mkffnl in which input parameters are passed via crystal_t and pseudopotential_type.
55 : !!
56 : !! INPUTS
57 : !! See mkffnl
58 : !!
59 : !! OUTPUT
60 : !! ffnl(npw,dimffnl,lmnmax,ntypat)=described below
61 : !! [request]=Used in conjunction with [comm] to perform non-blocking xmpi_isum_ip. Client code must
62 : !! wait on request before using ffnl. If not present, blocking API is used.
63 : !!
64 : !! SOURCE
65 :
66 158544 : subroutine mkffnl_objs(cryst, psps, dimffnl, ffnl, ider, idir, kg_k, kpg, kpt, nkpg, npw_k, ylm_k, ylm_gr_k, &
67 : comm, request) ! optional
68 :
69 : !Arguments ------------------------------------
70 : !scalars
71 : type(crystal_t),intent(in) :: cryst
72 : type(pseudopotential_type),intent(in) :: psps
73 : integer,intent(in) :: dimffnl, ider, idir, npw_k, nkpg
74 : integer,optional,intent(in) :: comm
75 : integer ABI_ASYNC, optional,intent(out):: request
76 : !arrays
77 : integer,intent(in) :: kg_k(3,npw_k)
78 : real(dp),intent(in) :: kpg(npw_k, nkpg), kpt(3), ylm_k(:,:), ylm_gr_k(:,:,:)
79 : real(dp),intent(out) :: ffnl(npw_k, dimffnl, psps%lmnmax, psps%ntypat)
80 : !
81 : !!Local variables-------------------------------
82 : integer :: my_comm
83 : ! *************************************************************************
84 :
85 158544 : my_comm = xmpi_comm_self; if (present(comm)) my_comm = comm
86 :
87 158544 : if (present(request)) then
88 : call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, &
89 : cryst%gmet, cryst%gprimd, ider, idir, psps%indlmn, kg_k, kpg, kpt, psps%lmnmax, &
90 : psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw_k, psps%ntypat, &
91 : psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm_k, ylm_gr_k, &
92 0 : comm=comm, request=request)
93 :
94 : else
95 : call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, &
96 : cryst%gmet, cryst%gprimd, ider, idir, psps%indlmn, kg_k, kpg, kpt, psps%lmnmax, &
97 : psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw_k, psps%ntypat, &
98 : psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm_k, ylm_gr_k, &
99 158544 : comm=comm)
100 : end if
101 :
102 158544 : end subroutine mkffnl_objs
103 : !!***
104 :
105 : !!****f* ABINIT/mkffnl
106 : !! NAME
107 : !! mkffnl
108 : !!
109 : !! FUNCTION
110 : !! Make FFNL, nonlocal form factors, for each type of atom up to ntypat
111 : !! and for each angular momentum.
112 : !! When Legendre polynomials are used in the application of the
113 : !! nonlocal operator, FFNLs depend on (l,n) components; in this
114 : !! case, form factors are real and divided by |k+G|^l;
115 : !! When spherical harmonics are used, FFNLs depend on (l,m,n)
116 : !! components; in this case, form factors are multiplied by Ylm(k+G).
117 : !!
118 : !! INPUTS
119 : !! dimekb=second dimension of ekb (see ekb)
120 : !! dimffnl=second dimension of ffnl (1+number of derivatives)
121 : !! ekb(dimekb,ntypat*(1-usepaw))=(Real) Kleinman-Bylander energies (hartree)
122 : !! ->NORM-CONSERVING PSPS ONLY
123 : !! ffspl(mqgrid,2,lnmax,ntypat)=form factors and spline fit to 2nd derivative
124 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
125 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
126 : !! ider=0=>no derivative wanted; 1=>1st derivative wanted; 2=>1st and 2nd derivatives wanted
127 : !! idir=ONLY WHEN YLMs ARE USED:
128 : !! When 1st derivative has to be computed: (see more info below)
129 : !! - Determine the direction(s) of the derivatives(s)
130 : !! - Determine the set of coordinates (reduced or cartesians)
131 : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,spin for i=ln (if useylm=0)
132 : !! or i=lmn (if useylm=1)
133 : !! [kinpw(npw)]=plane wave kinetic energy (useless here) and filter mask for dilatmx>1 (needed here)
134 : !! kg(3,npw)=integer coordinates of planewaves in basis sphere for this k point.
135 : !! kpg(npw,nkpg)= (k+G) components (only if useylm=1)
136 : !! kpt(3)=reduced coordinates of k point
137 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
138 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
139 : !! lnmax=max. number of (l,n) components over all type of psps
140 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
141 : !! mqgrid=size of q (or |G|) grid for f(q)
142 : !! nkpg=second dimension of kpg_k (0 if useylm=0)
143 : !! npw=number of planewaves in basis sphere
144 : !! ntypat=number of types of atoms
145 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
146 : !! pspso(ntypat)=spin-orbit characteristics for each atom type (1, 2, or 3)
147 : !! qgrid(mqgrid)=uniform grid of q values from 0 to qmax
148 : !! rmet(3,3)=real space metric (bohr**2)
149 : !! useylm=governs the way the nonlocal operator is to be applied:
150 : !! 1=using Ylm, 0=using Legendre polynomials
151 : !! ylm (npw,mpsang*mpsang*useylm)=real spherical harmonics for each G and k point
152 : !! ylm_gr(npw,3,mpsang*mpsang*useylm)=gradients of real spherical harmonics wrt (k+G)
153 : !! [comm]=MPI communicator. Default: xmpi_comm_self.
154 : !!
155 : !! OUTPUT
156 : !! ffnl(npw,dimffnl,lmnmax,ntypat)=described below
157 : !! [request]=Used in conjunction with [comm] to perform non-blocking xmpi_isum_ip. Client code must
158 : !! wait on request before using ffnl. If not present, blocking API is used.
159 : !!
160 : !! NOTES
161 : !! Uses spline fit ffspl provided by Numerical Recipes spline subroutine.
162 : !! Form factor $f_l(q)$ is defined by
163 : !! \begin{equation}
164 : !! \textrm{f}_l(q)=\frac{1}{dvrms} \int_0^\infty [j_l(2 \pi r q) u_l(r) dV(r) r dr]
165 : !! \end{equation}
166 : !! where u_l(r)=reference state wavefunction, dV(r)=nonlocal psp
167 : !! correction, j_l(arg)=spherical Bessel function for angular momentum l,
168 : !! and
169 : !! \begin{equation}
170 : !! \textrm{dvrms} = \int_0^\infty [(u_l(r) dV(r))^2 dr])^{1/2}
171 : !! \end{equation}
172 : !! which is square root of mean square dV, i.e.
173 : !! $ (\langle (dV)^2 \rangle)^{1/2} $ .
174 : !! This routine is passed f_l(q) in spline form in the array ffspl and then
175 : !! constructs the values of $f_l(q)$ on the relevant (k+G) in array ffnl.
176 : !! The evaluation of the integrals defining ffspl was done in mkkbff.
177 : !!
178 : !! Delivers the following (for each atom type t, or itypat):
179 : !! --------------------------
180 : !! Using Legendre polynomials in the application of nl operator:
181 : !! ffnl are real.
182 : !! ffnl(ig,1,(l,0,n),itypat) $= f_ln(k+G)/|k+G|^l $
183 : !! === if ider>=1
184 : !! ffnl(ig,2,(l,0,n),itypat) $=(fprime_ln(k+G)-l*f_ln(k+G)/|k+G|)/|k+G|^(l+1) $
185 : !! === if ider=2
186 : !! ffnl(ig,3,(l,0,n),itypat) $=(fprimeprime_ln(k+G)-(2l+1)*fprime_ln(k+G)/|k+G|
187 : !! +l(l+2)*f_ln(k+G)/|k+G|**2)/|k+G|^(l+2)
188 : !! --------------------------
189 : !! Using spherical harmonics in the application of nl operator:
190 : !! ffnl are real (we use REAL spherical harmonics).
191 : !! ffnl(ig,1,(l,m,n),itypat) = ffnl_1
192 : !! $= f_ln(k+G) * Y_lm(k+G) $
193 : !! === if ider>=1
194 : !! --if (idir==0)
195 : !! ffnl(ig,1+i,(l,m,n),itypat) = dffnl_i = 3 reduced coord. of d(ffnl_1)/dK^cart
196 : !! $= fprime_ln(k+G).Y_lm(k+G).(k+G)^red_i/|k+G|+f_ln(k+G).(dY_lm/dK^cart)^red_i $
197 : !! for i=1..3
198 : !! --if (1<=idir<=3)
199 : !! ffnl(ig,2,(l,m,n),itypat)= cart. coordinate idir of d(ffnl_1)/dK^red
200 : !! = Sum_(mu,nu) [ Gprim(mu,idir) Gprim(mu,nu) dffnl_nu ]
201 : !! --if (idir==4)
202 : !! ffnl(ig,1+i,(l,m,n),itypat)= 3 cart. coordinates of d(ffnl_1)/dK^red
203 : !! = Sum_(mu,nu) [ Gprim(mu,i) Gprim(mu,nu) dffnl_nu ]
204 : !! --if (-6<idir<-1)
205 : !! ffnl(ig,2,(l,m,n),itypat)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
206 : !! with d(ffnl)/dK^cart_i = Sum_nu [ Gprim(nu,i) dffnl_nu ]
207 : !! for |idir|->(mu,nu) (1->11,2->22,3->33,4->32,5->31,6->21)
208 : !! --if (idir==-7)
209 : !! ffnl(ig,2:7,(l,m,n),itypat)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
210 : !! with d(ffnl)/dK^cart_i = Sum_nu [ Gprim(nu,i) dffnl_nu ]
211 : !! for all (mu,nu) (6 independant terms)
212 : !! === if ider==2
213 : !! --if (idir==0)
214 : !! ffnl(ig,4+i,(l,m,n),itypat) = d2ffnl_mu,nu = 6 reduced coord. of d2(ffnl_1)/dK^cart.dK^cart
215 : !! for all i=(mu,nu) (6 independant terms)
216 : !! --if (idir==4)
217 : !! ffnl(ig,4+i,(l,m,n),itypat) = d2ffnl_i =6 cart. coordinates of d2(ffnl_1)/dK^red.dK^red
218 : !! for all i=(mu,nu) (6 independant terms)
219 : !! = Sum_(mu1,mu2,mu3,mu4) [ Gprim(mu1,mu) Gprim(mu2,nu) Gprim(mu1,mu3) Gprim(mu2,mu4) d2ffnl_mu3,mu4 ]
220 : !! --------------------------
221 : !!
222 : !! 1) l may be 0, 1, 2, or 3 in this version.
223 : !!
224 : !! 2) Norm-conserving psps: only FFNL for which ekb is not zero are calculated.
225 : !!
226 : !! 3) Each expression above approaches a constant as $|k+G| \rightarrow 0 $.
227 : !! In the cases where $|k+G|$ is in the denominator, there is always a
228 : !! factor of $(k+G)_mu$ multiplying the ffnl term where it is actually used,
229 : !! so that we may replace the ffnl term by any constant when $|k+G| = 0$.
230 : !! Below we replace 1/0 by 1/tol10, thus creating an arbitrary constant
231 : !! which will later be multiplied by 0.
232 : !!
233 : !! TODO
234 : !! Some parts can be rewritten with BLAS1 calls.
235 : !!
236 : !! SOURCE
237 :
238 4740506 : subroutine mkffnl(dimekb, dimffnl, ekb, ffnl, ffspl, gmet, gprimd, ider, idir, indlmn, &
239 4740506 : kg, kpg, kpt, lmnmax, lnmax, mpsang, mqgrid, nkpg, npw, ntypat, pspso, &
240 2370253 : qgrid, rmet, usepaw, useylm, ylm, ylm_gr, &
241 2370253 : comm, request, kinpw) ! optional
242 :
243 : !Arguments ------------------------------------
244 : !scalars
245 : integer,intent(in) :: dimekb,dimffnl,ider,idir,lmnmax,lnmax,mpsang,mqgrid,nkpg
246 : integer,intent(in) :: npw,ntypat,usepaw,useylm
247 : integer,optional,intent(in) :: comm
248 : real(dp),optional,intent(in) :: kinpw(:)
249 : integer ABI_ASYNC, optional,intent(out):: request
250 : !arrays
251 : integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg(3,npw),pspso(ntypat)
252 : real(dp),intent(in) :: ekb(dimekb,ntypat*(1-usepaw))
253 : real(dp),intent(in) :: ffspl(mqgrid,2,lnmax,ntypat),gmet(3,3),gprimd(3,3)
254 : real(dp),intent(in) :: kpg(npw,nkpg),kpt(3),qgrid(mqgrid),rmet(3,3)
255 : real(dp),intent(in) :: ylm(:,:),ylm_gr(:,:,:)
256 : real(dp),intent(out) :: ffnl(npw,dimffnl,lmnmax,ntypat)
257 : ! MG: Should be ABI_ASYNC due to optional non-Blocking API but NAG complains
258 : ! Error: m_d2frnl.F90, line 600: Array section FFNL_STR(:,:,:,:,MU) supplied for dummy FFNL (no. 4) of MKFFNL,
259 : ! the dummy is ASYNCHRONOUS but not assumed-shape
260 : ! so we declare request as ASYNCHRONOUS
261 :
262 : !Local variables-------------------------------
263 : !scalars
264 : integer :: ider_tmp,iffnl,ig,ig0,il,ilm,ilmn,iln,iln0,im,iylm,itypat,mu,mua,mub,nlmn,nu,nua,nub
265 : integer :: nprocs, my_rank, cnt, ierr
266 : real(dp),parameter :: renorm_factor=0.5d0/pi**2,tol_norm=tol10
267 : real(dp) :: ecut,ecutsm,effmass_free,fact,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3,rmetab,yp1
268 : logical :: testnl=.false.
269 : character(len=500) :: msg
270 : !arrays
271 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
272 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
273 : real(dp) :: rprimd(3,3),tsec(2)
274 2370253 : real(dp),allocatable :: dffnl_cart(:,:),dffnl_red(:,:),dffnl_tmp(:)
275 2370253 : real(dp),allocatable :: d2ffnl_cart(:,:),d2ffnl_red(:,:),d2ffnl_tmp(:)
276 2370253 : real(dp),allocatable :: kpgc(:,:),kpgn(:,:),kpgnorm(:),kpgnorm_inv(:)
277 2370253 : real(dp),allocatable :: wk_ffnl1(:),wk_ffnl2(:),wk_ffnl3(:),wk_ffspl(:,:)
278 :
279 : ! *************************************************************************
280 :
281 : ! Keep track of time spent in mkffnl
282 2370253 : call timab(16, 1, tsec)
283 :
284 2370253 : nprocs = 1; my_rank = 0
285 2370253 : if (present(comm)) then
286 157025 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
287 : end if
288 :
289 : ! Compatibility tests
290 : !if (mpsang>4) then
291 : ! write(msg,'(a,i0,a,a)')&
292 : ! 'Called with mpsang > 4, =',mpsang,ch10,&
293 : ! 'This subroutine will not accept lmax+1 > 4.'
294 : ! ABI_BUG(msg)
295 : !end if
296 2370253 : if (idir<-7.or.idir>4) then
297 0 : ABI_BUG('Called with idir<-6 or idir>4 !')
298 : end if
299 2370253 : if (useylm==0) then
300 1546807 : iffnl=1+ider
301 : else
302 823446 : iffnl=1
303 823446 : if (ider>=1) then
304 236428 : if (idir==0) iffnl=iffnl+3
305 236428 : if (idir/=0) iffnl=iffnl+1
306 236428 : if (idir==4) iffnl=iffnl+2
307 236428 : if (idir==-7) iffnl=iffnl+5
308 : end if
309 823446 : if (ider==2) then
310 67936 : if (idir==0) iffnl=iffnl+6
311 67936 : if (idir==4) iffnl=iffnl+6
312 : end if
313 : end if
314 2370253 : if (iffnl/=dimffnl) then
315 0 : write(msg,'(2(a,i1),a,i2)') 'Incompatibility between ider, idir and dimffnl : ider = ',ider,&
316 0 : ' idir = ',idir,' dimffnl = ',dimffnl
317 0 : ABI_BUG(msg)
318 : end if
319 2370253 : if (useylm==1) then
320 823446 : ABI_CHECK(size(ylm,1)==npw,'BUG: wrong ylm size (1)')
321 823446 : ABI_CHECK(size(ylm,2)==mpsang**2,'BUG: wrong ylm size (2)')
322 823446 : if(ider>0)then
323 236428 : ABI_CHECK(size(ylm_gr,1)==npw,'BUG: wrong ylm_gr size (1)')
324 236428 : ABI_CHECK(size(ylm_gr,2)>=3+6*(ider/2),'BUG: wrong ylm_gr size (2)')
325 236428 : ABI_CHECK(size(ylm_gr,3)==mpsang**2,'BUG: wrong ylm_gr size (3)')
326 : end if
327 : end if
328 :
329 : ! Get (k+G) and |k+G|
330 7110759 : ABI_MALLOC(kpgnorm,(npw))
331 4740506 : ABI_MALLOC(kpgnorm_inv,(npw))
332 :
333 2370253 : ig0=-1 ! index of |k+g|=0 vector
334 :
335 2370253 : if (useylm==1) then
336 2470338 : ABI_MALLOC(kpgc,(npw,3))
337 823446 : if (ider>=1) then
338 472856 : ABI_MALLOC(kpgn,(npw,3))
339 : end if
340 823446 : if (nkpg<3) then
341 : !$OMP PARALLEL DO PRIVATE(ig,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3)
342 85011489 : do ig=1,npw
343 84418196 : kpg1=kpt(1)+dble(kg(1,ig))
344 84418196 : kpg2=kpt(2)+dble(kg(2,ig))
345 84418196 : kpg3=kpt(3)+dble(kg(3,ig))
346 84418196 : kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
347 84418196 : kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
348 84418196 : kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
349 84418196 : kpgc(ig,1)=kpgc1
350 84418196 : kpgc(ig,2)=kpgc2
351 84418196 : kpgc(ig,3)=kpgc3
352 84418196 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
353 84418196 : if (kpgnorm(ig)<=tol_norm) ig0=ig
354 85011489 : if (ider>=1) then
355 25542469 : kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
356 25542469 : kpgn(ig,1)=kpg1*kpgnorm_inv(ig)
357 25542469 : kpgn(ig,2)=kpg2*kpgnorm_inv(ig)
358 25542469 : kpgn(ig,3)=kpg3*kpgnorm_inv(ig)
359 : end if
360 : end do
361 : else
362 : !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
363 25899681 : do ig=1,npw
364 25669528 : kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
365 25669528 : kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
366 25669528 : kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
367 25669528 : kpgc(ig,1)=kpgc1
368 25669528 : kpgc(ig,2)=kpgc2
369 25669528 : kpgc(ig,3)=kpgc3
370 25669528 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
371 25669528 : if (kpgnorm(ig)<=tol_norm) ig0=ig
372 25899681 : if (ider>=1) then
373 4166928 : kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
374 16667712 : kpgn(ig,1:3)=kpg(ig,1:3)*kpgnorm_inv(ig)
375 : end if
376 : end do
377 : end if
378 : else
379 1546807 : if (nkpg<3) then
380 1546807 : ecut=huge(zero)*0.1d0;ecutsm=zero;effmass_free=one
381 : ! Note that with ecutsm=0, the right kinetic energy is computed
382 1546807 : call mkkin(ecut,ecutsm,effmass_free,gmet,kg,kpgnorm,kpt,npw,0,0)
383 : !$OMP PARALLEL DO
384 369824234 : do ig=1,npw
385 368277427 : kpgnorm(ig)=sqrt(renorm_factor*kpgnorm(ig))
386 368277427 : kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
387 369824234 : if (kpgnorm(ig)<=tol_norm) ig0=ig
388 : end do
389 : else
390 : !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
391 0 : do ig=1,npw
392 0 : kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
393 0 : kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
394 0 : kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
395 0 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
396 0 : kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
397 0 : if (kpgnorm(ig)<=tol_norm) ig0=ig
398 : end do
399 : end if
400 : end if
401 :
402 : ! Treat dilatmx>1 (if kinpw is given)
403 2370253 : if (present(kinpw)) then
404 311235 : if (size(kinpw)/=npw) then
405 0 : ABI_ERROR("kinpw is not consistent with npw")
406 : end if
407 : end if
408 :
409 : ! Need rprimd in some cases
410 2370253 : if (ider>=1.and.useylm==1.and.ig0>0) then
411 5124 : do mu=1,3
412 16653 : do nu=1,3
413 15372 : rprimd(mu,nu)=gprimd(mu,1)*rmet(1,nu)+gprimd(mu,2)*rmet(2,nu)+gprimd(mu,3)*rmet(3,nu)
414 : end do
415 : end do
416 : end if
417 :
418 : ! Allocate several temporary arrays
419 4740506 : ABI_MALLOC(wk_ffnl1,(npw))
420 4740506 : ABI_MALLOC(wk_ffnl2,(npw))
421 4740506 : ABI_MALLOC(wk_ffnl3,(npw))
422 7110759 : ABI_MALLOC(wk_ffspl,(mqgrid,2))
423 :
424 2370253 : if (ider>=1.and.useylm==1) then
425 709284 : ABI_MALLOC(dffnl_red,(npw,3))
426 236428 : if (idir/=0) then
427 439324 : ABI_MALLOC(dffnl_cart,(npw,3))
428 : end if
429 236428 : if (idir>0) then
430 392554 : ABI_MALLOC(dffnl_tmp,(npw))
431 : end if
432 : end if
433 2370253 : if (ider>=2 .and. useylm==1) then
434 203808 : ABI_MALLOC(d2ffnl_red,(npw,6))
435 67936 : if (idir==4) then
436 133472 : ABI_MALLOC(d2ffnl_cart,(npw,6))
437 2436989 : ABI_MALLOC(d2ffnl_tmp,(npw))
438 : end if
439 : end if
440 :
441 : ! Loop over types of atoms
442 6172339042 : ffnl = zero; cnt = 0
443 6190430 : do itypat=1,ntypat
444 :
445 : ! Loop over (l,m,n) values
446 28012804 : iln0=0; nlmn=count(indlmn(3,:,itypat)>0)
447 :
448 26338418 : do ilmn=1,nlmn
449 20147988 : il=indlmn(1,ilmn,itypat)
450 20147988 : im=indlmn(2,ilmn,itypat)
451 20147988 : ilm =indlmn(4,ilmn,itypat)
452 20147988 : iln =indlmn(5,ilmn,itypat)
453 20147988 : iffnl=ilmn;if (useylm==0) iffnl=iln
454 :
455 : ! Special case: we enter the loop in case of spin-orbit calculation
456 : ! even if the psp has no spin-orbit component.
457 20147988 : if (indlmn(6,ilmn,itypat) ==1 .or. pspso(itypat) /=0) then
458 :
459 : ! Compute FFNL only if ekb>0 or paw
460 20090020 : if (usepaw==1) testnl=.true.
461 20090020 : if (usepaw==0) testnl=(abs(ekb(iln,itypat))>tol_norm)
462 :
463 20090020 : if (testnl) then
464 20089372 : cnt = cnt + 1
465 20089372 : if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism (optional)
466 : !
467 : ! Store form factors (from ffspl)
468 : ! -------------------------------
469 : ! MG: This part is an hotspot of in the EPH code due to the large number of k-points used
470 : ! To improve memory locality, I tried to:
471 : !
472 : ! 1) call a new version of splfit that operates on wk_ffspl with shape: (2,mqgrid)
473 : ! 2) pass a sorted kpgnorm array and then rearrange the output spline
474 : !
475 : ! but I didn't manage to make it significantly faster.
476 : ! For the time being, we rely on MPI-parallelism via the optional MPI communicator.
477 :
478 20089372 : if (iln > iln0) then
479 77676771200 : wk_ffspl(:,:)=ffspl(:,:,iln,itypat)
480 12892062 : ider_tmp = min(ider, 1)
481 12892062 : call splfit(qgrid,wk_ffnl2,wk_ffspl,ider_tmp,kpgnorm,wk_ffnl1,mqgrid,npw)
482 : ! Filter for dilatmx>1
483 12892062 : if (present(kinpw)) then
484 351734233 : do ig=1,npw
485 351734233 : if(kinpw(ig)>huge(zero)*1.d-11)then
486 11624266 : wk_ffnl1(ig) = zero
487 11624266 : wk_ffnl2(ig) = zero
488 : end if
489 : end do
490 : end if
491 12892062 : if (ider == 2) then
492 247620 : call splfit(qgrid,wk_ffnl3,wk_ffspl,ider,kpgnorm,wk_ffnl1,mqgrid,npw)
493 247620 : if (present(kinpw)) then
494 0 : do ig=1,npw
495 0 : if(kinpw(ig)>huge(zero)*1.d-11)then
496 0 : wk_ffnl3(ig) = zero
497 : end if
498 : end do
499 : end if
500 : end if
501 : end if
502 :
503 : ! Store FFNL and FFNL derivatives
504 : ! -------------------------------
505 :
506 : ! =========================================================================
507 : ! A-USE OF SPHER. HARMONICS IN APPLICATION OF NL OPERATOR:
508 : ! ffnl(K,l,m,n)=fnl(K).Ylm(K)
509 : ! --if (idir==0)
510 : ! ffnl_prime(K,1:3,l,m,n)=3 reduced coordinates of d(ffnl)/dK^cart
511 : ! =fnl_prime(K).Ylm(K).K^red_i/|K|+fnl(K).(dYlm/dK^cart)^red_i
512 : ! --if (0<idir<4)
513 : ! ffnl_prime(K,l,m,n)=cart. coordinate idir of d(ffnl)/dK^red
514 : ! --if (idir==4)
515 : ! ffnl_prime(K,l,m,n)=3 cart. coordinates of d(ffnl)/dK^red
516 : ! --if (-7<=idir<0) - |idir|=(mu,nu) (1->11,2->22,3->33,4->32,5->31,6->21)
517 : ! ffnl_prime(K,l,m,n)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
518 : ! ffnl_prime_prime(K,l,m,n)=6 reduced coordinates of d2(ffnl)/dK^cart.dK^cart
519 :
520 20089372 : if (useylm==1) then
521 11785117 : iylm = il**2 + il + 1 + im
522 : !$OMP PARALLEL DO
523 1667765652 : do ig=1,npw
524 1667765652 : ffnl(ig,1,iffnl,itypat)=ylm(ig,iylm)*wk_ffnl1(ig)
525 : end do
526 :
527 11785117 : if (ider>=1) then
528 : !$OMP PARALLEL DO COLLAPSE(2)
529 9292140 : do mu=1,3
530 982556583 : do ig=1,npw
531 980233548 : dffnl_red(ig,mu)=ylm(ig,iylm)*wk_ffnl2(ig)*kpgn(ig,mu)+ylm_gr(ig,mu,iylm)*wk_ffnl1(ig)
532 : end do
533 : end do
534 : ! Special cases |k+g|=0
535 2323035 : if (ig0>0) then
536 58752 : do mu=1,3
537 44064 : dffnl_red(ig0,mu)=zero
538 58752 : if (il==1) then
539 : !Retrieve 1st-deriv. of ffnl at q=zero according to spline routine
540 21825 : yp1=(wk_ffspl(2,1)-wk_ffspl(1,1))/qgrid(2)-sixth*qgrid(2)*(two*wk_ffspl(1,2)+wk_ffspl(2,2))
541 21825 : fact=yp1*sqrt(three/four_pi)
542 21825 : if (im==-1) dffnl_red(ig0,mu)=fact*rprimd(2,mu)
543 21825 : if (im== 0) dffnl_red(ig0,mu)=fact*rprimd(3,mu)
544 21825 : if (im==+1) dffnl_red(ig0,mu)=fact*rprimd(1,mu)
545 : end if
546 : end do
547 : end if
548 2323035 : if (idir==0) then
549 : !$OMP PARALLEL DO COLLAPSE(2)
550 1053512 : do mu=1,3
551 150553619 : do ig=1,npw
552 150290241 : ffnl(ig,1+mu,iffnl,itypat)=dffnl_red(ig,mu)
553 : end do
554 : end do
555 : else
556 832002964 : dffnl_cart=zero
557 : !$OMP PARALLEL DO COLLAPSE(2)
558 8238628 : do mu=1,3
559 832002964 : do ig=1,npw
560 3301236315 : do nu=1,3
561 3295057344 : dffnl_cart(ig,mu)=dffnl_cart(ig,mu)+dffnl_red(ig,nu)*gprimd(mu,nu)
562 : end do
563 : end do
564 : end do
565 2059657 : if (idir>=1.and.idir<=3) then
566 158985753 : dffnl_tmp=zero
567 : !$OMP PARALLEL PRIVATE(nu,ig)
568 : !$OMP DO
569 158985753 : do ig=1,npw
570 632578167 : do nu=1,3
571 631456552 : dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,idir)
572 : end do
573 : end do
574 : !$OMP END DO
575 : !$OMP WORKSHARE
576 158985753 : ffnl(:,2,iffnl,itypat)=dffnl_tmp(:)
577 : !$OMP END WORKSHARE
578 : !$OMP END PARALLEL
579 938042 : else if (idir==4) then
580 2658008 : do mu=1,3
581 : !$OMP PARALLEL PRIVATE(nu,ig)
582 : !$OMP WORKSHARE
583 215744808 : dffnl_tmp=zero
584 : !$OMP END WORKSHARE
585 : !$OMP DO
586 215744808 : do ig=1,npw
587 856998714 : do nu=1,3
588 855005208 : dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,mu)
589 : end do
590 : end do
591 : !$OMP END DO
592 : !$OMP WORKSHARE
593 216409310 : ffnl(:,1+mu,iffnl,itypat)=dffnl_tmp(:)
594 : !$OMP END WORKSHARE
595 : !$OMP END PARALLEL
596 : end do
597 273540 : else if (idir/=-7) then
598 223488 : mu=abs(idir);mua=alpha(mu);mub=beta(mu)
599 : !$OMP PARALLEL DO
600 37279888 : do ig=1,npw
601 37279888 : ffnl(ig,2,iffnl,itypat)=0.5d0* (dffnl_cart(ig,mua)*kpgc(ig,mub) + dffnl_cart(ig,mub)*kpgc(ig,mua))
602 : end do
603 : else if (idir==-7) then
604 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mua, mub)
605 350364 : do mu=1,6
606 50853204 : do ig=1,npw
607 50502840 : mua=alpha(mu);mub=beta(mu)
608 50803152 : ffnl(ig,1+mu,iffnl,itypat)=0.5d0 * (dffnl_cart(ig,mua)*kpgc(ig,mub) + dffnl_cart(ig,mub)*kpgc(ig,mua))
609 : end do
610 : end do
611 : end if
612 : end if
613 : end if
614 :
615 11785117 : if (ider==2) then
616 3468024 : do mu=1,6
617 2972592 : mua=alpha(mu);mub=beta(mu)
618 2972592 : rmetab=rmet(mua,mub)
619 : !$OMP PARALLEL DO
620 286909602 : do ig=1,npw
621 : d2ffnl_red(ig,mu)= &
622 : ylm_gr(ig,3+mu,iylm)*wk_ffnl1(ig) &
623 : + (rmetab-kpgn(ig,mua)*kpgn(ig,mub))*ylm(ig,iylm)*wk_ffnl2(ig)*kpgnorm_inv(ig) &
624 : + ylm(ig,iylm)*kpgn(ig,mua)*kpgn(ig,mub)*wk_ffnl3(ig) &
625 286909602 : + (ylm_gr(ig,mua,iylm)*kpgn(ig,mub)+ylm_gr(ig,mub,iylm)*kpgn(ig,mua))*wk_ffnl2(ig)
626 : end do
627 : ! Special cases |k+g|=0
628 3468024 : if (ig0>0) then
629 1452 : d2ffnl_red(ig0,mu)=zero
630 1452 : if (il==0) then
631 276 : d2ffnl_red(ig0,mu)=wk_ffspl(1,2)*rmetab/sqrt(four_pi)
632 : end if
633 1452 : if (il==2) then
634 600 : fact=wk_ffspl(1,2)*quarter*sqrt(15._dp/pi)
635 600 : if (im==-2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(2,mub)+rprimd(2,mua)*rprimd(1,mub))
636 600 : if (im==-1) d2ffnl_red(ig0,mu)=fact*(rprimd(2,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(2,mub))
637 600 : if (im==+1) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(1,mub))
638 600 : if (im==+2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
639 600 : if (im== 0) d2ffnl_red(ig0,mu)=(fact/sqrt3)*(two*rprimd(3,mua)*rprimd(3,mub) &
640 120 : -rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
641 : end if
642 : end if
643 : end do
644 495432 : if (idir==0) then
645 : !$OMP PARALLEL DO COLLAPSE(2)
646 113344 : do mu=1,6
647 11147056 : do ig=1,npw
648 11130864 : ffnl(ig,4+mu,iffnl,itypat)=d2ffnl_red(ig,mu)
649 : end do
650 : end do
651 479240 : else if (idir==4) then
652 276257978 : d2ffnl_cart=zero
653 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mu,mua,mub,ig,nu,nua,nub)
654 3354680 : do mu=1,6
655 276257978 : do ig=1,npw
656 272903298 : mua=alpha(mu);mub=beta(mu)
657 1094488632 : do nua=1,3
658 3547742874 : do nub=1,3
659 2456129682 : nu=gamma(nua,nub)
660 3274839576 : d2ffnl_cart(ig,mu)=d2ffnl_cart(ig,mu)+d2ffnl_red(ig,nu)*gprimd(mua,nua)*gprimd(mub,nub)
661 : end do
662 : end do
663 : end do
664 : end do
665 3354680 : do mu=1,6
666 2875440 : mua=alpha(mu);mub=beta(mu)
667 : !$OMP PARALLEL PRIVATE(nu,nua,nub,ig)
668 : !$OMP WORKSHARE
669 275778738 : d2ffnl_tmp=zero
670 : !$OMP END WORKSHARE
671 : !$OMP DO
672 275778738 : do ig=1,npw
673 1094488632 : do nua=1,3
674 3547742874 : do nub=1,3
675 2456129682 : nu=gamma(nua,nub)
676 3274839576 : d2ffnl_tmp(ig)=d2ffnl_tmp(ig)+d2ffnl_cart(ig,nu)*gprimd(nua,mua)*gprimd(nub,mub)
677 : end do
678 : end do
679 : end do
680 : !$OMP END DO
681 : !$OMP WORKSHARE
682 276257978 : ffnl(:,4+mu,iffnl,itypat)=d2ffnl_tmp(:)
683 : !$OMP END WORKSHARE
684 : !$OMP END PARALLEL
685 : end do
686 : end if
687 : end if
688 :
689 : ! =========================================================================
690 : ! B-USE OF LEGENDRE POLYNOMIAL IN APPLICATION OF NL OPERATOR:
691 : ! ffnl(K,l,n)=fnl(K)/|K|^l
692 : ! ffnl_prime(K,l,n)=(fnl_prime(K)-l*fnl(K)/|K|)/|K|^(l+1)
693 : ! ffnl_prime_prime(K,l,n)=(fnl_prime_prime(K)-(2*l+1)*fnl_prime(K)/|K|
694 : ! +l*(l+2)*fnl(K)/|K|^2)/|K|^(l+2)
695 8304255 : else if (iln>iln0) then
696 :
697 8304255 : if (il==0) then
698 : !$OMP PARALLEL DO
699 988533204 : do ig=1,npw
700 988533204 : ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)
701 : end do
702 : else
703 : !$OMP PARALLEL DO
704 1273196563 : do ig=1,npw
705 1273196563 : ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)*kpgnorm_inv(ig)**il
706 : end do
707 : end if
708 8304255 : if (ider>=1) then
709 : !$OMP PARALLEL DO
710 185045978 : do ig=1,npw
711 185045978 : ffnl(ig,2,iffnl,itypat)= (wk_ffnl2(ig)-dble(il)*wk_ffnl1(ig)*kpgnorm_inv(ig))*kpgnorm_inv(ig)**(il+1)
712 : end do
713 1743469 : if (ider==2) then
714 : !$OMP PARALLEL DO
715 4322143 : do ig=1,npw
716 : ffnl(ig,3,iffnl,itypat)= (wk_ffnl3(ig)- &
717 : dble(2*il+1)*wk_ffnl2(ig)*kpgnorm_inv(ig)+ &
718 4322143 : dble(il*(il+2))*wk_ffnl1(ig)*kpgnorm_inv(ig)**2)*kpgnorm_inv(ig)**(il+2)
719 : end do
720 : end if
721 : end if
722 :
723 : end if ! Use of Ylm or not
724 :
725 : else
726 : ! No NL part
727 : !$OMP PARALLEL DO COLLAPSE(2)
728 1541 : do mu=1,dimffnl
729 388740 : do ig=1,npw
730 330124 : ffnl(ig,mu,iffnl,itypat)=zero
731 : end do
732 : end do
733 :
734 : end if ! testnl (a nonlocal part exists)
735 : end if ! special case: spin orbit calc. & no spin-orbit psp
736 :
737 23968165 : if (iln > iln0) iln0 = iln
738 :
739 : end do ! loop over (l,m,n) values
740 : end do ! loop over atom types
741 :
742 2370253 : ABI_FREE(kpgnorm_inv)
743 2370253 : ABI_FREE(kpgnorm)
744 2370253 : ABI_FREE(wk_ffnl1)
745 2370253 : ABI_FREE(wk_ffnl2)
746 2370253 : ABI_FREE(wk_ffnl3)
747 2370253 : ABI_FREE(wk_ffspl)
748 :
749 : ! Optional deallocations.
750 2370253 : ABI_SFREE(kpgc)
751 2370253 : ABI_SFREE(kpgn)
752 2370253 : ABI_SFREE(dffnl_red)
753 2370253 : ABI_SFREE(d2ffnl_red)
754 2370253 : ABI_SFREE(dffnl_cart)
755 2370253 : ABI_SFREE(d2ffnl_cart)
756 2370253 : ABI_SFREE(dffnl_tmp)
757 2370253 : ABI_SFREE(d2ffnl_tmp)
758 :
759 2370253 : if (nprocs > 1) then
760 : ! Blocking/non-blocking depending on the presence of request.
761 0 : if (present(request)) then
762 0 : call xmpi_isum_ip(ffnl, comm, request, ierr)
763 : else
764 0 : call xmpi_sum(ffnl, comm, ierr)
765 : end if
766 : else
767 2370253 : if (present(request)) request = xmpi_request_null
768 : end if
769 :
770 2370253 : call timab(16, 2, tsec)
771 :
772 2370253 : end subroutine mkffnl
773 : !!***
774 :
775 : end module m_mkffnl
776 : !!***
|
