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 bra 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 bra center in nuclear attraction potential integrals 28 !> \author Bin Gao 29 !> \date 2011-10-18 30 !> \param orders_hgto_bra contains the minimum and maximum orders of HGTOs on bra center to return 31 !> \param orders_hgto_ket contains the minimum and maximum orders of HGTOs on ket center to return 32 !> \param orders_geo_pot contains the minimum and maximum orders of geometric derivatives to return 33 !> \param coord_ket contains the coordinates of bra center 34 !> \param exponent_ket is the exponent of HGTOs of bra center 35 !> \param coord_ket contains the coordinates of bra center 36 !> \param exponent_ket is the exponent of HGTOs of bra center 37 !> \param dim_hket is the dimension of HGTOs on ket center 38 !> \param dim_geo_hket is the dimension of geometric derivatives on nuclear potential origin 39 !> \param hket_pints contains the nuclear attraction integrals with zeroth order 40 !> HGTOs on bra center 41 !> \param dim_hgto_bra is the dimension of HGTOs of bra center afterwards 42 !> \param dim_hgto_ket is the dimension of HGTOs of ket center afterwards 43 !> \param dim_geo_hbra is the dimension of geometric derivatives afterwards 44 !> \return hbra_pints contains the integrals with required HGTOs on bra center 45 subroutine nucpot_hbra(orders_hgto_bra, orders_hgto_ket, orders_geo_pot, & 46 coord_bra, exponent_bra, coord_ket, exponent_ket, & 47 dim_hket, dim_geo_hket, hket_pints, dim_hgto_bra, & 48 dim_hgto_ket, dim_geo_hbra, hbra_pints) 49 use xkind 50 implicit none 51 integer, intent(in) :: orders_hgto_bra(2) 52 integer, intent(in) :: orders_hgto_ket(2) 53 integer, intent(in) :: orders_geo_pot(2) 54 real(REALK), intent(in) :: coord_bra(3) 55 real(REALK), intent(in) :: exponent_bra 56 real(REALK), intent(in) :: coord_ket(3) 57 real(REALK), intent(in) :: exponent_ket 58 integer, intent(in) :: dim_hket 59 integer, intent(in) :: dim_geo_hket 60 real(REALK), intent(in) :: hket_pints(dim_hket,dim_geo_hket) 61 integer, intent(in) :: dim_hgto_bra 62 integer, intent(in) :: dim_hgto_ket 63 integer, intent(in) :: dim_geo_hbra 64 real(REALK), intent(out) :: hbra_pints(dim_hgto_bra,dim_hgto_ket,dim_geo_hbra) 65!f2py intent(in) :: orders_hgto_bra 66!f2py intent(in) :: orders_hgto_ket 67!f2py intent(in) :: orders_geo_pot 68!f2py intent(in) :: coord_bra 69!f2py intent(in) :: exponent_bra 70!f2py intent(in) :: coord_ket 71!f2py intent(in) :: exponent_ket 72!f2py intent(hide) :: dim_hket 73!f2py intent(hide) :: dim_geo_hket 74!f2py intent(in) :: hket_pints 75!f2py intent(in) :: dim_hgto_bra 76!f2py intent(in) :: dim_hgto_ket 77!f2py intent(in) :: dim_geo_hbra 78!f2py intent(out) :: hbra_pints 79!f2py depend(dim_hgto_bra) :: hbra_pints 80!f2py depend(dim_hgto_ket) :: hbra_pints 81!f2py depend(dim_geo_hbra) :: hbra_pints 82 real(REALK) half_neg_rp !half of the negative reciprocal of total exponent \f$p_{ij}\f$ 83 real(REALK) ket_to_bra !ratio of exponent on ket center to that on bra center 84 real(REALK) cc_wrt_bra(3) !relative coordinates of center-of-charge w.r.t. bra center 85 integer max_low_geo !maximum of lower order of geometric derivatives 86 integer num_up_geo !number of higher order geometric derivatives 87 integer num_low_geo !number of lower order geometric derivatives 88 integer start_up_geo !start address of upper order geometric derivatives 89 integer end_up_geo !end address of upper order geometric derivatives 90 integer start_low_geo !start address of lower order geometric derivatives 91 integer end_low_geo !end address of lower order geometric derivatives 92 integer min_order_hket !minimum of current order of HGTOs on ket center 93 integer max_order_hbra !maximum of current order of HGTOs on bra center 94 integer dim_cur_hket !dimension of current order HGTOs on ket cetner of temporary integrals 95 integer dim_hket_zero !dimension of HGTOs on ket center with zeroth order HGTO on bra center 96 integer dim_hket_low !dimension of HGTOs on ket center with lower order HGTOs on bra center 97 integer dim_low_hbra !dimension of lower order HGTOs on bra center of temporary integrals 98 integer dim_up_hbra !dimension of upper order HGTOs on bra center of temporary integrals 99 integer order_geo !incremental recorder over orders of geometric derivatives 100 integer max_dim_hket !maximum dimension of HGTOs on ket center 101 integer size_low_hbra !size of temporary integrals of lower order HGTOs on bra center 102 integer size_up_hbra !size of temporary integrals of upper order HGTOs on bra center 103 integer size_bra_ket !size of HGTOs on bra and ket centers 104 integer low_hbra_int !pointer to temporary integrals of lower order HGTOs on bra center 105 integer up_hbra_int !pointer to temporary integrals of upper order HGTOs on bra center 106 integer dim_tmp !dimension of temporary integrals 107 real(REALK), allocatable :: tmp_ints(:,:) 108 !temporary integrals 109 integer offset_hbra_geo !offset of geometric derivatives in returned integrals 110 integer offset_hgto_ket !offset of HGTOs on ket center in integrals with zeroth 111 !order HGTO on bra center 112 logical zero_hbra !if returning zeroth order HGTO on bra center 113 integer offset_hket_low !offset of HGTOs on ket center in temporary integrals with 114 !lower order HGTOs on bra center 115 integer start_hket_up !start addresses of HGTOs on ket center in integrals with zeroth order 116 integer start_hket_low !HGTO on bra center, and upper/lower order geometric derivatives 117 integer ierr !error information 118#if defined(XTIME) 119 real(REALK) curr_time !current CPU time 120 ! sets current CPU time 121 call xtimer_set(curr_time) 122#endif 123 select case(orders_hgto_bra(2)) 124 ! returns s-shell HGTO 125 case(0) 126 hbra_pints(1,:,:) = hket_pints 127 ! maximum returned HGTOs are other shells, at least p-shell 128 case default 129 ! calculates the relative coordinates of center-of-charge w.r.t. bra center, 130 ! and the half of negative reciprocal of total exponent 131 half_neg_rp = 1.0_REALK/(exponent_bra+exponent_ket) 132 do ierr = 1, 3 133 cc_wrt_bra(ierr) = half_neg_rp*(exponent_bra*coord_bra(ierr) & 134 + exponent_ket*coord_ket(ierr))-coord_bra(ierr) 135 end do 136 half_neg_rp = -0.5_REALK*half_neg_rp 137 ! sets the ratio of exponent on ket center to that on bra center 138 ket_to_bra = exponent_ket/exponent_bra 139 ! maximum dimension of HGTOs on ket center 140 max_dim_hket = (orders_hgto_ket(2)+1)*(orders_hgto_ket(2)+2) & 141 * (orders_hgto_ket(2)+3)/6 142 ! sets the maximum of lower order of geometric derivatives 143 max_low_geo = orders_geo_pot(2)+orders_hgto_bra(2) 144 ! allocates memory for temporary integrals 145 call dim_nucpot_hbra(orders_hgto_ket(1), max_low_geo, orders_geo_pot(2), & 146 max_dim_hket, dim_tmp) 147 allocate(tmp_ints(dim_tmp,2), stat=ierr) 148 if (ierr/=0) & 149 call error_stop("nucpot_hbra", "failed to allocate tmp_ints", & 150 dim_tmp*2) 151 ! \var(max_low_geo)-1 order geometric derivatives 152 end_up_geo = dim_geo_hket !end address of upper order geometric derivatives 153 num_up_geo = (max_low_geo+1)*(max_low_geo+2)/2 !number of upper order geometric derivatives 154 end_low_geo = end_up_geo-num_up_geo !end address of lower order geometric derivatives 155 start_up_geo = end_low_geo+1 !start address of upper order geometric derivatives 156 max_low_geo = max_low_geo-1 157 num_low_geo = num_up_geo-(max_low_geo+2) !number of lower order geometric derivatives 158 start_low_geo = end_low_geo-num_low_geo+1 !end address of \var(orders_geo_pot(2)) order derivatives 159#if defined(DEBUG) 160 write(STDOUT,100) "GEO-s-HGTO/upper/start/end:", & 161 max_low_geo+1, start_up_geo, end_up_geo 162 write(STDOUT,100) "GEO-s-HGTO/upper/start/end:", & 163 max_low_geo, start_low_geo, end_low_geo 164#endif 165 ! sets the minimum of current order of HGTOs on ket center 166 min_order_hket = orders_hgto_ket(1) 167 ! sets the dimensions of HGTOs on ket center 168 if (min_order_hket==0) then 169 dim_cur_hket = max_dim_hket 170 dim_hket_zero = dim_cur_hket 171 else 172 ierr = min_order_hket*(min_order_hket+1)/2 173 dim_cur_hket = max_dim_hket-ierr*(min_order_hket+2)/3 174 dim_hket_zero = dim_cur_hket+ierr 175 end if 176 ! initializes the maximum of current order of HGTOs on bra center 177 max_order_hbra = 0 178 ! sets the dimensions and size of temporary integrals with lower order 179 ! HGTOs on bra center (not used in recurrence relations) 180 dim_low_hbra = 1 181 dim_hket_low = dim_cur_hket 182 size_low_hbra = dim_hket_low*num_up_geo 183 ! sets the dimension and size of temporary integrals with upper order HGTOs on bra center 184 dim_up_hbra = 3 185 size_bra_ket = dim_up_hbra*dim_cur_hket 186 size_up_hbra = size_bra_ket*num_low_geo 187 ! sets the start addresses of HGTOs on ket center in integrals with zeroth 188 ! order HGTO on bra center, and upper/lower order geometric derivatives 189 start_hket_up = dim_hket-dim_cur_hket+1 190 start_hket_low = dim_hket-dim_hket_zero+1 191 ! gets the temporary integrals for order of HGTOs on bra center up to 192 ! \var(max_order_hbra)+1 and order of geometric derivatives as \var(max_low_geo) 193 call sub_nucpot_hbra(max_low_geo, orders_hgto_ket, min_order_hket, & 194 max_order_hbra, cc_wrt_bra, half_neg_rp, & 195 ket_to_bra, dim_cur_hket, num_up_geo, & 196 hket_pints(start_hket_up:dim_hket,start_up_geo:end_up_geo), & 197 dim_hket_zero, num_low_geo, & 198 hket_pints(start_hket_low:dim_hket,start_low_geo:end_low_geo), & 199 0, dim_low_hbra, dim_hket_low, tmp_ints(1:size_low_hbra,1), & 200 dim_up_hbra, tmp_ints(1:size_up_hbra,2)) 201 ! initializes the pointers of HGTOs on bra center 202 low_hbra_int = 1 203 up_hbra_int = 2 204 ! loops over the orders of geometric derivatives not returned, the maximum 205 ! of current order of HGTOs \var(max_order_hbra) needs to update in the cycle 206 do order_geo = max_low_geo-1, orders_geo_pot(2), -1 207 ! updates the numbers of lower and upper order geometric derivatives 208 num_up_geo = num_low_geo 209 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 210 ! updates the start and end addresses of lower and upper order geometric derivatives 211 end_up_geo = end_low_geo 212 start_up_geo = start_low_geo 213 end_low_geo = start_low_geo-1 214 start_low_geo = end_low_geo-num_low_geo+1 215#if defined(DEBUG) 216 write(STDOUT,100) "GEO-HGTO/loop/1/upper/start/end:", & 217 order_geo+1, start_up_geo, end_up_geo 218 write(STDOUT,100) "GEO-HGTO/loop/1/lower/start/end:", & 219 order_geo, start_low_geo, end_low_geo 220#endif 221 ! updates maximum of current order of HGTOs on bra center 222 max_order_hbra = max_order_hbra+1 223 ! updates the dimensions of HGTOs on bra center 224 dim_low_hbra = dim_up_hbra 225 dim_up_hbra = dim_low_hbra+(max_order_hbra+2)*(max_order_hbra+3)/2 226 ! sets minimum of current order of HGTOs on ket center 227 min_order_hket = max(orders_hgto_ket(1)-max_order_hbra,0) 228 ! sets the dimensions of HGTOs on ket center 229 dim_hket_low = dim_cur_hket 230 if (min_order_hket==0) then 231 dim_cur_hket = max_dim_hket 232 dim_hket_zero = dim_cur_hket 233 else 234 ierr = min_order_hket*(min_order_hket+1)/2 235 dim_cur_hket = max_dim_hket-ierr*(min_order_hket+2)/3 236 dim_hket_zero = dim_cur_hket+ierr 237 end if 238 ! updates the sizes of temporary integrals 239 size_low_hbra = size_up_hbra 240 size_bra_ket = dim_up_hbra*dim_cur_hket 241 size_up_hbra = size_bra_ket*num_low_geo 242 ! switches the pointers 243 low_hbra_int = 3-low_hbra_int 244 up_hbra_int = 3-up_hbra_int 245 ! sets the start addresses of HGTOs on ket center in integrals with zeroth 246 ! order HGTO on bra center, and upper/lower order geometric derivatives 247 start_hket_up = dim_hket-dim_cur_hket+1 248 start_hket_low = dim_hket-dim_hket_zero+1 249 ! gets the temporary integrals for order of HGTOs on bra center up to 250 ! \var(max_order_hbra)+1 and order of geometric derivatives as \var(order_geo) 251 call sub_nucpot_hbra(order_geo, orders_hgto_ket, min_order_hket, & 252 max_order_hbra, cc_wrt_bra, half_neg_rp, & 253 ket_to_bra, dim_cur_hket, num_up_geo, & 254 hket_pints(start_hket_up:dim_hket,start_up_geo:end_up_geo), & 255 dim_hket_zero, num_low_geo, & 256 hket_pints(start_hket_low:dim_hket,start_low_geo:end_low_geo), & 257 0, dim_low_hbra, dim_hket_low, & 258 tmp_ints(1:size_low_hbra,low_hbra_int), & 259 dim_up_hbra, tmp_ints(1:size_up_hbra,up_hbra_int)) 260 end do 261 ! if returing zero order HGTO 262 zero_hbra = orders_hgto_bra(1)==0 263 ! sets the offsets of geometric derivatives and HGTOs 264 offset_hbra_geo = dim_geo_hbra-num_low_geo 265 offset_hgto_ket = dim_hket-dim_hgto_ket 266 ! assigns the returned integrals 267 call nucpot_hbra_assign(offset_hgto_ket, start_low_geo-1, dim_hket, & 268 dim_geo_hket, hket_pints, dim_up_hbra, & 269 dim_cur_hket, num_low_geo, & 270 tmp_ints(1:size_up_hbra,up_hbra_int), & 271 zero_hbra, offset_hbra_geo, dim_hgto_bra, & 272 dim_hgto_ket, dim_geo_hbra, hbra_pints) 273 ! sets offset of HGTOs on ket center in temporary integrals with lower order HGTOs on bra center 274 offset_hket_low = dim_cur_hket-dim_hket_low 275 ! loops over other returned orders of geometric derivatives, maximum of current 276 ! order of HGTOs on bra center \var(max_order_hbra) does not need to update 277 do order_geo = orders_geo_pot(2)-1, orders_geo_pot(1), -1 278 ! updates the numbers of lower and upper order geometric derivatives 279 num_up_geo = num_low_geo 280 num_low_geo = num_low_geo-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 281 ! updates the start and end addresses of lower and upper order geometric derivatives 282 end_up_geo = end_low_geo 283 start_up_geo = start_low_geo 284 end_low_geo = start_low_geo-1 285 start_low_geo = end_low_geo-num_low_geo+1 286#if defined(DEBUG) 287 write(STDOUT,100) "GEO-HGTO/loop/2/upper/start/end:", & 288 order_geo+1, start_up_geo, end_up_geo 289 write(STDOUT,100) "GEO-HGTO/loop/2/lower/start/end:", & 290 order_geo, start_low_geo, end_low_geo 291#endif 292 ! updates the sizes of temporary integrals 293 size_low_hbra = size_up_hbra 294 size_up_hbra = size_bra_ket*num_low_geo 295 ! switches the pointers 296 low_hbra_int = 3-low_hbra_int 297 up_hbra_int = 3-up_hbra_int 298 ! gets the temporary integrals for order of HGTOs on bra center up to 299 ! \var(max_order_hbra)+1 and order of geometric derivatives as \var(order_geo) 300 call sub_nucpot_hbra(order_geo, orders_hgto_ket, min_order_hket, & 301 max_order_hbra, cc_wrt_bra, half_neg_rp, & 302 ket_to_bra, dim_cur_hket, num_up_geo, & 303 hket_pints(start_hket_up:dim_hket,start_up_geo:end_up_geo), & 304 dim_hket_zero, num_low_geo, & 305 hket_pints(start_hket_low:dim_hket,start_low_geo:end_low_geo), & 306 offset_hket_low, dim_up_hbra, dim_cur_hket, & 307 tmp_ints(1:size_low_hbra,low_hbra_int), & 308 dim_up_hbra, tmp_ints(1:size_up_hbra,up_hbra_int)) 309 ! updates the offset of geometric derivatives 310 offset_hbra_geo = offset_hbra_geo-num_low_geo 311 ! assigns the returned integrals 312 call nucpot_hbra_assign(offset_hgto_ket, start_low_geo-1, dim_hket, & 313 dim_geo_hket, hket_pints, dim_up_hbra, & 314 dim_cur_hket, num_low_geo, & 315 tmp_ints(1:size_up_hbra,up_hbra_int), & 316 zero_hbra, offset_hbra_geo, dim_hgto_bra, & 317 dim_hgto_ket, dim_geo_hbra, hbra_pints) 318 end do 319 deallocate(tmp_ints) 320 end select 321#if defined(XTIME) 322 ! prints the CPU elapsed time 323 call xtimer_view(curr_time, "nucpot_hbra", STDOUT) 324#endif 325 return 326#if defined(DEBUG) 327100 format("nucpot_hbra>> ",A,I6,2I8) 328#endif 329 end subroutine nucpot_hbra 330 331 !> \brief gets the maximum dimension of temporary integrals used in recurrence relations 332 !> \author Bin Gao 333 !> \date 2012-03-05 334 !> \param min_order_hket is the minimum order of HGTOs on ket center 335 !> \param max_order_geo is the maximum order of geometric derivatives 336 !> \param min_order_geo is the minimum order of geometric derivatives 337 !> \param max_dim_hket is the maximum dimension of HGTOs on ket center 338 !> \return dim_ints is the maximum dimension of temporary integrals 339 subroutine dim_nucpot_hbra(min_order_hket, max_order_geo, min_order_geo, & 340 max_dim_hket, dim_ints) 341 use xkind 342 implicit none 343 integer, intent(in) :: min_order_hket 344 integer, intent(in) :: max_order_geo 345 integer, intent(in) :: min_order_geo 346 integer, intent(in) :: max_dim_hket 347 integer, intent(out) :: dim_ints 348!f2py intent(in) :: min_order_hket 349!f2py intent(in) :: max_order_geo 350!f2py intent(in) :: min_order_geo 351!f2py intent(in) :: max_dim_hket 352!f2py intent(out) :: dim_ints 353 integer num_geo_pot !number of upper order geometric derivatives 354 integer max_hgto_bra !maximum order of HGTOs on bra center 355 integer dim_hgto_bra !dimension of HGTOs on bra center 356 integer order_geo !incremental recorder over orders of geometric derivatives 357 integer min_hgto_ket !minimum order of HGTOs on ket center in the recurrence relations 358 integer dim_hgto_ket !dimension of HGTOs on ket center 359 integer dim_tmp !temporary result of dimension 360#if defined(XTIME) 361 real(REALK) curr_time !current CPU time 362 ! sets current CPU time 363 call xtimer_set(curr_time) 364#endif 365 ! initializes the return value 366 dim_ints = 0 367 ! initializes the number of geometric derivatives 368 num_geo_pot = (max_order_geo+1)*(max_order_geo+2)/2 369 ! initializes the maximum order and dimension of HGTOs on bra center 370 max_hgto_bra = 0 371 dim_hgto_bra = 0 372 ! loops over the orders of geometric derivatives not returned, the maximum 373 ! order of HGTOs \var(max_hgto_bra) needs to update in the cycle 374 do order_geo = max_order_geo-1, min_order_geo, -1 375 ! updates the number of geometric derivatives 376 num_geo_pot = num_geo_pot-(order_geo+2) !=(order_geo+1)*(order_geo+2)/2 377 ! updates maximum order of HGTOs on bra center 378 max_hgto_bra = max_hgto_bra+1 379 ! updates the dimension of HGTOs on bra center 380 dim_hgto_bra = dim_hgto_bra+(max_hgto_bra+1)*(max_hgto_bra+2)/2 381 ! sets minimum order of HGTOs on ket center 382 min_hgto_ket = max(min_order_hket-max_hgto_bra,0) 383 ! sets the dimension of HGTOs on ket center 384 if (min_hgto_ket==0) then 385 dim_hgto_ket = max_dim_hket 386 else 387 dim_hgto_ket = max_dim_hket-min_hgto_ket*(min_hgto_ket+1)*(min_hgto_ket+2)/6 388 end if 389 ! updates the maximum dimension 390 dim_tmp = dim_hgto_bra*dim_hgto_ket*num_geo_pot 391 if (dim_tmp>dim_ints) dim_ints = dim_tmp 392 end do 393#if defined(XTIME) 394 ! prints the CPU elapsed time 395 call xtimer_view(curr_time, "dim_nucpot_hbra", STDOUT) 396#endif 397 return 398 end subroutine dim_nucpot_hbra 399 400 !> \brief sub-recurrence relations by recovering upper order HGTOs on bra center 401 !> \author Bin Gao 402 !> \date 2011-10-18 403 !> \param cur_order_geo is current order of geometric derivatives 404 !> \param orders_hgto_ket contains the minimum and maximum orders HGTOs on ket center 405 !> \param min_order_hket is minimum of current order of HGTOs on ket center 406 !> \param max_order_hbra is maximum of current order of HGTOs on bra center 407 !> \param cc_wrt_bra contains the relative coordinates of center-of-charge w.r.t. bra center 408 !> \param half_neg_rp is the half of the negative reciprocal of total exponent \f$p_{ij}\f$ 409 !> \param ket_to_bra is the ratio of exponent on ket center to that on bra center 410 !> \param dim_cur_hket is the dimension of HGTOs on ket center for integrals with 411 !> upper order geometric derivatives and zeroth order HGTOs on bra center 412 !> \param num_up_geo is the number of upper order geometric derivatives 413 !> \param up_geo_zero contains the integrals with zeroth order HGTO on bra center and 414 !> upper order geometric derivatives 415 !> \param dim_hket_zero is the dimension of HGTOs on ket center for integrals with 416 !> lower order geometric derivatives and zeroth order HGTO on bra center 417 !> \param num_low_geo is the number lower order geometric derivatives 418 !> \param low_geo_zero contains the integrals with zeroth order HGTO on bra center and 419 !> lower order geometric derivatives 420 !> \param offset_hket_low is the offset of HGTOs on ket center in temporary integrals with 421 !> lower order HGTOs on bra center 422 !> \param dim_low_hbra is the dimension of lower order HGTOs on bra center 423 !> \param dim_hket_low is the dimension of HGTOs on ket center with upper order 424 !> geometric derivatives and lower order HGTOs on bra center 425 !> \param low_hbra_pints contains the integrals with lower order HGTOs on bra center and 426 !> upper order geometric derivatives 427 !> \param dim_up_hbra is the dimension of upper order HGTOs on bra center 428 !> \return up_hbra_pints contains the integrals with upper order HGTOs on bra center and 429 !> lower order geometric derivatives 430 subroutine sub_nucpot_hbra(cur_order_geo, orders_hgto_ket, min_order_hket, & 431 max_order_hbra, cc_wrt_bra, half_neg_rp, & 432 ket_to_bra, dim_cur_hket, num_up_geo, & 433 up_geo_zero, dim_hket_zero, num_low_geo, & 434 low_geo_zero, offset_hket_low, dim_low_hbra, & 435 dim_hket_low, low_hbra_pints, dim_up_hbra, up_hbra_pints) 436 use xkind 437 implicit none 438 integer, intent(in) :: cur_order_geo 439 integer, intent(in) :: orders_hgto_ket(2) 440 integer, intent(in) :: min_order_hket 441 integer, intent(in) :: max_order_hbra 442 real(REALK), intent(in) :: cc_wrt_bra(3) 443 real(REALK), intent(in) :: half_neg_rp 444 real(REALK), intent(in) :: ket_to_bra 445 integer, intent(in) :: dim_cur_hket 446 integer, intent(in) :: num_up_geo 447 real(REALK), intent(in) :: up_geo_zero(dim_cur_hket,num_up_geo) 448 integer, intent(in) :: dim_hket_zero 449 integer, intent(in) :: num_low_geo 450 real(REALK), intent(in) :: low_geo_zero(dim_hket_zero,num_low_geo) 451 integer, intent(in) :: offset_hket_low 452 integer, intent(in) :: dim_low_hbra 453 integer, intent(in) :: dim_hket_low 454 real(REALK), intent(in) :: low_hbra_pints(dim_low_hbra,dim_hket_low,num_up_geo) 455 integer, intent(in) :: dim_up_hbra 456 real(REALK), intent(out) :: up_hbra_pints(dim_up_hbra,dim_cur_hket,num_low_geo) 457!f2py intent(in) :: cur_order_geo 458!f2py intent(in) :: orders_hgto_ket 459!f2py intent(in) :: min_order_hket 460!f2py intent(in) :: max_order_hbra 461!f2py intent(in) :: cc_wrt_bra 462!f2py intent(in) :: half_neg_rp 463!f2py intent(in) :: ket_to_bra 464!f2py intent(hide) :: dim_cur_hket 465!f2py intent(hide) :: num_up_geo 466!f2py intent(in) :: up_geo_zero 467!f2py intent(hide) :: dim_hket_zero 468!f2py intent(hide) :: num_low_geo 469!f2py intent(in) :: low_geo_zero 470!f2py intent(in) :: offset_hket_low 471!f2py intent(hide) :: dim_low_hbra 472!f2py intent(hide) :: dim_hket_low 473!f2py intent(in) :: low_hbra_pints 474!f2py depend(num_up_geo) :: low_hbra_pints 475!f2py intent(in) :: dim_up_hbra 476!f2py intent(out) :: up_hbra_pints 477!f2py depend(dim_up_hbra) :: up_hbra_pints 478!f2py depend(dim_cur_hket) :: up_hbra_pints 479!f2py depend(num_low_geo) :: up_hbra_pints 480 integer addr_up_geo !addresses of upper order geometric derivatives 481 integer addr_up_geo_y 482 integer addr_up_geo_z 483 integer addr_cur_geo !address of current order geometric derivatives 484 integer igeo, jgeo !incremental recorders over geometric derivatives 485 integer min_cur_hket !minimum of current order of HGTOs on ket center 486 logical zero_cur_hket !if the minimum order of HGTOs on ket center is zeroth order 487 integer max_cur_hbra !maximum of current order of HGTOs on bra center for given order HGTOs on ket center 488 integer base_low_hket !base address of lower order HGTOs on ket center 489 integer base_hket_zero !base address of lower order HGTOs on ket center in integrals with zeroth 490 !order HGTO on bra center and lower order geometric derivatives 491 integer addr_cur_hket !address of current order HGTOs on ket center 492 integer addr_hket_zero !address of current order HGTOs on ket center in integrals with zeroth 493 !order HGTO on bra center and lower order geometric derivatives 494 integer addr_hket_low !address of current order HGTOs on ket center in integrals with lower 495 !order HGTOs on bra center and lower order geometric derivatives 496 integer addr_low_xket !addresses of lower order HGTOs on ket center 497 integer addr_low_yket 498 integer addr_low_zket 499 integer addr_xket_zero !addresses of lower order HGTOs on ket center in integrals with zeroth 500 integer addr_yket_zero !order HGTO on bra center and lower order geometric derivatives 501 integer addr_zket_zero 502 integer order_hket !order of HGTOs on ket center 503 integer order_xket !orders of xyz components of HGTOs on ket center 504 integer order_yket 505 integer order_zket 506 integer addr_up_hbra !address of upper order HGTOs on bra center 507 integer addr_cur_hbra !address of current order HGTOs on bra center 508 integer addr_low_hbra !address of lower order HGTOs on bra center 509 integer order_hbra !order of HGTOs on bra center 510 integer ibra, jbra !incremental recorder over HGTOs on bra center 511#if defined(XTIME) 512 real(REALK) curr_time !current CPU time 513 ! sets current CPU time 514 call xtimer_set(curr_time) 515#endif 516 if (min_order_hket==0) then 517 zero_cur_hket = .true. 518 min_cur_hket = 1 519 else 520 zero_cur_hket = .false. 521 min_cur_hket = min_order_hket 522 end if 523 addr_up_geo = 0 524 addr_cur_geo = 0 525 ! loops over xyz components of geometric derivatives 526 do igeo = cur_order_geo, 0, -1 527 do jgeo = 0, igeo 528 addr_up_geo = addr_up_geo+1 529 addr_cur_geo = addr_cur_geo+1 530 addr_cur_hket = 0 531 addr_hket_low = offset_hket_low 532 ! zeroth order HGTOs on ket center 533 if (zero_cur_hket) then 534 addr_cur_hket = addr_cur_hket+1 535 addr_hket_zero = 1 536 ! sets maximum of current order of HGTOs on bra center for zeroth order HGTOs on ket center 537 max_cur_hbra = max_order_hbra-orders_hgto_ket(1) 538 ! px HGTO on bra center 539 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 540 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 541 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo) 542 ! py HGTO on bra center 543 addr_up_geo_y = addr_up_geo+1 544 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 545 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 546 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_y) 547 ! pz HGTO on bra center 548 addr_up_geo_z = addr_up_geo+igeo+2 549 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 550 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 551 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_z) 552 ! d-shell on bra center 553 if (max_cur_hbra>0) then 554 addr_hket_low = addr_hket_low+1 555 ! dxx 556 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 557 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 558 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 559 + low_hbra_pints(1,addr_hket_low,addr_up_geo)) 560 ! dxy 561 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 562 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 563 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_y) 564 ! dyy 565 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 566 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 567 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 568 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y)) 569 ! dxz 570 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 571 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 572 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_z) 573 ! dyz 574 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 575 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 576 + half_neg_rp*low_hbra_pints(2,addr_hket_low,addr_up_geo_z) 577 ! dzz 578 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 579 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 580 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 581 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z)) 582 end if 583 if (max_cur_hbra>1) then 584 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 585 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 586 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 587 ! loops over other current order of HGTOs on bra center, starting from d-shell 588 do order_hbra = 2, max_cur_hbra 589 addr_up_hbra = addr_up_hbra+1 590 addr_cur_hbra = addr_cur_hbra+1 591#if defined(DEBUG) 592 write(STDOUT,100) "orders:", order_hbra, order_hket, cur_order_geo 593 write(STDOUT,100) "x-direction" 594 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 595 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 596 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 597 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo 598 write(STDOUT,100) "------------------------------------" 599#endif 600 ! x-direction 601 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 602 = cc_wrt_bra(1) & 603 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 604 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 605 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 606 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo)) 607 ! y-direction 608 addr_up_hbra = addr_up_hbra+1 609#if defined(DEBUG) 610 write(STDOUT,100) "y-direction" 611 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 612 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 613 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 614 write(STDOUT,100) "------------------------------------" 615#endif 616 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 617 = cc_wrt_bra(2) & 618 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 619 + half_neg_rp & 620 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) 621 do ibra = 1, order_hbra 622 addr_up_hbra = addr_up_hbra+1 623#if defined(DEBUG) 624 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 625 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 626 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 627 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 628 write(STDOUT,100) "------------------------------------" 629#endif 630 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 631 = cc_wrt_bra(2) & 632 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 633 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 634 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 635 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1)) 636 end do 637 ! z-direction 638#if defined(DEBUG) 639 write(STDOUT,100) "z-direction" 640#endif 641 addr_cur_hbra = addr_cur_hbra-1 642 do jbra = 0, order_hbra 643 addr_up_hbra = addr_up_hbra+1 644 addr_cur_hbra = addr_cur_hbra+1 645#if defined(DEBUG) 646 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 647 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 648 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 649 write(STDOUT,100) "------------------------------------" 650#endif 651 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 652 = cc_wrt_bra(3) & 653 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 654 + half_neg_rp & 655 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) 656 end do 657 do ibra = 1, order_hbra 658 do jbra = 0, order_hbra-ibra 659 addr_up_hbra = addr_up_hbra+1 660 addr_cur_hbra = addr_cur_hbra+1 661 addr_low_hbra = addr_low_hbra+1 662#if defined(DEBUG) 663 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 664 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 665 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 666 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 667 write(STDOUT,100) "------------------------------------" 668#endif 669 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 670 = cc_wrt_bra(3) & 671 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 672 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 673 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 674 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2)) 675 end do 676 end do 677 end do 678 end if 679 else 680 addr_hket_zero = min_order_hket*(min_order_hket+1)/2 681 end if 682 base_low_hket = 0 683 base_hket_zero = 0 684 ! loops over other current order HGTOs on ket center 685 do order_hket = min_cur_hket, orders_hgto_ket(2) 686 addr_cur_hket = addr_cur_hket+1 687 addr_hket_zero = addr_hket_zero+1 688 addr_xket_zero = base_hket_zero+1 689 ! sets maximum of current order of HGTOs on bra center for order \var(order_hket) 690 ! HGTOs on ket center 691 max_cur_hbra = min(max_order_hbra,max_order_hbra-orders_hgto_ket(1)+order_hket) 692 ! (1) x...x component of upper order HGTOs on ket center 693 ! 694 ! px HGTO on bra center 695 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 696 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 697 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo) & 698 - real(order_hket,REALK)*low_geo_zero(addr_xket_zero,addr_cur_geo)) 699 ! py HGTO on bra center 700 addr_up_geo_y = addr_up_geo+1 701 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 702 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 703 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_y) 704 ! pz HGTO on bra center 705 addr_up_geo_z = addr_up_geo+igeo+2 706 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 707 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 708 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_z) 709 ! d-shell on bra center 710 addr_low_xket = base_low_hket+1 !should put after if statement, but to avoid warning from compiler 711 if (max_cur_hbra>0) then 712 addr_hket_low = addr_hket_low+1 713 ! dxx 714 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 715 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 716 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 717 + low_hbra_pints(1,addr_hket_low,addr_up_geo) & 718 - real(order_hket,REALK)*up_hbra_pints(1,addr_low_xket,addr_cur_geo)) 719 ! dxy 720 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 721 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 722 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_y) 723 ! dyy 724 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 725 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 726 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 727 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y)) 728 ! dxz 729 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 730 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 731 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_z) 732 ! dyz 733 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 734 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 735 + half_neg_rp*low_hbra_pints(2,addr_hket_low,addr_up_geo_z) 736 ! dzz 737 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 738 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 739 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 740 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z)) 741 if (max_cur_hbra>1) then 742 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 743 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 744 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 745 ! loops over other current order of HGTOs on bra center, starting from d-shell 746 do order_hbra = 2, max_cur_hbra 747 addr_up_hbra = addr_up_hbra+1 748 addr_cur_hbra = addr_cur_hbra+1 749#if defined(DEBUG) 750 write(STDOUT,100) "orders:", order_hbra, order_hket, cur_order_geo 751 write(STDOUT,100) "x-direction" 752 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 753 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 754 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 755 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo 756 write(STDOUT,110) addr_up_hbra, addr_low_xket, addr_cur_geo 757 write(STDOUT,100) "------------------------------------" 758#endif 759 ! x-direction 760 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 761 = cc_wrt_bra(1) & 762 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 763 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 764 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 765 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo) & 766 - real(order_hket,REALK) & 767 * up_hbra_pints(addr_cur_hbra,addr_low_xket,addr_cur_geo)) 768 ! y-direction 769 addr_up_hbra = addr_up_hbra+1 770#if defined(DEBUG) 771 write(STDOUT,100) "y-direction" 772 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 773 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 774 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 775 write(STDOUT,100) "------------------------------------" 776#endif 777 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 778 = cc_wrt_bra(2) & 779 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 780 + half_neg_rp & 781 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) 782 do ibra = 1, order_hbra 783 addr_up_hbra = addr_up_hbra+1 784#if defined(DEBUG) 785 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 786 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 787 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 788 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 789 write(STDOUT,100) "------------------------------------" 790#endif 791 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 792 = cc_wrt_bra(2) & 793 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 794 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 795 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 796 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1)) 797 end do 798 ! z-direction 799#if defined(DEBUG) 800 write(STDOUT,100) "z-direction" 801#endif 802 addr_cur_hbra = addr_cur_hbra-1 803 do jbra = 0, order_hbra 804 addr_up_hbra = addr_up_hbra+1 805 addr_cur_hbra = addr_cur_hbra+1 806#if defined(DEBUG) 807 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 808 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 809 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 810 write(STDOUT,100) "------------------------------------" 811#endif 812 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 813 = cc_wrt_bra(3) & 814 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 815 + half_neg_rp & 816 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) 817 end do 818 do ibra = 1, order_hbra 819 do jbra = 0, order_hbra-ibra 820 addr_up_hbra = addr_up_hbra+1 821 addr_cur_hbra = addr_cur_hbra+1 822 addr_low_hbra = addr_low_hbra+1 823#if defined(DEBUG) 824 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 825 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 826 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 827 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 828 write(STDOUT,100) "------------------------------------" 829#endif 830 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 831 = cc_wrt_bra(3) & 832 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 833 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 834 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 835 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2)) 836 end do 837 end do 838 end do 839 end if 840 end if 841 ! (2) x...xy to xy...y components of upper order HGTOs on ket center 842 addr_low_yket = base_low_hket 843 addr_yket_zero = base_hket_zero 844 do order_yket = 1, order_hket-1 845 addr_cur_hket = addr_cur_hket+1 846 addr_hket_zero = addr_hket_zero+1 847 addr_xket_zero = addr_xket_zero+1 848 addr_yket_zero = addr_yket_zero+1 849 ! sets the order along x-direction 850 order_xket = order_hket-order_yket 851 ! px HGTO on bra center 852 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 853 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 854 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo) & 855 - real(order_xket,REALK)*low_geo_zero(addr_xket_zero,addr_cur_geo)) 856 ! py HGTO on bra center 857 addr_up_geo_y = addr_up_geo+1 858 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 859 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 860 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_y) & 861 - real(order_yket,REALK)*low_geo_zero(addr_yket_zero,addr_cur_geo)) 862 ! pz HGTO on bra center 863 addr_up_geo_z = addr_up_geo+igeo+2 864 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 865 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 866 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_z) 867 ! d-shell on bra center 868 if (max_cur_hbra>0) then 869 addr_low_xket = addr_low_xket+1 870 addr_low_yket = addr_low_yket+1 871 addr_hket_low = addr_hket_low+1 872 ! dxx 873 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 874 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 875 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 876 + low_hbra_pints(1,addr_hket_low,addr_up_geo) & 877 - real(order_xket,REALK)*up_hbra_pints(1,addr_low_xket,addr_cur_geo)) 878 ! dxy 879 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 880 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 881 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_y) & 882 - real(order_yket,REALK)*up_hbra_pints(1,addr_low_yket,addr_cur_geo)) 883 ! dyy 884 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 885 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 886 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 887 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y) & 888 - real(order_yket,REALK)*up_hbra_pints(2,addr_low_yket,addr_cur_geo)) 889 ! dxz 890 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 891 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 892 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_z) 893 ! dyz 894 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 895 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 896 + half_neg_rp*low_hbra_pints(2,addr_hket_low,addr_up_geo_z) 897 ! dzz 898 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 899 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 900 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 901 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z)) 902 if (max_cur_hbra>1) then 903 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 904 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 905 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 906 ! loops over other current order of HGTOs on bra center, starting from d-shell 907 do order_hbra = 2, max_cur_hbra 908 addr_up_hbra = addr_up_hbra+1 909 addr_cur_hbra = addr_cur_hbra+1 910#if defined(DEBUG) 911 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 912 write(STDOUT,100) "x-direction" 913 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 914 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 915 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 916 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo 917 write(STDOUT,110) addr_cur_hbra, addr_low_xket, addr_cur_geo 918 write(STDOUT,100) "------------------------------------" 919#endif 920 ! x-direction 921 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 922 = cc_wrt_bra(1) & 923 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 924 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 925 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 926 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo) & 927 - real(order_xket,REALK) & 928 * up_hbra_pints(addr_cur_hbra,addr_low_xket,addr_cur_geo)) 929 ! y-direction 930 addr_up_hbra = addr_up_hbra+1 931#if defined(DEBUG) 932 write(STDOUT,100) "y-direction" 933 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 934 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 935 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 936 write(STDOUT,110) addr_up_hbra, addr_low_yket, addr_cur_geo 937 write(STDOUT,100) "------------------------------------" 938#endif 939 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 940 = cc_wrt_bra(2) & 941 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 942 + half_neg_rp & 943 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) & 944 - real(order_yket,REALK) & 945 * up_hbra_pints(addr_cur_hbra,addr_low_yket,addr_cur_geo)) 946 do ibra = 1, order_hbra 947 addr_up_hbra = addr_up_hbra+1 948#if defined(DEBUG) 949 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 950 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 951 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 952 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 953 write(STDOUT,110) addr_cur_hbra+ibra, addr_low_yket, addr_cur_geo 954 write(STDOUT,100) "------------------------------------" 955#endif 956 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 957 = cc_wrt_bra(2) & 958 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 959 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 960 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 961 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1) & 962 - real(order_yket,REALK) & 963 * up_hbra_pints(addr_cur_hbra+ibra,addr_low_yket,addr_cur_geo)) 964 end do 965 ! z-direction 966#if defined(DEBUG) 967 write(STDOUT,100) "z-direction" 968#endif 969 addr_cur_hbra = addr_cur_hbra-1 970 do jbra = 0, order_hbra 971 addr_up_hbra = addr_up_hbra+1 972 addr_cur_hbra = addr_cur_hbra+1 973#if defined(DEBUG) 974 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 975 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 976 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 977 write(STDOUT,100) "------------------------------------" 978#endif 979 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 980 = cc_wrt_bra(3) & 981 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 982 + half_neg_rp & 983 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) 984 end do 985 do ibra = 1, order_hbra 986 do jbra = 0, order_hbra-ibra 987 addr_up_hbra = addr_up_hbra+1 988 addr_cur_hbra = addr_cur_hbra+1 989 addr_low_hbra = addr_low_hbra+1 990#if defined(DEBUG) 991 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 992 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 993 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 994 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 995 write(STDOUT,100) "------------------------------------" 996#endif 997 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 998 = cc_wrt_bra(3) & 999 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1000 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1001 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1002 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2)) 1003 end do 1004 end do 1005 end do 1006 end if 1007 end if 1008 end do 1009 ! (3) y...y component of upper order HGTOs on ket center 1010 addr_cur_hket = addr_cur_hket+1 1011 addr_hket_zero = addr_hket_zero+1 1012 addr_yket_zero = addr_yket_zero+1 1013 ! px HGTO on bra center 1014 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1015 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1016 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo) 1017 ! py HGTO on bra center 1018 addr_up_geo_y = addr_up_geo+1 1019 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1020 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1021 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_y) & 1022 - real(order_hket,REALK)*low_geo_zero(addr_yket_zero,addr_cur_geo)) 1023 ! pz HGTO on bra center 1024 addr_up_geo_z = addr_up_geo+igeo+2 1025 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1026 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1027 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_z) 1028 ! d-shell on bra center 1029 if (max_cur_hbra>0) then 1030 addr_low_yket = addr_low_yket+1 1031 addr_hket_low = addr_hket_low+1 1032 ! dxx 1033 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 1034 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1035 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1036 + low_hbra_pints(1,addr_hket_low,addr_up_geo)) 1037 ! dxy 1038 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 1039 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1040 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_y) & 1041 - real(order_hket,REALK)*up_hbra_pints(1,addr_low_yket,addr_cur_geo)) 1042 ! dyy 1043 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 1044 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1045 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1046 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y) & 1047 - real(order_hket,REALK)*up_hbra_pints(2,addr_low_yket,addr_cur_geo)) 1048 ! dxz 1049 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 1050 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1051 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_z) 1052 ! dyz 1053 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 1054 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1055 + half_neg_rp*low_hbra_pints(2,addr_hket_low,addr_up_geo_z) 1056 ! dzz 1057 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 1058 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1059 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1060 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z)) 1061 if (max_cur_hbra>1) then 1062 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 1063 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 1064 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 1065 ! loops over other current order of HGTOs on bra center, starting from d-shell 1066 do order_hbra = 2, max_cur_hbra 1067 addr_up_hbra = addr_up_hbra+1 1068 addr_cur_hbra = addr_cur_hbra+1 1069#if defined(DEBUG) 1070 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 1071 write(STDOUT,100) "x-direction" 1072 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1073 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1074 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 1075 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1076 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo 1077 write(STDOUT,100) "------------------------------------" 1078#endif 1079 ! x-direction 1080 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1081 = cc_wrt_bra(1) & 1082 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1083 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 1084 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 1085 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo)) 1086 ! y-direction 1087 addr_up_hbra = addr_up_hbra+1 1088#if defined(DEBUG) 1089 write(STDOUT,100) "y-direction" 1090 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1091 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1092 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 1093 write(STDOUT,110) addr_cur_hbra, addr_low_yket, addr_cur_geo 1094 write(STDOUT,100) "------------------------------------" 1095#endif 1096 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1097 = cc_wrt_bra(2) & 1098 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1099 + half_neg_rp & 1100 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) & 1101 - real(order_hket,REALK) & 1102 * up_hbra_pints(addr_cur_hbra,addr_low_yket,addr_cur_geo)) 1103 do ibra = 1, order_hbra 1104 addr_up_hbra = addr_up_hbra+1 1105#if defined(DEBUG) 1106 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1107 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 1108 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 1109 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 1110 write(STDOUT,110) addr_cur_hbra+ibra, addr_low_yket, addr_cur_geo 1111 write(STDOUT,100) "------------------------------------" 1112#endif 1113 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1114 = cc_wrt_bra(2) & 1115 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1116 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1117 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1118 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1) & 1119 - real(order_hket,REALK) & 1120 * up_hbra_pints(addr_cur_hbra+ibra,addr_low_yket,addr_cur_geo)) 1121 end do 1122 ! z-direction 1123#if defined(DEBUG) 1124 write(STDOUT,100) "z-direction" 1125#endif 1126 addr_cur_hbra = addr_cur_hbra-1 1127 do jbra = 0, order_hbra 1128 addr_up_hbra = addr_up_hbra+1 1129 addr_cur_hbra = addr_cur_hbra+1 1130#if defined(DEBUG) 1131 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1132 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1133 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1134 write(STDOUT,100) "------------------------------------" 1135#endif 1136 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1137 = cc_wrt_bra(3) & 1138 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1139 + half_neg_rp & 1140 * low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) 1141 end do 1142 do ibra = 1, order_hbra 1143 do jbra = 0, order_hbra-ibra 1144 addr_up_hbra = addr_up_hbra+1 1145 addr_cur_hbra = addr_cur_hbra+1 1146 addr_low_hbra = addr_low_hbra+1 1147#if defined(DEBUG) 1148 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1149 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1150 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 1151 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1152 write(STDOUT,100) "------------------------------------" 1153#endif 1154 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1155 = cc_wrt_bra(3) & 1156 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1157 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1158 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1159 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2)) 1160 end do 1161 end do 1162 end do 1163 end if 1164 end if 1165 ! (4) x...xz to yz...z components of upper order HGTOs on ket center 1166 addr_low_xket = base_low_hket+order_hket 1167 addr_low_zket = base_low_hket 1168 addr_xket_zero = base_hket_zero+order_hket 1169 addr_zket_zero = base_hket_zero 1170 do order_zket = 1, order_hket-1 1171 addr_cur_hket = addr_cur_hket+1 1172 addr_hket_zero = addr_hket_zero+1 1173 addr_xket_zero = addr_xket_zero+1 1174 addr_zket_zero = addr_zket_zero+1 1175 ! (4.1) x...xz...z component of upper order HGTOs on ket center 1176 order_xket = order_hket-order_zket 1177 ! px HGTO on bra center 1178 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1179 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1180 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo) & 1181 - real(order_xket,REALK)*low_geo_zero(addr_xket_zero,addr_cur_geo)) 1182 ! py HGTO on bra center 1183 addr_up_geo_y = addr_up_geo+1 1184 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1185 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1186 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_y) 1187 ! pz HGTO on bra center 1188 addr_up_geo_z = addr_up_geo+igeo+2 1189 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1190 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1191 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_z) & 1192 - real(order_zket,REALK)*low_geo_zero(addr_zket_zero,addr_cur_geo)) 1193 ! d-shell on bra center 1194 if (max_cur_hbra>0) then 1195 addr_low_xket = addr_low_xket+1 1196 addr_low_zket = addr_low_zket+1 1197 addr_hket_low = addr_hket_low+1 1198 ! dxx 1199 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 1200 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1201 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1202 + low_hbra_pints(1,addr_hket_low,addr_up_geo) & 1203 - real(order_xket,REALK)*up_hbra_pints(1,addr_low_xket,addr_cur_geo)) 1204 ! dxy 1205 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 1206 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1207 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_y) 1208 ! dyy 1209 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 1210 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1211 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1212 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y)) 1213 ! dxz 1214 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 1215 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1216 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_z) & 1217 - real(order_zket,REALK)*up_hbra_pints(1,addr_low_zket,addr_cur_geo)) 1218 ! dyz 1219 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 1220 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1221 + half_neg_rp*(low_hbra_pints(2,addr_hket_low,addr_up_geo_z) & 1222 - real(order_zket,REALK)*up_hbra_pints(2,addr_low_zket,addr_cur_geo)) 1223 ! dzz 1224 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 1225 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1226 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1227 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z) & 1228 - real(order_zket,REALK)*up_hbra_pints(3,addr_low_zket,addr_cur_geo)) 1229 if (max_cur_hbra>1) then 1230 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 1231 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 1232 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 1233 ! loops over other current order of HGTOs on bra center, starting from d-shell 1234 do order_hbra = 2, max_cur_hbra 1235 addr_up_hbra = addr_up_hbra+1 1236 addr_cur_hbra = addr_cur_hbra+1 1237#if defined(DEBUG) 1238 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 1239 write(STDOUT,100) "x-direction" 1240 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1241 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1242 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 1243 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1244 write(STDOUT,110) addr_cur_hbra, addr_low_xket, addr_up_geo 1245 write(STDOUT,100) "------------------------------------" 1246#endif 1247 ! x-direction 1248 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1249 = cc_wrt_bra(1) & 1250 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1251 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 1252 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 1253 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo) & 1254 - real(order_xket,REALK) & 1255 * up_hbra_pints(addr_cur_hbra,addr_low_xket,addr_cur_geo)) 1256 ! y-direction 1257 addr_up_hbra = addr_up_hbra+1 1258#if defined(DEBUG) 1259 write(STDOUT,100) "y-direction" 1260 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1261 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1262 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 1263 write(STDOUT,100) "------------------------------------" 1264#endif 1265 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1266 = cc_wrt_bra(2) & 1267 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1268 + half_neg_rp & 1269 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1)) 1270 do ibra = 1, order_hbra 1271 addr_up_hbra = addr_up_hbra+1 1272#if defined(DEBUG) 1273 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1274 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 1275 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 1276 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 1277 write(STDOUT,100) "------------------------------------" 1278#endif 1279 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1280 = cc_wrt_bra(2) & 1281 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1282 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1283 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1284 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1)) 1285 end do 1286 ! z-direction 1287#if defined(DEBUG) 1288 write(STDOUT,100) "z-direction" 1289#endif 1290 addr_cur_hbra = addr_cur_hbra-1 1291 do jbra = 0, order_hbra 1292 addr_up_hbra = addr_up_hbra+1 1293 addr_cur_hbra = addr_cur_hbra+1 1294#if defined(DEBUG) 1295 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1296 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1297 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1298 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1299 write(STDOUT,100) "------------------------------------" 1300#endif 1301 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1302 = cc_wrt_bra(3) & 1303 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1304 + half_neg_rp & 1305 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) & 1306 - real(order_zket,REALK) & 1307 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1308 end do 1309 do ibra = 1, order_hbra 1310 do jbra = 0, order_hbra-ibra 1311 addr_up_hbra = addr_up_hbra+1 1312 addr_cur_hbra = addr_cur_hbra+1 1313 addr_low_hbra = addr_low_hbra+1 1314#if defined(DEBUG) 1315 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1316 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1317 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 1318 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1319 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1320 write(STDOUT,100) "------------------------------------" 1321#endif 1322 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1323 = cc_wrt_bra(3) & 1324 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1325 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1326 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1327 + low_hbra_pints(addr_cur_hbra,addr_hket_low, & 1328 addr_up_geo+igeo+2) & 1329 - real(order_zket,REALK) & 1330 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1331 end do 1332 end do 1333 end do 1334 end if 1335 end if 1336 ! (4.2) x...xyz...z to xy...yz...z components of upper order HGTOs on ket center 1337 do order_yket = 1, order_hket-(order_zket+1) 1338 addr_cur_hket = addr_cur_hket+1 1339 addr_hket_zero = addr_hket_zero+1 1340 addr_xket_zero = addr_xket_zero+1 1341 addr_yket_zero = addr_yket_zero+1 1342 addr_zket_zero = addr_zket_zero+1 1343 order_xket = order_xket-1 1344 ! px HGTO on bra center 1345 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1346 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1347 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo) & 1348 - real(order_xket,REALK)*low_geo_zero(addr_xket_zero,addr_cur_geo)) 1349 ! py HGTO on bra center 1350 addr_up_geo_y = addr_up_geo+1 1351 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1352 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1353 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_y) & 1354 - real(order_yket,REALK)*low_geo_zero(addr_yket_zero,addr_cur_geo)) 1355 ! pz HGTO on bra center 1356 addr_up_geo_z = addr_up_geo+igeo+2 1357 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1358 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1359 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_z) & 1360 - real(order_zket,REALK)*low_geo_zero(addr_zket_zero,addr_cur_geo)) 1361 ! d-shell on bra center 1362 if (max_cur_hbra>0) then 1363 addr_low_xket = addr_low_xket+1 1364 addr_low_yket = addr_low_yket+1 1365 addr_low_zket = addr_low_zket+1 1366 addr_hket_low = addr_hket_low+1 1367 ! dxx 1368 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 1369 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1370 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1371 + low_hbra_pints(1,addr_hket_low,addr_up_geo) & 1372 - real(order_xket,REALK)*up_hbra_pints(1,addr_low_xket,addr_cur_geo)) 1373 ! dxy 1374 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 1375 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1376 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_y) & 1377 - real(order_yket,REALK)*up_hbra_pints(1,addr_low_yket,addr_cur_geo)) 1378 ! dyy 1379 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 1380 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1381 + half_neg_rp*(ket_to_bra & 1382 * low_geo_zero(addr_hket_zero,addr_cur_geo) & 1383 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y) & 1384 - real(order_yket,REALK)*up_hbra_pints(2,addr_low_yket,addr_cur_geo)) 1385 ! dxz 1386 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 1387 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1388 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_z) & 1389 - real(order_zket,REALK)*up_hbra_pints(1,addr_low_zket,addr_cur_geo)) 1390 ! dyz 1391 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 1392 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1393 + half_neg_rp*(low_hbra_pints(2,addr_hket_low,addr_up_geo_z) & 1394 - real(order_zket,REALK)*up_hbra_pints(2,addr_low_zket,addr_cur_geo)) 1395 ! dzz 1396 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 1397 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1398 + half_neg_rp*(ket_to_bra & 1399 * low_geo_zero(addr_hket_zero,addr_cur_geo) & 1400 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z) & 1401 - real(order_zket,REALK)*up_hbra_pints(3,addr_low_zket,addr_cur_geo)) 1402 if (max_cur_hbra>1) then 1403 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 1404 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 1405 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 1406 ! loops over other current order of HGTOs on bra center, starting from d-shell 1407 do order_hbra = 2, max_cur_hbra 1408 addr_up_hbra = addr_up_hbra+1 1409 addr_cur_hbra = addr_cur_hbra+1 1410#if defined(DEBUG) 1411 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 1412 write(STDOUT,100) "x-direction" 1413 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1414 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1415 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 1416 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1417 write(STDOUT,110) addr_cur_hbra, addr_low_xket, addr_up_geo 1418 write(STDOUT,100) "------------------------------------" 1419#endif 1420 ! x-direction 1421 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1422 = cc_wrt_bra(1) & 1423 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1424 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 1425 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 1426 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo) & 1427 - real(order_xket,REALK) & 1428 * up_hbra_pints(addr_cur_hbra,addr_low_xket,addr_cur_geo)) 1429 ! y-direction 1430 addr_up_hbra = addr_up_hbra+1 1431#if defined(DEBUG) 1432 write(STDOUT,100) "y-direction" 1433 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1434 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1435 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 1436 write(STDOUT,110) addr_cur_hbra, addr_low_yket, addr_cur_geo 1437 write(STDOUT,100) "------------------------------------" 1438#endif 1439 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1440 = cc_wrt_bra(2) & 1441 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1442 + half_neg_rp & 1443 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) & 1444 - real(order_yket,REALK) & 1445 * up_hbra_pints(addr_cur_hbra,addr_low_yket,addr_cur_geo)) 1446 do ibra = 1, order_hbra 1447 addr_up_hbra = addr_up_hbra+1 1448#if defined(DEBUG) 1449 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1450 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 1451 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 1452 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 1453 write(STDOUT,110) addr_cur_hbra+ibra, addr_low_yket, addr_cur_geo 1454 write(STDOUT,100) "------------------------------------" 1455#endif 1456 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1457 = cc_wrt_bra(2) & 1458 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket, & 1459 addr_cur_geo) & 1460 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1461 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket, & 1462 addr_cur_geo) & 1463 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low, & 1464 addr_up_geo+1) & 1465 - real(order_yket,REALK) & 1466 * up_hbra_pints(addr_cur_hbra+ibra,addr_low_yket,addr_cur_geo)) 1467 end do 1468 ! z-direction 1469#if defined(DEBUG) 1470 write(STDOUT,100) "z-direction" 1471#endif 1472 addr_cur_hbra = addr_cur_hbra-1 1473 do jbra = 0, order_hbra 1474 addr_up_hbra = addr_up_hbra+1 1475 addr_cur_hbra = addr_cur_hbra+1 1476#if defined(DEBUG) 1477 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1478 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1479 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1480 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1481 write(STDOUT,100) "------------------------------------" 1482#endif 1483 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1484 = cc_wrt_bra(3) & 1485 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1486 + half_neg_rp*(low_hbra_pints(addr_cur_hbra,addr_hket_low, & 1487 addr_up_geo+igeo+2) & 1488 - real(order_zket,REALK) & 1489 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1490 end do 1491 do ibra = 1, order_hbra 1492 do jbra = 0, order_hbra-ibra 1493 addr_up_hbra = addr_up_hbra+1 1494 addr_cur_hbra = addr_cur_hbra+1 1495 addr_low_hbra = addr_low_hbra+1 1496#if defined(DEBUG) 1497 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1498 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1499 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 1500 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1501 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1502 write(STDOUT,100) "------------------------------------" 1503#endif 1504 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1505 = cc_wrt_bra(3) & 1506 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1507 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1508 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1509 + low_hbra_pints(addr_cur_hbra,addr_hket_low, & 1510 addr_up_geo+igeo+2) & 1511 - real(order_zket,REALK) & 1512 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1513 end do 1514 end do 1515 end do 1516 end if 1517 end if 1518 end do 1519 ! (4.3) y...yz...z component of upper order HGTOs on ket center 1520 addr_cur_hket = addr_cur_hket+1 1521 addr_hket_zero = addr_hket_zero+1 1522 addr_yket_zero = addr_yket_zero+1 1523 addr_zket_zero = addr_zket_zero+1 1524 order_yket = order_hket-order_zket 1525 ! px HGTO on bra center 1526 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1527 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1528 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo) 1529 ! py HGTO on bra center 1530 addr_up_geo_y = addr_up_geo+1 1531 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1532 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1533 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_y) & 1534 - real(order_yket,REALK)*low_geo_zero(addr_yket_zero,addr_cur_geo)) 1535 ! pz HGTO on bra center 1536 addr_up_geo_z = addr_up_geo+igeo+2 1537 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1538 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1539 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_z) & 1540 - real(order_zket,REALK)*low_geo_zero(addr_zket_zero,addr_cur_geo)) 1541 ! d-shell on bra center 1542 if (max_cur_hbra>0) then 1543 addr_low_yket = addr_low_yket+1 1544 addr_low_zket = addr_low_zket+1 1545 addr_hket_low = addr_hket_low+1 1546 ! dxx 1547 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 1548 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1549 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1550 + low_hbra_pints(1,addr_hket_low,addr_up_geo)) 1551 ! dxy 1552 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 1553 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1554 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_y) & 1555 - real(order_yket,REALK)*up_hbra_pints(1,addr_low_yket,addr_cur_geo)) 1556 ! dyy 1557 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 1558 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1559 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1560 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y) & 1561 - real(order_yket,REALK)*up_hbra_pints(2,addr_low_yket,addr_cur_geo)) 1562 ! dxz 1563 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 1564 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1565 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_z) & 1566 - real(order_zket,REALK)*up_hbra_pints(1,addr_low_zket,addr_cur_geo)) 1567 ! dyz 1568 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 1569 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1570 + half_neg_rp*(low_hbra_pints(2,addr_hket_low,addr_up_geo_z) & 1571 - real(order_zket,REALK)*up_hbra_pints(2,addr_low_zket,addr_cur_geo)) 1572 ! dzz 1573 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 1574 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1575 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1576 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z) & 1577 - real(order_zket,REALK)*up_hbra_pints(3,addr_low_zket,addr_cur_geo)) 1578 if (max_cur_hbra>1) then 1579 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 1580 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 1581 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 1582 ! loops over other current order of HGTOs on bra center, starting from d-shell 1583 do order_hbra = 2, max_cur_hbra 1584 addr_up_hbra = addr_up_hbra+1 1585 addr_cur_hbra = addr_cur_hbra+1 1586#if defined(DEBUG) 1587 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 1588 write(STDOUT,100) "x-direction" 1589 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1590 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1591 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 1592 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1593 write(STDOUT,100) "------------------------------------" 1594#endif 1595 ! x-direction 1596 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1597 = cc_wrt_bra(1) & 1598 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1599 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 1600 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 1601 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo)) 1602 ! y-direction 1603 addr_up_hbra = addr_up_hbra+1 1604#if defined(DEBUG) 1605 write(STDOUT,100) "y-direction" 1606 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1607 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1608 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 1609 write(STDOUT,110) addr_cur_hbra, addr_low_yket, addr_cur_geo 1610 write(STDOUT,100) "------------------------------------" 1611#endif 1612 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1613 = cc_wrt_bra(2) & 1614 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1615 + half_neg_rp & 1616 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1) & 1617 - real(order_yket,REALK) & 1618 * up_hbra_pints(addr_cur_hbra,addr_low_yket,addr_cur_geo)) 1619 do ibra = 1, order_hbra 1620 addr_up_hbra = addr_up_hbra+1 1621#if defined(DEBUG) 1622 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1623 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 1624 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 1625 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 1626 write(STDOUT,110) addr_cur_hbra+ibra, addr_low_yket, addr_cur_geo 1627 write(STDOUT,100) "------------------------------------" 1628#endif 1629 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1630 = cc_wrt_bra(2) & 1631 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1632 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1633 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1634 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1) & 1635 - real(order_yket,REALK) & 1636 * up_hbra_pints(addr_cur_hbra+ibra,addr_low_yket,addr_cur_geo)) 1637 end do 1638 ! z-direction 1639#if defined(DEBUG) 1640 write(STDOUT,100) "z-direction" 1641#endif 1642 addr_cur_hbra = addr_cur_hbra-1 1643 do jbra = 0, order_hbra 1644 addr_up_hbra = addr_up_hbra+1 1645 addr_cur_hbra = addr_cur_hbra+1 1646#if defined(DEBUG) 1647 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1648 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1649 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1650 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1651 write(STDOUT,100) "------------------------------------" 1652#endif 1653 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1654 = cc_wrt_bra(3) & 1655 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1656 + half_neg_rp & 1657 * (low_hbra_pints(addr_cur_hbra,addr_hket_low, & 1658 addr_up_geo+igeo+2) & 1659 - real(order_zket,REALK) & 1660 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1661 end do 1662 do ibra = 1, order_hbra 1663 do jbra = 0, order_hbra-ibra 1664 addr_up_hbra = addr_up_hbra+1 1665 addr_cur_hbra = addr_cur_hbra+1 1666 addr_low_hbra = addr_low_hbra+1 1667#if defined(DEBUG) 1668 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1669 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1670 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 1671 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1672 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1673 write(STDOUT,100) "------------------------------------" 1674#endif 1675 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1676 = cc_wrt_bra(3) & 1677 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1678 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1679 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1680 + low_hbra_pints(addr_cur_hbra,addr_hket_low, & 1681 addr_up_geo+igeo+2) & 1682 - real(order_zket,REALK) & 1683 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1684 end do 1685 end do 1686 end do 1687 end if 1688 end if 1689 end do 1690 ! (5) z...z component of upper order HGTOs on ket center 1691 addr_cur_hket = addr_cur_hket+1 1692 addr_hket_zero = addr_hket_zero+1 1693 addr_zket_zero = addr_zket_zero+1 1694 ! px HGTO on bra center 1695 up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1696 = cc_wrt_bra(1)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1697 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo) 1698 ! py HGTO on bra center 1699 addr_up_geo_y = addr_up_geo+1 1700 up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1701 = cc_wrt_bra(2)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1702 + half_neg_rp*up_geo_zero(addr_cur_hket,addr_up_geo_y) 1703 ! pz HGTO on bra center 1704 addr_up_geo_z = addr_up_geo+igeo+2 1705 up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1706 = cc_wrt_bra(3)*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1707 + half_neg_rp*(up_geo_zero(addr_cur_hket,addr_up_geo_z) & 1708 - real(order_hket,REALK)*low_geo_zero(addr_zket_zero,addr_cur_geo)) 1709 ! d-shell on bra center 1710 if (max_cur_hbra>0) then 1711 addr_low_zket = addr_low_zket+1 1712 addr_hket_low = addr_hket_low+1 1713 ! dxx 1714 up_hbra_pints(4,addr_cur_hket,addr_cur_geo) & 1715 = cc_wrt_bra(1)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1716 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1717 + low_hbra_pints(1,addr_hket_low,addr_up_geo)) 1718 ! dxy 1719 up_hbra_pints(5,addr_cur_hket,addr_cur_geo) & 1720 = cc_wrt_bra(2)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1721 + half_neg_rp*low_hbra_pints(1,addr_hket_low,addr_up_geo_y) 1722 ! dyy 1723 up_hbra_pints(6,addr_cur_hket,addr_cur_geo) & 1724 = cc_wrt_bra(2)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1725 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1726 + low_hbra_pints(2,addr_hket_low,addr_up_geo_y)) 1727 ! dxz 1728 up_hbra_pints(7,addr_cur_hket,addr_cur_geo) & 1729 = cc_wrt_bra(3)*up_hbra_pints(1,addr_cur_hket,addr_cur_geo) & 1730 + half_neg_rp*(low_hbra_pints(1,addr_hket_low,addr_up_geo_z) & 1731 - real(order_hket,REALK)*up_hbra_pints(1,addr_low_zket,addr_cur_geo)) 1732 ! dyz 1733 up_hbra_pints(8,addr_cur_hket,addr_cur_geo) & 1734 = cc_wrt_bra(3)*up_hbra_pints(2,addr_cur_hket,addr_cur_geo) & 1735 + half_neg_rp*(low_hbra_pints(2,addr_hket_low,addr_up_geo_z) & 1736 - real(order_hket,REALK)*up_hbra_pints(2,addr_low_zket,addr_cur_geo)) 1737 ! dzz 1738 up_hbra_pints(9,addr_cur_hket,addr_cur_geo) & 1739 = cc_wrt_bra(3)*up_hbra_pints(3,addr_cur_hket,addr_cur_geo) & 1740 + half_neg_rp*(ket_to_bra*low_geo_zero(addr_hket_zero,addr_cur_geo) & 1741 + low_hbra_pints(3,addr_hket_low,addr_up_geo_z) & 1742 - real(order_hket,REALK)*up_hbra_pints(3,addr_low_zket,addr_cur_geo)) 1743 if (max_cur_hbra>1) then 1744 addr_up_hbra = 9 !base address of the f-shell HGTOs on bra center 1745 addr_cur_hbra = 3 !base address of the d-shell HGTOs on bra center 1746 addr_low_hbra = 0 !base address of the p-shell HGTOs on bra center 1747 ! loops over other current order of HGTOs on bra center, starting from d-shell 1748 do order_hbra = 2, max_cur_hbra 1749 addr_up_hbra = addr_up_hbra+1 1750 addr_cur_hbra = addr_cur_hbra+1 1751#if defined(DEBUG) 1752 write(STDOUT,100) "orders:", order_hbra, order_xket, cur_order_geo 1753 write(STDOUT,100) "x-direction" 1754 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1755 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1756 write(STDOUT,110) addr_low_hbra+1, addr_cur_hket, addr_cur_geo 1757 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1758 write(STDOUT,100) "------------------------------------" 1759#endif 1760 ! x-direction 1761 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1762 = cc_wrt_bra(1) & 1763 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1764 + half_neg_rp*(real(order_hbra,REALK)*ket_to_bra & 1765 * up_hbra_pints(addr_low_hbra+1,addr_cur_hket,addr_cur_geo) & 1766 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo)) 1767 ! y-direction 1768 addr_up_hbra = addr_up_hbra+1 1769#if defined(DEBUG) 1770 write(STDOUT,100) "y-direction" 1771 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1772 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1773 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+1 1774 write(STDOUT,100) "------------------------------------" 1775#endif 1776 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1777 = cc_wrt_bra(2) & 1778 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1779 + half_neg_rp & 1780 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+1)) 1781 do ibra = 1, order_hbra 1782 addr_up_hbra = addr_up_hbra+1 1783#if defined(DEBUG) 1784 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1785 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_cur_geo 1786 write(STDOUT,110) addr_low_hbra+ibra, addr_cur_hket, addr_cur_geo 1787 write(STDOUT,110) addr_cur_hbra+ibra, addr_cur_hket, addr_up_geo+1 1788 write(STDOUT,100) "------------------------------------" 1789#endif 1790 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1791 = cc_wrt_bra(2) & 1792 * up_hbra_pints(addr_cur_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1793 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1794 * up_hbra_pints(addr_low_hbra+ibra,addr_cur_hket,addr_cur_geo) & 1795 + low_hbra_pints(addr_cur_hbra+ibra,addr_hket_low,addr_up_geo+1)) 1796 end do 1797 ! z-direction 1798#if defined(DEBUG) 1799 write(STDOUT,100) "z-direction" 1800#endif 1801 addr_cur_hbra = addr_cur_hbra-1 1802 do jbra = 0, order_hbra 1803 addr_up_hbra = addr_up_hbra+1 1804 addr_cur_hbra = addr_cur_hbra+1 1805#if defined(DEBUG) 1806 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1807 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1808 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1809 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1810 write(STDOUT,100) "------------------------------------" 1811#endif 1812 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1813 = cc_wrt_bra(3) & 1814 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1815 + half_neg_rp & 1816 * (low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) & 1817 - real(order_hket,REALK) & 1818 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1819 end do 1820 do ibra = 1, order_hbra 1821 do jbra = 0, order_hbra-ibra 1822 addr_up_hbra = addr_up_hbra+1 1823 addr_cur_hbra = addr_cur_hbra+1 1824 addr_low_hbra = addr_low_hbra+1 1825#if defined(DEBUG) 1826 write(STDOUT,110) addr_up_hbra, addr_cur_hket, addr_cur_geo 1827 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_cur_geo 1828 write(STDOUT,110) addr_low_hbra, addr_cur_hket, addr_cur_geo 1829 write(STDOUT,110) addr_cur_hbra, addr_cur_hket, addr_up_geo+igeo+2 1830 write(STDOUT,110) addr_cur_hbra, addr_low_zket, addr_cur_geo 1831 write(STDOUT,100) "------------------------------------" 1832#endif 1833 up_hbra_pints(addr_up_hbra,addr_cur_hket,addr_cur_geo) & 1834 = cc_wrt_bra(3) & 1835 * up_hbra_pints(addr_cur_hbra,addr_cur_hket,addr_cur_geo) & 1836 + half_neg_rp*(real(ibra,REALK)*ket_to_bra & 1837 * up_hbra_pints(addr_low_hbra,addr_cur_hket,addr_cur_geo) & 1838 + low_hbra_pints(addr_cur_hbra,addr_hket_low,addr_up_geo+igeo+2) & 1839 - real(order_hket,REALK) & 1840 * up_hbra_pints(addr_cur_hbra,addr_low_zket,addr_cur_geo)) 1841 end do 1842 end do 1843 end do 1844 end if 1845 end if 1846 ! updates the base addresses of lower order HGTOs on ket center 1847 base_low_hket = addr_low_zket 1848 base_hket_zero = addr_zket_zero 1849 end do 1850 end do 1851 addr_up_geo = addr_up_geo+1 1852 end do 1853#if defined(XTIME) 1854 ! prints the CPU elapsed time 1855 call xtimer_view(curr_time, "sub_nucpot_hbra", STDOUT) 1856#endif 1857 return 1858#if defined(DEBUG) 1859100 format("sub_nucpot_hbra>> ",A,3I6) 1860110 format("sub_nucpot_hbra>> ","HBRA",I8,4X,"HKET",I8,4X,"DERV",I8) 1861#endif 1862 end subroutine sub_nucpot_hbra 1863 1864 !> \brief assigns the integrals with required HGTOs on bra center 1865 !> \author Bin Gao 1866 !> \date 2012-03-04 1867 !> \param offset_hgto_ket is the offset of HGTOs on ket center 1868 !> \param offset_hket_geo is the offset of geometric derivatives on nuclear potential origin 1869 !> \param dim_hket is the dimension of HGTOs on ket center 1870 !> \param dim_geo_hket is the dimension of geometric derivatives on nuclear potential origin 1871 !> \param hket_pints contains the nuclear attraction integrals with zeroth order 1872 !> HGTOs on bra center 1873 !> \param dim_up_hbra is the dimension of upper order HGTOs on bra center in temporary integrals 1874 !> \param dim_cur_hket is the dimension of current order HGTOs on ket center in temporary integrals 1875 !> \param num_low_geo is the number of lower order geometric derivatives in temporary integrals 1876 !> \param recur_pints contains the temporary integrals from recurrence relations 1877 !> \param zero_hbra indicates if zeroth order HGTO returned 1878 !> \param offset_hbra_geo is the offset of geometric derivatives on nuclear potential origin 1879 !> in returned integrals 1880 !> \param dim_hgto_bra is the dimension of HGTOs of bra center afterwards 1881 !> \param dim_hgto_ket is the dimension of HGTOs of ket center afterwards 1882 !> \param dim_geo_hbra is the dimension of geometric derivatives afterwards 1883 !> \return hbra_pints contains the integrals with required HGTOs on bra center 1884 subroutine nucpot_hbra_assign(offset_hgto_ket, offset_hket_geo, dim_hket, & 1885 dim_geo_hket, hket_pints, dim_up_hbra, & 1886 dim_cur_hket, num_low_geo, recur_pints, & 1887 zero_hbra, offset_hbra_geo, dim_hgto_bra, & 1888 dim_hgto_ket, dim_geo_hbra, hbra_pints) 1889 use xkind 1890 implicit none 1891 integer, intent(in) :: offset_hgto_ket 1892 integer, intent(in) :: offset_hket_geo 1893 integer, intent(in) :: dim_hket 1894 integer, intent(in) :: dim_geo_hket 1895 real(REALK), intent(in) :: hket_pints(dim_hket,dim_geo_hket) 1896 integer, intent(in) :: dim_up_hbra 1897 integer, intent(in) :: dim_cur_hket 1898 integer, intent(in) :: num_low_geo 1899 real(REALK), intent(in) :: recur_pints(dim_up_hbra,dim_cur_hket,num_low_geo) 1900 logical, intent(in) :: zero_hbra 1901 integer, intent(in) :: offset_hbra_geo 1902 integer, intent(in) :: dim_hgto_bra 1903 integer, intent(in) :: dim_hgto_ket 1904 integer, intent(in) :: dim_geo_hbra 1905 real(REALK), intent(inout) :: hbra_pints(dim_hgto_bra,dim_hgto_ket,dim_geo_hbra) 1906!f2py intent(in) :: offset_hgto_ket 1907!f2py intent(in) :: offset_hket_geo 1908!f2py intent(hide) :: dim_hket 1909!f2py intent(hide) :: dim_geo_hket 1910!f2py intent(in) :: hket_pints 1911!f2py intent(hide) :: dim_up_hbra 1912!f2py intent(hide) :: dim_cur_hket 1913!f2py intent(hide) :: num_low_geo 1914!f2py intent(in) :: recur_pints 1915!f2py intent(in) :: zero_hbra 1916!f2py intent(in) :: offset_hbra_geo 1917!f2py intent(hide) :: dim_hgto_bra 1918!f2py intent(hide) :: dim_hgto_ket 1919!f2py intent(hide) :: dim_geo_hbra 1920!f2py intent(inout) :: hbra_pints 1921 integer offset_cur_hket !offset of current order HGTOs on ket center in temporary integrals 1922 integer start_up_hbra !start address of upper order HGTOs on bra center in temporary integrals 1923 integer addr_hket_geo !address of geometric derivatives on nuclear potential origin 1924 integer addr_hbra_geo !address of geometric derivatives on nuclear potential origin in returned integrals 1925 integer igeo !incremental recorder over geometric derivatives 1926 integer iket !incremental recorder over HGTOs on ket center 1927#if defined(XTIME) 1928 real(REALK) curr_time !current CPU time 1929 ! sets current CPU time 1930 call xtimer_set(curr_time) 1931#endif 1932 ! sets the offset of current order HGTOs on ket center in temporary integrals 1933 offset_cur_hket = dim_cur_hket-dim_hgto_ket 1934 ! s-shell HGTO returned 1935 if (zero_hbra) then 1936 start_up_hbra = dim_up_hbra-dim_hgto_bra+2 1937 do igeo = 1, num_low_geo 1938 addr_hket_geo = offset_hket_geo+igeo 1939 addr_hbra_geo = offset_hbra_geo+igeo 1940 do iket = 1, dim_hgto_ket 1941 hbra_pints(1,iket,addr_hbra_geo) & 1942 = hket_pints(offset_hgto_ket+iket,addr_hket_geo) 1943 hbra_pints(2:dim_hgto_bra,iket,addr_hbra_geo) & 1944 = recur_pints(start_up_hbra:dim_up_hbra,offset_cur_hket+iket,igeo) 1945 end do 1946 end do 1947 else 1948 start_up_hbra = dim_up_hbra-dim_hgto_bra+1 1949 do igeo = 1, num_low_geo 1950 addr_hbra_geo = offset_hbra_geo+igeo 1951 do iket = 1, dim_hgto_ket 1952 hbra_pints(:,iket,addr_hbra_geo) & 1953 = recur_pints(start_up_hbra:dim_up_hbra,offset_cur_hket+iket,igeo) 1954 end do 1955 end do 1956 end if 1957#if defined(XTIME) 1958 ! prints the CPU elapsed time 1959 call xtimer_view(curr_time, "nucpot_hbra_assign", STDOUT) 1960#endif 1961 return 1962 end subroutine nucpot_hbra_assign 1963