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 nuclear attraction potential integrals. 18!! 19!! 2012-03-19, Bin Gao 20!! * first version 21 22#include "stdout.h" 23 24!> \brief recursive functions of nuclear attraction potential integrals 25!> \author Bin Gao 26!> \date 2012-03-19 27module recur_nucpot 28 use xkind 29 implicit none 30! constant Pi 31#include "private/pi.h" 32 real(REALK), private, save :: half_recip_p !half of the reciprocal of total exponent \f$p_{ij}\f$ 33 real(REALK), private, save :: half_neg_rp !half of the negative reciprocal of total exponent \f$p_{ij}\f$ 34 real(REALK), private, save :: bra_to_ket !ratio of exponent on bra center to that on ket center 35 real(REALK), private, save :: ket_to_bra !ratio of exponent on ket center to that on bra center 36 real(REALK), private, save :: cc_wrt_bra(3) !relative coordinates of center-of-charge w.r.t. bra center 37 real(REALK), private, save :: cc_wrt_ket(3) !relative coordinates of center-of-charge w.r.t. ket center 38 real(REALK), private, save :: nucorg_wrt_cc(3) !relative coordinates of nuclear potential origin ... 39 !w.r.t. center-of-charge 40 real(REALK), private, save :: cc_wrt_diporg(3) !relative coordinates of center-of-charge w.r.t. dipole origin 41 integer, private, save :: max_order_aux !maximum order of auxiliary integrals 42 real(REALK), private, allocatable, save :: zero_pint(:) 43 !Gauss-Weierstrass transform of nuclear attraction potential 44 45 public :: recur_nucpot_create 46 public :: recur_nucpot_destroy 47 public :: recur_nucpot_geom 48 public :: recur_nucpot_hket 49 public :: recur_nucpot_hbra 50 public :: recur_nucpot_moment 51 52 contains 53 54 !> \brief prepares the prerequisite quantities used in the following recursive functions 55 !> \author Bin Gao 56 !> \date 2011-07-27 57 !> \param coord_bra is the coordinates of bra center 58 !> \param exponent_bra is the exponent of primitive Gaussian of bra center 59 !> \param coord_ket is the coordinates of ket center 60 !> \param exponent_ket is the exponent of primitive Gaussian of ket center 61 !> \param nucpot_origin contains the coordinates of nuclear potential origin 62 !> \param dipole_origin contains the coordinates of dipole origin 63 !> \param scal_const is the scale constant for operator 64 !> \param order_elec is the order of electronic derivatives 65 !> \param order_aux is the maximum order of auxiliary integrals 66 !> \param gaupot_expt is the exponent used in the Gaussian broadening function of the charge 67 subroutine recur_nucpot_create(coord_bra, exponent_bra, coord_ket, exponent_ket, & 68 nucpot_origin, dipole_origin, scal_const, & 69 order_elec, order_aux, gaupot_expt) 70 real(REALK), intent(in) :: coord_bra(3) 71 real(REALK), intent(in) :: exponent_bra 72 real(REALK), intent(in) :: coord_ket(3) 73 real(REALK), intent(in) :: exponent_ket 74 real(REALK), intent(in) :: nucpot_origin(3) 75 real(REALK), intent(in) :: dipole_origin(3) 76 real(REALK), intent(in) :: scal_const 77 integer, intent(in) :: order_elec 78 integer, intent(in) :: order_aux 79 real(REALK), optional, intent(in) :: gaupot_expt 80 real(REALK) total_expnt !total exponent 81 real(REALK) recip_total_expnt !reciprocal of total exponent 82 real(REALK) coord_cc !xyz coordinates of center-of-charge 83 real(REALK) sd_bra_ket !square of the relative distance between bra and ket centers 84 real(REALK) reduced_expnt !reduced exponent 85 real(REALK) prefact_zero !prefactor of integrals of zero order GTOs and operators 86 real(REALK) arg_aux_fun !argument of auxiliary functions 87 real(REALK) scale_factor !scale factor 88 integer ixyz !incremental recorder over xyz components 89 integer iorder !incremental recorder over auxiliary functions 90 integer ierr !error information 91 ! computes the total exponent and its reciprocal 92 total_expnt = exponent_bra+exponent_ket 93 recip_total_expnt = 1.0_REALK/total_expnt 94 ! computes the half of the (negative) reciprocal of total exponent \f$p_{ij}\f$ 95 half_recip_p = 0.5_REALK*recip_total_expnt 96 half_neg_rp = -half_recip_p 97 ! computes the ratios of exponents 98 bra_to_ket = exponent_bra/exponent_ket 99 ket_to_bra = exponent_ket/exponent_bra 100 ! computes the relative coordinates and square of the relative distance between bra and ket centers 101 sd_bra_ket = 0.0_REALK 102 do ixyz = 1, 3 103 coord_cc = recip_total_expnt*(exponent_bra*coord_bra(ixyz) & 104 + exponent_ket*coord_ket(ixyz)) 105 cc_wrt_bra(ixyz) = coord_cc-coord_bra(ixyz) 106 cc_wrt_ket(ixyz) = coord_cc-coord_ket(ixyz) 107 nucorg_wrt_cc(ixyz) = nucpot_origin(ixyz)-coord_cc 108 cc_wrt_diporg(ixyz) = coord_cc-dipole_origin(ixyz) 109 sd_bra_ket = sd_bra_ket+(coord_ket(ixyz)-coord_bra(ixyz))**2 110 end do 111 ! computes the reduced exponent 112 reduced_expnt = exponent_bra*exponent_ket*recip_total_expnt 113 ! computes the prefactor of auxiliary integrals 114 prefact_zero = scal_const*(-exponent_ket-exponent_ket)**order_elec & 115 * exp(-reduced_expnt*sd_bra_ket) 116 ! computes the auxiliary integrals for Gaussian charge potential 117 if (present(gaupot_expt)) then 118 ! computes the argument of Boys functions 119 arg_aux_fun = 0.0_REALK 120 do ixyz = 1, 3 121 arg_aux_fun = arg_aux_fun+nucorg_wrt_cc(ixyz)**2 122 end do 123 scale_factor = gaupot_expt/(gaupot_expt+total_expnt) 124 total_expnt = scale_factor*total_expnt 125 arg_aux_fun = total_expnt*arg_aux_fun 126 ! sets the maximum order of Boys functions 127 max_order_aux = order_aux 128 allocate(zero_pint(0:max_order_aux), stat=ierr) 129 if (ierr/=0) & 130 call error_stop("recur_nucpot_create", & 131 "failed to allocate zero_pint", max_order_aux+1) 132 ! computes the Boys functions 133 call aux_boys_vec(0, max_order_aux, arg_aux_fun, zero_pint) 134 arg_aux_fun = prefact_zero*(PI+PI)*recip_total_expnt*sqrt(scale_factor) 135 total_expnt = -2.0_REALK*total_expnt 136 ! loops over the order of Boys functions 137 do iorder = 0, max_order_aux 138 zero_pint(iorder) = arg_aux_fun*zero_pint(iorder) 139 arg_aux_fun = arg_aux_fun*total_expnt 140 end do 141 ! computes the auxiliary integrals for nuclear attraction potential 142 else 143 ! computes the argument of Boys functions 144 arg_aux_fun = 0.0_REALK 145 do ixyz = 1, 3 146 arg_aux_fun = arg_aux_fun+nucorg_wrt_cc(ixyz)**2 147 end do 148 arg_aux_fun = total_expnt*arg_aux_fun 149 ! sets the maximum order of Boys functions 150 max_order_aux = order_aux 151 allocate(zero_pint(0:max_order_aux), stat=ierr) 152 if (ierr/=0) & 153 call error_stop("recur_nucpot_create", & 154 "failed to allocate zero_pint", max_order_aux+1) 155 ! computes the Boys functions 156 call aux_boys_vec(0, max_order_aux, arg_aux_fun, zero_pint) 157 arg_aux_fun = prefact_zero*(PI+PI)*recip_total_expnt 158 total_expnt = -2.0_REALK*total_expnt 159 ! loops over the order of Boys functions 160 do iorder = 0, max_order_aux 161 zero_pint(iorder) = arg_aux_fun*zero_pint(iorder) 162 arg_aux_fun = arg_aux_fun*total_expnt 163 end do 164 end if 165 end subroutine recur_nucpot_create 166 167 !> \brief cleans quantities used in the recursive functions 168 !> \author Bin Gao 169 !> \date 2011-07-27 170 subroutine recur_nucpot_destroy() 171 deallocate(zero_pint) 172 end subroutine recur_nucpot_destroy 173 174 !> \brief recursive function for \fn(nucpot_geom) 175 !> \author Bin Gao 176 !> \date 2011-07-27 177 !> \param order_geo contains the orders of xyz components of geometric derivative on 178 !> nuclear potential center 179 !> \param order_boys is the order of Boys function 180 !> \return geo_pint is the geometric derivative on nuclear potential center 181 recursive function recur_nucpot_geom(order_geo, order_boys) result(geo_pint) 182 real(REALK) geo_pint 183 integer, intent(in) :: order_geo(3) 184 integer, intent(in) :: order_boys 185 ! recurrence relation along x-direction 186 if (order_geo(1)>0) then 187 geo_pint = nucorg_wrt_cc(1) & 188 * recur_nucpot_geom((/order_geo(1)-1,order_geo(2:3)/), & 189 order_boys+1) & 190 + real(order_geo(1)-1,REALK) & 191 * recur_nucpot_geom((/order_geo(1)-2,order_geo(2:3)/), & 192 order_boys+1) 193 ! recurrence relation along y-direction 194 else if (order_geo(2)>0) then 195 geo_pint = nucorg_wrt_cc(2) & 196 * recur_nucpot_geom((/order_geo(1),order_geo(2)-1,order_geo(3)/), & 197 order_boys+1) & 198 + real(order_geo(2)-1,REALK) & 199 * recur_nucpot_geom((/order_geo(1),order_geo(2)-2,order_geo(3)/), & 200 order_boys+1) 201 ! recurrence relation along z-direction 202 else if (order_geo(3)>0) then 203 geo_pint = nucorg_wrt_cc(3) & 204 * recur_nucpot_geom((/order_geo(1:2),order_geo(3)-1/), & 205 order_boys+1) & 206 + real(order_geo(3)-1,REALK) & 207 * recur_nucpot_geom((/order_geo(1:2),order_geo(3)-2/), & 208 order_boys+1) 209 ! zero integral 210 else if (any(order_geo<0)) then 211 geo_pint = 0.0_REALK 212 ! auxiliary integral 213 else 214 if (order_boys<=max_order_aux) then 215 geo_pint = zero_pint(order_boys) 216 else 217 geo_pint = 0.0_REALK !gets rid of compiler's warning 218 call error_stop("recur_nucpot_geom", & 219 "no required order of auxiliary integrals", order_boys) 220 end if 221 end if 222 end function recur_nucpot_geom 223 224 !> \brief recursive function for \fn(nucpot_hket) 225 !> \author Bin Gao 226 !> \date 2011-10-18 227 !> \param order_hket contains the orders of xyz components of HGTO on ket center 228 !> \param order_geo contains the orders of xyz components of geometric derivatives 229 !> \return hket_pint is the primitive integral with requied order of HGTOs on ket center 230 recursive function recur_nucpot_hket(order_hket, order_geo) result(hket_pint) 231 real(REALK) hket_pint 232 integer, intent(in) :: order_hket(3) 233 integer, intent(in) :: order_geo(3) 234 ! recurrence relation along x-direction 235 if (order_hket(1)>0) then 236 hket_pint = cc_wrt_ket(1) & 237 * recur_nucpot_hket((/order_hket(1)-1,order_hket(2:3)/), & 238 order_geo) & 239 + half_neg_rp*(bra_to_ket*real(order_hket(1)-1,REALK) & 240 * recur_nucpot_hket((/order_hket(1)-2,order_hket(2:3)/), & 241 order_geo) & 242 + recur_nucpot_hket((/order_hket(1)-1,order_hket(2:3)/), & 243 (/order_geo(1)+1,order_geo(2:3)/))) 244 ! recurrence relation along y-direction 245 else if (order_hket(2)>0) then 246 hket_pint = cc_wrt_ket(2) & 247 * recur_nucpot_hket((/order_hket(1),order_hket(2)-1,order_hket(3)/), & 248 order_geo) & 249 + half_neg_rp*(bra_to_ket*real(order_hket(2)-1,REALK) & 250 * recur_nucpot_hket((/order_hket(1),order_hket(2)-2,order_hket(3)/), & 251 order_geo) & 252 + recur_nucpot_hket((/order_hket(1),order_hket(2)-1,order_hket(3)/), & 253 (/order_geo(1),order_geo(2)+1,order_geo(3)/))) 254 ! recurrence relation along z-direction 255 else if (order_hket(3)>0) then 256 hket_pint = cc_wrt_ket(3) & 257 * recur_nucpot_hket((/order_hket(1:2),order_hket(3)-1/), & 258 order_geo) & 259 + half_neg_rp*(bra_to_ket*real(order_hket(3)-1,REALK) & 260 * recur_nucpot_hket((/order_hket(1:2),order_hket(3)-2/), & 261 order_geo) & 262 + recur_nucpot_hket((/order_hket(1:2),order_hket(3)-1/), & 263 (/order_geo(1:2),order_geo(3)+1/))) 264 ! zero integral 265 else if (any(order_hket<0)) then 266 hket_pint = 0.0_REALK 267 ! recurrence relations on the order of geometric derivatives 268 else 269 hket_pint = recur_nucpot_geom(order_geo, 0) 270 end if 271 end function recur_nucpot_hket 272 273 !> \brief recursive function for \fn(nucpot_hbra) 274 !> \author Bin Gao 275 !> \date 2011-10-18 276 !> \param order_hbra contains the orders of xyz components of HGTO on bra center 277 !> \param order_hket contains the orders of xyz components of HGTO on ket center 278 !> \param order_geo contains the orders of xyz components of geometric derivatives 279 !> \return hbra_pint is the primitive integral with requied order of HGTOs on bra center 280 recursive function recur_nucpot_hbra(order_hbra, order_hket, order_geo) & 281 result(hbra_pint) 282 real(REALK) hbra_pint 283 integer, intent(in) :: order_hbra(3) 284 integer, intent(in) :: order_hket(3) 285 integer, intent(in) :: order_geo(3) 286 ! recurrence relation along x-direction 287 if (order_hbra(1)>0) then 288 hbra_pint = cc_wrt_bra(1) & 289 * recur_nucpot_hbra((/order_hbra(1)-1,order_hbra(2:3)/), & 290 order_hket, order_geo) & 291 + half_neg_rp*(ket_to_bra*real(order_hbra(1)-1,REALK) & 292 * recur_nucpot_hbra((/order_hbra(1)-2,order_hbra(2:3)/), & 293 order_hket, order_geo) & 294 + recur_nucpot_hbra((/order_hbra(1)-1,order_hbra(2:3)/), & 295 order_hket, & 296 (/order_geo(1)+1,order_geo(2:3)/)) & 297 - real(order_hket(1),REALK) & 298 * recur_nucpot_hbra((/order_hbra(1)-1,order_hbra(2:3)/), & 299 (/order_hket(1)-1,order_hket(2:3)/), & 300 order_geo)) 301 ! recurrence relation along y-direction 302 else if (order_hbra(2)>0) then 303 hbra_pint = cc_wrt_bra(2) & 304 * recur_nucpot_hbra((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 305 order_hket, order_geo) & 306 + half_neg_rp*(ket_to_bra*real(order_hbra(2)-1,REALK) & 307 * recur_nucpot_hbra((/order_hbra(1),order_hbra(2)-2,order_hbra(3)/), & 308 order_hket, order_geo) & 309 + recur_nucpot_hbra((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 310 order_hket, & 311 (/order_geo(1),order_geo(2)+1,order_geo(3)/)) & 312 - real(order_hket(2),REALK) & 313 * recur_nucpot_hbra((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 314 (/order_hket(1),order_hket(2)-1,order_hket(3)/), & 315 order_geo)) 316 ! recurrence relation along z-direction 317 else if (order_hbra(3)>0) then 318 hbra_pint = cc_wrt_bra(3) & 319 * recur_nucpot_hbra((/order_hbra(1:2),order_hbra(3)-1/), & 320 order_hket, order_geo) & 321 + half_neg_rp*(ket_to_bra*real(order_hbra(3)-1,REALK) & 322 * recur_nucpot_hbra((/order_hbra(1:2),order_hbra(3)-2/), & 323 order_hket, order_geo) & 324 + recur_nucpot_hbra((/order_hbra(1:2),order_hbra(3)-1/), & 325 order_hket, & 326 (/order_geo(1:2),order_geo(3)+1/)) & 327 - real(order_hket(3),REALK) & 328 * recur_nucpot_hbra((/order_hbra(1:2),order_hbra(3)-1/), & 329 (/order_hket(1:2),order_hket(3)-1/), & 330 order_geo)) 331 ! zero integral 332 else if (any(order_hbra<0) .or. any(order_hket<0)) then 333 hbra_pint = 0.0_REALK 334 ! recurrence relations on the order of HGTOs on ket center 335 else 336 hbra_pint = recur_nucpot_hket(order_hket, order_geo) 337 end if 338 end function recur_nucpot_hbra 339 340 !> \brief recursive function for \fn(nucpot_moment) 341 !> \author Bin Gao 342 !> \date 2011-10-18 343 !> \param order_hbra contains the orders of xyz components of HGTO on bra center 344 !> \param order_hket contains the orders of xyz components of HGTO on ket center 345 !> \param order_geo contains the orders of xyz components of geometric derivatives 346 !> \param order_mom contains the orders of xyz components of Cartesian multipole moment 347 !> \return hmom_pint is the nuclear attraction integral with required orders of Cartesian 348 !> multipole moment integral and HGTOs on bra and ket centers 349 recursive function recur_nucpot_moment(order_hbra, order_hket, order_mom, order_geo) & 350 result(hmom_pint) 351 real(REALK) hmom_pint 352 integer, intent(in) :: order_hbra(3) 353 integer, intent(in) :: order_hket(3) 354 integer, intent(in) :: order_mom(3) 355 integer, intent(in) :: order_geo(3) 356 ! recurrence relation along x-direction 357 if (order_mom(1)>0) then 358 hmom_pint = cc_wrt_diporg(1) & 359 * recur_nucpot_moment(order_hbra, order_hket, & 360 (/order_mom(1)-1,order_mom(2:3)/), & 361 order_geo) & 362 + half_recip_p*(real(order_hbra(1),REALK) & 363 * recur_nucpot_moment((/order_hbra(1)-1,order_hbra(2:3)/), & 364 order_hket, & 365 (/order_mom(1)-1,order_mom(2:3)/), & 366 order_geo) & 367 + real(order_hket(1),REALK) & 368 * recur_nucpot_moment(order_hbra, & 369 (/order_hket(1)-1,order_hket(2:3)/), & 370 (/order_mom(1)-1,order_mom(2:3)/), & 371 order_geo) & 372 + real(order_mom(1)-1,REALK) & 373 * recur_nucpot_moment(order_hbra, order_hket, & 374 (/order_mom(1)-2,order_mom(2:3)/), & 375 order_geo) & 376 - recur_nucpot_moment(order_hbra, order_hket, & 377 (/order_mom(1)-1,order_mom(2:3)/), & 378 (/order_geo(1)+1,order_geo(2:3)/))) 379 ! recurrence relation along y-direction 380 else if (order_mom(2)>0) then 381 hmom_pint = cc_wrt_diporg(2) & 382 * recur_nucpot_moment(order_hbra, order_hket, & 383 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 384 order_geo) & 385 + half_recip_p*(real(order_hbra(2),REALK) & 386 * recur_nucpot_moment((/order_hbra(1),order_hbra(2)-1,order_hbra(3)/), & 387 order_hket, & 388 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 389 order_geo) & 390 + real(order_hket(2),REALK) & 391 * recur_nucpot_moment(order_hbra, & 392 (/order_hket(1),order_hket(2)-1,order_hket(3)/), & 393 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 394 order_geo) & 395 + real(order_mom(2)-1,REALK) & 396 * recur_nucpot_moment(order_hbra, order_hket, & 397 (/order_mom(1),order_mom(2)-2,order_mom(3)/), & 398 order_geo) & 399 - recur_nucpot_moment(order_hbra, order_hket, & 400 (/order_mom(1),order_mom(2)-1,order_mom(3)/), & 401 (/order_geo(1),order_geo(2)+1,order_geo(3)/))) 402 ! recurrence relation along z-direction 403 else if (order_mom(3)>0) then 404 hmom_pint = cc_wrt_diporg(3) & 405 * recur_nucpot_moment(order_hbra, order_hket, & 406 (/order_mom(1:2),order_mom(3)-1/), & 407 order_geo) & 408 + half_recip_p*(real(order_hbra(3),REALK) & 409 * recur_nucpot_moment((/order_hbra(1:2),order_hbra(3)-1/), & 410 order_hket, & 411 (/order_mom(1:2),order_mom(3)-1/), & 412 order_geo) & 413 + real(order_hket(3),REALK) & 414 * recur_nucpot_moment(order_hbra, & 415 (/order_hket(1:2),order_hket(3)-1/), & 416 (/order_mom(1:2),order_mom(3)-1/), & 417 order_geo) & 418 + real(order_mom(3)-1,REALK) & 419 * recur_nucpot_moment(order_hbra, order_hket, & 420 (/order_mom(1:2),order_mom(3)-2/), & 421 order_geo) & 422 - recur_nucpot_moment(order_hbra, order_hket, & 423 (/order_mom(1:2),order_mom(3)-1/), & 424 (/order_geo(1:2),order_geo(3)+1/))) 425 ! zero integral 426 else if (any(order_hbra<0) .or. any(order_hket<0) .or. any(order_mom<0)) then 427 hmom_pint = 0.0_REALK 428 ! recurrence relations on the order of HGTOs on bra center 429 else 430 hmom_pint = recur_nucpot_hbra(order_hbra, order_hket, order_geo) 431 end if 432 end function recur_nucpot_moment 433 434end module recur_nucpot 435