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 recovers the HGTOs on ket center in nuclear attraction potential integrals. 18!! 19!! 2012-03-04, Bin Gao 20!! * rewrites to improve efficiency 21!! 22!! 2011-10-18, Bin Gao 23!! * first version 24 25#include "stdout.h" 26 27 !> \brief recovers the HGTOs on ket center in nuclear attraction potential integrals 28 !> \author Bin Gao 29 !> \date 2011-10-18 30 !> \param orders_hgto_ket contains the minimum and maximum orders of HGTOs to return 31 !> \param orders_geo_pot contains the minimum and maximum orders of geometric derivatives to return 32 !> \param coord_ket contains the coordinates of bra center 33 !> \param exponent_ket is the exponent of HGTOs of bra center 34 !> \param coord_ket contains the coordinates of ket center 35 !> \param exponent_ket is the exponent of HGTOs of ket center 36 !> \param dim_geo_pot is the dimension of geometric derivatives 37 !> \param geo_pot_pints contains the nuclear attraction integrals with zeroth order 38 !> HGTOs on ket center 39 !> \param dim_hgto_ket is the dimension of HGTOs of ket center afterwards 40 !> \param dim_geo_hket is the dimension of geometric derivatives afterwards 41 !> \return hket_pints contains the integrals with required HGTOs on ket center 42 subroutine nucpot_hket(orders_hgto_ket, orders_geo_pot, coord_bra, exponent_bra, & 43 coord_ket, exponent_ket, dim_geo_pot, geo_pot_pints, & 44 dim_hgto_ket, dim_geo_hket, hket_pints) 45 use xkind 46 implicit none 47 integer, intent(in) :: orders_hgto_ket(2) 48 integer, intent(in) :: orders_geo_pot(2) 49 real(REALK), intent(in) :: coord_bra(3) 50 real(REALK), intent(in) :: exponent_bra 51 real(REALK), intent(in) :: coord_ket(3) 52 real(REALK), intent(in) :: exponent_ket 53 integer, intent(in) :: dim_geo_pot 54 real(REALK), intent(in) :: geo_pot_pints(dim_geo_pot) 55 integer, intent(in) :: dim_hgto_ket 56 integer, intent(in) :: dim_geo_hket 57 real(REALK), intent(out) :: hket_pints(dim_hgto_ket,dim_geo_hket) 58!f2py intent(in) :: orders_hgto_ket 59!f2py intent(in) :: orders_geo_pot 60!f2py intent(in) :: coord_bra 61!f2py intent(in) :: exponent_bra 62!f2py intent(in) :: coord_ket 63!f2py intent(in) :: exponent_ket 64!f2py intent(hide) :: dim_geo_pot 65!f2py intent(in) :: geo_pot_pints 66!f2py intent(in) :: dim_hgto_ket 67!f2py intent(in) :: dim_geo_hket 68!f2py intent(out) :: hket_pints 69!f2py depend(dim_hgto_ket) :: hket_pints 70!f2py depend(dim_geo_hket) :: hket_pints 71 logical zero_hket !if returning zeroth order HGTO on ket center 72 integer max_low_geo !maximum of lower order of geometric derivatives 73 integer max_cur_hket !maximum of current order of HGTOs 74 real(REALK) half_neg_rp !half of the negative reciprocal of total exponent \f$p_{ij}\f$ 75 real(REALK) cc_wrt_ket(3) !relative coordinates of center-of-charge w.r.t. ket center 76 real(REALK) bra_to_ket !ratio of exponent on bra center to that on ket center 77 integer order_geo !incremental recorder over orders of geometric derivatives 78 integer num_up_geo !number of upper order geometric derivatives 79 integer num_low_geo !number of lower order geometric derivatives 80 integer start_up_geo !start address of upper order geometric derivatives 81 integer end_up_geo !end address of upper order geometric derivatives 82 integer start_low_geo !start address of lower order geometric derivatives 83 integer end_low_geo !end address of lower order geometric derivatives 84 integer dim_low_hket !dimension of lower order HGTOs on ket center 85 integer dim_up_hket !dimension of upper order HGTOs on ket center 86 integer size_low_hket !size of temporary integrals of lower order HGTOs on ket center 87 integer size_up_hket !size of temporary integrals of upper order HGTOs on ket center 88 integer low_hket_int !pointer to temporary integrals of lower order HGTOs on ket center 89 integer up_hket_int !pointer to temporary integrals of upper order HGTOs on ket center 90 integer dim_tmp !dimension of temporary integrals 91 real(REALK), allocatable :: tmp_ints(:,:) 92 !temporary integrals 93 integer offset_hket_geo !offset of geometric derivatives in returned integrals 94 integer ierr !error information 95#if defined(XTIME) 96 real(REALK) curr_time !current CPU time 97 ! sets current CPU time 98 call xtimer_set(curr_time) 99#endif 100 select case(orders_hgto_ket(2)) 101 ! returns s-shell HGTO 102 case(0) 103 hket_pints(1,:) = geo_pot_pints 104 ! maximum returned HGTOs are p-shell 105 case(1) 106 end_low_geo = dim_geo_pot !end address of lower order geometric derivatives 107 num_low_geo = (orders_geo_pot(2)+2) & !number of lower order geometric derivatives 108 * (orders_geo_pot(2)+3)/2 109 start_low_geo = end_low_geo-num_low_geo+1 !start address of lower order derivatives 110 ! allocates memory for temporary integrals 111 allocate(tmp_ints(3*(num_low_geo-(orders_geo_pot(2)+2)),1), stat=ierr) 112 if (ierr/=0) & 113 call error_stop("nucpot_hket", "failed to allocate tmp_ints", & 114 3*(num_low_geo-orders_geo_pot(2)-2)) 115 ! if returing zero order HGTO 116 zero_hket = orders_hgto_ket(1)==0 117 ! sets the offset of geometric derivatives in returned integrals 118 offset_hket_geo = dim_geo_hket 119 ! calculates the relative coordinates of center-of-charge w.r.t. ket center, 120 ! and the half of negative reciprocal of total exponent 121 half_neg_rp = 1.0_REALK/(exponent_bra+exponent_ket) 122 do ierr = 1, 3 123 cc_wrt_ket(ierr) = half_neg_rp*(exponent_bra*coord_bra(ierr) & 124 + exponent_ket*coord_ket(ierr))-coord_ket(ierr) 125 end do 126 half_neg_rp = -0.5_REALK*half_neg_rp 127 ! loops over orders of geometric derivatives 128 do order_geo = orders_geo_pot(2), orders_geo_pot(1), -1 129 ! updates the numbers of lower and upper order geometric derivatives 130 num_up_geo = num_low_geo 131 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 132 ! updates the start and end addresses of lower and upper order geometric derivatives 133 end_up_geo = end_low_geo 134 start_up_geo = start_low_geo 135 end_low_geo = start_low_geo-1 136 start_low_geo = end_low_geo-num_low_geo+1 137#if defined(DEBUG) 138 write(STDOUT,100) "GEO_HGTO/p/loop/upper/start/end:", & 139 order_geo+1, start_up_geo, end_up_geo 140 write(STDOUT,100) "GEO-HGTO/p/loop/lower/start/end:", & 141 order_geo, start_low_geo, end_low_geo 142#endif 143 ! sets the size of temporary integrals 144 size_up_hket = 3*num_low_geo 145 ! gets the p-shell HGTO and order \var(orders_geo_pot(2)) geometric derivatives integrals 146 call nucpot_hket_p(order_geo, cc_wrt_ket, half_neg_rp, & 147 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 148 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 149 3, tmp_ints(1:size_up_hket,1)) 150 ! updates the offset of geometric derivatives 151 offset_hket_geo = offset_hket_geo-num_low_geo 152 ! assigns the returned integrals 153 call nucpot_hket_assign(start_low_geo-1, dim_geo_pot, geo_pot_pints, & 154 3, num_low_geo, tmp_ints(1:size_up_hket,1), & 155 zero_hket, offset_hket_geo, dim_hgto_ket, & 156 dim_geo_hket, hket_pints) 157 end do 158 deallocate(tmp_ints) 159 ! maximum returned HGTOs are d-shell 160 case(2) 161 ! (1) \var(max_low_geo)-1 order geometric derivatives 162 max_low_geo = orders_geo_pot(2)+2 !maximum of lower order of geometric derivatives 163 end_up_geo = dim_geo_pot !end address of upper order geometric derivatives 164 num_up_geo = (max_low_geo+1)*(max_low_geo+2)/2 !number of upper order geometric derivatives 165 end_low_geo = end_up_geo-num_up_geo !end address of lower order geometric derivatives 166 start_up_geo = end_low_geo+1 !start address of upper order geometric derivatives 167 max_low_geo = max_low_geo-1 168 num_low_geo = num_up_geo-(max_low_geo+2) !number of lower order geometric derivatives 169 start_low_geo = end_low_geo-num_low_geo+1 !end address of \var(orders_geo_pot(2)) order derivatives 170 ! allocates memory for temporary integrals 171 allocate(tmp_ints(9*num_low_geo,2), stat=ierr) 172 if (ierr/=0) & 173 call error_stop("nucpot_hket", "failed to allocate tmp_ints", & 174 9*num_low_geo*2) 175#if defined(DEBUG) 176 write(STDOUT,100) "GEO-s-HGTO/d/upper/start/end:", & 177 max_low_geo+1, start_up_geo, end_up_geo 178 write(STDOUT,100) "GEO-s-HGTO/d/upper/start/end:", & 179 max_low_geo, start_low_geo, end_low_geo 180#endif 181 ! calculates the relative coordinates of center-of-charge w.r.t. ket center, 182 ! and the half of negative reciprocal of total exponent 183 half_neg_rp = 1.0_REALK/(exponent_bra+exponent_ket) 184 do ierr = 1, 3 185 cc_wrt_ket(ierr) = half_neg_rp*(exponent_bra*coord_bra(ierr) & 186 + exponent_ket*coord_ket(ierr))-coord_ket(ierr) 187 end do 188 half_neg_rp = -0.5_REALK*half_neg_rp 189 ! sets the dimension of lower order HGTOs and size of temporary integrals, p-shell 190 dim_low_hket = 3 191 size_low_hket = dim_low_hket*num_low_geo 192 ! gets the temporary p-shell HGTO integrals 193 call nucpot_hket_p(max_low_geo, cc_wrt_ket, half_neg_rp, & 194 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 195 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 196 dim_low_hket, tmp_ints(1:size_low_hket,1)) 197 ! (2) \var(max_low_geo)-2 order geometric derivatives 198 max_low_geo = max_low_geo-1 199 num_up_geo = num_low_geo 200 num_low_geo = num_low_geo-(max_low_geo+2) 201 end_up_geo = end_low_geo 202 start_up_geo = start_low_geo 203 end_low_geo = start_low_geo-1 204 start_low_geo = end_low_geo-num_low_geo+1 205#if defined(DEBUG) 206 write(STDOUT,100) "GEO-p-HGTO/d/upper/start/end:", & 207 max_low_geo+1, start_up_geo, end_up_geo 208 write(STDOUT,100) "GEO-p-HGTO/d/lower/start/end:", & 209 max_low_geo, start_low_geo, end_low_geo 210#endif 211 ! sets the dimension of upper order HGTOs and size of temporary integrals, p- and d-shell 212 dim_up_hket = 9 213 size_up_hket = dim_up_hket*num_low_geo 214 ! gets the temporary p-shell HGTO integrals 215 call nucpot_hket_p(max_low_geo, cc_wrt_ket, half_neg_rp, & 216 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 217 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 218 dim_up_hket, tmp_ints(1:size_up_hket,2)) 219 ! sets the ratio of exponent on bra center to that on ket center 220 bra_to_ket = exponent_bra/exponent_ket 221 ! recovers d-shell HGTOs 222 call nucpot_hket_d(max_low_geo, cc_wrt_ket, half_neg_rp, bra_to_ket, & 223 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 224 dim_low_hket, num_up_geo, tmp_ints(1:size_low_hket,1), & 225 dim_up_hket, tmp_ints(1:size_up_hket,2)) 226 ! if returing zero order HGTO 227 zero_hket = orders_hgto_ket(1)==0 228 ! sets the offset of geometric derivatives 229 offset_hket_geo = dim_geo_hket-num_low_geo 230 ! assigns the returned integrals 231 call nucpot_hket_assign(start_low_geo-1, dim_geo_pot, geo_pot_pints, & 232 dim_up_hket, num_low_geo, & 233 tmp_ints(1:size_up_hket,2), zero_hket, & 234 offset_hket_geo, dim_hgto_ket, dim_geo_hket, & 235 hket_pints) 236 ! initializes the pointers of HGTOs on ket center 237 low_hket_int = 1 238 up_hket_int = 2 239 ! updates the dimension of lower order HGTOs 240 dim_low_hket = dim_up_hket 241 ! loops over other returned orders of geometric derivatives 242 do order_geo = orders_geo_pot(2)-1, orders_geo_pot(1), -1 243 ! updates the numbers of lower and upper order geometric derivatives 244 num_up_geo = num_low_geo 245 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 246 ! updates the start and end addresses of lower and upper order geometric derivatives 247 end_up_geo = end_low_geo 248 start_up_geo = start_low_geo 249 end_low_geo = start_low_geo-1 250 start_low_geo = end_low_geo-num_low_geo+1 251#if defined(DEBUG) 252 write(STDOUT,100) "GEO-HGTO/d/loop/upper/start/end:", & 253 order_geo+1, start_up_geo, end_up_geo 254 write(STDOUT,100) "GEO-HGTO/d/loop/lower/start/end:", & 255 order_geo, start_low_geo, end_low_geo 256#endif 257 ! updates the sizes of temporary integrals 258 size_low_hket = size_up_hket 259 size_up_hket = dim_up_hket*num_low_geo 260 ! switches the pointers 261 low_hket_int = 3-low_hket_int 262 up_hket_int = 3-up_hket_int 263 ! gets the temporary p-shell HGTO integrals 264 call nucpot_hket_p(order_geo, cc_wrt_ket, half_neg_rp, & 265 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 266 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 267 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 268 ! gets the temporary d-shell HGTO integrals 269 call nucpot_hket_d(order_geo, cc_wrt_ket, half_neg_rp, & 270 bra_to_ket, num_low_geo, & 271 geo_pot_pints(start_low_geo:end_low_geo), & 272 dim_low_hket, num_up_geo, & 273 tmp_ints(1:size_low_hket,low_hket_int), & 274 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 275 ! updates the offset of geometric derivatives 276 offset_hket_geo = offset_hket_geo-num_low_geo 277 ! assigns the returned integrals 278 call nucpot_hket_assign(start_low_geo-1, dim_geo_pot, geo_pot_pints, & 279 dim_up_hket, num_low_geo, & 280 tmp_ints(1:size_up_hket,up_hket_int), & 281 zero_hket, offset_hket_geo, dim_hgto_ket, & 282 dim_geo_hket, hket_pints) 283 end do 284 deallocate(tmp_ints) 285 ! maximum returned HGTOs are other shells, at least f-shell 286 case default 287 ! (1) \var(max_low_geo)-1 order geometric derivatives 288 max_low_geo = orders_geo_pot(2) & !maximum of lower order of geometric derivatives 289 + orders_hgto_ket(2) 290 ! allocates memory for temporary integrals 291 call dim_nucpot_hket(max_low_geo, orders_geo_pot(2), dim_tmp) 292 allocate(tmp_ints(dim_tmp,2), stat=ierr) 293 if (ierr/=0) & 294 call error_stop("nucpot_hket", "failed to allocate tmp_ints", & 295 dim_tmp*2) 296 end_up_geo = dim_geo_pot !end address of upper order geometric derivatives 297 num_up_geo = (max_low_geo+1)*(max_low_geo+2)/2 !number of upper order geometric derivatives 298 end_low_geo = end_up_geo-num_up_geo !end address of lower order geometric derivatives 299 start_up_geo = end_low_geo+1 !start address of upper order geometric derivatives 300 max_low_geo = max_low_geo-1 301 num_low_geo = num_up_geo-(max_low_geo+2) !number of lower order geometric derivatives 302 start_low_geo = end_low_geo-num_low_geo+1 !end address of \var(orders_geo_pot(2)) order derivatives 303#if defined(DEBUG) 304 write(STDOUT,100) "GEO-s-HGTO/upper/start/end:", & 305 max_low_geo+1, start_up_geo, end_up_geo 306 write(STDOUT,100) "GEO-s-HGTO/upper/start/end:", & 307 max_low_geo, start_low_geo, end_low_geo 308#endif 309 ! calculates the relative coordinates of center-of-charge w.r.t. ket center, 310 ! and the half of negative reciprocal of total exponent 311 half_neg_rp = 1.0_REALK/(exponent_bra+exponent_ket) 312 do ierr = 1, 3 313 cc_wrt_ket(ierr) = half_neg_rp*(exponent_bra*coord_bra(ierr) & 314 + exponent_ket*coord_ket(ierr))-coord_ket(ierr) 315 end do 316 half_neg_rp = -0.5_REALK*half_neg_rp 317 ! sets the dimension of lower order HGTOs and size of temporary integrals, p-shell 318 dim_low_hket = 3 319 size_low_hket = dim_low_hket*num_low_geo 320 ! gets the temporary p-shell HGTO integrals 321 call nucpot_hket_p(max_low_geo, cc_wrt_ket, half_neg_rp, & 322 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 323 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 324 dim_low_hket, tmp_ints(1:size_low_hket,1)) 325 ! (2) \var(max_low_geo)-2 order geometric derivatives 326 max_low_geo = max_low_geo-1 327 num_up_geo = num_low_geo 328 num_low_geo = num_low_geo-(max_low_geo+2) 329 end_up_geo = end_low_geo 330 start_up_geo = start_low_geo 331 end_low_geo = start_low_geo-1 332 start_low_geo = end_low_geo-num_low_geo+1 333#if defined(DEBUG) 334 write(STDOUT,100) "GEO-p-HGTO/upper/start/end:", & 335 max_low_geo+1, start_up_geo, end_up_geo 336 write(STDOUT,100) "GEO-p-HGTO/lower/start/end:", & 337 max_low_geo, start_low_geo, end_low_geo 338#endif 339 ! sets the dimension of upper order HGTOs and size of temporary integrals, p- and d-shell 340 dim_up_hket = 9 341 size_up_hket = dim_up_hket*num_low_geo 342 ! gets the temporary p-shell HGTO integrals 343 call nucpot_hket_p(max_low_geo, cc_wrt_ket, half_neg_rp, & 344 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 345 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 346 dim_up_hket, tmp_ints(1:size_up_hket,2)) 347 ! sets the ratio of exponent on bra center to that on ket center 348 bra_to_ket = exponent_bra/exponent_ket 349 ! recovers d-shell HGTOs 350 call nucpot_hket_d(max_low_geo, cc_wrt_ket, half_neg_rp, bra_to_ket, & 351 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 352 dim_low_hket, num_up_geo, tmp_ints(1:size_low_hket,1), & 353 dim_up_hket, tmp_ints(1:size_up_hket,2)) 354 ! initializes the maximum of current order of HGTOs on ket center 355 max_cur_hket = 1 356 ! initializes the pointers of HGTOs on ket center 357 low_hket_int = 1 358 up_hket_int = 2 359 ! (3) loops over the orders of geometric derivatives till the maximum returned 360 ! order, the maximum of current order of HGTOs \var(max_cur_hket) needs to 361 ! update each iteration 362 do order_geo = max_low_geo-1, orders_geo_pot(2), -1 363 ! updates the numbers of lower and upper order geometric derivatives 364 num_up_geo = num_low_geo 365 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 366 ! updates the start and end addresses of lower and upper order geometric derivatives 367 end_up_geo = end_low_geo 368 start_up_geo = start_low_geo 369 end_low_geo = start_low_geo-1 370 start_low_geo = end_low_geo-num_low_geo+1 371#if defined(DEBUG) 372 write(STDOUT,100) "GEO-HGTO/loop/1/upper/start/end:", & 373 order_geo+1, start_up_geo, end_up_geo 374 write(STDOUT,100) "GEO-HGTO/loop/1/lower/start/end:", & 375 order_geo, start_low_geo, end_low_geo 376#endif 377 ! updates the maximum of current order of HGTOs 378 max_cur_hket = max_cur_hket+1 379 ! updates the dimensions of HGTOs on ket center 380 dim_low_hket = dim_up_hket 381 dim_up_hket = dim_low_hket+(max_cur_hket+2)*(max_cur_hket+3)/2 382 ! updates the sizes of temporary integrals 383 size_low_hket = size_up_hket 384 size_up_hket = dim_up_hket*num_low_geo 385 ! switches the pointers 386 low_hket_int = 3-low_hket_int 387 up_hket_int = 3-up_hket_int 388 ! gets the temporary p-shell HGTO integrals 389 call nucpot_hket_p(order_geo, cc_wrt_ket, half_neg_rp, & 390 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 391 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 392 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 393 ! gets the temporary d-shell HGTO integrals 394 call nucpot_hket_d(order_geo, cc_wrt_ket, half_neg_rp, & 395 bra_to_ket, num_low_geo, & 396 geo_pot_pints(start_low_geo:end_low_geo), & 397 dim_low_hket, num_up_geo, & 398 tmp_ints(1:size_low_hket,low_hket_int), & 399 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 400 ! gets the temporary integrals for other HGTOs (f, g, ...) 401 ! and \var(order_geo) order geometric derivatives 402 call sub_nucpot_hket(order_geo, max_cur_hket, cc_wrt_ket, half_neg_rp, & 403 bra_to_ket, dim_low_hket, num_up_geo, & 404 tmp_ints(1:size_low_hket,low_hket_int), dim_up_hket, & 405 num_low_geo, tmp_ints(1:size_up_hket,up_hket_int)) 406 end do 407 ! if returing zero order HGTO 408 zero_hket = orders_hgto_ket(1)==0 409 ! sets the offset of geometric derivatives 410 offset_hket_geo = dim_geo_hket-num_low_geo 411 ! assigns the returned integrals 412 call nucpot_hket_assign(start_low_geo-1, dim_geo_pot, geo_pot_pints, & 413 dim_up_hket, num_low_geo, & 414 tmp_ints(1:size_up_hket,up_hket_int), & 415 zero_hket, offset_hket_geo, dim_hgto_ket, & 416 dim_geo_hket, hket_pints) 417 ! updates the dimension of lower order HGTOs 418 dim_low_hket = dim_up_hket 419 ! loops over other returned orders of geometric derivatives, the maximum 420 ! of current order of HGTOs \var(max_cur_hket) does not need to update 421 do order_geo = orders_geo_pot(2)-1, orders_geo_pot(1), -1 422 ! updates the numbers of lower and upper order geometric derivatives 423 num_up_geo = num_low_geo 424 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 425 ! updates the start and end addresses of lower and upper order geometric derivatives 426 end_up_geo = end_low_geo 427 start_up_geo = start_low_geo 428 end_low_geo = start_low_geo-1 429 start_low_geo = end_low_geo-num_low_geo+1 430#if defined(DEBUG) 431 write(STDOUT,100) "GEO-HGTO/loop/2/upper/start/end:", & 432 order_geo+1, start_up_geo, end_up_geo 433 write(STDOUT,100) "GEO-HGTO/loop/2/lower/start/end:", & 434 order_geo, start_low_geo, end_low_geo 435#endif 436 ! updates the sizes of temporary integrals 437 size_low_hket = size_up_hket 438 size_up_hket = dim_up_hket*num_low_geo 439 ! switches the pointers 440 low_hket_int = 3-low_hket_int 441 up_hket_int = 3-up_hket_int 442 ! gets the temporary p-shell HGTO integrals 443 call nucpot_hket_p(order_geo, cc_wrt_ket, half_neg_rp, & 444 num_up_geo, geo_pot_pints(start_up_geo:end_up_geo), & 445 num_low_geo, geo_pot_pints(start_low_geo:end_low_geo), & 446 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 447 ! gets the temporary d-shell HGTO integrals 448 call nucpot_hket_d(order_geo, cc_wrt_ket, half_neg_rp, & 449 bra_to_ket, num_low_geo, & 450 geo_pot_pints(start_low_geo:end_low_geo), & 451 dim_low_hket, num_up_geo, & 452 tmp_ints(1:size_low_hket,low_hket_int), & 453 dim_up_hket, tmp_ints(1:size_up_hket,up_hket_int)) 454 ! gets the temporary integrals for other HGTOs (f, g, ...) 455 ! and \var(order_geo) order geometric derivatives 456 call sub_nucpot_hket(order_geo, max_cur_hket, cc_wrt_ket, half_neg_rp, & 457 bra_to_ket, dim_low_hket, num_up_geo, & 458 tmp_ints(1:size_low_hket,low_hket_int), dim_up_hket, & 459 num_low_geo, tmp_ints(1:size_up_hket,up_hket_int)) 460 ! updates the offset of geometric derivatives 461 offset_hket_geo = offset_hket_geo-num_low_geo 462 ! assigns the returned integrals 463 call nucpot_hket_assign(start_low_geo-1, dim_geo_pot, geo_pot_pints, & 464 dim_up_hket, num_low_geo, & 465 tmp_ints(1:size_up_hket,up_hket_int), & 466 zero_hket, offset_hket_geo, dim_hgto_ket, & 467 dim_geo_hket, hket_pints) 468 end do 469 deallocate(tmp_ints) 470 end select 471#if defined(XTIME) 472 ! prints the CPU elapsed time 473 call xtimer_view(curr_time, "nucpot_hket", STDOUT) 474#endif 475 return 476#if defined(DEBUG) 477100 format("nucpot_hket>> ",A,I6,2I8) 478#endif 479 end subroutine nucpot_hket 480 481 !> \brief gets the maximum dimension of temporary integrals used in recurrence relations 482 !> \author Bin Gao 483 !> \date 2012-03-05 484 !> \param max_order_geo is the maximum order of geometric derivatives 485 !> \param min_order_geo is the minimum order of geometric derivatives 486 !> \return dim_ints is the maximum dimension of temporary integrals 487 subroutine dim_nucpot_hket(max_order_geo, min_order_geo, dim_ints) 488 use xkind 489 implicit none 490 integer, intent(in) :: max_order_geo 491 integer, intent(in) :: min_order_geo 492 integer, intent(out) :: dim_ints 493!f2py intent(in) :: max_order_geo 494!f2py intent(in) :: min_order_geo 495!f2py intent(out) :: dim_ints 496 integer num_geo_pot !number of upper order geometric derivatives 497 integer max_hgto_ket !maximum order of HGTOs on ket center 498 integer dim_hgto_ket !dimension of HGTOs on ket center 499 integer order_geo !incremental recorder over orders of geometric derivatives 500 integer dim_tmp !temporary result of dimension 501#if defined(XTIME) 502 real(REALK) curr_time !current CPU time 503 ! sets current CPU time 504 call xtimer_set(curr_time) 505#endif 506 ! initializes the return value 507 dim_ints = 0 508 ! initializes the number of geometric derivatives 509 num_geo_pot = (max_order_geo+1)*(max_order_geo+2)/2 510 ! initializes the maximum order and dimension of HGTOs on ket center 511 max_hgto_ket = 0 512 dim_hgto_ket = 0 513 ! loops over the orders of geometric derivatives till the maximum order 514 ! of returned geometric derivatives, the maximum order of HGTOs needs to 515 ! update each iteration 516 do order_geo = max_order_geo-1, min_order_geo, -1 517 ! updates the numbers of geometric derivatives 518 num_geo_pot = num_geo_pot-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 519 ! updates the maximum order of HGTOs 520 max_hgto_ket = max_hgto_ket+1 521 ! updates the dimension of HGTOs on ket center 522 dim_hgto_ket = dim_hgto_ket+(max_hgto_ket+1)*(max_hgto_ket+2)/2 523 ! updates the maximum dimension 524 dim_tmp = dim_hgto_ket*num_geo_pot 525 if (dim_tmp>dim_ints) dim_ints = dim_tmp 526 end do 527#if defined(XTIME) 528 ! prints the CPU elapsed time 529 call xtimer_view(curr_time, "dim_nucpot_hket", STDOUT) 530#endif 531 return 532 end subroutine dim_nucpot_hket 533 534 !> \brief recovers p-shell HGTOs 535 !> \author Bin Gao 536 !> \date 2011-10-18 537 !> \param cur_order_geo is current order of geometric derivatives 538 !> \param cc_wrt_ket contains the relative coordinates of center-of-charge w.r.t. ket center 539 !> \param half_neg_rp is the half of the negative reciprocal of total exponent \f$p_{ij}\f$ 540 !> \param num_up_geo is the number of upper order geometric derivatives 541 !> \param up_geo_pints contains the integrals with s-shell HGTO and upper order 542 !> geometric derivatives 543 !> \param num_low_geo is the number lower order geometric derivatives 544 !> \param low_geo_pints contains the integrals with s-shell HGTO and lower order 545 !> geometric derivatives 546 !> \param dim_up_hket is the dimension of upper order HGTOs 547 !> \return up_hket_pints contains the integrals of p-shell HGTOs and lower order 548 !> geometric derivatives 549 subroutine nucpot_hket_p(cur_order_geo, cc_wrt_ket, half_neg_rp, & 550 num_up_geo, up_geo_pints, num_low_geo, & 551 low_geo_pints, dim_up_hket, up_hket_pints) 552 use xkind 553 implicit none 554 integer, intent(in) :: cur_order_geo 555 real(REALK), intent(in) :: cc_wrt_ket(3) 556 real(REALK), intent(in) :: half_neg_rp 557 integer, intent(in) :: num_up_geo 558 real(REALK), intent(in) :: up_geo_pints(num_up_geo) 559 integer, intent(in) :: num_low_geo 560 real(REALK), intent(in) :: low_geo_pints(num_low_geo) 561 integer, intent(in) :: dim_up_hket 562 real(REALK), intent(inout) :: up_hket_pints(dim_up_hket,num_low_geo) 563!f2py intent(in) :: cur_order_geo 564!f2py intent(in) :: cc_wrt_ket 565!f2py intent(in) :: half_neg_rp 566!f2py intent(hide) :: num_up_geo 567!f2py intent(in) :: up_geo_pints 568!f2py intent(hide) :: num_low_geo 569!f2py intent(in) :: low_geo_pints 570!f2py intent(hide) :: dim_up_hket 571!f2py intent(inout) :: up_hket_pints 572!f2py depend(num_low_geo) :: up_hket_pints 573 integer addr_up_geo !address of upper order geometric derivatives 574 integer addr_cur_geo !address of current order geometric derivatives 575 integer igeo, jgeo !incremental recorders over geometric derivatives 576#if defined(XTIME) 577 real(REALK) curr_time !current CPU time 578 ! sets current CPU time 579 call xtimer_set(curr_time) 580#endif 581 addr_up_geo = 0 582 addr_cur_geo = 0 583 do igeo = cur_order_geo, 0, -1 584 do jgeo = 0, igeo 585 addr_up_geo = addr_up_geo+1 586 addr_cur_geo = addr_cur_geo+1 587 ! px 588 up_hket_pints(1,addr_cur_geo) & 589 = cc_wrt_ket(1)*low_geo_pints(addr_cur_geo) & 590 + half_neg_rp*up_geo_pints(addr_up_geo) 591 ! py 592 up_hket_pints(2,addr_cur_geo) & 593 = cc_wrt_ket(2)*low_geo_pints(addr_cur_geo) & 594 + half_neg_rp*up_geo_pints(addr_up_geo+1) 595 ! pz 596 up_hket_pints(3,addr_cur_geo) & 597 = cc_wrt_ket(3)*low_geo_pints(addr_cur_geo) & 598 + half_neg_rp*up_geo_pints(addr_up_geo+igeo+2) 599 end do 600 addr_up_geo = addr_up_geo+1 601 end do 602#if defined(XTIME) 603 ! prints the CPU elapsed time 604 call xtimer_view(curr_time, "nucpot_hket_p", STDOUT) 605#endif 606 return 607 end subroutine nucpot_hket_p 608 609 !> \brief recovers d-shell HGTOs 610 !> \author Bin Gao 611 !> \date 2011-10-18 612 !> \param cur_order_geo is current order of geometric derivatives 613 !> \param cc_wrt_ket contains the relative coordinates of center-of-charge w.r.t. ket center 614 !> \param half_neg_rp is the half of the negative reciprocal of total exponent \f$p_{ij}\f$ 615 !> \param bra_to_ket is the ratio of exponent on bra center to that on ket center 616 !> \param low_geo_pints contains the integrals with s-shell HGTO and lower order 617 !> geometric derivatives 618 !> \param dim_low_hket is the dimension of lower order HGTOs 619 !> \param num_up_geo is the number of upper order geometric derivatives 620 !> \param low_hket_pints contains the integrals with p-shell HGTOs and upper order 621 !> geometric derivatives 622 !> \param dim_up_hket is the dimension of upper order HGTOs 623 !> \return up_hket_pints contains the integrals of d-shell HGTOs and lower order 624 !> geometric derivatives 625 subroutine nucpot_hket_d(cur_order_geo, cc_wrt_ket, half_neg_rp, & 626 bra_to_ket, num_low_geo, low_geo_pints, & 627 dim_low_hket, num_up_geo, low_hket_pints, & 628 dim_up_hket, up_hket_pints) 629 use xkind 630 implicit none 631 integer, intent(in) :: cur_order_geo 632 real(REALK), intent(in) :: cc_wrt_ket(3) 633 real(REALK), intent(in) :: half_neg_rp 634 real(REALK), intent(in) :: bra_to_ket 635 integer, intent(in) :: num_low_geo 636 real(REALK), intent(in) :: low_geo_pints(num_low_geo) 637 integer, intent(in) :: dim_low_hket 638 integer, intent(in) :: num_up_geo 639 real(REALK), intent(in) :: low_hket_pints(dim_low_hket,num_up_geo) 640 integer, intent(in) :: dim_up_hket 641 real(REALK), intent(inout) :: up_hket_pints(dim_up_hket,num_low_geo) 642!f2py intent(in) :: cur_order_geo 643!f2py intent(in) :: cc_wrt_ket 644!f2py intent(in) :: half_neg_rp 645!f2py intent(in) :: bra_to_ket 646!f2py intent(hide) :: num_low_geo 647!f2py intent(in) :: low_geo_pints 648!f2py intent(hide) :: dim_low_hket 649!f2py intent(hide) :: num_up_geo 650!f2py intent(in) :: low_hket_pints 651!f2py intent(hide) :: dim_up_hket 652!f2py intent(inout) :: up_hket_pints 653!f2py depend(num_low_geo) :: up_hket_pints 654 integer addr_up_geo !addresses of upper order geometric derivatives 655 integer addr_up_geo_y 656 integer addr_up_geo_z 657 integer addr_cur_geo !address of current order geometric derivatives 658 integer igeo, jgeo !incremental recorders over geometric derivatives 659#if defined(XTIME) 660 real(REALK) curr_time !current CPU time 661 ! sets current CPU time 662 call xtimer_set(curr_time) 663#endif 664 addr_up_geo = 0 665 addr_cur_geo = 0 666 do igeo = cur_order_geo, 0, -1 667 do jgeo = 0, igeo 668 addr_up_geo = addr_up_geo+1 669 addr_cur_geo = addr_cur_geo+1 670 ! dxx 671 up_hket_pints(4,addr_cur_geo) & 672 = cc_wrt_ket(1)*up_hket_pints(1,addr_cur_geo) & 673 + half_neg_rp*(bra_to_ket*low_geo_pints(addr_cur_geo) & 674 + low_hket_pints(1,addr_up_geo)) 675 addr_up_geo_y = addr_up_geo+1 676 ! dxy 677 up_hket_pints(5,addr_cur_geo) & 678 = cc_wrt_ket(2)*up_hket_pints(1,addr_cur_geo) & 679 + half_neg_rp*low_hket_pints(1,addr_up_geo_y) 680 ! dyy 681 up_hket_pints(6,addr_cur_geo) & 682 = cc_wrt_ket(2)*up_hket_pints(2,addr_cur_geo) & 683 + half_neg_rp*(bra_to_ket*low_geo_pints(addr_cur_geo) & 684 + low_hket_pints(2,addr_up_geo_y)) 685 addr_up_geo_z = addr_up_geo+igeo+2 686 ! dxz 687 up_hket_pints(7,addr_cur_geo) & 688 = cc_wrt_ket(3)*up_hket_pints(1,addr_cur_geo) & 689 + half_neg_rp*low_hket_pints(1,addr_up_geo_z) 690 ! dyz 691 up_hket_pints(8,addr_cur_geo) & 692 = cc_wrt_ket(3)*up_hket_pints(2,addr_cur_geo) & 693 + half_neg_rp*low_hket_pints(2,addr_up_geo_z) 694 ! dzz 695 up_hket_pints(9,addr_cur_geo) & 696 = cc_wrt_ket(3)*up_hket_pints(3,addr_cur_geo) & 697 + half_neg_rp*(bra_to_ket*low_geo_pints(addr_cur_geo) & 698 + low_hket_pints(3,addr_up_geo_z)) 699 end do 700 addr_up_geo = addr_up_geo+1 701 end do 702#if defined(XTIME) 703 ! prints the CPU elapsed time 704 call xtimer_view(curr_time, "nucpot_hket_d", STDOUT) 705#endif 706 return 707 end subroutine nucpot_hket_d 708 709 !> \brief sub-recurrence relations by recovering upper order HGTOs on ket center 710 !> \author Bin Gao 711 !> \date 2011-10-18 712 !> \param cur_order_geo is current order of geometric derivatives 713 !> \param max_cur_hket is maximum of current order of HGTOs on ket center 714 !> \param cc_wrt_ket contains the relative coordinates of center-of-charge w.r.t. ket center 715 !> \param half_neg_rp is the half of the negative reciprocal of total exponent \f$p_{ij}\f$ 716 !> \param bra_to_ket is the ratio of exponent on bra center to that on ket center 717 !> \param dim_low_hket is the dimension of lower order HGTOs 718 !> \param num_up_geo is the number of upper order geometric derivatives 719 !> \param low_hket_pints contains the integrals with lower order HGTOs and upper order 720 !> geometric derivatives 721 !> \param dim_up_hket is the dimension of upper order HGTOs 722 !> \param num_low_geo is the number lower order geometric derivatives 723 !> \return up_hket_pints contains the integrals of upper order HGTOs and lower order 724 !> geometric derivatives 725 subroutine sub_nucpot_hket(cur_order_geo, max_cur_hket, cc_wrt_ket, & 726 half_neg_rp, bra_to_ket, dim_low_hket, & 727 num_up_geo, low_hket_pints, dim_up_hket, & 728 num_low_geo, up_hket_pints) 729 use xkind 730 implicit none 731 integer, intent(in) :: cur_order_geo 732 integer, intent(in) :: max_cur_hket 733 real(REALK), intent(in) :: cc_wrt_ket(3) 734 real(REALK), intent(in) :: half_neg_rp 735 real(REALK), intent(in) :: bra_to_ket 736 integer, intent(in) :: dim_low_hket 737 integer, intent(in) :: num_up_geo 738 real(REALK), intent(in) :: low_hket_pints(dim_low_hket,num_up_geo) 739 integer, intent(in) :: dim_up_hket 740 integer, intent(in) :: num_low_geo 741 real(REALK), intent(inout) :: up_hket_pints(dim_up_hket,num_low_geo) 742!f2py intent(in) :: cur_order_geo 743!f2py intent(in) :: max_cur_hket 744!f2py intent(in) :: cc_wrt_ket 745!f2py intent(in) :: half_neg_rp 746!f2py intent(in) :: bra_to_ket 747!f2py intent(hide) :: dim_low_hket 748!f2py intent(hide) :: num_up_geo 749!f2py intent(in) :: low_hket_pints 750!f2py intent(hide) :: dim_up_hket 751!f2py intent(hide) :: num_low_geo 752!f2py intent(inout) :: up_hket_pints 753 integer addr_up_geo !address of upper order geometric derivatives 754 integer addr_cur_geo !address of current order geometric derivatives 755 integer igeo, jgeo !incremental recorders over geometric derivatives 756 integer order_hket !order of HGTOs on ket center 757 integer addr_up_hket !address of upper order HGTOs 758 integer addr_cur_hket !address of current order HGTOs 759 integer addr_low_hket !address of lower order HGTOs 760 integer iket, jket !incremental recorder over HGTOs 761#if defined(XTIME) 762 real(REALK) curr_time !current CPU time 763 ! sets current CPU time 764 call xtimer_set(curr_time) 765#endif 766 addr_up_geo = 0 767 addr_cur_geo = 0 768 ! loops over xyz components of geometric derivatives 769 do igeo = cur_order_geo, 0, -1 770 do jgeo = 0, igeo 771 addr_up_geo = addr_up_geo+1 772 addr_cur_geo = addr_cur_geo+1 773 addr_up_hket = 9 !base address of the f-shell HGTOs 774 addr_cur_hket = 3 !base address of the d-shell HGTOs 775 addr_low_hket = 0 !base address of the p-shell HGTOs 776 ! loops over current order of HGTOs, starting from d-shell 777 do order_hket = 2, max_cur_hket 778 addr_up_hket = addr_up_hket+1 779 addr_cur_hket = addr_cur_hket+1 780#if defined(DEBUG) 781 write(STDOUT,100) "orders:", order_hket, cur_order_geo 782 write(STDOUT,100) "x-direction" 783 write(STDOUT,110) addr_up_hket, addr_cur_geo 784 write(STDOUT,110) addr_cur_hket, addr_cur_geo 785 write(STDOUT,110) addr_low_hket+1, addr_cur_geo 786 write(STDOUT,110) addr_cur_hket, addr_up_geo 787 write(STDOUT,100) "------------------------------------" 788#endif 789 ! x-direction 790 up_hket_pints(addr_up_hket,addr_cur_geo) & 791 = cc_wrt_ket(1)*up_hket_pints(addr_cur_hket,addr_cur_geo) & 792 + half_neg_rp*(real(order_hket,REALK)*bra_to_ket & 793 * up_hket_pints(addr_low_hket+1,addr_cur_geo) & 794 + low_hket_pints(addr_cur_hket,addr_up_geo)) 795 ! y-direction 796 addr_up_hket = addr_up_hket+1 797#if defined(DEBUG) 798 write(STDOUT,100) "y-direction" 799 write(STDOUT,110) addr_up_hket, addr_cur_geo 800 write(STDOUT,110) addr_cur_hket, addr_cur_geo 801 write(STDOUT,110) addr_cur_hket, addr_up_geo+1 802 write(STDOUT,100) "------------------------------------" 803#endif 804 up_hket_pints(addr_up_hket,addr_cur_geo) & 805 = cc_wrt_ket(2)*up_hket_pints(addr_cur_hket,addr_cur_geo) & 806 + half_neg_rp*low_hket_pints(addr_cur_hket,addr_up_geo+1) 807 do iket = 1, order_hket 808 addr_up_hket = addr_up_hket+1 809#if defined(DEBUG) 810 write(STDOUT,110) addr_up_hket, addr_cur_geo 811 write(STDOUT,110) addr_cur_hket+iket, addr_cur_geo 812 write(STDOUT,110) addr_low_hket+iket, addr_cur_geo 813 write(STDOUT,110) addr_cur_hket+iket, addr_up_geo+1 814 write(STDOUT,100) "------------------------------------" 815#endif 816 up_hket_pints(addr_up_hket,addr_cur_geo) & 817 = cc_wrt_ket(2)*up_hket_pints(addr_cur_hket+iket,addr_cur_geo) & 818 + half_neg_rp*(real(iket,REALK)*bra_to_ket & 819 * up_hket_pints(addr_low_hket+iket,addr_cur_geo) & 820 + low_hket_pints(addr_cur_hket+iket,addr_up_geo+1)) 821 end do 822 ! z-direction 823#if defined(DEBUG) 824 write(STDOUT,100) "z-direction" 825#endif 826 addr_cur_hket = addr_cur_hket-1 827 do jket = 0, order_hket 828 addr_up_hket = addr_up_hket+1 829 addr_cur_hket = addr_cur_hket+1 830#if defined(DEBUG) 831 write(STDOUT,110) addr_up_hket, addr_cur_geo 832 write(STDOUT,110) addr_cur_hket, addr_cur_geo 833 write(STDOUT,110) addr_cur_hket, addr_up_geo+igeo+2 834 write(STDOUT,100) "------------------------------------" 835#endif 836 up_hket_pints(addr_up_hket,addr_cur_geo) & 837 = cc_wrt_ket(3)*up_hket_pints(addr_cur_hket,addr_cur_geo) & 838 + half_neg_rp*low_hket_pints(addr_cur_hket,addr_up_geo+igeo+2) 839 end do 840 do iket = 1, order_hket 841 do jket = 0, order_hket-iket 842 addr_up_hket = addr_up_hket+1 843 addr_cur_hket = addr_cur_hket+1 844 addr_low_hket = addr_low_hket+1 845#if defined(DEBUG) 846 write(STDOUT,110) addr_up_hket, addr_cur_geo 847 write(STDOUT,110) addr_cur_hket, addr_cur_geo 848 write(STDOUT,110) addr_low_hket, addr_cur_geo 849 write(STDOUT,110) addr_cur_hket, addr_up_geo+igeo+2 850 write(STDOUT,100) "------------------------------------" 851#endif 852 up_hket_pints(addr_up_hket,addr_cur_geo) & 853 = cc_wrt_ket(3)*up_hket_pints(addr_cur_hket,addr_cur_geo) & 854 + half_neg_rp*(real(iket,REALK)*bra_to_ket & 855 * up_hket_pints(addr_low_hket,addr_cur_geo) & 856 + low_hket_pints(addr_cur_hket,addr_up_geo+igeo+2)) 857 end do 858 end do 859 end do 860 end do 861 addr_up_geo = addr_up_geo+1 862 end do 863#if defined(XTIME) 864 ! prints the CPU elapsed time 865 call xtimer_view(curr_time, "sub_nucpot_hket", STDOUT) 866#endif 867 return 868#if defined(DEBUG) 869100 format("sub_nucpot_hket>> ",A,2I6) 870110 format("sub_nucpot_hket>> ","HGTO",I8,4X,"GEO",I8) 871#endif 872 end subroutine sub_nucpot_hket 873 874 !> \brief assigns the integrals with required HGTOs on ket center 875 !> \author Bin Gao 876 !> \date 2012-03-04 877 !> \param offset_geo_pot is the offset of geometric derivatives on nuclear potential origin 878 !> \param dim_geo_pot is the dimension of geometric derivatives 879 !> \param geo_pot_pints contains the nuclear attraction integrals with zeroth order 880 !> HGTOs on ket center 881 !> \param dim_up_hket is the dimension of upper order HGTOs in temporary integrals 882 !> \param num_low_geo is the number of lower order geometric derivatives in temporary integrals 883 !> \param recur_pints contains the temporary integrals from recurrence relations 884 !> \param zero_hket indicates if zeroth order HGTO returned 885 !> \param offset_hket_geo is the offset of geometric derivatives on nuclear potential origin 886 !> in returned integrals 887 !> \param dim_hgto_ket is the dimension of HGTOs of ket center afterwards 888 !> \param dim_geo_hket is the dimension of geometric derivatives afterwards 889 !> \return hket_pints contains the integrals with required HGTOs on ket center 890 subroutine nucpot_hket_assign(offset_geo_pot, dim_geo_pot, geo_pot_pints, & 891 dim_up_hket, num_low_geo, recur_pints, & 892 zero_hket, offset_hket_geo, dim_hgto_ket, & 893 dim_geo_hket, hket_pints) 894 use xkind 895 implicit none 896 integer, intent(in) :: offset_geo_pot 897 integer, intent(in) :: dim_geo_pot 898 real(REALK), intent(in) :: geo_pot_pints(dim_geo_pot) 899 integer, intent(in) :: dim_up_hket 900 integer, intent(in) :: num_low_geo 901 real(REALK), intent(in) :: recur_pints(dim_up_hket,num_low_geo) 902 logical, intent(in) :: zero_hket 903 integer, intent(in) :: offset_hket_geo 904 integer, intent(in) :: dim_hgto_ket 905 integer, intent(in) :: dim_geo_hket 906 real(REALK), intent(inout) :: hket_pints(dim_hgto_ket,dim_geo_hket) 907!f2py intent(in) :: offset_geo_pot 908!f2py intent(hide) :: dim_geo_pot 909!f2py intent(in) :: geo_pot_pints 910!f2py intent(hide) :: dim_up_hket 911!f2py intent(hide) :: num_low_geo 912!f2py intent(in) :: recur_pints 913!f2py intent(in) :: zero_hket 914!f2py intent(in) :: offset_hket_geo 915!f2py intent(hide) :: dim_hgto_ket 916!f2py intent(hide) :: dim_geo_hket 917!f2py intent(inout) :: hket_pints 918 integer start_up_hket !start address of upper order HGTOs in temporary integrals 919 integer addr_geo_pot !address of geometric derivatives on nuclear potential origin 920 integer addr_hket_geo !address of geometric derivatives on nuclear potential origin in returned integrals 921 integer igeo !incremental recorder over geometric derivatives 922#if defined(XTIME) 923 real(REALK) curr_time !current CPU time 924 ! sets current CPU time 925 call xtimer_set(curr_time) 926#endif 927 ! s-shell HGTO returned 928 if (zero_hket) then 929 start_up_hket = dim_up_hket-dim_hgto_ket+2 930 do igeo = 1, num_low_geo 931 addr_geo_pot = offset_geo_pot+igeo 932 addr_hket_geo = offset_hket_geo+igeo 933 hket_pints(1,addr_hket_geo) = geo_pot_pints(addr_geo_pot) 934 hket_pints(2:dim_hgto_ket,addr_hket_geo) & 935 = recur_pints(start_up_hket:dim_up_hket,igeo) 936 end do 937 else 938 start_up_hket = dim_up_hket-dim_hgto_ket+1 939 do igeo = 1, num_low_geo 940 addr_hket_geo = offset_hket_geo+igeo 941 hket_pints(:,addr_hket_geo) & 942 = recur_pints(start_up_hket:dim_up_hket,igeo) 943 end do 944 end if 945#if defined(XTIME) 946 ! prints the CPU elapsed time 947 call xtimer_view(curr_time, "nucpot_hket_assign", STDOUT) 948#endif 949 return 950 end subroutine nucpot_hket_assign 951