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 returns the zeroth order Cartesian multipole moment integrals 18!! with specific orders of HGTOs on bra center. 19!! 20!! 2012-02-12, Bin Gao: 21!! * first version 22 23#include "stdout.h" 24 25 !> \brief computes the zeroth order Cartesian multipole moment integrals 26 !> with specific orders of HGTOs on bra center 27 !> \author Bin Gao 28 !> \date 2012-02-12 29 !> \param orders_hgto_bra contains the minimum and maximum orders of HGTOs on bra center 30 !> \param coord_bra contains the coordinates of bra center 31 !> \param exponent_bra is the exponent of primitive HGTO of bra center 32 !> \param coord_ket contains the coordinates of ket center 33 !> \param exponent_ket is the exponent of primitive HGTO of ket center 34 !> \param order_elec is the order of electronic derivatives 35 !> \param scal_const is the scale constant for Cartesian multipole moments 36 !> \param dim_hgto_bra is the dimension of integrals 37 !> \return hbra_pints contains the primitive HGTO integrals 38 subroutine carmom_hbra(orders_hgto_bra, coord_bra, exponent_bra, & 39 coord_ket, exponent_ket, order_elec, & 40 scal_const, dim_hgto_bra, hbra_pints) 41 use xkind 42 implicit none 43 integer, intent(in) :: orders_hgto_bra(2) 44 real(REALK), intent(in) :: coord_bra(3) 45 real(REALK), intent(in) :: exponent_bra 46 real(REALK), intent(in) :: coord_ket(3) 47 real(REALK), intent(in) :: exponent_ket 48 integer, intent(in) :: order_elec 49 real(REALK), intent(in) :: scal_const 50 integer, intent(in) :: dim_hgto_bra 51 real(REALK), intent(out) :: hbra_pints(dim_hgto_bra) 52!f2py intent(in) :: orders_hgto_bra 53!f2py intent(in) :: coord_bra 54!f2py intent(in) :: exponent_bra 55!f2py intent(in) :: coord_ket 56!f2py intent(in) :: exponent_ket 57!f2py intent(in) :: order_elec 58!f2py intent(in) :: scal_const 59!f2py intent(in) :: dim_hgto_bra 60!f2py intent(out) :: hbra_pints 61!f2py depend(dim_hgto_bra) :: hbra_pints 62#include "private/pi.h" 63 real(REALK) recip_total_expnt !reciprocal of total exponent 64 real(REALK) reduced_expnt !reduced exponent 65 real(REALK) sd_bra_ket !square of the relative distance between bra and ket centers 66 real(REALK) cc_wrt_bra(3) !relative coordinates of center-of-charge w.r.t. bra center 67 real(REALK) half_recip_nup !\f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 68 integer base_low_hgto !base address of lower order HGTOs 69 integer base_cur_hgto !base address of current order HGTOs 70 integer base_up_hgto !base address of upper order HGTOs 71 integer order_hgto, iorder, jorder !incremental recorders over orders of HGTOs 72 integer addr_low_hgto !address of lower order HGTOs 73 integer addr_cur_hgto !address of current order HGTOs 74 integer addr_up_hgto !address of upper order HGTOs 75 integer num_min_hgto !number of minimum returned order HGTOs 76 integer dim_tmp !dimension of temporary integrals 77 real(REALK), allocatable :: tmp_ints(:) !temporary integrals 78 real(REALK) zero_pint !integral with zeroth order HGTO 79 integer ierr !error information 80#if defined(XTIME) 81 real(REALK) curr_time !current CPU time 82 ! sets current CPU time 83 call xtimer_set(curr_time) 84#endif 85 ! computes the reciprocal of total exponent 86 recip_total_expnt = 1.0_REALK/(exponent_bra+exponent_ket) 87 ! computes the reduced exponent 88 reduced_expnt = exponent_bra*exponent_ket*recip_total_expnt 89 ! computes the square of the relative distance between bra and ket centers 90 sd_bra_ket = 0.0_REALK 91 do iorder = 1, 3 92 sd_bra_ket = sd_bra_ket+(coord_ket(iorder)-coord_bra(iorder))**2 93 end do 94 select case(orders_hgto_bra(1)) 95 ! zeroth order HGTOs returned 96 case(0) 97 select case(orders_hgto_bra(2)) 98 ! only zeroth order HGTO returned 99 case(0) 100 hbra_pints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 101 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 102 ! zeroth and first order HGTOs returned 103 case(1) 104 ! zeroth order 105 hbra_pints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 106 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 107 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 108 half_recip_nup = -exponent_ket*recip_total_expnt 109 ! computes the relative coordinates of center-of-charge w.r.t. bra center 110 do iorder = 1, 3 111 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 112 end do 113 ! first order 114 hbra_pints(2) = cc_wrt_bra(1)*hbra_pints(1) 115 hbra_pints(3) = cc_wrt_bra(2)*hbra_pints(1) 116 hbra_pints(4) = cc_wrt_bra(3)*hbra_pints(1) 117 ! up to, at least, second order HGTOs returned 118 case default 119 ! zeroth order 120 hbra_pints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 121 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 122 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 123 half_recip_nup = -exponent_ket*recip_total_expnt 124 ! computes the relative coordinates of center-of-charge w.r.t. bra center 125 do iorder = 1, 3 126 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 127 end do 128 ! first order 129 hbra_pints(2) = cc_wrt_bra(1)*hbra_pints(1) 130 hbra_pints(3) = cc_wrt_bra(2)*hbra_pints(1) 131 hbra_pints(4) = cc_wrt_bra(3)*hbra_pints(1) 132 ! computes parameter \f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 133 half_recip_nup = 0.5_REALK*half_recip_nup/exponent_bra 134 ! initializes base addresses 135 base_low_hgto = 0 136 base_cur_hgto = 1 137 base_up_hgto = 4 138 ! left returned orders up to \var(orders_hgto_bra(2)) 139 do order_hgto = 1, orders_hgto_bra(2)-1 140 ! recurrence relations along x-direction 141 addr_up_hgto = base_up_hgto+1 142 addr_cur_hgto = base_cur_hgto+1 143 hbra_pints(addr_up_hgto) & 144 = cc_wrt_bra(1)*hbra_pints(addr_cur_hgto) & 145 + real(order_hgto,REALK)*half_recip_nup*hbra_pints(base_low_hgto+1) 146 ! recurrence relations along y-direction 147 addr_up_hgto = addr_up_hgto+1 148 hbra_pints(addr_up_hgto) = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) 149 do iorder = 1, order_hgto 150 addr_up_hgto = addr_up_hgto+1 151 addr_cur_hgto = addr_cur_hgto+1 152 hbra_pints(addr_up_hgto) & 153 = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) & 154 + real(iorder,REALK)*half_recip_nup*hbra_pints(base_low_hgto+iorder) 155 end do 156 ! recurrence relations along z-direction 157 addr_cur_hgto = base_cur_hgto 158 do jorder = 0, order_hgto 159 addr_up_hgto = addr_up_hgto+1 160 addr_cur_hgto = addr_cur_hgto+1 161 hbra_pints(addr_up_hgto) = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) 162 end do 163 addr_low_hgto = base_low_hgto 164 do iorder = 1, order_hgto 165 do jorder = 0, order_hgto-iorder 166 addr_up_hgto = addr_up_hgto+1 167 addr_cur_hgto = addr_cur_hgto+1 168 addr_low_hgto = addr_low_hgto+1 169 hbra_pints(addr_up_hgto) & 170 = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) & 171 + real(iorder,REALK)*half_recip_nup*hbra_pints(addr_low_hgto) 172 end do 173 end do 174 ! updates base addresses 175 base_low_hgto = base_cur_hgto 176 base_cur_hgto = base_up_hgto 177 base_up_hgto = addr_up_hgto 178 end do 179 end select 180 ! first order HGTOs returned 181 case (1) 182 select case(orders_hgto_bra(2)) 183 ! only first order HGTOs returned 184 case(1) 185 ! zeroth order 186 hbra_pints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 187 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 188 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 189 half_recip_nup = -exponent_ket*recip_total_expnt 190 ! computes the relative coordinates of center-of-charge w.r.t. bra center 191 do iorder = 1, 3 192 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 193 end do 194 ! first order 195 hbra_pints(2) = cc_wrt_bra(2)*hbra_pints(1) 196 hbra_pints(3) = cc_wrt_bra(3)*hbra_pints(1) 197 hbra_pints(1) = cc_wrt_bra(1)*hbra_pints(1) 198 ! first and second orders HGTOs returned 199 case(2) 200 ! computes the integral with zeroth order HGTO 201 zero_pint = scal_const*(-exponent_ket-exponent_ket)**order_elec & 202 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 203 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 204 half_recip_nup = -exponent_ket*recip_total_expnt 205 ! computes the relative coordinates of center-of-charge w.r.t. bra center 206 do iorder = 1, 3 207 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 208 end do 209 ! computes parameter \f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 210 half_recip_nup = 0.5_REALK*half_recip_nup/exponent_bra 211 ! first order 212 hbra_pints(1) = cc_wrt_bra(1)*zero_pint 213 hbra_pints(2) = cc_wrt_bra(2)*zero_pint 214 hbra_pints(3) = cc_wrt_bra(3)*zero_pint 215 ! second order 216 hbra_pints(4) = cc_wrt_bra(1)*hbra_pints(1)+half_recip_nup*zero_pint 217 hbra_pints(5) = cc_wrt_bra(2)*hbra_pints(1) 218 hbra_pints(6) = cc_wrt_bra(2)*hbra_pints(2)+half_recip_nup*zero_pint 219 hbra_pints(7) = cc_wrt_bra(3)*hbra_pints(1) 220 hbra_pints(8) = cc_wrt_bra(3)*hbra_pints(2) 221 hbra_pints(9) = cc_wrt_bra(3)*hbra_pints(3)+half_recip_nup*zero_pint 222 ! up to, at least, third order HGTOs returned 223 case default 224 ! computes the integral with zeroth order HGTO 225 zero_pint = scal_const*(-exponent_ket-exponent_ket)**order_elec & 226 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 227 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 228 half_recip_nup = -exponent_ket*recip_total_expnt 229 ! computes the relative coordinates of center-of-charge w.r.t. bra center 230 do iorder = 1, 3 231 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 232 end do 233 ! computes parameter \f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 234 half_recip_nup = 0.5_REALK*half_recip_nup/exponent_bra 235 ! first order 236 hbra_pints(1) = cc_wrt_bra(1)*zero_pint 237 hbra_pints(2) = cc_wrt_bra(2)*zero_pint 238 hbra_pints(3) = cc_wrt_bra(3)*zero_pint 239 ! second order 240 hbra_pints(4) = cc_wrt_bra(1)*hbra_pints(1)+half_recip_nup*zero_pint 241 hbra_pints(5) = cc_wrt_bra(2)*hbra_pints(1) 242 hbra_pints(6) = cc_wrt_bra(2)*hbra_pints(2)+half_recip_nup*zero_pint 243 hbra_pints(7) = cc_wrt_bra(3)*hbra_pints(1) 244 hbra_pints(8) = cc_wrt_bra(3)*hbra_pints(2) 245 hbra_pints(9) = cc_wrt_bra(3)*hbra_pints(3)+half_recip_nup*zero_pint 246 ! initializes base addresses 247 base_low_hgto = 0 248 base_cur_hgto = 3 249 base_up_hgto = 9 250 ! left returned orders up to \var(orders_hgto_bra(2)) 251 do order_hgto = 2, orders_hgto_bra(2)-1 252 ! recurrence relations along x-direction 253 addr_up_hgto = base_up_hgto+1 254 addr_cur_hgto = base_cur_hgto+1 255 hbra_pints(addr_up_hgto) & 256 = cc_wrt_bra(1)*hbra_pints(addr_cur_hgto) & 257 + real(order_hgto,REALK)*half_recip_nup*hbra_pints(base_low_hgto+1) 258 ! recurrence relations along y-direction 259 addr_up_hgto = addr_up_hgto+1 260 hbra_pints(addr_up_hgto) = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) 261 do iorder = 1, order_hgto 262 addr_up_hgto = addr_up_hgto+1 263 addr_cur_hgto = addr_cur_hgto+1 264 hbra_pints(addr_up_hgto) & 265 = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) & 266 + real(iorder,REALK)*half_recip_nup*hbra_pints(base_low_hgto+iorder) 267 end do 268 ! recurrence relations along z-direction 269 addr_cur_hgto = base_cur_hgto 270 do jorder = 0, order_hgto 271 addr_up_hgto = addr_up_hgto+1 272 addr_cur_hgto = addr_cur_hgto+1 273 hbra_pints(addr_up_hgto) = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) 274 end do 275 addr_low_hgto = base_low_hgto 276 do iorder = 1, order_hgto 277 do jorder = 0, order_hgto-iorder 278 addr_up_hgto = addr_up_hgto+1 279 addr_cur_hgto = addr_cur_hgto+1 280 addr_low_hgto = addr_low_hgto+1 281 hbra_pints(addr_up_hgto) & 282 = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) & 283 + real(iorder,REALK)*half_recip_nup*hbra_pints(addr_low_hgto) 284 end do 285 end do 286 ! updates base addresses 287 base_low_hgto = base_cur_hgto 288 base_cur_hgto = base_up_hgto 289 base_up_hgto = addr_up_hgto 290 end do 291 end select 292 ! high order HGTOs (>1) returned 293 case default 294 ! only order \var(orders_hgto_bra(1)) HGTOs returned 295 if (orders_hgto_bra(1)==orders_hgto_bra(2)) then 296 dim_tmp = (orders_hgto_bra(1)+1)*(orders_hgto_bra(1)+2)/2 297 allocate(tmp_ints(3*dim_tmp), stat=ierr) 298 if (ierr/=0) & 299 call error_stop("carmom_hbra", & 300 "failed to allocate tmp_ints", & 301 3*dim_tmp) 302 ! computes the integral with zeroth order HGTO 303 tmp_ints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 304 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 305 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 306 half_recip_nup = -exponent_ket*recip_total_expnt 307 ! computes the relative coordinates of center-of-charge w.r.t. bra center 308 do iorder = 1, 3 309 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 310 end do 311 ! computes parameter \f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 312 half_recip_nup = 0.5_REALK*half_recip_nup/exponent_bra 313 ! first order 314 tmp_ints(dim_tmp+1) = cc_wrt_bra(1)*tmp_ints(1) 315 tmp_ints(dim_tmp+2) = cc_wrt_bra(2)*tmp_ints(1) 316 tmp_ints(dim_tmp+3) = cc_wrt_bra(3)*tmp_ints(1) 317 ! initializes base addresses 318 base_low_hgto = 0 319 base_cur_hgto = dim_tmp 320 base_up_hgto = 2*dim_tmp 321 do order_hgto = 1, orders_hgto_bra(1)-1 322 ! recurrence relations along x-direction 323 addr_up_hgto = base_up_hgto+1 324 addr_cur_hgto = base_cur_hgto+1 325 tmp_ints(addr_up_hgto) & 326 = cc_wrt_bra(1)*tmp_ints(addr_cur_hgto) & 327 + real(order_hgto,REALK)*half_recip_nup*tmp_ints(base_low_hgto+1) 328 ! recurrence relations along y-direction 329 addr_up_hgto = addr_up_hgto+1 330 tmp_ints(addr_up_hgto) = cc_wrt_bra(2)*tmp_ints(addr_cur_hgto) 331 do iorder = 1, order_hgto 332 addr_up_hgto = addr_up_hgto+1 333 addr_cur_hgto = addr_cur_hgto+1 334 tmp_ints(addr_up_hgto) & 335 = cc_wrt_bra(2)*tmp_ints(addr_cur_hgto) & 336 + real(iorder,REALK)*half_recip_nup*tmp_ints(base_low_hgto+iorder) 337 end do 338 ! recurrence relations along z-direction 339 addr_cur_hgto = base_cur_hgto 340 do jorder = 0, order_hgto 341 addr_up_hgto = addr_up_hgto+1 342 addr_cur_hgto = addr_cur_hgto+1 343 tmp_ints(addr_up_hgto) = cc_wrt_bra(3)*tmp_ints(addr_cur_hgto) 344 end do 345 addr_low_hgto = base_low_hgto 346 do iorder = 1, order_hgto 347 do jorder = 0, order_hgto-iorder 348 addr_up_hgto = addr_up_hgto+1 349 addr_cur_hgto = addr_cur_hgto+1 350 addr_low_hgto = addr_low_hgto+1 351 tmp_ints(addr_up_hgto) & 352 = cc_wrt_bra(3)*tmp_ints(addr_cur_hgto) & 353 + real(iorder,REALK)*half_recip_nup*tmp_ints(addr_low_hgto) 354 end do 355 end do 356 ! updates base addresses 357 addr_up_hgto = base_low_hgto 358 base_low_hgto = base_cur_hgto 359 base_cur_hgto = base_up_hgto 360 base_up_hgto = addr_up_hgto 361 end do 362 ! assigns order \var(orders_hgto_bra(1)) 363 hbra_pints(:) = tmp_ints(base_cur_hgto+1:base_cur_hgto+dim_tmp) 364 deallocate(tmp_ints) 365 else 366 num_min_hgto = (orders_hgto_bra(1)+1)*(orders_hgto_bra(1)+2)/2 367 dim_tmp = num_min_hgto+orders_hgto_bra(1)+2 368 allocate(tmp_ints(3*dim_tmp), stat=ierr) 369 if (ierr/=0) & 370 call error_stop("carmom_hbra", & 371 "failed to allocate tmp_ints", & 372 3*dim_tmp) 373 ! computes the integral with zeroth order HGTO 374 tmp_ints(1) = scal_const*(-exponent_ket-exponent_ket)**order_elec & 375 * exp(-reduced_expnt*sd_bra_ket)*sqrt((PI*recip_total_expnt)**3) 376 ! computes parameter \f$-\frac{b_{j\lambda}}{p_{ij}}\f$ 377 half_recip_nup = -exponent_ket*recip_total_expnt 378 ! computes the relative coordinates of center-of-charge w.r.t. bra center 379 do iorder = 1, 3 380 cc_wrt_bra(iorder) = half_recip_nup*(coord_bra(iorder)-coord_ket(iorder)) 381 end do 382 ! computes parameter \f$-\frac{1}{2p_{ij}}\frac{b_{j\lambda}}{a_{i\kappa}}\f$ 383 half_recip_nup = 0.5_REALK*half_recip_nup/exponent_bra 384 ! first order 385 tmp_ints(dim_tmp+1) = cc_wrt_bra(1)*tmp_ints(1) 386 tmp_ints(dim_tmp+2) = cc_wrt_bra(2)*tmp_ints(1) 387 tmp_ints(dim_tmp+3) = cc_wrt_bra(3)*tmp_ints(1) 388 ! initializes base addresses 389 base_low_hgto = 0 390 base_cur_hgto = dim_tmp 391 base_up_hgto = 2*dim_tmp 392 ! orders up to \var(orders_hgto_bra(1))+1 393 do order_hgto = 1, orders_hgto_bra(1) 394 ! recurrence relations along x-direction 395 addr_up_hgto = base_up_hgto+1 396 addr_cur_hgto = base_cur_hgto+1 397 tmp_ints(addr_up_hgto) & 398 = cc_wrt_bra(1)*tmp_ints(addr_cur_hgto) & 399 + real(order_hgto,REALK)*half_recip_nup*tmp_ints(base_low_hgto+1) 400 ! recurrence relations along y-direction 401 addr_up_hgto = addr_up_hgto+1 402 tmp_ints(addr_up_hgto) = cc_wrt_bra(2)*tmp_ints(addr_cur_hgto) 403 do iorder = 1, order_hgto 404 addr_up_hgto = addr_up_hgto+1 405 addr_cur_hgto = addr_cur_hgto+1 406 tmp_ints(addr_up_hgto) & 407 = cc_wrt_bra(2)*tmp_ints(addr_cur_hgto) & 408 + real(iorder,REALK)*half_recip_nup*tmp_ints(base_low_hgto+iorder) 409 end do 410 ! recurrence relations along z-direction 411 addr_cur_hgto = base_cur_hgto 412 do jorder = 0, order_hgto 413 addr_up_hgto = addr_up_hgto+1 414 addr_cur_hgto = addr_cur_hgto+1 415 tmp_ints(addr_up_hgto) = cc_wrt_bra(3)*tmp_ints(addr_cur_hgto) 416 end do 417 addr_low_hgto = base_low_hgto 418 do iorder = 1, order_hgto 419 do jorder = 0, order_hgto-iorder 420 addr_up_hgto = addr_up_hgto+1 421 addr_cur_hgto = addr_cur_hgto+1 422 addr_low_hgto = addr_low_hgto+1 423 tmp_ints(addr_up_hgto) & 424 = cc_wrt_bra(3)*tmp_ints(addr_cur_hgto) & 425 + real(iorder,REALK)*half_recip_nup*tmp_ints(addr_low_hgto) 426 end do 427 end do 428 ! updates base addresses 429 addr_up_hgto = base_low_hgto 430 base_low_hgto = base_cur_hgto 431 base_cur_hgto = base_up_hgto 432 base_up_hgto = addr_up_hgto 433 end do 434 ! assigns orders \var(orders_hgto_bra(1)) and \var(orders_hgto_bra(1))+1, 435 ! and initializes base addresses 436 hbra_pints(1:num_min_hgto) & 437 = tmp_ints(base_low_hgto+1:base_low_hgto+num_min_hgto) 438 base_up_hgto = num_min_hgto+dim_tmp 439 hbra_pints(num_min_hgto+1:base_up_hgto) & 440 = tmp_ints(base_cur_hgto+1:base_cur_hgto+dim_tmp) 441 deallocate(tmp_ints) 442 base_low_hgto = 0 443 base_cur_hgto = num_min_hgto 444 ! left returned orders up to \var(orders_hgto_bra(2)) 445 do order_hgto = orders_hgto_bra(1)+1, orders_hgto_bra(2)-1 446 ! recurrence relations along x-direction 447 addr_up_hgto = base_up_hgto+1 448 addr_cur_hgto = base_cur_hgto+1 449 hbra_pints(addr_up_hgto) & 450 = cc_wrt_bra(1)*hbra_pints(addr_cur_hgto) & 451 + real(order_hgto,REALK)*half_recip_nup*hbra_pints(base_low_hgto+1) 452 ! recurrence relations along y-direction 453 addr_up_hgto = addr_up_hgto+1 454 hbra_pints(addr_up_hgto) = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) 455 do iorder = 1, order_hgto 456 addr_up_hgto = addr_up_hgto+1 457 addr_cur_hgto = addr_cur_hgto+1 458 hbra_pints(addr_up_hgto) & 459 = cc_wrt_bra(2)*hbra_pints(addr_cur_hgto) & 460 + real(iorder,REALK)*half_recip_nup*hbra_pints(base_low_hgto+iorder) 461 end do 462 ! recurrence relations along z-direction 463 addr_cur_hgto = base_cur_hgto 464 do jorder = 0, order_hgto 465 addr_up_hgto = addr_up_hgto+1 466 addr_cur_hgto = addr_cur_hgto+1 467 hbra_pints(addr_up_hgto) = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) 468 end do 469 addr_low_hgto = base_low_hgto 470 do iorder = 1, order_hgto 471 do jorder = 0, order_hgto-iorder 472 addr_up_hgto = addr_up_hgto+1 473 addr_cur_hgto = addr_cur_hgto+1 474 addr_low_hgto = addr_low_hgto+1 475 hbra_pints(addr_up_hgto) & 476 = cc_wrt_bra(3)*hbra_pints(addr_cur_hgto) & 477 + real(iorder,REALK)*half_recip_nup*hbra_pints(addr_low_hgto) 478 end do 479 end do 480 ! updates base addresses 481 base_low_hgto = base_cur_hgto 482 base_cur_hgto = base_up_hgto 483 base_up_hgto = addr_up_hgto 484 end do 485 end if 486 end select 487#if defined(XTIME) 488 ! prints the CPU elapsed time 489 call xtimer_view(curr_time, "carmom_hbra", STDOUT) 490#endif 491 return 492 end subroutine carmom_hbra 493