TABLE OF CONTENTS
- ABINIT/m_pawrad
- m_pawrad/bound_deriv
- m_pawrad/calc_slatradl
- m_pawrad/nderiv_gen
- m_pawrad/nderiv_lin
- m_pawrad/pawrad_bcast
- m_pawrad/pawrad_copy
- m_pawrad/pawrad_deducer0
- m_pawrad/pawrad_free_0D
- m_pawrad/pawrad_free_1D
- m_pawrad/pawrad_ifromr
- m_pawrad/pawrad_init
- m_pawrad/pawrad_isame
- m_pawrad/pawrad_print
- m_pawrad/pawrad_type
- m_pawrad/poisson
- m_pawrad/screened_coul_kernel
- m_pawrad/simp_gen
ABINIT/m_pawrad [ Modules ]
NAME
m_pawrad
FUNCTION
Module containing all the functions related to the PAW radial meshes
COPYRIGHT
Copyright (C) 2013-2024 ABINIT group (MT,FJ,MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
* Routines tagged with "@type_name" are strongly connected to the definition of the data type.
Strongly connected means that the proper functioning of the implementation relies on the
assumption that the tagged procedure is consistent with the type declaration.
Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
that all the strongly connected routines are changed accordingly to accommodate the modification of the data type
Typical examples of strongly connected routines are creation, destruction or reset methods.
* FOR DEVELOPERS: in order to preserve the portability of libPAW library,
please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
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
m_pawrad/bound_deriv [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
bound_deriv
FUNCTION
Computes derivatives of a function a boundaries of interval (first and last derivative)
INPUTS
func(n)= array containing function mesh <type(pawrad_type)>= radial mesh and related data nn= size of intervall
OUTPUT
yp1,ypn= derivatives of func at r(1) and r(n)
SOURCE
1114 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 if (mesh%radfact(1)>zero) then 1125 yp1=1._dp/12._dp/mesh%stepint/mesh%radfact(1) & 1126 & *(-25._dp*func(1)+48._dp*func(2)-36._dp*func(3)+16._dp*func(4)-3._dp*func(5)) 1127 else 1128 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 & +25._dp*func(nn))/mesh%radfact(nn) 1133 1134 end subroutine bound_deriv
m_pawrad/calc_slatradl [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
calc_slatradl
FUNCTION
Calculate the radial part of Slater integrals. See below.
INPUTS
ll= l quantum number in the expansion of the Coulomb term. mesh_size=Number of points on the radial mesh. ff1(radmesh), ff2(radmesh)= The two functions to be integrated. Pawrad <type(pawrad_type)>=Structure containing radial grid information.
OUTPUT
integral =
$ \dfrac{4\pi}{2L+1} \int ff1(r1) \dfrac{r_<^L}{r_>^{L+1}} ff2(r2) dr1 dr2 $
where $r_< = min(r1,r2)$ and $r_> = Max(r1,r2)$.
SOURCE
1702 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 real(dp),allocatable :: hh(:),gg(:) 1717 1718 ! ************************************************************************* 1719 1720 if (mesh_size > Pawrad%mesh_size) then 1721 msg='mesh_size > pawrad%mesh_size!' 1722 LIBPAW_BUG(msg) 1723 end if 1724 1725 LIBPAW_ALLOCATE(hh,(mesh_size)) 1726 LIBPAW_ALLOCATE(gg,(mesh_size)) 1727 !hh = zero 1728 !gg = zero 1729 1730 int_meshsz=Pawrad%int_meshsz 1731 !$int_meshsz=Pawrad%mesh_size 1732 1733 ! the line below requires hh as work array. 1734 hh = ff2 1735 1736 ! TODO find where int_meshsz is calculated and if it can affects the results. 1737 if (int_meshsz<mesh_size) hh(int_meshsz+1:mesh_size)=zero 1738 1739 call poisson(hh,ll,Pawrad,gg) 1740 1741 gg(2:mesh_size) = gg(2:mesh_size)/Pawrad%rad(2:mesh_size) 1742 1743 hh(1) = zero 1744 hh(2:mesh_size)= ff1(2:mesh_size) * gg(2:mesh_size) 1745 LIBPAW_DEALLOCATE(gg) 1746 1747 call simp_gen(integral,hh,Pawrad) 1748 integral = four_pi * integral 1749 1750 LIBPAW_DEALLOCATE(hh) 1751 1752 end subroutine calc_slatradl
m_pawrad/nderiv_gen [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
nderiv_gen
FUNCTION
Do corrected first (and -if requested- second) derivation on a given (generalized) grid. This routine interfaces nderiv_lin (derivation on a linear grid).
INPUTS
func(:)=input function radmesh <type(pawrad_type)>=data containing radial grid information
OUTPUT
der(:)= 1st derivative of input function [der2(:)]= -- optional -- 2nd derivative of input function
NOTES
Possible mesh types (radmesh%mesh_type) mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
SOURCE
931 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 msz=size(func) 950 if (size(der)/=msz.or.msz>radmesh%mesh_size) then 951 msg='wrong sizes for in/out arrays!' 952 LIBPAW_BUG(msg) 953 end if 954 955 compute_2der=(present(der2)) 956 957 if (radmesh%mesh_type==1) then 958 959 call nderiv_lin(radmesh%rstep,func,der,msz,1) 960 if (compute_2der) then 961 call nderiv_lin(radmesh%rstep,func,der2,msz,2) 962 end if 963 964 else if (radmesh%mesh_type==2) then 965 966 call nderiv_lin(radmesh%lstep,func,der,msz,1) 967 der(1:msz)=der(1:msz)/radmesh%radfact(1:msz) 968 if (compute_2der)then 969 call nderiv_lin(radmesh%lstep,func,der2,msz,2) 970 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 call nderiv_lin(radmesh%lstep,func(2:msz),der(2:msz),msz-1,1) 976 der(2:msz)=der(2:msz)/radmesh%radfact(2:msz) 977 call pawrad_deducer0(der,msz,radmesh) 978 if (compute_2der)then 979 call nderiv_lin(radmesh%lstep,func(2:msz),der2(2:msz),msz-1,2) 980 der2(2:msz)=(der2(2:msz)/radmesh%radfact(2:msz)-der(2:msz))/radmesh%radfact(2:msz) 981 call pawrad_deducer0(der2,msz,radmesh) 982 end if 983 984 else if (radmesh%mesh_type==4) then 985 986 call nderiv_lin(radmesh%lstep,func,der,msz,1) 987 der(1:msz)=der(1:msz)/radmesh%radfact(1:msz) 988 if (compute_2der)then 989 call nderiv_lin(radmesh%lstep,func,der2,msz,2) 990 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 call nderiv_lin(one,func,der,msz,1) 996 der(1:msz)=der(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep) 997 if (compute_2der)then 998 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 & (radmesh%rstep+radmesh%rad(1:msz)) 1001 end if 1002 end if 1003 1004 end subroutine nderiv_gen
m_pawrad/nderiv_lin [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
nderiv_lin
FUNCTION
Do corrected first (and -if requested- second) derivation on a given LINEAR grid.
INPUTS
hh= radial step ndim= radial mesh size yy(ndim)= input function norder= order of derivation (1 or 2)
OUTPUT
zz(ndim)= first or second derivative of y
SOURCE
1027 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 h1=one/(12.d0*hh) 1046 y1=yy(ndim-4) 1047 1048 !FIRST DERIVATIVE 1049 !================ 1050 if (norder==1) then 1051 1052 ! Prepare differentiation loop 1053 bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5)) 1054 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 do ii=5,ndim 1057 aa=bb;bb=cc 1058 cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3))) 1059 zz(ii-4)=aa 1060 end do 1061 ! Normal exit 1062 ier=0 1063 aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim)) 1064 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 zz(ndim-1)=aa 1066 zz(ndim-2)=cc 1067 zz(ndim-3)=bb 1068 1069 ! SECOND DERIVATIVE 1070 ! ================= 1071 else 1072 h1=h1/hh 1073 ! Prepare differentiation loop 1074 bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5)) 1075 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 do ii=5,ndim 1078 aa=bb;bb=cc 1079 cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2)) 1080 zz(ii-4)=aa 1081 end do 1082 ! Normal exit 1083 ier=0 1084 aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim)) 1085 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 zz(ndim-1)=aa 1087 zz(ndim-2)=cc 1088 zz(ndim-3)=bb 1089 1090 end if !norder 1091 1092 end subroutine nderiv_lin
m_pawrad/pawrad_bcast [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_bcast
FUNCTION
Communicate pawrad data over a given MPI communicator
INPUTS
comm_mpi= communicator used to broadcast data
SIDE EFFECTS
pawrad=<type pawrad_type>= a radial mesh datastructure for PAW
SOURCE
648 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 integer,allocatable :: list_int(:) 662 real(dp),allocatable :: list_dpr(:) 663 664 !************************************************************************* 665 666 me=xmpi_comm_rank(comm_mpi) 667 668 !Initializations 669 if_rad=0; if_radfact=0; if_simfact=0 670 671 !calculate the size of the reals 672 if(me==0) then 673 if (allocated(pawrad%rad)) then 674 if_rad=1 !communicate rad 675 isz1=size(pawrad%rad) 676 if(isz1/=pawrad%mesh_size) then 677 msg='rad: sz1 /= pawrad%mesh_size (1)' 678 LIBPAW_BUG(msg) 679 end if 680 end if 681 if (allocated(pawrad%radfact)) then 682 if_radfact=1 !communicate radfact 683 isz1=size(pawrad%radfact) 684 if(isz1/=pawrad%mesh_size) then 685 msg='radfact: sz1 /= pawrad%mesh_size (2)' 686 LIBPAW_BUG(msg) 687 end if 688 end if 689 if (allocated(pawrad%simfact)) then 690 if_simfact=1 !communicate simfact 691 isz1=size(pawrad%simfact) 692 if(isz1/=pawrad%mesh_size) then 693 msg='simfact: sz1 /= pawrad%mesh_size (3)' 694 LIBPAW_BUG(msg) 695 end if 696 end if 697 end if 698 699 !Brodcast the integers 700 LIBPAW_ALLOCATE(list_int,(6)) 701 if(me==0) then 702 list_int(1)=pawrad%int_meshsz 703 list_int(2)=pawrad%mesh_size 704 list_int(3)=pawrad%mesh_type 705 list_int(4)=if_rad 706 list_int(5)=if_radfact 707 list_int(6)=if_simfact 708 end if 709 call xmpi_bcast(list_int,0,comm_mpi,ierr) 710 if(me/=0) then 711 pawrad%int_meshsz =list_int(1) 712 pawrad%mesh_size =list_int(2) 713 pawrad%mesh_type =list_int(3) 714 if_rad=list_int(4) 715 if_radfact=list_int(5) 716 if_simfact=list_int(6) 717 end if 718 LIBPAW_DEALLOCATE(list_int) 719 720 !Broadcast the reals 721 nn=4+pawrad%mesh_size*(if_rad+if_radfact+if_simfact) 722 LIBPAW_ALLOCATE(list_dpr,(nn)) 723 if(me==0) then 724 list_dpr(1)=pawrad%lstep 725 list_dpr(2)=pawrad%rmax 726 list_dpr(3)=pawrad%rstep 727 list_dpr(4)=pawrad%stepint 728 indx=5 729 if (if_rad==1) then 730 isz1=pawrad%mesh_size 731 list_dpr(indx:indx+isz1-1)=pawrad%rad(1:isz1) 732 indx=indx+isz1 733 end if 734 if (if_radfact==1) then 735 isz1=pawrad%mesh_size 736 list_dpr(indx:indx+isz1-1)=pawrad%radfact(1:isz1) 737 indx=indx+isz1 738 end if 739 if (if_simfact==1) then 740 isz1=pawrad%mesh_size 741 list_dpr(indx:indx+isz1-1)=pawrad%simfact(1:isz1) 742 indx=indx+isz1 743 end if 744 end if 745 call xmpi_bcast(list_dpr,0,comm_mpi,ierr) 746 if(me/=0) then 747 pawrad%lstep=list_dpr(1) 748 pawrad%rmax=list_dpr(2) 749 pawrad%rstep=list_dpr(3) 750 pawrad%stepint=list_dpr(4) 751 indx=5 752 ! Deallocate all arrays: 753 if (allocated(pawrad%rad)) then 754 LIBPAW_DEALLOCATE(pawrad%rad) 755 end if 756 if (allocated(pawrad%radfact)) then 757 LIBPAW_DEALLOCATE(pawrad%radfact) 758 end if 759 if (allocated(pawrad%simfact)) then 760 LIBPAW_DEALLOCATE(pawrad%simfact) 761 end if 762 ! Communicate if flag is set to 1: 763 if(if_rad==1) then 764 isz1=pawrad%mesh_size 765 LIBPAW_ALLOCATE(pawrad%rad,(isz1)) 766 pawrad%rad(1:isz1)=list_dpr(indx:indx+isz1-1) 767 indx=indx+isz1 768 end if 769 if(if_radfact==1) then 770 isz1=pawrad%mesh_size 771 LIBPAW_ALLOCATE(pawrad%radfact,(isz1)) 772 pawrad%radfact(1:isz1)=list_dpr(indx:indx+isz1-1) 773 indx=indx+isz1 774 end if 775 if(if_simfact==1) then 776 isz1=pawrad%mesh_size 777 LIBPAW_ALLOCATE(pawrad%simfact,(isz1)) 778 pawrad%simfact(1:isz1)=list_dpr(indx:indx+isz1-1) 779 indx=indx+isz1 780 end if 781 end if 782 LIBPAW_DEALLOCATE(list_dpr) 783 784 end subroutine pawrad_bcast
m_pawrad/pawrad_copy [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_copy
FUNCTION
Copy one radial mesh (in a generalized format) to another
INPUTS
mesh1 <type(pawrad_type)>=data containing radial grid information of input mesh
OUTPUT
mesh2 <type(pawrad_type)>=data containing radial grid information of output mesh
NOTES
Possible mesh types (mesh%mesh_type) mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
SOURCE
559 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 mesh2%mesh_type =mesh1%mesh_type 573 mesh2%mesh_size =mesh1%mesh_size 574 mesh2%int_meshsz=mesh1%int_meshsz 575 mesh2%lstep =mesh1%lstep 576 mesh2%rstep =mesh1%rstep 577 mesh2%stepint =mesh1%stepint 578 mesh2%rmax =mesh1%rmax 579 580 LIBPAW_ALLOCATE(mesh2%rad,(mesh1%mesh_size)) 581 LIBPAW_ALLOCATE(mesh2%radfact,(mesh1%mesh_size)) 582 LIBPAW_ALLOCATE(mesh2%simfact,(mesh1%mesh_size)) 583 do ir=1,mesh1%mesh_size 584 mesh2%rad(ir) =mesh1%rad(ir) 585 mesh2%radfact(ir)=mesh1%radfact(ir) 586 mesh2%simfact(ir)=mesh1%simfact(ir) 587 end do 588 589 end subroutine pawrad_copy
m_pawrad/pawrad_deducer0 [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_deducer0
FUNCTION
Extrapolate r=0 value of a function from values near r=0 using a 3 points formula
INPUTS
funcsz=size of array func radmesh <type(pawrad_type)>=data containing radial grid information
SIDE EFFECTS
func(funcsz)=array containing values of function to extrapolate
SOURCE
611 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 if (radmesh%mesh_type==1.or.radmesh%mesh_type==2.or.radmesh%mesh_type==4.or.radmesh%mesh_type==5) then 623 func(1)=func(4)+3*(func(2)-func(3)) 624 else if (radmesh%mesh_type==3) then 625 func(1)=func(4)+exp(two*radmesh%lstep)/(exp(radmesh%lstep)-one)*(func(2)-func(3)) 626 end if 627 628 end subroutine pawrad_deducer0
m_pawrad/pawrad_free_0D [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_free_0D
FUNCTION
Frees all memory allocated in the object
SOURCE
312 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 if (allocated(Rmesh%rad )) then 325 LIBPAW_DEALLOCATE(Rmesh%rad) 326 end if 327 if (allocated(Rmesh%radfact)) then 328 LIBPAW_DEALLOCATE(Rmesh%radfact) 329 end if 330 if (allocated(Rmesh%simfact)) then 331 LIBPAW_DEALLOCATE(Rmesh%simfact) 332 end if 333 Rmesh%int_meshsz=0 334 Rmesh%mesh_size=0 335 Rmesh%mesh_type=-1 336 337 end subroutine pawrad_free_0D
m_pawrad/pawrad_free_1D [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_free_1D
FUNCTION
Destroy all objects in an array of pawrad data structures
SOURCE
351 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 nn=size(Rmesh) 364 if (nn==0) return 365 366 do ii=1,nn 367 call pawrad_free_0D(Rmesh(ii)) 368 end do 369 370 end subroutine pawrad_free_1D
m_pawrad/pawrad_ifromr [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_ifromr
FUNCTION
Retreive Index FROM a given R value in a radial grid Grid can be regular or logarithimc
INPUTS
rr=given input r value radmesh <type(pawrad_type)>=data containing radial grid information
OUTPUT
pawrad_ifromr=index of rr in radial grid
NOTES
Possible mesh types (radmesh%mesh_type) mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
SOURCE
1406 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 if (radmesh%mesh_type==1) then 1420 pawrad_ifromr=int(tol8+rr/radmesh%rstep)+1 1421 else if (radmesh%mesh_type==2) then 1422 pawrad_ifromr=int(tol8+log(1.d0+rr/radmesh%rstep)/radmesh%lstep)+1 1423 else if (radmesh%mesh_type==3) then 1424 if (rr<radmesh%rstep) then 1425 pawrad_ifromr=1 1426 else 1427 pawrad_ifromr=int(tol8+log(rr/radmesh%rstep)/radmesh%lstep)+2 1428 end if 1429 else if (radmesh%mesh_type==4) then 1430 pawrad_ifromr=int(tol8+(1.d0-exp(-rr/radmesh%rstep))/radmesh%lstep)+1 1431 else if (radmesh%mesh_type==5) then 1432 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 write(msg,'(a,i0)')" Unknown value of %mesh_type ",radmesh%mesh_type 1436 LIBPAW_ERROR(msg) 1437 end if 1438 1439 end function pawrad_ifromr
m_pawrad/pawrad_init [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_init
FUNCTION
Creation method for radial meshes. Compute all points (and related weights) of a radial mesh. Grid can be regular or logarithimc.
INPUTS
[mesh_size]=Dimension of the radial mesh
If not present, take mesh%mesh_size
[mesh_type]=Type of mesh
If not present, take mesh%mesh_type
[rstep]=Radial step of the mesh (AA parameter above)
If not present, take mesh%rstep
[lstep]=Exponential step of the mesh (BB parameter above)
If not present, take mesh%lstep
Needed only if mesh type is logarithmic.
[r_for_intg]=Mesh size used in integrals computation
If not present, take mesh%r_for_intg
Integrals will be computed up to rr(r_for_intg)
(can be negative for an integration over the whole grid)
OUTPUT
SIDE EFFECTS
mesh<pawrad_type>=The object completely initialized (containing radial grid information). The following quantities are calculated inside the routine: %stepint = Radial step used to convert any function from the present grid onto a regular grid in order to integrate it using trapeze method %rad(mesh_size) = Coordinates of all the points of the mesh. %radfact(mesh_size) = Factors used to compute radial integrals. %int_meshsz = Integrals will be computed up to r(int_meshsz) %simfact(mesh_size) = Factor used to compute radial integrals by the a Simpson scheme Integral[f] = Sum_i [simfact(i)*f(i)] %rmax = Max. value of r = rad(mesh_size)
NOTES
Possible mesh types (mesh%mesh_type) mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
SOURCE
190 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 mesh_size_ = mesh%mesh_size ; if (present(mesh_size)) mesh_size_ = mesh_size 211 mesh_type_ = mesh%mesh_type ; if (present(mesh_type)) mesh_type_ = mesh_type 212 rstep_ = mesh%rstep ; if (present(rstep)) rstep_ = rstep 213 lstep_ = mesh%lstep ; if (present(lstep)) lstep_ = lstep 214 215 r_for_intg_=-1._dp;if (present(r_for_intg)) r_for_intg_=r_for_intg 216 217 mesh%mesh_size = mesh_size_ 218 mesh%mesh_type = mesh_type_ 219 mesh%rstep = rstep_ 220 mesh%lstep = lstep_ 221 LIBPAW_ALLOCATE(mesh%rad ,(mesh%mesh_size)) 222 LIBPAW_ALLOCATE(mesh%radfact,(mesh%mesh_size)) 223 LIBPAW_ALLOCATE(mesh%simfact,(mesh%mesh_size)) 224 mesh%simfact=zero 225 if (mesh%mesh_type==1) then 226 isim=3 227 mesh%stepint=mesh%rstep 228 mesh%rad(1)=zero;mesh%radfact(1)=one 229 do ir=2,mesh%mesh_size 230 mesh%rad(ir) =mesh%rstep*dble(ir-1) 231 mesh%radfact(ir)=one 232 end do 233 else if (mesh%mesh_type==2) then 234 isim=3 235 mesh%stepint=mesh%lstep 236 mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep 237 do ir=2,mesh%mesh_size 238 mesh%rad(ir) =mesh%rstep*(exp(mesh%lstep*dble(ir-1))-one) 239 mesh%radfact(ir)=mesh%rad(ir)+mesh%rstep 240 end do 241 else if (mesh%mesh_type==3) then 242 isim=4 243 mesh%stepint=mesh%lstep 244 mesh%rad(1)=zero;mesh%radfact(1)=zero 245 do ir=2,mesh%mesh_size 246 mesh%rad(ir) =mesh%rstep*exp(mesh%lstep*dble(ir-2)) 247 mesh%radfact(ir)=mesh%rad(ir) 248 end do 249 else if (mesh%mesh_type==4) then 250 isim=3 251 mesh%lstep=one/dble(mesh%mesh_size) 252 mesh%stepint=mesh%lstep 253 mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep 254 do ir=2,mesh%mesh_size 255 mesh%rad(ir) =-mesh%rstep*log(one-mesh%lstep*dble(ir-1)) 256 mesh%radfact(ir)=mesh%rstep/(one-mesh%lstep*dble(ir-1)) 257 end do 258 else if (mesh%mesh_type==5) then 259 isim=3 260 mesh%stepint=mesh%rstep 261 mesh%rad(1)=zero;mesh%radfact(1)=1/mesh%lstep 262 do ir=2,mesh%mesh_size 263 mesh%rad(ir) =mesh%rstep*dble(ir-1)/(mesh%lstep-dble(ir-1)) 264 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 write(msg,'(a,i0)')" Unknown value of mesh_type: ",mesh%mesh_type 269 LIBPAW_ERROR(msg) 270 end if 271 272 mesh%int_meshsz=mesh%mesh_size 273 if (r_for_intg_>0.d0) then 274 ir=min(pawrad_ifromr(mesh,r_for_intg_),mesh%mesh_size) 275 if (ir<mesh%mesh_size) then 276 if (abs(mesh%rad(ir+1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir+1 277 end if 278 if (ir>1) then 279 if (abs(mesh%rad(ir-1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir-1 280 end if 281 mesh%int_meshsz=ir 282 end if 283 284 hh=mesh%stepint/3.d0 285 mesh%simfact(mesh%int_meshsz)=hh*mesh%radfact(mesh%int_meshsz) 286 mesh%simfact(1:isim-2)=zero 287 ir_last=1 288 do ir=mesh%int_meshsz,isim,-2 289 mesh%simfact(ir-1)=4.d0*hh*mesh%radfact(ir-1) 290 mesh%simfact(ir-2)=2.d0*hh*mesh%radfact(ir-2) 291 ir_last=ir-2 292 end do 293 mesh%simfact(ir_last)=half*mesh%simfact(ir_last) 294 if (mesh%int_meshsz<mesh%mesh_size) mesh%simfact(mesh%int_meshsz+1:mesh%mesh_size)=zero 295 296 mesh%rmax=mesh%rad(mesh%mesh_size) 297 298 end subroutine pawrad_init
m_pawrad/pawrad_isame [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_isame
FUNCTION
Check two radial meshes, returns a logical flag defining whether the meshes have the same equation and an integer flag
INPUTS
Rmesh1,Rmesh2<pawrad_type>=The two radial meshes.
OUTPUT
hasameq=.true. if the two meshes are defined by the same equation.
whichdenser=
* 0 if meshes are not compatible
* 1 if Rmesh1 is denser than Rmesh2
* 2 if Rmesh2 is denser than Rmesh1
SOURCE
483 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 whichdenser =0 ; hasameq=.FALSE. 499 500 if (Rmesh1%mesh_type /= Rmesh2%mesh_type) RETURN 501 502 SELECT CASE (Rmesh1%mesh_type) 503 504 CASE (RMESH_LINEAR) !check for linear meshes 505 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 & .and.Rmesh1%lstep == Rmesh2%lstep ) 513 514 CASE (RMESH_NL) !check for linear meshes 515 hasameq = (Rmesh1%rstep == Rmesh2%rstep) 516 517 CASE DEFAULT 518 msg='Unknown mesh type' 519 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 if (hasameq) then 526 whichdenser= 1 527 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
m_pawrad/pawrad_print [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
pawrad_print
FUNCTION
Reports basic info on the object.
INPUTS
Rmesh<pawrad_type>=Object defining the radial mesh header=String for the header provided by the user. [unit]=Unit number for output, defaults to std_out [prtvol]=Verbosity level, minimal if not specified. [mode_paral]=Either "COLL" or "PERS". Passed to wrtout. Defaults to "COLL"
OUTPUT
Only writing.
SOURCE
394 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 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 412 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 413 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 414 415 msg=ch10//' ==== Info on the Radial Mesh ==== ' 416 if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== ' 417 call wrtout(my_unt,msg,my_mode) 418 419 SELECT CASE (Rmesh%mesh_type) 420 421 CASE (RMESH_LINEAR) 422 write(msg,'(a,i4,a,g12.5)')& 423 & ' - 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 & ' - 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 & ' - 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 & ' - 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 & ' - Non-linear mesh: r(i)=-AA*i/(n-i), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep 440 441 CASE DEFAULT 442 msg = ' Unknown mesh type! Action : check your pseudopotential or input file.' 443 LIBPAW_ERROR(msg) 444 END SELECT 445 446 call wrtout(my_unt,msg,my_mode) 447 448 if (my_prtvol>1) then 449 write(msg,'(a,i4)')' Mesh size for integrals = ',Rmesh%int_meshsz 450 call wrtout(my_unt,msg,my_mode) 451 write(msg,'(a,g12.5)')' rmax=rad(mesh_size) = ',Rmesh%rmax 452 call wrtout(my_unt,msg,my_mode) 453 write(msg,'(a,g12.5)')' Value of stepint = ',Rmesh%stepint 454 call wrtout(my_unt,msg,my_mode) 455 end if 456 457 end subroutine pawrad_print
m_pawrad/pawrad_type [ Types ]
[ Top ] [ m_pawrad ] [ Types ]
NAME
pawrad_type
FUNCTION
For PAW, RADial mesh discretization and related data
SOURCE
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
m_pawrad/poisson [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
poisson
FUNCTION
Solve poisson equation for angularly dependent charge distribution of angular momentum l Densities and potentials are given on a (generalized) radial grid
INPUTS
den(:)= electron density * (4*pi*r**2) appropriate for l ll= l quantum number radmesh <type(pawrad_type)>=data containing radial grid information [screened_sr_separation]= --optional-- separation for screened short-range kernel no screening by default
OUTPUT
[qq]= --optional-- lth moment of the charge; not compatible with screened Coulomb interaction
rv(:)= electrostatic potential * r in (Hartree*Bohr) units
where v(r)=\frac{1}{2l+1}(\frac{int[(r''^(l+2))g(r'')dr'']} {r^(l+1)}
+(r^l) int[r''^(1-l)g(r'')dr''])
SOURCE
1164 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 real(dp),allocatable :: aa(:),bb(:),cc(:),dd(:),radl(:),radl1(:) 1184 1185 ! *************************************************************************** 1186 1187 mesh_size=size(den) 1188 if (size(rv)/=mesh_size.or.mesh_size>radmesh%mesh_size) then 1189 LIBPAW_BUG('wrong sizes!') 1190 end if 1191 1192 use_numerov=(radmesh%mesh_type==1) 1193 1194 use_screened=.false.;sr_omega=zero 1195 if (present(screened_sr_separation)) then 1196 sr_omega=screened_sr_separation 1197 use_screened=(sr_omega>tol8) 1198 ! Numerov method not coded for screened Coulomb interaction 1199 if (use_numerov.and.use_screened) use_numerov=.false. 1200 end if 1201 1202 !============================================== 1203 !UNIFORM GRID - NUMEROV ALGORITHM 1204 !============================================== 1205 if (use_numerov) then 1206 hh=radmesh%rstep 1207 nn=radmesh%int_meshsz-1 1208 LIBPAW_ALLOCATE(aa,(nn)) 1209 LIBPAW_ALLOCATE(bb,(nn)) 1210 LIBPAW_ALLOCATE(cc,(nn+1)) 1211 do ir=1,nn 1212 aa(ir)=two*hh*den(ir+1)/(ir) 1213 bb(ir)=den(ir+1)*((ir*hh)**ll) 1214 end do 1215 cc(1)=zero 1216 cc(2:nn+1)=bb(1:nn) 1217 call simp_gen(qq_,cc,radmesh) 1218 qq_=qq_/dble(2*ll+1) 1219 rv(1)=aa(1)+0.1_dp*aa(2) 1220 do ir=2,nn-1 1221 rv(ir)=aa(ir)+0.1_dp*(aa(ir+1)+aa(ir-1)) 1222 end do 1223 rv(nn)=aa(nn)+0.1_dp*aa(nn-1) 1224 angm=dble(ll*(ll+1)) 1225 rr=(nn+1)*hh 1226 rv(nn)=rv(nn)+(2.4_dp-0.2_dp*angm/((nn+1)**2))*qq_/(rr**ll) 1227 do ir=1,nn 1228 aa(ir)=angm/(ir*ir) 1229 bb(ir)=2.4_dp+aa(ir) 1230 end do 1231 do ir=1,nn-1 1232 cc(ir)=-1.2_dp+0.1_dp*aa(ir+1) 1233 end do 1234 do ir=nn,2,-1 1235 aa(ir)=-1.2_dp+0.1_dp*aa(ir-1) 1236 end do 1237 if (nn.eq.1) then 1238 rv(2)=rv(1)/bb(1) 1239 rv(1)=zero 1240 else 1241 do ir=2,nn 1242 bb(ir)=bb(ir)-aa(ir)*cc(ir-1)/bb(ir-1) 1243 end do 1244 rv(1)=rv(1)/bb(1) 1245 do ir=2,nn 1246 rv(ir)=(rv(ir)-aa(ir)*rv(ir-1))/bb(ir) 1247 end do 1248 do ir=nn-1,1,-1 1249 rv(ir)=rv(ir)-cc(ir)*rv(ir+1)/bb(ir) 1250 end do 1251 do ir=nn+1,2,-1 1252 rv(ir)=rv(ir-1) 1253 end do 1254 rv(1)=zero 1255 if (nn+1<mesh_size) rv(nn+2:mesh_size)=zero 1256 end if 1257 rv(:)=half*rv(:) 1258 if (present(qq)) qq=qq_ 1259 LIBPAW_DEALLOCATE(aa) 1260 LIBPAW_DEALLOCATE(bb) 1261 LIBPAW_DEALLOCATE(cc) 1262 1263 ! ============================================== 1264 ! ANY OTHER GRID - SIMPSON ALGORITHM 1265 ! ============================================== 1266 else 1267 nn=mesh_size;hh=third*radmesh%stepint 1268 do while (abs(den(nn))<tol16.and.nn>radmesh%int_meshsz) 1269 nn=nn-1 1270 end do 1271 mm=nn;if (radmesh%mesh_type==3) mm=mm-1 1272 LIBPAW_ALLOCATE(aa,(nn)) 1273 LIBPAW_ALLOCATE(bb,(nn)) 1274 LIBPAW_ALLOCATE(cc,(nn)) 1275 LIBPAW_ALLOCATE(dd,(nn)) 1276 1277 if (.not.use_screened) then 1278 1279 ! Standard Coulomb integral 1280 LIBPAW_ALLOCATE(radl,(nn)) 1281 LIBPAW_ALLOCATE(radl1,(nn)) 1282 do jr=nn,2,-1 1283 ir=nn-jr+1 1284 radl(jr) =radmesh%rad(jr)**ll 1285 radl1(jr)=radmesh%rad(jr)*radl(jr) 1286 aa(ir)=den(jr)*radmesh%radfact(jr)*radl(jr) 1287 bb(ir)=den(jr)*radmesh%radfact(jr)/radl1(jr) 1288 end do 1289 radl(1)=zero;radl1(1)=zero 1290 ee(2)=aa(nn-1);ee(3)=aa(nn-2);ee(4)=aa(nn-3) 1291 call pawrad_deducer0(ee,4,radmesh) 1292 aa(nn)=ee(1) 1293 ee(2)=bb(nn-1);ee(3)=bb(nn-2);ee(4)=bb(nn-3) 1294 call pawrad_deducer0(ee,4,radmesh) 1295 bb(nn)=ee(1) 1296 cc(1)=zero;dd(1)=zero 1297 do ir=3,mm,2 1298 cc(ir) =cc(ir-2)+hh*(aa(ir-2)+four*aa(ir-1)+aa(ir)) 1299 cc(ir-1)=cc(ir-2)+hh*(1.25_dp*aa(ir-2)+two*aa(ir-1)-quarter*aa(ir)) 1300 dd(ir) =dd(ir-2)+hh*(bb(ir-2)+four*bb(ir-1)+bb(ir)) 1301 dd(ir-1)=dd(ir-2)+hh*(1.25_dp*bb(ir-2)+two*bb(ir-1)-quarter*bb(ir)) 1302 end do 1303 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 cc(mm)=cc(mm-1)+hh*(1.25_dp*aa(mm-2)+two*aa(mm-1)-quarter*aa(mm)) 1307 dd(mm)=dd(mm-1)+hh*(1.25_dp*bb(mm-2)+two*bb(mm-1)-quarter*bb(mm)) 1308 end if 1309 if (mm<nn) then 1310 cc(nn)=cc(mm)+half*(aa(mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1)) 1311 dd(nn)=dd(mm)+half*(bb(mm)+bb(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1)) 1312 end if 1313 rv(1)=zero 1314 do ir=2,nn 1315 jr=nn-ir+1 1316 rv(ir)=(dd(jr)*radl1(ir)+(cc(nn)-cc(jr))/radl(ir))/(two*ll+one) 1317 end do 1318 if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn) 1319 if (present(qq)) qq=cc(nn)/(two*ll+one) 1320 LIBPAW_DEALLOCATE(radl) 1321 LIBPAW_DEALLOCATE(radl1) 1322 else 1323 1324 ! Short-range screened Coulomb integral 1325 rv(1)=zero 1326 do ir=2,nn 1327 ri=sr_omega*radmesh%rad(ir) 1328 do jr=2,nn 1329 rj=sr_omega*radmesh%rad(jr) 1330 aa(jr)=den(jr)*radmesh%radfact(jr)*sr_omega*screened_coul_kernel(ll,ri,rj) 1331 end do 1332 call pawrad_deducer0(aa(1:4),4,radmesh) 1333 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 if (ir>2) then 1337 do jr=ir-2,1,-2 1338 intg=intg+hh*(aa(jr)+four*aa(jr+1)+aa(jr+2)) 1339 end do 1340 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 if (nn >2) intg=intg+hh*(+1.25_dp*aa(1)+two*aa(2)-quarter*aa(3)) 1343 if (nn==2) intg=intg+(aa(1)+aa(2))*hh*half 1344 end if 1345 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 if (ir<nn-2) then 1348 do jr=ir+2,nn,2 1349 intg=intg+hh*(aa(jr-2)+four*aa(jr-1)+aa(jr)) 1350 end do 1351 if (mod(nn-ir,2)==1) then 1352 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 & +106._dp*aa(nn-3)-19._dp*aa(nn-4))/240._dp 1356 else 1357 intg=intg+hh*(1.25_dp*aa(nn)+two*aa(nn-1)-quarter*aa(nn-2)) 1358 endif 1359 end if 1360 else if (ir==nn-1) then 1361 intg=intg+(aa(nn-1)+aa(nn))*hh*half 1362 end if 1363 rv(ir)=intg/(two*ll+one)*radmesh%rad(ir) 1364 end do 1365 if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn) 1366 if (present(qq)) qq=zero ! Not relevant 1367 1368 end if 1369 LIBPAW_DEALLOCATE(aa) 1370 LIBPAW_DEALLOCATE(bb) 1371 LIBPAW_DEALLOCATE(cc) 1372 LIBPAW_DEALLOCATE(dd) 1373 1374 end if 1375 1376 end subroutine poisson
m_pawrad/screened_coul_kernel [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
screened_coul_kernel
FUNCTION
Evaluates the kernel function used to compute the integral of a screened Coulomb potential (with erfc function) See Angyan, Gerber, Marsman, J. Phys. A: Math. Gen. 39, 8613 (2006) [[cite:Angyan2006]]
INPUTS
[formula]=optional; used to force formula (1=full; 2=dev. near r2=0; 3=dev near r1*r2=0) order= order of the function (typically l quantum number) r1,r2=input arguments
OUTPUT
screened_coul_kernel=output radial function
SOURCE
1463 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 if (order>6) then 1479 msg='PAW screened exchange not coded for l>2!' 1480 LIBPAW_ERROR(msg) 1481 end if 1482 1483 !Use max and min of arguments 1484 xx=max(r1,r2) ; yy=min(r1,r2) 1485 1486 !Unscreened Coulomb interaction for very small (or negative) arguments 1487 if (xx<tol8) then 1488 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 formula_=1;if (xx<0.25_dp.and.yy<0.25_dp) formula_=2 1494 if (present(formula)) formula_=formula 1495 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 xy=xx*yy ; sqrtpi=sqrt(pi) 1503 sumexp=exp(-(xx+yy)**2)+exp(-(xx-yy)**2) 1504 difexp=exp(-(xx+yy)**2)-exp(-(xx-yy)**2) 1505 erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy) 1506 if (order==0) then 1507 f0=-difexp/(two*sqrtpi*xy) 1508 hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy) 1509 screened_coul_kernel = hh + f0 1510 else if(order==1) then 1511 f0=-difexp/(two*sqrtpi*xy) 1512 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1513 hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2) 1514 screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy 1515 else if (order==2) then 1516 f0=-difexp/(two*sqrtpi*xy) 1517 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1518 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1519 hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3) 1520 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 f0=-difexp/(two*sqrtpi*xy) 1523 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1524 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1525 f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4) 1526 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 & + f0*(xx**6+yy**6)/xy**3 1529 else if (order==4) then 1530 f0=-difexp/(two*sqrtpi*xy) 1531 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1532 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1533 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 & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5) 1536 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 & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4 1539 else if (order==5) then 1540 f0=-difexp/(two*sqrtpi*xy) 1541 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1542 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1543 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 & +(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 & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6) 1548 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 & + f0*(xx**10+yy**10)/xy**5 1552 else if (order==6) then 1553 f0=-difexp/(two*sqrtpi*xy) 1554 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1555 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1556 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 & +(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 & +(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 & (672._dp*xy**5+10080._dp*xy**3+20790._dp*xy)*sumexp)/(sqrtpi*(two*xy)**7) 1563 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 & + 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 sqrtpi=sqrt(pi) 1574 if (order==0) then 1575 screened_coul_kernel = paw_derfc(xx)/xx + exp(-xx**2)/sqrtpi * & 1576 & (two/three *yy**2 & 1577 & +(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 & +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 & +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 & +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 & +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 & +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 & +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 xy=xx*yy ; sqrtpi=sqrt(pi) 1616 sumexp=exp(-(xx**2+yy**2)) 1617 erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy) 1618 if (order==0) then 1619 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1620 hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy) 1621 screened_coul_kernel = hh + f0 1622 else if(order==1) then 1623 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1624 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1625 hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2) 1626 screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy 1627 else if (order==2) then 1628 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1629 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1630 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1631 hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3) 1632 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 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1635 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1636 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1637 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1638 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 & + f0*(xx**6+yy**6)/xy**3 1641 else if (order==4) then 1642 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1643 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1644 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1645 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1646 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1647 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 & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4 1650 else if (order==5) then 1651 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1652 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1653 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1654 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1655 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1656 f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp 1657 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 & + f0*(xx**10+yy**10)/xy**5 1661 else if (order==6) then 1662 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1663 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1664 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1665 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1666 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1667 f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp 1668 f6 = 128._dp*xy**2*(15._dp+two*xy**2)/(1027025._dp*sqrtpi)*sumexp 1669 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 & + 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
m_pawrad/simp_gen [ Functions ]
[ Top ] [ m_pawrad ] [ Functions ]
NAME
simp_gen
FUNCTION
Do integral on a given (generalized) grid using Simpson rule
INPUTS
radmesh <type(pawrad_type)>=data containing radial grid information func(:)=integrand values r_for_intg=upper boundary for (future) integration over the radial grid (can be negative for an integration over the whole grid)
OUTPUT
intg=resulting integral by Simpson rule
NOTES
Possible mesh types (radmesh%mesh_type) mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
SOURCE
815 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 real(dp),allocatable :: simfact(:) 834 character(len=500) :: msg 835 836 ! ************************************************************************* 837 838 if (present(r_for_intg)) then 839 if (r_for_intg>0.d0) then 840 ir=min(pawrad_ifromr(radmesh,r_for_intg),radmesh%mesh_size) 841 if (ir<radmesh%mesh_size) then 842 if (abs(radmesh%rad(ir+1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir+1 843 end if 844 if (ir>1) then 845 if (abs(radmesh%rad(ir-1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir-1 846 end if 847 int_meshsz=ir 848 else 849 int_meshsz=radmesh%mesh_size 850 end if 851 if (int_meshsz>radmesh%mesh_size.or.int_meshsz>size(func)) then 852 write(msg,'(3(a,i4))')"int_meshsz= ",int_meshsz," > mesh_size=",radmesh%mesh_size,& 853 & ", size(func)=",size(func) 854 LIBPAW_BUG(msg) 855 end if 856 isim=3; if (radmesh%mesh_type==3)isim=4 857 LIBPAW_ALLOCATE(simfact,(radmesh%mesh_size)) 858 hh=radmesh%stepint/3.d0 859 simfact(int_meshsz)=hh*radmesh%radfact(int_meshsz) 860 simfact(1:isim-2)=zero 861 ir_last=1 862 do ir=int_meshsz,isim,-2 863 simfact(ir-1)=4.d0*hh*radmesh%radfact(ir-1) 864 simfact(ir-2)=2.d0*hh*radmesh%radfact(ir-2) 865 ir_last=ir-2 866 end do 867 simfact(ir_last)=half*simfact(ir_last) 868 if (int_meshsz<radmesh%mesh_size) simfact(int_meshsz+1:radmesh%mesh_size)=zero 869 870 nn=int_meshsz 871 simp=zero 872 do ii=1,nn 873 simp=simp+func(ii)*simfact(ii) 874 end do 875 LIBPAW_DEALLOCATE(simfact) 876 877 else 878 if (radmesh%int_meshsz>size(func)) then 879 write(msg,'(2(a,i4))')"int_meshsz= ",int_meshsz," > size(func)=",size(func) 880 LIBPAW_BUG(msg) 881 end if 882 nn=radmesh%int_meshsz 883 simp=zero 884 do ii=1,nn 885 simp=simp+func(ii)*radmesh%simfact(ii) 886 end do 887 end if 888 889 resid=zero 890 if (radmesh%mesh_type==3) then 891 resid=half*(func(2)+func(1))*(radmesh%rad(2)-radmesh%rad(1)) 892 if (mod(nn,2)==1) resid=resid+radmesh%stepint/3.d0*(1.25d0*func(2)*radmesh%radfact(2) & 893 & +2.d0*func(3)*radmesh%radfact(3)-0.25d0*func(4)*radmesh%radfact(4)) 894 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 & -0.25d0*func(3)*radmesh%radfact(3)) 897 end if 898 899 intg=simp+resid 900 901 end subroutine simp_gen