Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_sphharm
2 : : !! NAME
3 : : !! m_paw_sphharm
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains a set of routines to compute the complex (resp. real)
7 : : !! spherical harmonics Ylm (resp. Slm) (and gradients).
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2013-2026 ABINIT group (MT, FJ, NH, TRangel)
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 : : !! NOTES
16 : : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
17 : : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
18 : : !!
19 : : !! SOURCE
20 : :
21 : : #include "libpaw.h"
22 : :
23 : : #define MATPACK(row,col) (MAX(row,col)*(MAX(row,col)-1)/2 + MIN(row,col))
24 : : #define LMPACK(lqn,mqn) (lqn*lqn+lqn+mqn+1)
25 : :
26 : : MODULE m_paw_sphharm
27 : :
28 : : USE_DEFS
29 : : USE_MSG_HANDLING
30 : : USE_MEMORY_PROFILING
31 : :
32 : : implicit none
33 : :
34 : : private
35 : :
36 : : !Public procedures.
37 : : public :: ylmc ! Complex Spherical harmonics for l<=3.
38 : : public :: ylmcd ! First derivative of complex Ylm wrt theta and phi up to l<=3
39 : : public :: ylm_cmplx ! All (complex) spherical harmonics for lx<=4
40 : : public :: initylmr ! Real Spherical Harmonics on a set of vectors
41 : : public :: ys ! Matrix element <Yl'm'|Slm>
42 : : public :: lxyz ! Matrix element <Yl'm'|L_idir|Ylm>
43 : : public :: slxyzs ! Matrix element <Sl'm'|L_idir|Slm>
44 : : public :: lsylm ! Compute the LS operator in the real spherical harmonics basis
45 : : public :: plm_coeff ! Coefficients depending on Plm used to compute the 2nd der of Ylm
46 : : public :: ass_leg_pol ! Associated Legendre Polynomial Plm(x)
47 : : public :: plm_dphi ! m*P_lm(x)/sqrt((1-x^2) (P_lm= associatedLegendre polynomial)
48 : : public :: plm_dtheta ! -(1-x^2)^1/2*d/dx{P_lm(x)} (P_lm= associated Legendre polynomial)
49 : : public :: plm_d2theta ! d2(Plm (cos(theta)))/d(theta)2 (P_lm= associated Legendre polynomial)
50 : : public :: pl_deriv ! d2(Pl (x)))/d(x)2 where P_l is a Legendre polynomial
51 : : public :: ylm_angular_mesh ! Build (theta, phi) angular mesh
52 : : public :: mat_mlms2jmj ! Change a matrix from the Ylm basis to the J,M_J basis
53 : : public :: mat_slm2ylm ! Change a matrix from the Slm to the Ylm basis or from Ylm to Slm
54 : : public :: setsym_ylm ! Compute rotation matrices expressed in the basis of real spherical harmonics
55 : : public :: setnabla_ylm ! Evaluate several integrals involving spherical harmonics and their gradient
56 : : public :: gaunt ! Gaunt coeffients for complex Yml
57 : : public :: realgaunt ! Compute "real Gaunt coefficients" with "real spherical harmonics"
58 : : public :: nablarealgaunt ! Compute the integrals Grad(Slimi).Grad(Sjmj) Slkmk of real spherical harmonics
59 : : public :: make_dyadic ! Compute integrals of S_limi(alpha*1 - \beta*\hat{r}\hat{r})S_ljmj
60 : :
61 : : !Private functions
62 : : private :: create_slm2ylm ! For a given angular momentum lcor, compute slm2ylm
63 : : private :: create_mlms2jmj ! For a given angular momentum lcor, give the rotation matrix msml2jmj
64 : : private :: mkeuler ! For a given symmetry operation, determines the corresponding Euler angles
65 : : private :: dbeta ! Calculate the rotation matrix d^l_{m{\prim}m}(beta)
66 : : private :: phim ! Computes Phi_m[t]=Sqrt2.cos[m.t] (m>0), Sqrt2.sin[|m|.t] (m<0), 1 (m=0)
67 : : private :: gauleg ! Compute the coefficients for Gauss-Legendre integration
68 : : private :: perms ! Returns N!/(N-k)! if N>=0 and N>k
69 : : private :: rfactorial ! Calculates N! as a double precision real
70 : : !!***
71 : :
72 : : CONTAINS
73 : :
74 : : !===========================================================
75 : : !!***
76 : :
77 : : !!****f* m_paw_sphharm/ylmc
78 : : !! NAME
79 : : !! ylmc
80 : : !!
81 : : !! FUNCTION
82 : : !! Return a complex spherical harmonic with l <= 3
83 : : !!
84 : : !! INPUTS
85 : : !! il=angular quantum number
86 : : !! im=magnetic quantum number
87 : : !! kcart=vector in cartesian coordinates defining the value of \theta and \psi
88 : : !! where calculate the spherical harmonic
89 : : !!
90 : : !! OUTPUT
91 : : !! ylm= spherical harmonic
92 : : !!
93 : : !! NOTES
94 : : !! Note the use of double precision complex.
95 : : !! Case l>3 not implemented.
96 : : !!
97 : : !! SOURCE
98 : :
99 : 191291224 : function ylmc(il,im,kcart)
100 : :
101 : : !Arguments ------------------------------------
102 : : !scalars
103 : : integer,intent(in) :: il,im
104 : : complex(dp) :: ylmc
105 : : !arrays
106 : : real(dp),intent(in) :: kcart(3)
107 : :
108 : : !Local variables-------------------------------
109 : : !scalars
110 : : integer,parameter :: LMAX=3
111 : : real(dp),parameter :: PPAD=tol8
112 : : real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi
113 : : !complex(dp) :: new_ylmc
114 : : character(len=500) :: msg
115 : : complex(dp) :: ctmp
116 : : ! *************************************************************************
117 : :
118 [ - + ]: 191291224 : if (ABS(im)>ABS(il)) then
119 : 0 : write(msg,'(3(a,i0))') 'm is,',im,' however it should be between ',-il,' and ',il
120 : 0 : LIBPAW_ERROR(msg)
121 : : end if
122 : :
123 : 191291224 : ylmc = czero
124 : :
125 : 191291224 : r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
126 [ - + ]: 191291224 : if (r<PPAD) r=r+PPAD
127 : : !$if (r<tol10) RETURN
128 : :
129 : 191291224 : rxy=SQRT(kcart(1)**2+kcart(2)**2)
130 [ - + ]: 191291224 : if (rxy<PPAD)rxy=r+PPAD
131 : : !
132 : : ! Determine theta and phi
133 : 191291224 : costh= kcart(3)/r
134 : :
135 : : #if 1
136 : : ! old buggy coding
137 : 191291224 : sinth= rxy/r
138 : 191291224 : cosphi= kcart(1)/rxy
139 : 191291224 : sinphi= kcart(2)/rxy
140 : : #else
141 : : sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg
142 : : cosphi=one
143 : : sinphi=zero
144 : : if (sinth>tol10) then
145 : : cosphi=kcart(1)/(r*sinth)
146 : : sinphi=kcart(2)/(r*sinth)
147 : : end if
148 : : #endif
149 : :
150 : 191291224 : costwophi= two*cosphi**2 - one
151 : 191291224 : sintwophi= two*sinphi*cosphi
152 : 191291224 : costhreephi=cosphi*costwophi-sinphi*sintwophi
153 : 191291224 : sinthreephi=cosphi*sintwophi+sinphi*costwophi
154 : :
155 : : select case (il)
156 : :
157 : : case (0)
158 : 271666888 : ylmc= one/SQRT(four_pi)
159 : :
160 : : case (1)
161 [ + + ]: 80375664 : if (ABS(im)==0) then
162 : 26791888 : ylmc = SQRT(three/(four_pi))*costh
163 [ + - ]: 53583776 : else if (ABS(im)==1) then
164 : 53583776 : ylmc = -SQRT(three/(8._dp*pi))*sinth*CMPLX(cosphi,sinphi)
165 : : else
166 : 0 : msg='wrong im'
167 : 0 : LIBPAW_ERROR(msg)
168 : : end if
169 : :
170 : : case (2)
171 [ + + ]: 83963160 : if (ABS(im)==0) then
172 : 16792632 : ylmc = SQRT(5.d0/(16.d0*pi))*(three*costh**2-one)
173 [ + + ]: 67170528 : else if (ABS(im)==1) then
174 : 33585264 : ylmc = -SQRT(15.d0/(8.d0*pi))*sinth*costh*cmplx(cosphi,sinphi)
175 [ + - ]: 33585264 : else if (ABS(im)==2) then
176 : 33585264 : ylmc = SQRT(15.d0/(32.d0*pi))*(sinth)**2*CMPLX(costwophi,sintwophi)
177 : : else
178 : 0 : msg='wrong im'
179 : 0 : LIBPAW_ERROR(msg)
180 : : end if
181 : :
182 : : case (3)
183 [ - + + - : 191291224 : if (ABS(im)==0) then
+ - - ]
184 : 0 : ylmc= SQRT(7.d0/(16.d0*pi))*(5.d0*costh**3 -3.d0*costh)
185 [ # # ]: 0 : else if (ABS(im)==1) then
186 : 0 : ylmc= -SQRT(21.d0/(64.d0*pi))*sinth*(5.d0*costh**2-one)*CMPLX(cosphi,sinphi)
187 [ # # ]: 0 : else if (ABS(im)==2) then
188 : 0 : ylmc= SQRT(105.d0/(32.d0*pi))*sinth**2*costh*CMPLX(costwophi,sintwophi)
189 [ # # ]: 0 : else if (ABS(im)==3) then
190 : 0 : ylmc=-SQRT(35.d0/(64.d0*pi))*sinth**3*CMPLX(costhreephi,sinthreephi)
191 : : else
192 : 0 : msg='wrong im'
193 : 0 : LIBPAW_ERROR(msg)
194 : : end if
195 : :
196 : : case default
197 : : !write(msg,'(a,i6,a,i6)')' The maximum allowed value for l is,',LMAX,' however l=',il
198 : : !LIBPAW_ERROR(msg)
199 : : end select
200 : : !
201 : : !=== Treat the case im < 0 ===
202 [ + + ]: 191291224 : if (im < 0) then
203 : 60377152 : ctmp = (-one)**(im)*CONJG(ylmc)
204 : 60377152 : ylmc = ctmp
205 : : end if
206 : :
207 : : ! FIXME: Use the piece of code below as it works for arbitrary (l,m)
208 : : ! the implementation above is buggy when the vector is along z!
209 : : !
210 : : #if 0
211 : : ! Remember the expression of complex spherical harmonics:
212 : : ! $Y_{lm}(\theta,\phi)=sqrt{{(2l+1) over (4\pi)} {fact(l-m)/fact(l+m)} } P_l^m(cos(\theta)) e^{i m\phi}$
213 : : new_ylmc = SQRT((2*il+1)*rfactorial(il-ABS(im))/(rfactorial(il+ABS(im))*four_pi)) * &
214 : : & ass_leg_pol(il,ABS(im),costh) * CMPLX(cosphi,sinphi)**ABS(im)
215 : : if (im<0) new_ylmc=(-one)**(im)*CONJG(new_ylmc)
216 : :
217 : : if (ABS(new_ylmc-ylmc)>tol6) then
218 : : !LIBPAW_WARNING("Check new_ylmc")
219 : : !write(std_out,*)"il,im,new_ylmc, ylmc",il,im,new_ylmc,ylmc
220 : : !write(std_out,*)"fact",SQRT((2*il+1)*rfactorial(il-ABS(im))/(rfactorial(il+ABS(im))*four_pi))
221 : : !write(std_out,*)"costh,sinth,ass_leg_pol",costh,sinth,ass_leg_pol(il,ABS(im),costh)
222 : : !write(std_out,*)"cosphi,sinphi,e^{imphi}",cosphi,sinphi,CMPLX(cosphi,sinphi)**ABS(im)
223 : : end if
224 : : ylmc = new_ylmc
225 : : #endif
226 : :
227 : 191291224 : end function ylmc
228 : : !!***
229 : :
230 : : !----------------------------------------------------------------------
231 : :
232 : : !!****f* m_paw_sphharm/ylmcd
233 : : !! NAME
234 : : !! ylmcd
235 : : !!
236 : : !! FUNCTION
237 : : !! Computes dth and dphi, the first derivatives of complex Ylm as a function of
238 : : !! th and phi (the angles of the spherical coordinates)
239 : : !! It works for all spherical harmonics with l <= 3
240 : : !!
241 : : !! INPUTS
242 : : !! il=angular quantum number
243 : : !! im=magnetic quantum number
244 : : !! kcart=cartesian coordinates of the vector where the first derivatives of Ylm are evaluated
245 : : !!
246 : : !! OUTPUT
247 : : !! dth =derivative of Y_lm with respect to \theta
248 : : !! dphi=derivative of Y_lm with respect to \phi
249 : : !!
250 : : !! NOTES
251 : : !! Note the use of double precision complex.
252 : : !! Case l>3 not implemented.
253 : : !!
254 : : !! SOURCE
255 : :
256 : 47822806 : subroutine ylmcd(il,im,kcart,dth,dphi)
257 : :
258 : : !Arguments ------------------------------------
259 : : !scalars
260 : : integer,intent(in) :: il,im
261 : : complex(dp),intent(out) :: dphi,dth
262 : : !arrays
263 : : real(dp),intent(in) :: kcart(3)
264 : :
265 : : !Local variables-------------------------------
266 : : !scalars
267 : : integer,parameter :: LMAX=3
268 : : real(dp),parameter :: PPAD=tol8
269 : : real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi,c
270 : : character(len=500) :: msg
271 : : complex(dp) :: ctmp
272 : : ! *************************************************************************
273 : :
274 [ - + ]: 47822806 : if (ABS(im)>ABS(il))then
275 : 0 : write(msg,'(3(a,i0))')' m is,',im,' however it should be between ',-il,' and ',il
276 : 0 : LIBPAW_ERROR(msg)
277 : : end if
278 : :
279 : 47822806 : dphi=czero; dth=czero
280 : :
281 : 47822806 : r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
282 [ - + ]: 47822806 : if (r<PPAD) r=r+PPAD
283 : : !$if (r<tol10) RETURN
284 : :
285 : 47822806 : rxy=SQRT(kcart(1)**2+kcart(2)**2)
286 [ - + ]: 47822806 : if (rxy<PPAD) rxy=r+PPAD
287 : :
288 : : ! Determine theta and phi
289 : 47822806 : costh= kcart(3)/r
290 : : #if 1
291 : : ! old buggy coding
292 : 47822806 : sinth= rxy/r
293 : 47822806 : cosphi= kcart(1)/rxy
294 : 47822806 : sinphi= kcart(2)/rxy
295 : : #else
296 : : sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg
297 : : cosphi=one
298 : : sinphi=zero
299 : : if (sinth>tol10) then
300 : : cosphi=kcart(1)/(r*sinth)
301 : : sinphi=kcart(2)/(r*sinth)
302 : : end if
303 : : #endif
304 : :
305 : 47822806 : costwophi= two*cosphi**2 - one
306 : 47822806 : sintwophi= two*sinphi*cosphi
307 : 47822806 : costhreephi=cosphi*costwophi-sinphi*sintwophi
308 : 47822806 : sinthreephi=cosphi*sintwophi+sinphi*costwophi
309 : :
310 : : select case (il)
311 : :
312 : : case (0)
313 : : dth = czero
314 : 20093916 : dphi = czero
315 : :
316 : : case (1)
317 [ + + ]: 20093916 : if (ABS(im)==0) then
318 : 6697972 : dth= -SQRT(three/(four_pi))*sinth
319 : : dphi= czero
320 [ + - ]: 13395944 : else if (abs(im)==1) then
321 : 13395944 : dth= -SQRT(3.d0/(8.d0*pi))*costh*CMPLX(cosphi,sinphi)
322 : 13395944 : dphi=-SQRT(3.d0/(8.d0*pi))*sinth*CMPLX(-sinphi,cosphi)
323 : : end if
324 : :
325 : : case (2)
326 [ + + ]: 20990790 : if (ABS(im)==0) then
327 : 4198158 : dth= -SQRT(5.d0/(16.d0*pi))*6.d0*costh*sinth
328 : : dphi= czero
329 [ + + ]: 16792632 : else if (ABS(im)==1) then
330 : 8396316 : dth= -SQRT(15.d0/(8.d0*pi))*(costh**2-sinth**2)*CMPLX(cosphi,sinphi)
331 : 8396316 : dphi= -SQRT(15.d0/(8.d0*pi))*costh*sinth*(0.d0,1.d0)*CMPLX(cosphi,sinphi)
332 [ + - ]: 8396316 : else if (abs(im)==2) then
333 : 8396316 : dth = SQRT(15.d0/(32.d0*pi))*2.d0*costh*sinth*CMPLX(costwophi,sintwophi)
334 : 8396316 : dphi = SQRT(15.d0/(32.d0*pi))*sinth**2*(0.d0,2.d0)*CMPLX(costwophi,sintwophi)
335 : : end if
336 : :
337 : : case (3)
338 [ # # ]: 0 : if (ABS(im)==0) then
339 : 0 : dth = SQRT(7.d0/(16*pi))*(-15.d0*costh**2*sinth + 3.d0**sinth)
340 : : dphi= czero
341 [ # # ]: 0 : else if (ABS(im)==1) then
342 : 0 : c = SQRT(21.d0/(64.d0*pi))
343 : 0 : dth= -c* (15.d0*costh**3-11.d0*costh)* CMPLX(cosphi,sinphi)
344 : 0 : dphi=-c*sinth*( 5.d0*costh**2-1 )*(0.d0,1.d0)*CMPLX(cosphi,sinphi)
345 [ # # ]: 0 : else if (ABS(im)==2) then
346 : 0 : c = SQRT(105.d0/(32.d0*pi))
347 : 0 : dth =c*(2.d0*sinth*costh**2-sinth**3) *CMPLX(costwophi,sintwophi)
348 : 0 : dphi=c*(2.d0*sinth**2*costh)*(0.d0,1.d0)*CMPLX(costwophi,sintwophi)
349 [ # # ]: 0 : else if (abs(im)==3) then
350 : 0 : dth =-SQRT(35.d0/(64.d0*pi))*3.d0*sinth**2*costh*CMPLX(costhreephi,sinthreephi)
351 : 0 : dphi=-SQRT(35.d0/(64.d0*pi))*sinth**3*(0.d0,3.d0)*CMPLX(costhreephi,sinthreephi)
352 : : end if
353 : :
354 : : case default
355 : 0 : write(msg,'(2(a,i0))')' The maximum allowed value for l is,',LMAX,' however, l=',il
356 [ + + - - : 47822806 : LIBPAW_ERROR(msg)
+ ]
357 : : end select
358 : : !
359 : : !=== Treat the case im < 0 ===
360 [ + + ]: 47822806 : if (im<0) then
361 : 15094288 : ctmp = (-one)**(im)*CONJG(dth)
362 : 15094288 : dth = ctmp
363 : 15094288 : ctmp= (-one)**(im)*CONJG(dphi)
364 : 15094288 : dphi= ctmp
365 : : end if
366 : :
367 : 47822806 : end subroutine ylmcd
368 : : !!***
369 : :
370 : : !----------------------------------------------------------------------
371 : :
372 : : !!****f* m_paw_sphharm/ylm_cmplx
373 : : !! NAME
374 : : !! ylm_cmplx
375 : : !!
376 : : !! FUNCTION
377 : : !! Calculate all (complex) spherical harmonics for lx<=4
378 : : !!
379 : : !! INPUTS
380 : : !! lx= quantum numbers.
381 : : !! xx= cartesian coordinate in the x direction
382 : : !! yy= cartesian coordinate in the y direction
383 : : !! zz= cartesian coordinate in the z direction
384 : : !!
385 : : !! cartesian coordinates
386 : : !! OUTPUT
387 : : !! ylm((lx+1)*(lx+1)) complex spherical harmonics for all l<=lx and all
388 : : !! possible values of m.
389 : : !!
390 : : !! NOTES
391 : : !! We are supressing the so-called Condon-Shortley phase
392 : : !!
393 : : !! SOURCE
394 : :
395 : 134873 : subroutine ylm_cmplx(lx,ylm,xx,yy,zz)
396 : :
397 : : !Arguments ------------------------------------
398 : : !scalars
399 : : integer,intent(in) :: lx
400 : : real(dp),intent(in) :: xx,yy,zz
401 : : !arrays
402 : : complex(dp),intent(out) :: ylm((lx+1)*(lx+1))
403 : :
404 : : !Local variables-------------------------------
405 : : !scalars
406 : : integer :: ii,l1,m1,nc,nn
407 : : real(dp) :: dc,dl,dm,ds,rr,rrs,rs,sq2,w,x,xs,ya,yi,yr
408 : : !arrays
409 : 269746 : real(dp) :: cosa(lx+1),fact(2*(lx+1)),plm(lx+2,lx+2),qlm(lx+2,lx+2),sgn(lx+1)
410 : 269746 : real(dp) :: sina(lx+1)
411 : : ! *************************************************************************
412 : :
413 : : !normalization coefficients
414 : 134873 : sq2=sqrt(2.0d0)
415 : 134873 : fact(1)=1.0d0
416 [ + + ]: 475363 : do ii=2,2*lx+1
417 : 475363 : fact(ii)=(ii-1)*fact(ii-1)
418 : : end do
419 [ + + ]: 439991 : do l1=1,lx+1
420 : 305118 : sgn(l1)=(-1.d0)**(l1-1)
421 [ + + ]: 950726 : do m1=1,l1
422 : : qlm(l1,m1)=sqrt((2*l1-1)*fact(l1-m1+1)/&
423 : 815853 : & (four_pi*fact(l1+m1-1)))
424 : : end do
425 : : end do
426 : :
427 : : !legendre polynomials
428 : 134873 : rs=xx**2 + yy**2 + zz**2
429 [ + + ]: 134873 : if(rs > tol8) then
430 : 134861 : xs=zz**2/rs
431 : 134861 : x=zz/sqrt(rs)
432 : 134861 : w=sqrt(abs(1.0d0 - xs))
433 : : else
434 : : x=0.0d0
435 : :
436 : : w=1.0d0
437 : : end if
438 : 134873 : plm(1,1)=1.0d0
439 : 134873 : plm(2,1)=x
440 : 134873 : plm(2,2)=w
441 : 134873 : plm(3,2)=3.0d0*x*w
442 [ + + ]: 305118 : do m1=1,lx
443 : 170245 : dm=m1-1
444 [ + + ]: 170245 : if(m1 > 1) then
445 : 35372 : plm(m1+1,m1)=x*plm(m1,m1) + 2*dm*w*plm(m1,m1-1)
446 : : end if
447 [ + + ]: 170245 : if(m1 < lx) then
448 [ + + ]: 70744 : do l1=m1+2,lx+1
449 : 35372 : dl=l1-1
450 : : plm(l1,m1)=((2*dl-1)*x*plm(l1-1,m1)&
451 : 70744 : & - (dl+dm-1)*plm(l1-2,m1))/(dl-dm)
452 : : end do
453 : : end if
454 : 305118 : plm(m1+1,m1+1)=(2*dm+1)*w*plm(m1,m1)
455 : : end do
456 : :
457 : : !azimuthal angle phase factors
458 : 134873 : rrs=xx**2 + yy**2
459 [ + + ]: 134873 : if(rrs > tol8) then
460 : 134623 : rr=sqrt(rrs)
461 : 134623 : dc=xx/rr
462 : 134623 : ds=yy/rr
463 : : else
464 : : dc=1.0d0
465 : : ds=0.0d0
466 : : end if
467 : 134873 : cosa(1)=1.0d0
468 : 134873 : sina(1)=0.0d0
469 [ + + ]: 305118 : do m1=2,lx+1
470 : 170245 : cosa(m1)=dc*cosa(m1-1) - ds*sina(m1-1)
471 : 305118 : sina(m1)=ds*cosa(m1-1) + dc*sina(m1-1)
472 : : end do
473 : :
474 : : !combine factors
475 [ + + ]: 439991 : do l1=1,lx+1
476 [ + + ]: 645608 : do m1=2,l1
477 : 205617 : nn=(l1-1)**2 + (l1-1) + (m1-1) + 1
478 : 205617 : nc=(l1-1)**2 + (l1-1) - (m1-1) + 1
479 : : ! note that we are supressing the so-called Condon-Shortley phase
480 : : ! ya=sgn(m1)*qlm(l1,m1)*plm(l1,m1)
481 : 205617 : ya=qlm(l1,m1)*plm(l1,m1)
482 : 205617 : yr=ya*cosa(m1)
483 : 205617 : yi=ya*sina(m1)
484 : 205617 : ylm(nc)=sgn(m1)*cmplx(yr,-yi)
485 : 510735 : ylm(nn)=cmplx(yr,yi)
486 : : end do
487 : : end do
488 [ + + ]: 439991 : do l1=1,lx+1
489 : 305118 : nn=(l1-1)**2 + (l1-1) + 1
490 : 305118 : ya=qlm(l1,1)*plm(l1,1)
491 : 439991 : ylm(nn)=cmplx(ya,0.d0)
492 : : end do
493 : :
494 : 134873 : end subroutine ylm_cmplx
495 : : !!***
496 : :
497 : : !----------------------------------------------------------------------
498 : :
499 : : !!****f* m_paw_sphharm/initylmr
500 : : !! NAME
501 : : !! initylmr
502 : : !!
503 : : !! FUNCTION
504 : : !! Calculate the real spherical harmonics Ylm (and gradients)
505 : : !! over a set of (r) vectors given in Cartesian coordinates.
506 : : !!
507 : : !! INPUTS
508 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotential
509 : : !! normchoice=0 the input rr vectors are normalized
510 : : !! =1 the norm of the input vector is in nrm() array
511 : : !! nrm(npts) = Depending of normchoice, this array contains
512 : : !! either the weight of the point or the norm of rr.
513 : : !! npts = number of rr vectors
514 : : !! option= 1=compute Ylm(R), 2=compute Ylm(R) and dYlm/dRi (cartesian derivatives),
515 : : !! 3=compute Ylm(R), dYlm/dRi and d2Ylm/dRidRj (cartesian derivatives)
516 : : !! rr(3,npts)= vectors for which ylmr have to be calculated
517 : : !! For each point of the spherical mesh, gives the
518 : : !! Cartesian coordinates of the corresponding point.
519 : : !!
520 : : !! OUTPUT
521 : : !! if (option=1, 2 or 3)
522 : : !! ylm(mpsang*mpsang,npts) = real spherical harmonics for each r point
523 : : !! if (option=2 or 3)
524 : : !! ylmr_gr(1:3,mpsang*mpsang,npts)= gradients of real spherical harmonics
525 : : !! if (option=3)
526 : : !! ylmr_gr(4:9,mpsang*mpsang,npts)= first and second gradients of real spherical harmonics
527 : : !!
528 : : !! NOTES
529 : : !! Remember the expression of complex spherical harmonics:
530 : : !! $Y_{lm}(%theta ,%phi)=sqrt{{(2l+1) over (4 %pi)} {fact(l-m) over fact(l+m)} } P_l^m(cos(%theta)) func e^{i m %phi}$
531 : : !! Remember the expression of real spherical harmonics as linear combination of imaginary spherical harmonics:
532 : : !! $Yr_{lm}(%theta ,%phi)=(Re{Y_{l-m}}+(-1)^m Re{Y_{lm}})/sqrt{2}
533 : : !! $Yr_{l-m}(%theta ,%phi)=(Im{Y_{l-m}}-(-1)^m Im{Y_{lm}})/sqrt{2}
534 : : !!
535 : : !! SOURCE
536 : :
537 : 435067 : subroutine initylmr(mpsang,normchoice,npts,nrm,option,rr,ylmr,ylmr_gr)
538 : :
539 : : !Arguments ------------------------------------
540 : : !scalars
541 : : integer,intent(in) :: mpsang,normchoice,npts,option
542 : : !arrays
543 : : real(dp),intent(in) :: nrm(npts),rr(3,npts)
544 : : real(dp),intent(out) :: ylmr(mpsang*mpsang,npts)
545 : : real(dp),optional,intent(out) :: ylmr_gr(3*(option/2)+6*(option/3),mpsang*mpsang,npts)
546 : :
547 : : !Local variables ------------------------------
548 : : !scalars
549 : : integer :: dimgr,ilang,inpt,l0,ll,mm
550 : : real(dp) :: cphi,ctheta,fact,onem,rnorm,sphi,stheta,work1,work2,ylmcst,ylmcst2
551 : : logical :: compute_ylm,compute_ylm2gr,compute_ylmgr
552 : : !arrays
553 : 940036 : real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),rphase(mpsang-1)
554 : 470018 : real(dp),allocatable :: blm(:,:)
555 : :
556 : : !************************************************************************
557 : :
558 : : !What has to be computed ?
559 : 470018 : compute_ylm = (option==1.or.option==2.or.option==3)
560 [ + + - + ]: 470018 : compute_ylmgr =(( option==2.or.option==3).and.present(ylmr_gr))
561 [ + + - + ]: 470018 : compute_ylm2gr=(( option==3).and.present(ylmr_gr))
562 : 470018 : dimgr=3*(option/2)+6*(option/3)
563 : :
564 : : !Initialisation of spherical harmonics
565 [ + - + + : 527248414 : if (compute_ylm ) ylmr (: ,1:npts)=zero
+ + ]
566 [ + + + + : 665806039 : if (compute_ylmgr) ylmr_gr(:,:,1:npts)=zero
+ + + + ]
567 : :
568 : : !Special case for l=0
569 [ + - + + ]: 42044984 : if (compute_ylm ) ylmr(1,1:npts)=1._dp/sqrt(four_pi)
570 [ + + + + : 45551226 : if (compute_ylmgr) ylmr_gr(1:dimgr,1,1:npts)=zero
+ + ]
571 [ + + ]: 470018 : if (mpsang>1) then
572 : :
573 : : ! Loop over all rr
574 [ + + ]: 41212419 : do inpt=1,npts
575 : :
576 : : ! Load module of rr
577 : 40744711 : rnorm=one
578 [ + + ]: 40744711 : if (normchoice==1) rnorm=nrm(inpt)
579 : :
580 : : ! Continue only for r<>0
581 : :
582 [ + + ]: 40773775 : if (rnorm>tol10) then
583 : :
584 : : ! Determine theta and phi
585 : 40707744 : cphi=one
586 : 40707744 : sphi=zero
587 : 40707744 : ctheta=rr(3,inpt)/rnorm
588 : : ! LB-2025/01:
589 : : ! If ctheta is too close to 1 (or -1), then stheta is small with poor accuracy,
590 : : ! which leads to a numerical instability visible (but negligible) in stress and forces components.
591 : : ! Here we have stheta>1e-6 with a minimum of 2 correct digits when close to 1e-6.
592 [ + + ]: 40707744 : if (abs(one-ctheta)<tol12) then
593 : 117289 : ctheta=one
594 : 117289 : stheta=zero
595 [ + + ]: 40590455 : else if (abs(one+ctheta)<tol12) then
596 : 117407 : ctheta=-one
597 : 117407 : stheta=zero
598 : : else
599 : 40473048 : stheta=sqrt(one-ctheta*ctheta)
600 : 40473048 : cphi=rr(1,inpt)/(rnorm*stheta)
601 : 40473048 : sphi=rr(2,inpt)/(rnorm*stheta)
602 : : end if
603 [ + + ]: 135644594 : do mm=1,mpsang-1
604 : 94936850 : rphase(mm)=dreal(dcmplx(cphi,sphi)**mm)
605 : 135644594 : iphase(mm)=aimag(dcmplx(cphi,sphi)**mm)
606 : : end do
607 : :
608 : : ! Determine gradients of theta and phi
609 [ + + ]: 40707744 : if (compute_ylmgr) then
610 : 10468806 : dtheta(1)=ctheta*cphi
611 : 10468806 : dtheta(2)=ctheta*sphi
612 : 10468806 : dtheta(3)=-stheta
613 : 10468806 : dphi(1)=-sphi
614 : 10468806 : dphi(2)=cphi
615 : 10468806 : dphi(3)=zero
616 : : end if
617 : :
618 : : ! COMPUTE Ylm(R)
619 [ + - ]: 40707744 : if (compute_ylm) then
620 : : ! Loop over angular momentum l
621 [ + + ]: 135644594 : do ilang=2,mpsang
622 : 94936850 : ll=ilang-1
623 : 94936850 : l0=ll**2+ll+1
624 : 94936850 : fact=1._dp/real(ll*(ll+1),dp)
625 : 94936850 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)
626 : : ! Special case m=0
627 : 94936850 : ylmr(l0,inpt)=ylmcst*ass_leg_pol(ll,0,ctheta)
628 : : ! Compute for m>0
629 : 94936850 : onem=one
630 [ + + ]: 309834480 : do mm=1,ll
631 : 174189886 : onem=-onem
632 : 174189886 : work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp)
633 : 174189886 : ylmr(l0+mm,inpt)=work1*rphase(mm)
634 : 174189886 : ylmr(l0-mm,inpt)=work1*iphase(mm)
635 [ + + ]: 269126736 : if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
636 : : end do ! End loop over m
637 : : end do ! End loop over l
638 : : end if
639 : :
640 : : ! COMPUTE dYlm/dRi
641 [ + + ]: 40707744 : if (compute_ylmgr) then
642 : : ! Loop over angular momentum l
643 [ + + ]: 38793150 : do ilang=2,mpsang
644 : 28324344 : ll=ilang-1
645 : 28324344 : l0=ll**2+ll+1
646 : 28324344 : fact=1._dp/real(ll*(ll+1),dp)
647 : 28324344 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rnorm
648 : : ! Special case m=0
649 : 28324344 : work1=ylmcst*plm_dtheta(ll,0,ctheta)
650 [ + + ]: 113297376 : ylmr_gr(1:3,l0,inpt)=work1*dtheta(1:3)
651 : : ! Compute for m>0
652 : 28324344 : onem=one
653 [ + + ]: 98424084 : do mm=1,ll
654 : 59630934 : onem=-onem
655 : 59630934 : work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp)
656 : 59630934 : work2=ylmcst*sqrt(fact)*onem*plm_dphi (ll,mm,ctheta)*sqrt(2._dp)
657 [ + + ]: 238523736 : ylmr_gr(1:3,l0+mm,inpt)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3)
658 [ + + ]: 238523736 : ylmr_gr(1:3,l0-mm,inpt)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3)
659 [ + + ]: 87955278 : if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
660 : : end do ! End loop over m
661 : : end do ! End loop over l
662 : : end if
663 : :
664 : : ! COMPUTE d2Ylm/dRidRj
665 [ + + ]: 40707744 : if (compute_ylm2gr) then
666 [ - + + - : 591351 : LIBPAW_ALLOCATE(blm,(5,mpsang*mpsang))
+ - - + -
+ ]
667 : 197117 : call plm_coeff(blm,mpsang,ctheta)
668 : :
669 : : ! Loop over angular momentum l
670 [ + + ]: 784651 : do ilang=2,mpsang
671 : 587534 : ll=ilang-1
672 : 587534 : l0=ll**2+ll+1
673 : 587534 : fact=1._dp/real(ll*(ll+1),dp)
674 : 587534 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rnorm**2)
675 : : ! Special case m=0
676 : 587534 : ylmr_gr(4,l0,inpt)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi)
677 : 587534 : ylmr_gr(5,l0,inpt)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi)
678 : 587534 : ylmr_gr(6,l0,inpt)=ylmcst*blm(1,l0)
679 : 587534 : ylmr_gr(7,l0,inpt)=ylmcst*blm(2,l0)*sphi
680 : 587534 : ylmr_gr(8,l0,inpt)=ylmcst*blm(2,l0)*cphi
681 : 587534 : ylmr_gr(9,l0,inpt)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi
682 : : ! Compute for m>0
683 : 587534 : onem=one
684 [ + + ]: 2052396 : do mm=1,ll
685 : 1267745 : onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two)
686 : : ylmr_gr(4,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-&
687 : 1267745 : & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
688 : : ylmr_gr(4,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+&
689 : 1267745 : & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
690 : : ylmr_gr(5,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+&
691 : 1267745 : & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
692 : : ylmr_gr(5,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-&
693 : 1267745 : & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
694 : 1267745 : ylmr_gr(6,l0+mm,inpt)=ylmcst2*blm(1,l0+mm)*rphase(mm)
695 : 1267745 : ylmr_gr(6,l0-mm,inpt)=ylmcst2*blm(1,l0+mm)*iphase(mm)
696 : : ylmr_gr(7,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+&
697 : 1267745 : & mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
698 : : ylmr_gr(7,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-&
699 : 1267745 : & mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
700 : : ylmr_gr(8,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-&
701 : 1267745 : & mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
702 : : ylmr_gr(8,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+&
703 : 1267745 : & mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
704 : : ylmr_gr(9,l0+mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-&
705 : 1267745 : & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm))
706 : : ylmr_gr(9,l0-mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+&
707 : 1267745 : & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm))
708 [ + + ]: 1855279 : if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
709 : : end do ! End loop over m
710 : : end do ! End loop over l
711 [ - + ]: 197117 : LIBPAW_DEALLOCATE(blm)
712 : : end if
713 : :
714 : : ! End condition r<>0
715 : : end if
716 : :
717 : : ! End loop over rr
718 : : end do
719 : :
720 : : ! End condition l<>0
721 : : end if
722 : :
723 [ + + ]: 470018 : end subroutine initylmr
724 : : !!***
725 : :
726 : : !----------------------------------------------------------------------
727 : :
728 : : !!****f* m_paw_sphharm/ys
729 : : !! NAME
730 : : !! ys
731 : : !!
732 : : !! FUNCTION
733 : : !! Computes the matrix element <Y_(l2,m2)|S_(l1,m1)>
734 : : !!
735 : : !! INPUTS
736 : : !! integer :: l2,m2,l1,m1
737 : : !!
738 : : !! OUTPUT
739 : : !! complex(dp) :: ys_val
740 : : !!
741 : : !! NOTES
742 : : !! Ylm is the standard complex-valued spherical harmonic, Slm is the real spherical harmonic
743 : : !! used througout abinit.
744 : : !!
745 : : !! SOURCE
746 : :
747 : 533124 : subroutine ys(l2,m2,l1,m1,ys_val)
748 : :
749 : : !Arguments ---------------------------------------------
750 : : !scalars
751 : : integer,intent(in) :: l1,l2,m1,m2
752 : : complex(dp),intent(out) :: ys_val
753 : :
754 : : !Local variables ---------------------------------------
755 : : !scalars
756 : : integer :: mp1
757 : :
758 : : ! *********************************************************************
759 : :
760 : : ! See Blanco et al., J. Mol Struct. 419, 19-27 (1997) Eq. 19
761 : : ! <Y_l2,m2|S_l1,m1> is given by C^l_{m1,m2} where
762 : : ! l1 == l2 and |m1| == |m2|, 0 otherwise
763 : :
764 : 533124 : ys_val = czero
765 : :
766 [ + - ]: 533124 : if ( l2 /= l1 ) return
767 [ + + ]: 533124 : if ( abs(m2) /= abs(m1) ) return
768 : :
769 : 393780 : mp1=(-1)**abs(m1)
770 : :
771 [ + + ]: 393780 : if(m1.EQ.0) then
772 : 25116 : ys_val=cone
773 [ + + + + ]: 368664 : else if((m1.GT.0).AND.(m2.GT.0)) then
774 : 95106 : ys_val=mp1*sqrthalf
775 [ + + + - ]: 273558 : else if((m1.GT.0).AND.(m2.LT.0)) then
776 : 95106 : ys_val=sqrthalf
777 [ + - + + ]: 178452 : else if((m1.LT.0).AND.(m2.GT.0)) then
778 : 89226 : ys_val=-j_dpc*mp1*sqrthalf
779 [ + - + - ]: 89226 : else if((m1.LT.0).AND.(m2.LT.0)) then
780 : 89226 : ys_val=j_dpc*sqrthalf
781 : : else
782 : : ys_val=czero
783 : : end if
784 : :
785 : : end subroutine ys
786 : : !!***
787 : :
788 : : !----------------------------------------------------------------------
789 : :
790 : : !!****f* m_paw_sphharm/lxyz.F90
791 : : !! NAME
792 : : !! lxyz
793 : : !!
794 : : !! FUNCTION
795 : : !! Computes the matrix element <Yl'm'|L_idir|Ylm>
796 : : !!
797 : : !! INPUTS
798 : : !! integer :: lp,mp,idir,ll,mm
799 : : !!
800 : : !! OUTPUT
801 : : !! complex(dp) :: lidir
802 : : !!
803 : : !! NOTES
804 : : !! Ylm is the standard complex-valued spherical harmonic,
805 : : !! idir is the direction in space of L
806 : : !!
807 : : !! SOURCE
808 : :
809 : 215163 : subroutine lxyz(lp,mp,idir,ll,mm,lidir)
810 : :
811 : : !Arguments ---------------------------------------------
812 : : !scalars
813 : : integer,intent(in) :: idir,ll,lp,mm,mp
814 : : complex(dp),intent(out) :: lidir
815 : :
816 : : !Local variables ---------------------------------------
817 : : !scalars
818 : : complex(dp) :: jme, jmme, jpme
819 : : ! *********************************************************************
820 : :
821 : 215163 : lidir = czero
822 [ + - ]: 215163 : if ( lp /= ll ) return
823 : :
824 : 215163 : jpme=czero; jmme=czero; jme=czero
825 [ + + ]: 215163 : if (mp==mm) then
826 : 62721 : jme=cone*mm
827 [ + + ]: 152442 : else if (mp==mm+1) then
828 : 43614 : jpme=-cone*sqrt(half*((ll*(ll+1))-mm*(mm+1)))
829 [ + + ]: 108828 : else if (mp==mm-1) then
830 : 43614 : jmme= cone*sqrt(half*((ll*(ll+1))-mm*(mm-1)))
831 : : end if
832 : :
833 : 319959 : select case (idir)
834 : : case (1) ! Lx
835 : 104796 : lidir = -sqrthalf*(jpme - jmme)
836 : : case (2) ! Ly
837 : 46584 : lidir = j_dpc*sqrthalf*(jpme + jmme)
838 : : case (3) ! Lz
839 [ + + + - ]: 215163 : lidir = jme
840 : : end select
841 : :
842 : : end subroutine lxyz
843 : : !!***
844 : :
845 : : !----------------------------------------------------------------------
846 : :
847 : : !!****f* m_paw_sphharm/slxyzs
848 : : !! NAME
849 : : !! slxyzs
850 : : !!
851 : : !! FUNCTION
852 : : !! computes the matrix element <Sl'm'|L_idir|Slm>
853 : : !!
854 : : !! INPUTS
855 : : !! integer :: lp,mp,idir,ll,mm
856 : : !!
857 : : !! OUTPUT
858 : : !! complex(dp) :: sls_val
859 : : !!
860 : : !! NOTES
861 : : !! Slm is the real spherical harmonic used througout abinit,
862 : : !! L_idir is a component of the angular momentum operator.
863 : : !! The subroutine computes <S_l'm'|L_idir|S_lm>
864 : : !!
865 : : !! SOURCE
866 : :
867 : 19107 : subroutine slxyzs(lp,mp,idir,ll,mm,sls_val)
868 : :
869 : : !Arguments ---------------------------------------------
870 : : !scalars
871 : : integer,intent(in) :: idir,ll,lp,mm,mp
872 : : complex(dp),intent(out) :: sls_val
873 : :
874 : : !Local variables ---------------------------------------
875 : : !scalars
876 : : integer :: mpp,mppp
877 : : complex(dp) :: lidir,sy_val,ys_val
878 : : ! *********************************************************************
879 : :
880 : 19107 : sls_val = czero
881 : :
882 [ - + ]: 19107 : if ( lp /= ll ) return
883 : :
884 [ + + ]: 81828 : do mpp = -ll, ll
885 : 62721 : call ys(ll,mpp,ll,mp,sy_val)
886 [ + + ]: 296991 : do mppp = -ll, ll
887 : 215163 : call lxyz(ll,mpp,idir,ll,mppp,lidir)
888 : 215163 : call ys(ll,mppp,ll,mm,ys_val)
889 : 277884 : sls_val = sls_val + conjg(sy_val)*lidir*ys_val
890 : : end do
891 : : end do
892 : :
893 : : end subroutine slxyzs
894 : : !!***
895 : :
896 : : !----------------------------------------------------------------------
897 : :
898 : : !!****f* m_paw_sphharm/lsylm
899 : : !! NAME
900 : : !! lsylm
901 : : !!
902 : : !! FUNCTION
903 : : !! Compute the LS operator in the real spherical harmonics basis
904 : : !! ls_ylm(ilm1,ilm2,ispin)= <sigma, S_lm1| L.S |S_lm2, sigma_prime>
905 : : !! ilm,1m2=(l,m1,m2) with -l<=m1<=l, -l<=m2<=l and 0<l<=lmax
906 : : !! ispin=(sigma,sigma_prime) 1=(up,up), 2=(up,dn), 3=(dn,up), 4=(dn,dn)
907 : : !!
908 : : !! INPUTS
909 : : !! lmax= max. value of angular momentum l
910 : : !!
911 : : !! OUTPUT
912 : : !! ls_ylm(2,l_max**2*(l_max**2+1)/2,2)=LS operator in the real spherical harmonics basis
913 : : !! ls_ylm(:,:,1)=<up, S_lm1| L.S |S_lm2, up>
914 : : !! ls_ylm(:,:,2)=<up, S_lm1| L.S |S_lm2, down>
915 : : !! One can deduce:
916 : : !! <down, S_lm1| L.S |S_lm2, down>=-<up, S_lm1| L.S |S_lm2, up>
917 : : !! <down, S_lm1| L.S |S_lm2, up> =-Conjg[<up, S_lm1| L.S |S_lm2, down>]
918 : : !! Also, only ilm1<=ilm2 terms are stored, because:
919 : : !! <sigma, S_lm1| L.S |S_lm2, sigma_prime>=-<sigma_prime, S_lm1| L.S |S_lm2, sigma>
920 : : !!
921 : : !! SOURCE
922 : :
923 : 67 : subroutine lsylm(ls_ylm,lmax)
924 : :
925 : : !Arguments ---------------------------------------------
926 : : !scalars
927 : : integer,intent(in) :: lmax
928 : : !arrays
929 : : real(dp),allocatable :: ls_ylm(:,:,:)
930 : :
931 : : !Local variables ---------------------------------------
932 : : !scalars
933 : : integer :: ii,ilm,im,j0lm,jj,jlm,jm,klm,ll,lm0,mm,ispden
934 : : real(dp),parameter :: invsqrt2=one/sqrt2
935 : : real(dp) :: onem
936 : : character(len=500) :: msg
937 : : logical,parameter :: tso=.false. ! use true to Test Spin Orbit and
938 : : ! write the matrix of L.S in different basis
939 : : !arrays
940 : : complex(dp) :: tmp(2)
941 : 67 : complex(dp),allocatable :: ls_cplx(:,:,:),slm2ylm(:,:)
942 : 67 : complex(dp),allocatable :: mat_inp_c(:,:,:),mat_out_c(:,:,:)
943 : 67 : complex(dp),allocatable :: mat_ls_ylm(:,:,:),mat_jmj(:,:)
944 : : character(len=9),parameter :: dspin2(2)=(/"up-up ","up-dn "/)
945 : : character(len=9),parameter :: dspin6(6)=(/"dn ","up ","dn-dn ","up-up ","dn-up ","up-dn "/)
946 : : character(len=9),parameter :: dspinm(6)=(/"dn ","up ","n ","mx ","my ","mz "/)
947 : : ! *************************************************************************
948 : :
949 [ - + ]: 67 : if (.not.allocated(ls_ylm)) then
950 : 0 : msg='ls_ylm is not allocated!'
951 : 0 : LIBPAW_BUG(msg)
952 : : end if
953 [ + + - + ]: 268 : if ( size(ls_ylm) < 2*(lmax+1)**2 * ((lmax+1)**2+1) ) then
954 : 0 : msg='wrong size for ls_ylm!'
955 : 0 : LIBPAW_BUG(msg)
956 : : end if
957 : :
958 : : !Initialization
959 [ + + + + : 14217 : ls_ylm=zero
+ + ]
960 : :
961 : : !Nothing to do if lmax=0
962 [ - + ]: 67 : if (lmax<=0) return
963 : :
964 : : !Loop on l quantum number
965 [ + + ]: 174 : do ll=1,lmax
966 : :
967 : : ! Transformation matrixes: real->complex spherical harmonics
968 [ - + + - : 428 : LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1))
- + + - +
- - + -
+ ]
969 [ + + + + ]: 2223 : slm2ylm=czero
970 [ + + ]: 516 : do im=1,2*ll+1
971 : 409 : mm=im-ll-1;jm=-mm+ll+1
972 : 409 : onem=dble((-1)**mm)
973 [ + + ]: 409 : if (mm> 0) then
974 : 151 : slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
975 : 151 : slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
976 : : end if
977 [ + + ]: 409 : if (mm==0) then
978 : 107 : slm2ylm(im,im)=cone
979 : : end if
980 [ + + ]: 516 : if (mm< 0) then
981 : 151 : slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
982 : 151 : slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
983 : : end if
984 : : end do
985 : :
986 : : ! Compute <sigma, Y_lm1|L.S|Y_lm2, sigma_prime> (Y_lm=complex spherical harmonics)
987 : : ! 1= <up|L.S|up> ; 2= <up|L.S|dn>
988 [ - + + - : 535 : LIBPAW_ALLOCATE(ls_cplx,(2*ll+1,2*ll+1,2))
- + + - +
- + - - +
- + ]
989 [ + + + + : 4553 : ls_cplx=czero
+ + ]
990 : : if(tso) then
991 : : LIBPAW_ALLOCATE(mat_ls_ylm,(2*ll+1,2*ll+1,4))
992 : : if(tso) mat_ls_ylm=czero
993 : : end if
994 : : if(tso) then
995 : : LIBPAW_ALLOCATE(mat_jmj,(2*(2*ll+1),2*(2*ll+1)))
996 : : if(tso) mat_jmj=czero
997 : : end if
998 [ + + ]: 516 : do im=1,2*ll+1
999 : 409 : mm=im-ll-1
1000 : 409 : ls_cplx(im,im,1)=half*mm
1001 : : if(tso) mat_ls_ylm(im,im,1)=-half*mm ! dn dn
1002 : : if(tso) mat_ls_ylm(im,im,2)=half*mm ! up up
1003 [ + + ]: 409 : if ((mm+1)<= ll) then
1004 : 302 : ls_cplx(im,im+1,2)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp))
1005 : : if(tso) mat_ls_ylm(im,im+1,4)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp)) ! up dn
1006 : : if(tso) mat_ls_ylm(im+1,im,3)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp)) ! dn up
1007 : : end if
1008 [ + + ]: 516 : if ((mm-1)>=-ll) then
1009 : 302 : ls_cplx(im-1,im,2)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp))
1010 : : if(tso) mat_ls_ylm(im-1,im,4)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp)) ! up dn
1011 : : if(tso) mat_ls_ylm(im,im-1,3)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp)) ! dn up
1012 : : end if
1013 : : end do
1014 : :
1015 : : ! test : print LS in J,M_J basis
1016 : : if(tso) then
1017 : : do ispden=1,4
1018 : : write(msg,'(3a)') ch10,"value of LS in the Ylm basis for " ,trim(dspin6(ispden+2*(4/4)))
1019 : : call wrtout(std_out,msg,'COLL')
1020 : : do im=1,ll*2+1
1021 : : write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1)
1022 : : call wrtout(std_out,msg,'COLL')
1023 : : end do
1024 : : end do
1025 : : call mat_mlms2jmj(ll,mat_ls_ylm,mat_jmj,4,1,2,3,std_out,'COLL') ! optspin=2 : dn spin are first
1026 : : end if
1027 : :
1028 : : ! Compute <sigma, S_lm1|L.S|S_lm2, sigma_prime> (S_lm=real spherical harmonics)
1029 : : ! 1= <up|L.S|up> ; 2= <up|L.S|dn>
1030 : : if(tso) then
1031 : : LIBPAW_ALLOCATE(mat_inp_c,(2*ll+1,2*ll+1,4))
1032 : : LIBPAW_ALLOCATE(mat_out_c,(2*ll+1,2*ll+1,4))
1033 : : end if
1034 : 107 : lm0=ll**2
1035 [ + + ]: 516 : do jm=1,2*ll+1
1036 : 409 : jlm=lm0+jm;j0lm=jlm*(jlm-1)/2
1037 [ + + ]: 1574 : do im=1,jm
1038 : 1058 : ilm=lm0+im;klm=j0lm+ilm
1039 : 1058 : tmp(:)=czero
1040 [ + + ]: 5836 : do ii=1,2*ll+1
1041 [ + + ]: 29718 : do jj=1,2*ll+1
1042 [ + + ]: 76424 : tmp(:)=tmp(:)+ls_cplx(ii,jj,:)*CONJG(slm2ylm(ii,im))*slm2ylm(jj,jm)
1043 : : end do
1044 : : end do
1045 [ + + ]: 3174 : ls_ylm(1,klm,:)=REAL(tmp(:),kind=dp)
1046 [ + + ]: 3583 : ls_ylm(2,klm,:)=AIMAG(tmp(:))
1047 : : end do
1048 : : end do
1049 : :
1050 : : ! Test: print LS in Slm basis
1051 : : if(tso) then
1052 : : call mat_slm2ylm(ll,mat_ls_ylm,mat_inp_c,4,2,2,3,std_out,'COLL') ! from Ylm to Slm, and dn spin are first
1053 : : do ispden=1,4
1054 : : write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspin6(ispden+2*(4/4)))
1055 : : call wrtout(std_out,msg,'COLL')
1056 : : do im=1,ll*2+1
1057 : : write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_inp_c(im,jm,ispden),jm=1,ll*2+1)
1058 : : call wrtout(std_out,msg,'COLL')
1059 : : end do
1060 : : end do
1061 : : ! change into n,m basis
1062 : : mat_ls_ylm(:,:,1)=(mat_inp_c(:,:,1)+mat_inp_c(:,:,2))
1063 : : mat_ls_ylm(:,:,2)=(mat_inp_c(:,:,3)+mat_inp_c(:,:,4))
1064 : : mat_ls_ylm(:,:,3)=-cmplx(0.d0,1.d0)*(mat_inp_c(:,:,4)-mat_inp_c(:,:,3))
1065 : : mat_ls_ylm(:,:,4)=(mat_inp_c(:,:,1)-mat_inp_c(:,:,2))
1066 : : do ispden=1,4
1067 : : write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspinm(ispden+2*(4/4)))
1068 : : call wrtout(std_out,msg,'COLL')
1069 : : do im=1,ll*2+1
1070 : : write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1)
1071 : : call wrtout(std_out,msg,'COLL')
1072 : : end do
1073 : : end do
1074 : : LIBPAW_DEALLOCATE(mat_inp_c)
1075 : : LIBPAW_DEALLOCATE(mat_ls_ylm)
1076 : : LIBPAW_DEALLOCATE(mat_jmj)
1077 : : LIBPAW_DEALLOCATE(mat_out_c)
1078 : : end if ! tso
1079 : :
1080 [ - + ]: 107 : LIBPAW_DEALLOCATE(ls_cplx)
1081 [ - + ]: 174 : LIBPAW_DEALLOCATE(slm2ylm)
1082 : :
1083 : : ! End loop on l
1084 : : end do
1085 : :
1086 : : end subroutine lsylm
1087 : : !!***
1088 : :
1089 : : !----------------------------------------------------------------------
1090 : :
1091 : : !!****f* m_paw_sphharm/plm_coeff
1092 : : !! NAME
1093 : : !! plm_coeff
1094 : : !!
1095 : : !! FUNCTION
1096 : : !! Compute coefficients depending on Plm and its derivatives where P_lm is a legendre polynomial.
1097 : : !! They are used to compute the second derivatives of spherical harmonics
1098 : : !!
1099 : : !! INPUTS
1100 : : !! mpsang=1+ maximum l quantum number
1101 : : !! xx= input value
1102 : : !!
1103 : : !! OUTPUT
1104 : : !! blm(5,mpsang*mpsang)=coefficients depending on Plm and its derivatives where P_lm is a legendre polynome
1105 : : !!
1106 : : !! SOURCE
1107 : :
1108 : 4222762 : subroutine plm_coeff(blm,mpsang,xx)
1109 : :
1110 : : !Arguments ---------------------------------------------
1111 : : !scalars
1112 : : integer,intent(in) :: mpsang
1113 : : real(dp),intent(in) :: xx
1114 : : !arrays
1115 : : real(dp),intent(out) :: blm(5,mpsang*mpsang)
1116 : :
1117 : : !Local variables ---------------------------------------
1118 : : !scalars
1119 : : integer :: il,ilm,ilm0,ilm1,im
1120 : : real(dp) :: dplm_dt,d2plm_dt2,llp1,onemx2,plm,sqrx,xsqrx,xx2,yy
1121 : : logical :: is_one
1122 : : character(len=500) :: msg
1123 : : !arrays
1124 : 8445524 : real(dp) :: pl_d2(mpsang),plm_d2t(mpsang*mpsang)
1125 : :
1126 : : !************************************************************************
1127 : :
1128 [ - + ]: 4222762 : if (abs(xx).gt.1.d0) then
1129 : 0 : msg = ' plm_coeff : xx > 1 !'
1130 : 0 : LIBPAW_ERROR(msg)
1131 : : end if
1132 : :
1133 [ + + + + ]: 238089142 : blm=zero
1134 : 4222762 : is_one=(abs(abs(xx)-one)<=tol12)
1135 : 4222762 : xx2=xx**2
1136 : 4222762 : onemx2=abs(one-xx2)
1137 : 4222762 : sqrx=sqrt(onemx2)
1138 : 4222762 : xsqrx=xx*sqrt(onemx2)
1139 : :
1140 : 4222762 : call plm_d2theta(mpsang,plm_d2t,xx)
1141 [ + + ]: 4222762 : if (is_one) then
1142 : 1583 : yy=sign(one,xx)
1143 : 1583 : call pl_deriv(mpsang,pl_d2,yy)
1144 : : end if
1145 : :
1146 [ + + ]: 16837556 : do il=0,mpsang-1
1147 : 12614794 : llp1=dble(il*(il+1))
1148 : 12614794 : ilm0=il*il+il+1
1149 [ + + ]: 42633818 : do im=0,il
1150 : 25796262 : ilm=ilm0+im;ilm1=ilm0-im
1151 : :
1152 : 25796262 : plm =(-1)**im*ass_leg_pol(il,im,xx)
1153 : 25796262 : dplm_dt =(-1)**im*plm_dtheta(il,im,xx)
1154 : 25796262 : d2plm_dt2= plm_d2t(ilm)
1155 : :
1156 : 25796262 : blm(1,ilm)= two*xsqrx *dplm_dt+onemx2*d2plm_dt2
1157 : 25796262 : blm(2,ilm)= (one-two*xx2)*dplm_dt-xsqrx *d2plm_dt2
1158 : 25796262 : blm(3,ilm)=llp1*plm+ d2plm_dt2
1159 : 25796262 : blm(4,ilm)= -two*xsqrx *dplm_dt+xx2 *d2plm_dt2
1160 : :
1161 : :
1162 [ + + ]: 25796262 : if (is_one) then
1163 [ + + ]: 9204 : if (im==1) then
1164 : 2934 : blm(5,ilm)=llp1*plm+d2plm_dt2
1165 : : end if
1166 [ + + ]: 9204 : if (im==2) then
1167 : 1351 : blm(5,ilm)=d2plm_dt2-three*pl_d2(il+1)
1168 : : end if
1169 : : else
1170 [ + + ]: 25787058 : if(im>0) then
1171 : 13176781 : blm(5,ilm)=plm/onemx2-dplm_dt*xx/sqrx
1172 : : end if
1173 : : end if
1174 : :
1175 [ + + ]: 38411056 : if (im>0) then
1176 : 13181468 : blm(1,ilm1)=blm(1,ilm)
1177 : 13181468 : blm(2,ilm1)=blm(2,ilm)
1178 : 13181468 : blm(3,ilm1)=blm(3,ilm)
1179 : 13181468 : blm(4,ilm1)=blm(4,ilm)
1180 : 13181468 : blm(5,ilm1)=blm(5,ilm)
1181 : : end if
1182 : :
1183 : : end do
1184 : : end do
1185 : :
1186 : 4222762 : end subroutine plm_coeff
1187 : : !!***
1188 : :
1189 : : !----------------------------------------------------------------------
1190 : :
1191 : : !!****f* m_paw_sphharm/ass_leg_pol
1192 : : !! NAME
1193 : : !! ass_leg_pol
1194 : : !!
1195 : : !! FUNCTION
1196 : : !! Compute the associated Legendre Polynomial Plm(x),
1197 : : !! using a stable recursion formula.
1198 : : !! Here m and l are integers satisfying 0<=m<=l,
1199 : : !! while x lies in the range -1<=x<=1
1200 : : !!
1201 : : !! INPUTS
1202 : : !! l,m= l,m numbers
1203 : : !! xarg=argument of the polynom
1204 : : !!
1205 : : !! OUTPUT
1206 : : !!
1207 : : !! SOURCE
1208 : :
1209 : 394520503 : function ass_leg_pol(l,m,xarg)
1210 : :
1211 : : !Arguments ------------------------------------
1212 : : !scalars
1213 : : integer, intent(in) :: l,m
1214 : : real(dp), intent(in) :: xarg
1215 : : real(dp) :: ass_leg_pol
1216 : :
1217 : : !Local variables-------------------------------
1218 : : !scalars
1219 : : integer :: i,ll
1220 : : real(dp) :: pll,polmm,tmp1,sqrx,x
1221 : : character(len=100) :: msg
1222 : :
1223 : : ! *************************************************************************
1224 : :
1225 : 394520503 : x=xarg
1226 [ + - + - : 394520503 : if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0) then
- + ]
1227 [ # # # # : 0 : if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0+1.d-10) then
# # ]
1228 : 0 : msg='Bad choice of l, m or x !'
1229 : 0 : LIBPAW_BUG(msg)
1230 : : endif
1231 : : x=1.d0
1232 : : endif
1233 : :
1234 : 394520503 : polmm=1.d0
1235 [ + + ]: 394520503 : if (m>0) then
1236 : 243233327 : sqrx=sqrt(abs((1.d0-x)*(1.d0+x)))
1237 [ + + ]: 632652314 : do i=1,m
1238 : 632652314 : polmm=polmm*(1.0d0-2.0d0*i)*sqrx
1239 : : enddo
1240 : : endif
1241 : :
1242 [ + + ]: 394520503 : if (l==m) then
1243 : : ass_leg_pol=polmm
1244 : : else
1245 : 239064057 : tmp1=x*(2.0d0*m+1.0d0)*polmm
1246 [ + + ]: 239064057 : if (l==(m+1)) then
1247 : : ass_leg_pol=tmp1
1248 : : else
1249 [ + + ]: 246010661 : do ll=m+2,l
1250 : 145565494 : pll=(x*(2.0d0*ll-1.0d0)*tmp1-(ll+m-1.0d0)*polmm)/dble(ll-m)
1251 : 145565494 : polmm=tmp1
1252 : 246010661 : tmp1=pll
1253 : : enddo
1254 : : ass_leg_pol=pll
1255 : : endif
1256 : : endif
1257 : :
1258 : 394520503 : end function ass_leg_pol
1259 : : !!***
1260 : :
1261 : : !----------------------------------------------------------------------
1262 : :
1263 : : !!****f* m_paw_sphharm/plm_dphi
1264 : : !! NAME
1265 : : !! plm_dphi
1266 : : !!
1267 : : !! FUNCTION
1268 : : !! Compute m*P_lm(x)/sqrt((1-x^2)where P_lm is a legendre polynome
1269 : : !!
1270 : : !! INPUTS
1271 : : !! ll= l quantum number
1272 : : !! mm= m quantum number
1273 : : !! xx= input value
1274 : : !!
1275 : : !! OUTPUT
1276 : : !! plm_dphi(xx)
1277 : : !!
1278 : : !! NOTES
1279 : : !! This routine comes from Function Der_Phi_P(L,m,x)
1280 : : !! (pwpaw code from N. Holzwarth, implemented by Y. Abraham))
1281 : : !!
1282 : : !! SOURCE
1283 : :
1284 : 81783805 : function plm_dphi(ll,mm,xx)
1285 : :
1286 : : !Arguments ---------------------------------------------
1287 : : !scalars
1288 : : integer,intent(in) :: ll,mm
1289 : : real(dp) :: plm_dphi
1290 : : real(dp),intent(in) :: xx
1291 : :
1292 : : !Local variables ---------------------------------------
1293 : : !scalars
1294 : : integer :: il,im
1295 : : real(dp) :: dosomx2,fact,pll,pmm,pmmp1,somx2
1296 : : character(len=500) :: msg
1297 : :
1298 : : ! *********************************************************************
1299 : :
1300 [ + - + - : 81783805 : if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then
- + ]
1301 : 0 : msg = 'plm_dphi : mm < 0 or mm > ll or xx > 1 !'
1302 : 0 : LIBPAW_ERROR(msg)
1303 : : end if
1304 : :
1305 : 81783805 : plm_dphi=zero
1306 [ + - ]: 81783805 : if (mm==0) return
1307 : :
1308 : 81783805 : pmm=one
1309 : 81783805 : dosomx2=one
1310 [ + - ]: 81783805 : if (mm > 0) then
1311 : 81783805 : somx2=sqrt((1-xx)*(1+xx))
1312 : 81783805 : fact=one
1313 [ + + ]: 223829778 : do im=1,mm
1314 : 142045973 : pmm=-pmm*fact
1315 : 223829778 : fact=fact+2
1316 : : end do
1317 [ + + ]: 81783805 : if (mm > 1) then
1318 [ + + ]: 98089053 : do im=2,mm
1319 : 98089053 : dosomx2=somx2*dosomx2
1320 : : end do
1321 : : end if
1322 : 81783805 : pmm=pmm*dosomx2 !due to one more term (-1^M)
1323 : : end if
1324 [ + + ]: 81783805 : if(ll==mm) then
1325 : 43956920 : plm_dphi=pmm*mm
1326 : : else
1327 : 37826885 : pmmp1=xx*(2*mm+1)*pmm
1328 [ + + ]: 37826885 : if(ll==mm+1) then
1329 : 23911136 : plm_dphi=pmmp1*mm
1330 [ + - ]: 13915749 : else if(ll>=mm+2) then
1331 [ + + ]: 36351032 : do il=mm+2,ll
1332 : 22435283 : pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1333 : 22435283 : pmm=pmmp1
1334 : 36351032 : pmmp1=pll
1335 : : end do
1336 : 13915749 : plm_dphi=pll*mm
1337 : : end if
1338 : : end if
1339 : :
1340 : : end function plm_dphi
1341 : : !!***
1342 : :
1343 : : !----------------------------------------------------------------------
1344 : :
1345 : : !!****f* m_paw_sphharm/plm_dtheta
1346 : : !! NAME
1347 : : !! plm_dtheta
1348 : : !!
1349 : : !! FUNCTION
1350 : : !! Compute -(1-x^2)^1/2*d/dx{P_lm(x)} where P_lm is a legendre polynome
1351 : : !!
1352 : : !! INPUTS
1353 : : !! ll= l quantum number
1354 : : !! mm= m quantum number
1355 : : !! xx= input value
1356 : : !!
1357 : : !! OUTPUT
1358 : : !! plm_dtheta(xx)
1359 : : !!
1360 : : !! NOTES
1361 : : !! This routine comes from Function Der_Theta_P(L,m,x)
1362 : : !! (pwpaw code from N. Holzwarth, implemented by Y. Abraham))
1363 : : !!
1364 : : !! SOURCE
1365 : :
1366 : 225836359 : function plm_dtheta(ll,mm,xx)
1367 : :
1368 : : !Arguments ---------------------------------------------
1369 : : !scalars
1370 : : integer,intent(in) :: ll,mm
1371 : : real(dp) :: plm_dtheta
1372 : : real(dp),intent(in) :: xx
1373 : :
1374 : : !Local variables ---------------------------------------
1375 : : !scalars
1376 : : integer :: il,im
1377 : : real(dp) :: dosomx2,dpll,dpmm,dpmmp1,fact,pll,pmm,pmmp1,somx2
1378 : : character(len=500) :: msg
1379 : :
1380 : : ! *********************************************************************
1381 : :
1382 [ + - + - : 225836359 : if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then
- + ]
1383 : 0 : msg = 'plm_dtheta : mm < 0 or mm > ll or xx > 1 !'
1384 : 0 : LIBPAW_ERROR(msg)
1385 : : end if
1386 : :
1387 : 225836359 : plm_dtheta=zero
1388 : 225836359 : pmm=one
1389 : 225836359 : dpmm=one
1390 : 225836359 : dosomx2=one
1391 : 225836359 : somx2=sqrt((1-xx)*(1+xx))
1392 [ + + ]: 225836359 : if(mm==0)then
1393 : : dpmm=zero
1394 [ + - ]: 156649851 : elseif (mm > 0) then
1395 : : fact=one
1396 [ + + ]: 402429176 : do im=1,mm
1397 : 245779325 : pmm=-pmm*fact*somx2
1398 : 245779325 : dpmm=-dpmm*fact
1399 : 402429176 : fact=fact+2
1400 : : end do
1401 [ + + ]: 156649851 : if(mm>1)then
1402 [ + + ]: 152047377 : do im=2,mm
1403 : 152047377 : dosomx2=dosomx2*somx2
1404 : : end do
1405 : : end if
1406 : 156649851 : dpmm= dpmm*mm*xx*dosomx2
1407 : : end if
1408 [ + + ]: 225836359 : if(ll==mm)then
1409 : : plm_dtheta=dpmm
1410 : : else
1411 : 118912453 : pmmp1=xx*(2*mm+1)*pmm
1412 : 118912453 : dpmmp1=-(2*mm+1)*somx2*pmm+xx*(2*mm+1)*dpmm
1413 [ + + ]: 118912453 : if(ll==mm+1) then
1414 : : plm_dtheta=dpmmp1
1415 [ + - ]: 45717151 : else if(ll>=mm+2)then
1416 [ + + ]: 115068715 : do il=mm+2,ll
1417 : 69351564 : pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1418 : 69351564 : dpll=(-somx2*(2*il-1)*pmmp1+(xx*(2*il-1)*dpmmp1-(il+mm-1)*dpmm))/(il-mm)
1419 : 69351564 : pmm=pmmp1
1420 : 69351564 : pmmp1=pll
1421 : 69351564 : dpmm=dpmmp1
1422 : 115068715 : dpmmp1=dpll
1423 : : end do
1424 : : plm_dtheta=dpll
1425 : : end if
1426 : : end if
1427 : :
1428 : 225836359 : end function plm_dtheta
1429 : : !!***
1430 : :
1431 : : !----------------------------------------------------------------------
1432 : :
1433 : : !!****f* m_paw_sphharm/plm_d2theta
1434 : : !! NAME
1435 : : !! plm_d2theta
1436 : : !!
1437 : : !! FUNCTION
1438 : : !! Compute d2(Plm (cos(theta)))/d(theta)2 where P_lm is a legendre polynome
1439 : : !!
1440 : : !! INPUTS
1441 : : !! mpsang=1+ maximum l quantum number
1442 : : !! xx= input value
1443 : : !!
1444 : : !! OUTPUT
1445 : : !! plm_d2t(mpsang*mpsang)
1446 : : !!
1447 : : !! SOURCE
1448 : :
1449 : 4222762 : subroutine plm_d2theta(mpsang,plm_d2t,xx)
1450 : :
1451 : : !Arguments ---------------------------------------------
1452 : : !scalars
1453 : : integer,intent(in) :: mpsang
1454 : : real(dp),intent(in) :: xx
1455 : : !arrays
1456 : : real(dp),intent(out) :: plm_d2t(mpsang*mpsang)
1457 : :
1458 : : !Local variables ---------------------------------------
1459 : : !scalars
1460 : : integer :: il,ilm,ilmm1,ilmm2,im
1461 : : real(dp) :: sqrx
1462 : : character(len=500) :: msg
1463 : :
1464 : : !************************************************************************
1465 [ - + ]: 4222762 : if (abs(xx).gt.1.d0) then
1466 : 0 : msg = 'plm_d2theta : xx > 1 !'
1467 : 0 : LIBPAW_ERROR(msg)
1468 : : end if
1469 : :
1470 [ + + ]: 43200492 : plm_d2t=zero
1471 [ + - ]: 4222762 : if (mpsang>1) then
1472 : 4222762 : sqrx=sqrt(abs((1.d0-xx)*(1.d0+xx)))
1473 : :
1474 [ + + ]: 12614794 : do il=1,mpsang-1
1475 : 8392032 : ilm=il*il+2*il+1
1476 : 8392032 : ilmm1=(il-1)*(il-1)+2*(il-1)+1
1477 : : ! terme d2(Pll)/dtet2
1478 : : plm_d2t(ilm)=(2*il-1)*(sqrx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))+&
1479 : 8392032 : & 2.d0*xx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1480 : 8392032 : plm_d2t(ilm-2*il)=plm_d2t(ilm)
1481 : : ! terme d2(Pl(l-1))/dtet2
1482 : : plm_d2t(ilm-1)=(2*il-1)*(xx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))-&
1483 : 8392032 : & 2.d0*sqrx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1484 [ + + ]: 12614794 : if(il>1) plm_d2t(il*il+2)=plm_d2t(ilm-1)
1485 : : end do
1486 : : ! terme d2(Plm)/dtet2
1487 [ + + ]: 4222762 : if(mpsang>2) then
1488 [ + + ]: 7814868 : do il=2,mpsang-1
1489 [ + + ]: 12604304 : do im=0,il-2
1490 : 4789436 : ilm=il*il+il+1+im
1491 : 4789436 : ilmm1=(il-1)*(il-1)+il+im
1492 : 4789436 : ilmm2=(il-2)*(il-2)+il-1+im
1493 : : plm_d2t(ilm)=dble(2*il-1)/dble(il-im)*(xx*(plm_d2t(ilmm1)-(-1)**im*ass_leg_pol(il-1,im,xx))-&
1494 : : & 2.d0*sqrx*(-1)**im*plm_dtheta(il-1,im,xx))-&
1495 : 4789436 : & dble(il+im-1)/dble(il-im)*plm_d2t(ilmm2)
1496 : 8958706 : plm_d2t(ilm-2*im)=plm_d2t(ilm)
1497 : : end do
1498 : : end do
1499 : : end if
1500 : : end if
1501 : :
1502 : 4222762 : end subroutine plm_d2theta
1503 : : !!***
1504 : :
1505 : : !----------------------------------------------------------------------
1506 : :
1507 : : !!****f* m_paw_sphharm/pl_deriv
1508 : : !! NAME
1509 : : !! pl_deriv
1510 : : !!
1511 : : !! FUNCTION
1512 : : !! Compute d2(Pl (x)))/d(x)2 where P_l is a legendre polynomial
1513 : : !!
1514 : : !! INPUTS
1515 : : !! mpsang=1+ maximum l quantum number
1516 : : !! xx= input value
1517 : : !!
1518 : : !! OUTPUT
1519 : : !! pl_d2(mpsang*mpsang)
1520 : : !!
1521 : : !! SOURCE
1522 : :
1523 : 1583 : subroutine pl_deriv(mpsang,pl_d2,xx)
1524 : :
1525 : : !Arguments ---------------------------------------------
1526 : : !scalars
1527 : : integer,intent(in) :: mpsang
1528 : : real(dp),intent(in) :: xx
1529 : : !arrays
1530 : : real(dp),intent(out) :: pl_d2(mpsang)
1531 : :
1532 : : !Local variables ---------------------------------------
1533 : : !scalars
1534 : : integer :: il,ilm
1535 : : real(dp) :: il_,il_m1,il_2m1
1536 : : character(len=500) :: msg
1537 : : !arrays
1538 : 3166 : real(dp) :: pl(mpsang),pl_d1(mpsang)
1539 : :
1540 : : ! *********************************************************************
1541 : :
1542 [ - + ]: 1583 : if (abs(xx).gt.1.d0) then
1543 : 0 : msg = 'pl_deriv : xx > 1 !'
1544 : 0 : LIBPAW_ERROR(msg)
1545 : : end if
1546 : :
1547 [ + + + + : 18300 : pl_d2=zero; pl_d1=zero; pl=zero
+ + ]
1548 : 1583 : pl(1)=one; pl(2)=xx
1549 : 1583 : pl_d1(1)=zero; pl_d1(2)=one
1550 : 1583 : pl_d2(1)=zero; pl_d2(2)=zero
1551 [ + + ]: 1583 : if (mpsang>2) then
1552 [ + + ]: 2434 : do il=2,mpsang-1
1553 : 1351 : il_=dble(il);il_m1=dble(il-1);il_2m1=dble(2*il-1)
1554 : 1351 : ilm=il+1
1555 : 1351 : pl(ilm)=(il_2m1*xx*pl(ilm-1)-il_m1*pl(ilm-2))/il_
1556 : 1351 : pl_d1(ilm)=(il_2m1*(xx*pl_d1(ilm-1)+pl(ilm-1))-il_m1*pl_d1(ilm-2))/il_
1557 : 2434 : pl_d2(ilm)=(il_2m1*(xx*pl_d2(ilm-1)+two*pl_d1(ilm-1))-il_m1*pl_d2(ilm-2))/il_
1558 : : end do
1559 : : end if
1560 : :
1561 : 1583 : end subroutine pl_deriv
1562 : : !!***
1563 : :
1564 : : !----------------------------------------------------------------------
1565 : :
1566 : : !!****f* m_paw_sphharm/ylm_angular_mesh
1567 : : !! NAME
1568 : : !! ylm_angular_mesh
1569 : : !!
1570 : : !! FUNCTION
1571 : : !! Build (theta, phi) angular mesh from (ntheta, nphi)
1572 : : !!
1573 : : !! INPUTS
1574 : : !! ntheta= number of sample points in the theta dir
1575 : : !! nphi= number of sample points in the phi dir
1576 : : !!
1577 : : !! OUTPUT
1578 : : !! angl_size= total number of sample points in the angular mesh, i.e. (ntheta * nphi)
1579 : : !! cart_coord(3, angl_size)= for each point of the angular mesh, gives the Cartesian coordinates
1580 : : !! of the corresponding point on an unitary sphere.
1581 : : !! ang_wgth(angl_size)= for each point of the angular mesh, gives the weight
1582 : : !! of the corresponding point on an unitary sphere.
1583 : : !!
1584 : : !! NOTE
1585 : : !! Summing over f * angwgth gives the spherical average 1/(4pi) \int domega f(omega)
1586 : : !!
1587 : : !! SOURCE
1588 : :
1589 : 63 : subroutine ylm_angular_mesh(ntheta, nphi, angl_size, cart_coord, ang_wgth)
1590 : :
1591 : : !Arguments ------------------------------------
1592 : : integer,intent(in) :: ntheta, nphi
1593 : : integer,intent(out) :: angl_size
1594 : : real(dp),allocatable,intent(out) :: cart_coord(:,:)
1595 : : real(dp),allocatable,intent(out) :: ang_wgth(:)
1596 : :
1597 : : !Local variables ------------------------------
1598 : : !scalars
1599 : : integer :: it, ip, npoints
1600 : : real(dp) :: ang, con, cos_phi, cos_theta, sin_phi, sin_theta
1601 : : character(len=500) :: msg
1602 : : !arrays
1603 : 63 : real(dp),allocatable :: th(:),wth(:)
1604 : :
1605 : : ! *************************************************************************
1606 : :
1607 [ - + + - : 189 : LIBPAW_ALLOCATE(th, (ntheta))
+ - + - -
+ - + ]
1608 [ - + + - : 126 : LIBPAW_ALLOCATE(wth, (ntheta))
- + - + ]
1609 : :
1610 : 63 : con = two_pi / nphi
1611 : 63 : call gauleg(-one, one, th, wth, ntheta)
1612 : :
1613 : 63 : angl_size = ntheta * nphi
1614 [ - + + - : 189 : LIBPAW_ALLOCATE(cart_coord, (3, angl_size))
+ - + - -
+ - + -
+ ]
1615 [ - + + - : 189 : LIBPAW_ALLOCATE(ang_wgth, (angl_size))
+ - + - -
+ - + -
+ ]
1616 : 63 : npoints = 0
1617 [ + + ]: 1099 : do it = 1, ntheta
1618 : 1036 : cos_theta = th(it)
1619 : 1036 : sin_theta = sqrt(one - cos_theta*cos_theta)
1620 [ + + ]: 32125 : do ip = 1, nphi
1621 : 31026 : ang = con * (ip-1)
1622 : 31026 : cos_phi = cos(ang); sin_phi = sin(ang)
1623 : 31026 : npoints = npoints + 1
1624 : 31026 : cart_coord(1, npoints) = sin_theta * cos_phi
1625 : 31026 : cart_coord(2, npoints) = sin_theta * sin_phi
1626 : 31026 : cart_coord(3, npoints) = cos_theta
1627 : : ! Normalization required
1628 : 32062 : ang_wgth(npoints) = wth(it) / (two * nphi)
1629 : : end do
1630 : : end do
1631 : :
1632 [ - + ]: 63 : LIBPAW_DEALLOCATE(th)
1633 [ - + ]: 63 : LIBPAW_DEALLOCATE(wth)
1634 : :
1635 : : !Error if npoints exceeds angl_size
1636 [ - + ]: 63 : if (npoints > angl_size) then
1637 : 0 : write(msg, '(a,i4,a,a,i4)' ) 'npoints =',npoints,ch10,&
1638 : 0 : & 'angl_size =',angl_size
1639 : 0 : LIBPAW_BUG(msg)
1640 : : end if
1641 : :
1642 : 63 : end subroutine ylm_angular_mesh
1643 : : !!***
1644 : :
1645 : : !----------------------------------------------------------------------
1646 : :
1647 : : !!****f* m_paw_sphharm/mat_mlms2jmj
1648 : : !! NAME
1649 : : !! mat_mlms2jmj
1650 : : !!
1651 : : !! FUNCTION
1652 : : !! For a given angular momentum lcor, change a matrix of dimension 2(2*lcor+1)
1653 : : !! from the Ylm basis to the J,M_J basis if option==1
1654 : : !!
1655 : : !! INPUTS
1656 : : !! lcor= angular momentum
1657 : : !! ndij= ndij = 4
1658 : : !! option= 1 matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis
1659 : : !! 2 matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis
1660 : : !! optspin= 1 Spin up are first
1661 : : !! 2 Spin dn are first
1662 : : !! prtvol=printing volume
1663 : : !! unitfi=printing file unit ; -1 for no printing
1664 : : !! wrt_mode=printing mode in parallel ('COLL' or 'PERS')
1665 : : !!
1666 : : !! SIDE EFFECTS
1667 : : !! mat_mlms= Input/Output matrix in the Ylm basis, size of the matrix is (2*lcor+1,2*lcor+1,ndij)
1668 : : !! mat_jmj= Input/Output matrix in the J,M_J basis, size is 2*(2*lcor+1),2*(2*lcor+1)
1669 : : !!
1670 : : !! NOTES
1671 : : !! usefull only in ndij==4
1672 : : !!
1673 : : !! SOURCE
1674 : :
1675 : 863 : subroutine mat_mlms2jmj(lcor,mat_mlms,mat_jmj,ndij,option,optspin,prtvol,unitfi,wrt_mode)
1676 : :
1677 : : !Arguments ---------------------------------------------
1678 : : !scalars
1679 : : integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi
1680 : : character(len=4),intent(in) :: wrt_mode
1681 : : !arrays
1682 : : complex(dp),intent(inout) :: mat_mlms(2*lcor+1,2*lcor+1,ndij)
1683 : : complex(dp),intent(inout) :: mat_jmj(2*(2*lcor+1),2*(2*lcor+1))
1684 : :
1685 : : !Local variables ---------------------------------------
1686 : : !scalars
1687 : : integer :: ii,im,im1,im2,ispden,jc1,jc2,jj,jm,ll,ml1,ml2,ms1,ms2
1688 : : real(dp),parameter :: invsqrt2=one/sqrt2
1689 : : real(dp) :: invsqrt2lp1,xj,xmj
1690 : : complex(dp) :: mat_tmp,tmp2
1691 : : character(len=9),parameter :: dspinold(6)=(/"up ","down ","up-up ","down-down","up-dn ","dn-up "/)
1692 : : character(len=9),parameter :: dspin(6)=(/"dn ","up ","dn-dn ","up-up ","dn-up ","up-dn "/)
1693 : : character(len=500) :: msg
1694 : : !arrays
1695 : 863 : integer, allocatable :: ind_msml(:,:)
1696 : 863 : complex(dp),allocatable :: mat_mlms2(:,:),mlms2jmj(:,:)
1697 : : !*********************************************************************
1698 : :
1699 [ - + ]: 863 : if(ndij/=4) then
1700 : 0 : msg=" ndij/=4 !"
1701 : 0 : LIBPAW_BUG(msg)
1702 : : end if
1703 [ - + ]: 863 : if (option/=1.and.option/=2) then
1704 : 0 : msg=' option=/1 and =/2 !'
1705 : 0 : LIBPAW_BUG(msg)
1706 : : end if
1707 [ - + ]: 863 : if (optspin/=1.and.optspin/=2) then
1708 : 0 : msg=' optspin=/1 and =/2 !'
1709 : 0 : LIBPAW_BUG(msg)
1710 : : end if
1711 : :
1712 [ + + ]: 863 : if (unitfi/=-1) then
1713 [ + - ]: 801 : if(option==1) then
1714 : 801 : write(msg,'(3a)') ch10,&
1715 : 1602 : & "matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis"
1716 : 801 : call wrtout(unitfi,msg,wrt_mode)
1717 [ # # ]: 0 : else if(option==2) then
1718 : 0 : write(msg,'(3a)') ch10,&
1719 : 0 : & "matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis"
1720 : 0 : call wrtout(unitfi,msg,wrt_mode)
1721 : : end if
1722 : : end if
1723 : :
1724 [ + - ]: 863 : if(option==1) then
1725 [ - + ]: 863 : if(optspin==2) then
1726 [ # # # # ]: 0 : if(abs(prtvol)>2.and.unitfi/=-1)&
1727 : 0 : & write(msg,'(3a)') ch10,"assume spin dn is the first in the array"
1728 [ + - ]: 863 : else if (optspin==1) then
1729 [ + + + + ]: 863 : if(abs(prtvol)>2.and.unitfi/=-1)&
1730 : 783 : & write(msg,'(3a)') ch10,"change array in order that spin dn is the first in the array"
1731 [ + + ]: 5250 : do ii=1,2*lcor+1
1732 [ + + ]: 27689 : do jj=1,2*lcor+1
1733 : 22439 : mat_tmp=mat_mlms(ii,jj,2)
1734 : 22439 : mat_mlms(ii,jj,2)=mat_mlms(ii,jj,1)
1735 : 22439 : mat_mlms(ii,jj,1)=mat_tmp
1736 : 22439 : mat_tmp=mat_mlms(ii,jj,4)
1737 : 22439 : mat_mlms(ii,jj,4)=mat_mlms(ii,jj,3)
1738 : 26826 : mat_mlms(ii,jj,3)=mat_tmp
1739 : : end do
1740 : : end do
1741 : : ! mat_tmp(:,:,1)=mat_mlms(:,:,2);mat_tmp(:,:,2)=mat_mlms(:,:,1)
1742 : : ! mat_tmp(:,:,3)=mat_mlms(:,:,4);mat_tmp(:,:,4)=mat_mlms(:,:,3)
1743 : : ! mat_mlms(:,:,:)=mat_tmp(:,:,:)
1744 : : end if
1745 [ + + + + ]: 863 : if(abs(prtvol)>2.and.unitfi/=-1) then
1746 : 783 : call wrtout(unitfi,msg,wrt_mode)
1747 : : end if
1748 : : end if
1749 : :
1750 [ + - + + : 863 : if(option==1.and.abs(prtvol)>2.and.unitfi/=-1) then
+ + ]
1751 [ + + ]: 3915 : do ispden=1,ndij
1752 : 3132 : write(msg,'(3a)') ch10,&
1753 [ + - ]: 6264 : & "Input matrix in the Ylm basis for component ",trim(dspin(ispden+2*(ndij/4)))
1754 : 3132 : call wrtout(unitfi,msg,wrt_mode)
1755 [ + + ]: 19831 : do im1=1,lcor*2+1
1756 : : write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
1757 : 15916 : & (mat_mlms(im1,im2,ispden),im2=1,lcor*2+1)
1758 : 19048 : call wrtout(unitfi,msg,wrt_mode)
1759 : : end do
1760 : : end do
1761 : : end if ! option==1
1762 : :
1763 : : !--------------- Built indices + allocations
1764 : 863 : ll=lcor
1765 [ + - - + : 4315 : LIBPAW_ALLOCATE(mlms2jmj,(2*(2*ll+1),2*(2*ll+1)))
+ - - + +
- + - + -
- + - + ]
1766 [ + + + + ]: 99393 : mlms2jmj=czero
1767 [ - + + - : 2589 : LIBPAW_BOUND2_ALLOCATE(ind_msml,BOUNDS(1,2),BOUNDS(-ll,ll))
+ - + - -
+ - + ]
1768 [ - + + - : 2589 : LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1)))
- + + - -
+ - + ]
1769 [ + + + + ]: 99393 : mlms2jmj=czero
1770 : : jc1=0
1771 [ + + ]: 2589 : do ms1=1,2
1772 [ + + ]: 11363 : do ml1=-ll,ll
1773 : 8774 : jc1=jc1+1
1774 : 10500 : ind_msml(ms1,ml1)=jc1
1775 : : end do
1776 : : end do
1777 : : !--------------- Change representation of input matrix for ndij==4
1778 [ + - ]: 863 : if(option==1) then
1779 : : jc1=0
1780 [ + + ]: 2589 : do ms1=1,2
1781 [ + + ]: 11363 : do ml1=1,2*ll+1
1782 : 8774 : jc1=jc1+1
1783 : 8774 : jc2=0
1784 [ + + ]: 28048 : do ms2=1,2
1785 [ + + ]: 116078 : do ml2=1,2*ll+1
1786 : 89756 : jc2=jc2+1
1787 [ + + ]: 89756 : if(ms1==ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,ms1)
1788 [ + + ]: 89756 : if(ms1<ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,3)
1789 [ + + ]: 107304 : if(ms1>ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,4)
1790 : : end do
1791 : : end do
1792 : : end do
1793 : : end do
1794 [ + + + + ]: 863 : if(abs(prtvol)>1.and.unitfi/=-1) then
1795 : 783 : write(msg,'(3a)') ch10,"Input matrix in the lms basis for all component"
1796 : 783 : call wrtout(unitfi,msg,wrt_mode)
1797 [ + + ]: 8741 : do im1=1,2*(lcor*2+1)
1798 : : write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
1799 : 7958 : & (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1))
1800 : 8741 : call wrtout(unitfi,msg,wrt_mode)
1801 : : end do
1802 : : end if
1803 : : end if ! option==1
1804 : :
1805 : : !--------------- built mlms2jmj
1806 : : !do jj=ll,ll+1 ! the physical value of j are ll-0.5,ll+0.5
1807 : : !xj(jj)=jj-0.5
1808 [ - + ]: 863 : if(ll==0)then
1809 : 0 : msg=' ll should not be equal to zero !'
1810 : 0 : LIBPAW_BUG(msg)
1811 : : end if
1812 : 863 : jc1=0
1813 : 863 : invsqrt2lp1=one/sqrt(float(2*lcor+1))
1814 [ + + ]: 2589 : do jj=ll,ll+1
1815 : 1726 : xj=float(jj)-half
1816 [ + + ]: 11363 : do jm=-jj,jj-1
1817 : 8774 : xmj=float(jm)+half
1818 : 8774 : jc1=jc1+1
1819 [ + + ]: 8774 : if(nint(xj+0.5)==ll+1) then
1820 [ + + ]: 5250 : if(nint(xmj+0.5)==ll+1) then
1821 : 863 : mlms2jmj(ind_msml(2,ll),jc1)=1.0 ! J=L+0.5 and m_J=L+0.5
1822 [ + + ]: 4387 : else if(nint(xmj-0.5)==-ll-1) then
1823 : 863 : mlms2jmj(ind_msml(1,-ll),jc1)=1.0 ! J=L+0.5 and m_J=-L-0.5
1824 : : else
1825 : 3524 : mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1826 : 3524 : mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1827 : : end if
1828 : : end if
1829 [ + + ]: 10500 : if(nint(xj+0.5)==ll) then
1830 : 3524 : mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1831 : 3524 : mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1832 : : end if
1833 : : end do
1834 : : end do
1835 [ + + + + ]: 863 : if(abs(prtvol)>2.and.unitfi/=-1) then
1836 : 783 : write(msg,'(3a)') ch10,"Matrix to go from |M_L,M_S> to |J,M_J>"
1837 : 783 : call wrtout(unitfi,msg,wrt_mode)
1838 [ + + ]: 8741 : do im1=1,2*(lcor*2+1)
1839 : 7958 : write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mlms2jmj(im1,im2),im2=1,2*(lcor*2+1))
1840 : 8741 : call wrtout(unitfi,msg,wrt_mode)
1841 : : end do
1842 : : end if
1843 : :
1844 [ + + ]: 9637 : do jm=1,2*(2*ll+1)
1845 [ + + ]: 99393 : do im=1,2*(2*ll+1)
1846 : : tmp2=czero
1847 [ + + ]: 1015540 : do ii=1,2*(2*ll+1)
1848 [ + + ]: 10668516 : do jj=1,2*(2*ll+1)
1849 [ + - ]: 10578760 : if(option==1) then
1850 : 9652976 : tmp2=tmp2+mat_mlms2(ii,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))
1851 [ # # ]: 0 : else if(option==2) then
1852 : 0 : tmp2=tmp2+mat_jmj(ii,jj)*(mlms2jmj(im,ii))*CONJG(mlms2jmj(jm,jj)) ! inv=t*
1853 : : end if
1854 : : end do
1855 : : end do
1856 [ + - ]: 98530 : if(option==1) then
1857 : 89756 : mat_jmj(im,jm)=tmp2
1858 [ # # ]: 0 : else if(option==2) then
1859 : 0 : mat_mlms2(im,jm)=tmp2
1860 : : end if
1861 : : end do
1862 : : end do
1863 [ + - ]: 863 : if(option==1) then
1864 [ + + + + ]: 863 : if (abs(prtvol)>=1.and.unitfi/=-1) then
1865 : 783 : write(msg,'(3a)') ch10," Matrix in the J,M_J basis"
1866 : 783 : call wrtout(unitfi,msg,wrt_mode)
1867 [ + + ]: 8741 : do im1=1,2*(lcor*2+1)
1868 : 7958 : write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_jmj(im1,im2),im2=1,2*(lcor*2+1))
1869 : 8741 : call wrtout(unitfi,msg,wrt_mode)
1870 : : end do
1871 : : end if
1872 [ # # ]: 0 : else if(option==2) then
1873 [ # # # # ]: 0 : if (abs(prtvol)>=1.and.unitfi/=-1) then
1874 : 0 : write(msg,'(3a)') ch10," Matrix in the m_s m_l basis"
1875 : 0 : call wrtout(unitfi,msg,wrt_mode)
1876 [ # # ]: 0 : do im1=1,2*(lcor*2+1)
1877 : 0 : write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1))
1878 : 0 : call wrtout(unitfi,msg,wrt_mode)
1879 : : end do
1880 : : end if
1881 : : jc1=0
1882 [ # # ]: 0 : do ms1=1,2
1883 [ # # ]: 0 : do ml1=1,2*ll+1
1884 : 0 : jc1=jc1+1
1885 : 0 : jc2=0
1886 [ # # ]: 0 : do ms2=1,2
1887 [ # # ]: 0 : do ml2=1,2*ll+1
1888 : 0 : jc2=jc2+1
1889 [ # # ]: 0 : if(ms1==ms2) mat_mlms(ml1,ml2,ms1)=mat_mlms2(jc1,jc2)
1890 [ # # ]: 0 : if(ms1<ms2) mat_mlms(ml1,ml2,3)=mat_mlms2(jc1,jc2)
1891 [ # # ]: 0 : if(ms1>ms2) mat_mlms(ml1,ml2,4)=mat_mlms2(jc1,jc2)
1892 : : end do
1893 : : end do
1894 : : end do
1895 : : end do
1896 : : end if
1897 [ - + ]: 863 : LIBPAW_DEALLOCATE(mlms2jmj)
1898 [ - + ]: 863 : LIBPAW_DEALLOCATE(mat_mlms2)
1899 [ - + ]: 863 : LIBPAW_DEALLOCATE(ind_msml)
1900 : :
1901 : 863 : end subroutine mat_mlms2jmj
1902 : : !!***
1903 : :
1904 : : !----------------------------------------------------------------------
1905 : :
1906 : : !!****f* m_paw_sphharm/mat_slm2ylm
1907 : : !! NAME
1908 : : !! mat_slm2ylm
1909 : : !!
1910 : : !! FUNCTION
1911 : : !! For a given angular momentum lcor, change a matrix of dimension (2*lcor+1)
1912 : : !! from the Slm to the Ylm basis if option==1 or from Ylm to Slm if !option==2
1913 : : !!
1914 : : !! INPUTS
1915 : : !! lcor= angular momentum, size of the matrix is 2(2*lcor+1)
1916 : : !! mat_inp_c= Input matrix
1917 : : !! ndij= ndij = 4
1918 : : !! option= -1 Change matrix from Slm to Ylm basis
1919 : : !! 1 Change matrix from Ylm to Slm basis
1920 : : !! optspin= 1 Spin up are first
1921 : : !! 2 Spin dn are first
1922 : : !! prtvol=printing volume
1923 : : !! unitfi=printing file unit ; -1 for no printing
1924 : : !! wrt_mode=printing mode in parallel ('COLL' or 'PERS')
1925 : : !!
1926 : : !! OUTPUT
1927 : : !! mat_inp_c= Output matrix in Ylm or Slm basis according to option
1928 : : !!
1929 : : !! NOTES
1930 : : !! usefull only in ndij==4
1931 : : !!
1932 : : !! SOURCE
1933 : :
1934 : 871 : subroutine mat_slm2ylm(lcor,mat_inp_c,mat_out_c,ndij,option,optspin,prtvol,unitfi,wrt_mode)
1935 : :
1936 : : !Arguments ---------------------------------------------
1937 : : !scalars
1938 : : integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi
1939 : : character(len=4),intent(in) :: wrt_mode
1940 : : !arrays
1941 : 1742 : complex(dp) :: mat_inp_c(2*lcor+1,2*lcor+1,ndij),mat_out(2*lcor+1,2*lcor+1,ndij)
1942 : : complex(dp) :: mat_out_c(2*lcor+1,2*lcor+1,ndij)
1943 : :
1944 : : !Local variables ---------------------------------------
1945 : : !scalars
1946 : : integer :: jm,ii,jj,ll,mm,ispden,im,im1,im2
1947 : : real(dp),parameter :: invsqrt2=one/sqrt2
1948 : : real(dp) :: onem
1949 : : complex(dp) :: tmp2
1950 : : character(len=9),parameter :: dspinc(6)=(/"up ","down ","up-up ","down-down","up-dn ","dn-up "/)! optspin 1
1951 : : character(len=9),parameter :: dspinc2(6)=(/"up ","down ","dn-dn ","up-up ","dn-up ","up-dn "/)! optspin 2
1952 : : character(len=500) :: msg
1953 : : !arrays
1954 : 871 : complex(dp),allocatable :: slm2ylm(:,:)
1955 : : ! *********************************************************************
1956 : :
1957 [ - + ]: 871 : if(ndij/=4) then
1958 : 0 : msg=' ndij:=4 !'
1959 : 0 : LIBPAW_BUG(msg)
1960 : : end if
1961 [ - + ]: 871 : if (option/=1.and.option/=2.and.option/=3.and.option/=4) then
1962 : 0 : msg=' option=/1 or 2 or 3 or 4 !'
1963 : 0 : LIBPAW_BUG(msg)
1964 : : end if
1965 : :
1966 [ + + + + ]: 871 : if(abs(prtvol)>2.and.unitfi/=-1) then
1967 : 816 : write(msg,'(3a)') ch10, " mat_slm2ylm"
1968 : 816 : call wrtout(unitfi,msg,wrt_mode)
1969 : : end if
1970 : :
1971 [ + + + + ]: 871 : if(abs(prtvol)>2.and.unitfi/=-1) then
1972 [ + - ]: 816 : if(option==1.or.option==3) then
1973 : 816 : write(msg,'(3a)') ch10,"matrix in Slm basis is changed into Ylm basis"
1974 : 816 : call wrtout(unitfi,msg,wrt_mode)
1975 [ # # ]: 0 : else if(option==2.or.option==4) then
1976 : 0 : write(msg,'(3a)') ch10,"matrix in Ylm basis is changed into Slm basis"
1977 : 0 : call wrtout(unitfi,msg,wrt_mode)
1978 : : end if
1979 : : end if
1980 : :
1981 : 871 : ll=lcor
1982 [ + - - + : 4355 : LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1))
+ - - + +
- + - + -
- + + - ]
1983 [ + + + + ]: 27937 : slm2ylm=czero
1984 [ + + + + : 112619 : mat_out=zero
+ + ]
1985 [ + + + + : 112619 : mat_out_c=czero
+ + ]
1986 [ + + ]: 5298 : do im=1,2*ll+1
1987 : 4427 : mm=im-ll-1;jm=-mm+ll+1
1988 : 4427 : onem=dble((-1)**mm)
1989 [ + + ]: 4427 : if (mm> 0) then
1990 : 1778 : slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
1991 : 1778 : slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
1992 : : end if
1993 [ + + ]: 4427 : if (mm==0) then
1994 : 871 : slm2ylm(im,im)=cone
1995 : : end if
1996 [ + + ]: 5298 : if (mm< 0) then
1997 : 1778 : slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
1998 : 1778 : slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
1999 : : end if
2000 : : end do
2001 [ + + + + ]: 871 : if(abs(prtvol)>2.and.unitfi/=-1) then
2002 [ + + ]: 4080 : do ispden=1,ndij
2003 [ + - ]: 3264 : if(optspin==1) then
2004 [ + - ]: 3264 : if(option==1.or.option==3)&
2005 : 3264 : & write(msg,'(3a)') ch10,&
2006 [ + - ]: 6528 : & "Input matrix in the Slm basis for component ",trim(dspinc(ispden+2*(ndij/4)))
2007 [ - + ]: 3264 : if(option==2.or.option==3)&
2008 : 0 : & write(msg,'(3a)') ch10,&
2009 [ # # ]: 0 : & "Input matrix in the Ylm basis for component ",trim(dspinc(ispden+2*(ndij/4)))
2010 : : else
2011 [ # # ]: 0 : if(option==1.or.option==3)&
2012 : 0 : & write(msg,'(3a)') ch10,&
2013 [ # # ]: 0 : & "Input matrix in the Slm basis for component ",trim(dspinc2(ispden+2*(ndij/4)))
2014 [ # # ]: 0 : if(option==2.or.option==3)&
2015 : 0 : & write(msg,'(3a)') ch10,&
2016 [ # # ]: 0 : & "Input matrix in the Ylm basis for component ",trim(dspinc2(ispden+2*(ndij/4)))
2017 : : end if
2018 : 3264 : call wrtout(unitfi,msg,wrt_mode)
2019 [ + + ]: 20672 : do im1=1,lcor*2+1
2020 : : write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2021 : 16592 : & (mat_inp_c(im1,im2,ispden),im2=1,lcor*2+1)
2022 : 19856 : call wrtout(unitfi,msg,wrt_mode)
2023 : : end do
2024 : : end do
2025 : : end if
2026 [ + + ]: 4355 : do ispden=1,ndij
2027 [ + + ]: 22063 : do jm=1,2*ll+1
2028 [ + + ]: 111748 : do im=1,2*ll+1
2029 : : tmp2=czero
2030 [ + + ]: 557448 : do ii=1,2*ll+1
2031 [ + + ]: 2990692 : do jj=1,2*ll+1
2032 [ + - ]: 2900136 : if(option==1) then
2033 : 2433244 : tmp2=tmp2+mat_inp_c(ii,jj,ispden)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))
2034 [ # # ]: 0 : else if(option==2) then
2035 : 0 : tmp2=tmp2+mat_inp_c(ii,jj,ispden)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))
2036 : : end if
2037 : : end do
2038 : : end do
2039 : 108264 : mat_out_c(im,jm,ispden)=tmp2
2040 : : end do
2041 : : end do
2042 : : end do ! ispden
2043 [ + + ]: 5298 : do ii=1,2*ll+1
2044 [ + + ]: 27937 : do jj=1,2*ll+1
2045 : 22639 : mat_out(ii,jj,1)=real(mat_out_c(ii,jj,1))
2046 : 22639 : mat_out(ii,jj,2)=real(mat_out_c(ii,jj,2))
2047 : 22639 : mat_out(ii,jj,3)=real(mat_out_c(ii,jj,3))
2048 : 22639 : mat_out(ii,jj,4)=aimag(mat_out_c(ii,jj,3))
2049 : : ! check that n_{m,m'}^{alpha,beta}=conjg(n_{m',m"}^{beta,alpha}).
2050 [ + - - + ]: 22639 : if((abs(aimag(mat_out_c(ii,jj,3))+aimag(mat_out_c(jj,ii,4))).ge.0.0001).or. &
2051 : 4427 : & (abs(real(mat_out_c(ii,jj,3))-real(mat_out_c(jj,ii,4))).ge.0.0001)) then
2052 : : write(msg,'(a,4f10.4)') &
2053 : 0 : & ' prb with mat_out_c ',mat_out_c(ii,jj,3),mat_out_c(ii,jj,4)
2054 : 0 : LIBPAW_BUG(msg)
2055 : : end if
2056 : : end do
2057 : : end do
2058 : :
2059 [ - + ]: 871 : LIBPAW_DEALLOCATE(slm2ylm)
2060 : :
2061 : 871 : end subroutine mat_slm2ylm
2062 : : !!***
2063 : :
2064 : : !----------------------------------------------------------------------
2065 : :
2066 : : !!****f* m_paw_sphharm/setsym_ylm
2067 : : !! NAME
2068 : : !! setsym_ylm
2069 : : !!
2070 : : !! FUNCTION
2071 : : !! Compute rotation matrices expressed in the basis of real spherical harmonics
2072 : : !! This coefficients are used later to symmetrize PAW on-site quantities (rhoij, dij, ...).
2073 : : !!
2074 : : !! INPUTS
2075 : : !! gprimd(3,3)==dimensional primitive translations for reciprocal space (bohr^-1)
2076 : : !! lmax=value of lmax mentioned at the second line of the psp file
2077 : : !! nsym=number of symmetry elements in space group
2078 : : !! pawprtvol=control print volume and debugging output
2079 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2080 : : !! sym(3,3,nsym)=symmetries of group in terms of operations on primitive translations
2081 : : !!
2082 : : !! OUTPUT
2083 : : !! zarot(2*lmax+1,2*lmax+1,lmax+1,nsym)=coefficients of the
2084 : : !! transformation of real spherical harmonics
2085 : : !! under the symmetry operations
2086 : : !!
2087 : : !! NOTES
2088 : : !! Typical use: sym(:,:,:) is symrec(:,:,:) (rotations in reciprocal space)
2089 : : !! because we need symrel^-1 (=transpose[symrec])
2090 : : !! to symmetrize quantities.
2091 : : !!
2092 : : !! - This file comes from the file crystal_symmetry.f
2093 : : !! by N.A.W. Holzwarth and A. Tackett for the code pwpaw
2094 : : !! - Uses sign & phase convension of M. E. Rose, Elementary Theory of Angular
2095 : : !! Momentum, John Wiley & Sons,. inc. 1957)
2096 : : !! zalpha = exp(-i*alpha) zgamma = exp (-i*gamma)
2097 : : !! - Assumes each transformation can be expressed in terms of 3 Euler
2098 : : !! angles with or without inversion
2099 : : !!
2100 : : !! Reference for evaluation of rotation matrices in the basis of real SH:
2101 : : !! Blanco M.A., Florez M. and Bermejo M.
2102 : : !! Journal of Molecular Structure: THEOCHEM, Volume 419, Number 1, 8 December 1997 , pp. 19-27(9)
2103 : : !! http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf
2104 : : !!
2105 : : !! SOURCE
2106 : :
2107 : 1979 : subroutine setsym_ylm(gprimd,lmax,nsym,pawprtvol,rprimd,sym,zarot)
2108 : :
2109 : : !Arguments ---------------------------------------------
2110 : : !scalars
2111 : : integer,intent(in) :: lmax,nsym,pawprtvol
2112 : : !arrays
2113 : : integer,intent(in) :: sym(3,3,nsym)
2114 : : real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
2115 : : real(dp),intent(out) :: zarot(2*lmax+1,2*lmax+1,lmax+1,nsym)
2116 : :
2117 : : !Local variables ------------------------------
2118 : : !scalars
2119 : : integer :: i1,ii,il,irot,isn,j1,jj,k1,ll,mm,mp
2120 : : real(dp) :: cosalp,cosbeta,cosgam,sinalp,sinbeta,singam
2121 : : character(len=1000) :: msg
2122 : : !arrays
2123 : : real(dp) :: prod(3,3),rot(3,3)
2124 : :
2125 : : !************************************************************************
2126 : :
2127 [ + + ]: 1979 : if (abs(pawprtvol)>=3) then
2128 : 101 : write(msg,'(8a,i4)') ch10,&
2129 : 101 : & ' PAW TEST:',ch10,&
2130 : 101 : & ' ==== setsym_ylm: rotation matrices in the basis ============',ch10,&
2131 : 101 : & ' ==== of real spherical harmonics ============',ch10,&
2132 : 202 : & ' > Number of symmetries (nsym)=',nsym
2133 : 101 : call wrtout(std_out,msg,'COLL')
2134 : : end if
2135 : :
2136 [ + + + + : 1697687 : zarot=zero
+ + + + ]
2137 : :
2138 [ + + ]: 30098 : do irot=1,nsym
2139 : :
2140 [ + + ]: 28119 : if (abs(pawprtvol)>=3) then
2141 : 3620 : write(msg,'(a,i2,a,9i2,a)') ' >For symmetry ',irot,' (',sym(:,:,irot),')'
2142 : 3620 : call wrtout(std_out,msg,'COLL')
2143 : : end if
2144 : :
2145 : : ! === l=0 case ===
2146 : 28119 : zarot(1,1,1,irot)=one
2147 : :
2148 : : ! === l>0 case ===
2149 [ + + ]: 28119 : if (lmax>0) then
2150 : : ! Calculate the rotations in the cartesian basis
2151 : 27816 : rot=zero;prod=zero
2152 [ + + ]: 111264 : do k1=1,3
2153 [ + + ]: 361608 : do j1=1,3
2154 [ + + ]: 1084824 : do i1=1,3
2155 : 1001376 : prod(i1,j1)=prod(i1,j1)+sym(i1,k1,irot)*rprimd(j1,k1)
2156 : : end do
2157 : : end do
2158 : : end do
2159 [ + + ]: 111264 : do j1=1,3
2160 [ + + ]: 361608 : do i1=1,3
2161 [ + + ]: 1001376 : do k1=1,3
2162 : 1001376 : rot(i1,j1)=rot(i1,j1)+gprimd(i1,k1)*prod(k1,j1)
2163 : : end do
2164 [ + + ]: 333792 : if(abs(rot(i1,j1))<tol10) rot(i1,j1)=zero
2165 : : end do
2166 : : end do
2167 : 27816 : call mkeuler(rot,cosbeta,sinbeta,cosalp,sinalp,cosgam,singam,isn)
2168 [ + + ]: 68688 : do ll=1,lmax
2169 : 40872 : il=(isn)**ll
2170 [ + + ]: 219420 : do mp=-ll,ll
2171 : 150732 : jj=mp+ll+1
2172 [ + + ]: 792588 : do mm=-ll,ll
2173 : 600984 : ii=mm+ll+1
2174 : :
2175 : : ! Formula (47) from the paper of Blanco et al
2176 : : zarot(ii,jj,ll+1,irot)=il&
2177 : : & *(phim(cosalp,sinalp,mm)*phim(cosgam,singam,mp)*sign(1,mp)&
2178 : : *(dbeta(cosbeta,sinbeta,ll,abs(mp),abs(mm))&
2179 : : & +(-1._dp)**mm*dbeta(cosbeta,sinbeta,ll,abs(mm),-abs(mp)))*half&
2180 : : & -phim(cosalp,sinalp,-mm)*phim(cosgam,singam,-mp)*sign(1,mm)&
2181 : : *(dbeta(cosbeta,sinbeta,ll,abs(mp),abs(mm))&
2182 : 751716 : & -(-1._dp)**mm*dbeta(cosbeta,sinbeta,ll,abs(mm),-abs(mp)))*half)
2183 : : end do
2184 : : end do
2185 : : end do
2186 : : end if ! lmax case
2187 : :
2188 [ + + ]: 30098 : if (abs(pawprtvol)>=3) then
2189 [ + - ]: 3620 : if(lmax>0) then
2190 : : write(msg,'(2a,3(3(2x,f7.3),a))') &
2191 : 3620 : & ' Rotation matrice for l=1:',ch10,&
2192 : 3620 : & (zarot(1,jj,2,irot),jj=1,3),ch10,&
2193 : 3620 : & (zarot(2,jj,2,irot),jj=1,3),ch10,&
2194 : 7240 : & (zarot(3,jj,2,irot),jj=1,3)
2195 : 3620 : call wrtout(std_out,msg,'COLL')
2196 : : end if
2197 [ + + ]: 3620 : if(lmax>1) then
2198 : : write(msg,'(2a,5(5(2x,f7.3),a))') &
2199 : 3380 : & ' Rotation matrice for l=2:',ch10,&
2200 : 3380 : & (zarot(1,jj,3,irot),jj=1,5),ch10,&
2201 : 3380 : & (zarot(2,jj,3,irot),jj=1,5),ch10,&
2202 : 3380 : & (zarot(3,jj,3,irot),jj=1,5),ch10,&
2203 : 3380 : & (zarot(4,jj,3,irot),jj=1,5),ch10,&
2204 : 6760 : & (zarot(5,jj,3,irot),jj=1,5)
2205 : 3380 : call wrtout(std_out,msg,'COLL')
2206 : : end if
2207 [ + + ]: 3620 : if(lmax>2) then
2208 : : write(msg,'(2a,7(7(2x,f7.3),a))') &
2209 : 576 : & ' Rotation matrice for l=3:',ch10,&
2210 : 576 : & (zarot(1,jj,4,irot),jj=1,7),ch10,&
2211 : 576 : & (zarot(2,jj,4,irot),jj=1,7),ch10,&
2212 : 576 : & (zarot(3,jj,4,irot),jj=1,7),ch10,&
2213 : 576 : & (zarot(4,jj,4,irot),jj=1,7),ch10,&
2214 : 576 : & (zarot(5,jj,4,irot),jj=1,7),ch10,&
2215 : 576 : & (zarot(6,jj,4,irot),jj=1,7),ch10,&
2216 : 1152 : & (zarot(7,jj,4,irot),jj=1,7)
2217 : 576 : call wrtout(std_out,msg,'COLL')
2218 : : end if
2219 : : end if
2220 : :
2221 : : end do ! isym loop
2222 : :
2223 : 1979 : end subroutine setsym_ylm
2224 : : !!***
2225 : :
2226 : : !----------------------------------------------------------------------
2227 : :
2228 : : !!****f* m_paw_sphharm/setnabla_ylm
2229 : : !! NAME
2230 : : !! setnabla_ylm
2231 : : !!
2232 : : !! FUNCTION
2233 : : !! Evaluate several inegrals involving spherical harmonics and their gradient.
2234 : : !! These integrals are angular part for <phi_i|nabla|phi_j> and <tphi_i|nabla|tphi_j>.
2235 : : !!
2236 : : !! INPUTS
2237 : : !! mpsang=1+ max. angular momentum
2238 : : !!
2239 : : !! OUTPUT
2240 : : !! ang_phipphj :: angular part for <phi_i|nabla|phi_j> and <tphi_i|nabla|tphi_j>
2241 : : !! ang_phipphj(i,j,1)=\int sin\theta cos\phi Si Sj d\omega
2242 : : !! ang_phipphj(i,j,2)=\int cos\theta cos\phi Si \frac{d}{d\theta}Sj d\Omega
2243 : : !! ang_phipphj(i,j,3)=\int -sin\phi Si \frac{d}{d\phi}Sj d\Omega
2244 : : !! ang_phipphj(i,j,4)=\int sin\theta sin\phi Si Sj d\Omega
2245 : : !! ang_phipphj(i,j,5)=\int cos\theta sin\phi Si \frac{d}{d\theta}Sj d\Omega
2246 : : !! ang_phipphj(i,j,6)=\int cos\phi Si \frac{d}{d\phi}Sj d\Omega
2247 : : !! ang_phipphj(i,j,7)=\int cos\theta Si Sj d\Omega
2248 : : !! ang_phipphj(i,j,8)=\int -sin\theta Si \frac{d}{d\theta}Sj d\Omega
2249 : : !!
2250 : : !! NOTES
2251 : : !! See : Mazevet, S., Torrent, M., Recoules, V. and Jollet, F., High Energy Density Physics, 6, 84-88 (2010)
2252 : : !! Calculations of the Transport Properties within the PAW Formalism
2253 : : !! SOURCE
2254 : :
2255 : 1158 : subroutine setnabla_ylm(ang_phipphj,mpsang)
2256 : :
2257 : : !Arguments ------------------------------------
2258 : : !scalars
2259 : : integer,intent(in) :: mpsang
2260 : : !arrays
2261 : : real(dp),intent(out) :: ang_phipphj(mpsang**2,mpsang**2,8)
2262 : :
2263 : : !Local variables-------------------------------
2264 : : character(len=500) :: msg
2265 : : real(dp) :: ang_phipphj_tmp(16,16,8)
2266 : :
2267 : : ! ************************************************************************
2268 : :
2269 [ - + ]: 1158 : if (mpsang>4) then
2270 : 0 : msg=' Not designed for angular momentum greater than 3!'
2271 : 0 : LIBPAW_ERROR(msg)
2272 : : end if
2273 : :
2274 : : !8 angular integrals for l=0..3, m=-l..+l
2275 : : !ang_phipphj(1,4,1)=\frac{1}{\sqrt{3}}
2276 : : !ang_phipphj(2,5,1)=\frac{1}{\sqrt{5}}
2277 : : !ang_phipphj(3,8,1)=\frac{1}{\sqrt{5}}
2278 : : !ang_phipphj(4,1,1)=\frac{1}{\sqrt{3}}
2279 : : !ang_phipphj(4,7,1)=-\frac{1}{\sqrt{15}}
2280 : : !ang_phipphj(4,9,1)=\frac{1}{\sqrt{5}}
2281 : : !ang_phipphj(5,2,1)=\frac{1}{\sqrt{5}}
2282 : : !ang_phipphj(5,10,1)=\sqrt{\frac{3}{14}}
2283 : : !ang_phipphj(5,12,1)=-\frac{1}{\sqrt{70}}
2284 : : !ang_phipphj(6,11,1)=\frac{1}{\sqrt{7}}
2285 : : !ang_phipphj(7,4,1)=-\frac{1}{\sqrt{15}}
2286 : : !ang_phipphj(7,14,1)=\sqrt{\frac{6}{35}}
2287 : : !ang_phipphj(8,3,1)=\frac{1}{\sqrt{5}}
2288 : : !ang_phipphj(8,13,1)=-\sqrt{\frac{3}{35}}
2289 : : !ang_phipphj(8,15,1)=\frac{1}{\sqrt{7}}
2290 : : !ang_phipphj(9,4,1)=\frac{1}{\sqrt{5}}
2291 : : !ang_phipphj(9,14,1)=-\frac{1}{\sqrt{70}}
2292 : : !ang_phipphj(9,16,1)=\sqrt{\frac{3}{14}}
2293 : : !ang_phipphj(10,5,1)=\sqrt{\frac{3}{14}}
2294 : : !ang_phipphj(11,6,1)=\frac{1}{\sqrt{7}}
2295 : : !ang_phipphj(12,5,1)=-\frac{1}{\sqrt{70}}
2296 : : !ang_phipphj(13,8,1)=-\sqrt{\frac{3}{35}}
2297 : : !ang_phipphj(14,7,1)=\sqrt{\frac{6}{35}}
2298 : : !ang_phipphj(14,9,1)=-\frac{1}{\sqrt{70}}
2299 : : !ang_phipphj(15,8,1)=\frac{1}{\sqrt{7}}
2300 : : !ang_phipphj(16,9,1)=\sqrt{\frac{3}{14}}
2301 : : !ang_phipphj(1,4,2)=\frac{1}{2 \sqrt{3}}
2302 : : !ang_phipphj(1,14,2)=-\frac{\sqrt{\frac{7}{6}}}{2}
2303 : : !ang_phipphj(2,5,2)=\frac{1}{2 \sqrt{5}}
2304 : : !ang_phipphj(3,8,2)=\frac{1}{2 \sqrt{5}}
2305 : : !ang_phipphj(4,7,2)=-\sqrt{\frac{3}{5}}
2306 : : !ang_phipphj(4,9,2)=\frac{1}{2 \sqrt{5}}
2307 : : !ang_phipphj(5,2,2)=\frac{1}{4 \sqrt{5}}
2308 : : !ang_phipphj(5,10,2)=\frac{\sqrt{\frac{3}{14}}}{2}
2309 : : !ang_phipphj(5,12,2)=-2 \sqrt{\frac{2}{35}}
2310 : : !ang_phipphj(6,11,2)=\frac{1}{2 \sqrt{7}}
2311 : : !ang_phipphj(7,4,2)=\frac{1}{\sqrt{15}}
2312 : : !ang_phipphj(7,14,2)=\frac{13}{2 \sqrt{210}}
2313 : : !ang_phipphj(8,3,2)=-\frac{1}{\sqrt{5}}
2314 : : !ang_phipphj(8,13,2)=-4 \sqrt{\frac{3}{35}}
2315 : : !ang_phipphj(8,15,2)=\frac{1}{2 \sqrt{7}}
2316 : : !ang_phipphj(9,4,2)=\frac{1}{4 \sqrt{5}}
2317 : : !ang_phipphj(9,14,2)=-2 \sqrt{\frac{2}{35}}
2318 : : !ang_phipphj(9,16,2)=\frac{\sqrt{\frac{3}{14}}}{2}
2319 : : !ang_phipphj(10,5,2)=\frac{1}{\sqrt{42}}
2320 : : !ang_phipphj(11,6,2)=-\frac{1}{4 \sqrt{7}}
2321 : : !ang_phipphj(12,5,2)=\sqrt{\frac{2}{35}}
2322 : : !ang_phipphj(13,8,2)=2 \sqrt{\frac{3}{35}}
2323 : : !ang_phipphj(14,7,2)=-2 \sqrt{\frac{6}{35}}
2324 : : !ang_phipphj(14,9,2)=\sqrt{\frac{2}{35}}
2325 : : !ang_phipphj(15,8,2)=-\frac{1}{4 \sqrt{7}}
2326 : : !ang_phipphj(16,9,2)=\frac{1}{\sqrt{42}}
2327 : : !ang_phipphj(1,4,3)=\frac{\sqrt{3}}{2}
2328 : : !ang_phipphj(1,14,3)=\frac{\sqrt{\frac{7}{6}}}{2}
2329 : : !ang_phipphj(2,5,3)=\frac{\sqrt{5}}{2}
2330 : : !ang_phipphj(3,8,3)=\frac{\sqrt{5}}{2}
2331 : : !ang_phipphj(4,9,3)=\frac{\sqrt{5}}{2}
2332 : : !ang_phipphj(5,2,3)=-\frac{\sqrt{5}}{4}
2333 : : !ang_phipphj(5,10,3)=\frac{\sqrt{\frac{21}{2}}}{2}
2334 : : !ang_phipphj(6,11,3)=\frac{\sqrt{7}}{2}
2335 : : !ang_phipphj(7,14,3)=\frac{\sqrt{\frac{35}{6}}}{2}
2336 : : !ang_phipphj(8,15,3)=\frac{\sqrt{7}}{2}
2337 : : !ang_phipphj(9,4,3)=-\frac{\sqrt{5}}{4}
2338 : : !ang_phipphj(9,16,3)=\frac{\sqrt{\frac{21}{2}}}{2}
2339 : : !ang_phipphj(10,5,3)=-\sqrt{\frac{7}{6}}
2340 : : !ang_phipphj(11,6,3)=-\frac{\sqrt{7}}{4}
2341 : : !ang_phipphj(15,8,3)=-\frac{\sqrt{7}}{4}
2342 : : !ang_phipphj(16,9,3)=-\sqrt{\frac{7}{6}}
2343 : : !ang_phipphj(1,2,4)=\frac{1}{\sqrt{3}}
2344 : : !ang_phipphj(2,1,4)=\frac{1}{\sqrt{3}}
2345 : : !ang_phipphj(2,7,4)=-\frac{1}{\sqrt{15}}
2346 : : !ang_phipphj(2,9,4)=-\frac{1}{\sqrt{5}}
2347 : : !ang_phipphj(3,6,4)=\frac{1}{\sqrt{5}}
2348 : : !ang_phipphj(4,5,4)=\frac{1}{\sqrt{5}}
2349 : : !ang_phipphj(5,4,4)=\frac{1}{\sqrt{5}}
2350 : : !ang_phipphj(5,14,4)=-\frac{1}{\sqrt{70}}
2351 : : !ang_phipphj(5,16,4)=-\sqrt{\frac{3}{14}}
2352 : : !ang_phipphj(6,3,4)=\frac{1}{\sqrt{5}}
2353 : : !ang_phipphj(6,13,4)=-\sqrt{\frac{3}{35}}
2354 : : !ang_phipphj(6,15,4)=-\frac{1}{\sqrt{7}}
2355 : : !ang_phipphj(7,2,4)=-\frac{1}{\sqrt{15}}
2356 : : !ang_phipphj(7,12,4)=\sqrt{\frac{6}{35}}
2357 : : !ang_phipphj(8,11,4)=\frac{1}{\sqrt{7}}
2358 : : !ang_phipphj(9,2,4)=-\frac{1}{\sqrt{5}}
2359 : : !ang_phipphj(9,10,4)=\sqrt{\frac{3}{14}}
2360 : : !ang_phipphj(9,12,4)=\frac{1}{\sqrt{70}}
2361 : : !ang_phipphj(10,9,4)=\sqrt{\frac{3}{14}}
2362 : : !ang_phipphj(11,8,4)=\frac{1}{\sqrt{7}}
2363 : : !ang_phipphj(12,7,4)=\sqrt{\frac{6}{35}}
2364 : : !ang_phipphj(12,9,4)=\frac{1}{\sqrt{70}}
2365 : : !ang_phipphj(13,6,4)=-\sqrt{\frac{3}{35}}
2366 : : !ang_phipphj(14,5,4)=-\frac{1}{\sqrt{70}}
2367 : : !ang_phipphj(15,6,4)=-\frac{1}{\sqrt{7}}
2368 : : !ang_phipphj(16,5,4)=-\sqrt{\frac{3}{14}}
2369 : : !ang_phipphj(1,2,5)=\frac{1}{2 \sqrt{3}}
2370 : : !ang_phipphj(1,12,5)=-\frac{\sqrt{\frac{7}{6}}}{2}
2371 : : !ang_phipphj(2,7,5)=-\sqrt{\frac{3}{5}}
2372 : : !ang_phipphj(2,9,5)=-\frac{1}{2 \sqrt{5}}
2373 : : !ang_phipphj(3,6,5)=\frac{1}{2 \sqrt{5}}
2374 : : !ang_phipphj(4,5,5)=\frac{1}{2 \sqrt{5}}
2375 : : !ang_phipphj(5,4,5)=\frac{1}{4 \sqrt{5}}
2376 : : !ang_phipphj(5,14,5)=-2 \sqrt{\frac{2}{35}}
2377 : : !ang_phipphj(5,16,5)=-\frac{\sqrt{\frac{3}{14}}}{2}
2378 : : !ang_phipphj(6,3,5)=-\frac{1}{\sqrt{5}}
2379 : : !ang_phipphj(6,13,5)=-4 \sqrt{\frac{3}{35}}
2380 : : !ang_phipphj(6,15,5)=-\frac{1}{2 \sqrt{7}}
2381 : : !ang_phipphj(7,2,5)=\frac{1}{\sqrt{15}}
2382 : : !ang_phipphj(7,12,5)=\frac{13}{2 \sqrt{210}}
2383 : : !ang_phipphj(8,11,5)=\frac{1}{2 \sqrt{7}}
2384 : : !ang_phipphj(9,2,5)=-\frac{1}{4 \sqrt{5}}
2385 : : !ang_phipphj(9,10,5)=\frac{\sqrt{\frac{3}{14}}}{2}
2386 : : !ang_phipphj(9,12,5)=2 \sqrt{\frac{2}{35}}
2387 : : !ang_phipphj(10,9,5)=\frac{1}{\sqrt{42}}
2388 : : !ang_phipphj(11,8,5)=-\frac{1}{4 \sqrt{7}}
2389 : : !ang_phipphj(12,7,5)=-2 \sqrt{\frac{6}{35}}
2390 : : !ang_phipphj(12,9,5)=-\sqrt{\frac{2}{35}}
2391 : : !ang_phipphj(13,6,5)=2 \sqrt{\frac{3}{35}}
2392 : : !ang_phipphj(14,5,5)=\sqrt{\frac{2}{35}}
2393 : : !ang_phipphj(15,6,5)=\frac{1}{4 \sqrt{7}}
2394 : : !ang_phipphj(16,5,5)=-\frac{1}{\sqrt{42}}
2395 : : !ang_phipphj(1,2,6)=\frac{\sqrt{3}}{2}
2396 : : !ang_phipphj(1,12,6)=\frac{\sqrt{\frac{7}{6}}}{2}
2397 : : !ang_phipphj(2,9,6)=-\frac{\sqrt{5}}{2}
2398 : : !ang_phipphj(3,6,6)=\frac{\sqrt{5}}{2}
2399 : : !ang_phipphj(4,5,6)=\frac{\sqrt{5}}{2}
2400 : : !ang_phipphj(5,4,6)=-\frac{\sqrt{5}}{4}
2401 : : !ang_phipphj(5,16,6)=-\frac{\sqrt{\frac{21}{2}}}{2}
2402 : : !ang_phipphj(6,15,6)=-\frac{\sqrt{7}}{2}
2403 : : !ang_phipphj(7,12,6)=\frac{\sqrt{\frac{35}{6}}}{2}
2404 : : !ang_phipphj(8,11,6)=\frac{\sqrt{7}}{2}
2405 : : !ang_phipphj(9,2,6)=\frac{\sqrt{5}}{4}
2406 : : !ang_phipphj(9,10,6)=\frac{\sqrt{\frac{21}{2}}}{2}
2407 : : !ang_phipphj(10,9,6)=-\sqrt{\frac{7}{6}}
2408 : : !ang_phipphj(11,8,6)=-\frac{\sqrt{7}}{4}
2409 : : !ang_phipphj(15,6,6)=\frac{\sqrt{7}}{4}
2410 : : !ang_phipphj(16,5,6)=\sqrt{\frac{7}{6}}
2411 : : !ang_phipphj(1,3,7)=\frac{1}{\sqrt{3}}
2412 : : !ang_phipphj(2,6,7)=\frac{1}{\sqrt{5}}
2413 : : !ang_phipphj(3,1,7)=\frac{1}{\sqrt{3}}
2414 : : !ang_phipphj(3,7,7)=\frac{2}{\sqrt{15}}
2415 : : !ang_phipphj(4,8,7)=\frac{1}{\sqrt{5}}
2416 : : !ang_phipphj(5,11,7)=\frac{1}{\sqrt{7}}
2417 : : !ang_phipphj(6,2,7)=\frac{1}{\sqrt{5}}
2418 : : !ang_phipphj(6,12,7)=2 \sqrt{\frac{2}{35}}
2419 : : !ang_phipphj(7,3,7)=\frac{2}{\sqrt{15}}
2420 : : !ang_phipphj(7,13,7)=\frac{3}{\sqrt{35}}
2421 : : !ang_phipphj(8,4,7)=\frac{1}{\sqrt{5}}
2422 : : !ang_phipphj(8,14,7)=2 \sqrt{\frac{2}{35}}
2423 : : !ang_phipphj(9,15,7)=\frac{1}{\sqrt{7}}
2424 : : !ang_phipphj(11,5,7)=\frac{1}{\sqrt{7}}
2425 : : !ang_phipphj(12,6,7)=2 \sqrt{\frac{2}{35}}
2426 : : !ang_phipphj(13,7,7)=\frac{3}{\sqrt{35}}
2427 : : !ang_phipphj(14,8,7)=2 \sqrt{\frac{2}{35}}
2428 : : !ang_phipphj(15,9,7)=\frac{1}{\sqrt{7}}
2429 : : !ang_phipphj(1,3,8)=\frac{2}{\sqrt{3}}
2430 : : !ang_phipphj(2,6,8)=\frac{3}{\sqrt{5}}
2431 : : !ang_phipphj(3,7,8)=2 \sqrt{\frac{3}{5}}
2432 : : !ang_phipphj(4,8,8)=\frac{3}{\sqrt{5}}
2433 : : !ang_phipphj(5,11,8)=\frac{4}{\sqrt{7}}
2434 : : !ang_phipphj(6,2,8)=-\frac{1}{\sqrt{5}}
2435 : : !ang_phipphj(6,12,8)=8 \sqrt{\frac{2}{35}}
2436 : : !ang_phipphj(7,3,8)=-\frac{2}{\sqrt{15}}
2437 : : !ang_phipphj(7,13,8)=\frac{12}{\sqrt{35}}
2438 : : !ang_phipphj(8,4,8)=-\frac{1}{\sqrt{5}}
2439 : : !ang_phipphj(8,14,8)=8 \sqrt{\frac{2}{35}}
2440 : : !ang_phipphj(9,15,8)=\frac{4}{\sqrt{7}}
2441 : : !ang_phipphj(11,5,8)=-\frac{2}{\sqrt{7}}
2442 : : !ang_phipphj(12,6,8)=-4 \sqrt{\frac{2}{35}}
2443 : : !ang_phipphj(13,7,8)=-\frac{6}{\sqrt{35}}
2444 : : !ang_phipphj(14,8,8)=-4 \sqrt{\frac{2}{35}}
2445 : : !ang_phipphj(15,9,8)=-\frac{2}{\sqrt{7}}
2446 : :
2447 : :
2448 : 1158 : ang_phipphj_tmp=zero
2449 : : !
2450 : 1158 : ang_phipphj_tmp(1,4,1)=0.57735026918962576451_dp
2451 : 1158 : ang_phipphj_tmp(2,5,1)=0.44721359549995793928_dp
2452 : 1158 : ang_phipphj_tmp(3,8,1)=0.44721359549995793928_dp
2453 : 1158 : ang_phipphj_tmp(4,1,1)=0.57735026918962576451_dp
2454 : 1158 : ang_phipphj_tmp(4,7,1)=-0.25819888974716112568_dp
2455 : 1158 : ang_phipphj_tmp(4,9,1)=0.44721359549995793928_dp
2456 : 1158 : ang_phipphj_tmp(5,2,1)=0.44721359549995793928_dp
2457 : 1158 : ang_phipphj_tmp(5,10,1)=0.46291004988627573078_dp
2458 : 1158 : ang_phipphj_tmp(5,12,1)=-0.11952286093343936400_dp
2459 : 1158 : ang_phipphj_tmp(6,11,1)=0.37796447300922722721_dp
2460 : 1158 : ang_phipphj_tmp(7,4,1)=-0.25819888974716112568_dp
2461 : 1158 : ang_phipphj_tmp(7,14,1)=0.41403933560541253068_dp
2462 : 1158 : ang_phipphj_tmp(8,3,1)=0.44721359549995793928_dp
2463 : 1158 : ang_phipphj_tmp(8,13,1)=-0.29277002188455995381_dp
2464 : 1158 : ang_phipphj_tmp(8,15,1)=0.37796447300922722721_dp
2465 : 1158 : ang_phipphj_tmp(9,4,1)=0.44721359549995793928_dp
2466 : 1158 : ang_phipphj_tmp(9,14,1)=-0.11952286093343936400_dp
2467 : 1158 : ang_phipphj_tmp(9,16,1)=0.46291004988627573078_dp
2468 : 1158 : ang_phipphj_tmp(10,5,1)=0.46291004988627573078_dp
2469 : 1158 : ang_phipphj_tmp(11,6,1)=0.37796447300922722721_dp
2470 : 1158 : ang_phipphj_tmp(12,5,1)=-0.11952286093343936400_dp
2471 : 1158 : ang_phipphj_tmp(13,8,1)=-0.29277002188455995381_dp
2472 : 1158 : ang_phipphj_tmp(14,7,1)=0.41403933560541253068_dp
2473 : 1158 : ang_phipphj_tmp(14,9,1)=-0.11952286093343936400_dp
2474 : 1158 : ang_phipphj_tmp(15,8,1)=0.37796447300922722721_dp
2475 : 1158 : ang_phipphj_tmp(16,9,1)=0.46291004988627573078_dp
2476 : : !
2477 : 1158 : ang_phipphj_tmp(1,4,2)=0.28867513459481288225_dp
2478 : 1158 : ang_phipphj_tmp(1,14,2)=-0.54006172486732168591_dp
2479 : 1158 : ang_phipphj_tmp(2,5,2)=0.22360679774997896964_dp
2480 : 1158 : ang_phipphj_tmp(3,8,2)=0.22360679774997896964_dp
2481 : 1158 : ang_phipphj_tmp(4,7,2)=-0.77459666924148337704_dp
2482 : 1158 : ang_phipphj_tmp(4,9,2)=0.22360679774997896964_dp
2483 : 1158 : ang_phipphj_tmp(5,2,2)=0.11180339887498948482_dp
2484 : 1158 : ang_phipphj_tmp(5,10,2)=0.23145502494313786539_dp
2485 : 1158 : ang_phipphj_tmp(5,12,2)=-0.47809144373375745599_dp
2486 : 1158 : ang_phipphj_tmp(6,11,2)=0.18898223650461361361_dp
2487 : 1158 : ang_phipphj_tmp(7,4,2)=0.25819888974716112568_dp
2488 : 1158 : ang_phipphj_tmp(7,14,2)=0.44854261357253024157_dp
2489 : 1158 : ang_phipphj_tmp(8,3,2)=-0.44721359549995793928_dp
2490 : 1158 : ang_phipphj_tmp(8,13,2)=-1.1710800875382398152_dp
2491 : 1158 : ang_phipphj_tmp(8,15,2)=0.18898223650461361361_dp
2492 : 1158 : ang_phipphj_tmp(9,4,2)=0.11180339887498948482_dp
2493 : 1158 : ang_phipphj_tmp(9,14,2)=-0.47809144373375745599_dp
2494 : 1158 : ang_phipphj_tmp(9,16,2)=0.23145502494313786539_dp
2495 : 1158 : ang_phipphj_tmp(10,5,2)=0.15430334996209191026_dp
2496 : 1158 : ang_phipphj_tmp(11,6,2)=-0.094491118252306806804_dp
2497 : 1158 : ang_phipphj_tmp(12,5,2)=0.23904572186687872799_dp
2498 : 1158 : ang_phipphj_tmp(13,8,2)=0.58554004376911990761_dp
2499 : 1158 : ang_phipphj_tmp(14,7,2)=-0.82807867121082506136_dp
2500 : 1158 : ang_phipphj_tmp(14,9,2)=0.23904572186687872799_dp
2501 : 1158 : ang_phipphj_tmp(15,8,2)=-0.094491118252306806804_dp
2502 : 1158 : ang_phipphj_tmp(16,9,2)=0.15430334996209191026_dp
2503 : : !
2504 : 1158 : ang_phipphj_tmp(1,4,3)=0.86602540378443864676_dp
2505 : 1158 : ang_phipphj_tmp(1,14,3)=0.54006172486732168591_dp
2506 : 1158 : ang_phipphj_tmp(2,5,3)=1.1180339887498948482_dp
2507 : 1158 : ang_phipphj_tmp(3,8,3)=1.1180339887498948482_dp
2508 : 1158 : ang_phipphj_tmp(4,9,3)=1.1180339887498948482_dp
2509 : 1158 : ang_phipphj_tmp(5,2,3)=-0.55901699437494742410_dp
2510 : 1158 : ang_phipphj_tmp(5,10,3)=1.6201851746019650577_dp
2511 : 1158 : ang_phipphj_tmp(6,11,3)=1.3228756555322952953_dp
2512 : 1158 : ang_phipphj_tmp(7,14,3)=1.2076147288491198811_dp
2513 : 1158 : ang_phipphj_tmp(8,15,3)=1.3228756555322952953_dp
2514 : 1158 : ang_phipphj_tmp(9,4,3)=-0.55901699437494742410_dp
2515 : 1158 : ang_phipphj_tmp(9,16,3)=1.6201851746019650577_dp
2516 : 1158 : ang_phipphj_tmp(10,5,3)=-1.0801234497346433718_dp
2517 : 1158 : ang_phipphj_tmp(11,6,3)=-0.66143782776614764763_dp
2518 : 1158 : ang_phipphj_tmp(15,8,3)=-0.66143782776614764763_dp
2519 : 1158 : ang_phipphj_tmp(16,9,3)=-1.0801234497346433718_dp
2520 : : !
2521 : 1158 : ang_phipphj_tmp(1,2,4)=0.57735026918962576451_dp
2522 : 1158 : ang_phipphj_tmp(2,1,4)=0.57735026918962576451_dp
2523 : 1158 : ang_phipphj_tmp(2,7,4)=-0.25819888974716112568_dp
2524 : 1158 : ang_phipphj_tmp(2,9,4)=-0.44721359549995793928_dp
2525 : 1158 : ang_phipphj_tmp(3,6,4)=0.44721359549995793928_dp
2526 : 1158 : ang_phipphj_tmp(4,5,4)=0.44721359549995793928_dp
2527 : 1158 : ang_phipphj_tmp(5,4,4)=0.44721359549995793928_dp
2528 : 1158 : ang_phipphj_tmp(5,14,4)=-0.11952286093343936400_dp
2529 : 1158 : ang_phipphj_tmp(5,16,4)=-0.46291004988627573078_dp
2530 : 1158 : ang_phipphj_tmp(6,3,4)=0.44721359549995793928_dp
2531 : 1158 : ang_phipphj_tmp(6,13,4)=-0.29277002188455995381_dp
2532 : 1158 : ang_phipphj_tmp(6,15,4)=-0.37796447300922722721_dp
2533 : 1158 : ang_phipphj_tmp(7,2,4)=-0.25819888974716112568_dp
2534 : 1158 : ang_phipphj_tmp(7,12,4)=0.41403933560541253068_dp
2535 : 1158 : ang_phipphj_tmp(8,11,4)=0.37796447300922722721_dp
2536 : 1158 : ang_phipphj_tmp(9,2,4)=-0.44721359549995793928_dp
2537 : 1158 : ang_phipphj_tmp(9,10,4)=0.46291004988627573078_dp
2538 : 1158 : ang_phipphj_tmp(9,12,4)=0.11952286093343936400_dp
2539 : 1158 : ang_phipphj_tmp(10,9,4)=0.46291004988627573078_dp
2540 : 1158 : ang_phipphj_tmp(11,8,4)=0.37796447300922722721_dp
2541 : 1158 : ang_phipphj_tmp(12,7,4)=0.41403933560541253068_dp
2542 : 1158 : ang_phipphj_tmp(12,9,4)=0.11952286093343936400_dp
2543 : 1158 : ang_phipphj_tmp(13,6,4)=-0.29277002188455995381_dp
2544 : 1158 : ang_phipphj_tmp(14,5,4)=-0.11952286093343936400_dp
2545 : 1158 : ang_phipphj_tmp(15,6,4)=-0.37796447300922722721_dp
2546 : 1158 : ang_phipphj_tmp(16,5,4)=-0.46291004988627573078_dp
2547 : : !
2548 : 1158 : ang_phipphj_tmp(1,2,5)=0.28867513459481288225_dp
2549 : 1158 : ang_phipphj_tmp(1,12,5)=-0.54006172486732168591_dp
2550 : 1158 : ang_phipphj_tmp(2,7,5)=-0.77459666924148337704_dp
2551 : 1158 : ang_phipphj_tmp(2,9,5)=-0.22360679774997896964_dp
2552 : 1158 : ang_phipphj_tmp(3,6,5)=0.22360679774997896964_dp
2553 : 1158 : ang_phipphj_tmp(4,5,5)=0.22360679774997896964_dp
2554 : 1158 : ang_phipphj_tmp(5,4,5)=0.11180339887498948482_dp
2555 : 1158 : ang_phipphj_tmp(5,14,5)=-0.47809144373375745599_dp
2556 : 1158 : ang_phipphj_tmp(5,16,5)=-0.23145502494313786539_dp
2557 : 1158 : ang_phipphj_tmp(6,3,5)=-0.44721359549995793928_dp
2558 : 1158 : ang_phipphj_tmp(6,13,5)=-1.1710800875382398152_dp
2559 : 1158 : ang_phipphj_tmp(6,15,5)=-0.18898223650461361361_dp
2560 : 1158 : ang_phipphj_tmp(7,2,5)=0.25819888974716112568_dp
2561 : 1158 : ang_phipphj_tmp(7,12,5)=0.44854261357253024157_dp
2562 : 1158 : ang_phipphj_tmp(8,11,5)=0.18898223650461361361_dp
2563 : 1158 : ang_phipphj_tmp(9,2,5)=-0.11180339887498948482_dp
2564 : 1158 : ang_phipphj_tmp(9,10,5)=0.23145502494313786539_dp
2565 : 1158 : ang_phipphj_tmp(9,12,5)=0.47809144373375745599_dp
2566 : 1158 : ang_phipphj_tmp(10,9,5)=0.15430334996209191026_dp
2567 : 1158 : ang_phipphj_tmp(11,8,5)=-0.094491118252306806804_dp
2568 : 1158 : ang_phipphj_tmp(12,7,5)=-0.82807867121082506136_dp
2569 : 1158 : ang_phipphj_tmp(12,9,5)=-0.23904572186687872799_dp
2570 : 1158 : ang_phipphj_tmp(13,6,5)=0.58554004376911990761_dp
2571 : 1158 : ang_phipphj_tmp(14,5,5)=0.23904572186687872799_dp
2572 : 1158 : ang_phipphj_tmp(15,6,5)=0.094491118252306806804_dp
2573 : 1158 : ang_phipphj_tmp(16,5,5)=-0.15430334996209191026_dp
2574 : : !
2575 : 1158 : ang_phipphj_tmp(1,2,6)=0.86602540378443864676_dp
2576 : 1158 : ang_phipphj_tmp(1,12,6)=0.54006172486732168591_dp
2577 : 1158 : ang_phipphj_tmp(2,9,6)=-1.1180339887498948482_dp
2578 : 1158 : ang_phipphj_tmp(3,6,6)=1.1180339887498948482_dp
2579 : 1158 : ang_phipphj_tmp(4,5,6)=1.1180339887498948482_dp
2580 : 1158 : ang_phipphj_tmp(5,4,6)=-0.55901699437494742410_dp
2581 : 1158 : ang_phipphj_tmp(5,16,6)=-1.6201851746019650577_dp
2582 : 1158 : ang_phipphj_tmp(6,15,6)=-1.3228756555322952953_dp
2583 : 1158 : ang_phipphj_tmp(7,12,6)=1.2076147288491198811_dp
2584 : 1158 : ang_phipphj_tmp(8,11,6)=1.3228756555322952953_dp
2585 : 1158 : ang_phipphj_tmp(9,2,6)=0.55901699437494742410_dp
2586 : 1158 : ang_phipphj_tmp(9,10,6)=1.6201851746019650577_dp
2587 : 1158 : ang_phipphj_tmp(10,9,6)=-1.0801234497346433718_dp
2588 : 1158 : ang_phipphj_tmp(11,8,6)=-0.66143782776614764763_dp
2589 : 1158 : ang_phipphj_tmp(15,6,6)=0.66143782776614764763_dp
2590 : 1158 : ang_phipphj_tmp(16,5,6)=1.0801234497346433718_dp
2591 : : !
2592 : 1158 : ang_phipphj_tmp(1,3,7)=0.57735026918962576451_dp
2593 : 1158 : ang_phipphj_tmp(2,6,7)=0.44721359549995793928_dp
2594 : 1158 : ang_phipphj_tmp(3,1,7)=0.57735026918962576451_dp
2595 : 1158 : ang_phipphj_tmp(3,7,7)=0.51639777949432225136_dp
2596 : 1158 : ang_phipphj_tmp(4,8,7)=0.44721359549995793928_dp
2597 : 1158 : ang_phipphj_tmp(5,11,7)=0.37796447300922722721_dp
2598 : 1158 : ang_phipphj_tmp(6,2,7)=0.44721359549995793928_dp
2599 : 1158 : ang_phipphj_tmp(6,12,7)=0.47809144373375745599_dp
2600 : 1158 : ang_phipphj_tmp(7,3,7)=0.51639777949432225136_dp
2601 : 1158 : ang_phipphj_tmp(7,13,7)=0.50709255283710994651_dp
2602 : 1158 : ang_phipphj_tmp(8,4,7)=0.44721359549995793928_dp
2603 : 1158 : ang_phipphj_tmp(8,14,7)=0.47809144373375745599_dp
2604 : 1158 : ang_phipphj_tmp(9,15,7)=0.37796447300922722721_dp
2605 : 1158 : ang_phipphj_tmp(11,5,7)=0.37796447300922722721_dp
2606 : 1158 : ang_phipphj_tmp(12,6,7)=0.47809144373375745599_dp
2607 : 1158 : ang_phipphj_tmp(13,7,7)=0.50709255283710994651_dp
2608 : 1158 : ang_phipphj_tmp(14,8,7)=0.47809144373375745599_dp
2609 : 1158 : ang_phipphj_tmp(15,9,7)=0.37796447300922722721_dp
2610 : : !
2611 : 1158 : ang_phipphj_tmp(1,3,8)=1.1547005383792515290_dp
2612 : 1158 : ang_phipphj_tmp(2,6,8)=1.3416407864998738178_dp
2613 : 1158 : ang_phipphj_tmp(3,7,8)=1.5491933384829667541_dp
2614 : 1158 : ang_phipphj_tmp(4,8,8)=1.3416407864998738178_dp
2615 : 1158 : ang_phipphj_tmp(5,11,8)=1.5118578920369089089_dp
2616 : 1158 : ang_phipphj_tmp(6,2,8)=-0.44721359549995793928_dp
2617 : 1158 : ang_phipphj_tmp(6,12,8)=1.9123657749350298240_dp
2618 : 1158 : ang_phipphj_tmp(7,3,8)=-0.51639777949432225136_dp
2619 : 1158 : ang_phipphj_tmp(7,13,8)=2.0283702113484397860_dp
2620 : 1158 : ang_phipphj_tmp(8,4,8)=-0.44721359549995793928_dp
2621 : 1158 : ang_phipphj_tmp(8,14,8)=1.9123657749350298240_dp
2622 : 1158 : ang_phipphj_tmp(9,15,8)=1.5118578920369089089_dp
2623 : 1158 : ang_phipphj_tmp(11,5,8)=-0.75592894601845445443_dp
2624 : 1158 : ang_phipphj_tmp(12,6,8)=-0.95618288746751491198_dp
2625 : 1158 : ang_phipphj_tmp(13,7,8)=-1.0141851056742198930_dp
2626 : 1158 : ang_phipphj_tmp(14,8,8)=-0.95618288746751491198_dp
2627 : 1158 : ang_phipphj_tmp(15,9,8)=-0.75592894601845445443_dp
2628 : :
2629 [ + + + + : 198502 : ang_phipphj(:,:,:)=ang_phipphj_tmp(1:mpsang**2,1:mpsang**2,:)
+ + ]
2630 : :
2631 : 1158 : end subroutine setnabla_ylm
2632 : : !!***
2633 : :
2634 : : !----------------------------------------------------------------------
2635 : :
2636 : : !!****f* m_paw_sphharm/gaunt
2637 : : !! NAME
2638 : : !! gaunt
2639 : : !!
2640 : : !! FUNCTION
2641 : : !! Returns gaunt coefficient, i.e.
2642 : : !! the integral of Sqrt[4 \pi] Y*(l_i,m_i) Y*(ll,mm) Y(l_j,m_j)
2643 : : !! See the 3-j and 6-j symbols by Rotenberg, etc., (Technology Press, 1959), pg.5.
2644 : : !!
2645 : : !! INPUTS
2646 : : !! ll,mm,l1,l2,m1,m2= six quantum numbers defining the Gaunt coef.
2647 : : !!
2648 : : !! OUTPUT
2649 : : !! gaunt(ll,mm,l1,l2,m1,m2)=the value of the integral
2650 : : !!
2651 : : !! SOURCE
2652 : :
2653 : 315928 : function gaunt(ll,mm,l1,m1,l2,m2)
2654 : :
2655 : : !Arguments ---------------------------------------------
2656 : : !scalars
2657 : : integer,intent(in) :: l1,l2,ll,m1,m2,mm
2658 : : real(dp) :: gaunt
2659 : :
2660 : : !Local variables ------------------------------
2661 : : !scalars
2662 : : integer :: i1,i2,j1,j1half,j2,j2half,j3,j3half,j_half,jj,k1,k2,n1,n2
2663 : : real(dp) :: argument,sign,sum,xx,yy
2664 : : logical :: ok
2665 : :
2666 : : !************************************************************************
2667 : :
2668 : 315928 : gaunt=zero;sum=zero;ok =.true.
2669 : :
2670 [ - + ]: 315928 : if((-m1-mm+m2) /= 0) ok = .false.
2671 [ - + ]: 315928 : if(abs(m1) > l1) ok = .false.
2672 [ - + ]: 315928 : if(abs(mm) > ll) ok = .false.
2673 [ - + ]: 315928 : if(abs(m2) > l2) ok = .false.
2674 : :
2675 : 315928 : jj = l1 + ll + l2
2676 [ - + ]: 315928 : if (mod(jj,2)/=0) ok = .false.
2677 : 315928 : j1 = jj-2*l2
2678 : 315928 : j2 = jj-2*ll
2679 : 315928 : j3 = jj-2*l1
2680 : :
2681 [ + - + - ]: 315928 : if (j1<0 .or. j2<0 .or. j3<0) ok = .false.
2682 : :
2683 [ + - ]: 315928 : if (ok) then
2684 : :
2685 : 315928 : xx = (2 * l1 + 1) * (2 * ll + 1) * (2 * l2 + 1)
2686 : :
2687 : 315928 : j1half = j1/2
2688 : 315928 : j2half = j2/2
2689 : 315928 : j3half = j3/2
2690 : 315928 : j_half = jj/2
2691 : :
2692 : 315928 : gaunt = (-1)**j1half * sqrt(xx)
2693 : 1263712 : gaunt = gaunt * rfactorial(j2)*rfactorial(j3)/rfactorial(jj+1)
2694 : : gaunt = gaunt * rfactorial(j_half)/(rfactorial(j1half)&
2695 : 1307312 : & * rfactorial(j2half)*rfactorial(j3half))
2696 : :
2697 : 1459041 : yy = rfactorial(l2 + m2) * rfactorial(l2 - m2)
2698 : :
2699 [ + + ]: 315928 : if (mm>=0) then
2700 : 194884 : yy = yy * perms(ll+mm,2*mm)
2701 : : else
2702 : 121044 : yy = yy / perms(ll-mm,-2*mm)
2703 : : end if
2704 : :
2705 [ + + ]: 315928 : if (m1>=0) then
2706 : 190982 : yy = yy / perms(l1+m1,2*m1)
2707 : : else
2708 : 124946 : yy = yy * perms(l1-m1,-2*m1)
2709 : : end if
2710 : :
2711 : 315928 : gaunt = gaunt * sqrt(yy)
2712 : :
2713 : 315928 : i1 = l2 - ll - m1
2714 : 315928 : i2 = l2 - l1 + mm
2715 : 315928 : k1 = -min(0, i1, i2)
2716 : 315928 : n1 = l1 + m1
2717 : 315928 : n2 = ll - mm
2718 : 315928 : k2 = min(j1, n1, n2)
2719 : :
2720 : 315928 : sign = 1._dp
2721 [ + + ]: 315928 : if(k1>0) sign = (-1._dp)**k1
2722 : :
2723 : 631856 : argument = sign * perms(n1,k1)/rfactorial(k1)
2724 : 631856 : argument = argument * perms(n2,k1)/rfactorial(i1 + k1)
2725 : 631856 : argument = argument * perms(j1,k1)/rfactorial(i2 + k1)
2726 : 315928 : sum = sum + argument
2727 : :
2728 : 315928 : sign = -sign
2729 : 315928 : k1 = k1 + 1
2730 [ + + ]: 359290 : do while(k1 <= k2)
2731 : 86724 : argument = sign * perms(n1, k1)/rfactorial(k1)
2732 : 86724 : argument = argument * perms(n2, k1)/rfactorial(i1 + k1)
2733 : 86724 : argument = argument * perms(j1, k1)/rfactorial(i2 + k1)
2734 : 43362 : sum = sum + argument
2735 : 43362 : sign = -sign
2736 : 43362 : k1 = k1 + 1
2737 : : end do
2738 : :
2739 : : end if
2740 : :
2741 : 315928 : gaunt = gaunt * sum
2742 : :
2743 : 315928 : end function gaunt
2744 : : !!***
2745 : :
2746 : : !----------------------------------------------------------------------
2747 : :
2748 : : !!****f* m_pawsphharm/realgaunt
2749 : : !! NAME
2750 : : !! realgaunt
2751 : : !!
2752 : : !! FUNCTION
2753 : : !! This routine compute "real Gaunt coefficients", i.e. gaunt
2754 : : !! coefficients according to "real spherical harmonics"
2755 : : !! RealGaunt(ilm,ilm_i,ilm_j) = Int[ S_lm Slm_i Slm_j]
2756 : : !!
2757 : : !! INPUTS
2758 : : !! l_max= max. value of ang. momentum l+1; Gaunt coeffs up to
2759 : : !! [(2*l_max-1,m),(l_max,m),(l_max,m)] are computed
2760 : : !!
2761 : : !! OUTPUT
2762 : : !! gntselect((2*l_max-1)**2,l_max**2*(l_max**2+1)/2)=
2763 : : !! selection rules for Gaunt coefficients
2764 : : !! if Gaunt coeff. is zero, gntselect=0
2765 : : !! if Gaunt coeff. is non-zero, gntselect is the index of
2766 : : !! the coeff. in realgnt(:) array
2767 : : !! ngnt= number of non-zero Gaunt coefficients
2768 : : !! realgnt((2*l_max-1)**2*l_max**4)= non-zero real Gaunt coefficients
2769 : : !!
2770 : : !! NOTE
2771 : : !! Second index of gntselect is in "upper triangle" format.
2772 : : !! Its formula is klm_ij = ilm_i*(ilm_i-1)/2 + ilm_j,
2773 : : !! corresponding to the two index pairs: (ilm_i,ilm_j) and (ilm_j,ilmj)
2774 : : !!
2775 : : !! SOURCE
2776 : :
2777 [ - + - + ]: 2318 : subroutine realgaunt(l_max,ngnt,gntselect,realgnt)
2778 : :
2779 : : !Arguments ---------------------------------------------
2780 : : !scalars
2781 : : integer,intent(in) :: l_max
2782 : : integer,intent(out) :: ngnt
2783 : : !arrays
2784 : : integer,intent(out) :: gntselect(:,:)
2785 : : real(dp),intent(out) :: realgnt(:)
2786 : :
2787 : : !Local variables ------------------------------
2788 : : !scalars
2789 : : integer :: ilm1,ilm2,ilmp1,k0lm1,klm1,l1,l2,ll,lp1,m1,m2,mm,mm1,mm2,mm3,mp1
2790 : : real(dp) :: c11,c12,c21,c22,c31,c32,fact,realgnt_tmp
2791 : : character(len=500) :: msg
2792 : : !arrays
2793 : 2318 : integer,allocatable :: ssgn(:)
2794 : 2318 : type(coeff3_type), allocatable :: coeff(:)
2795 : :
2796 : : !************************************************************************
2797 : :
2798 [ + + + - : 6954 : if ( size(gntselect)<(2*l_max-1)**2*(l_max**2*(l_max**2+1))/2 .or. &
- + ]
2799 : : & size(realgnt) <(2*l_max-1)**2*(l_max**2*(l_max**2+1))/2 ) then
2800 : 0 : msg='Too small sizes for gntselect/realgnt!'
2801 : 0 : LIBPAW_BUG(msg)
2802 : : end if
2803 : :
2804 : : !Initialize output arrays with zeros.
2805 [ + + + + : 6519833 : gntselect = 0; realgnt = zero
+ + ]
2806 : :
2807 : : !Compute matrix cc where Sl=cc*Yl (Sl=real sph. harm.)
2808 : : !------------------------------------------------
2809 [ + - - + : 21748 : LIBPAW_DATATYPE_ALLOCATE(coeff,(4*l_max-3))
+ - + - +
- - + - +
+ + ]
2810 [ + + ]: 17112 : do ll=1,4*l_max-3
2811 [ - + + - : 59176 : LIBPAW_ALLOCATE(coeff(ll)%value,(2,2*ll-1,2*ll-1))
- + + - +
- - + - +
- + ]
2812 [ + + + + : 4235198 : coeff(ll)%value(:,:,:)=zero
+ + ]
2813 : 14794 : coeff(ll)%value(1,ll,ll)=one
2814 [ + + ]: 65838 : do mm=1,ll-1
2815 : 48726 : coeff(ll)%value(1,ll+mm,ll+mm)= (-1._dp)**mm/sqrt(2._dp)
2816 : 48726 : coeff(ll)%value(1,ll-mm,ll+mm)= ( 1._dp) /sqrt(2._dp)
2817 : 48726 : coeff(ll)%value(2,ll+mm,ll-mm)=-(-1._dp)**mm/sqrt(2._dp)
2818 : 63520 : coeff(ll)%value(2,ll-mm,ll-mm)= ( 1._dp) /sqrt(2._dp)
2819 : : end do
2820 : : end do
2821 : :
2822 [ - + + - : 6954 : LIBPAW_ALLOCATE(ssgn,(l_max**2))
+ - + - -
+ - + ]
2823 [ + + ]: 16185 : ssgn(:)=1
2824 [ + - ]: 2318 : if (l_max>0) then
2825 [ + + ]: 5437 : do l1=1,l_max-1
2826 : 3119 : ilm1=1+l1**2+l1
2827 [ + + ]: 9652 : do m1=-l1,-1
2828 : 7334 : ssgn(ilm1+m1)=-1
2829 : : end do
2830 : : end do
2831 : : end if
2832 : :
2833 : 2318 : ngnt=0
2834 : :
2835 : : !Loop on (lp1,mp1)
2836 : : !------------------------------------------------
2837 [ + + ]: 7755 : do lp1=0,l_max-1
2838 [ + + ]: 21622 : do mp1=-lp1,lp1
2839 : 13867 : ilmp1=1+lp1**2+lp1+mp1
2840 : 13867 : k0lm1=ilmp1*(ilmp1-1)/2
2841 : :
2842 : : ! Loop on (l1,m1)<=(lp1,mp1)
2843 : : ! ------------------------------------------------
2844 [ + + ]: 58413 : do l1=0,l_max-1
2845 [ + + ]: 175227 : do m1=-l1,l1
2846 : 122251 : ilm1=1+l1**2+l1+m1
2847 : :
2848 [ + + ]: 161360 : if (ilm1<=ilmp1) then
2849 : :
2850 : 68059 : klm1=k0lm1+ilm1
2851 [ + + ]: 2528462 : gntselect(:,klm1)=0
2852 : :
2853 : : ! Loop on (l2,m2)
2854 : : ! ------------------------------------------------
2855 [ - + + + ]: 226047 : do l2=abs(l1-lp1),l1+lp1,2
2856 [ + + ]: 1110735 : do m2=-l2,l2
2857 : 884688 : ilm2=1+l2**2+l2+m2
2858 : :
2859 : : ! Real Gaunt coeffs selection rules
2860 : : ! ------------------------------------------------
2861 [ + + + + : 884688 : if ((l2<=l1+lp1).and.&
+ + + + +
+ + + + +
+ + ]
2862 : : & (((m1== mp1).and.((m2==0).or.(m2==2*abs(mp1)))).or.&
2863 : : & ((m1==-mp1).and.(m2==-abs(m1)-abs(mp1))).or.&
2864 : : & ((abs(m1)/=(abs(mp1)).and.&
2865 : : & ((m2==ssgn(ilm1)*ssgn(ilmp1)* (abs(m1)+abs(mp1))).or.&
2866 : : & (m2==ssgn(ilm1)*ssgn(ilmp1)*abs(abs(m1)-abs(mp1)))&
2867 [ + - ]: 1042676 : ))))) then
2868 : :
2869 : : ! Compute selected real Gaunt coefficient
2870 : : ! ------------------------------------------------
2871 : 165596 : realgnt_tmp=zero
2872 [ + + ]: 905536 : do mm1=-l1,l1
2873 : 739940 : c11=coeff(l1+1)%value(1,l1+mm1+1,l1+m1+1)
2874 : 739940 : c12=coeff(l1+1)%value(2,l1+mm1+1,l1+m1+1)
2875 [ + + ]: 5527020 : do mm2= -lp1,lp1
2876 : 4621484 : c21=coeff(lp1+1)%value(1,lp1+mm2+1,lp1+mp1+1)
2877 : 4621484 : c22=coeff(lp1+1)%value(2,lp1+mm2+1,lp1+mp1+1)
2878 [ + + ]: 43348876 : do mm3= -l2,l2
2879 : 37987452 : c31=coeff(l2+1)%value(1,l2+mm3+1,l2+m2+1)
2880 : 37987452 : c32=coeff(l2+1)%value(2,l2+mm3+1,l2+m2+1)
2881 : : fact=c11*c21*c31 - c12*c22*c31&
2882 : 37987452 : & -c11*c22*c32 - c12*c21*c32
2883 [ + + + + ]: 37987452 : if((abs(fact)>=tol12).and.(mm3==-mm2-mm1)) &
2884 : : & realgnt_tmp=realgnt_tmp+fact*(-1)**mm2 &
2885 : 4937412 : & *gaunt(l2,mm3,l1,mm1,lp1,-mm2)
2886 : : end do
2887 : : end do
2888 : : end do
2889 : :
2890 : : ! Count additional non-zero real Gaunt coeffs
2891 : : ! ------------------------------------------------
2892 [ + + ]: 165596 : if (abs(realgnt_tmp)>=tol12) then
2893 : 164512 : ngnt=ngnt+1
2894 : 164512 : gntselect(ilm2,klm1)=ngnt
2895 : 164512 : realgnt(ngnt)=realgnt_tmp/sqrt(four_pi)
2896 : : end if
2897 : :
2898 : : ! End loops
2899 : : ! ------------------------------------------------
2900 : : end if
2901 : : end do
2902 : : end do
2903 : : end if
2904 : : end do
2905 : : end do
2906 : : end do
2907 : : end do
2908 : :
2909 : : !Deallocate memory
2910 : : !------------------------------------------------
2911 [ + + ]: 17112 : do ll=1,4*l_max-3
2912 [ - + ]: 17112 : LIBPAW_DEALLOCATE(coeff(ll)%value)
2913 : : end do
2914 [ + - + + : 17112 : LIBPAW_DATATYPE_DEALLOCATE(coeff)
- + - + ]
2915 [ - + ]: 2318 : LIBPAW_DEALLOCATE(ssgn)
2916 : :
2917 : 2318 : end subroutine realgaunt
2918 : : !!***
2919 : :
2920 : : !----------------------------------------------------------------------
2921 : :
2922 : : !!****f* m_pawsphharm/make_dyadic
2923 : : !! NAME
2924 : : !! make_dyadic
2925 : : !!
2926 : : !! FUNCTION
2927 : : !! compute integrals S_limi (\alpha*1-\beta\hat{r}\hat{r}) S_ljmj
2928 : : !! for the general dyadic \alpha*1-\beta\hat{r}\hat{r}
2929 : : !!
2930 : : !! INPUTS
2931 : : !!
2932 : : !! OUTPUT
2933 : : !! dyadic(3,3,lm2_size)=dyadic term for each S_iS_j pair
2934 : : !!
2935 : : !! NOTE
2936 : : !! dyadic is indexed in voigt-style format:
2937 : : !! (1,1,klm) : xx
2938 : : !! (2,2,klm) : xx
2939 : : !! (3,3,klm) : xx
2940 : : !! (2,3,klm) : yz
2941 : : !! (1,3,klm) : xz
2942 : : !! (1,2,klm) : xy
2943 : : !! and symmetric around diagonal
2944 : : !!
2945 : : !! SOURCE
2946 : :
2947 : 620 : subroutine make_dyadic(alpha,beta,dyadic,gntselect,gs1,gs2,lm2_size,ngnt,realgnt)
2948 : :
2949 : : !Arguments ---------------------------------------------
2950 : : !scalars
2951 : : integer,intent(in) :: gs1,gs2,lm2_size,ngnt
2952 : : real(dp),intent(in) :: alpha,beta
2953 : : !arrays
2954 : : integer,intent(in) :: gntselect(gs1,gs2)
2955 : : real(dp),intent(in) :: realgnt(ngnt)
2956 : : real(dp),intent(out) :: dyadic(3,3,lm2_size)
2957 : :
2958 : : !Local variables ------------------------------
2959 : : !scalars
2960 : : integer :: dlm,ignt,klm
2961 : : real(dp) :: afact,rgnt
2962 : : real(dp),parameter :: c1=sqrt(four_pi/five)
2963 : : real(dp),parameter :: c2=one/sqrt(three)
2964 : : real(dp),parameter :: c3=sqrt(four_pi)
2965 : : !arrays
2966 : :
2967 : : !************************************************************************
2968 : :
2969 [ + + + + : 90320 : dyadic = zero
+ + ]
2970 [ + + ]: 7520 : do klm = 1, lm2_size
2971 [ + + ]: 69620 : do dlm = 1, 9
2972 [ + + ]: 62100 : if ( (dlm>1) .AND. (dlm<5) ) cycle ! no L=1 contributions
2973 : 41400 : ignt = gntselect(dlm,klm)
2974 [ + + ]: 41400 : if (ignt<1) cycle
2975 : 7960 : rgnt=realgnt(ignt)
2976 : 2580 : select case(dlm)
2977 : :
2978 : : ! S00, appears in xx, yy, zz
2979 : : case (1)
2980 : 2580 : afact = rgnt*c3*(alpha-beta/three)
2981 : 2580 : dyadic(1,1,klm) = dyadic(1,1,klm)+afact
2982 : 2580 : dyadic(2,2,klm) = dyadic(2,2,klm)+afact
2983 : 2580 : dyadic(3,3,klm) = dyadic(3,3,klm)+afact
2984 : :
2985 : : ! S2,-2 appears in xy
2986 : : case (5)
2987 : 680 : dyadic(1,2,klm) = dyadic(1,2,klm)-rgnt*beta*c1*c2
2988 : :
2989 : : ! S2,-1 appears yz
2990 : : case (6)
2991 : 700 : dyadic(2,3,klm) = dyadic(2,3,klm)-rgnt*beta*c1*c2
2992 : :
2993 : : ! S2,0 appears in xx, yy, zz
2994 : : case (7)
2995 : 1980 : afact = rgnt*beta*c1/three
2996 : 1980 : dyadic(1,1,klm) = dyadic(1,1,klm) + afact
2997 : 1980 : dyadic(2,2,klm) = dyadic(2,2,klm) + afact
2998 : 1980 : dyadic(3,3,klm) = dyadic(3,3,klm) - two*afact
2999 : :
3000 : : ! S2,1 appears in xz
3001 : : case (8)
3002 : 700 : dyadic(1,3,klm) = dyadic(1,3,klm)-rgnt*beta*c1*c2
3003 : :
3004 : : ! S2,2 appears in xx, yy
3005 : : case (9)
3006 : 1320 : afact=rgnt*beta*c1*c2
3007 : 1320 : dyadic(1,1,klm) = dyadic(1,1,klm)-afact
3008 [ + + + + : 7960 : dyadic(2,2,klm) = dyadic(2,2,klm)+afact
+ + - ]
3009 : : case default
3010 : : end select
3011 : :
3012 : 7960 : dyadic(3,1,klm) = dyadic(1,3,klm)
3013 : 7960 : dyadic(2,1,klm) = dyadic(1,2,klm)
3014 : 69000 : dyadic(3,2,klm) = dyadic(2,3,klm)
3015 : :
3016 : : end do
3017 : : end do
3018 : :
3019 : 620 : end subroutine make_dyadic
3020 : : !!***
3021 : :
3022 : : !----------------------------------------------------------------------
3023 : :
3024 : :
3025 : : !!****f* m_paw_sphharm/nablarealgaunt
3026 : : !! NAME
3027 : : !! nablarealgaunt
3028 : : !!
3029 : : !! FUNCTION
3030 : : !! Evaluate integrals involving spherical harmonics and their gradient.
3031 : : !! These integrals are angular part for <nablaphi|nablaphj> and <tnablaphi|tnablaphj>
3032 : : !! Nabla_RealGaunt(ilm,ilm_i,ilm_j) = Int[ Slm Grad(Slm_i).Grad(Slm_j)]
3033 : : !!
3034 : : !! INPUTS
3035 : : !! l_max = 1 + max. l value for Slm (see description above)
3036 : : !! l_max_ij = 1 + max. l value for Slm_i and Slm_j (see description above)
3037 : : !!
3038 : : !! OUTPUT
3039 : : !! nnablagnt= number of non-zero integrals
3040 : : !! nabgauntselect(l_max**2,l_max_ij**2,l_max_ij**2)= indexes of the non-zero integrals
3041 : : !! nablagaunt(l_max**2*l_max_ij**4)= values of the integrals
3042 : : !!
3043 : : !! SOURCE
3044 : :
3045 [ - + - + ]: 16 : subroutine nablarealgaunt(l_max,l_max_ij,nnablagnt,nabgauntselect,nablagaunt)
3046 : :
3047 : : !Arguments ---------------------------------------------
3048 : : !scalars
3049 : : integer, intent(in) :: l_max,l_max_ij
3050 : : integer, intent(out) :: nnablagnt
3051 : : !array
3052 : : integer,intent(out) :: nabgauntselect(:,:,:)
3053 : : real(dp),intent(out) :: nablagaunt(:)
3054 : :
3055 : : !Local variables ---------------------------------------
3056 : : logical,parameter :: debug=.false.
3057 : : integer :: angl_size,ii,ilm,ilm_i,ilm_j,ipt,mpsang,ntheta,nphi,ylm_size
3058 : : real(dp) :: nabla_rg, yylmgr
3059 : : character(len=500) :: msg
3060 : 16 : real(dp),allocatable :: ang_wgth(:),cart_coord(:,:),ylmr(:,:),ylmrgr(:,:,:)
3061 : :
3062 : : !************************************************************************
3063 : :
3064 [ + + + - : 64 : if ( size(nabgauntselect)< (l_max**2)*(l_max_ij**4) .or. &
- + ]
3065 : : & size(nablagaunt) < (l_max**2)*(l_max_ij**4) ) then
3066 : 0 : msg='Too small sizes for nabgauntselect/nablagaunt!'
3067 : 0 : LIBPAW_BUG(msg)
3068 : : end if
3069 : :
3070 [ + + + + : 22150 : nabgauntselect(:,:,:)=-1
+ + ]
3071 [ + + ]: 21130 : nablagaunt(:)=zero
3072 : :
3073 : 16 : ii=0
3074 [ + - ]: 16 : if (l_max>1) then
3075 [ + - ]: 16 : if (l_max_ij>=1) then
3076 : 16 : ii=ii+1 ; nabgauntselect(1,2,2)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi))
3077 : 16 : ii=ii+1 ; nabgauntselect(1,3,3)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi))
3078 : 16 : ii=ii+1 ; nabgauntselect(1,4,4)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi))
3079 : : end if
3080 [ + + ]: 16 : if (l_max_ij>2) then
3081 : 10 : ii=ii+1 ; nabgauntselect(1,5,5)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
3082 : 10 : ii=ii+1 ; nabgauntselect(1,6,6)=ii ; nablagaunt(ii)=1.692568750643269_dp !\dfrac{3}{\sqrt{\pi}}
3083 : 10 : ii=ii+1 ; nabgauntselect(1,7,7)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
3084 : 10 : ii=ii+1 ; nabgauntselect(1,8,8)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
3085 : 10 : ii=ii+1 ; nabgauntselect(1,9,9)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
3086 : 10 : ii=ii+1 ; nabgauntselect(2,2,7)=ii ; nablagaunt(ii)=-0.37846987830302403_dp !\frac{-3}{2\sqrt{5\pi}}
3087 : 10 : ii=ii+1 ; nabgauntselect(2,2,9)=ii ; nablagaunt(ii)=-0.6555290583552474_dp !\frac{-1.5\sqrt{3}}{\sqrt{5\pi}}
3088 : 10 : ii=ii+1 ; nabgauntselect(2,3,6)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3089 : 10 : ii=ii+1 ; nabgauntselect(2,4,5)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3090 : 10 : ii=ii+1 ; nabgauntselect(2,5,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3091 : 10 : ii=ii+1 ; nabgauntselect(2,6,3)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3092 : 10 : ii=ii+1 ; nabgauntselect(2,7,2)=ii ; nablagaunt(ii)=-0.37846987830302403_dp!\frac{-3}{2\sqrt{5\pi}}
3093 : 10 : ii=ii+1 ; nabgauntselect(2,9,2)=ii ; nablagaunt(ii)=-0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
3094 : 10 : ii=ii+1 ; nabgauntselect(3,2,6)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3095 : 10 : ii=ii+1 ; nabgauntselect(3,3,7)=ii ; nablagaunt(ii)=0.75693974607408354_dp !\frac{3}{\sqrt{5\pi}
3096 : 10 : ii=ii+1 ; nabgauntselect(3,4,8)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3097 : 10 : ii=ii+1 ; nabgauntselect(3,6,2)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3098 : 10 : ii=ii+1 ; nabgauntselect(3,7,3)=ii ; nablagaunt(ii)=0.7569397566060481_dp !\frac{3}{\sqrt{5\pi}}
3099 : 10 : ii=ii+1 ; nabgauntselect(3,8,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
3100 : 10 : ii=ii+1 ; nabgauntselect(4,2,5)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3101 : 10 : ii=ii+1 ; nabgauntselect(4,3,8)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3102 : 10 : ii=ii+1 ; nabgauntselect(4,4,7)=ii ; nablagaunt(ii)=-0.37846987830302403_dp !\frac{-3}{2\sqrt{5\pi}}
3103 : 10 : ii=ii+1 ; nabgauntselect(4,4,9)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3104 : 10 : ii=ii+1 ; nabgauntselect(4,5,2)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
3105 : 10 : ii=ii+1 ; nabgauntselect(4,7,4)=ii ; nablagaunt(ii)=-0.37846987830302403_dp!\frac{-3}{2\sqrt{5\pi}}
3106 : 10 : ii=ii+1 ; nabgauntselect(4,8,3)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
3107 : 10 : ii=ii+1 ; nabgauntselect(4,9,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
3108 : : end if
3109 : : end if
3110 : :
3111 [ + - ]: 16 : if (l_max>2) then
3112 [ + - ]: 16 : if (l_max_ij>1) then
3113 : 16 : ii=ii+1 ; nabgauntselect(5,2,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3114 : 16 : ii=ii+1 ; nabgauntselect(5,4,2)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3115 : 16 : ii=ii+1 ; nabgauntselect(6,2,3)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3116 : 16 : ii=ii+1 ; nabgauntselect(6,3,2)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3117 : 16 : ii=ii+1 ; nabgauntselect(7,2,2)=ii ; nablagaunt(ii)=0.126156626101008_dp !\frac{1}{2\sqrt{5\pi}}
3118 : 16 : ii=ii+1 ; nabgauntselect(7,3,3)=ii ; nablagaunt(ii)=-0.252313252202016_dp !-\frac{1}{\sqrt{5\pi}}
3119 : 16 : ii=ii+1 ; nabgauntselect(7,4,4)=ii ; nablagaunt(ii)=0.126156626101008_dp !\frac{1}{2\sqrt{5\pi}}
3120 : 16 : ii=ii+1 ; nabgauntselect(8,3,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3121 : 16 : ii=ii+1 ; nabgauntselect(8,4,3)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3122 : 16 : ii=ii+1 ; nabgauntselect(9,2,2)=ii ; nablagaunt(ii)=0.2185096861184158_dp !\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3123 : 16 : ii=ii+1 ; nabgauntselect(9,4,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3124 : : end if
3125 [ + + ]: 16 : if (l_max_ij>2) then
3126 : 10 : ii=ii+1 ; nabgauntselect(5,5,7)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3127 : 10 : ii=ii+1 ; nabgauntselect(5,6,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3128 : 10 : ii=ii+1 ; nabgauntselect(5,7,5)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3129 : 10 : ii=ii+1 ; nabgauntselect(5,8,6)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3130 : 10 : ii=ii+1 ; nabgauntselect(6,5,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3131 : 10 : ii=ii+1 ; nabgauntselect(6,6,7)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3132 : 10 : ii=ii+1 ; nabgauntselect(6,6,9)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !-\frac{3}{14}\sqrt{\frac{15}{\pi}}
3133 : 10 : ii=ii+1 ; nabgauntselect(6,7,6)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3134 : 10 : ii=ii+1 ; nabgauntselect(6,8,5)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3135 : 10 : ii=ii+1 ; nabgauntselect(6,9,6)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !-\frac{3}{14}\sqrt{\frac{15}{\pi}}
3136 : 10 : ii=ii+1 ; nabgauntselect(7,5,5)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3137 : 10 : ii=ii+1 ; nabgauntselect(7,6,6)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3138 : 10 : ii=ii+1 ; nabgauntselect(7,7,7)=ii ; nablagaunt(ii)=0.5406712547186058_dp !\frac{3}{7}\sqrt{\frac{5}{\pi}}
3139 : 10 : ii=ii+1 ; nabgauntselect(7,8,8)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3140 : 10 : ii=ii+1 ; nabgauntselect(7,9,9)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3141 : 10 : ii=ii+1 ; nabgauntselect(8,5,6)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3142 : 10 : ii=ii+1 ; nabgauntselect(8,6,5)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3143 : 10 : ii=ii+1 ; nabgauntselect(8,7,8)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3144 : 10 : ii=ii+1 ; nabgauntselect(8,8,7)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3145 : 10 : ii=ii+1 ; nabgauntselect(8,8,9)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3146 : 10 : ii=ii+1 ; nabgauntselect(8,9,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3147 : 10 : ii=ii+1 ; nabgauntselect(9,6,6)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3148 : 10 : ii=ii+1 ; nabgauntselect(9,7,9)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3149 : 10 : ii=ii+1 ; nabgauntselect(9,8,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3150 : 10 : ii=ii+1 ; nabgauntselect(9,9,7)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3151 : : end if
3152 : : end if
3153 : :
3154 : 16 : nnablagnt=ii
3155 : :
3156 : : !If not tabulated, compute the integrals
3157 [ + + - + ]: 16 : if (l_max>3.or.l_max_ij>3) then
3158 : :
3159 : 10 : ntheta=25 ; nphi=25
3160 : 10 : call ylm_angular_mesh(ntheta,nphi,angl_size,cart_coord,ang_wgth)
3161 : :
3162 : 10 : mpsang=1+max(l_max,l_max_ij)
3163 : 10 : ylm_size=mpsang**2
3164 [ - + + - : 40 : LIBPAW_ALLOCATE(ylmr,(ylm_size,angl_size))
- + + - +
- + - - +
- + ]
3165 [ - + + - : 40 : LIBPAW_ALLOCATE(ylmrgr,(3,ylm_size,angl_size))
- + + - +
- + - - +
- + ]
3166 : 10 : call initylmr(mpsang,0,angl_size,ang_wgth,2,cart_coord,ylmr,ylmr_gr=ylmrgr)
3167 : :
3168 : : if (debug) open(unit=111,file='nablarealgaunt.dat',form='formatted')
3169 : :
3170 [ + + ]: 260 : do ilm=1,l_max**2
3171 [ + + ]: 2510 : do ilm_i=1,l_max_ij**2
3172 [ + + ]: 22750 : do ilm_j=1,l_max_ij**2
3173 : :
3174 [ + + - + ]: 20250 : if (ilm<10.and.ilm_i<10.and.ilm_j<10) cycle ! Already stored (tabulated)
3175 : :
3176 : : ! Compute integral
3177 : : nabla_rg=zero
3178 [ + + ]: 8112960 : do ipt=1,angl_size
3179 : : yylmgr=ylmrgr(1,ilm_i,ipt)*ylmrgr(1,ilm_j,ipt) &
3180 : : & +ylmrgr(2,ilm_i,ipt)*ylmrgr(2,ilm_j,ipt) &
3181 : 8100000 : & +ylmrgr(3,ilm_i,ipt)*ylmrgr(3,ilm_j,ipt)
3182 : 8112960 : nabla_rg=nabla_rg+ang_wgth(ipt)*ylmr(ilm,ipt)*yylmgr
3183 : : end do
3184 : 12960 : nabla_rg=four_pi*nabla_rg
3185 : :
3186 : : ! Store it if non-zero
3187 [ + + ]: 15210 : if (abs(nabla_rg)>tol12) then
3188 : : if (debug) then
3189 : : write(111,'(5x,a,i2,a,i2,a,i2,a,f19.15,a)') &
3190 : : & "ii=ii+1 ; nabgauntselect(",ilm,",",ilm_i,",",ilm_j, &
3191 : : & ")=ii ; nablagaunt(ii)=",nabla_rg,"_dp"
3192 : : end if
3193 : 790 : nnablagnt=nnablagnt+1
3194 : 790 : nabgauntselect(ilm,ilm_i,ilm_j)=nnablagnt
3195 : 790 : nablagaunt(nnablagnt)=nabla_rg
3196 : : end if
3197 : :
3198 : : end do ! ilm_j
3199 : : end do ! ilm_i
3200 : : end do !ilm
3201 : :
3202 : : if (debug) close(111)
3203 [ - + ]: 10 : LIBPAW_DEALLOCATE(ylmr)
3204 [ - + ]: 10 : LIBPAW_DEALLOCATE(ylmrgr)
3205 [ - + ]: 10 : LIBPAW_DEALLOCATE(cart_coord)
3206 [ - + ]: 10 : LIBPAW_DEALLOCATE(ang_wgth)
3207 : : end if
3208 : :
3209 [ - + - + ]: 16 : end subroutine nablarealgaunt
3210 : : !!***
3211 : :
3212 : : !----------------------------------------------------------------------
3213 : :
3214 : : !!****f* m_paw_sphharm/create_slm2ylm
3215 : : !! NAME
3216 : : !! create_slm2ylm
3217 : : !!
3218 : : !! FUNCTION
3219 : : !! For a given angular momentum lcor, compute slm2ylm.
3220 : : !!
3221 : : !! INPUTS
3222 : : !! lcor= angular momentum, size of the matrix is 2(2*lcor+1)
3223 : : !!
3224 : : !! OUTPUT
3225 : : !! slm2ylm(2lcor+1,2lcor+1) = rotation matrix.
3226 : : !!
3227 : : !! NOTES
3228 : : !! useful only in ndij==4
3229 : : !!
3230 : : !! SOURCE
3231 : :
3232 : : subroutine create_slm2ylm(lcor,slmtwoylm)
3233 : :
3234 : : !Arguments ---------------------------------------------
3235 : : !scalars
3236 : : integer,intent(in) :: lcor
3237 : : !arrays
3238 : : complex(dp),intent(out) :: slmtwoylm(2*lcor+1,2*lcor+1)
3239 : :
3240 : : !Local variables ---------------------------------------
3241 : : !scalars
3242 : : integer :: jm,ll,mm,im
3243 : : real(dp),parameter :: invsqrt2=one/sqrt2
3244 : : real(dp) :: onem
3245 : : ! *********************************************************************
3246 : :
3247 : : ll=lcor
3248 : : slmtwoylm=czero
3249 : : do im=1,2*ll+1
3250 : : mm=im-ll-1;jm=-mm+ll+1
3251 : : onem=dble((-1)**mm)
3252 : : if (mm> 0) then
3253 : : slmtwoylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
3254 : : slmtwoylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
3255 : : end if
3256 : : if (mm==0) then
3257 : : slmtwoylm(im,im)=cone
3258 : : end if
3259 : : if (mm< 0) then
3260 : : slmtwoylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
3261 : : slmtwoylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
3262 : : end if
3263 : : end do
3264 : :
3265 : : end subroutine create_slm2ylm
3266 : : !!***
3267 : :
3268 : : !----------------------------------------------------------------------
3269 : :
3270 : : !!****f* m_paw_sphharm/create_mlms2jmj
3271 : : !! NAME
3272 : : !! create_mlms2jmj
3273 : : !!
3274 : : !! FUNCTION
3275 : : !! For a given angular momentum lcor, give the rotation matrix msml2jmj
3276 : : !!
3277 : : !! INPUTS
3278 : : !! lcor= angular momentum
3279 : : !!
3280 : : !! SIDE EFFECTS
3281 : : !! mlms2jmj= rotation matrix
3282 : : !!
3283 : : !! SOURCE
3284 : :
3285 : : subroutine create_mlms2jmj(lcor,mlmstwojmj)
3286 : :
3287 : : !Arguments ---------------------------------------------
3288 : : !scalars
3289 : : integer,intent(in) :: lcor
3290 : : !arrays
3291 : : complex(dp),intent(out) :: mlmstwojmj(2*(2*lcor+1),2*(2*lcor+1))
3292 : :
3293 : : !Local variables ---------------------------------------
3294 : : !scalars
3295 : : integer :: jc1,jj,jm,ll,ml1,ms1
3296 : : real(dp) :: invsqrt2lp1,xj,xmj
3297 : : character(len=500) :: msg
3298 : : !arrays
3299 : : integer, allocatable :: ind_msml(:,:)
3300 : : complex(dp),allocatable :: mat_mlms2(:,:)
3301 : : !*********************************************************************
3302 : :
3303 : : !--------------- Built indices + allocations
3304 : : ll=lcor
3305 : : mlmstwojmj=czero
3306 : : LIBPAW_BOUND2_ALLOCATE(ind_msml,BOUNDS(1,2),BOUNDS(-ll,ll))
3307 : : LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1)))
3308 : : mlmstwojmj=czero
3309 : : jc1=0
3310 : : do ms1=1,2
3311 : : do ml1=-ll,ll
3312 : : jc1=jc1+1
3313 : : ind_msml(ms1,ml1)=jc1
3314 : : end do
3315 : : end do
3316 : :
3317 : : !--------------- built mlmstwojmj
3318 : : !do jj=ll,ll+1 ! the physical value of j are ll-0.5,ll+0.5
3319 : : !xj(jj)=jj-0.5
3320 : : if(ll==0)then
3321 : : msg=' ll should not be equal to zero !'
3322 : : LIBPAW_BUG(msg)
3323 : : end if
3324 : : jc1=0
3325 : : invsqrt2lp1=one/sqrt(float(2*lcor+1))
3326 : : do jj=ll,ll+1
3327 : : xj=float(jj)-half
3328 : : do jm=-jj,jj-1
3329 : : xmj=float(jm)+half
3330 : : jc1=jc1+1
3331 : : if(nint(xj+0.5)==ll+1) then
3332 : : if(nint(xmj+0.5)==ll+1) then
3333 : : mlmstwojmj(ind_msml(2,ll),jc1)=1.0 ! J=L+0.5 and m_J=L+0.5
3334 : : else if(nint(xmj-0.5)==-ll-1) then
3335 : : mlmstwojmj(ind_msml(1,-ll),jc1)=1.0 ! J=L+0.5 and m_J=-L-0.5
3336 : : else
3337 : : mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
3338 : : mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
3339 : : end if
3340 : : end if
3341 : : if(nint(xj+0.5)==ll) then
3342 : : mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
3343 : : mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
3344 : : end if
3345 : : end do
3346 : : end do
3347 : :
3348 : : LIBPAW_DEALLOCATE(ind_msml)
3349 : : LIBPAW_DEALLOCATE(mat_mlms2)
3350 : :
3351 : : end subroutine create_mlms2jmj
3352 : : !!***
3353 : :
3354 : : !----------------------------------------------------------------------
3355 : :
3356 : : !!****f* m_paw_sphharm/mkeuler
3357 : : !! NAME
3358 : : !! mkeuler
3359 : : !!
3360 : : !! FUNCTION
3361 : : !! Private function
3362 : : !! For a given symmetry operation, determines the corresponding Euler angles
3363 : : !!
3364 : : !! INPUTS
3365 : : !! rot(3,3)= symmetry matrix
3366 : : !!
3367 : : !! OUTPUT
3368 : : !! cosalp= cos(alpha) with alpha=Euler angle 1
3369 : : !! cosbeta= cos(beta) with beta =Euler angle 2
3370 : : !! cosgam= cos(gamma) with gamma=Euler angle 3
3371 : : !! isn= error code (0 if the routine exit normally)
3372 : : !! sinalp= sin(alpha) with alpha=Euler angle 1
3373 : : !! sinbeta= sin(beta) with beta=Euler angle 2
3374 : : !! singam= sin(gamma) with gamma=Euler angle 3
3375 : : !!
3376 : : !! NOTES
3377 : : !! This file comes from the file crystal_symmetry.f
3378 : : !! by N.A.W. Holzwarth and A. Tackett for the code pwpaw
3379 : : !! XG20200718 However, this routine was not accurate in the determination
3380 : : !! of beta when cosbeta was close to one (indeed this is a special case).
3381 : : !! This has been corrected. Moreover, sinbeta has been made an output in order
3382 : : !! to allow accurate calculations in dbeta. Also, tolerances have been made consistent.
3383 : : !!
3384 : : !! SOURCE
3385 : :
3386 : 27816 : subroutine mkeuler(rot,cosbeta,sinbeta,cosalp,sinalp,cosgam,singam,isn)
3387 : :
3388 : : !Arguments ---------------------------------------------
3389 : : !scalars
3390 : : integer,intent(out) :: isn
3391 : : real(dp),intent(out) :: cosalp,cosbeta,cosgam,sinalp,sinbeta,singam
3392 : : !arrays
3393 : : real(dp),intent(in) :: rot(3,3)
3394 : :
3395 : : !Local variables ---------------------------------------
3396 : : !scalars
3397 : : integer :: ier
3398 : : real(dp) :: check,sinbeta2
3399 : : character(len=500) :: msg
3400 : :
3401 : : ! *********************************************************************
3402 : :
3403 [ + - ]: 42473 : do isn= -1,1,2
3404 : :
3405 : : !Old coding, inaccurate
3406 : : ! cosbeta=real(isn)*rot(3,3)
3407 : : ! if(abs(1._dp-cosbeta*cosbeta)<tol10) then
3408 : : ! sinbeta=zero
3409 : : ! else
3410 : : ! sinbeta=sqrt(1._dp-cosbeta*cosbeta)
3411 : : ! end if
3412 : : ! if (abs(sinbeta).gt.tol10) then
3413 : : ! cosalp=isn*rot(3,1)/sinbeta
3414 : : ! sinalp=isn*rot(3,2)/sinbeta
3415 : : ! cosgam=-isn*rot(1,3)/sinbeta
3416 : : ! singam=isn*rot(2,3)/sinbeta
3417 : : ! else
3418 : : ! cosalp=isn*rot(1,1)/cosbeta
3419 : : ! sinalp=isn*rot(1,2)/cosbeta
3420 : : ! cosgam=one
3421 : : ! singam=zero
3422 : : ! end if
3423 : :
3424 : : !New coding, more accurate
3425 : 42473 : cosbeta=real(isn)*rot(3,3)
3426 : 42473 : sinbeta2=rot(1,3)**2+rot(2,3)**2
3427 [ + + ]: 42473 : if(sinbeta2<tol8**2)then
3428 : 17044 : sinbeta=zero
3429 : 17044 : cosalp=isn*rot(1,1)/cosbeta
3430 : 17044 : sinalp=isn*rot(1,2)/cosbeta
3431 : 17044 : cosgam=one
3432 : 17044 : singam=zero
3433 : : else
3434 : 25429 : sinbeta=sqrt(sinbeta2)
3435 : 25429 : cosalp=isn*rot(3,1)/sinbeta
3436 : 25429 : sinalp=isn*rot(3,2)/sinbeta
3437 : 25429 : cosgam=-isn*rot(1,3)/sinbeta
3438 : 25429 : singam=isn*rot(2,3)/sinbeta
3439 : : end if
3440 : : !
3441 : :
3442 : : ! Check matrix:
3443 : 42473 : ier=0
3444 : 42473 : check=cosalp*cosbeta*cosgam-sinalp*singam
3445 [ + + ]: 42473 : if (abs(check-isn*rot(1,1))>tol8) ier=ier+1
3446 : 42473 : check=sinalp*cosbeta*cosgam+cosalp*singam
3447 [ + + ]: 42473 : if (abs(check-isn*rot(1,2))>tol8) ier=ier+1
3448 : 42473 : check=-sinbeta*cosgam
3449 [ - + ]: 42473 : if (abs(check-isn*rot(1,3))>tol8) ier=ier+1
3450 : 42473 : check=-cosalp*cosbeta*singam-sinalp*cosgam
3451 [ + + ]: 42473 : if (abs(check-isn*rot(2,1))>tol8) ier=ier+1
3452 : 42473 : check=-sinalp*cosbeta*singam+cosalp*cosgam
3453 [ + + ]: 42473 : if (abs(check-isn*rot(2,2))>tol8) ier=ier+1
3454 : 42473 : check=sinbeta*singam
3455 [ - + ]: 42473 : if (abs(check-isn*rot(2,3))>tol8) ier=ier+1
3456 : 42473 : check=cosalp*sinbeta
3457 [ - + ]: 42473 : if (abs(check-isn*rot(3,1))>tol8) ier=ier+1
3458 : 42473 : check=sinalp*sinbeta
3459 [ - + ]: 42473 : if (abs(check-isn*rot(3,2))>tol8) ier=ier+1
3460 [ + + ]: 42473 : if (ier.eq.0) return
3461 : : end do
3462 : :
3463 : 0 : isn=0
3464 : : write(msg, '(7a)' )&
3465 : 0 : & 'Error during determination of symetries!',ch10,&
3466 : 0 : & 'Action: check your input file:',ch10,&
3467 : 0 : & 'unit cell vectors and/or atoms positions',ch10,&
3468 : 0 : & 'have to be given with a better precision.'
3469 : 0 : LIBPAW_ERROR(msg)
3470 : :
3471 : : end subroutine mkeuler
3472 : : !!***
3473 : :
3474 : : !----------------------------------------------------------------------
3475 : :
3476 : : !!****f* m_paw_sphharm/phim
3477 : : !! NAME
3478 : : !! phim
3479 : : !!
3480 : : !! FUNCTION
3481 : : !! Computes Phi_m[theta]=Sqrt[2] cos[m theta], if m>0
3482 : : !! Sqrt[2] sin[Abs(m) theta], if m<0
3483 : : !! 1 , if m=0
3484 : : !!
3485 : : !! INPUTS
3486 : : !! costeta= cos(theta) (theta= input angle)
3487 : : !! mm = index m
3488 : : !! sinteta= sin(theta) (theta= input angle)
3489 : : !!
3490 : : !! OUTPUT
3491 : : !! phim= Phi_m(theta) (see above)
3492 : : !!
3493 : : !! NOTES
3494 : : !! - This file comes from the file crystal_symmetry.f
3495 : : !! by N.A.W. Holzwarth and A. Tackett for the code pwpaw
3496 : : !!
3497 : : !! SOURCE
3498 : :
3499 : 2403936 : pure function phim(costheta,sintheta,mm)
3500 : :
3501 : : !Arguments ---------------------------------------------
3502 : : !scalars
3503 : : integer,intent(in) :: mm
3504 : : real(dp) :: phim
3505 : : real(dp),intent(in) :: costheta,sintheta
3506 : :
3507 : : ! *********************************************************************
3508 : :
3509 : : if (mm==0) phim=one
3510 [ + + ]: 2403936 : if (mm==1) phim=sqrt2*costheta
3511 [ + + ]: 2403936 : if (mm==-1) phim=sqrt2*sintheta
3512 [ + + ]: 2403936 : if (mm==2) phim=sqrt2*(costheta*costheta-sintheta*sintheta)
3513 [ + + ]: 2403936 : if (mm==-2) phim=sqrt2*two*sintheta*costheta
3514 [ + + ]: 2403936 : if (mm==3) phim=sqrt2*&
3515 : : & (costheta*(costheta*costheta-sintheta*sintheta)&
3516 : 27576 : & -sintheta*two*sintheta*costheta)
3517 [ + + ]: 2403936 : if (mm==-3) phim=sqrt2*&
3518 : : & (sintheta*(costheta*costheta-sintheta*sintheta)&
3519 : 27576 : & +costheta*two*sintheta*costheta)
3520 : :
3521 : 2403936 : end function phim
3522 : : !!***
3523 : :
3524 : : !----------------------------------------------------------------------
3525 : :
3526 : : !!****f* m_paw_sphharm/dbeta
3527 : : !! NAME
3528 : : !! dbeta
3529 : : !!
3530 : : !! FUNCTION
3531 : : !! Calculate the rotation matrix d^l_{m{\prim}m}(beta) using Eq. 4.14 of
3532 : : !! M.E. Rose, Elementary Theory of Angular Momentum,
3533 : : !! John Wiley & Sons, New-York, 1957
3534 : : !!
3535 : : !! INPUTS
3536 : : !! cosbeta= cosinus of beta (=Euler angle)
3537 : : !! sinbeta= sinus of beta (=Euler angle)
3538 : : !! ll= index l
3539 : : !! mm= index m
3540 : : !! mp= index m_prime
3541 : : !!
3542 : : !! OUTPUT
3543 : : !! dbeta= rotation matrix
3544 : : !!
3545 : : !! NOTES
3546 : : !! - This file comes from the file crystal_symmetry.f
3547 : : !! by N.A.W. Holzwarth and A. Tackett for the code pwpaw
3548 : : !! - Assume l relatively small so that factorials do not cause
3549 : : !! roundoff error
3550 : : !! - XG20200718 This routine was inaccurate when cosbeta was close to one or minus one.
3551 : : !! This has been fixed by adding sinbeta argument obtained from mkeuler.
3552 : : !! Tolerances have been adjusted as well.
3553 : : !!
3554 : : !! SOURCE
3555 : :
3556 : 2403936 : function dbeta(cosbeta,sinbeta,ll,mp,mm)
3557 : :
3558 : : !Arguments ---------------------------------------------
3559 : : !scalars
3560 : : integer,intent(in) :: ll,mm,mp
3561 : : real(dp) :: dbeta
3562 : : real(dp),intent(in) :: cosbeta,sinbeta
3563 : :
3564 : : !Local variables ------------------------------
3565 : : !scalars
3566 : : integer,parameter :: mxterms=200
3567 : : integer :: ii,ina,inb,inc,ml,ms
3568 : : real(dp) :: arg,cosbetab2,pref,sinbetab2,sum,tt
3569 : :
3570 : : !************************************************************************
3571 : 2403936 : dbeta=zero
3572 : :
3573 : : !Special cases
3574 [ + + ]: 2403936 : if (abs(cosbeta-1._dp).lt.tol10) then
3575 [ + + ]: 460920 : if (mp.eq.mm) dbeta=1
3576 [ + + ]: 1943016 : else if (abs(cosbeta+1._dp).lt.tol10) then
3577 [ + + ]: 424800 : if (mp.eq.-mm) dbeta=(-1)**(ll+mm)
3578 : : else
3579 : : ! General case
3580 : :
3581 : : !!!!! Old coding
3582 : : !! This is inaccurate when cosbeta is close to -1
3583 : : ! cosbetab2=sqrt((1+cosbeta)*0.5_dp)
3584 : : !! This is inaccurate when cosbeta is close to +1
3585 : : ! sinbetab2=sqrt((1-cosbeta)*0.5_dp)
3586 : : !!!!! End old coding, begin new coding
3587 [ + + ]: 1518216 : if(cosbeta>-tol8)then
3588 : : !If cosbeta is positive, cosbeta2 is positive with value >0.7, so one can divide by cosbetab2
3589 : 1496456 : cosbetab2=sqrt((1+cosbeta)*half)
3590 : 1496456 : sinbetab2=sinbeta*half/cosbetab2
3591 : : else
3592 : : !If cosbeta is negative, sinbeta2 is positive with value >0.7, so one can divide by sinbetab2
3593 : 21760 : sinbetab2=sqrt((1-cosbeta)*half)
3594 : 21760 : cosbetab2=sinbeta*half/sinbetab2
3595 : : endif
3596 : : !!!!! End of new coding
3597 : :
3598 : 1518216 : ml=max(mp,mm)
3599 : 1518216 : ms=min(mp,mm)
3600 [ + + ]: 1518216 : if (ml.ne.mp) sinbetab2=-sinbetab2
3601 : 1518216 : tt=-(sinbetab2/cosbetab2)**2
3602 : : pref=sqrt((rfactorial(ll-ms)*rfactorial(ll+ml))&
3603 : : & /(rfactorial(ll+ms)*rfactorial(ll-ml)))&
3604 : : & /rfactorial(ml-ms)*(cosbetab2**(2*ll+ms-ml))&
3605 : 12128744 : & *((-sinbetab2)**(ml-ms))
3606 : 1518216 : sum=1._dp
3607 : 1518216 : arg=1._dp
3608 : 1518216 : ina=ml-ll
3609 : 1518216 : inb=-ms-ll
3610 : 1518216 : inc=ml-ms+1
3611 [ + - ]: 1992392 : do ii=1,mxterms
3612 [ + + ]: 1992392 : if (ina.eq.0.or.inb.eq.0) exit
3613 : 474176 : arg=(arg*ina*inb*tt)/(ii*inc)
3614 : 474176 : sum=sum+arg
3615 : 474176 : ina=ina+1
3616 : 474176 : inb=inb+1
3617 : 1992392 : inc=inc+1
3618 : : end do
3619 : 1518216 : dbeta=pref*sum
3620 : : end if
3621 : :
3622 : 2403936 : end function dbeta
3623 : : !!***
3624 : :
3625 : : !----------------------------------------------------------------------
3626 : :
3627 : : !!****f* m_paw_sphharm/gauleg
3628 : : !! NAME
3629 : : !! gauleg
3630 : : !!
3631 : : !! FUNCTION
3632 : : !! Private function
3633 : : !! Compute the coefficients (supports and weights) for Gauss-Legendre integration
3634 : : !!
3635 : : !! INPUTS
3636 : : !! xmin=lower bound of integration
3637 : : !! xmax=upper bound of integration
3638 : : !! nn=order of integration
3639 : : !!
3640 : : !! OUTPUT
3641 : : !! x(nn)=array of support points
3642 : : !! weights(n)=array of integration weights
3643 : : !!
3644 : : !! SOURCE
3645 : :
3646 : 63 : subroutine gauleg(xmin,xmax,x,weights,nn)
3647 : :
3648 : : !Arguments ---------------------------------------------
3649 : : !scalars
3650 : : integer,intent(in) :: nn
3651 : : real(dp),intent(in) :: xmax,xmin
3652 : : !arrays
3653 : : real(dp),intent(out) :: weights(nn),x(nn)
3654 : :
3655 : : !Local variables ------------------------------
3656 : : !scalars
3657 : : integer :: ii,jj
3658 : : real(dp),parameter :: tol=1.d-13
3659 : : real(dp) :: p1,p2,p3,pi,xl,pp,xmean,z,z1
3660 : :
3661 : : !************************************************************************
3662 : :
3663 : 63 : pi=4._dp*atan(1._dp)
3664 : 63 : xl=(xmax-xmin)*0.5_dp
3665 : 63 : xmean=(xmax+xmin)*0.5_dp
3666 : :
3667 [ + + ]: 588 : do ii=1,(nn+1)/2
3668 : 525 : z=cos(pi*(ii-0.25_dp)/(nn+0.5_dp))
3669 : : do
3670 : 1732 : p1=1._dp
3671 : 1732 : p2=0._dp
3672 [ + + ]: 36608 : do jj=1,nn
3673 : 34876 : p3=p2
3674 : 34876 : p2=p1
3675 : 36608 : p1=((2._dp*jj-1._dp)*z*p2-(jj-1._dp)*p3)/jj
3676 : : end do
3677 : 1732 : pp=nn*(p2-z*p1)/(1._dp-z**2)
3678 : 1732 : z1=z
3679 : 1732 : z=z1-p1/pp
3680 [ + + ]: 1732 : if(abs(z-z1) < tol) exit
3681 : : end do
3682 : 525 : x(ii)=xmean-xl*z
3683 : 525 : x(nn+1-ii)=xmean+xl*z
3684 : 525 : weights(ii)=2._dp*xl/((1._dp-z**2)*pp**2)
3685 : 588 : weights(nn+1-ii)=weights(ii)
3686 : : end do
3687 : :
3688 : 63 : end subroutine gauleg
3689 : : !!***
3690 : :
3691 : : !----------------------------------------------------------------------
3692 : :
3693 : : !!****f* m_paw_sphharm/perms
3694 : : !! NAME
3695 : : !! perms
3696 : : !!
3697 : : !! FUNCTION
3698 : : !! Private function
3699 : : !! Returns N!/(N-k)! if N>=0 and N>k ; otherwise 0 is returned
3700 : : !!
3701 : : !! INPUTS
3702 : : !! kk=number k to use
3703 : : !! nn=number N to use
3704 : : !!
3705 : : !! OUTPUT
3706 : : !! perms= n!/(n-k)!
3707 : : !!
3708 : : !! SOURCE
3709 : :
3710 : 1709726 : function perms(nn,kk)
3711 : :
3712 : : !Arguments ---------------------------------------------
3713 : : !scalars
3714 : : integer,intent(in) :: kk,nn
3715 : : real(dp) :: perms
3716 : :
3717 : : !Local variables ---------------------------------------
3718 : : !scalars
3719 : : integer :: ii
3720 : : real(dp) :: pp
3721 : :
3722 : : ! *********************************************************************
3723 : :
3724 [ + - + - ]: 1709726 : if (nn>=0.and.nn>=kk) then
3725 : 1709726 : pp=1._dp
3726 [ + + ]: 4710492 : do ii=nn-kk+1,nn
3727 : 4710492 : pp=pp*ii
3728 : : end do
3729 : : else
3730 : : pp=0._dp
3731 : : end if
3732 : :
3733 : 1709726 : perms=pp
3734 : :
3735 : 1709726 : end function perms
3736 : : !!***
3737 : :
3738 : : !----------------------------------------------------------------------
3739 : :
3740 : : !!****f* m_paw_sphharm/rfactorial
3741 : : !! NAME
3742 : : !! rfactorial
3743 : : !!
3744 : : !! FUNCTION
3745 : : !! Private function
3746 : : !! Calculates N! as a double precision real.
3747 : : !!
3748 : : !! INPUTS
3749 : : !! nn=number to use
3750 : : !!
3751 : : !! OUTPUT
3752 : : !! factorial= n! (real)
3753 : : !!
3754 : : !! SOURCE
3755 : :
3756 : 7844014 : elemental function rfactorial(nn)
3757 : :
3758 : : !Arguments ---------------------------------------------
3759 : : !scalars
3760 : : integer,intent(in) :: nn
3761 : : real(dp) :: rfactorial
3762 : :
3763 : : !Local variables ---------------------------------------
3764 : : !scalars
3765 : : integer :: ii
3766 : :
3767 : : ! *********************************************************************
3768 : :
3769 : 7844014 : rfactorial=one
3770 [ + + + + : 24606512 : do ii=2,nn
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ ]
3771 : 24606512 : rfactorial=rfactorial*ii
3772 : : end do
3773 : :
3774 : : end function rfactorial
3775 : : !!***
3776 : :
3777 : : END MODULE m_paw_sphharm
3778 : : !!***
|
