Line data Source code
1 : !!****m* ABINIT/m_pawrad
2 : !! NAME
3 : !! m_pawrad
4 : !!
5 : !! FUNCTION
6 : !! Module containing all the functions related to the PAW radial meshes
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2013-2026 ABINIT group (MT,FJ,MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! NOTES
15 : !! * Routines tagged with "@type_name" are strongly connected to the definition of the data type.
16 : !! Strongly connected means that the proper functioning of the implementation relies on the
17 : !! assumption that the tagged procedure is consistent with the type declaration.
18 : !! Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
19 : !! that all the strongly connected routines are changed accordingly to accommodate the modification of the data type
20 : !! Typical examples of strongly connected routines are creation, destruction or reset methods.
21 : !!
22 : !! * FOR DEVELOPERS: in order to preserve the portability of libPAW library,
23 : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
24 : !!
25 : !! SOURCE
26 :
27 : #include "libpaw.h"
28 :
29 : MODULE m_pawrad
30 :
31 : USE_DEFS
32 : USE_MSG_HANDLING
33 : USE_MPI_WRAPPERS
34 : USE_MEMORY_PROFILING
35 :
36 : use m_paw_numeric, only : paw_derfc
37 :
38 : implicit none
39 :
40 : private
41 :
42 : !public procedures.
43 : public :: pawrad_init ! Main creation method
44 : public :: pawrad_free ! Free the allocated memory
45 : public :: pawrad_print ! Printout of the basic info
46 : public :: pawrad_isame ! Checks whether two meshes are equivalent or have the same equation.
47 : public :: pawrad_copy ! Returns a copy of the mesh.
48 : public :: pawrad_ifromr ! Retrieve the Index FROM a given R value in a radial grid.
49 : public :: pawrad_deducer0 ! Extrapolate r=0 value of a function from values near r=0.
50 : public :: pawrad_bcast ! Broadcast pawrad datastructure over a given MPI communicator
51 : public :: simp_gen ! Performs integral on a given (generalized) grid using Simpson rule.
52 : public :: nderiv_gen ! Do corrected first (and 2nd) derivation on a given (generalized) grid.
53 : public :: nderiv_lin ! Do corrected first (and 2nd) derivation on a given linear grid.
54 : public :: bound_deriv ! Computes derivatives at boundaries of the mesh
55 : public :: poisson ! Solves Poisson eq. for angularly dependent charge distribution of angular momentum l
56 : public :: screened_coul_kernel ! Kernel used to compute short-range screened Coulomb integrals
57 : public :: calc_slatradl ! Calculates the radial part of Slater integrals.
58 :
59 : interface pawrad_free
60 : module procedure pawrad_free_0D
61 : module procedure pawrad_free_1D
62 : end interface pawrad_free
63 :
64 : ! TODO: Might use bit flags, but all radmesh stuff should be encapsulated here!
65 : integer,private,parameter :: RMESH_LINEAR = 1
66 : integer,private,parameter :: RMESH_LOG1 = 2
67 : integer,private,parameter :: RMESH_LOG2 = 3
68 : integer,private,parameter :: RMESH_LOG3 = 4
69 : integer,private,parameter :: RMESH_NL = 5
70 : !!***
71 :
72 : !-------------------------------------------------------------------------
73 :
74 : !!****t* m_pawrad/pawrad_type
75 : !! NAME
76 : !! pawrad_type
77 : !!
78 : !! FUNCTION
79 : !! For PAW, RADial mesh discretization and related data
80 : !!
81 : !! SOURCE
82 :
83 : type, public :: pawrad_type
84 :
85 : !Integer scalars
86 :
87 : integer :: int_meshsz=0
88 : ! Mesh size used in integrals computation
89 : ! Integrals will be computed up to r(int_meshsz)
90 :
91 : integer :: mesh_size=0
92 : ! Dimension of radial mesh
93 :
94 : integer :: mesh_type=-1
95 : ! Type of mesh
96 : ! 1=regular grid: r(i)=(i-1)*AA
97 : ! 2=logarithmic grid: r(i)=AA*(exp[BB*(i-1)]-1)
98 : ! 3=logarithmic grid: r(i>1)=AA*exp[BB*(i-1)] and r(1)=0
99 : ! 4=logarithmic grid: r(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
100 :
101 : !Real (real(dp)) scalars
102 :
103 : real(dp) :: lstep=zero
104 : ! Exponential step of the mesh (BB parameter above)
105 : ! Defined only if mesh type is logarithmic
106 :
107 : real(dp) :: rmax=zero
108 : ! Max. value of r = rad(mesh_size)
109 :
110 : real(dp) :: rstep=zero
111 : ! Radial step of the mesh (AA parameter above)
112 :
113 : real(dp) :: stepint=zero
114 : ! Radial step used to convert any function from the
115 : ! present grid onto a regular grid in order to
116 : ! integrate it using trapeze method
117 :
118 : !Real (real(dp)) arrays
119 :
120 : real(dp), allocatable :: rad(:)
121 : ! rad(mesh_size)
122 : ! Coordinates of all the points of the mesh
123 :
124 : real(dp), allocatable :: radfact(:)
125 : ! radfact(mesh_size)
126 : ! Factor used to compute radial integrals
127 : ! Before being integrated on the present mesh,
128 : ! any function is multiplied by this factor
129 :
130 : real(dp), allocatable :: simfact(:)
131 : ! simfact(mesh_size)
132 : ! Factor used to compute radial integrals by the a Simpson scheme
133 : ! Integral[f] = Sum_i [simfact(i)*f(i)]
134 :
135 : end type pawrad_type
136 : !!***
137 :
138 : CONTAINS
139 : !===========================================================
140 : !!***
141 :
142 : !!****f* m_pawrad/pawrad_init
143 : !! NAME
144 : !! pawrad_init
145 : !!
146 : !! FUNCTION
147 : !! Creation method for radial meshes.
148 : !! Compute all points (and related weights) of a radial mesh.
149 : !! Grid can be regular or logarithimc.
150 : !!
151 : !! INPUTS
152 : !! [mesh_size]=Dimension of the radial mesh
153 : !! If not present, take mesh%mesh_size
154 : !! [mesh_type]=Type of mesh
155 : !! If not present, take mesh%mesh_type
156 : !! [rstep]=Radial step of the mesh (AA parameter above)
157 : !! If not present, take mesh%rstep
158 : !! [lstep]=Exponential step of the mesh (BB parameter above)
159 : !! If not present, take mesh%lstep
160 : !! Needed only if mesh type is logarithmic.
161 : !! [r_for_intg]=Mesh size used in integrals computation
162 : !! If not present, take mesh%r_for_intg
163 : !! Integrals will be computed up to rr(r_for_intg)
164 : !! (can be negative for an integration over the whole grid)
165 : !!
166 : !! OUTPUT
167 : !!
168 : !! SIDE EFFECTS
169 : !! mesh<pawrad_type>=The object completely initialized (containing radial grid information).
170 : !! The following quantities are calculated inside the routine:
171 : !! %stepint = Radial step used to convert any function from the
172 : !! present grid onto a regular grid in order to integrate it using trapeze method
173 : !! %rad(mesh_size) = Coordinates of all the points of the mesh.
174 : !! %radfact(mesh_size) = Factors used to compute radial integrals.
175 : !! %int_meshsz = Integrals will be computed up to r(int_meshsz)
176 : !! %simfact(mesh_size) = Factor used to compute radial integrals by the a Simpson scheme
177 : !! Integral[f] = Sum_i [simfact(i)*f(i)]
178 : !! %rmax = Max. value of r = rad(mesh_size)
179 : !!
180 : !! NOTES
181 : !! Possible mesh types (mesh%mesh_type)
182 : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
183 : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
184 : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
185 : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
186 : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
187 : !!
188 : !! SOURCE
189 :
190 5156 : subroutine pawrad_init(mesh,mesh_size,mesh_type,rstep,lstep,r_for_intg)
191 :
192 : !Arguments ------------------------------------
193 : !scalars
194 : integer,intent(in),optional :: mesh_size,mesh_type
195 : real(dp),intent(in),optional :: rstep,lstep
196 : real(dp),intent(in),optional :: r_for_intg
197 : type(pawrad_type),intent(inout) :: mesh
198 :
199 : !Local variables-------------------------------
200 : !scalars
201 : integer :: ir,ir_last,isim,mesh_size_,mesh_type_
202 : real(dp) :: hh,lstep_,rstep_,r_for_intg_
203 : character(len=500) :: msg
204 :
205 : ! *************************************************************************
206 :
207 : !@pawrad_type
208 :
209 : !Retrieve mesh data
210 5156 : mesh_size_ = mesh%mesh_size ; if (present(mesh_size)) mesh_size_ = mesh_size
211 5156 : mesh_type_ = mesh%mesh_type ; if (present(mesh_type)) mesh_type_ = mesh_type
212 5156 : rstep_ = mesh%rstep ; if (present(rstep)) rstep_ = rstep
213 5156 : lstep_ = mesh%lstep ; if (present(lstep)) lstep_ = lstep
214 :
215 5156 : r_for_intg_=-1._dp;if (present(r_for_intg)) r_for_intg_=r_for_intg
216 :
217 5156 : mesh%mesh_size = mesh_size_
218 5156 : mesh%mesh_type = mesh_type_
219 5156 : mesh%rstep = rstep_
220 5156 : mesh%lstep = lstep_
221 15468 : LIBPAW_ALLOCATE(mesh%rad ,(mesh%mesh_size))
222 10312 : LIBPAW_ALLOCATE(mesh%radfact,(mesh%mesh_size))
223 10312 : LIBPAW_ALLOCATE(mesh%simfact,(mesh%mesh_size))
224 7829511 : mesh%simfact=zero
225 : if (mesh%mesh_type==1) then
226 1636 : isim=3
227 1636 : mesh%stepint=mesh%rstep
228 1636 : mesh%rad(1)=zero;mesh%radfact(1)=one
229 3038160 : do ir=2,mesh%mesh_size
230 3036524 : mesh%rad(ir) =mesh%rstep*dble(ir-1)
231 3038160 : mesh%radfact(ir)=one
232 : end do
233 : else if (mesh%mesh_type==2) then
234 3371 : isim=3
235 3371 : mesh%stepint=mesh%lstep
236 3371 : mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep
237 4693128 : do ir=2,mesh%mesh_size
238 4689757 : mesh%rad(ir) =mesh%rstep*(exp(mesh%lstep*dble(ir-1))-one)
239 4693128 : mesh%radfact(ir)=mesh%rad(ir)+mesh%rstep
240 : end do
241 : else if (mesh%mesh_type==3) then
242 133 : isim=4
243 133 : mesh%stepint=mesh%lstep
244 133 : mesh%rad(1)=zero;mesh%radfact(1)=zero
245 85883 : do ir=2,mesh%mesh_size
246 85750 : mesh%rad(ir) =mesh%rstep*exp(mesh%lstep*dble(ir-2))
247 85883 : mesh%radfact(ir)=mesh%rad(ir)
248 : end do
249 : else if (mesh%mesh_type==4) then
250 0 : isim=3
251 0 : mesh%lstep=one/dble(mesh%mesh_size)
252 0 : mesh%stepint=mesh%lstep
253 0 : mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep
254 0 : do ir=2,mesh%mesh_size
255 0 : mesh%rad(ir) =-mesh%rstep*log(one-mesh%lstep*dble(ir-1))
256 0 : mesh%radfact(ir)=mesh%rstep/(one-mesh%lstep*dble(ir-1))
257 : end do
258 : else if (mesh%mesh_type==5) then
259 16 : isim=3
260 16 : mesh%stepint=mesh%rstep
261 16 : mesh%rad(1)=zero;mesh%radfact(1)=1/mesh%lstep
262 7184 : do ir=2,mesh%mesh_size
263 7168 : mesh%rad(ir) =mesh%rstep*dble(ir-1)/(mesh%lstep-dble(ir-1))
264 7184 : mesh%radfact(ir)=(mesh%rstep+mesh%rad(ir))/(mesh%lstep-dble(ir-1))/mesh%rstep
265 : end do
266 :
267 : else ! Other values of mesh_type are not allowed (see psp7in.F90)
268 0 : write(msg,'(a,i0)')" Unknown value of mesh_type: ",mesh%mesh_type
269 0 : LIBPAW_ERROR(msg)
270 : end if
271 :
272 5156 : mesh%int_meshsz=mesh%mesh_size
273 5156 : if (r_for_intg_>0.d0) then
274 2108 : ir=min(pawrad_ifromr(mesh,r_for_intg_),mesh%mesh_size)
275 2108 : if (ir<mesh%mesh_size) then
276 1786 : if (abs(mesh%rad(ir+1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir+1
277 : end if
278 2108 : if (ir>1) then
279 2108 : if (abs(mesh%rad(ir-1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir-1
280 : end if
281 2108 : mesh%int_meshsz=ir
282 : end if
283 :
284 5156 : hh=mesh%stepint/3.d0
285 5156 : mesh%simfact(mesh%int_meshsz)=hh*mesh%radfact(mesh%int_meshsz)
286 10445 : mesh%simfact(1:isim-2)=zero
287 5156 : ir_last=1
288 5156 : do ir=mesh%int_meshsz,isim,-2
289 3487353 : mesh%simfact(ir-1)=4.d0*hh*mesh%radfact(ir-1)
290 3487353 : mesh%simfact(ir-2)=2.d0*hh*mesh%radfact(ir-2)
291 3487353 : ir_last=ir-2
292 : end do
293 5156 : mesh%simfact(ir_last)=half*mesh%simfact(ir_last)
294 847601 : if (mesh%int_meshsz<mesh%mesh_size) mesh%simfact(mesh%int_meshsz+1:mesh%mesh_size)=zero
295 :
296 5156 : mesh%rmax=mesh%rad(mesh%mesh_size)
297 :
298 5156 : end subroutine pawrad_init
299 : !!***
300 :
301 : !----------------------------------------------------------------------
302 :
303 : !!****f* m_pawrad/pawrad_free_0D
304 : !! NAME
305 : !! pawrad_free_0D
306 : !!
307 : !! FUNCTION
308 : !! Frees all memory allocated in the object
309 : !!
310 : !! SOURCE
311 :
312 9633 : subroutine pawrad_free_0D(Rmesh)
313 :
314 : !Arguments ------------------------------------
315 : !arrays
316 : type(pawrad_type),intent(inout) :: Rmesh
317 :
318 : !Local variables-------------------------------
319 :
320 : ! *************************************************************************
321 :
322 : !@Pawrad_type
323 :
324 9633 : if (allocated(Rmesh%rad )) then
325 8690 : LIBPAW_DEALLOCATE(Rmesh%rad)
326 : end if
327 9633 : if (allocated(Rmesh%radfact)) then
328 8690 : LIBPAW_DEALLOCATE(Rmesh%radfact)
329 : end if
330 9633 : if (allocated(Rmesh%simfact)) then
331 8690 : LIBPAW_DEALLOCATE(Rmesh%simfact)
332 : end if
333 9633 : Rmesh%int_meshsz=0
334 9633 : Rmesh%mesh_size=0
335 9633 : Rmesh%mesh_type=-1
336 :
337 9633 : end subroutine pawrad_free_0D
338 : !!***
339 :
340 : !----------------------------------------------------------------------
341 :
342 : !!****f* m_pawrad/pawrad_free_1D
343 : !! NAME
344 : !! pawrad_free_1D
345 : !!
346 : !! FUNCTION
347 : !! Destroy all objects in an array of pawrad data structures
348 : !!
349 : !! SOURCE
350 :
351 1979 : subroutine pawrad_free_1D(Rmesh)
352 :
353 : !Arguments ------------------------------------
354 : type(pawrad_type),intent(inout) :: Rmesh(:)
355 :
356 : !Local variables-------------------------------
357 : integer :: ii,nn
358 :
359 : ! *************************************************************************
360 :
361 : !@pawrad_type
362 :
363 1979 : nn=size(Rmesh)
364 1979 : if (nn==0) return
365 :
366 2918 : do ii=1,nn
367 2918 : call pawrad_free_0D(Rmesh(ii))
368 : end do
369 :
370 : end subroutine pawrad_free_1D
371 : !!***
372 :
373 : !----------------------------------------------------------------------
374 :
375 : !!****f* m_pawrad/pawrad_print
376 : !! NAME
377 : !! pawrad_print
378 : !!
379 : !! FUNCTION
380 : !! Reports basic info on the object.
381 : !!
382 : !! INPUTS
383 : !! Rmesh<pawrad_type>=Object defining the radial mesh
384 : !! header=String for the header provided by the user.
385 : !! [unit]=Unit number for output, defaults to std_out
386 : !! [prtvol]=Verbosity level, minimal if not specified.
387 : !! [mode_paral]=Either "COLL" or "PERS". Passed to wrtout. Defaults to "COLL"
388 : !!
389 : !! OUTPUT
390 : !! Only writing.
391 : !!
392 : !! SOURCE
393 :
394 0 : subroutine pawrad_print(Rmesh,header,unit,prtvol,mode_paral)
395 :
396 : !Arguments ------------------------------------
397 : integer,intent(in),optional :: prtvol,unit
398 : character(len=4),intent(in),optional :: mode_paral
399 : character(len=*),intent(in),optional :: header
400 : type(pawrad_type),intent(in) :: Rmesh
401 :
402 : !Local variables-------------------------------
403 : !scalars
404 : integer :: my_unt,my_prtvol
405 : character(len=4) :: my_mode
406 : character(len=500) :: msg
407 :
408 : ! *************************************************************************
409 :
410 : !@pawrad_type
411 0 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
412 0 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
413 0 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
414 :
415 0 : msg=ch10//' ==== Info on the Radial Mesh ==== '
416 0 : if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== '
417 0 : call wrtout(my_unt,msg,my_mode)
418 :
419 0 : SELECT CASE (Rmesh%mesh_type)
420 :
421 : CASE (RMESH_LINEAR)
422 : write(msg,'(a,i4,a,g12.5)')&
423 0 : & ' - Linear mesh: r(i)=step*(i-1), size=',Rmesh%mesh_size,', step=',Rmesh%rstep
424 :
425 : CASE (RMESH_LOG1)
426 : write(msg,'(a,i4,2(a,g12.5))')&
427 0 : & ' - Logarithimc mesh: r(i)=AA*[exp(BB*(i-1))-1], size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep
428 :
429 : CASE (RMESH_LOG2)
430 : write(msg,'(a,i4,2(a,g12.5))')&
431 0 : & ' - Logarithimc mesh: r(i)=AA*exp(BB*(i-2)), size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep
432 :
433 : CASE (RMESH_LOG3)
434 : write(msg,'(a,i1,a,i4,a,g12.5)')&
435 0 : & ' - Logarithimc mesh: r(i)=-AA*ln(1-(i-1)/n), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep
436 :
437 : CASE (RMESH_NL)
438 : write(msg,'(a,i1,a,i4,a,g12.5)')&
439 0 : & ' - Non-linear mesh: r(i)=-AA*i/(n-i), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep
440 :
441 : CASE DEFAULT
442 0 : msg = ' Unknown mesh type! Action : check your pseudopotential or input file.'
443 0 : LIBPAW_ERROR(msg)
444 : END SELECT
445 :
446 0 : call wrtout(my_unt,msg,my_mode)
447 :
448 0 : if (my_prtvol>1) then
449 0 : write(msg,'(a,i4)')' Mesh size for integrals = ',Rmesh%int_meshsz
450 0 : call wrtout(my_unt,msg,my_mode)
451 0 : write(msg,'(a,g12.5)')' rmax=rad(mesh_size) = ',Rmesh%rmax
452 0 : call wrtout(my_unt,msg,my_mode)
453 0 : write(msg,'(a,g12.5)')' Value of stepint = ',Rmesh%stepint
454 0 : call wrtout(my_unt,msg,my_mode)
455 : end if
456 :
457 0 : end subroutine pawrad_print
458 : !!***
459 :
460 : !----------------------------------------------------------------------
461 :
462 : !!****f* m_pawrad/pawrad_isame
463 : !! NAME
464 : !! pawrad_isame
465 : !!
466 : !! FUNCTION
467 : !! Check two radial meshes, returns a logical flag defining
468 : !! whether the meshes have the same equation and an integer
469 : !! flag
470 : !!
471 : !! INPUTS
472 : !! Rmesh1,Rmesh2<pawrad_type>=The two radial meshes.
473 : !!
474 : !! OUTPUT
475 : !! hasameq=.true. if the two meshes are defined by the same equation.
476 : !! whichdenser=
477 : !! * 0 if meshes are not compatible
478 : !! * 1 if Rmesh1 is denser than Rmesh2
479 : !! * 2 if Rmesh2 is denser than Rmesh1
480 : !!
481 : !! SOURCE
482 :
483 2 : subroutine pawrad_isame(Rmesh1,Rmesh2,hasameq,whichdenser)
484 :
485 : !Arguments ------------------------------------
486 : integer,intent(out) :: whichdenser
487 : logical,intent(out) :: hasameq
488 : type(pawrad_type),intent(in) :: Rmesh1
489 : type(pawrad_type),intent(in) :: Rmesh2
490 :
491 : !Local variables-------------------------------
492 : character(len=50) :: msg
493 :
494 : ! *************************************************************************
495 :
496 : !@pawrad_type
497 :
498 2 : whichdenser =0 ; hasameq=.FALSE.
499 :
500 2 : if (Rmesh1%mesh_type /= Rmesh2%mesh_type) RETURN
501 :
502 0 : SELECT CASE (Rmesh1%mesh_type)
503 :
504 : CASE (RMESH_LINEAR) !check for linear meshes
505 0 : hasameq = (Rmesh1%rstep == Rmesh2%rstep)
506 :
507 : CASE (RMESH_LOG1,& !check for logarithmic meshes
508 : & RMESH_LOG2,&
509 : & RMESH_LOG3)
510 :
511 : hasameq = ( Rmesh1%rstep == Rmesh2%rstep &
512 4 : & .and.Rmesh1%lstep == Rmesh2%lstep )
513 :
514 : CASE (RMESH_NL) !check for linear meshes
515 0 : hasameq = (Rmesh1%rstep == Rmesh2%rstep)
516 :
517 : CASE DEFAULT
518 0 : msg='Unknown mesh type'
519 2 : LIBPAW_ERROR(msg)
520 :
521 : END SELECT
522 :
523 : ! === If meshes have same equation, check whether they are equal ===
524 : ! * Note that also int_meshsz must be equal
525 2 : if (hasameq) then
526 0 : whichdenser= 1
527 0 : if (Rmesh2%mesh_size > Rmesh1%mesh_size ) whichdenser = 2
528 : !if (Rmesh1%int_meshsz == Rmesh2%int_meshsz) whichdenser = 2
529 : end if
530 :
531 : end subroutine pawrad_isame
532 : !!***
533 :
534 : !----------------------------------------------------------------------
535 :
536 : !!****f* m_pawrad/pawrad_copy
537 : !! NAME
538 : !! pawrad_copy
539 : !!
540 : !! FUNCTION
541 : !! Copy one radial mesh (in a generalized format) to another
542 : !!
543 : !! INPUTS
544 : !! mesh1 <type(pawrad_type)>=data containing radial grid information of input mesh
545 : !!
546 : !! OUTPUT
547 : !! mesh2 <type(pawrad_type)>=data containing radial grid information of output mesh
548 : !!
549 : !! NOTES
550 : !! Possible mesh types (mesh%mesh_type)
551 : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
552 : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
553 : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
554 : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
555 : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
556 : !!
557 : !! SOURCE
558 :
559 3350 : subroutine pawrad_copy(mesh1,mesh2)
560 :
561 : !Arguments ------------------------------------
562 : !scalars
563 : type(pawrad_type),intent(in) :: mesh1
564 : type(pawrad_type),intent(out) :: mesh2
565 :
566 : !Local variables-------------------------------
567 : !scalars
568 : integer :: ir
569 :
570 : ! *************************************************************************
571 :
572 3350 : mesh2%mesh_type =mesh1%mesh_type
573 3350 : mesh2%mesh_size =mesh1%mesh_size
574 3350 : mesh2%int_meshsz=mesh1%int_meshsz
575 3350 : mesh2%lstep =mesh1%lstep
576 3350 : mesh2%rstep =mesh1%rstep
577 3350 : mesh2%stepint =mesh1%stepint
578 3350 : mesh2%rmax =mesh1%rmax
579 :
580 10050 : LIBPAW_ALLOCATE(mesh2%rad,(mesh1%mesh_size))
581 6700 : LIBPAW_ALLOCATE(mesh2%radfact,(mesh1%mesh_size))
582 6700 : LIBPAW_ALLOCATE(mesh2%simfact,(mesh1%mesh_size))
583 3741489 : do ir=1,mesh1%mesh_size
584 3738139 : mesh2%rad(ir) =mesh1%rad(ir)
585 3738139 : mesh2%radfact(ir)=mesh1%radfact(ir)
586 3741489 : mesh2%simfact(ir)=mesh1%simfact(ir)
587 : end do
588 :
589 3350 : end subroutine pawrad_copy
590 : !!***
591 :
592 : !----------------------------------------------------------------------
593 :
594 : !!****f* m_pawrad/pawrad_deducer0
595 : !! NAME
596 : !! pawrad_deducer0
597 : !!
598 : !! FUNCTION
599 : !! Extrapolate r=0 value of a function from values near r=0
600 : !! using a 3 points formula
601 : !!
602 : !! INPUTS
603 : !! funcsz=size of array func
604 : !! radmesh <type(pawrad_type)>=data containing radial grid information
605 : !!
606 : !! SIDE EFFECTS
607 : !! func(funcsz)=array containing values of function to extrapolate
608 : !!
609 : !! SOURCE
610 :
611 11374689 : subroutine pawrad_deducer0(func,funcsz,radmesh)
612 :
613 : !Arguments ------------------------------------
614 : !scalars
615 : integer,intent(in) :: funcsz
616 : type(pawrad_type),intent(in) :: radmesh
617 : !arrays
618 : real(dp),intent(inout) :: func(funcsz)
619 :
620 : ! *************************************************************************
621 :
622 11374689 : if (radmesh%mesh_type==1.or.radmesh%mesh_type==2.or.radmesh%mesh_type==4.or.radmesh%mesh_type==5) then
623 11265318 : func(1)=func(4)+3*(func(2)-func(3))
624 109371 : else if (radmesh%mesh_type==3) then
625 109371 : func(1)=func(4)+exp(two*radmesh%lstep)/(exp(radmesh%lstep)-one)*(func(2)-func(3))
626 : end if
627 :
628 11374689 : end subroutine pawrad_deducer0
629 : !!***
630 :
631 : !----------------------------------------------------------------------
632 :
633 : !!****f* m_pawrad/pawrad_bcast
634 : !! NAME
635 : !! pawrad_bcast
636 : !!
637 : !! FUNCTION
638 : !! Communicate pawrad data over a given MPI communicator
639 : !!
640 : !! INPUTS
641 : !! comm_mpi= communicator used to broadcast data
642 : !!
643 : !! SIDE EFFECTS
644 : !! pawrad=<type pawrad_type>= a radial mesh datastructure for PAW
645 : !!
646 : !! SOURCE
647 :
648 276 : subroutine pawrad_bcast(pawrad,comm_mpi)
649 :
650 : !Arguments ------------------------------------
651 : !scalars
652 : integer,intent(in) :: comm_mpi
653 : type(pawrad_type),intent(inout) :: pawrad
654 :
655 : !Local variables-------------------------------
656 : !scalars
657 : integer :: ierr,indx,me,nn,isz1
658 : integer :: if_rad,if_radfact,if_simfact !flags used to communicate
659 : character(len=500) :: msg
660 : !arrays
661 276 : integer,allocatable :: list_int(:)
662 276 : real(dp),allocatable :: list_dpr(:)
663 :
664 : !*************************************************************************
665 :
666 276 : me=xmpi_comm_rank(comm_mpi)
667 :
668 : !Initializations
669 276 : if_rad=0; if_radfact=0; if_simfact=0
670 :
671 : !calculate the size of the reals
672 276 : if(me==0) then
673 70 : if (allocated(pawrad%rad)) then
674 70 : if_rad=1 !communicate rad
675 70 : isz1=size(pawrad%rad)
676 70 : if(isz1/=pawrad%mesh_size) then
677 0 : msg='rad: sz1 /= pawrad%mesh_size (1)'
678 0 : LIBPAW_BUG(msg)
679 : end if
680 : end if
681 70 : if (allocated(pawrad%radfact)) then
682 70 : if_radfact=1 !communicate radfact
683 70 : isz1=size(pawrad%radfact)
684 70 : if(isz1/=pawrad%mesh_size) then
685 0 : msg='radfact: sz1 /= pawrad%mesh_size (2)'
686 0 : LIBPAW_BUG(msg)
687 : end if
688 : end if
689 70 : if (allocated(pawrad%simfact)) then
690 70 : if_simfact=1 !communicate simfact
691 70 : isz1=size(pawrad%simfact)
692 70 : if(isz1/=pawrad%mesh_size) then
693 0 : msg='simfact: sz1 /= pawrad%mesh_size (3)'
694 0 : LIBPAW_BUG(msg)
695 : end if
696 : end if
697 : end if
698 :
699 : !Brodcast the integers
700 276 : LIBPAW_ALLOCATE(list_int,(6))
701 276 : if(me==0) then
702 70 : list_int(1)=pawrad%int_meshsz
703 70 : list_int(2)=pawrad%mesh_size
704 70 : list_int(3)=pawrad%mesh_type
705 70 : list_int(4)=if_rad
706 70 : list_int(5)=if_radfact
707 70 : list_int(6)=if_simfact
708 : end if
709 276 : call xmpi_bcast(list_int,0,comm_mpi,ierr)
710 276 : if(me/=0) then
711 206 : pawrad%int_meshsz =list_int(1)
712 206 : pawrad%mesh_size =list_int(2)
713 206 : pawrad%mesh_type =list_int(3)
714 206 : if_rad=list_int(4)
715 206 : if_radfact=list_int(5)
716 206 : if_simfact=list_int(6)
717 : end if
718 276 : LIBPAW_DEALLOCATE(list_int)
719 :
720 : !Broadcast the reals
721 276 : nn=4+pawrad%mesh_size*(if_rad+if_radfact+if_simfact)
722 828 : LIBPAW_ALLOCATE(list_dpr,(nn))
723 276 : if(me==0) then
724 70 : list_dpr(1)=pawrad%lstep
725 70 : list_dpr(2)=pawrad%rmax
726 70 : list_dpr(3)=pawrad%rstep
727 70 : list_dpr(4)=pawrad%stepint
728 70 : indx=5
729 70 : if (if_rad==1) then
730 70 : isz1=pawrad%mesh_size
731 79770 : list_dpr(indx:indx+isz1-1)=pawrad%rad(1:isz1)
732 70 : indx=indx+isz1
733 : end if
734 70 : if (if_radfact==1) then
735 70 : isz1=pawrad%mesh_size
736 79770 : list_dpr(indx:indx+isz1-1)=pawrad%radfact(1:isz1)
737 : indx=indx+isz1
738 : end if
739 70 : if (if_simfact==1) then
740 70 : isz1=pawrad%mesh_size
741 79770 : list_dpr(indx:indx+isz1-1)=pawrad%simfact(1:isz1)
742 276 : indx=indx+isz1
743 : end if
744 : end if
745 276 : call xmpi_bcast(list_dpr,0,comm_mpi,ierr)
746 276 : if(me/=0) then
747 206 : pawrad%lstep=list_dpr(1)
748 206 : pawrad%rmax=list_dpr(2)
749 206 : pawrad%rstep=list_dpr(3)
750 206 : pawrad%stepint=list_dpr(4)
751 206 : indx=5
752 : ! Deallocate all arrays:
753 206 : if (allocated(pawrad%rad)) then
754 21 : LIBPAW_DEALLOCATE(pawrad%rad)
755 : end if
756 206 : if (allocated(pawrad%radfact)) then
757 21 : LIBPAW_DEALLOCATE(pawrad%radfact)
758 : end if
759 206 : if (allocated(pawrad%simfact)) then
760 21 : LIBPAW_DEALLOCATE(pawrad%simfact)
761 : end if
762 : ! Communicate if flag is set to 1:
763 206 : if(if_rad==1) then
764 206 : isz1=pawrad%mesh_size
765 618 : LIBPAW_ALLOCATE(pawrad%rad,(isz1))
766 243438 : pawrad%rad(1:isz1)=list_dpr(indx:indx+isz1-1)
767 206 : indx=indx+isz1
768 : end if
769 206 : if(if_radfact==1) then
770 206 : isz1=pawrad%mesh_size
771 618 : LIBPAW_ALLOCATE(pawrad%radfact,(isz1))
772 243438 : pawrad%radfact(1:isz1)=list_dpr(indx:indx+isz1-1)
773 : indx=indx+isz1
774 : end if
775 206 : if(if_simfact==1) then
776 206 : isz1=pawrad%mesh_size
777 618 : LIBPAW_ALLOCATE(pawrad%simfact,(isz1))
778 243438 : pawrad%simfact(1:isz1)=list_dpr(indx:indx+isz1-1)
779 276 : indx=indx+isz1
780 : end if
781 : end if
782 276 : LIBPAW_DEALLOCATE(list_dpr)
783 :
784 276 : end subroutine pawrad_bcast
785 : !!***
786 :
787 : !----------------------------------------------------------------------
788 :
789 : !!****f* m_pawrad/simp_gen
790 : !! NAME
791 : !! simp_gen
792 : !!
793 : !! FUNCTION
794 : !! Do integral on a given (generalized) grid using Simpson rule
795 : !!
796 : !! INPUTS
797 : !! radmesh <type(pawrad_type)>=data containing radial grid information
798 : !! func(:)=integrand values
799 : !! r_for_intg=upper boundary for (future) integration over the radial grid
800 : !! (can be negative for an integration over the whole grid)
801 : !!
802 : !! OUTPUT
803 : !! intg=resulting integral by Simpson rule
804 : !!
805 : !! NOTES
806 : !! Possible mesh types (radmesh%mesh_type)
807 : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
808 : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
809 : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
810 : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
811 : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
812 : !!
813 : !! SOURCE
814 :
815 89191545 : subroutine simp_gen(intg,func,radmesh,r_for_intg)
816 :
817 : #if defined HAVE_AVX_SAFE_MODE
818 : !DEC$ NOOPTIMIZE
819 : #endif
820 :
821 : !Arguments ------------------------------------
822 : !scalars
823 : real(dp),intent(out) :: intg
824 : real(dp),intent(in),optional :: r_for_intg
825 : type(pawrad_type),intent(in) :: radmesh
826 : !arrays
827 : real(dp),intent(in) :: func(:)
828 :
829 : !Local variables-------------------------------
830 : !scalars
831 : integer :: ii,int_meshsz,ir,ir_last,isim,nn
832 : real(dp) :: hh,resid,simp
833 89191545 : real(dp),allocatable :: simfact(:)
834 : character(len=500) :: msg
835 :
836 : ! *************************************************************************
837 :
838 89191545 : if (present(r_for_intg)) then
839 3275 : if (r_for_intg>0.d0) then
840 3275 : ir=min(pawrad_ifromr(radmesh,r_for_intg),radmesh%mesh_size)
841 3275 : if (ir<radmesh%mesh_size) then
842 2650 : if (abs(radmesh%rad(ir+1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir+1
843 : end if
844 3275 : if (ir>1) then
845 3275 : if (abs(radmesh%rad(ir-1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir-1
846 : end if
847 3275 : int_meshsz=ir
848 : else
849 0 : int_meshsz=radmesh%mesh_size
850 : end if
851 3275 : if (int_meshsz>radmesh%mesh_size.or.int_meshsz>size(func)) then
852 0 : write(msg,'(3(a,i4))')"int_meshsz= ",int_meshsz," > mesh_size=",radmesh%mesh_size,&
853 0 : & ", size(func)=",size(func)
854 0 : LIBPAW_BUG(msg)
855 : end if
856 3275 : isim=3; if (radmesh%mesh_type==3)isim=4
857 9825 : LIBPAW_ALLOCATE(simfact,(radmesh%mesh_size))
858 3275 : hh=radmesh%stepint/3.d0
859 3275 : simfact(int_meshsz)=hh*radmesh%radfact(int_meshsz)
860 6743 : simfact(1:isim-2)=zero
861 3275 : ir_last=1
862 3275 : do ir=int_meshsz,isim,-2
863 1537138 : simfact(ir-1)=4.d0*hh*radmesh%radfact(ir-1)
864 1537138 : simfact(ir-2)=2.d0*hh*radmesh%radfact(ir-2)
865 1537138 : ir_last=ir-2
866 : end do
867 3275 : simfact(ir_last)=half*simfact(ir_last)
868 998441 : if (int_meshsz<radmesh%mesh_size) simfact(int_meshsz+1:radmesh%mesh_size)=zero
869 :
870 3275 : nn=int_meshsz
871 3275 : simp=zero
872 3082713 : do ii=1,nn
873 3082713 : simp=simp+func(ii)*simfact(ii)
874 : end do
875 3275 : LIBPAW_DEALLOCATE(simfact)
876 :
877 : else
878 89188270 : if (radmesh%int_meshsz>size(func)) then
879 0 : write(msg,'(2(a,i4))')"int_meshsz= ",int_meshsz," > size(func)=",size(func)
880 0 : LIBPAW_BUG(msg)
881 : end if
882 89188270 : nn=radmesh%int_meshsz
883 89188270 : simp=zero
884 >10995*10^7 : do ii=1,nn
885 >10995*10^7 : simp=simp+func(ii)*radmesh%simfact(ii)
886 : end do
887 : end if
888 :
889 89191545 : resid=zero
890 89191545 : if (radmesh%mesh_type==3) then
891 1302684 : resid=half*(func(2)+func(1))*(radmesh%rad(2)-radmesh%rad(1))
892 1302684 : if (mod(nn,2)==1) resid=resid+radmesh%stepint/3.d0*(1.25d0*func(2)*radmesh%radfact(2) &
893 99204 : & +2.d0*func(3)*radmesh%radfact(3)-0.25d0*func(4)*radmesh%radfact(4))
894 87888861 : else if (mod(nn,2)==0) then
895 : resid=radmesh%stepint/3.d0*(1.25d0*func(1)*radmesh%radfact(1)+2.d0*func(2)*radmesh%radfact(2) &
896 51853799 : & -0.25d0*func(3)*radmesh%radfact(3))
897 : end if
898 :
899 89191545 : intg=simp+resid
900 :
901 89191545 : end subroutine simp_gen
902 : !!***
903 :
904 : !----------------------------------------------------------------------
905 :
906 : !!****f* m_pawrad/nderiv_gen
907 : !! NAME
908 : !! nderiv_gen
909 : !!
910 : !! FUNCTION
911 : !! Do corrected first (and -if requested- second) derivation on a given (generalized) grid.
912 : !! This routine interfaces nderiv_lin (derivation on a linear grid).
913 : !!
914 : !! INPUTS
915 : !! func(:)=input function
916 : !! radmesh <type(pawrad_type)>=data containing radial grid information
917 : !!
918 : !! OUTPUT
919 : !! der(:)= 1st derivative of input function
920 : !! [der2(:)]= -- optional -- 2nd derivative of input function
921 : !!
922 : !! NOTES
923 : !! Possible mesh types (radmesh%mesh_type)
924 : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
925 : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
926 : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
927 : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
928 : !!
929 : !! SOURCE
930 :
931 7645232 : subroutine nderiv_gen(der,func,radmesh,der2)
932 :
933 : !Arguments ------------------------------------
934 : !scalars
935 : type(pawrad_type),intent(in) :: radmesh
936 : !arrays
937 : real(dp),intent(in) :: func(:)
938 : real(dp),intent(out) :: der(:)
939 : real(dp),optional,intent(out) :: der2(:)
940 :
941 : !Local variables-------------------------------
942 : !scalars
943 : integer :: msz
944 : logical :: compute_2der
945 : character(len=500) :: msg
946 :
947 : ! *************************************************************************
948 :
949 7645232 : msz=size(func)
950 7645232 : if (size(der)/=msz.or.msz>radmesh%mesh_size) then
951 0 : msg='wrong sizes for in/out arrays!'
952 0 : LIBPAW_BUG(msg)
953 : end if
954 :
955 7645232 : compute_2der=(present(der2))
956 :
957 7645232 : if (radmesh%mesh_type==1) then
958 :
959 138 : call nderiv_lin(radmesh%rstep,func,der,msz,1)
960 138 : if (compute_2der) then
961 0 : call nderiv_lin(radmesh%rstep,func,der2,msz,2)
962 : end if
963 :
964 : else if (radmesh%mesh_type==2) then
965 :
966 7644958 : call nderiv_lin(radmesh%lstep,func,der,msz,1)
967 9622438658 : der(1:msz)=der(1:msz)/radmesh%radfact(1:msz)
968 7644958 : if (compute_2der)then
969 0 : call nderiv_lin(radmesh%lstep,func,der2,msz,2)
970 0 : der2(1:msz)=(der2(1:msz)/radmesh%radfact(1:msz)-der(1:msz))/radmesh%radfact(1:msz)
971 : end if
972 :
973 : else if (radmesh%mesh_type==3) then
974 :
975 136 : call nderiv_lin(radmesh%lstep,func(2:msz),der(2:msz),msz-1,1)
976 47600 : der(2:msz)=der(2:msz)/radmesh%radfact(2:msz)
977 136 : call pawrad_deducer0(der,msz,radmesh)
978 136 : if (compute_2der)then
979 0 : call nderiv_lin(radmesh%lstep,func(2:msz),der2(2:msz),msz-1,2)
980 0 : der2(2:msz)=(der2(2:msz)/radmesh%radfact(2:msz)-der(2:msz))/radmesh%radfact(2:msz)
981 0 : call pawrad_deducer0(der2,msz,radmesh)
982 : end if
983 :
984 : else if (radmesh%mesh_type==4) then
985 :
986 0 : call nderiv_lin(radmesh%lstep,func,der,msz,1)
987 0 : der(1:msz)=der(1:msz)/radmesh%radfact(1:msz)
988 0 : if (compute_2der)then
989 0 : call nderiv_lin(radmesh%lstep,func,der2,msz,2)
990 0 : der2(1:msz)=der2(1:msz)/radmesh%radfact(1:msz)**2-der(1:msz)/radmesh%rstep
991 : end if
992 :
993 : else if (radmesh%mesh_type==5) then
994 :
995 0 : call nderiv_lin(one,func,der,msz,1)
996 0 : der(1:msz)=der(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep)
997 0 : if (compute_2der)then
998 0 : call nderiv_lin(one,func,der2,msz,2)
999 : der2(1:msz)=der2(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep)**2-two*der(1:msz)/&
1000 0 : & (radmesh%rstep+radmesh%rad(1:msz))
1001 : end if
1002 : end if
1003 :
1004 7645232 : end subroutine nderiv_gen
1005 : !!***
1006 :
1007 : !----------------------------------------------------------------------
1008 :
1009 : !!****f* m_pawrad/nderiv_lin
1010 : !! NAME
1011 : !! nderiv_lin
1012 : !!
1013 : !! FUNCTION
1014 : !! Do corrected first (and -if requested- second) derivation on a given LINEAR grid.
1015 : !!
1016 : !! INPUTS
1017 : !! hh= radial step
1018 : !! ndim= radial mesh size
1019 : !! yy(ndim)= input function
1020 : !! norder= order of derivation (1 or 2)
1021 : !!
1022 : !! OUTPUT
1023 : !! zz(ndim)= first or second derivative of y
1024 : !!
1025 : !! SOURCE
1026 :
1027 7645232 : subroutine nderiv_lin(hh,yy,zz,ndim,norder)
1028 :
1029 : !Arguments ---------------------------------------------
1030 : !scalars
1031 : integer,intent(in) :: ndim,norder
1032 : real(dp),intent(in) :: hh
1033 : !arrays
1034 : real(dp),intent(in) :: yy(ndim)
1035 : real(dp),intent(out) :: zz(ndim)
1036 :
1037 : !Local variables ---------------------------------------
1038 : !scalars
1039 : integer :: ier,ii
1040 : real(dp) :: aa,bb,cc,h1,y1
1041 :
1042 : ! *************************************************************************
1043 :
1044 : !Initialization (common to 1st and 2nd derivative)
1045 7645232 : h1=one/(12.d0*hh)
1046 7645232 : y1=yy(ndim-4)
1047 :
1048 : !FIRST DERIVATIVE
1049 : !================
1050 7645232 : if (norder==1) then
1051 :
1052 : ! Prepare differentiation loop
1053 7645232 : bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
1054 7645232 : cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
1055 : ! Start differentiation loop
1056 9591973696 : do ii=5,ndim
1057 9584328464 : aa=bb;bb=cc
1058 9584328464 : cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
1059 9591973696 : zz(ii-4)=aa
1060 : end do
1061 : ! Normal exit
1062 7645232 : ier=0
1063 7645232 : aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
1064 7645232 : zz(ndim)=h1*(3.d0*y1-16.d0*yy(ndim-3)+36.d0*yy(ndim-2) -48.d0*yy(ndim-1)+25.d0*yy(ndim))
1065 7645232 : zz(ndim-1)=aa
1066 7645232 : zz(ndim-2)=cc
1067 7645232 : zz(ndim-3)=bb
1068 :
1069 : ! SECOND DERIVATIVE
1070 : ! =================
1071 : else
1072 0 : h1=h1/hh
1073 : ! Prepare differentiation loop
1074 0 : bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
1075 0 : cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
1076 : ! Start differentiation loop
1077 0 : do ii=5,ndim
1078 0 : aa=bb;bb=cc
1079 0 : cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
1080 0 : zz(ii-4)=aa
1081 : end do
1082 : ! Normal exit
1083 0 : ier=0
1084 0 : aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
1085 0 : zz(ndim)=h1*(11.d0*y1-56.d0*yy(ndim-3)+114.d0*yy(ndim-2) -104.d0*yy(ndim-1)+35.d0*yy(ndim))
1086 0 : zz(ndim-1)=aa
1087 0 : zz(ndim-2)=cc
1088 0 : zz(ndim-3)=bb
1089 :
1090 : end if !norder
1091 :
1092 7645232 : end subroutine nderiv_lin
1093 : !!***
1094 :
1095 : !----------------------------------------------------------------------
1096 :
1097 : !!****f* m_pawrad/bound_deriv
1098 : !! NAME
1099 : !! bound_deriv
1100 : !!
1101 : !! FUNCTION
1102 : !! Computes derivatives of a function a boundaries of interval (first and last derivative)
1103 : !!
1104 : !! INPUTS
1105 : !! func(n)= array containing function
1106 : !! mesh <type(pawrad_type)>= radial mesh and related data
1107 : !! nn= size of intervall
1108 : !!
1109 : !! OUTPUT
1110 : !! yp1,ypn= derivatives of func at r(1) and r(n)
1111 : !!
1112 : !! SOURCE
1113 :
1114 104 : subroutine bound_deriv(func,mesh,nn,yp1,ypn)
1115 :
1116 : !Arguments----------------------
1117 : integer, intent(in) :: nn
1118 : real(dp), intent(in) :: func(nn)
1119 : real(dp), intent(out) :: yp1,ypn
1120 : type(pawrad_type),intent(in) :: mesh
1121 :
1122 : !*************************************************************************
1123 :
1124 104 : if (mesh%radfact(1)>zero) then
1125 : yp1=1._dp/12._dp/mesh%stepint/mesh%radfact(1) &
1126 104 : & *(-25._dp*func(1)+48._dp*func(2)-36._dp*func(3)+16._dp*func(4)-3._dp*func(5))
1127 : else
1128 0 : yp1=(func(2)-func(1))/(mesh%rad(2)-mesh%rad(1))
1129 : end if
1130 : ypn=1._dp/12._dp/mesh%stepint &
1131 : & *( 3._dp*func(nn-4)-16._dp*func(nn-3)+36._dp*func(nn-2)-48._dp*func(nn-1) &
1132 104 : & +25._dp*func(nn))/mesh%radfact(nn)
1133 :
1134 104 : end subroutine bound_deriv
1135 :
1136 : !!***
1137 :
1138 : !----------------------------------------------------------------------
1139 :
1140 : !!****f* m_pawrad/poisson
1141 : !! NAME
1142 : !! poisson
1143 : !!
1144 : !! FUNCTION
1145 : !! Solve poisson equation for angularly dependent charge
1146 : !! distribution of angular momentum l
1147 : !! Densities and potentials are given on a (generalized) radial grid
1148 : !!
1149 : !! INPUTS
1150 : !! den(:)= electron density * (4*pi*r**2) appropriate for l
1151 : !! ll= l quantum number
1152 : !! radmesh <type(pawrad_type)>=data containing radial grid information
1153 : !! [screened_sr_separation]= --optional-- separation for screened short-range kernel
1154 : !! no screening by default
1155 : !!
1156 : !! OUTPUT
1157 : !! [qq]= --optional-- lth moment of the charge; not compatible with screened Coulomb interaction
1158 : !! rv(:)= electrostatic potential * r in (Hartree*Bohr) units
1159 : !! where v(r)=\frac{1}{2l+1}(\frac{int[(r''^(l+2))g(r'')dr'']} {r^(l+1)}
1160 : !! +(r^l) int[r''^(1-l)g(r'')dr''])
1161 : !!
1162 : !! SOURCE
1163 :
1164 296755 : subroutine poisson(den,ll,radmesh,rv,screened_sr_separation,qq)
1165 :
1166 : !Arguments ---------------------------------------------
1167 : !scalars
1168 : integer,intent(in) :: ll
1169 : real(dp),intent(out),optional :: qq
1170 : real(dp),intent(in),optional :: screened_sr_separation
1171 : type(pawrad_type),intent(in) :: radmesh
1172 : !arrays
1173 : real(dp),intent(in) :: den(:)
1174 : real(dp),intent(out) :: rv(:)
1175 :
1176 : !Local variables ---------------------------------------
1177 : !scalars
1178 : integer :: ir,jr,mesh_size,mm,nn
1179 : logical :: use_numerov,use_screened
1180 : real(dp) :: angm,hh,intg,qq_,ri,rj,rr,sr_omega
1181 : !arrays
1182 : real(dp) :: ee(4)
1183 296755 : real(dp),allocatable :: aa(:),bb(:),cc(:),dd(:),radl(:),radl1(:)
1184 :
1185 : ! ***************************************************************************
1186 :
1187 296755 : mesh_size=size(den)
1188 296755 : if (size(rv)/=mesh_size.or.mesh_size>radmesh%mesh_size) then
1189 0 : LIBPAW_BUG('wrong sizes!')
1190 : end if
1191 :
1192 296755 : use_numerov=(radmesh%mesh_type==1)
1193 :
1194 296755 : use_screened=.false.;sr_omega=zero
1195 296755 : if (present(screened_sr_separation)) then
1196 358 : sr_omega=screened_sr_separation
1197 358 : use_screened=(sr_omega>tol8)
1198 : ! Numerov method not coded for screened Coulomb interaction
1199 358 : if (use_numerov.and.use_screened) use_numerov=.false.
1200 : end if
1201 :
1202 : !==============================================
1203 : !UNIFORM GRID - NUMEROV ALGORITHM
1204 : !==============================================
1205 296755 : if (use_numerov) then
1206 7449 : hh=radmesh%rstep
1207 7449 : nn=radmesh%int_meshsz-1
1208 22347 : LIBPAW_ALLOCATE(aa,(nn))
1209 14898 : LIBPAW_ALLOCATE(bb,(nn))
1210 29796 : LIBPAW_ALLOCATE(cc,(nn+1))
1211 5207879 : do ir=1,nn
1212 5200430 : aa(ir)=two*hh*den(ir+1)/(ir)
1213 5207879 : bb(ir)=den(ir+1)*((ir*hh)**ll)
1214 : end do
1215 7449 : cc(1)=zero
1216 5207879 : cc(2:nn+1)=bb(1:nn)
1217 7449 : call simp_gen(qq_,cc,radmesh)
1218 7449 : qq_=qq_/dble(2*ll+1)
1219 7449 : rv(1)=aa(1)+0.1_dp*aa(2)
1220 5192981 : do ir=2,nn-1
1221 5192981 : rv(ir)=aa(ir)+0.1_dp*(aa(ir+1)+aa(ir-1))
1222 : end do
1223 7449 : rv(nn)=aa(nn)+0.1_dp*aa(nn-1)
1224 7449 : angm=dble(ll*(ll+1))
1225 7449 : rr=(nn+1)*hh
1226 7449 : rv(nn)=rv(nn)+(2.4_dp-0.2_dp*angm/((nn+1)**2))*qq_/(rr**ll)
1227 5207879 : do ir=1,nn
1228 5200430 : aa(ir)=angm/(ir*ir)
1229 5207879 : bb(ir)=2.4_dp+aa(ir)
1230 : end do
1231 5200430 : do ir=1,nn-1
1232 5200430 : cc(ir)=-1.2_dp+0.1_dp*aa(ir+1)
1233 : end do
1234 5200430 : do ir=nn,2,-1
1235 5200430 : aa(ir)=-1.2_dp+0.1_dp*aa(ir-1)
1236 : end do
1237 7449 : if (nn.eq.1) then
1238 0 : rv(2)=rv(1)/bb(1)
1239 0 : rv(1)=zero
1240 : else
1241 5200430 : do ir=2,nn
1242 5200430 : bb(ir)=bb(ir)-aa(ir)*cc(ir-1)/bb(ir-1)
1243 : end do
1244 7449 : rv(1)=rv(1)/bb(1)
1245 5200430 : do ir=2,nn
1246 5200430 : rv(ir)=(rv(ir)-aa(ir)*rv(ir-1))/bb(ir)
1247 : end do
1248 5200430 : do ir=nn-1,1,-1
1249 5200430 : rv(ir)=rv(ir)-cc(ir)*rv(ir+1)/bb(ir)
1250 : end do
1251 5207879 : do ir=nn+1,2,-1
1252 5207879 : rv(ir)=rv(ir-1)
1253 : end do
1254 7449 : rv(1)=zero
1255 44544 : if (nn+1<mesh_size) rv(nn+2:mesh_size)=zero
1256 : end if
1257 5252423 : rv(:)=half*rv(:)
1258 7449 : if (present(qq)) qq=qq_
1259 7449 : LIBPAW_DEALLOCATE(aa)
1260 7449 : LIBPAW_DEALLOCATE(bb)
1261 14898 : LIBPAW_DEALLOCATE(cc)
1262 :
1263 : ! ==============================================
1264 : ! ANY OTHER GRID - SIMPSON ALGORITHM
1265 : ! ==============================================
1266 : else
1267 289306 : nn=mesh_size;hh=third*radmesh%stepint
1268 1477793 : do while (abs(den(nn))<tol16.and.nn>radmesh%int_meshsz)
1269 1192746 : nn=nn-1
1270 : end do
1271 289306 : mm=nn;if (radmesh%mesh_type==3) mm=mm-1
1272 867918 : LIBPAW_ALLOCATE(aa,(nn))
1273 578612 : LIBPAW_ALLOCATE(bb,(nn))
1274 578612 : LIBPAW_ALLOCATE(cc,(nn))
1275 578612 : LIBPAW_ALLOCATE(dd,(nn))
1276 :
1277 289306 : if (.not.use_screened) then
1278 :
1279 : ! Standard Coulomb integral
1280 577896 : LIBPAW_ALLOCATE(radl,(nn))
1281 577896 : LIBPAW_ALLOCATE(radl1,(nn))
1282 278857225 : do jr=nn,2,-1
1283 278568277 : ir=nn-jr+1
1284 278568277 : radl(jr) =radmesh%rad(jr)**ll
1285 278568277 : radl1(jr)=radmesh%rad(jr)*radl(jr)
1286 278568277 : aa(ir)=den(jr)*radmesh%radfact(jr)*radl(jr)
1287 278857225 : bb(ir)=den(jr)*radmesh%radfact(jr)/radl1(jr)
1288 : end do
1289 288948 : radl(1)=zero;radl1(1)=zero
1290 288948 : ee(2)=aa(nn-1);ee(3)=aa(nn-2);ee(4)=aa(nn-3)
1291 288948 : call pawrad_deducer0(ee,4,radmesh)
1292 288948 : aa(nn)=ee(1)
1293 288948 : ee(2)=bb(nn-1);ee(3)=bb(nn-2);ee(4)=bb(nn-3)
1294 288948 : call pawrad_deducer0(ee,4,radmesh)
1295 288948 : bb(nn)=ee(1)
1296 288948 : cc(1)=zero;dd(1)=zero
1297 288948 : do ir=3,mm,2
1298 139194772 : cc(ir) =cc(ir-2)+hh*(aa(ir-2)+four*aa(ir-1)+aa(ir))
1299 139194772 : cc(ir-1)=cc(ir-2)+hh*(1.25_dp*aa(ir-2)+two*aa(ir-1)-quarter*aa(ir))
1300 139194772 : dd(ir) =dd(ir-2)+hh*(bb(ir-2)+four*bb(ir-1)+bb(ir))
1301 139194772 : dd(ir-1)=dd(ir-2)+hh*(1.25_dp*bb(ir-2)+two*bb(ir-1)-quarter*bb(ir))
1302 : end do
1303 288948 : if (mod(mm,2)==0) then
1304 : ! cc(mm)=cc(mm-2)+hh*(aa(mm-2)+four*aa(mm-1)+aa(mm))
1305 : ! dd(mm)=dd(mm-2)+hh*(bb(mm-2)+four*bb(mm-1)+bb(mm))
1306 167704 : cc(mm)=cc(mm-1)+hh*(1.25_dp*aa(mm-2)+two*aa(mm-1)-quarter*aa(mm))
1307 167704 : dd(mm)=dd(mm-1)+hh*(1.25_dp*bb(mm-2)+two*bb(mm-1)-quarter*bb(mm))
1308 : end if
1309 288948 : if (mm<nn) then
1310 11029 : cc(nn)=cc(mm)+half*(aa(mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1311 11029 : dd(nn)=dd(mm)+half*(bb(mm)+bb(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1312 : end if
1313 288948 : rv(1)=zero
1314 278857225 : do ir=2,nn
1315 278568277 : jr=nn-ir+1
1316 278857225 : rv(ir)=(dd(jr)*radl1(ir)+(cc(nn)-cc(jr))/radl(ir))/(two*ll+one)
1317 : end do
1318 1475645 : if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn)
1319 288948 : if (present(qq)) qq=cc(nn)/(two*ll+one)
1320 288948 : LIBPAW_DEALLOCATE(radl)
1321 288948 : LIBPAW_DEALLOCATE(radl1)
1322 : else
1323 :
1324 : ! Short-range screened Coulomb integral
1325 358 : rv(1)=zero
1326 278133 : do ir=2,nn
1327 277775 : ri=sr_omega*radmesh%rad(ir)
1328 264402876 : do jr=2,nn
1329 264125101 : rj=sr_omega*radmesh%rad(jr)
1330 264402876 : aa(jr)=den(jr)*radmesh%radfact(jr)*sr_omega*screened_coul_kernel(ll,ri,rj)
1331 : end do
1332 277775 : call pawrad_deducer0(aa(1:4),4,radmesh)
1333 277775 : intg=zero
1334 : ! Compute the integral in 2 parts because the function is not differentiable at r(ir)
1335 : ! Integral from zero to r(ir)
1336 277775 : if (ir>2) then
1337 277417 : do jr=ir-2,1,-2
1338 66031246 : intg=intg+hh*(aa(jr)+four*aa(jr+1)+aa(jr+2))
1339 : end do
1340 277417 : if (mod(ir,2)==0) intg=intg+hh*(-quarter*aa(1)+two*aa(2)+1.25_dp*aa(3))
1341 : else if (ir==2) then
1342 358 : if (nn >2) intg=intg+hh*(+1.25_dp*aa(1)+two*aa(2)-quarter*aa(3))
1343 358 : if (nn==2) intg=intg+(aa(1)+aa(2))*hh*half
1344 : end if
1345 277775 : if (mm<nn) intg=intg+half*(aa(1+nn-mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1346 : ! Integral from r(ir) to rmax
1347 277775 : if (ir<nn-2) then
1348 276701 : do jr=ir+2,nn,2
1349 65892059 : intg=intg+hh*(aa(jr-2)+four*aa(jr-1)+aa(jr))
1350 : end do
1351 276701 : if (mod(nn-ir,2)==1) then
1352 138471 : if (ir<=nn-4) then
1353 : !Use a high-order formula because points are spaced
1354 : intg=intg+hh*(251._dp*aa(nn)+646._dp*aa(nn-1)-264._dp*aa(nn-2) &
1355 138113 : & +106._dp*aa(nn-3)-19._dp*aa(nn-4))/240._dp
1356 : else
1357 358 : intg=intg+hh*(1.25_dp*aa(nn)+two*aa(nn-1)-quarter*aa(nn-2))
1358 : endif
1359 : end if
1360 1074 : else if (ir==nn-1) then
1361 358 : intg=intg+(aa(nn-1)+aa(nn))*hh*half
1362 : end if
1363 278133 : rv(ir)=intg/(two*ll+one)*radmesh%rad(ir)
1364 : end do
1365 2148 : if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn)
1366 358 : if (present(qq)) qq=zero ! Not relevant
1367 :
1368 : end if
1369 289306 : LIBPAW_DEALLOCATE(aa)
1370 289306 : LIBPAW_DEALLOCATE(bb)
1371 289306 : LIBPAW_DEALLOCATE(cc)
1372 289306 : LIBPAW_DEALLOCATE(dd)
1373 :
1374 : end if
1375 :
1376 296755 : end subroutine poisson
1377 : !!***
1378 :
1379 : !----------------------------------------------------------------------
1380 :
1381 : !!****f* m_pawrad/pawrad_ifromr
1382 : !! NAME
1383 : !! pawrad_ifromr
1384 : !!
1385 : !! FUNCTION
1386 : !! Retreive Index FROM a given R value in a radial grid
1387 : !! Grid can be regular or logarithimc
1388 : !!
1389 : !! INPUTS
1390 : !! rr=given input r value
1391 : !! radmesh <type(pawrad_type)>=data containing radial grid information
1392 : !!
1393 : !! OUTPUT
1394 : !! pawrad_ifromr=index of rr in radial grid
1395 : !!
1396 : !! NOTES
1397 : !! Possible mesh types (radmesh%mesh_type)
1398 : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
1399 : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
1400 : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
1401 : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
1402 : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
1403 : !!
1404 : !! SOURCE
1405 :
1406 11145 : function pawrad_ifromr(radmesh,rr)
1407 :
1408 : !Arguments ------------------------------------
1409 : !scalars
1410 : integer :: pawrad_ifromr
1411 : real(dp),intent(in) :: rr
1412 : type(pawrad_type),intent(in) :: radmesh
1413 :
1414 : !Local variables-------------------------------
1415 : character(len=500) :: msg
1416 :
1417 : ! *************************************************************************
1418 :
1419 11145 : if (radmesh%mesh_type==1) then
1420 158 : pawrad_ifromr=int(tol8+rr/radmesh%rstep)+1
1421 : else if (radmesh%mesh_type==2) then
1422 10463 : pawrad_ifromr=int(tol8+log(1.d0+rr/radmesh%rstep)/radmesh%lstep)+1
1423 : else if (radmesh%mesh_type==3) then
1424 484 : if (rr<radmesh%rstep) then
1425 : pawrad_ifromr=1
1426 : else
1427 484 : pawrad_ifromr=int(tol8+log(rr/radmesh%rstep)/radmesh%lstep)+2
1428 : end if
1429 : else if (radmesh%mesh_type==4) then
1430 0 : pawrad_ifromr=int(tol8+(1.d0-exp(-rr/radmesh%rstep))/radmesh%lstep)+1
1431 : else if (radmesh%mesh_type==5) then
1432 40 : pawrad_ifromr=int(tol8+(radmesh%lstep*rr)/(radmesh%rstep+rr))+1
1433 : else
1434 : ! Other values of mesh_type are not allowed (see psp7in.F90)
1435 0 : write(msg,'(a,i0)')" Unknown value of %mesh_type ",radmesh%mesh_type
1436 0 : LIBPAW_ERROR(msg)
1437 : end if
1438 :
1439 11145 : end function pawrad_ifromr
1440 : !!***
1441 :
1442 : !----------------------------------------------------------------------
1443 :
1444 : !!****f* m_pawrad/screened_coul_kernel
1445 : !! NAME
1446 : !! screened_coul_kernel
1447 : !!
1448 : !! FUNCTION
1449 : !! Evaluates the kernel function used to compute the integral
1450 : !! of a screened Coulomb potential (with erfc function)
1451 : !! See Angyan, Gerber, Marsman, J. Phys. A: Math. Gen. 39, 8613 (2006) [[cite:Angyan2006]]
1452 : !!
1453 : !! INPUTS
1454 : !! [formula]=optional; used to force formula (1=full; 2=dev. near r2=0; 3=dev near r1*r2=0)
1455 : !! order= order of the function (typically l quantum number)
1456 : !! r1,r2=input arguments
1457 : !!
1458 : !! OUTPUT
1459 : !! screened_coul_kernel=output radial function
1460 : !!
1461 : !! SOURCE
1462 :
1463 264125101 : function screened_coul_kernel(order,r1,r2,formula)
1464 :
1465 : !Arguments ------------------------------------
1466 : !scalars
1467 : integer,intent(in) :: order
1468 : integer,optional :: formula
1469 : real(dp),intent(in) :: r1,r2
1470 : real(dp) :: screened_coul_kernel
1471 : !Local variables-------------------------------
1472 : integer :: formula_
1473 : real(dp) :: difexp,erfcp,erfcm,f0,f1,f2,f3,f4,f5,f6,hh,sqrtpi,sumexp,xx,xy,yy
1474 : character(len=100) :: msg
1475 :
1476 : !******************************************************************
1477 :
1478 264125101 : if (order>6) then
1479 0 : msg='PAW screened exchange not coded for l>2!'
1480 0 : LIBPAW_ERROR(msg)
1481 : end if
1482 :
1483 : !Use max and min of arguments
1484 264125101 : xx=max(r1,r2) ; yy=min(r1,r2)
1485 :
1486 : !Unscreened Coulomb interaction for very small (or negative) arguments
1487 264125101 : if (xx<tol8) then
1488 0 : screened_coul_kernel=yy**order/xx**(order+1) ; return
1489 : end if
1490 :
1491 : !Choice of formula
1492 : !Empirical criterion, inspired by J. Phys. A 39, pp8624 [[cite:Angyan2006]] and adjusted
1493 264125101 : formula_=1;if (xx<0.25_dp.and.yy<0.25_dp) formula_=2
1494 264125101 : if (present(formula)) formula_=formula
1495 0 : select case (formula_)
1496 :
1497 : !Full formula
1498 : !J. Phys. A 39, 8613 (2006) - Eq. (21), (22), (24) [[cite:Angyan2006]]
1499 : ! Note: typo in the paper: in eq (22), the sum begins at p=0
1500 : !------------------------------------------------------------------
1501 : case(1)
1502 0 : xy=xx*yy ; sqrtpi=sqrt(pi)
1503 0 : sumexp=exp(-(xx+yy)**2)+exp(-(xx-yy)**2)
1504 0 : difexp=exp(-(xx+yy)**2)-exp(-(xx-yy)**2)
1505 0 : erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy)
1506 264125101 : if (order==0) then
1507 0 : f0=-difexp/(two*sqrtpi*xy)
1508 0 : hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy)
1509 0 : screened_coul_kernel = hh + f0
1510 : else if(order==1) then
1511 0 : f0=-difexp/(two*sqrtpi*xy)
1512 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1513 0 : hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2)
1514 0 : screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy
1515 : else if (order==2) then
1516 0 : f0=-difexp/(two*sqrtpi*xy)
1517 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1518 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1519 0 : hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3)
1520 0 : screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2
1521 : else if (order==3) then
1522 0 : f0=-difexp/(two*sqrtpi*xy)
1523 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1524 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1525 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1526 0 : hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4)
1527 : screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 &
1528 0 : & + f0*(xx**6+yy**6)/xy**3
1529 : else if (order==4) then
1530 0 : f0=-difexp/(two*sqrtpi*xy)
1531 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1532 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1533 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1534 : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1535 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1536 0 : hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5)
1537 : screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 &
1538 0 : & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4
1539 : else if (order==5) then
1540 0 : f0=-difexp/(two*sqrtpi*xy)
1541 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1542 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1543 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1544 : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1545 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1546 : f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp &
1547 0 : & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6)
1548 0 : hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6)
1549 : screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 &
1550 : & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 &
1551 0 : & + f0*(xx**10+yy**10)/xy**5
1552 : else if (order==6) then
1553 0 : f0=-difexp/(two*sqrtpi*xy)
1554 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1555 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1556 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1557 : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1558 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1559 : f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp &
1560 0 : & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6)
1561 : f6=-((64._dp*xy**6+3360._dp*xy**4+18900._dp*xy**2+10395._dp)*difexp+ &
1562 0 : & (672._dp*xy**5+10080._dp*xy**3+20790._dp*xy)*sumexp)/(sqrtpi*(two*xy)**7)
1563 0 : hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7)
1564 : screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 &
1565 : & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 &
1566 0 : & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6
1567 : end if
1568 :
1569 : !Development for yy->0
1570 : !J. Phys. A 39, 8613 (2006) - Eq. (28), (29), (30) [[cite:Angyan2006]]
1571 : !------------------------------------------------------------------
1572 : case(2)
1573 264125101 : sqrtpi=sqrt(pi)
1574 0 : if (order==0) then
1575 : screened_coul_kernel = paw_derfc(xx)/xx + exp(-xx**2)/sqrtpi * &
1576 : & (two/three *yy**2 &
1577 113142729 : & +(two*xx**2-three)/15._dp*yy**4)
1578 : else if(order==1) then
1579 : screened_coul_kernel = paw_derfc(xx)*yy/xx**2 + exp(-xx**2)/sqrtpi * &
1580 : & (two/xx *yy &
1581 : & +four*xx/5._dp *yy**3 &
1582 55508225 : & +two*xx*(two*xx**2-5._dp)/35._dp*yy**5)
1583 : else if (order==2) then
1584 : screened_coul_kernel = paw_derfc(xx)*yy**2/xx**3 + exp(-xx**2)/sqrtpi * &
1585 : & (two*(two*xx**2+three)/(three*xx**2) *yy**2 &
1586 : & +8._dp*xx**2/21._dp *yy**4 &
1587 95474147 : & +four*xx**2*(two*xx**2-7._dp)/189._dp*yy**6)
1588 : else if (order==3) then
1589 : screened_coul_kernel = paw_derfc(xx)*yy**3/xx**4 + exp(-xx**2)/sqrtpi * &
1590 : & (two*(four*xx**4+10._dp*xx**2+15._dp)/(15._dp*xx**3)*yy**3 &
1591 : & +16._dp*xx**3/135._dp *yy**5 &
1592 0 : & +8._dp*xx**3*(two*xx**2-9._dp)/1485._dp *yy**7)
1593 : else if (order==4) then
1594 : screened_coul_kernel = paw_derfc(xx)*yy**4/xx**5 + exp(-xx**2)/sqrtpi * &
1595 : & (two*(8._dp*xx**6+28._dp*xx**4+70._dp*xx**2+105._dp)/(105._dp*xx**4)*yy**4 &
1596 : & +32._dp*xx**4/1155._dp *yy**6 &
1597 0 : & +16._dp*xx**4*(two*xx**2-11._dp)/15015._dp *yy**8)
1598 : else if (order==5) then
1599 : screened_coul_kernel = paw_derfc(xx)*yy**5/xx**6 + exp(-xx**2)/sqrtpi * &
1600 : & (two*(16._dp*xx**8+72._dp*xx**6+252._dp*xx**4+630._dp*xx**2+945._dp)/(945._dp*xx**5)*yy**5 &
1601 : & +64._dp*xx**5/12285._dp *yy**7 &
1602 0 : & +32._dp*xx**5*(two*xx**2-13._dp)/184275._dp *yy**9)
1603 : else if (order==6) then
1604 : screened_coul_kernel = paw_derfc(xx)*yy**6/xx**7 + exp(-xx**2)/sqrtpi * &
1605 : & (two*(32._dp*xx**10+176._dp*xx**8+792._dp*xx**6+2772._dp*xx**4+6930._dp*xx**2+10395._dp)/(10395._dp*xx**6)*yy**6 &
1606 : & +128._dp*xx**6/155925._dp *yy**8 &
1607 0 : & +64._dp*xx**6*(two*xx**2-15._dp)/2650725._dp *yy**10)
1608 : end if
1609 :
1610 : !Development for xx*yy->0
1611 : !J. Phys. A 39, 8613 (2006) - Eq. (24), (26) [[cite:Angyan2006]]
1612 : ! Note: typo in the paper: (2n+3)! should be (2n+3)!!
1613 : !------------------------------------------------------------------
1614 : case(3)
1615 0 : xy=xx*yy ; sqrtpi=sqrt(pi)
1616 0 : sumexp=exp(-(xx**2+yy**2))
1617 0 : erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy)
1618 264125101 : if (order==0) then
1619 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1620 0 : hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy)
1621 0 : screened_coul_kernel = hh + f0
1622 : else if(order==1) then
1623 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1624 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1625 0 : hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2)
1626 0 : screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy
1627 : else if (order==2) then
1628 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1629 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1630 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1631 0 : hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3)
1632 0 : screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2
1633 : else if (order==3) then
1634 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1635 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1636 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1637 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1638 0 : hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4)
1639 : screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 &
1640 0 : & + f0*(xx**6+yy**6)/xy**3
1641 : else if (order==4) then
1642 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1643 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1644 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1645 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1646 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1647 0 : hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5)
1648 : screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 &
1649 0 : & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4
1650 : else if (order==5) then
1651 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1652 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1653 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1654 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1655 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1656 0 : f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp
1657 0 : hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6)
1658 : screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 &
1659 : & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 &
1660 0 : & + f0*(xx**10+yy**10)/xy**5
1661 : else if (order==6) then
1662 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1663 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1664 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1665 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1666 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1667 0 : f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp
1668 0 : f6 = 128._dp*xy**2*(15._dp+two*xy**2)/(1027025._dp*sqrtpi)*sumexp
1669 0 : hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7)
1670 : screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 &
1671 : & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 &
1672 0 : & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6
1673 : end if
1674 :
1675 : end select
1676 :
1677 : end function screened_coul_kernel
1678 : !!***
1679 :
1680 : !----------------------------------------------------------------------
1681 :
1682 : !!****f* m_pawrad/calc_slatradl
1683 : !! NAME
1684 : !! calc_slatradl
1685 : !!
1686 : !! FUNCTION
1687 : !! Calculate the radial part of Slater integrals. See below.
1688 : !!
1689 : !! INPUTS
1690 : !! ll= l quantum number in the expansion of the Coulomb term.
1691 : !! mesh_size=Number of points on the radial mesh.
1692 : !! ff1(radmesh), ff2(radmesh)= The two functions to be integrated.
1693 : !! Pawrad <type(pawrad_type)>=Structure containing radial grid information.
1694 : !!
1695 : !! OUTPUT
1696 : !! integral =
1697 : !! $ \dfrac{4\pi}{2L+1} \int ff1(r1) \dfrac{r_<^L}{r_>^{L+1}} ff2(r2) dr1 dr2 $
1698 : !! where $r_< = min(r1,r2)$ and $r_> = Max(r1,r2)$.
1699 : !!
1700 : !! SOURCE
1701 :
1702 162 : subroutine calc_slatradl(ll,mesh_size,ff1,ff2,Pawrad,integral)
1703 :
1704 : !scalars
1705 : integer,intent(in) :: mesh_size,ll
1706 : real(dp),intent(out) :: integral
1707 : !arrays
1708 : real(dp),intent(in) :: ff1(mesh_size),ff2(mesh_size)
1709 : type(pawrad_type),intent(in) :: Pawrad
1710 :
1711 : !Local variables ---------------------------------------
1712 : !scalars
1713 : integer :: int_meshsz
1714 : character(len=100) :: msg
1715 : !arrays
1716 162 : real(dp),allocatable :: hh(:),gg(:)
1717 :
1718 : ! *************************************************************************
1719 :
1720 162 : if (mesh_size > Pawrad%mesh_size) then
1721 0 : msg='mesh_size > pawrad%mesh_size!'
1722 0 : LIBPAW_BUG(msg)
1723 : end if
1724 :
1725 486 : LIBPAW_ALLOCATE(hh,(mesh_size))
1726 324 : LIBPAW_ALLOCATE(gg,(mesh_size))
1727 : !hh = zero
1728 : !gg = zero
1729 :
1730 162 : int_meshsz=Pawrad%int_meshsz
1731 : !$int_meshsz=Pawrad%mesh_size
1732 :
1733 : ! the line below requires hh as work array.
1734 224208 : hh = ff2
1735 :
1736 : ! TODO find where int_meshsz is calculated and if it can affects the results.
1737 162 : if (int_meshsz<mesh_size) hh(int_meshsz+1:mesh_size)=zero
1738 :
1739 162 : call poisson(hh,ll,Pawrad,gg)
1740 :
1741 223884 : gg(2:mesh_size) = gg(2:mesh_size)/Pawrad%rad(2:mesh_size)
1742 :
1743 162 : hh(1) = zero
1744 223884 : hh(2:mesh_size)= ff1(2:mesh_size) * gg(2:mesh_size)
1745 162 : LIBPAW_DEALLOCATE(gg)
1746 :
1747 162 : call simp_gen(integral,hh,Pawrad)
1748 162 : integral = four_pi * integral
1749 :
1750 162 : LIBPAW_DEALLOCATE(hh)
1751 :
1752 162 : end subroutine calc_slatradl
1753 : !!***
1754 :
1755 : !----------------------------------------------------------------------
1756 :
1757 0 : END MODULE m_pawrad
1758 : !!***
|
