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-2018 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
28 #include "libpaw.h" 29 30 MODULE m_pawrad 31 32 USE_DEFS 33 USE_MSG_HANDLING 34 USE_MPI_WRAPPERS 35 USE_MEMORY_PROFILING 36 37 use m_paw_numeric, only : paw_derfc 38 39 implicit none 40 41 private 42 43 !public procedures. 44 public :: pawrad_init ! Main creation method 45 public :: pawrad_free ! Free the allocated memory 46 public :: pawrad_print ! Printout of the basic info 47 public :: pawrad_isame ! Checks whether two meshes are equivalent or have the same equation. 48 public :: pawrad_copy ! Returns a copy of the mesh. 49 public :: pawrad_ifromr ! Retrieve the Index FROM a given R value in a radial grid. 50 public :: pawrad_deducer0 ! Extrapolate r=0 value of a function from values near r=0. 51 public :: pawrad_bcast ! Broadcast pawrad datastructure over a given MPI communicator 52 public :: simp_gen ! Performs integral on a given (generalized) grid using Simpson rule. 53 public :: nderiv_gen ! Do corrected first (and 2nd) derivation on a given (generalized) grid. 54 public :: nderiv_lin ! Do corrected first (and 2nd) derivation on a given linear grid. 55 public :: bound_deriv ! Computes derivatives at boundaries of the mesh 56 public :: poisson ! Solves Poisson eq. for angularly dependent charge distribution of angular momentum l 57 public :: screened_coul_kernel ! Kernel used to compute short-range screened Coulomb integrals 58 public :: calc_slatradl ! Calculates the radial part of Slater integrals. 59 60 interface pawrad_free 61 module procedure pawrad_free_0D 62 module procedure pawrad_free_1D 63 end interface pawrad_free 64 65 ! TODO: Might use bit flags, but all radmesh stuff should be encapsulated here! 66 integer,private,parameter :: RMESH_LINEAR = 1 67 integer,private,parameter :: RMESH_LOG1 = 2 68 integer,private,parameter :: RMESH_LOG2 = 3 69 integer,private,parameter :: RMESH_LOG3 = 4 70 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)
PARENTS
m_paw_atom,m_pawpsp,m_pawxmlps,outscfcv
CHILDREN
poisson,simp_gen
SOURCE
1292 subroutine bound_deriv(func,mesh,nn,yp1,ypn) 1293 1294 1295 !This section has been created automatically by the script Abilint (TD). 1296 !Do not modify the following lines by hand. 1297 #undef ABI_FUNC 1298 #define ABI_FUNC 'bound_deriv' 1299 !End of the abilint section 1300 1301 implicit none 1302 1303 !Arguments---------------------- 1304 integer, intent(in) :: nn 1305 real(dp), intent(in) :: func(nn) 1306 real(dp), intent(out) :: yp1,ypn 1307 type(pawrad_type),intent(in) :: mesh 1308 1309 !************************************************************************* 1310 1311 if (mesh%radfact(1)>zero) then 1312 yp1=1._dp/12._dp/mesh%stepint/mesh%radfact(1) & 1313 & *(-25._dp*func(1)+48._dp*func(2)-36._dp*func(3)+16._dp*func(4)-3._dp*func(5)) 1314 else 1315 yp1=(func(2)-func(1))/(mesh%rad(2)-mesh%rad(1)) 1316 end if 1317 ypn=1._dp/12._dp/mesh%stepint & 1318 & *( 3._dp*func(nn-4)-16._dp*func(nn-3)+36._dp*func(nn-2)-48._dp*func(nn-1) & 1319 & +25._dp*func(nn))/mesh%radfact(nn) 1320 1321 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)$.
PARENTS
m_paw_slater
CHILDREN
poisson,simp_gen
SOURCE
1936 subroutine calc_slatradl(ll,mesh_size,ff1,ff2,Pawrad,integral) 1937 1938 1939 !This section has been created automatically by the script Abilint (TD). 1940 !Do not modify the following lines by hand. 1941 #undef ABI_FUNC 1942 #define ABI_FUNC 'calc_slatradl' 1943 !End of the abilint section 1944 1945 implicit none 1946 1947 !scalars 1948 integer,intent(in) :: mesh_size,ll 1949 real(dp),intent(out) :: integral 1950 !arrays 1951 real(dp),intent(in) :: ff1(mesh_size),ff2(mesh_size) 1952 type(pawrad_type),intent(in) :: Pawrad 1953 1954 !Local variables --------------------------------------- 1955 !scalars 1956 integer :: int_meshsz 1957 character(len=100) :: msg 1958 !arrays 1959 real(dp),allocatable :: hh(:),gg(:) 1960 1961 ! ************************************************************************* 1962 1963 if (mesh_size > Pawrad%mesh_size) then 1964 msg='mesh_size > pawrad%mesh_size!' 1965 MSG_BUG(msg) 1966 end if 1967 1968 LIBPAW_ALLOCATE(hh,(mesh_size)) 1969 LIBPAW_ALLOCATE(gg,(mesh_size)) 1970 !hh = zero 1971 !gg = zero 1972 1973 int_meshsz=Pawrad%int_meshsz 1974 !$int_meshsz=Pawrad%mesh_size 1975 1976 ! the line below requires hh as work array. 1977 hh = ff2 1978 1979 ! TODO find where int_meshsz is calculated and if it can affects the results. 1980 if (int_meshsz<mesh_size) hh(int_meshsz+1:mesh_size)=zero 1981 1982 call poisson(hh,ll,Pawrad,gg) 1983 1984 gg(2:mesh_size) = gg(2:mesh_size)/Pawrad%rad(2:mesh_size) 1985 1986 hh(1) = zero 1987 hh(2:mesh_size)= ff1(2:mesh_size) * gg(2:mesh_size) 1988 LIBPAW_DEALLOCATE(gg) 1989 1990 call simp_gen(integral,hh,Pawrad) 1991 integral = four_pi * integral 1992 1993 LIBPAW_DEALLOCATE(hh) 1994 1995 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
PARENTS
m_paw_pwaves_lmn,m_pawdij,m_pawpsp,m_pawxc,optics_paw,optics_paw_core pawinit,pawnabla_init,poslifetime,posratecore,spline_paw_fncs
CHILDREN
poisson,simp_gen
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
1079 subroutine nderiv_gen(der,func,radmesh,der2) 1080 1081 1082 !This section has been created automatically by the script Abilint (TD). 1083 !Do not modify the following lines by hand. 1084 #undef ABI_FUNC 1085 #define ABI_FUNC 'nderiv_gen' 1086 !End of the abilint section 1087 1088 implicit none 1089 1090 !Arguments ------------------------------------ 1091 !scalars 1092 type(pawrad_type),intent(in) :: radmesh 1093 !arrays 1094 real(dp),intent(in) :: func(:) 1095 real(dp),intent(out) :: der(:) 1096 real(dp),optional,intent(out) :: der2(:) 1097 1098 !Local variables------------------------------- 1099 !scalars 1100 integer :: msz 1101 logical :: compute_2der 1102 character(len=500) :: msg 1103 1104 ! ************************************************************************* 1105 1106 msz=size(func) 1107 if (size(der)/=msz.or.msz>radmesh%mesh_size) then 1108 msg='wrong sizes for in/out arrays!' 1109 MSG_BUG(msg) 1110 end if 1111 1112 compute_2der=(present(der2)) 1113 1114 if (radmesh%mesh_type==1) then 1115 1116 call nderiv_lin(radmesh%rstep,func,der,msz,1) 1117 if (compute_2der) then 1118 call nderiv_lin(radmesh%rstep,func,der2,msz,2) 1119 end if 1120 1121 else if (radmesh%mesh_type==2) then 1122 1123 call nderiv_lin(radmesh%lstep,func,der,msz,1) 1124 der(1:msz)=der(1:msz)/radmesh%radfact(1:msz) 1125 if (compute_2der)then 1126 call nderiv_lin(radmesh%lstep,func,der2,msz,2) 1127 der2(1:msz)=(der2(1:msz)/radmesh%radfact(1:msz)-der(1:msz))/radmesh%radfact(1:msz) 1128 end if 1129 1130 else if (radmesh%mesh_type==3) then 1131 1132 call nderiv_lin(radmesh%lstep,func(2:msz),der(2:msz),msz-1,1) 1133 der(2:msz)=der(2:msz)/radmesh%radfact(2:msz) 1134 call pawrad_deducer0(der,msz,radmesh) 1135 if (compute_2der)then 1136 call nderiv_lin(radmesh%lstep,func(2:msz),der2(2:msz),msz-1,2) 1137 der2(2:msz)=(der2(2:msz)/radmesh%radfact(2:msz)-der(2:msz))/radmesh%radfact(2:msz) 1138 call pawrad_deducer0(der2,msz,radmesh) 1139 end if 1140 1141 else if (radmesh%mesh_type==4) then 1142 1143 call nderiv_lin(radmesh%lstep,func,der,msz,1) 1144 der(1:msz)=der(1:msz)/radmesh%radfact(1:msz) 1145 if (compute_2der)then 1146 call nderiv_lin(radmesh%lstep,func,der2,msz,2) 1147 der2(1:msz)=der2(1:msz)/radmesh%radfact(1:msz)**2-der(1:msz)/radmesh%rstep 1148 end if 1149 1150 else if (radmesh%mesh_type==5) then 1151 1152 call nderiv_lin(one,func,der,msz,1) 1153 der(1:msz)=der(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep) 1154 if (compute_2der)then 1155 call nderiv_lin(one,func,der2,msz,2) 1156 der2(1:msz)=der2(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep)**2-two*der(1:msz)/& 1157 & (radmesh%rstep+radmesh%rad(1:msz)) 1158 end if 1159 end if 1160 1161 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
PARENTS
m_pawrad
CHILDREN
poisson,simp_gen
SOURCE
1190 subroutine nderiv_lin(hh,yy,zz,ndim,norder) 1191 1192 1193 !This section has been created automatically by the script Abilint (TD). 1194 !Do not modify the following lines by hand. 1195 #undef ABI_FUNC 1196 #define ABI_FUNC 'nderiv_lin' 1197 !End of the abilint section 1198 1199 implicit none 1200 1201 !Arguments --------------------------------------------- 1202 !scalars 1203 integer,intent(in) :: ndim,norder 1204 real(dp),intent(in) :: hh 1205 !arrays 1206 real(dp),intent(in) :: yy(ndim) 1207 real(dp),intent(out) :: zz(ndim) 1208 1209 !Local variables --------------------------------------- 1210 !scalars 1211 integer :: ier,ii 1212 real(dp) :: aa,bb,cc,h1,y1 1213 1214 ! ************************************************************************* 1215 1216 !Initialization (common to 1st and 2nd derivative) 1217 h1=one/(12.d0*hh) 1218 y1=yy(ndim-4) 1219 1220 !FIRST DERIVATIVE 1221 !================ 1222 if (norder==1) then 1223 1224 ! Prepare differentiation loop 1225 bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5)) 1226 cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5)) 1227 ! Start differentiation loop 1228 do ii=5,ndim 1229 aa=bb;bb=cc 1230 cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3))) 1231 zz(ii-4)=aa 1232 end do 1233 ! Normal exit 1234 ier=0 1235 aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim)) 1236 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)) 1237 zz(ndim-1)=aa 1238 zz(ndim-2)=cc 1239 zz(ndim-3)=bb 1240 1241 ! SECOND DERIVATIVE 1242 ! ================= 1243 else 1244 h1=h1/hh 1245 ! Prepare differentiation loop 1246 bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5)) 1247 cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5)) 1248 ! Start differentiation loop 1249 do ii=5,ndim 1250 aa=bb;bb=cc 1251 cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2)) 1252 zz(ii-4)=aa 1253 end do 1254 ! Normal exit 1255 ier=0 1256 aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim)) 1257 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)) 1258 zz(ndim-1)=aa 1259 zz(ndim-2)=cc 1260 zz(ndim-3)=bb 1261 1262 end if !norder 1263 1264 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
PARENTS
m_pawpsp
CHILDREN
poisson,simp_gen
SOURCE
762 subroutine pawrad_bcast(pawrad,comm_mpi) 763 764 765 !This section has been created automatically by the script Abilint (TD). 766 !Do not modify the following lines by hand. 767 #undef ABI_FUNC 768 #define ABI_FUNC 'pawrad_bcast' 769 !End of the abilint section 770 771 implicit none 772 773 !Arguments ------------------------------------ 774 !scalars 775 integer,intent(in) :: comm_mpi 776 type(pawrad_type),intent(inout) :: pawrad 777 778 !Local variables------------------------------- 779 !scalars 780 integer :: ierr,indx,me,nn,isz1 781 integer :: if_rad,if_radfact,if_simfact !flags used to communicate 782 character(len=500) :: msg 783 !arrays 784 integer,allocatable :: list_int(:) 785 real(dp),allocatable :: list_dpr(:) 786 787 !************************************************************************* 788 789 me=xmpi_comm_rank(comm_mpi) 790 791 !Initializations 792 if_rad=0; if_radfact=0; if_simfact=0 793 794 !calculate the size of the reals 795 if(me==0) then 796 if (allocated(pawrad%rad)) then 797 if_rad=1 !communicate rad 798 isz1=size(pawrad%rad) 799 if(isz1/=pawrad%mesh_size) then 800 msg='rad: sz1 /= pawrad%mesh_size (1)' 801 MSG_BUG(msg) 802 end if 803 end if 804 if (allocated(pawrad%radfact)) then 805 if_radfact=1 !communicate radfact 806 isz1=size(pawrad%radfact) 807 if(isz1/=pawrad%mesh_size) then 808 msg='radfact: sz1 /= pawrad%mesh_size (2)' 809 MSG_BUG(msg) 810 end if 811 end if 812 if (allocated(pawrad%simfact)) then 813 if_simfact=1 !communicate simfact 814 isz1=size(pawrad%simfact) 815 if(isz1/=pawrad%mesh_size) then 816 msg='simfact: sz1 /= pawrad%mesh_size (3)' 817 MSG_BUG(msg) 818 end if 819 end if 820 end if 821 822 !Brodcast the integers 823 LIBPAW_ALLOCATE(list_int,(6)) 824 if(me==0) then 825 list_int(1)=pawrad%int_meshsz 826 list_int(2)=pawrad%mesh_size 827 list_int(3)=pawrad%mesh_type 828 list_int(4)=if_rad 829 list_int(5)=if_radfact 830 list_int(6)=if_simfact 831 end if 832 call xmpi_bcast(list_int,0,comm_mpi,ierr) 833 if(me/=0) then 834 pawrad%int_meshsz =list_int(1) 835 pawrad%mesh_size =list_int(2) 836 pawrad%mesh_type =list_int(3) 837 if_rad=list_int(4) 838 if_radfact=list_int(5) 839 if_simfact=list_int(6) 840 end if 841 LIBPAW_DEALLOCATE(list_int) 842 843 !Broadcast the reals 844 nn=4+pawrad%mesh_size*(if_rad+if_radfact+if_simfact) 845 LIBPAW_ALLOCATE(list_dpr,(nn)) 846 if(me==0) then 847 list_dpr(1)=pawrad%lstep 848 list_dpr(2)=pawrad%rmax 849 list_dpr(3)=pawrad%rstep 850 list_dpr(4)=pawrad%stepint 851 indx=5 852 if (if_rad==1) then 853 isz1=pawrad%mesh_size 854 list_dpr(indx:indx+isz1-1)=pawrad%rad(1:isz1) 855 indx=indx+isz1 856 end if 857 if (if_radfact==1) then 858 isz1=pawrad%mesh_size 859 list_dpr(indx:indx+isz1-1)=pawrad%radfact(1:isz1) 860 indx=indx+isz1 861 end if 862 if (if_simfact==1) then 863 isz1=pawrad%mesh_size 864 list_dpr(indx:indx+isz1-1)=pawrad%simfact(1:isz1) 865 indx=indx+isz1 866 end if 867 end if 868 call xmpi_bcast(list_dpr,0,comm_mpi,ierr) 869 if(me/=0) then 870 pawrad%lstep=list_dpr(1) 871 pawrad%rmax=list_dpr(2) 872 pawrad%rstep=list_dpr(3) 873 pawrad%stepint=list_dpr(4) 874 indx=5 875 ! Deallocate all arrays: 876 if (allocated(pawrad%rad)) then 877 LIBPAW_DEALLOCATE(pawrad%rad) 878 end if 879 if (allocated(pawrad%radfact)) then 880 LIBPAW_DEALLOCATE(pawrad%radfact) 881 end if 882 if (allocated(pawrad%simfact)) then 883 LIBPAW_DEALLOCATE(pawrad%simfact) 884 end if 885 ! Communicate if flag is set to 1: 886 if(if_rad==1) then 887 isz1=pawrad%mesh_size 888 LIBPAW_ALLOCATE(pawrad%rad,(isz1)) 889 pawrad%rad(1:isz1)=list_dpr(indx:indx+isz1-1) 890 indx=indx+isz1 891 end if 892 if(if_radfact==1) then 893 isz1=pawrad%mesh_size 894 LIBPAW_ALLOCATE(pawrad%radfact,(isz1)) 895 pawrad%radfact(1:isz1)=list_dpr(indx:indx+isz1-1) 896 indx=indx+isz1 897 end if 898 if(if_simfact==1) then 899 isz1=pawrad%mesh_size 900 LIBPAW_ALLOCATE(pawrad%simfact,(isz1)) 901 pawrad%simfact(1:isz1)=list_dpr(indx:indx+isz1-1) 902 indx=indx+isz1 903 end if 904 end if 905 LIBPAW_DEALLOCATE(list_dpr) 906 907 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
PARENTS
m_pawpsp,m_pawpwij
CHILDREN
poisson,simp_gen
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
641 subroutine pawrad_copy(mesh1,mesh2) 642 643 644 !This section has been created automatically by the script Abilint (TD). 645 !Do not modify the following lines by hand. 646 #undef ABI_FUNC 647 #define ABI_FUNC 'pawrad_copy' 648 !End of the abilint section 649 650 implicit none 651 652 !Arguments ------------------------------------ 653 !scalars 654 type(pawrad_type),intent(in) :: mesh1 655 type(pawrad_type),intent(out) :: mesh2 656 657 !Local variables------------------------------- 658 !scalars 659 integer :: ir 660 661 ! ************************************************************************* 662 663 mesh2%mesh_type =mesh1%mesh_type 664 mesh2%mesh_size =mesh1%mesh_size 665 mesh2%int_meshsz=mesh1%int_meshsz 666 mesh2%lstep =mesh1%lstep 667 mesh2%rstep =mesh1%rstep 668 mesh2%stepint =mesh1%stepint 669 mesh2%rmax =mesh1%rmax 670 671 LIBPAW_ALLOCATE(mesh2%rad,(mesh1%mesh_size)) 672 LIBPAW_ALLOCATE(mesh2%radfact,(mesh1%mesh_size)) 673 LIBPAW_ALLOCATE(mesh2%simfact,(mesh1%mesh_size)) 674 do ir=1,mesh1%mesh_size 675 mesh2%rad(ir) =mesh1%rad(ir) 676 mesh2%radfact(ir)=mesh1%radfact(ir) 677 mesh2%simfact(ir)=mesh1%simfact(ir) 678 end do 679 680 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
PARENTS
denfgr,m_paw_atom,m_paw_gaussfit,m_paw_pwaves_lmn,m_paw_slater,m_pawdij m_pawpsp,m_pawrad,m_pawxc,make_efg_onsite,optics_paw,optics_paw_core pawdenpot,pawdensities,pawnabla_init
CHILDREN
poisson,simp_gen
SOURCE
710 subroutine pawrad_deducer0(func,funcsz,radmesh) 711 712 713 !This section has been created automatically by the script Abilint (TD). 714 !Do not modify the following lines by hand. 715 #undef ABI_FUNC 716 #define ABI_FUNC 'pawrad_deducer0' 717 !End of the abilint section 718 719 implicit none 720 721 !Arguments ------------------------------------ 722 !scalars 723 integer,intent(in) :: funcsz 724 type(pawrad_type),intent(in) :: radmesh 725 !arrays 726 real(dp),intent(inout) :: func(funcsz) 727 728 ! ************************************************************************* 729 730 if (radmesh%mesh_type==1.or.radmesh%mesh_type==2.or.radmesh%mesh_type==4.or.radmesh%mesh_type==5) then 731 func(1)=func(4)+3*(func(2)-func(3)) 732 else if (radmesh%mesh_type==3) then 733 func(1)=func(4)+exp(two*radmesh%lstep)/(exp(radmesh%lstep)-one)*(func(2)-func(3)) 734 end if 735 736 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
PARENTS
m_pawrad
CHILDREN
poisson,simp_gen
SOURCE
335 subroutine pawrad_free_0D(Rmesh) 336 337 338 !This section has been created automatically by the script Abilint (TD). 339 !Do not modify the following lines by hand. 340 #undef ABI_FUNC 341 #define ABI_FUNC 'pawrad_free_0D' 342 !End of the abilint section 343 344 implicit none 345 346 !Arguments ------------------------------------ 347 !arrays 348 type(pawrad_type),intent(inout) :: Rmesh 349 350 !Local variables------------------------------- 351 352 ! ************************************************************************* 353 354 !@Pawrad_type 355 356 if (allocated(Rmesh%rad )) then 357 LIBPAW_DEALLOCATE(Rmesh%rad) 358 end if 359 if (allocated(Rmesh%radfact)) then 360 LIBPAW_DEALLOCATE(Rmesh%radfact) 361 end if 362 if (allocated(Rmesh%simfact)) then 363 LIBPAW_DEALLOCATE(Rmesh%simfact) 364 end if 365 Rmesh%int_meshsz=0 366 Rmesh%mesh_size=0 367 Rmesh%mesh_type=-1 368 369 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
PARENTS
CHILDREN
poisson,simp_gen
SOURCE
388 subroutine pawrad_free_1D(Rmesh) 389 390 391 !This section has been created automatically by the script Abilint (TD). 392 !Do not modify the following lines by hand. 393 #undef ABI_FUNC 394 #define ABI_FUNC 'pawrad_free_1D' 395 !End of the abilint section 396 397 implicit none 398 399 !Arguments ------------------------------------ 400 type(pawrad_type),intent(inout) :: Rmesh(:) 401 402 !Local variables------------------------------- 403 integer :: ii,nn 404 405 ! ************************************************************************* 406 407 !@pawrad_type 408 409 nn=size(Rmesh) 410 if (nn==0) return 411 412 do ii=1,nn 413 call pawrad_free_0D(Rmesh(ii)) 414 end do 415 416 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
PARENTS
CHILDREN
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
1612 function pawrad_ifromr(radmesh,rr) 1613 1614 1615 !This section has been created automatically by the script Abilint (TD). 1616 !Do not modify the following lines by hand. 1617 #undef ABI_FUNC 1618 #define ABI_FUNC 'pawrad_ifromr' 1619 !End of the abilint section 1620 1621 implicit none 1622 1623 !Arguments ------------------------------------ 1624 !scalars 1625 integer :: pawrad_ifromr 1626 real(dp),intent(in) :: rr 1627 type(pawrad_type),intent(in) :: radmesh 1628 1629 !Local variables------------------------------- 1630 character(len=500) :: msg 1631 1632 ! ************************************************************************* 1633 1634 if (radmesh%mesh_type==1) then 1635 pawrad_ifromr=int(tol8+rr/radmesh%rstep)+1 1636 else if (radmesh%mesh_type==2) then 1637 pawrad_ifromr=int(tol8+log(1.d0+rr/radmesh%rstep)/radmesh%lstep)+1 1638 else if (radmesh%mesh_type==3) then 1639 if (rr<radmesh%rstep) then 1640 pawrad_ifromr=1 1641 else 1642 pawrad_ifromr=int(tol8+log(rr/radmesh%rstep)/radmesh%lstep)+2 1643 end if 1644 else if (radmesh%mesh_type==4) then 1645 pawrad_ifromr=int(tol8+(1.d0-exp(-rr/radmesh%rstep))/radmesh%lstep)+1 1646 else if (radmesh%mesh_type==5) then 1647 pawrad_ifromr=int(tol8+(radmesh%lstep*rr)/(radmesh%rstep+rr))+1 1648 else 1649 ! Other values of mesh_type are not allowed (see psp7in.F90) 1650 write(msg,'(a,i0)')" Unknown value of %mesh_type ",radmesh%mesh_type 1651 MSG_ERROR(msg) 1652 end if 1653 1654 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)
PARENTS
dfpt_eltfrxc,m_atom,m_paw_gaussfit,m_pawpsp,m_pawpwij,m_pawxmlps,m_psps mkcore,mkcore_paw,mkcore_wvl,psp8in,psp9in,wvl_initro
CHILDREN
poisson,simp_gen
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
198 subroutine pawrad_init(mesh,mesh_size,mesh_type,rstep,lstep,r_for_intg) 199 200 201 !This section has been created automatically by the script Abilint (TD). 202 !Do not modify the following lines by hand. 203 #undef ABI_FUNC 204 #define ABI_FUNC 'pawrad_init' 205 !End of the abilint section 206 207 implicit none 208 209 !Arguments ------------------------------------ 210 !scalars 211 integer,intent(in),optional :: mesh_size,mesh_type 212 real(dp),intent(in),optional :: rstep,lstep 213 real(dp),intent(in),optional :: r_for_intg 214 type(pawrad_type),intent(inout) :: mesh 215 216 !Local variables------------------------------- 217 !scalars 218 integer :: ir,ir_last,isim,mesh_size_,mesh_type_ 219 real(dp) :: hh,lstep_,rstep_,r_for_intg_ 220 character(len=500) :: msg 221 222 ! ************************************************************************* 223 224 !@pawrad_type 225 226 !Retrieve mesh data 227 mesh_size_ = mesh%mesh_size ; if (present(mesh_size)) mesh_size_ = mesh_size 228 mesh_type_ = mesh%mesh_type ; if (present(mesh_type)) mesh_type_ = mesh_type 229 rstep_ = mesh%rstep ; if (present(rstep)) rstep_ = rstep 230 lstep_ = mesh%lstep ; if (present(lstep)) lstep_ = lstep 231 232 r_for_intg_=-1._dp;if (present(r_for_intg)) r_for_intg_=r_for_intg 233 234 mesh%mesh_size = mesh_size_ 235 mesh%mesh_type = mesh_type_ 236 mesh%rstep = rstep_ 237 mesh%lstep = lstep_ 238 LIBPAW_ALLOCATE(mesh%rad ,(mesh%mesh_size)) 239 LIBPAW_ALLOCATE(mesh%radfact,(mesh%mesh_size)) 240 LIBPAW_ALLOCATE(mesh%simfact,(mesh%mesh_size)) 241 mesh%simfact=zero 242 if (mesh%mesh_type==1) then 243 isim=3 244 mesh%stepint=mesh%rstep 245 mesh%rad(1)=zero;mesh%radfact(1)=one 246 do ir=2,mesh%mesh_size 247 mesh%rad(ir) =mesh%rstep*dble(ir-1) 248 mesh%radfact(ir)=one 249 end do 250 else if (mesh%mesh_type==2) then 251 isim=3 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*(exp(mesh%lstep*dble(ir-1))-one) 256 mesh%radfact(ir)=mesh%rad(ir)+mesh%rstep 257 end do 258 else if (mesh%mesh_type==3) then 259 isim=4 260 mesh%stepint=mesh%lstep 261 mesh%rad(1)=zero;mesh%radfact(1)=zero 262 do ir=2,mesh%mesh_size 263 mesh%rad(ir) =mesh%rstep*exp(mesh%lstep*dble(ir-2)) 264 mesh%radfact(ir)=mesh%rad(ir) 265 end do 266 else if (mesh%mesh_type==4) then 267 isim=3 268 mesh%lstep=one/dble(mesh%mesh_size) 269 mesh%stepint=mesh%lstep 270 mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep 271 do ir=2,mesh%mesh_size 272 mesh%rad(ir) =-mesh%rstep*log(one-mesh%lstep*dble(ir-1)) 273 mesh%radfact(ir)=mesh%rstep/(one-mesh%lstep*dble(ir-1)) 274 end do 275 else if (mesh%mesh_type==5) then 276 isim=3 277 mesh%stepint=mesh%rstep 278 mesh%rad(1)=zero;mesh%radfact(1)=1/mesh%lstep 279 do ir=2,mesh%mesh_size 280 mesh%rad(ir) =mesh%rstep*dble(ir-1)/(mesh%lstep-dble(ir-1)) 281 mesh%radfact(ir)=(mesh%rstep+mesh%rad(ir))/(mesh%lstep-dble(ir-1))/mesh%rstep 282 end do 283 284 else ! Other values of mesh_type are not allowed (see psp7in.F90) 285 write(msg,'(a,i0)')" Unknown value of mesh_type: ",mesh%mesh_type 286 MSG_ERROR(msg) 287 end if 288 289 mesh%int_meshsz=mesh%mesh_size 290 if (r_for_intg_>0.d0) then 291 ir=min(pawrad_ifromr(mesh,r_for_intg_),mesh%mesh_size) 292 if (ir<mesh%mesh_size) then 293 if (abs(mesh%rad(ir+1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir+1 294 end if 295 if (ir>1) then 296 if (abs(mesh%rad(ir-1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir-1 297 end if 298 mesh%int_meshsz=ir 299 end if 300 301 hh=mesh%stepint/3.d0 302 mesh%simfact(mesh%int_meshsz)=hh*mesh%radfact(mesh%int_meshsz) 303 mesh%simfact(1:isim-2)=zero 304 ir_last=1 305 do ir=mesh%int_meshsz,isim,-2 306 mesh%simfact(ir-1)=4.d0*hh*mesh%radfact(ir-1) 307 mesh%simfact(ir-2)=2.d0*hh*mesh%radfact(ir-2) 308 ir_last=ir-2 309 end do 310 mesh%simfact(ir_last)=half*mesh%simfact(ir_last) 311 if (mesh%int_meshsz<mesh%mesh_size) mesh%simfact(mesh%int_meshsz+1:mesh%mesh_size)=zero 312 313 mesh%rmax=mesh%rad(mesh%mesh_size) 314 315 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
PARENTS
m_atom,m_paw_slater
CHILDREN
poisson,simp_gen
SOURCE
550 subroutine pawrad_isame(Rmesh1,Rmesh2,hasameq,whichdenser) 551 552 553 !This section has been created automatically by the script Abilint (TD). 554 !Do not modify the following lines by hand. 555 #undef ABI_FUNC 556 #define ABI_FUNC 'pawrad_isame' 557 !End of the abilint section 558 559 implicit none 560 561 !Arguments ------------------------------------ 562 integer,intent(out) :: whichdenser 563 logical,intent(out) :: hasameq 564 type(pawrad_type),intent(in) :: Rmesh1 565 type(pawrad_type),intent(in) :: Rmesh2 566 567 !Local variables------------------------------- 568 character(len=50) :: msg 569 570 ! ************************************************************************* 571 572 !@pawrad_type 573 574 whichdenser =0 ; hasameq=.FALSE. 575 576 if (Rmesh1%mesh_type /= Rmesh2%mesh_type) RETURN 577 578 SELECT CASE (Rmesh1%mesh_type) 579 580 CASE (RMESH_LINEAR) !check for linear meshes 581 hasameq = (Rmesh1%rstep == Rmesh2%rstep) 582 583 CASE (RMESH_LOG1,& !check for logarithmic meshes 584 & RMESH_LOG2,& 585 & RMESH_LOG3) 586 587 hasameq = ( Rmesh1%rstep == Rmesh2%rstep & 588 & .and.Rmesh1%lstep == Rmesh2%lstep ) 589 590 CASE (RMESH_NL) !check for linear meshes 591 hasameq = (Rmesh1%rstep == Rmesh2%rstep) 592 593 CASE DEFAULT 594 msg='Unknown mesh type' 595 MSG_ERROR(msg) 596 597 END SELECT 598 599 ! === If meshes have same equation, check whether they are equal === 600 ! * Note that also int_meshsz must be equal 601 if (hasameq) then 602 whichdenser= 1 603 if (Rmesh2%mesh_size > Rmesh1%mesh_size ) whichdenser = 2 604 !if (Rmesh1%int_meshsz == Rmesh2%int_meshsz) whichdenser = 2 605 end if 606 607 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.
PARENTS
m_atom
CHILDREN
poisson,simp_gen
SOURCE
446 subroutine pawrad_print(Rmesh,header,unit,prtvol,mode_paral) 447 448 449 !This section has been created automatically by the script Abilint (TD). 450 !Do not modify the following lines by hand. 451 #undef ABI_FUNC 452 #define ABI_FUNC 'pawrad_print' 453 !End of the abilint section 454 455 implicit none 456 457 !Arguments ------------------------------------ 458 integer,intent(in),optional :: prtvol,unit 459 character(len=4),intent(in),optional :: mode_paral 460 character(len=*),intent(in),optional :: header 461 type(pawrad_type),intent(in) :: Rmesh 462 463 !Local variables------------------------------- 464 !scalars 465 integer :: my_unt,my_prtvol 466 character(len=4) :: my_mode 467 character(len=500) :: msg 468 469 ! ************************************************************************* 470 471 !@pawrad_type 472 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 473 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 474 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 475 476 msg=ch10//' ==== Info on the Radial Mesh ==== ' 477 if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== ' 478 call wrtout(my_unt,msg,my_mode) 479 480 SELECT CASE (Rmesh%mesh_type) 481 482 CASE (RMESH_LINEAR) 483 write(msg,'(a,i4,a,g12.5)')& 484 & ' - Linear mesh: r(i)=step*(i-1), size=',Rmesh%mesh_size,', step=',Rmesh%rstep 485 486 CASE (RMESH_LOG1) 487 write(msg,'(a,i4,2(a,g12.5))')& 488 & ' - Logarithimc mesh: r(i)=AA*[exp(BB*(i-1))-1], size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep 489 490 CASE (RMESH_LOG2) 491 write(msg,'(a,i4,2(a,g12.5))')& 492 & ' - Logarithimc mesh: r(i)=AA*exp(BB*(i-2)), size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep 493 494 CASE (RMESH_LOG3) 495 write(msg,'(a,i1,a,i4,a,g12.5)')& 496 & ' - Logarithimc mesh: r(i)=-AA*ln(1-(i-1)/n), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep 497 498 CASE (RMESH_NL) 499 write(msg,'(a,i1,a,i4,a,g12.5)')& 500 & ' - Non-linear mesh: r(i)=-AA*i/(n-i), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep 501 502 CASE DEFAULT 503 msg = ' Unknown mesh type! Action : check your pseudopotential or input file.' 504 MSG_ERROR(msg) 505 END SELECT 506 507 call wrtout(my_unt,msg,my_mode) 508 509 if (my_prtvol>1) then 510 write(msg,'(a,i4)')' Mesh size for integrals = ',Rmesh%int_meshsz 511 call wrtout(my_unt,msg,my_mode) 512 write(msg,'(a,g12.5)')' rmax=rad(mesh_size) = ',Rmesh%rmax 513 call wrtout(my_unt,msg,my_mode) 514 write(msg,'(a,g12.5)')' Value of stepint = ',Rmesh%stepint 515 call wrtout(my_unt,msg,my_mode) 516 end if 517 518 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
84 type, public :: pawrad_type 85 86 !Integer scalars 87 88 integer :: int_meshsz=0 89 ! Mesh size used in integrals computation 90 ! Integrals will be computed up to r(int_meshsz) 91 92 integer :: mesh_size=0 93 ! Dimension of radial mesh 94 95 integer :: mesh_type=-1 96 ! Type of mesh 97 ! 1=regular grid: r(i)=(i-1)*AA 98 ! 2=logarithmic grid: r(i)=AA*(exp[BB*(i-1)]-1) 99 ! 3=logarithmic grid: r(i>1)=AA*exp[BB*(i-1)] and r(1)=0 100 ! 4=logarithmic grid: r(i)=-AA*ln[1-BB*(i-1)] with BB=1/n 101 102 !Real (real(dp)) scalars 103 104 real(dp) :: lstep=zero 105 ! Exponential step of the mesh (BB parameter above) 106 ! Defined only if mesh type is logarithmic 107 108 real(dp) :: rmax=zero 109 ! Max. value of r = rad(mesh_size) 110 111 real(dp) :: rstep=zero 112 ! Radial step of the mesh (AA parameter above) 113 114 real(dp) :: stepint=zero 115 ! Radial step used to convert any function from the 116 ! present grid onto a regular grid in order to 117 ! integrate it using trapeze method 118 119 !Real (real(dp)) arrays 120 121 real(dp), allocatable :: rad(:) 122 ! rad(mesh_size) 123 ! Coordinates of all the points of the mesh 124 125 real(dp), allocatable :: radfact(:) 126 ! radfact(mesh_size) 127 ! Factor used to compute radial integrals 128 ! Before being integrated on the present mesh, 129 ! any function is multiplied by this factor 130 131 real(dp), allocatable :: simfact(:) 132 ! simfact(mesh_size) 133 ! Factor used to compute radial integrals by the a Simpson scheme 134 ! Integral[f] = Sum_i [simfact(i)*f(i)] 135 136 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''])
PARENTS
calc_ubare,m_paw_atom,m_pawpsp,m_pawrad,pawdenpot,pawinit,pawpuxinit
CHILDREN
poisson,simp_gen
SOURCE
1357 subroutine poisson(den,ll,radmesh,rv,screened_sr_separation,qq) 1358 1359 1360 !This section has been created automatically by the script Abilint (TD). 1361 !Do not modify the following lines by hand. 1362 #undef ABI_FUNC 1363 #define ABI_FUNC 'poisson' 1364 !End of the abilint section 1365 1366 implicit none 1367 1368 !Arguments --------------------------------------------- 1369 !scalars 1370 integer,intent(in) :: ll 1371 real(dp),intent(out),optional :: qq 1372 real(dp),intent(in),optional :: screened_sr_separation 1373 type(pawrad_type),intent(in) :: radmesh 1374 !arrays 1375 real(dp),intent(in) :: den(:) 1376 real(dp),intent(out) :: rv(:) 1377 1378 !Local variables --------------------------------------- 1379 !scalars 1380 integer :: ir,jr,mesh_size,mm,nn 1381 logical :: use_numerov,use_screened 1382 real(dp) :: angm,hh,intg,qq_,ri,rj,rr,sr_omega 1383 !arrays 1384 real(dp) :: ee(4) 1385 real(dp),allocatable :: aa(:),bb(:),cc(:),dd(:),radl(:),radl1(:) 1386 1387 ! *************************************************************************** 1388 1389 mesh_size=size(den) 1390 if (size(rv)/=mesh_size.or.mesh_size>radmesh%mesh_size) then 1391 MSG_BUG('wrong sizes!') 1392 end if 1393 1394 use_numerov=(radmesh%mesh_type==1) 1395 1396 use_screened=.false.;sr_omega=zero 1397 if (present(screened_sr_separation)) then 1398 sr_omega=screened_sr_separation 1399 use_screened=(sr_omega>tol8) 1400 ! Numerov method not coded for screened Coulomb interaction 1401 if (use_numerov.and.use_screened) use_numerov=.false. 1402 end if 1403 1404 !============================================== 1405 !UNIFORM GRID - NUMEROV ALGORITHM 1406 !============================================== 1407 if (use_numerov) then 1408 hh=radmesh%rstep 1409 nn=radmesh%int_meshsz-1 1410 LIBPAW_ALLOCATE(aa,(nn)) 1411 LIBPAW_ALLOCATE(bb,(nn)) 1412 LIBPAW_ALLOCATE(cc,(nn+1)) 1413 do ir=1,nn 1414 aa(ir)=two*hh*den(ir+1)/(ir) 1415 bb(ir)=den(ir+1)*((ir*hh)**ll) 1416 end do 1417 cc(1)=zero 1418 cc(2:nn+1)=bb(1:nn) 1419 call simp_gen(qq_,cc,radmesh) 1420 qq_=qq_/dble(2*ll+1) 1421 rv(1)=aa(1)+0.1_dp*aa(2) 1422 do ir=2,nn-1 1423 rv(ir)=aa(ir)+0.1_dp*(aa(ir+1)+aa(ir-1)) 1424 end do 1425 rv(nn)=aa(nn)+0.1_dp*aa(nn-1) 1426 angm=dble(ll*(ll+1)) 1427 rr=(nn+1)*hh 1428 rv(nn)=rv(nn)+(2.4_dp-0.2_dp*angm/((nn+1)**2))*qq_/(rr**ll) 1429 do ir=1,nn 1430 aa(ir)=angm/(ir*ir) 1431 bb(ir)=2.4_dp+aa(ir) 1432 end do 1433 do ir=1,nn-1 1434 cc(ir)=-1.2_dp+0.1_dp*aa(ir+1) 1435 end do 1436 do ir=nn,2,-1 1437 aa(ir)=-1.2_dp+0.1_dp*aa(ir-1) 1438 end do 1439 if (nn.eq.1) then 1440 rv(2)=rv(1)/bb(1) 1441 rv(1)=zero 1442 else 1443 do ir=2,nn 1444 bb(ir)=bb(ir)-aa(ir)*cc(ir-1)/bb(ir-1) 1445 end do 1446 rv(1)=rv(1)/bb(1) 1447 do ir=2,nn 1448 rv(ir)=(rv(ir)-aa(ir)*rv(ir-1))/bb(ir) 1449 end do 1450 do ir=nn-1,1,-1 1451 rv(ir)=rv(ir)-cc(ir)*rv(ir+1)/bb(ir) 1452 end do 1453 do ir=nn+1,2,-1 1454 rv(ir)=rv(ir-1) 1455 end do 1456 rv(1)=zero 1457 if (nn+1<mesh_size) rv(nn+2:mesh_size)=zero 1458 end if 1459 rv(:)=half*rv(:) 1460 if (present(qq)) qq=qq_ 1461 LIBPAW_DEALLOCATE(aa) 1462 LIBPAW_DEALLOCATE(bb) 1463 LIBPAW_DEALLOCATE(cc) 1464 1465 ! ============================================== 1466 ! ANY OTHER GRID - SIMPSON ALGORITHM 1467 ! ============================================== 1468 else 1469 nn=mesh_size;hh=third*radmesh%stepint 1470 do while (abs(den(nn))<tol16.and.nn>radmesh%int_meshsz) 1471 nn=nn-1 1472 end do 1473 mm=nn;if (radmesh%mesh_type==3) mm=mm-1 1474 LIBPAW_ALLOCATE(aa,(nn)) 1475 LIBPAW_ALLOCATE(bb,(nn)) 1476 LIBPAW_ALLOCATE(cc,(nn)) 1477 LIBPAW_ALLOCATE(dd,(nn)) 1478 1479 if (.not.use_screened) then 1480 1481 ! Standard Coulomb integral 1482 LIBPAW_ALLOCATE(radl,(nn)) 1483 LIBPAW_ALLOCATE(radl1,(nn)) 1484 do jr=nn,2,-1 1485 ir=nn-jr+1 1486 radl(jr) =radmesh%rad(jr)**ll 1487 radl1(jr)=radmesh%rad(jr)*radl(jr) 1488 aa(ir)=den(jr)*radmesh%radfact(jr)*radl(jr) 1489 bb(ir)=den(jr)*radmesh%radfact(jr)/radl1(jr) 1490 end do 1491 radl(1)=zero;radl1(1)=zero 1492 ee(2)=aa(nn-1);ee(3)=aa(nn-2);ee(4)=aa(nn-3) 1493 call pawrad_deducer0(ee,4,radmesh) 1494 aa(nn)=ee(1) 1495 ee(2)=bb(nn-1);ee(3)=bb(nn-2);ee(4)=bb(nn-3) 1496 call pawrad_deducer0(ee,4,radmesh) 1497 bb(nn)=ee(1) 1498 cc(1)=zero;dd(1)=zero 1499 do ir=3,mm,2 1500 cc(ir) =cc(ir-2)+hh*(aa(ir-2)+four*aa(ir-1)+aa(ir)) 1501 cc(ir-1)=cc(ir-2)+hh*(1.25_dp*aa(ir-2)+two*aa(ir-1)-quarter*aa(ir)) 1502 dd(ir) =dd(ir-2)+hh*(bb(ir-2)+four*bb(ir-1)+bb(ir)) 1503 dd(ir-1)=dd(ir-2)+hh*(1.25_dp*bb(ir-2)+two*bb(ir-1)-quarter*bb(ir)) 1504 end do 1505 if (mod(mm,2)==0) then 1506 ! cc(mm)=cc(mm-2)+hh*(aa(mm-2)+four*aa(mm-1)+aa(mm)) 1507 ! dd(mm)=dd(mm-2)+hh*(bb(mm-2)+four*bb(mm-1)+bb(mm)) 1508 cc(mm)=cc(mm-1)+hh*(1.25_dp*aa(mm-2)+two*aa(mm-1)-quarter*aa(mm)) 1509 dd(mm)=dd(mm-1)+hh*(1.25_dp*bb(mm-2)+two*bb(mm-1)-quarter*bb(mm)) 1510 end if 1511 if (mm<nn) then 1512 cc(nn)=cc(mm)+half*(aa(mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1)) 1513 dd(nn)=dd(mm)+half*(bb(mm)+bb(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1)) 1514 end if 1515 rv(1)=zero 1516 do ir=2,nn 1517 jr=nn-ir+1 1518 rv(ir)=(dd(jr)*radl1(ir)+(cc(nn)-cc(jr))/radl(ir))/(two*ll+one) 1519 end do 1520 if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn) 1521 if (present(qq)) qq=cc(nn)/(two*ll+one) 1522 LIBPAW_DEALLOCATE(radl) 1523 LIBPAW_DEALLOCATE(radl1) 1524 else 1525 1526 ! Short-range screened Coulomb integral 1527 rv(1)=zero 1528 do ir=2,nn 1529 ri=sr_omega*radmesh%rad(ir) 1530 do jr=2,nn 1531 rj=sr_omega*radmesh%rad(jr) 1532 aa(jr)=den(jr)*radmesh%radfact(jr)*sr_omega*screened_coul_kernel(ll,ri,rj) 1533 end do 1534 call pawrad_deducer0(aa(1:4),4,radmesh) 1535 intg=zero 1536 ! Compute the integral in 2 parts because the function is not differentiable at r(ir) 1537 ! Integral from zero to r(ir) 1538 if (ir>2) then 1539 do jr=ir-2,1,-2 1540 intg=intg+hh*(aa(jr)+four*aa(jr+1)+aa(jr+2)) 1541 end do 1542 if (mod(ir,2)==0) intg=intg+hh*(-quarter*aa(1)+two*aa(2)+1.25_dp*aa(3)) 1543 else if (ir==2) then 1544 if (nn >2) intg=intg+hh*(+1.25_dp*aa(1)+two*aa(2)-quarter*aa(3)) 1545 if (nn==2) intg=intg+(aa(1)+aa(2))*hh*half 1546 end if 1547 if (mm<nn) intg=intg+half*(aa(1+nn-mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1)) 1548 ! Integral from r(ir) to rmax 1549 if (ir<nn-2) then 1550 do jr=ir+2,nn,2 1551 intg=intg+hh*(aa(jr-2)+four*aa(jr-1)+aa(jr)) 1552 end do 1553 if (mod(nn-ir,2)==1) then 1554 if (ir<=nn-4) then 1555 !Use a high-order formula because points are spaced 1556 intg=intg+hh*(251._dp*aa(nn)+646._dp*aa(nn-1)-264._dp*aa(nn-2) & 1557 & +106._dp*aa(nn-3)-19._dp*aa(nn-4))/240._dp 1558 else 1559 intg=intg+hh*(1.25_dp*aa(nn)+two*aa(nn-1)-quarter*aa(nn-2)) 1560 endif 1561 end if 1562 else if (ir==nn-1) then 1563 intg=intg+(aa(nn-1)+aa(nn))*hh*half 1564 end if 1565 rv(ir)=intg/(two*ll+one)*radmesh%rad(ir) 1566 end do 1567 if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn) 1568 if (present(qq)) qq=zero ! Not relevant 1569 1570 end if 1571 LIBPAW_DEALLOCATE(aa) 1572 LIBPAW_DEALLOCATE(bb) 1573 LIBPAW_DEALLOCATE(cc) 1574 LIBPAW_DEALLOCATE(dd) 1575 1576 end if 1577 1578 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
PARENTS
CHILDREN
SOURCE
1682 function screened_coul_kernel(order,r1,r2,formula) 1683 1684 1685 !This section has been created automatically by the script Abilint (TD). 1686 !Do not modify the following lines by hand. 1687 #undef ABI_FUNC 1688 #define ABI_FUNC 'screened_coul_kernel' 1689 !End of the abilint section 1690 1691 implicit none 1692 1693 !Arguments ------------------------------------ 1694 !scalars 1695 integer,intent(in) :: order 1696 integer,optional :: formula 1697 real(dp),intent(in) :: r1,r2 1698 real(dp) :: screened_coul_kernel 1699 !Local variables------------------------------- 1700 integer :: formula_ 1701 real(dp) :: difexp,erfcp,erfcm,f0,f1,f2,f3,f4,f5,f6,hh,sqrtpi,sumexp,xx,xy,yy 1702 character(len=100) :: msg 1703 1704 !****************************************************************** 1705 1706 if (order>6) then 1707 msg='PAW screened exchange not coded for l>2!' 1708 MSG_ERROR(msg) 1709 end if 1710 1711 !Use max and min of arguments 1712 xx=max(r1,r2) ; yy=min(r1,r2) 1713 1714 !Unscreened Coulomb interaction for very small (or negative) arguments 1715 if (xx<tol8) then 1716 screened_coul_kernel=yy**order/xx**(order+1) ; return 1717 end if 1718 1719 !Choice of formula 1720 !Empirical criterion, inspired by J. Phys. A 39, pp8624 [[cite:Angyan2006]] and adjusted 1721 formula_=1;if (xx<0.25_dp.and.yy<0.25_dp) formula_=2 1722 if (present(formula)) formula_=formula 1723 select case (formula_) 1724 1725 !Full formula 1726 !J. Phys. A 39, 8613 (2006) - Eq. (21), (22), (24) [[cite:Angyan2006]] 1727 ! Note: typo in the paper: in eq (22), the sum begins at p=0 1728 !------------------------------------------------------------------ 1729 case(1) 1730 xy=xx*yy ; sqrtpi=sqrt(pi) 1731 sumexp=exp(-(xx+yy)**2)+exp(-(xx-yy)**2) 1732 difexp=exp(-(xx+yy)**2)-exp(-(xx-yy)**2) 1733 erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy) 1734 if (order==0) then 1735 f0=-difexp/(two*sqrtpi*xy) 1736 hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy) 1737 screened_coul_kernel = hh + f0 1738 else if(order==1) then 1739 f0=-difexp/(two*sqrtpi*xy) 1740 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1741 hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2) 1742 screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy 1743 else if (order==2) then 1744 f0=-difexp/(two*sqrtpi*xy) 1745 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1746 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1747 hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3) 1748 screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2 1749 else if (order==3) then 1750 f0=-difexp/(two*sqrtpi*xy) 1751 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1752 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1753 f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4) 1754 hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4) 1755 screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 & 1756 & + f0*(xx**6+yy**6)/xy**3 1757 else if (order==4) then 1758 f0=-difexp/(two*sqrtpi*xy) 1759 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1760 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1761 f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4) 1762 f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp & 1763 & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5) 1764 hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5) 1765 screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 & 1766 & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4 1767 else if (order==5) then 1768 f0=-difexp/(two*sqrtpi*xy) 1769 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1770 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1771 f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4) 1772 f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp & 1773 & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5) 1774 f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp & 1775 & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6) 1776 hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6) 1777 screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 & 1778 & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 & 1779 & + f0*(xx**10+yy**10)/xy**5 1780 else if (order==6) then 1781 f0=-difexp/(two*sqrtpi*xy) 1782 f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2) 1783 f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3) 1784 f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4) 1785 f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp & 1786 & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5) 1787 f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp & 1788 & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6) 1789 f6=-((64._dp*xy**6+3360._dp*xy**4+18900._dp*xy**2+10395._dp)*difexp+ & 1790 & (672._dp*xy**5+10080._dp*xy**3+20790._dp*xy)*sumexp)/(sqrtpi*(two*xy)**7) 1791 hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7) 1792 screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 & 1793 & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 & 1794 & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6 1795 end if 1796 1797 !Development for yy->0 1798 !J. Phys. A 39, 8613 (2006) - Eq. (28), (29), (30) [[cite:Angyan2006]] 1799 !------------------------------------------------------------------ 1800 case(2) 1801 sqrtpi=sqrt(pi) 1802 if (order==0) then 1803 screened_coul_kernel = paw_derfc(xx)/xx + exp(-xx**2)/sqrtpi * & 1804 & (two/three *yy**2 & 1805 & +(two*xx**2-three)/15._dp*yy**4) 1806 else if(order==1) then 1807 screened_coul_kernel = paw_derfc(xx)*yy/xx**2 + exp(-xx**2)/sqrtpi * & 1808 & (two/xx *yy & 1809 & +four*xx/5._dp *yy**3 & 1810 & +two*xx*(two*xx**2-5._dp)/35._dp*yy**5) 1811 else if (order==2) then 1812 screened_coul_kernel = paw_derfc(xx)*yy**2/xx**3 + exp(-xx**2)/sqrtpi * & 1813 & (two*(two*xx**2+three)/(three*xx**2) *yy**2 & 1814 & +8._dp*xx**2/21._dp *yy**4 & 1815 & +four*xx**2*(two*xx**2-7._dp)/189._dp*yy**6) 1816 else if (order==3) then 1817 screened_coul_kernel = paw_derfc(xx)*yy**3/xx**4 + exp(-xx**2)/sqrtpi * & 1818 & (two*(four*xx**4+10._dp*xx**2+15._dp)/(15._dp*xx**3)*yy**3 & 1819 & +16._dp*xx**3/135._dp *yy**5 & 1820 & +8._dp*xx**3*(two*xx**2-9._dp)/1485._dp *yy**7) 1821 else if (order==4) then 1822 screened_coul_kernel = paw_derfc(xx)*yy**4/xx**5 + exp(-xx**2)/sqrtpi * & 1823 & (two*(8._dp*xx**6+28._dp*xx**4+70._dp*xx**2+105._dp)/(105._dp*xx**4)*yy**4 & 1824 & +32._dp*xx**4/1155._dp *yy**6 & 1825 & +16._dp*xx**4*(two*xx**2-11._dp)/15015._dp *yy**8) 1826 else if (order==5) then 1827 screened_coul_kernel = paw_derfc(xx)*yy**5/xx**6 + exp(-xx**2)/sqrtpi * & 1828 & (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 & 1829 & +64._dp*xx**5/12285._dp *yy**7 & 1830 & +32._dp*xx**5*(two*xx**2-13._dp)/184275._dp *yy**9) 1831 else if (order==6) then 1832 screened_coul_kernel = paw_derfc(xx)*yy**6/xx**7 + exp(-xx**2)/sqrtpi * & 1833 & (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 & 1834 & +128._dp*xx**6/155925._dp *yy**8 & 1835 & +64._dp*xx**6*(two*xx**2-15._dp)/2650725._dp *yy**10) 1836 end if 1837 1838 !Development for xx*yy->0 1839 !J. Phys. A 39, 8613 (2006) - Eq. (24), (26) [[cite:Angyan2006]] 1840 ! Note: typo in the paper: (2n+3)! should be (2n+3)!! 1841 !------------------------------------------------------------------ 1842 case(3) 1843 xy=xx*yy ; sqrtpi=sqrt(pi) 1844 sumexp=exp(-(xx**2+yy**2)) 1845 erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy) 1846 if (order==0) then 1847 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1848 hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy) 1849 screened_coul_kernel = hh + f0 1850 else if(order==1) then 1851 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1852 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1853 hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2) 1854 screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy 1855 else if (order==2) then 1856 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1857 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1858 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1859 hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3) 1860 screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2 1861 else if (order==3) then 1862 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1863 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1864 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1865 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1866 hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4) 1867 screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 & 1868 & + f0*(xx**6+yy**6)/xy**3 1869 else if (order==4) then 1870 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1871 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1872 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1873 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1874 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1875 hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5) 1876 screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 & 1877 & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4 1878 else if (order==5) then 1879 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1880 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1881 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1882 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1883 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1884 f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp 1885 hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6) 1886 screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 & 1887 & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 & 1888 & + f0*(xx**10+yy**10)/xy**5 1889 else if (order==6) then 1890 f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp 1891 f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp 1892 f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp 1893 f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp 1894 f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp 1895 f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp 1896 f6 = 128._dp*xy**2*(15._dp+two*xy**2)/(1027025._dp*sqrtpi)*sumexp 1897 hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7) 1898 screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 & 1899 & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 & 1900 & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6 1901 end if 1902 1903 end select 1904 1905 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
PARENTS
calc_ubare,m_atom,m_paw_atom,m_paw_slater,m_pawdij,m_pawhr,m_pawpsp m_pawpwij,m_pawrad,m_pawxc,m_plowannier,m_psps,make_efg_onsite mlwfovlp_projpaw,optics_paw,optics_paw_core,partial_dos_fractions_paw pawdensities,pawinit,pawnabla_init,pawpuxinit,posdoppler,poslifetime posratecore,qijb_kk,smatrix_pawinit
CHILDREN
poisson,simp_gen
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
948 subroutine simp_gen(intg,func,radmesh,r_for_intg) 949 950 #if defined HAVE_AVX_SAFE_MODE 951 !DEC$ NOOPTIMIZE 952 #endif 953 954 !This section has been created automatically by the script Abilint (TD). 955 !Do not modify the following lines by hand. 956 #undef ABI_FUNC 957 #define ABI_FUNC 'simp_gen' 958 !End of the abilint section 959 960 implicit none 961 962 !Arguments ------------------------------------ 963 !scalars 964 real(dp),intent(out) :: intg 965 real(dp),intent(in),optional :: r_for_intg 966 type(pawrad_type),intent(in) :: radmesh 967 !arrays 968 real(dp),intent(in) :: func(:) 969 970 !Local variables------------------------------- 971 !scalars 972 integer :: ii,int_meshsz,ir,ir_last,isim,nn 973 real(dp) :: hh,resid,simp 974 real(dp),allocatable :: simfact(:) 975 character(len=500) :: msg 976 977 ! ************************************************************************* 978 979 if (present(r_for_intg)) then 980 if (r_for_intg>0.d0) then 981 ir=min(pawrad_ifromr(radmesh,r_for_intg),radmesh%mesh_size) 982 if (ir<radmesh%mesh_size) then 983 if (abs(radmesh%rad(ir+1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir+1 984 end if 985 if (ir>1) then 986 if (abs(radmesh%rad(ir-1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir-1 987 end if 988 int_meshsz=ir 989 else 990 int_meshsz=radmesh%mesh_size 991 end if 992 if (int_meshsz>radmesh%mesh_size.or.int_meshsz>size(func)) then 993 write(msg,'(3(a,i4))')"int_meshsz= ",int_meshsz," > mesh_size=",radmesh%mesh_size,& 994 & ", size(func)=",size(func) 995 MSG_BUG(msg) 996 end if 997 isim=3; if (radmesh%mesh_type==3)isim=4 998 LIBPAW_ALLOCATE(simfact,(radmesh%mesh_size)) 999 hh=radmesh%stepint/3.d0 1000 simfact(int_meshsz)=hh*radmesh%radfact(int_meshsz) 1001 simfact(1:isim-2)=zero 1002 ir_last=1 1003 do ir=int_meshsz,isim,-2 1004 simfact(ir-1)=4.d0*hh*radmesh%radfact(ir-1) 1005 simfact(ir-2)=2.d0*hh*radmesh%radfact(ir-2) 1006 ir_last=ir-2 1007 end do 1008 simfact(ir_last)=half*simfact(ir_last) 1009 if (int_meshsz<radmesh%mesh_size) simfact(int_meshsz+1:radmesh%mesh_size)=zero 1010 1011 nn=int_meshsz 1012 simp=zero 1013 do ii=1,nn 1014 simp=simp+func(ii)*simfact(ii) 1015 end do 1016 LIBPAW_DEALLOCATE(simfact) 1017 1018 else 1019 if (radmesh%int_meshsz>size(func)) then 1020 write(msg,'(2(a,i4))')"int_meshsz= ",int_meshsz," > size(func)=",size(func) 1021 MSG_BUG(msg) 1022 end if 1023 nn=radmesh%int_meshsz 1024 simp=zero 1025 do ii=1,nn 1026 simp=simp+func(ii)*radmesh%simfact(ii) 1027 end do 1028 end if 1029 1030 resid=zero 1031 if (radmesh%mesh_type==3) then 1032 resid=half*(func(2)+func(1))*(radmesh%rad(2)-radmesh%rad(1)) 1033 if (mod(nn,2)==1) resid=resid+radmesh%stepint/3.d0*(1.25d0*func(2)*radmesh%radfact(2) & 1034 & +2.d0*func(3)*radmesh%radfact(3)-0.25d0*func(4)*radmesh%radfact(4)) 1035 else if (mod(nn,2)==0) then 1036 resid=radmesh%stepint/3.d0*(1.25d0*func(1)*radmesh%radfact(1)+2.d0*func(2)*radmesh%radfact(2) & 1037 & -0.25d0*func(3)*radmesh%radfact(3)) 1038 end if 1039 1040 intg=simp+resid 1041 1042 end subroutine simp_gen