1!! gen1int: compute one-electron integrals using rotational London atomic-orbitals 2!! Copyright 2009-2012 Bin Gao, and Andreas Thorvaldsen 3!! 4!! gen1int is free software: you can redistribute it and/or modify 5!! it under the terms of the GNU Lesser General Public License as published by 6!! the Free Software Foundation, either version 3 of the License, or 7!! (at your option) any later version. 8!! 9!! gen1int is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU Lesser General Public License for more details. 13!! 14!! You should have received a copy of the GNU Lesser General Public License 15!! along with gen1int. If not, see <http://www.gnu.org/licenses/>. 16!! 17!! Recursive functions of Dirac delta function integrals. 18!! 19!! 2012-03-18, Bin Gao 20!! * first version 21 22#include "stdout.h" 23 24!> \brief recursive functions of Dirac delta function integrals 25!> \author Bin Gao 26!> \date 2012-03-18 27module recur_delta 28 use xkind 29 implicit none 30 real(REALK), private, save :: zero_pint !Gauss-Weierstrass transform of Dirac delta function 31 real(REALK), private, save :: half_nr_a !half of the negative reciprocal of exponent on bra center 32 real(REALK), private, save :: half_nr_b !half of the negative reciprocal of exponent on ket center 33 real(REALK), private, save :: neg_twice_p !negative of the twice total exponent 34 real(REALK), private, save :: delta_wrt_bra(3) !relative coordinates of Dirac delta function origin ... 35 !w.r.t. bra center 36 real(REALK), private, save :: delta_wrt_ket(3) !relative coordinates of Dirac delta function origin ... 37 !w.r.t. ket center 38 real(REALK), private, save :: delta_wrt_cc(3) !relative coordinates of Dirac delta function origin ... 39 !w.r.t. center-of-charge 40 real(REALK), private, save :: delta_wrt_diporg(3) !relative coordinates of Dirac delta function origin ... 41 !w.r.t. dipole origin 42 43 public :: recur_delta_create 44 public :: recur_delta_geom 45 public :: recur_delta_hbra 46 public :: recur_delta_hket 47 public :: recur_delta_moment 48 49 contains 50 51 !> \brief prepares the prerequisite quantities used in the following recursive functions 52 !> \author Bin Gao 53 !> \date 2011-07-27 54 !> \param coord_bra is the coordinates of bra center 55 !> \param exponent_bra is the exponent of primitive Gaussian of bra center 56 !> \param coord_ket is the coordinates of ket center 57 !> \param exponent_ket is the exponent of primitive Gaussian of ket center 58 !> \param delta_origin contains the coordinates of Dirac delta function origin 59 !> \param dipole_origin contains the coordinates of dipole origin 60 !> \param scal_const is the scale constant 61 !> \param order_elec is the order of electronic derivatives 62 subroutine recur_delta_create(coord_bra, exponent_bra, coord_ket, exponent_ket, & 63 delta_origin, dipole_origin, scal_const, order_elec) 64 real(REALK), intent(in) :: coord_bra(3) 65 real(REALK), intent(in) :: exponent_bra 66 real(REALK), intent(in) :: coord_ket(3) 67 real(REALK), intent(in) :: exponent_ket 68 real(REALK), intent(in) :: delta_origin(3) 69 real(REALK), intent(in) :: dipole_origin(3) 70 real(REALK), intent(in) :: scal_const 71 integer, intent(in) :: order_elec 72 real(REALK) neg_total_expnt !negative total exponent 73 real(REALK) recip_total_expnt !reciprocal of total exponent 74 real(REALK) neg_reduced_expnt !negative reduced exponent 75 real(REALK) sd_bra_ket !square of the relative distance between bra and ket centers 76 real(REALK) sd_delta_cc !square of the relative distance between delta function ... 77 !and center-of-charge 78 integer ixyz !incremental recorder over xyz components 79 ! sets the half of the negative reciprocal of exponents 80 half_nr_a = -0.5_REALK/exponent_bra 81 half_nr_b = -0.5_REALK/exponent_ket 82 ! computes the negative total exponent, reciprocal of total exponent and negative reduced exponent 83 neg_total_expnt = -(exponent_bra+exponent_ket) 84 recip_total_expnt = 1.0_REALK/neg_total_expnt 85 neg_reduced_expnt = exponent_bra*exponent_ket*recip_total_expnt 86 recip_total_expnt = -recip_total_expnt 87 ! sets the negative of twice the total exponents 88 neg_twice_p = 2.0_REALK*neg_total_expnt 89 ! computes the relative coordinates, and squares of the relative distance 90 ! between bra and ket centers, delta function and center-of-charge 91 sd_bra_ket = 0.0_REALK 92 sd_delta_cc = 0.0_REALK 93 do ixyz = 1, 3 94 delta_wrt_bra(ixyz) = delta_origin(ixyz)-coord_bra(ixyz) 95 delta_wrt_ket(ixyz) = delta_origin(ixyz)-coord_ket(ixyz) 96 delta_wrt_cc(ixyz) = delta_origin(ixyz)-(exponent_bra*coord_bra(ixyz) & 97 + exponent_ket*coord_ket(ixyz))*recip_total_expnt 98 delta_wrt_diporg(ixyz) = delta_origin(ixyz)-dipole_origin(ixyz) 99 sd_bra_ket = sd_bra_ket+(coord_ket(ixyz)-coord_bra(ixyz))**2 100 sd_delta_cc = sd_delta_cc+delta_wrt_cc(ixyz)**2 101 end do 102 ! computes the Gauss-Weierstrass transform of Dirac delta function 103 zero_pint = scal_const*(-exponent_ket-exponent_ket)**order_elec & 104 * exp(neg_reduced_expnt*sd_bra_ket+neg_total_expnt*sd_delta_cc) 105 end subroutine recur_delta_create 106 107 !> \brief recursive function for \fn(delta_geom) 108 !> \author Bin Gao 109 !> \date 2012-02-12 110 !> \param order_geo contains the orders of xyz components of geometric derivatives 111 !> \return geo_pint is the geometric derivative on Dirac delta function 112 recursive function recur_delta_geom(order_geo) result(geo_pint) 113 real(REALK) geo_pint 114 integer, intent(in) :: order_geo(3) 115 ! recurrence relation along x-direction 116 if (order_geo(1)>0) then 117 geo_pint = neg_twice_p*(delta_wrt_cc(1) & 118 * recur_delta_geom((/order_geo(1)-1,order_geo(2:3)/)) & 119 + real(order_geo(1)-1,REALK) & 120 * recur_delta_geom((/order_geo(1)-2,order_geo(2:3)/))) 121 ! recurrence relation along y-direction 122 else if (order_geo(2)>0) then 123 geo_pint = neg_twice_p*(delta_wrt_cc(2) & 124 * recur_delta_geom((/order_geo(1),order_geo(2)-1,order_geo(3)/)) & 125 + real(order_geo(2)-1,REALK) & 126 * recur_delta_geom((/order_geo(1),order_geo(2)-2,order_geo(3)/))) 127 ! recurrence relation along z-direction 128 else if (order_geo(3)>0) then 129 geo_pint = neg_twice_p*(delta_wrt_cc(3) & 130 * recur_delta_geom((/order_geo(1:2),order_geo(3)-1/)) & 131 + real(order_geo(3)-1,REALK) & 132 * recur_delta_geom((/order_geo(1:2),order_geo(3)-2/))) 133 ! zero integral 134 else if (any(order_geo<0)) then 135 geo_pint = 0.0_REALK 136 ! zeroth order geometric derivative 137 else 138 geo_pint = zero_pint 139 end if 140 end function recur_delta_geom 141 142 !> \brief recursive function for \fn(delta_hket) but on bra center 143 !> \author Bin Gao 144 !> \date 2012-03-16 145 !> \param order_hbra contains the orders of xyz components of HGTO on bra center 146 !> \param order_geo contains the orders of xyz components of geometric derivatives 147 !> \return hbra_pint is the primitive integral with requied order of HGTOs on bra center 148 recursive function recur_delta_hbra(order_hbra, order_geo) result(hbra_pint) 149 real(REALK) hbra_pint 150 integer, intent(in) :: order_hbra(3) 151 integer, intent(in) :: order_geo(3) 152 ! recurrence relation along x-direction 153 if (order_hbra(1)>0) then 154 hbra_pint = delta_wrt_bra(1) & 155 * recur_delta_hbra((/order_hbra(1)-1,order_hbra(2:3)/), order_geo) & 156 + real(order_geo(1),REALK) & 157 * recur_delta_hbra((/order_hbra(1)-1,order_hbra(2:3)/), & 158 (/order_geo(1)-1,order_geo(2:3)/)) & 159 + real(order_hbra(1)-1,REALK)*half_nr_a & 160 * recur_delta_hbra((/order_hbra(1)-2,order_hbra(2:3)/), order_geo) 161 ! recurrence relation along y-direction 162 else if (order_hbra(2)>0) then 163 hbra_pint = delta_wrt_bra(2) & 164 * recur_delta_hbra((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 165 order_geo) & 166 + real(order_geo(2),REALK) & 167 * recur_delta_hbra((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 168 (/order_geo(1),order_geo(2)-1,order_geo(3)/)) & 169 + real(order_hbra(2)-1,REALK)*half_nr_a & 170 * recur_delta_hbra((/order_hbra(1),order_hbra(2)-2,order_hbra(3)/), & 171 order_geo) 172 ! recurrence relation along z-direction 173 else if (order_hbra(3)>0) then 174 hbra_pint = delta_wrt_bra(3) & 175 * recur_delta_hbra((/order_hbra(1:2),order_hbra(3)-1/), order_geo) & 176 + real(order_geo(3),REALK) & 177 * recur_delta_hbra((/order_hbra(1:2),order_hbra(3)-1/), & 178 (/order_geo(1:2),order_geo(3)-1/)) & 179 + real(order_hbra(3)-1,REALK)*half_nr_a & 180 * recur_delta_hbra((/order_hbra(1:2),order_hbra(3)-2/), order_geo) 181 ! zero integral 182 else if (any(order_hbra<0) .or. any(order_geo<0)) then 183 hbra_pint = 0.0_REALK 184 ! recurrence relations on the order of geometric derivatives 185 else 186 hbra_pint = recur_delta_geom(order_geo) 187 end if 188 end function recur_delta_hbra 189 190 !> \brief recursive function for \fn(delta_hket) 191 !> \author Bin Gao 192 !> \date 2012-03-16 193 !> \param order_hbra contains the orders of xyz components of HGTO on bra center 194 !> \param order_hket contains the orders of xyz components of HGTO on ket center 195 !> \param order_geo contains the orders of xyz components of geometric derivatives 196 !> \return hket_pint is the primitive integral with requied order of HGTOs on ket center 197 recursive function recur_delta_hket(order_hbra, order_hket, order_geo) & 198 result(hket_pint) 199 real(REALK) hket_pint 200 integer, intent(in) :: order_hbra(3) 201 integer, intent(in) :: order_hket(3) 202 integer, intent(in) :: order_geo(3) 203 ! recurrence relation along x-direction 204 if (order_hket(1)>0) then 205 hket_pint = delta_wrt_ket(1) & 206 * recur_delta_hket(order_hbra, (/order_hket(1)-1,order_hket(2:3)/), & 207 order_geo) & 208 + real(order_geo(1),REALK) & 209 * recur_delta_hket(order_hbra, (/order_hket(1)-1,order_hket(2:3)/), & 210 (/order_geo(1)-1,order_geo(2:3)/)) & 211 + real(order_hket(1)-1,REALK)*half_nr_b & 212 * recur_delta_hket(order_hbra, (/order_hket(1)-2,order_hket(2:3)/), & 213 order_geo) 214 ! recurrence relation along y-direction 215 else if (order_hket(2)>0) then 216 hket_pint = delta_wrt_ket(2) & 217 * recur_delta_hket(order_hbra, & 218 (/order_hket(1),order_hket(2)-1,order_hket(3)/), & 219 order_geo) & 220 + real(order_geo(2),REALK) & 221 * recur_delta_hket(order_hbra, & 222 (/order_hket(1),order_hket(2)-1,order_hket(3)/), & 223 (/order_geo(1),order_geo(2)-1,order_geo(3)/)) & 224 + real(order_hket(2)-1,REALK)*half_nr_b & 225 * recur_delta_hket(order_hbra, & 226 (/order_hket(1),order_hket(2)-2,order_hket(3)/), & 227 order_geo) 228 ! recurrence relation along z-direction 229 else if (order_hket(3)>0) then 230 hket_pint = delta_wrt_ket(3) & 231 * recur_delta_hket(order_hbra, (/order_hket(1:2),order_hket(3)-1/), & 232 order_geo) & 233 + real(order_geo(3),REALK) & 234 * recur_delta_hket(order_hbra, (/order_hket(1:2),order_hket(3)-1/), & 235 (/order_geo(1:2),order_geo(3)-1/)) & 236 + real(order_hket(3)-1,REALK)*half_nr_b & 237 * recur_delta_hket(order_hbra, (/order_hket(1:2),order_hket(3)-2/), & 238 order_geo) 239 ! zero integral 240 else if (any(order_hket<0) .or. any(order_geo<0)) then 241 hket_pint = 0.0_REALK 242 ! recurrence relations on the order of HGTOs on bra center 243 else 244 hket_pint = recur_delta_hbra(order_hbra, order_geo) 245 end if 246 end function recur_delta_hket 247 248 !> \brief recursive function for \fn(delta_moment) 249 !> \author Bin Gao 250 !> \date 2012-03-16 251 !> \param order_hbra contains the orders of xyz components of HGTO on bra center 252 !> \param order_hket contains the orders of xyz components of HGTO on ket center 253 !> \param order_mom contains the orders of xyz components of Cartesian multipole moment 254 !> \param order_geo contains the orders of xyz components of geometric derivatives 255 !> \return hmom_pint is the primitive Cartesian multipole moment integral with requied 256 !> order of geometric derivatives on Dirac delta function 257 recursive function recur_delta_moment(order_hbra, order_hket, order_mom, order_geo) & 258 result(hmom_pint) 259 real(REALK) hmom_pint 260 integer, intent(in) :: order_hbra(3) 261 integer, intent(in) :: order_hket(3) 262 integer, intent(in) :: order_mom(3) 263 integer, intent(in) :: order_geo(3) 264 ! recurrence relation along x-direction 265 if (order_mom(1)>0) then 266 hmom_pint = delta_wrt_diporg(1) & 267 * recur_delta_moment(order_hbra, order_hket, & 268 (/order_mom(1)-1,order_mom(2:3)/), & 269 order_geo) & 270 + real(order_geo(1),REALK) & 271 * recur_delta_moment(order_hbra, order_hket, & 272 (/order_mom(1)-1,order_mom(2:3)/), & 273 (/order_geo(1)-1,order_geo(2:3)/)) 274 ! recurrence relation along y-direction 275 else if (order_mom(2)>0) then 276 hmom_pint = delta_wrt_diporg(2) & 277 * recur_delta_moment(order_hbra, order_hket, & 278 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 279 order_geo) & 280 + real(order_geo(2),REALK) & 281 * recur_delta_moment(order_hbra, order_hket, & 282 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 283 (/order_geo(1),order_geo(2)-1,order_geo(3)/)) 284 ! recurrence relation along z-direction 285 else if (order_mom(3)>0) then 286 hmom_pint = delta_wrt_diporg(3) & 287 * recur_delta_moment(order_hbra, order_hket, & 288 (/order_mom(1:2),order_mom(3)-1/), & 289 order_geo) & 290 + real(order_geo(3),REALK) & 291 * recur_delta_moment(order_hbra, order_hket, & 292 (/order_mom(1:2),order_mom(3)-1/), & 293 (/order_geo(1:2),order_geo(3)-1/)) 294 ! zero integral 295 else if (any(order_mom<0) .or. any(order_geo<0)) then 296 hmom_pint = 0.0_REALK 297 ! recurrence relations on the order of HGTOs on ket center 298 else 299 hmom_pint = recur_delta_hket(order_hbra, order_hket, order_geo) 300 end if 301 end function recur_delta_moment 302 303end module recur_delta 304