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