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 191287784 : 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 191287784 : 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 191287784 : ylmc = czero
124 :
125 191287784 : r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
126 191287784 : if (r<PPAD) r=r+PPAD
127 : !$if (r<tol10) RETURN
128 :
129 191287784 : rxy=SQRT(kcart(1)**2+kcart(2)**2)
130 191287784 : if (rxy<PPAD)rxy=r+PPAD
131 : !
132 : ! Determine theta and phi
133 191287784 : costh= kcart(3)/r
134 :
135 : #if 1
136 : ! old buggy coding
137 191287784 : sinth= rxy/r
138 191287784 : cosphi= kcart(1)/rxy
139 191287784 : 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 191287784 : costwophi= two*cosphi**2 - one
151 191287784 : sintwophi= two*sinphi*cosphi
152 191287784 : costhreephi=cosphi*costwophi-sinphi*sintwophi
153 191287784 : sinthreephi=cosphi*sintwophi+sinphi*costwophi
154 :
155 : select case (il)
156 :
157 : case (0)
158 271661384 : ylmc= one/SQRT(four_pi)
159 :
160 : case (1)
161 80373600 : if (ABS(im)==0) then
162 26791200 : ylmc = SQRT(three/(four_pi))*costh
163 53582400 : else if (ABS(im)==1) then
164 53582400 : 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 191287784 : 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 191287784 : if (im < 0) then
203 60376464 : ctmp = (-one)**(im)*CONJG(ylmc)
204 60376464 : 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 191287784 : 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 47821946 : 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 47821946 : 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 47821946 : dphi=czero; dth=czero
280 :
281 47821946 : r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
282 47821946 : if (r<PPAD) r=r+PPAD
283 : !$if (r<tol10) RETURN
284 :
285 47821946 : rxy=SQRT(kcart(1)**2+kcart(2)**2)
286 47821946 : if (rxy<PPAD) rxy=r+PPAD
287 :
288 : ! Determine theta and phi
289 47821946 : costh= kcart(3)/r
290 : #if 1
291 : ! old buggy coding
292 47821946 : sinth= rxy/r
293 47821946 : cosphi= kcart(1)/rxy
294 47821946 : 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 47821946 : costwophi= two*cosphi**2 - one
306 47821946 : sintwophi= two*sinphi*cosphi
307 47821946 : costhreephi=cosphi*costwophi-sinphi*sintwophi
308 47821946 : sinthreephi=cosphi*sintwophi+sinphi*costwophi
309 :
310 : select case (il)
311 :
312 : case (0)
313 : dth = czero
314 20093400 : dphi = czero
315 :
316 : case (1)
317 20093400 : if (ABS(im)==0) then
318 6697800 : dth= -SQRT(three/(four_pi))*sinth
319 : dphi= czero
320 13395600 : else if (abs(im)==1) then
321 13395600 : dth= -SQRT(3.d0/(8.d0*pi))*costh*CMPLX(cosphi,sinphi)
322 13395600 : 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 47821946 : LIBPAW_ERROR(msg)
357 : end select
358 : !
359 : !=== Treat the case im < 0 ===
360 47821946 : if (im<0) then
361 15094116 : ctmp = (-one)**(im)*CONJG(dth)
362 15094116 : dth = ctmp
363 15094116 : ctmp= (-one)**(im)*CONJG(dphi)
364 15094116 : dphi= ctmp
365 : end if
366 :
367 47821946 : 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 434995 : 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 939828 : real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),rphase(mpsang-1)
554 469914 : real(dp),allocatable :: blm(:,:)
555 :
556 : !************************************************************************
557 :
558 : !What has to be computed ?
559 469914 : compute_ylm = (option==1.or.option==2.or.option==3)
560 469914 : compute_ylmgr =(( option==2.or.option==3).and.present(ylmr_gr))
561 469914 : compute_ylm2gr=(( option==3).and.present(ylmr_gr))
562 469914 : dimgr=3*(option/2)+6*(option/3)
563 :
564 : !Initialisation of spherical harmonics
565 514243186 : if (compute_ylm ) ylmr (: ,1:npts)=zero
566 616232963 : if (compute_ylmgr) ylmr_gr(:,:,1:npts)=zero
567 :
568 : !Special case for l=0
569 41440750 : if (compute_ylm ) ylmr(1,1:npts)=1._dp/sqrt(four_pi)
570 43240226 : if (compute_ylmgr) ylmr_gr(1:dimgr,1,1:npts)=zero
571 469914 : if (mpsang>1) then
572 :
573 : ! Loop over all rr
574 40608185 : do inpt=1,npts
575 :
576 : ! Load module of rr
577 40140581 : rnorm=one
578 40140581 : if (normchoice==1) rnorm=nrm(inpt)
579 :
580 : ! Continue only for r<>0
581 :
582 40169541 : if (rnorm>tol10) then
583 :
584 : ! Determine theta and phi
585 40103686 : cphi=one
586 40103686 : sphi=zero
587 40103686 : 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 40103686 : if (abs(one-ctheta)<tol12) then
593 116914 : ctheta=one
594 116914 : stheta=zero
595 39986772 : else if (abs(one+ctheta)<tol12) then
596 117032 : ctheta=-one
597 117032 : stheta=zero
598 : else
599 39869740 : stheta=sqrt(one-ctheta*ctheta)
600 39869740 : cphi=rr(1,inpt)/(rnorm*stheta)
601 39869740 : sphi=rr(2,inpt)/(rnorm*stheta)
602 : end if
603 132961988 : do mm=1,mpsang-1
604 92858302 : rphase(mm)=dreal(dcmplx(cphi,sphi)**mm)
605 132961988 : iphase(mm)=aimag(dcmplx(cphi,sphi)**mm)
606 : end do
607 :
608 : ! Determine gradients of theta and phi
609 40103686 : if (compute_ylmgr) then
610 9976190 : dtheta(1)=ctheta*cphi
611 9976190 : dtheta(2)=ctheta*sphi
612 9976190 : dtheta(3)=-stheta
613 9976190 : dphi(1)=-sphi
614 9976190 : dphi(2)=cphi
615 9976190 : dphi(3)=zero
616 : end if
617 :
618 : ! COMPUTE Ylm(R)
619 40103686 : if (compute_ylm) then
620 : ! Loop over angular momentum l
621 132961988 : do ilang=2,mpsang
622 92858302 : ll=ilang-1
623 92858302 : l0=ll**2+ll+1
624 92858302 : fact=1._dp/real(ll*(ll+1),dp)
625 92858302 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)
626 : ! Special case m=0
627 92858302 : ylmr(l0,inpt)=ylmcst*ass_leg_pol(ll,0,ctheta)
628 : ! Compute for m>0
629 92858302 : onem=one
630 302293188 : do mm=1,ll
631 169331200 : onem=-onem
632 169331200 : work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp)
633 169331200 : ylmr(l0+mm,inpt)=work1*rphase(mm)
634 169331200 : ylmr(l0-mm,inpt)=work1*iphase(mm)
635 262189502 : 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 40103686 : if (compute_ylmgr) then
642 : ! Loop over angular momentum l
643 36569774 : do ilang=2,mpsang
644 26593584 : ll=ilang-1
645 26593584 : l0=ll**2+ll+1
646 26593584 : fact=1._dp/real(ll*(ll+1),dp)
647 26593584 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rnorm
648 : ! Special case m=0
649 26593584 : work1=ylmcst*plm_dtheta(ll,0,ctheta)
650 106374336 : ylmr_gr(1:3,l0,inpt)=work1*dtheta(1:3)
651 : ! Compute for m>0
652 26593584 : onem=one
653 92113512 : do mm=1,ll
654 55543738 : onem=-onem
655 55543738 : work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp)
656 55543738 : work2=ylmcst*sqrt(fact)*onem*plm_dphi (ll,mm,ctheta)*sqrt(2._dp)
657 222174952 : ylmr_gr(1:3,l0+mm,inpt)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3)
658 222174952 : ylmr_gr(1:3,l0-mm,inpt)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3)
659 82137322 : 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 40103686 : if (compute_ylm2gr) then
666 421251 : LIBPAW_ALLOCATE(blm,(5,mpsang*mpsang))
667 140417 : call plm_coeff(blm,mpsang,ctheta)
668 :
669 : ! Loop over angular momentum l
670 522383 : do ilang=2,mpsang
671 381966 : ll=ilang-1
672 381966 : l0=ll**2+ll+1
673 381966 : fact=1._dp/real(ll*(ll+1),dp)
674 381966 : ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rnorm**2)
675 : ! Special case m=0
676 381966 : ylmr_gr(4,l0,inpt)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi)
677 381966 : ylmr_gr(5,l0,inpt)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi)
678 381966 : ylmr_gr(6,l0,inpt)=ylmcst*blm(1,l0)
679 381966 : ylmr_gr(7,l0,inpt)=ylmcst*blm(2,l0)*sphi
680 381966 : ylmr_gr(8,l0,inpt)=ylmcst*blm(2,l0)*cphi
681 381966 : ylmr_gr(9,l0,inpt)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi
682 : ! Compute for m>0
683 381966 : onem=one
684 1297440 : do mm=1,ll
685 775057 : 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 775057 : & 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 775057 : & 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 775057 : & 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 775057 : & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
694 775057 : ylmr_gr(6,l0+mm,inpt)=ylmcst2*blm(1,l0+mm)*rphase(mm)
695 775057 : 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 775057 : & 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 775057 : & 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 775057 : & 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 775057 : & 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 775057 : & 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 775057 : & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm))
708 1157023 : 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 140417 : 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 469914 : 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 503892 : 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 503892 : ys_val = czero
765 :
766 503892 : if ( l2 /= l1 ) return
767 503892 : if ( abs(m2) /= abs(m1) ) return
768 :
769 377540 : mp1=(-1)**abs(m1)
770 :
771 377540 : if(m1.EQ.0) then
772 21868 : ys_val=cone
773 355672 : else if((m1.GT.0).AND.(m2.GT.0)) then
774 91394 : ys_val=mp1*sqrthalf
775 264278 : else if((m1.GT.0).AND.(m2.LT.0)) then
776 91394 : ys_val=sqrthalf
777 172884 : else if((m1.LT.0).AND.(m2.GT.0)) then
778 86442 : ys_val=-j_dpc*mp1*sqrthalf
779 86442 : else if((m1.LT.0).AND.(m2.LT.0)) then
780 86442 : 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 193239 : 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 193239 : lidir = czero
822 193239 : if ( lp /= ll ) return
823 :
824 193239 : jpme=czero; jmme=czero; jme=czero
825 193239 : if (mp==mm) then
826 55413 : jme=cone*mm
827 137826 : else if (mp==mm+1) then
828 38742 : jpme=-cone*sqrt(half*((ll*(ll+1))-mm*(mm+1)))
829 99084 : else if (mp==mm-1) then
830 38742 : jmme= cone*sqrt(half*((ll*(ll+1))-mm*(mm-1)))
831 : end if
832 :
833 277623 : select case (idir)
834 : case (1) ! Lx
835 84384 : lidir = -sqrthalf*(jpme - jmme)
836 : case (2) ! Ly
837 45828 : lidir = j_dpc*sqrthalf*(jpme + jmme)
838 : case (3) ! Lz
839 193239 : 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 16671 : 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 16671 : sls_val = czero
881 :
882 16671 : if ( lp /= ll ) return
883 :
884 72084 : do mpp = -ll, ll
885 55413 : call ys(ll,mpp,ll,mp,sy_val)
886 265323 : do mppp = -ll, ll
887 193239 : call lxyz(ll,mpp,idir,ll,mppp,lidir)
888 193239 : call ys(ll,mppp,ll,mm,ys_val)
889 248652 : 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 66 : 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 66 : complex(dp),allocatable :: ls_cplx(:,:,:),slm2ylm(:,:)
942 66 : complex(dp),allocatable :: mat_inp_c(:,:,:),mat_out_c(:,:,:)
943 66 : 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 66 : if (.not.allocated(ls_ylm)) then
950 0 : msg='ls_ylm is not allocated!'
951 0 : LIBPAW_BUG(msg)
952 : end if
953 264 : 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 14154 : ls_ylm=zero
960 :
961 : !Nothing to do if lmax=0
962 66 : if (lmax<=0) return
963 :
964 : !Loop on l quantum number
965 172 : do ll=1,lmax
966 :
967 : ! Transformation matrixes: real->complex spherical harmonics
968 424 : LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1))
969 2210 : slm2ylm=czero
970 512 : do im=1,2*ll+1
971 406 : mm=im-ll-1;jm=-mm+ll+1
972 406 : onem=dble((-1)**mm)
973 406 : if (mm> 0) then
974 150 : slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
975 150 : slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
976 : end if
977 406 : if (mm==0) then
978 106 : slm2ylm(im,im)=cone
979 : end if
980 512 : if (mm< 0) then
981 150 : slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
982 150 : 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 530 : LIBPAW_ALLOCATE(ls_cplx,(2*ll+1,2*ll+1,2))
989 4526 : 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 512 : do im=1,2*ll+1
999 406 : mm=im-ll-1
1000 406 : 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 406 : if ((mm+1)<= ll) then
1004 300 : 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 512 : if ((mm-1)>=-ll) then
1009 300 : 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 106 : lm0=ll**2
1035 512 : do jm=1,2*ll+1
1036 406 : jlm=lm0+jm;j0lm=jlm*(jlm-1)/2
1037 1564 : do im=1,jm
1038 1052 : ilm=lm0+im;klm=j0lm+ilm
1039 1052 : tmp(:)=czero
1040 5812 : do ii=1,2*ll+1
1041 29640 : do jj=1,2*ll+1
1042 76244 : tmp(:)=tmp(:)+ls_cplx(ii,jj,:)*CONJG(slm2ylm(ii,im))*slm2ylm(jj,jm)
1043 : end do
1044 : end do
1045 3156 : ls_ylm(1,klm,:)=REAL(tmp(:),kind=dp)
1046 3562 : 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 106 : LIBPAW_DEALLOCATE(ls_cplx)
1081 172 : 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 4162200 : 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 8324400 : real(dp) :: pl_d2(mpsang),plm_d2t(mpsang*mpsang)
1125 :
1126 : !************************************************************************
1127 :
1128 4162200 : if (abs(xx).gt.1.d0) then
1129 0 : msg = ' plm_coeff : xx > 1 !'
1130 0 : LIBPAW_ERROR(msg)
1131 : end if
1132 :
1133 230334168 : blm=zero
1134 4162200 : is_one=(abs(abs(xx)-one)<=tol12)
1135 4162200 : xx2=xx**2
1136 4162200 : onemx2=abs(one-xx2)
1137 4162200 : sqrx=sqrt(onemx2)
1138 4162200 : xsqrx=xx*sqrt(onemx2)
1139 :
1140 4162200 : call plm_d2theta(mpsang,plm_d2t,xx)
1141 4162200 : if (is_one) then
1142 1527 : yy=sign(one,xx)
1143 1527 : call pl_deriv(mpsang,pl_d2,yy)
1144 : end if
1145 :
1146 16503140 : do il=0,mpsang-1
1147 12340940 : llp1=dble(il*(il+1))
1148 12340940 : ilm0=il*il+il+1
1149 41521274 : do im=0,il
1150 25018134 : ilm=ilm0+im;ilm1=ilm0-im
1151 :
1152 25018134 : plm =(-1)**im*ass_leg_pol(il,im,xx)
1153 25018134 : dplm_dt =(-1)**im*plm_dtheta(il,im,xx)
1154 25018134 : d2plm_dt2= plm_d2t(ilm)
1155 :
1156 25018134 : blm(1,ilm)= two*xsqrx *dplm_dt+onemx2*d2plm_dt2
1157 25018134 : blm(2,ilm)= (one-two*xx2)*dplm_dt-xsqrx *d2plm_dt2
1158 25018134 : blm(3,ilm)=llp1*plm+ d2plm_dt2
1159 25018134 : blm(4,ilm)= -two*xsqrx *dplm_dt+xx2 *d2plm_dt2
1160 :
1161 :
1162 25018134 : if (is_one) then
1163 8616 : if (im==1) then
1164 2766 : blm(5,ilm)=llp1*plm+d2plm_dt2
1165 : end if
1166 8616 : if (im==2) then
1167 1239 : blm(5,ilm)=d2plm_dt2-three*pl_d2(il+1)
1168 : end if
1169 : else
1170 25009518 : if(im>0) then
1171 12672871 : blm(5,ilm)=plm/onemx2-dplm_dt*xx/sqrx
1172 : end if
1173 : end if
1174 :
1175 37359074 : if (im>0) then
1176 12677194 : blm(1,ilm1)=blm(1,ilm)
1177 12677194 : blm(2,ilm1)=blm(2,ilm)
1178 12677194 : blm(3,ilm1)=blm(3,ilm)
1179 12677194 : blm(4,ilm1)=blm(4,ilm)
1180 12677194 : blm(5,ilm1)=blm(5,ilm)
1181 : end if
1182 :
1183 : end do
1184 : end do
1185 :
1186 4162200 : 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 383940126 : 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 383940126 : x=xarg
1226 383940126 : 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 383940126 : polmm=1.d0
1235 383940126 : if (m>0) then
1236 236273372 : sqrx=sqrt(abs((1.d0-x)*(1.d0+x)))
1237 613247643 : do i=1,m
1238 613247643 : polmm=polmm*(1.0d0-2.0d0*i)*sqrx
1239 : enddo
1240 : endif
1241 :
1242 383940126 : if (l==m) then
1243 : ass_leg_pol=polmm
1244 : else
1245 232256832 : tmp1=x*(2.0d0*m+1.0d0)*polmm
1246 232256832 : if (l==(m+1)) then
1247 : ass_leg_pol=tmp1
1248 : else
1249 237295663 : do ll=m+2,l
1250 140218985 : pll=(x*(2.0d0*ll-1.0d0)*tmp1-(ll+m-1.0d0)*polmm)/dble(ll-m)
1251 140218985 : polmm=tmp1
1252 237295663 : tmp1=pll
1253 : enddo
1254 : ass_leg_pol=pll
1255 : endif
1256 : endif
1257 :
1258 383940126 : 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 77064454 : 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 77064454 : 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 77064454 : plm_dphi=zero
1306 77064454 : if (mm==0) return
1307 :
1308 77064454 : pmm=one
1309 77064454 : dosomx2=one
1310 77064454 : if (mm > 0) then
1311 77064454 : somx2=sqrt((1-xx)*(1+xx))
1312 77064454 : fact=one
1313 210523059 : do im=1,mm
1314 133458605 : pmm=-pmm*fact
1315 210523059 : fact=fact+2
1316 : end do
1317 77064454 : if (mm > 1) then
1318 91844075 : do im=2,mm
1319 91844075 : dosomx2=somx2*dosomx2
1320 : end do
1321 : end if
1322 77064454 : pmm=pmm*dosomx2 !due to one more term (-1^M)
1323 : end if
1324 77064454 : if(ll==mm) then
1325 41614530 : plm_dphi=pmm*mm
1326 : else
1327 35449924 : pmmp1=xx*(2*mm+1)*pmm
1328 35449924 : if(ll==mm+1) then
1329 22652467 : plm_dphi=pmmp1*mm
1330 12797457 : else if(ll>=mm+2) then
1331 33741684 : do il=mm+2,ll
1332 20944227 : pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1333 20944227 : pmm=pmmp1
1334 33741684 : pmmp1=pll
1335 : end do
1336 12797457 : 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 215261828 : 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 215261828 : 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 215261828 : plm_dtheta=zero
1388 215261828 : pmm=one
1389 215261828 : dpmm=one
1390 215261828 : dosomx2=one
1391 215261828 : somx2=sqrt((1-xx)*(1+xx))
1392 215261828 : if(mm==0)then
1393 : dpmm=zero
1394 148965418 : elseif (mm > 0) then
1395 : fact=one
1396 380493115 : do im=1,mm
1397 231527697 : pmm=-pmm*fact*somx2
1398 231527697 : dpmm=-dpmm*fact
1399 380493115 : fact=fact+2
1400 : end do
1401 148965418 : if(mm>1)then
1402 141417891 : do im=2,mm
1403 141417891 : dosomx2=dosomx2*somx2
1404 : end do
1405 : end if
1406 148965418 : dpmm= dpmm*mm*xx*dosomx2
1407 : end if
1408 215261828 : if(ll==mm)then
1409 : plm_dtheta=dpmm
1410 : else
1411 112233918 : pmmp1=xx*(2*mm+1)*pmm
1412 112233918 : dpmmp1=-(2*mm+1)*somx2*pmm+xx*(2*mm+1)*dpmm
1413 112233918 : if(ll==mm+1) then
1414 : plm_dtheta=dpmmp1
1415 42357948 : else if(ll>=mm+2)then
1416 106444497 : do il=mm+2,ll
1417 64086549 : pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1418 64086549 : dpll=(-somx2*(2*il-1)*pmmp1+(xx*(2*il-1)*dpmmp1-(il+mm-1)*dpmm))/(il-mm)
1419 64086549 : pmm=pmmp1
1420 64086549 : pmmp1=pll
1421 64086549 : dpmm=dpmmp1
1422 106444497 : dpmmp1=dpll
1423 : end do
1424 : plm_dtheta=dpll
1425 : end if
1426 : end if
1427 :
1428 215261828 : 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 4162200 : 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 4162200 : if (abs(xx).gt.1.d0) then
1466 0 : msg = 'plm_d2theta : xx > 1 !'
1467 0 : LIBPAW_ERROR(msg)
1468 : end if
1469 :
1470 41857528 : plm_d2t=zero
1471 4162200 : if (mpsang>1) then
1472 4162200 : sqrx=sqrt(abs((1.d0-xx)*(1.d0+xx)))
1473 :
1474 12340940 : do il=1,mpsang-1
1475 8178740 : ilm=il*il+2*il+1
1476 8178740 : 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 8178740 : & 2.d0*xx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1480 8178740 : 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 8178740 : & 2.d0*sqrx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1484 12340940 : if(il>1) plm_d2t(il*il+2)=plm_d2t(ilm-1)
1485 : end do
1486 : ! terme d2(Plm)/dtet2
1487 4162200 : if(mpsang>2) then
1488 7601576 : do il=2,mpsang-1
1489 12100030 : do im=0,il-2
1490 4498454 : ilm=il*il+il+1+im
1491 4498454 : ilmm1=(il-1)*(il-1)+il+im
1492 4498454 : 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 4498454 : & dble(il+im-1)/dble(il-im)*plm_d2t(ilmm2)
1496 8514994 : plm_d2t(ilm-2*im)=plm_d2t(ilm)
1497 : end do
1498 : end do
1499 : end if
1500 : end if
1501 :
1502 4162200 : 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 1527 : 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 3054 : real(dp) :: pl(mpsang),pl_d1(mpsang)
1539 :
1540 : ! *********************************************************************
1541 :
1542 1527 : if (abs(xx).gt.1.d0) then
1543 0 : msg = 'pl_deriv : xx > 1 !'
1544 0 : LIBPAW_ERROR(msg)
1545 : end if
1546 :
1547 17460 : pl_d2=zero; pl_d1=zero; pl=zero
1548 1527 : pl(1)=one; pl(2)=xx
1549 1527 : pl_d1(1)=zero; pl_d1(2)=one
1550 1527 : pl_d2(1)=zero; pl_d2(2)=zero
1551 1527 : if (mpsang>2) then
1552 2266 : do il=2,mpsang-1
1553 1239 : il_=dble(il);il_m1=dble(il-1);il_2m1=dble(2*il-1)
1554 1239 : ilm=il+1
1555 1239 : pl(ilm)=(il_2m1*xx*pl(ilm-1)-il_m1*pl(ilm-2))/il_
1556 1239 : pl_d1(ilm)=(il_2m1*(xx*pl_d1(ilm-1)+pl(ilm-1))-il_m1*pl_d1(ilm-2))/il_
1557 2266 : 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 1527 : 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 1923 : 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 1923 : 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 1675118 : zarot=zero
2137 :
2138 29635 : do irot=1,nsym
2139 :
2140 27712 : 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 27712 : zarot(1,1,1,irot)=one
2147 :
2148 : ! === l>0 case ===
2149 27712 : if (lmax>0) then
2150 : ! Calculate the rotations in the cartesian basis
2151 27409 : rot=zero;prod=zero
2152 109636 : do k1=1,3
2153 356317 : do j1=1,3
2154 1068951 : do i1=1,3
2155 986724 : prod(i1,j1)=prod(i1,j1)+sym(i1,k1,irot)*rprimd(j1,k1)
2156 : end do
2157 : end do
2158 : end do
2159 109636 : do j1=1,3
2160 356317 : do i1=1,3
2161 986724 : do k1=1,3
2162 986724 : rot(i1,j1)=rot(i1,j1)+gprimd(i1,k1)*prod(k1,j1)
2163 : end do
2164 328908 : if(abs(rot(i1,j1))<tol10) rot(i1,j1)=zero
2165 : end do
2166 : end do
2167 27409 : call mkeuler(rot,cosbeta,sinbeta,cosalp,sinalp,cosgam,singam,isn)
2168 67702 : do ll=1,lmax
2169 40293 : il=(isn)**ll
2170 216353 : do mp=-ll,ll
2171 148651 : jj=mp+ll+1
2172 781965 : do mm=-ll,ll
2173 593021 : 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 741672 : & -(-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 29635 : 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 1923 : 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 1157 : 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 1157 : 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 1157 : ang_phipphj_tmp=zero
2449 : !
2450 1157 : ang_phipphj_tmp(1,4,1)=0.57735026918962576451_dp
2451 1157 : ang_phipphj_tmp(2,5,1)=0.44721359549995793928_dp
2452 1157 : ang_phipphj_tmp(3,8,1)=0.44721359549995793928_dp
2453 1157 : ang_phipphj_tmp(4,1,1)=0.57735026918962576451_dp
2454 1157 : ang_phipphj_tmp(4,7,1)=-0.25819888974716112568_dp
2455 1157 : ang_phipphj_tmp(4,9,1)=0.44721359549995793928_dp
2456 1157 : ang_phipphj_tmp(5,2,1)=0.44721359549995793928_dp
2457 1157 : ang_phipphj_tmp(5,10,1)=0.46291004988627573078_dp
2458 1157 : ang_phipphj_tmp(5,12,1)=-0.11952286093343936400_dp
2459 1157 : ang_phipphj_tmp(6,11,1)=0.37796447300922722721_dp
2460 1157 : ang_phipphj_tmp(7,4,1)=-0.25819888974716112568_dp
2461 1157 : ang_phipphj_tmp(7,14,1)=0.41403933560541253068_dp
2462 1157 : ang_phipphj_tmp(8,3,1)=0.44721359549995793928_dp
2463 1157 : ang_phipphj_tmp(8,13,1)=-0.29277002188455995381_dp
2464 1157 : ang_phipphj_tmp(8,15,1)=0.37796447300922722721_dp
2465 1157 : ang_phipphj_tmp(9,4,1)=0.44721359549995793928_dp
2466 1157 : ang_phipphj_tmp(9,14,1)=-0.11952286093343936400_dp
2467 1157 : ang_phipphj_tmp(9,16,1)=0.46291004988627573078_dp
2468 1157 : ang_phipphj_tmp(10,5,1)=0.46291004988627573078_dp
2469 1157 : ang_phipphj_tmp(11,6,1)=0.37796447300922722721_dp
2470 1157 : ang_phipphj_tmp(12,5,1)=-0.11952286093343936400_dp
2471 1157 : ang_phipphj_tmp(13,8,1)=-0.29277002188455995381_dp
2472 1157 : ang_phipphj_tmp(14,7,1)=0.41403933560541253068_dp
2473 1157 : ang_phipphj_tmp(14,9,1)=-0.11952286093343936400_dp
2474 1157 : ang_phipphj_tmp(15,8,1)=0.37796447300922722721_dp
2475 1157 : ang_phipphj_tmp(16,9,1)=0.46291004988627573078_dp
2476 : !
2477 1157 : ang_phipphj_tmp(1,4,2)=0.28867513459481288225_dp
2478 1157 : ang_phipphj_tmp(1,14,2)=-0.54006172486732168591_dp
2479 1157 : ang_phipphj_tmp(2,5,2)=0.22360679774997896964_dp
2480 1157 : ang_phipphj_tmp(3,8,2)=0.22360679774997896964_dp
2481 1157 : ang_phipphj_tmp(4,7,2)=-0.77459666924148337704_dp
2482 1157 : ang_phipphj_tmp(4,9,2)=0.22360679774997896964_dp
2483 1157 : ang_phipphj_tmp(5,2,2)=0.11180339887498948482_dp
2484 1157 : ang_phipphj_tmp(5,10,2)=0.23145502494313786539_dp
2485 1157 : ang_phipphj_tmp(5,12,2)=-0.47809144373375745599_dp
2486 1157 : ang_phipphj_tmp(6,11,2)=0.18898223650461361361_dp
2487 1157 : ang_phipphj_tmp(7,4,2)=0.25819888974716112568_dp
2488 1157 : ang_phipphj_tmp(7,14,2)=0.44854261357253024157_dp
2489 1157 : ang_phipphj_tmp(8,3,2)=-0.44721359549995793928_dp
2490 1157 : ang_phipphj_tmp(8,13,2)=-1.1710800875382398152_dp
2491 1157 : ang_phipphj_tmp(8,15,2)=0.18898223650461361361_dp
2492 1157 : ang_phipphj_tmp(9,4,2)=0.11180339887498948482_dp
2493 1157 : ang_phipphj_tmp(9,14,2)=-0.47809144373375745599_dp
2494 1157 : ang_phipphj_tmp(9,16,2)=0.23145502494313786539_dp
2495 1157 : ang_phipphj_tmp(10,5,2)=0.15430334996209191026_dp
2496 1157 : ang_phipphj_tmp(11,6,2)=-0.094491118252306806804_dp
2497 1157 : ang_phipphj_tmp(12,5,2)=0.23904572186687872799_dp
2498 1157 : ang_phipphj_tmp(13,8,2)=0.58554004376911990761_dp
2499 1157 : ang_phipphj_tmp(14,7,2)=-0.82807867121082506136_dp
2500 1157 : ang_phipphj_tmp(14,9,2)=0.23904572186687872799_dp
2501 1157 : ang_phipphj_tmp(15,8,2)=-0.094491118252306806804_dp
2502 1157 : ang_phipphj_tmp(16,9,2)=0.15430334996209191026_dp
2503 : !
2504 1157 : ang_phipphj_tmp(1,4,3)=0.86602540378443864676_dp
2505 1157 : ang_phipphj_tmp(1,14,3)=0.54006172486732168591_dp
2506 1157 : ang_phipphj_tmp(2,5,3)=1.1180339887498948482_dp
2507 1157 : ang_phipphj_tmp(3,8,3)=1.1180339887498948482_dp
2508 1157 : ang_phipphj_tmp(4,9,3)=1.1180339887498948482_dp
2509 1157 : ang_phipphj_tmp(5,2,3)=-0.55901699437494742410_dp
2510 1157 : ang_phipphj_tmp(5,10,3)=1.6201851746019650577_dp
2511 1157 : ang_phipphj_tmp(6,11,3)=1.3228756555322952953_dp
2512 1157 : ang_phipphj_tmp(7,14,3)=1.2076147288491198811_dp
2513 1157 : ang_phipphj_tmp(8,15,3)=1.3228756555322952953_dp
2514 1157 : ang_phipphj_tmp(9,4,3)=-0.55901699437494742410_dp
2515 1157 : ang_phipphj_tmp(9,16,3)=1.6201851746019650577_dp
2516 1157 : ang_phipphj_tmp(10,5,3)=-1.0801234497346433718_dp
2517 1157 : ang_phipphj_tmp(11,6,3)=-0.66143782776614764763_dp
2518 1157 : ang_phipphj_tmp(15,8,3)=-0.66143782776614764763_dp
2519 1157 : ang_phipphj_tmp(16,9,3)=-1.0801234497346433718_dp
2520 : !
2521 1157 : ang_phipphj_tmp(1,2,4)=0.57735026918962576451_dp
2522 1157 : ang_phipphj_tmp(2,1,4)=0.57735026918962576451_dp
2523 1157 : ang_phipphj_tmp(2,7,4)=-0.25819888974716112568_dp
2524 1157 : ang_phipphj_tmp(2,9,4)=-0.44721359549995793928_dp
2525 1157 : ang_phipphj_tmp(3,6,4)=0.44721359549995793928_dp
2526 1157 : ang_phipphj_tmp(4,5,4)=0.44721359549995793928_dp
2527 1157 : ang_phipphj_tmp(5,4,4)=0.44721359549995793928_dp
2528 1157 : ang_phipphj_tmp(5,14,4)=-0.11952286093343936400_dp
2529 1157 : ang_phipphj_tmp(5,16,4)=-0.46291004988627573078_dp
2530 1157 : ang_phipphj_tmp(6,3,4)=0.44721359549995793928_dp
2531 1157 : ang_phipphj_tmp(6,13,4)=-0.29277002188455995381_dp
2532 1157 : ang_phipphj_tmp(6,15,4)=-0.37796447300922722721_dp
2533 1157 : ang_phipphj_tmp(7,2,4)=-0.25819888974716112568_dp
2534 1157 : ang_phipphj_tmp(7,12,4)=0.41403933560541253068_dp
2535 1157 : ang_phipphj_tmp(8,11,4)=0.37796447300922722721_dp
2536 1157 : ang_phipphj_tmp(9,2,4)=-0.44721359549995793928_dp
2537 1157 : ang_phipphj_tmp(9,10,4)=0.46291004988627573078_dp
2538 1157 : ang_phipphj_tmp(9,12,4)=0.11952286093343936400_dp
2539 1157 : ang_phipphj_tmp(10,9,4)=0.46291004988627573078_dp
2540 1157 : ang_phipphj_tmp(11,8,4)=0.37796447300922722721_dp
2541 1157 : ang_phipphj_tmp(12,7,4)=0.41403933560541253068_dp
2542 1157 : ang_phipphj_tmp(12,9,4)=0.11952286093343936400_dp
2543 1157 : ang_phipphj_tmp(13,6,4)=-0.29277002188455995381_dp
2544 1157 : ang_phipphj_tmp(14,5,4)=-0.11952286093343936400_dp
2545 1157 : ang_phipphj_tmp(15,6,4)=-0.37796447300922722721_dp
2546 1157 : ang_phipphj_tmp(16,5,4)=-0.46291004988627573078_dp
2547 : !
2548 1157 : ang_phipphj_tmp(1,2,5)=0.28867513459481288225_dp
2549 1157 : ang_phipphj_tmp(1,12,5)=-0.54006172486732168591_dp
2550 1157 : ang_phipphj_tmp(2,7,5)=-0.77459666924148337704_dp
2551 1157 : ang_phipphj_tmp(2,9,5)=-0.22360679774997896964_dp
2552 1157 : ang_phipphj_tmp(3,6,5)=0.22360679774997896964_dp
2553 1157 : ang_phipphj_tmp(4,5,5)=0.22360679774997896964_dp
2554 1157 : ang_phipphj_tmp(5,4,5)=0.11180339887498948482_dp
2555 1157 : ang_phipphj_tmp(5,14,5)=-0.47809144373375745599_dp
2556 1157 : ang_phipphj_tmp(5,16,5)=-0.23145502494313786539_dp
2557 1157 : ang_phipphj_tmp(6,3,5)=-0.44721359549995793928_dp
2558 1157 : ang_phipphj_tmp(6,13,5)=-1.1710800875382398152_dp
2559 1157 : ang_phipphj_tmp(6,15,5)=-0.18898223650461361361_dp
2560 1157 : ang_phipphj_tmp(7,2,5)=0.25819888974716112568_dp
2561 1157 : ang_phipphj_tmp(7,12,5)=0.44854261357253024157_dp
2562 1157 : ang_phipphj_tmp(8,11,5)=0.18898223650461361361_dp
2563 1157 : ang_phipphj_tmp(9,2,5)=-0.11180339887498948482_dp
2564 1157 : ang_phipphj_tmp(9,10,5)=0.23145502494313786539_dp
2565 1157 : ang_phipphj_tmp(9,12,5)=0.47809144373375745599_dp
2566 1157 : ang_phipphj_tmp(10,9,5)=0.15430334996209191026_dp
2567 1157 : ang_phipphj_tmp(11,8,5)=-0.094491118252306806804_dp
2568 1157 : ang_phipphj_tmp(12,7,5)=-0.82807867121082506136_dp
2569 1157 : ang_phipphj_tmp(12,9,5)=-0.23904572186687872799_dp
2570 1157 : ang_phipphj_tmp(13,6,5)=0.58554004376911990761_dp
2571 1157 : ang_phipphj_tmp(14,5,5)=0.23904572186687872799_dp
2572 1157 : ang_phipphj_tmp(15,6,5)=0.094491118252306806804_dp
2573 1157 : ang_phipphj_tmp(16,5,5)=-0.15430334996209191026_dp
2574 : !
2575 1157 : ang_phipphj_tmp(1,2,6)=0.86602540378443864676_dp
2576 1157 : ang_phipphj_tmp(1,12,6)=0.54006172486732168591_dp
2577 1157 : ang_phipphj_tmp(2,9,6)=-1.1180339887498948482_dp
2578 1157 : ang_phipphj_tmp(3,6,6)=1.1180339887498948482_dp
2579 1157 : ang_phipphj_tmp(4,5,6)=1.1180339887498948482_dp
2580 1157 : ang_phipphj_tmp(5,4,6)=-0.55901699437494742410_dp
2581 1157 : ang_phipphj_tmp(5,16,6)=-1.6201851746019650577_dp
2582 1157 : ang_phipphj_tmp(6,15,6)=-1.3228756555322952953_dp
2583 1157 : ang_phipphj_tmp(7,12,6)=1.2076147288491198811_dp
2584 1157 : ang_phipphj_tmp(8,11,6)=1.3228756555322952953_dp
2585 1157 : ang_phipphj_tmp(9,2,6)=0.55901699437494742410_dp
2586 1157 : ang_phipphj_tmp(9,10,6)=1.6201851746019650577_dp
2587 1157 : ang_phipphj_tmp(10,9,6)=-1.0801234497346433718_dp
2588 1157 : ang_phipphj_tmp(11,8,6)=-0.66143782776614764763_dp
2589 1157 : ang_phipphj_tmp(15,6,6)=0.66143782776614764763_dp
2590 1157 : ang_phipphj_tmp(16,5,6)=1.0801234497346433718_dp
2591 : !
2592 1157 : ang_phipphj_tmp(1,3,7)=0.57735026918962576451_dp
2593 1157 : ang_phipphj_tmp(2,6,7)=0.44721359549995793928_dp
2594 1157 : ang_phipphj_tmp(3,1,7)=0.57735026918962576451_dp
2595 1157 : ang_phipphj_tmp(3,7,7)=0.51639777949432225136_dp
2596 1157 : ang_phipphj_tmp(4,8,7)=0.44721359549995793928_dp
2597 1157 : ang_phipphj_tmp(5,11,7)=0.37796447300922722721_dp
2598 1157 : ang_phipphj_tmp(6,2,7)=0.44721359549995793928_dp
2599 1157 : ang_phipphj_tmp(6,12,7)=0.47809144373375745599_dp
2600 1157 : ang_phipphj_tmp(7,3,7)=0.51639777949432225136_dp
2601 1157 : ang_phipphj_tmp(7,13,7)=0.50709255283710994651_dp
2602 1157 : ang_phipphj_tmp(8,4,7)=0.44721359549995793928_dp
2603 1157 : ang_phipphj_tmp(8,14,7)=0.47809144373375745599_dp
2604 1157 : ang_phipphj_tmp(9,15,7)=0.37796447300922722721_dp
2605 1157 : ang_phipphj_tmp(11,5,7)=0.37796447300922722721_dp
2606 1157 : ang_phipphj_tmp(12,6,7)=0.47809144373375745599_dp
2607 1157 : ang_phipphj_tmp(13,7,7)=0.50709255283710994651_dp
2608 1157 : ang_phipphj_tmp(14,8,7)=0.47809144373375745599_dp
2609 1157 : ang_phipphj_tmp(15,9,7)=0.37796447300922722721_dp
2610 : !
2611 1157 : ang_phipphj_tmp(1,3,8)=1.1547005383792515290_dp
2612 1157 : ang_phipphj_tmp(2,6,8)=1.3416407864998738178_dp
2613 1157 : ang_phipphj_tmp(3,7,8)=1.5491933384829667541_dp
2614 1157 : ang_phipphj_tmp(4,8,8)=1.3416407864998738178_dp
2615 1157 : ang_phipphj_tmp(5,11,8)=1.5118578920369089089_dp
2616 1157 : ang_phipphj_tmp(6,2,8)=-0.44721359549995793928_dp
2617 1157 : ang_phipphj_tmp(6,12,8)=1.9123657749350298240_dp
2618 1157 : ang_phipphj_tmp(7,3,8)=-0.51639777949432225136_dp
2619 1157 : ang_phipphj_tmp(7,13,8)=2.0283702113484397860_dp
2620 1157 : ang_phipphj_tmp(8,4,8)=-0.44721359549995793928_dp
2621 1157 : ang_phipphj_tmp(8,14,8)=1.9123657749350298240_dp
2622 1157 : ang_phipphj_tmp(9,15,8)=1.5118578920369089089_dp
2623 1157 : ang_phipphj_tmp(11,5,8)=-0.75592894601845445443_dp
2624 1157 : ang_phipphj_tmp(12,6,8)=-0.95618288746751491198_dp
2625 1157 : ang_phipphj_tmp(13,7,8)=-1.0141851056742198930_dp
2626 1157 : ang_phipphj_tmp(14,8,8)=-0.95618288746751491198_dp
2627 1157 : ang_phipphj_tmp(15,9,8)=-0.75592894601845445443_dp
2628 :
2629 197773 : ang_phipphj(:,:,:)=ang_phipphj_tmp(1:mpsang**2,1:mpsang**2,:)
2630 :
2631 1157 : 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 311951 : 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 311951 : gaunt=zero;sum=zero;ok =.true.
2669 :
2670 311951 : if((-m1-mm+m2) /= 0) ok = .false.
2671 311951 : if(abs(m1) > l1) ok = .false.
2672 311951 : if(abs(mm) > ll) ok = .false.
2673 311951 : if(abs(m2) > l2) ok = .false.
2674 :
2675 311951 : jj = l1 + ll + l2
2676 311951 : if (mod(jj,2)/=0) ok = .false.
2677 311951 : j1 = jj-2*l2
2678 311951 : j2 = jj-2*ll
2679 311951 : j3 = jj-2*l1
2680 :
2681 311951 : if (j1<0 .or. j2<0 .or. j3<0) ok = .false.
2682 :
2683 311951 : if (ok) then
2684 :
2685 311951 : xx = (2 * l1 + 1) * (2 * ll + 1) * (2 * l2 + 1)
2686 :
2687 311951 : j1half = j1/2
2688 311951 : j2half = j2/2
2689 311951 : j3half = j3/2
2690 311951 : j_half = jj/2
2691 :
2692 311951 : gaunt = (-1)**j1half * sqrt(xx)
2693 1247804 : gaunt = gaunt * rfactorial(j2)*rfactorial(j3)/rfactorial(jj+1)
2694 : gaunt = gaunt * rfactorial(j_half)/(rfactorial(j1half)&
2695 1290584 : & * rfactorial(j2half)*rfactorial(j3half))
2696 :
2697 1439254 : yy = rfactorial(l2 + m2) * rfactorial(l2 - m2)
2698 :
2699 311951 : if (mm>=0) then
2700 192479 : yy = yy * perms(ll+mm,2*mm)
2701 : else
2702 119472 : yy = yy / perms(ll-mm,-2*mm)
2703 : end if
2704 :
2705 311951 : if (m1>=0) then
2706 188617 : yy = yy / perms(l1+m1,2*m1)
2707 : else
2708 123334 : yy = yy * perms(l1-m1,-2*m1)
2709 : end if
2710 :
2711 311951 : gaunt = gaunt * sqrt(yy)
2712 :
2713 311951 : i1 = l2 - ll - m1
2714 311951 : i2 = l2 - l1 + mm
2715 311951 : k1 = -min(0, i1, i2)
2716 311951 : n1 = l1 + m1
2717 311951 : n2 = ll - mm
2718 311951 : k2 = min(j1, n1, n2)
2719 :
2720 311951 : sign = 1._dp
2721 311951 : if(k1>0) sign = (-1._dp)**k1
2722 :
2723 623902 : argument = sign * perms(n1,k1)/rfactorial(k1)
2724 623902 : argument = argument * perms(n2,k1)/rfactorial(i1 + k1)
2725 623902 : argument = argument * perms(j1,k1)/rfactorial(i2 + k1)
2726 311951 : sum = sum + argument
2727 :
2728 311951 : sign = -sign
2729 311951 : k1 = k1 + 1
2730 354549 : do while(k1 <= k2)
2731 85196 : argument = sign * perms(n1, k1)/rfactorial(k1)
2732 85196 : argument = argument * perms(n2, k1)/rfactorial(i1 + k1)
2733 85196 : argument = argument * perms(j1, k1)/rfactorial(i2 + k1)
2734 42598 : sum = sum + argument
2735 42598 : sign = -sign
2736 42598 : k1 = k1 + 1
2737 : end do
2738 :
2739 : end if
2740 :
2741 311951 : gaunt = gaunt * sum
2742 :
2743 311951 : 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 2296 : 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 2296 : integer,allocatable :: ssgn(:)
2794 2296 : type(coeff3_type), allocatable :: coeff(:)
2795 :
2796 : !************************************************************************
2797 :
2798 6888 : 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 6443557 : gntselect = 0; realgnt = zero
2806 :
2807 : !Compute matrix cc where Sl=cc*Yl (Sl=real sph. harm.)
2808 : !------------------------------------------------
2809 21532 : LIBPAW_DATATYPE_ALLOCATE(coeff,(4*l_max-3))
2810 16940 : do ll=1,4*l_max-3
2811 58576 : LIBPAW_ALLOCATE(coeff(ll)%value,(2,2*ll-1,2*ll-1))
2812 4186928 : coeff(ll)%value(:,:,:)=zero
2813 14644 : coeff(ll)%value(1,ll,ll)=one
2814 65138 : do mm=1,ll-1
2815 48198 : coeff(ll)%value(1,ll+mm,ll+mm)= (-1._dp)**mm/sqrt(2._dp)
2816 48198 : coeff(ll)%value(1,ll-mm,ll+mm)= ( 1._dp) /sqrt(2._dp)
2817 48198 : coeff(ll)%value(2,ll+mm,ll-mm)=-(-1._dp)**mm/sqrt(2._dp)
2818 62842 : coeff(ll)%value(2,ll-mm,ll-mm)= ( 1._dp) /sqrt(2._dp)
2819 : end do
2820 : end do
2821 :
2822 6888 : LIBPAW_ALLOCATE(ssgn,(l_max**2))
2823 16019 : ssgn(:)=1
2824 2296 : if (l_max>0) then
2825 5383 : do l1=1,l_max-1
2826 3087 : ilm1=1+l1**2+l1
2827 9553 : do m1=-l1,-1
2828 7257 : ssgn(ilm1+m1)=-1
2829 : end do
2830 : end do
2831 : end if
2832 :
2833 2296 : ngnt=0
2834 :
2835 : !Loop on (lp1,mp1)
2836 : !------------------------------------------------
2837 7679 : do lp1=0,l_max-1
2838 21402 : do mp1=-lp1,lp1
2839 13723 : ilmp1=1+lp1**2+lp1+mp1
2840 13723 : k0lm1=ilmp1*(ilmp1-1)/2
2841 :
2842 : ! Loop on (l1,m1)<=(lp1,mp1)
2843 : ! ------------------------------------------------
2844 57789 : do l1=0,l_max-1
2845 173241 : do m1=-l1,l1
2846 120835 : ilm1=1+l1**2+l1+m1
2847 :
2848 159518 : if (ilm1<=ilmp1) then
2849 :
2850 67279 : klm1=k0lm1+ilm1
2851 2492222 : gntselect(:,klm1)=0
2852 :
2853 : ! Loop on (l2,m2)
2854 : ! ------------------------------------------------
2855 223334 : do l2=abs(l1-lp1),l1+lp1,2
2856 1095873 : do m2=-l2,l2
2857 872539 : ilm2=1+l2**2+l2+m2
2858 :
2859 : ! Real Gaunt coeffs selection rules
2860 : ! ------------------------------------------------
2861 872539 : 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 1028594 : ))))) then
2868 :
2869 : ! Compute selected real Gaunt coefficient
2870 : ! ------------------------------------------------
2871 163527 : realgnt_tmp=zero
2872 893126 : do mm1=-l1,l1
2873 729599 : c11=coeff(l1+1)%value(1,l1+mm1+1,l1+m1+1)
2874 729599 : c12=coeff(l1+1)%value(2,l1+mm1+1,l1+m1+1)
2875 5439909 : do mm2= -lp1,lp1
2876 4546783 : c21=coeff(lp1+1)%value(1,lp1+mm2+1,lp1+mp1+1)
2877 4546783 : c22=coeff(lp1+1)%value(2,lp1+mm2+1,lp1+mp1+1)
2878 42549445 : do mm3= -l2,l2
2879 37273063 : c31=coeff(l2+1)%value(1,l2+mm3+1,l2+m2+1)
2880 37273063 : c32=coeff(l2+1)%value(2,l2+mm3+1,l2+m2+1)
2881 : fact=c11*c21*c31 - c12*c22*c31&
2882 37273063 : & -c11*c22*c32 - c12*c21*c32
2883 37273063 : if((abs(fact)>=tol12).and.(mm3==-mm2-mm1)) &
2884 : & realgnt_tmp=realgnt_tmp+fact*(-1)**mm2 &
2885 4858734 : & *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 163527 : if (abs(realgnt_tmp)>=tol12) then
2893 162449 : ngnt=ngnt+1
2894 162449 : gntselect(ilm2,klm1)=ngnt
2895 162449 : 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 16940 : do ll=1,4*l_max-3
2912 16940 : LIBPAW_DEALLOCATE(coeff(ll)%value)
2913 : end do
2914 16940 : LIBPAW_DATATYPE_DEALLOCATE(coeff)
2915 2296 : LIBPAW_DEALLOCATE(ssgn)
2916 :
2917 2296 : 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 616 : 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 87976 : dyadic = zero
2970 7336 : do klm = 1, lm2_size
2971 67816 : do dlm = 1, 9
2972 60480 : if ( (dlm>1) .AND. (dlm<5) ) cycle ! no L=1 contributions
2973 40320 : ignt = gntselect(dlm,klm)
2974 40320 : if (ignt<1) cycle
2975 7808 : rgnt=realgnt(ignt)
2976 2544 : select case(dlm)
2977 :
2978 : ! S00, appears in xx, yy, zz
2979 : case (1)
2980 2544 : afact = rgnt*c3*(alpha-beta/three)
2981 2544 : dyadic(1,1,klm) = dyadic(1,1,klm)+afact
2982 2544 : dyadic(2,2,klm) = dyadic(2,2,klm)+afact
2983 2544 : dyadic(3,3,klm) = dyadic(3,3,klm)+afact
2984 :
2985 : ! S2,-2 appears in xy
2986 : case (5)
2987 664 : dyadic(1,2,klm) = dyadic(1,2,klm)-rgnt*beta*c1*c2
2988 :
2989 : ! S2,-1 appears yz
2990 : case (6)
2991 680 : 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 1944 : afact = rgnt*beta*c1/three
2996 1944 : dyadic(1,1,klm) = dyadic(1,1,klm) + afact
2997 1944 : dyadic(2,2,klm) = dyadic(2,2,klm) + afact
2998 1944 : dyadic(3,3,klm) = dyadic(3,3,klm) - two*afact
2999 :
3000 : ! S2,1 appears in xz
3001 : case (8)
3002 680 : dyadic(1,3,klm) = dyadic(1,3,klm)-rgnt*beta*c1*c2
3003 :
3004 : ! S2,2 appears in xx, yy
3005 : case (9)
3006 1296 : afact=rgnt*beta*c1*c2
3007 1296 : dyadic(1,1,klm) = dyadic(1,1,klm)-afact
3008 7808 : dyadic(2,2,klm) = dyadic(2,2,klm)+afact
3009 : case default
3010 : end select
3011 :
3012 7808 : dyadic(3,1,klm) = dyadic(1,3,klm)
3013 7808 : dyadic(2,1,klm) = dyadic(1,2,klm)
3014 67200 : dyadic(3,2,klm) = dyadic(2,3,klm)
3015 :
3016 : end do
3017 : end do
3018 :
3019 616 : 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 27409 : 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 41841 : 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 41841 : cosbeta=real(isn)*rot(3,3)
3426 41841 : sinbeta2=rot(1,3)**2+rot(2,3)**2
3427 41841 : if(sinbeta2<tol8**2)then
3428 16774 : sinbeta=zero
3429 16774 : cosalp=isn*rot(1,1)/cosbeta
3430 16774 : sinalp=isn*rot(1,2)/cosbeta
3431 16774 : cosgam=one
3432 16774 : singam=zero
3433 : else
3434 25067 : sinbeta=sqrt(sinbeta2)
3435 25067 : cosalp=isn*rot(3,1)/sinbeta
3436 25067 : sinalp=isn*rot(3,2)/sinbeta
3437 25067 : cosgam=-isn*rot(1,3)/sinbeta
3438 25067 : singam=isn*rot(2,3)/sinbeta
3439 : end if
3440 : !
3441 :
3442 : ! Check matrix:
3443 41841 : ier=0
3444 41841 : check=cosalp*cosbeta*cosgam-sinalp*singam
3445 41841 : if (abs(check-isn*rot(1,1))>tol8) ier=ier+1
3446 41841 : check=sinalp*cosbeta*cosgam+cosalp*singam
3447 41841 : if (abs(check-isn*rot(1,2))>tol8) ier=ier+1
3448 41841 : check=-sinbeta*cosgam
3449 41841 : if (abs(check-isn*rot(1,3))>tol8) ier=ier+1
3450 41841 : check=-cosalp*cosbeta*singam-sinalp*cosgam
3451 41841 : if (abs(check-isn*rot(2,1))>tol8) ier=ier+1
3452 41841 : check=-sinalp*cosbeta*singam+cosalp*cosgam
3453 41841 : if (abs(check-isn*rot(2,2))>tol8) ier=ier+1
3454 41841 : check=sinbeta*singam
3455 41841 : if (abs(check-isn*rot(2,3))>tol8) ier=ier+1
3456 41841 : check=cosalp*sinbeta
3457 41841 : if (abs(check-isn*rot(3,1))>tol8) ier=ier+1
3458 41841 : check=sinalp*sinbeta
3459 41841 : if (abs(check-isn*rot(3,2))>tol8) ier=ier+1
3460 41841 : 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 2372084 : 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 2372084 : if (mm==1) phim=sqrt2*costheta
3511 2372084 : if (mm==-1) phim=sqrt2*sintheta
3512 2372084 : if (mm==2) phim=sqrt2*(costheta*costheta-sintheta*sintheta)
3513 2372084 : if (mm==-2) phim=sqrt2*two*sintheta*costheta
3514 2372084 : if (mm==3) phim=sqrt2*&
3515 : & (costheta*(costheta*costheta-sintheta*sintheta)&
3516 27576 : & -sintheta*two*sintheta*costheta)
3517 2372084 : if (mm==-3) phim=sqrt2*&
3518 : & (sintheta*(costheta*costheta-sintheta*sintheta)&
3519 27576 : & +costheta*two*sintheta*costheta)
3520 :
3521 2372084 : 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 2372084 : 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 2372084 : dbeta=zero
3572 :
3573 : !Special cases
3574 2372084 : if (abs(cosbeta-1._dp).lt.tol10) then
3575 454140 : if (mp.eq.mm) dbeta=1
3576 1917944 : else if (abs(cosbeta+1._dp).lt.tol10) then
3577 419840 : 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 1498104 : if(cosbeta>-tol8)then
3588 : !If cosbeta is positive, cosbeta2 is positive with value >0.7, so one can divide by cosbetab2
3589 1476344 : cosbetab2=sqrt((1+cosbeta)*half)
3590 1476344 : 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 1498104 : ml=max(mp,mm)
3599 1498104 : ms=min(mp,mm)
3600 1498104 : if (ml.ne.mp) sinbetab2=-sinbetab2
3601 1498104 : 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 11972568 : & *((-sinbetab2)**(ml-ms))
3606 1498104 : sum=1._dp
3607 1498104 : arg=1._dp
3608 1498104 : ina=ml-ll
3609 1498104 : inb=-ms-ll
3610 1498104 : inc=ml-ms+1
3611 1966752 : do ii=1,mxterms
3612 1966752 : if (ina.eq.0.or.inb.eq.0) exit
3613 468648 : arg=(arg*ina*inb*tt)/(ii*inc)
3614 468648 : sum=sum+arg
3615 468648 : ina=ina+1
3616 468648 : inb=inb+1
3617 1966752 : inc=inc+1
3618 : end do
3619 1498104 : dbeta=pref*sum
3620 : end if
3621 :
3622 2372084 : 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 1687549 : 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 1687549 : if (nn>=0.and.nn>=kk) then
3725 1687549 : pp=1._dp
3726 4644559 : do ii=nn-kk+1,nn
3727 4644559 : pp=pp*ii
3728 : end do
3729 : else
3730 : pp=0._dp
3731 : end if
3732 :
3733 1687549 : perms=pp
3734 :
3735 1687549 : 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 7741616 : 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 7741616 : rfactorial=one
3770 24280433 : do ii=2,nn
3771 24280433 : rfactorial=rfactorial*ii
3772 : end do
3773 :
3774 : end function rfactorial
3775 : !!***
3776 :
3777 : END MODULE m_paw_sphharm
3778 : !!***
|
