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!! This file calculates the Cartesian multipole moment integrals 18!! using contracted London Cartesian Gaussians at zero fields. 19!! 20!! 2011-03-24, Bin Gao: 21!! * first version 22 23 !> \brief calculates the Cartesian multipole moment integrals using 24 !> contracted London Cartesian Gaussians at zero fields 25 !> \author Bin Gao 26 !> \date 2011-03-24 27 !> \param idx_bra is the atomic index of bra center 28 !> \param coord_bra contains the coordinates of bra center 29 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 30 !> \param num_prim_bra is the number of primitive Gaussians of bra center 31 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 32 !> \param num_contr_bra is the number of contractions of bra center 33 !> \param contr_coef_bra contains the contraction coefficients of bra center 34 !> \param idx_ket is the atomic index of ket center 35 !> \param coord_ket contains the coordinates of ket center 36 !> \param angular_ket is the angular number of ket center 37 !> \param num_prim_ket is the number of primitive Gaussians of ket center 38 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 39 !> \param num_contr_ket is the number of contractions of ket center 40 !> \param contr_coef_ket contains the contraction coefficients of ket center 41 !> \param dipole_origin contains the coordinates of dipole origin 42 !> \param order_mom is the order of Cartesian multipole moments 43 !> \param order_elec is the order of electronic derivatives 44 !> \param order_mag is the order of magnetic derivatives 45 !> \param num_cents is the number of differentiated centers 46 !> \param idx_cent contains the indices of atomic centers 47 !> \param order_cent contains the order of derivatives of the corresponding atomic centers 48 !> \param num_cart_bra is the number of Cartesian GTOs on bra, 49 !> equals to \f$(\var(angular_bra)+1)(\var(angular_bra)+2)/2\f$ 50 !> \param num_cart_ket is the number of Cartesian GTOs on ket, 51 !> equals to \f$(\var(angular_ket)+1)(\var(angular_ket)+2)/2\f$ 52 !> \param num_opt is the number of operators including derivatives 53 !> \return contr_ints contains the contracted Cartesian multipole moment integrals 54 subroutine lcgto_zero_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 55 exponent_bra, num_contr_bra, contr_coef_bra, & 56 idx_ket, coord_ket, angular_ket, num_prim_ket, & 57 exponent_ket, num_contr_ket, contr_coef_ket, & 58 order_mag_bra, order_mag_ket, order_mag_total, & 59 order_ram_bra, order_ram_ket, order_ram_total, & 60 order_geo_bra, order_geo_ket, & 61 num_cents, idx_cent, order_cent, & 62 idx_diporg, dipole_origin, scal_const, & 63 order_geo_mom, order_mom, order_elec, & 64 num_cart_bra, num_cart_ket, num_opt, contr_ints) 65 implicit none 66 integer, intent(in) :: idx_bra 67 real(8), intent(in) :: coord_bra(3) 68 integer, intent(in) :: angular_bra 69 integer, intent(in) :: num_prim_bra 70 real(8), intent(in) :: exponent_bra(num_prim_bra) 71 integer, intent(in) :: num_contr_bra 72 real(8), intent(in) :: contr_coef_bra(num_prim_bra,num_contr_bra) 73 integer, intent(in) :: idx_ket 74 real(8), intent(in) :: coord_ket(3) 75 integer, intent(in) :: angular_ket 76 integer, intent(in) :: num_prim_ket 77 real(8), intent(in) :: exponent_ket(num_prim_ket) 78 integer, intent(in) :: num_contr_ket 79 real(8), intent(in) :: contr_coef_ket(num_prim_ket,num_contr_ket) 80 real(8), intent(in) :: dipole_origin(3) 81 integer, intent(in) :: order_mom 82 integer, intent(in) :: order_elec 83 integer, intent(in) :: order_mag 84 integer, intent(in) :: num_cents 85 integer, intent(in) :: idx_cent(2) 86 integer, intent(in) :: order_cent(2) 87 integer, intent(in) :: num_cints 88 real(8), intent(out) :: contr_ints(num_cints) 89!f2py intent(in) :: idx_bra 90!f2py intent(in) :: coord_bra 91!f2py intent(in) :: angular_bra 92!f2py intent(hide) :: num_prim_bra 93!f2py intent(in) :: exponent_bra 94!f2py intent(hide) :: num_contr_bra 95!f2py intent(in) :: contr_coef_bra 96!f2py depend(num_prim_bra) :: contr_coef_bra 97!f2py intent(in) :: idx_ket 98!f2py intent(in) :: coord_ket 99!f2py intent(in) :: angular_ket 100!f2py intent(hide) :: num_prim_ket 101!f2py intent(in) :: exponent_ket 102!f2py intent(hide) :: num_contr_ket 103!f2py intent(in) :: contr_coef_ket 104!f2py depend(num_prim_ket) :: contr_coef_ket 105!f2py intent(in) :: dipole_origin 106!f2py intent(in) :: order_mom 107!f2py intent(in) :: order_elec 108!f2py intent(in) :: order_mag 109!f2py intent(in) :: num_cents 110!f2py intent(in) :: idx_cent 111!f2py intent(in) :: order_cent 112!f2py intent(in) :: num_cints 113!f2py intent(out) :: contr_ints 114!f2py depend(num_cints) :: contr_ints 115#if defined(DEBUG) 116 ! dumps the contracted GTOs, operators and derivatives to check 117 call dump_int_info("LCGTO", idx_bra, coord_bra, angular_bra, & 118 num_prim_bra, exponent_bra, & 119 num_contr_bra, contr_coef_bra, & 120 "LCGTO", idx_ket, coord_ket, angular_ket, & 121 num_prim_ket, exponent_ket, & 122 num_contr_ket, contr_coef_ket, & 123 "carmom", 1, (/0/), dipole_origin, & 124 (/order_mom/), order_elec, order_mag, & 125 2, num_cents, idx_cent, order_cent) 126#endif 127 ! too many differentiated centers, or the centers of bra and ket coincide 128 if (num_cents>2 .or. (num_cents>0 .and. idx_bra==idx_ket)) then 129 contr_ints = 0.0D+00 130 else 131 select case(num_cents) 132 ! no geometric derivatives 133 case(0) 134 call lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 135 num_prim_bra, exponent_bra, & 136 num_contr_bra, contr_coef_bra, & 137 coord_ket, angular_ket, & 138 num_prim_ket, exponent_ket, & 139 num_contr_ket, contr_coef_ket, & 140 dipole_origin, order_mom, & 141 order_elec, order_mag, & 142 0, 0, num_cints, contr_ints) 143 ! one-center geometric derivatives 144 case(1) 145 if (idx_cent(1)==idx_bra) then 146 call lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 147 num_prim_bra, exponent_bra, & 148 num_contr_bra, contr_coef_bra, & 149 coord_ket, angular_ket, & 150 num_prim_ket, exponent_ket, & 151 num_contr_ket, contr_coef_ket, & 152 dipole_origin, order_mom, & 153 order_elec, order_mag, & 154 order_cent(1), 0, num_cints, contr_ints) 155 else if (idx_cent(1)==idx_ket) then 156 call lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 157 num_prim_bra, exponent_bra, & 158 num_contr_bra, contr_coef_bra, & 159 coord_ket, angular_ket, & 160 num_prim_ket, exponent_ket, & 161 num_contr_ket, contr_coef_ket, & 162 dipole_origin, order_mom, & 163 order_elec, order_mag, & 164 0, order_cent(1), num_cints, contr_ints) 165 else 166 contr_ints = 0.0D+00 167 end if 168 ! two-center geometric derivatives 169 case(2) 170 if (idx_cent(1)==idx_bra .and. idx_cent(2)==idx_ket) then 171 call lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 172 num_prim_bra, exponent_bra, & 173 num_contr_bra, contr_coef_bra, & 174 coord_ket, angular_ket, & 175 num_prim_ket, exponent_ket, & 176 num_contr_ket, contr_coef_ket, & 177 dipole_origin, order_mom, & 178 order_elec, order_mag, & 179 order_cent(1), order_cent(2), & 180 num_cints, contr_ints) 181 else if (idx_cent(2)==idx_bra .and. idx_cent(1)==idx_ket) then 182 call lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 183 num_prim_bra, exponent_bra, & 184 num_contr_bra, contr_coef_bra, & 185 coord_ket, angular_ket, & 186 num_prim_ket, exponent_ket, & 187 num_contr_ket, contr_coef_ket, & 188 dipole_origin, order_mom, & 189 order_elec, order_mag, & 190 order_cent(2), order_cent(1), & 191 num_cints, contr_ints) 192 else 193 contr_ints = 0.0D+00 194 end if 195 end select 196 end if 197 return 198 end subroutine lcgto_zero_carmom 199 200 !> \brief recurrence relations of Cartesian multipole moment integrals using 201 !> contracted London Cartesian Gaussians at zero field 202 !> \author Bin Gao 203 !> \date 2011-03-24 204 !> \param coord_bra contains the coordinates of bra center 205 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 206 !> \param num_prim_bra is the number of primitive Gaussians of bra center 207 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 208 !> \param num_contr_bra is the number of contractions of bra center 209 !> \param contr_coef_bra contains the contraction coefficients of bra center 210 !> \param coord_ket contains the coordinates of ket center 211 !> \param angular_ket is the angular number of ket center 212 !> \param num_prim_ket is the number of primitive Gaussians of ket center 213 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 214 !> \param num_contr_ket is the number of contractions of ket center 215 !> \param contr_coef_ket contains the contraction coefficients of ket center 216 !> \param dipole_origin contains the coordinates of dipole origin 217 !> \param order_mom is the order of Cartesian multipole moments 218 !> \param order_elec is the order of electronic derivatives 219 !> \param order_mag is the order of magnetic derivatives 220 !> \param order_geo_bra is the order of geometric derivatives with respect to the bra center 221 !> \param order_geo_ket is the order of geometric derivatives with respect to the ket center 222 !> \param num_cints is the number of contracted Cartesian multipole moment integrals 223 !> \return contr_ints contains the Cartesian multipole moment integrals 224 subroutine lcgto_zero_carmom_recurr(coord_bra, angular_bra, & 225 num_prim_bra, exponent_bra, & 226 num_contr_bra, contr_coef_bra, & 227 coord_ket, angular_ket, & 228 num_prim_ket, exponent_ket, & 229 num_contr_ket, contr_coef_ket, & 230 dipole_origin, order_mom, & 231 order_elec, & 232 order_geo_bra, order_geo_ket, & 233 order_mag_bra, order_mag_ket, order_ram_bra, order_ram_ket, 234 num_cints, contr_ints) 235 implicit none 236 real(8), intent(in) :: coord_bra(3) 237 integer, intent(in) :: angular_bra 238 integer, intent(in) :: num_prim_bra 239 real(8), intent(in) :: exponent_bra(num_prim_bra) 240 integer, intent(in) :: num_contr_bra 241 real(8), intent(in) :: contr_coef_bra(num_prim_bra,num_contr_bra) 242 real(8), intent(in) :: coord_ket(3) 243 integer, intent(in) :: angular_ket 244 integer, intent(in) :: num_prim_ket 245 real(8), intent(in) :: exponent_ket(num_prim_ket) 246 integer, intent(in) :: num_contr_ket 247 real(8), intent(in) :: contr_coef_ket(num_prim_ket,num_contr_ket) 248 real(8), intent(in) :: dipole_origin(3) 249 integer, intent(in) :: order_mom 250 integer, intent(in) :: order_elec 251 integer, intent(in) :: order_geo_bra 252 integer, intent(in) :: order_geo_ket 253 integer, intent(in) :: order_mag 254 integer, intent(in) :: order_mag_bra 255 integer, intent(in) :: order_mag_ket 256 integer, intent(in) :: order_ram 257 integer, intent(in) :: order_ram_bra 258 integer, intent(in) :: order_ram_ket 259 integer, intent(in) :: num_cints 260 real(8), intent(out) :: contr_ints(num_cints) 261!f2py intent(in) :: coord_bra 262!f2py intent(in) :: angular_bra 263!f2py intent(hide) :: num_prim_bra 264!f2py intent(in) :: exponent_bra 265!f2py intent(hide) :: num_contr_bra 266!f2py intent(in) :: contr_coef_bra 267!f2py depend(num_prim_bra) :: contr_coef_bra 268!f2py intent(in) :: coord_ket 269!f2py intent(in) :: angular_ket 270!f2py intent(hide) :: num_prim_ket 271!f2py intent(in) :: exponent_ket 272!f2py intent(hide) :: num_contr_ket 273!f2py intent(in) :: contr_coef_ket 274!f2py depend(num_prim_ket) :: contr_coef_ket 275!f2py intent(in) :: dipole_origin 276!f2py intent(in) :: order_mom 277!f2py intent(in) :: order_elec 278!f2py intent(in) :: order_mag 279!f2py intent(in) :: order_geo_bra 280!f2py intent(in) :: order_geo_ket 281!f2py intent(in) :: num_cints 282!f2py intent(out) :: contr_ints 283!f2py depend(num_cints) :: contr_ints 284 integer iprim, jprim !incremental recorders over primitives 285 integer icontr, jcontr !incremental recorders over contractions 286 integer ierr !error information 287 if (order_mom>0) then 288 if (order_ram>0) then 289 else 290 end if 291 else 292 if (order_ram>0) then 293 do iram = 0, order_ram 294 end do 295 else 296 call contr_cgto_carmom_recurr(coord_bra, angular_bra, & 297 num_prim_bra, exponent_bra, & 298 num_contr_bra, contr_coef_bra, & 299 coord_ket, angular_ket, & 300 num_prim_ket, exponent_ket, & 301 num_contr_ket, contr_coef_ket, & 302 order_elec, dipole_origin, & 303 order_mom, order_geo_mom, & 304 order_geo_bra, order_geo_ket, & 305 num_cart_bra, num_cart_ket, & 306 num_opt, contr_ints) 307 end if 308 end if 309 return 310 end subroutine lcgto_zero_carmom_recurr 311 312 !> \brief recurrence relations of Cartesian multipole moment integrals using 313 !> contracted Cartesian Gaussians 314 !> \author Bin Gao 315 !> \date 2011-03-24 316 !> \param coord_bra contains the coordinates of bra center 317 !> \param angular_bra contains the minimum and maximum angular numbers of bra center (s=0, p=1, d=2, ...) 318 !> \param num_prim_bra is the number of primitive Gaussians of bra center 319 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 320 !> \param num_contr_bra is the number of contractions of bra center 321 !> \param contr_coef_bra contains the contraction coefficients of bra center 322 !> \param coord_ket contains the coordinates of ket center 323 !> \param angular_ket contains the minimum and maximum angular number of ket center 324 !> \param num_prim_ket is the number of primitive Gaussians of ket center 325 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 326 !> \param num_contr_ket is the number of contractions of ket center 327 !> \param contr_coef_ket contains the contraction coefficients of ket center 328 !> \param order_elec is the order of electronic derivatives 329 !> \param dipole_origin contains the coordinates of dipole origin 330 !> \param order_mom is the order of Cartesian multipole moments 331 !> \param order_geo_mom is the order of geometric derivatives on Cartesian multipole moments 332 !> \param order_geo_bra contains the minimum and maximum orders of geometric derivatives 333 !> with respect to the bra center 334 !> \param order_geo_ket contains the minimum and maximum orders of geometric derivatives 335 !> with respect to the ket center 336 !> \param dim_cart_bra is the number of Cartesian GTOs on bra 337 !> \param dim_cart_ket is the number of Cartesian GTOs on ket 338 !> \param num_opt is the number of operators including derivatives 339 !> \return contr_ints contains the contracted Cartesian multipole moment integrals 340 subroutine contr_cgto_carmom_recurr(coord_bra, angular_bra, & 341 num_prim_bra, exponent_bra, & 342 num_contr_bra, contr_coef_bra, & 343 coord_ket, angular_ket, & 344 num_prim_ket, exponent_ket, & 345 num_contr_ket, contr_coef_ket, & 346 order_elec, dipole_origin, & 347 order_mom, order_geo_mom, & 348 order_geo_bra, order_geo_ket, & 349 dim_cart_bra, dim_cart_ket, & 350 num_opt, contr_ints) 351 implicit none 352 real(8), intent(in) :: coord_bra(3) 353 integer, intent(in) :: angular_bra(2) 354 integer, intent(in) :: num_prim_bra 355 real(8), intent(in) :: exponent_bra(num_prim_bra) 356 integer, intent(in) :: num_contr_bra 357 real(8), intent(in) :: contr_coef_bra(num_prim_bra,num_contr_bra) 358 real(8), intent(in) :: coord_ket(3) 359 integer, intent(in) :: angular_ket(2) 360 integer, intent(in) :: num_prim_ket 361 real(8), intent(in) :: exponent_ket(num_prim_ket) 362 integer, intent(in) :: num_contr_ket 363 real(8), intent(in) :: contr_coef_ket(num_prim_ket,num_contr_ket) 364 integer, intent(in) :: order_elec 365 real(8), intent(in) :: dipole_origin(3) 366 integer, intent(in) :: order_mom 367 integer, intent(in) :: order_geo_mom 368 integer, intent(in) :: order_geo_bra(2) 369 integer, intent(in) :: order_geo_ket(2) 370 integer, intent(in) :: dim_cart_bra 371 integer, intent(in) :: dim_cart_ket 372 integer, intent(in) :: num_opt 373 real(8), intent(out) :: contr_ints(dim_cart_bra,num_contr_bra, & 374 dim_cart_ket,num_contr_ket,num_opt) 375!f2py intent(in) :: coord_bra 376!f2py intent(in) :: angular_bra 377!f2py intent(hide) :: num_prim_bra 378!f2py intent(in) :: exponent_bra 379!f2py intent(hide) :: num_contr_bra 380!f2py intent(in) :: contr_coef_bra 381!f2py depend(num_prim_bra) :: contr_coef_bra 382!f2py intent(in) :: coord_ket 383!f2py intent(in) :: angular_ket 384!f2py intent(hide) :: num_prim_ket 385!f2py intent(in) :: exponent_ket 386!f2py intent(hide) :: num_contr_ket 387!f2py intent(in) :: contr_coef_ket 388!f2py depend(num_prim_ket) :: contr_coef_ket 389!f2py intent(in) :: order_elec 390!f2py intent(in) :: dipole_origin 391!f2py intent(in) :: order_mom 392!f2py intent(in) :: order_geo_mom 393!f2py intent(in) :: order_geo_bra 394!f2py intent(in) :: order_geo_ket 395!f2py intent(in) :: dim_cart_bra 396!f2py intent(in) :: dim_cart_ket 397!f2py intent(in) :: num_opt 398!f2py intent(out) :: contr_ints 399!f2py depend(dim_cart_bra) :: contr_ints 400!f2py depend(num_contr_bra) :: contr_ints 401!f2py depend(dim_cart_ket) :: contr_ints 402!f2py depend(num_contr_ket) :: contr_ints 403!f2py depend(num_opt) :: contr_ints 404 integer, parameter :: STDOUT = 6 !IO of standard output 405 integer order_herm_bra(2) !range of orders of Hermite Gaussians on bra 406 integer order_herm_ket(2) !range of orders of Hermite Gaussians on ket 407 integer low_order_mom !order of Cartesian multipole moment by considering its derivative 408 integer dim_herm_bra !number of Hermite Gaussians on bra 409 integer dim_herm_ket !number of Hermite Gaussians on ket 410 integer num_elec !number of xyz components of electronic derivatives 411 integer num_low_mom !number of xyz components of lower order Cartesian multipole moment 412 integer num_uniq_ops !number of unique operators after considering the geometric 413 !derivative of Cartesian multipole moment 414 integer num_geo_bra !number of geometric derivatives on bra 415 integer num_geo_ket !number of geometric derivatives on ket 416 real(8), allocatable :: herm_pints(:,:,:) !primitive Hermite integrals 417 real(8), allocatable :: herm_ket_pints(:,:,:,:) !primitive Hermite (ket) integrals 418 integer prod_op_geo !product of number of derivatives and operators 419 real(8), allocatable :: cart_pints(:,:,:,:,:) !primitive Cartesian integrals 420 integer dim_contr_cgto !dimension of contracted CGTOs 421 integer start_opt !start address of the first operator 422 integer iprim, jprim !incremental recorders over primitives 423 integer icontr, jcontr !incremental recorders over contractions 424 integer icart, jcart !incremental recorders over Cartesian Gaussians 425 integer iop !incremental recorder over operators 426 integer ierr !error information 427 ! computes the minimum and maximum orders of Hermite Gaussians on bra and ket 428 order_herm_bra(1) = order_geo_bra(1) !FIXME +mod(angular_bra,2) 429 order_herm_bra(2) = order_geo_bra(2)+angular_bra(2) 430 order_herm_ket(1) = order_geo_ket(1) !FIXME +mod(angular_ket,2) 431 order_herm_ket(2) = order_geo_ket(2)+angular_ket(2) 432 ! computes the order of Cartesian multipole moment after considering its geometric derivative 433 low_order_mom = order_mom-order_geo_mom 434 ! computes the number of primitive Hermite Gaussians used in recurrence relations 435 if (order_herm_bra(1)==0) then 436 dim_herm_bra = (order_herm_bra(2)+1)*(order_herm_bra(2)+2)*(order_herm_bra(2)+3)/6 437 else 438 dim_herm_bra = ((order_herm_bra(2)+1)*(order_herm_bra(2)+2)*(order_herm_bra(2)+3) & 439 - order_herm_bra(1)*(order_herm_bra(1)+1)*(order_herm_bra(1)+2))/6 440 end if 441 if (order_herm_ket(1)==0) then 442 dim_herm_ket = (order_herm_ket(2)+1)*(order_herm_ket(2)+2)*(order_herm_ket(2)+3)/6 443 else 444 dim_herm_ket = ((order_herm_ket(2)+1)*(order_herm_ket(2)+2)*(order_herm_ket(2)+3) & 445 - order_herm_ket(1)*(order_herm_ket(1)+1)*(order_herm_ket(1)+2))/6 446 end if 447 ! computes the number of operators, which is more efficient than using \var(TRI_SIZE) 448 num_elec = (order_elec+1)*(order_elec+2)/2 449 num_low_mom = (low_order_mom+1)*(low_order_mom+2)/2 450 num_uniq_ops = num_elec*num_low_mom 451 ! computes the number of geometric derivatives 452 num_geo_bra = (order_geo_bra+1)*(order_geo_bra+2)/2 453 num_geo_ket = (order_geo_ket+1)*(order_geo_ket+2)/2 454 ! allocates the memory for primitive integrals 455 allocate(herm_pints(dim_herm_bra,dim_herm_ket,num_uniq_ops), stat=ierr) 456 if (ierr/=0) then 457 write(STDOUT,100) "contr_cgto_carmom_recurr: failed to allocate", & 458 dim_herm_bra*dim_herm_ket*num_uniq_ops, & 459 "primitive Hermite integrals!" 460 stop 461 end if 462 allocate(herm_ket_pints(dim_cart_bra,dim_herm_ket,num_uniq_ops,num_geo_bra), & 463 stat=ierr) 464 if (ierr/=0) then 465 write(STDOUT,100) "contr_cgto_carmom_recurr: failed to allocate", & 466 dim_cart_bra*dim_herm_ket*num_uniq_ops*num_geo_bra, & 467 "primitive Cartesian-Hermite integrals!" 468 stop 469 end if 470 prod_op_geo = num_uniq_ops*num_geo_bra*num_geo_ket 471 ! computes the start address of the first operator 472 if (order_geo_mom==0) then 473 start_opt = 1 474 else 475 dim_contr_cgto = dim_cart_bra*num_contr_bra*dim_cart_ket*num_contr_ket 476 start_opt = num_opt-prod_op_geo+1 477 end if 478 allocate(cart_pints(dim_cart_bra,dim_cart_ket,start_opt:num_opt, & 479 num_prim_bra,num_prim_ket), stat=ierr) 480 if (ierr/=0) then 481 write(STDOUT,100) "contr_cgto_carmom_recurr: failed to allocate", & 482 dim_cart_bra*dim_cart_ket*prod_op_geo*num_prim_bra*num_prim_ket, & 483 "primitive Cartesian integrals!" 484 stop 485 end if 486#if defined(DEBUG) 487 write(STDOUT,100) "Number of primitive Hermite integrals:", & 488 dim_herm_bra*dim_herm_ket*num_uniq_ops 489 write(STDOUT,100) "Number of primitive Cartesian-Hermite integrals:", & 490 dim_cart_bra*dim_herm_ket*num_uniq_ops*num_geo_bra 491 write(STDOUT,100) "Number of primitive Cartesian integrals:", & 492 dim_cart_bra*dim_cart_ket*prod_op_geo*num_prim_bra*num_prim_ket 493#endif 494 ! computes the primitive integrals 495 do jprim = 1, num_prim_ket 496 do iprim = 1, num_prim_bra 497 ! computes the primitive Hermite integrals 498 call prim_hgto_carmom(order_herm_bra, coord_bra, exponent_bra(iprim), & 499 order_herm_ket, coord_ket, exponent_ket(jprim), & 500 order_elec, dipole_origin, low_order_mom, & 501 dim_herm_bra, dim_herm_ket, num_elec, num_low_mom, & 502 herm_pints) 503 ! transforms Hermite Gaussians on bra to Cartesian ones first 504 call hgto_to_lcgto_ket(angular_bra, & 505 (/order_geo_bra,order_geo_bra/), & 506 exponent_bra(iprim), 1, dim_herm_bra, & 507 num_uniq_ops, dim_herm_ket, herm_pints, & 508 dim_cart_bra, num_geo_bra, herm_ket_pints) 509 ! then Hermite Gaussians on ket 510 call hgto_to_lcgto_ket(angular_ket, & 511 (/order_geo_ket,order_geo_ket/), & 512 exponent_ket(jprim), dim_cart_bra, & 513 dim_herm_ket, num_uniq_ops, num_geo_bra, & 514 herm_ket_pints, dim_cart_ket, num_geo_ket, & 515 cart_pints(:,:,:,iprim,jprim)) 516 end do 517 end do 518 ! cleans 519 deallocate(herm_pints) 520 deallocate(herm_ket_pints) 521 ! constructs contracted integrals (without geometric derivative on Cartesian multipole moment), 522 ! in the order of contr_ints(dim_cart_bra,num_contr_bra,dim_cart_ket,num_contr_ket, & 523 ! num_elec,num_low_mom,num_geo_bra,num_geo_ket) 524 contr_ints(:,:,:,:,start_opt:num_opt) = 0.0 525 do jprim = 1, num_prim_ket 526 do iprim = 1, num_prim_bra 527 do iop = start_opt, num_opt 528 do jcontr = 1, num_contr_ket 529 do jcart = 1, dim_cart_ket 530 do icontr = 1, num_contr_bra 531 do icart = 1, dim_cart_bra 532 !FIXME: uses intermediate results to reduce the number of loops 533 contr_ints(icart,icontr,jcart,jcontr,iop) & 534 = contr_ints(icart,icontr,jcart,jcontr,iop) & 535 + contr_coef_bra(iprim,icontr) & !FIXME 536 * contr_coef_ket(jprim,jcontr) & !FIXME 537 * cart_pints(icart,jcart,iop,iprim,jprim) 538 end do 539 end do 540 end do 541 end do 542 end do 543 end do 544 end do 545 ! cleans 546 deallocate(cart_pints) 547 ! constructs the geometric derivatives of Cartesian multipole moments 548 ! in the order of contr_ints(dim_cart_bra,num_contr_bra,dim_cart_ket,num_contr_ket, & 549 ! num_elec,num_mom,num_geo_mom,num_geo_bra,num_geo_ket) 550 if (order_geo_mom>0) then 551 call carmom_deriv(order_mom, low_order_mom, order_geo_mom, & 552 .false., dim_contr_cgto*num_elec, & 553 num_low_mom, num_geo_bra*num_geo_ket, & 554 contr_ints(:,:,:,:,start_opt:num_opt), & 555 (order_mom+1)*(order_mom+2)/2, & 556 (order_geo_mom+1)*(order_geo_mom+2)/2, & 557 contr_ints(:,:,:,:,1:start_opt-1)) 558 end if 559 return 560100 format("contr_cgto_carmom_recurr>> ",A,I8,1X,A) 561 end subroutine contr_cgto_carmom_recurr 562