1!! gen1int: compute derivatives of one-electron integrals using Hermite Gaussians 2!! Copyright 2009-2012 Bin Gao, and Andreas Thorvaldsen 3!! 4!! This file is part of gen1int. 5!! 6!! gen1int is free software: you can redistribute it and/or modify 7!! it under the terms of the GNU Lesser General Public License as published by 8!! the Free Software Foundation, either version 3 of the License, or 9!! (at your option) any later version. 10!! 11!! gen1int is distributed in the hope that it will be useful, 12!! but WITHOUT ANY WARRANTY; without even the implied warranty of 13!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14!! GNU Lesser General Public License for more details. 15!! 16!! You should have received a copy of the GNU Lesser General Public License 17!! along with gen1int. If not, see <http://www.gnu.org/licenses/>. 18!! 19!! This file contains subroutines for auxiliary functions. 20!! 21!! 2009-09-22, Bin Gao: 22!! * change subroutine "gen1int_aux_order" due to the modifications of types "g1prop", 23!! "g1int1" and "g1opt" 24!! 25!! 2009-08-24, Bin Gao: 26!! * finish subroutine of calculations of G function; however, our previous defined 27!! G function ( \int_0^1\exp(-T(1-u^2))(1-u^2)^{n}du ) does not have a stable upward 28!! recursion relation, so we have to give it up, and choose another G function 29!! \int_0^1\exp(-T(1-u^2))u^{2n}du, which is much easier to evaluate compared with 30!! our previous G function; but it requires different recursion relations for 31!! calculating integrals 32!! * finish the tests of the accuray of auxiliary functions, *NOTE* that the referenced 33!! results may not be true sometimes, like for G function, its downward recursion 34!! relation is not stable for large argument, so if the tests involve higher order 35!! G function, the referenced result of 0th order is not always accurate 36!! 37!! 2009-08-23, Bin Gao: 38!! * change subroutine "gen1int_aux_order" with operators and geometric derivatives 39!! as arguments 40!! * add "max_order_aux" (maximum order of auxiliary functions) 41!! * finish subroutine of calculations of Boys function 42!! * add tests of the accuray of auxiliary functions (not finished) 43!! 44!! 2009-08-19, Bin Gao: 45!! * add determination of orders of auxiliary functions from the *input* generalized 46!! one-electron operator and/or *calculated* property integrals in subroutine 47!! gen1int_aux_order 48!! 49!! 2009-08-17, Bin Gao: 50!! * make changes because the modification of *input* generalized one-electron operator 51!! 52!! 2009-08-10, Bin Gao: 53!! * move the definitions of variables related to auxiliary functions from gen1int_data.F90. 54!! 55!! 2009-08-03, Bin Gao: 56!! * make changes according to the modifications of operators. 57!! 58!! 2009-07-23, Bin Gao: 59!! * remove the determination of orders of auxiliary functions from the 60!! property integrals in subroutine gen1int_aux_order, since this 61!! is taken care in gen1int_input_process. 62!! 63!! 2009-06-27, Bin Gao: 64!! * first version 65 66module gen1int_rfunction 67 ! precision 68 use xprecision 69 ! some constants 70 use xconst 71 ! some basic variables & subroutines 72 use xsystem 73 ! basis sets 74 use xbasis, only : max_angular_num 75 ! some basic variables of GEN1INT 76 use gen1int_data 77 implicit none 78 ! maximum order of auxiliary functions, for safety, and also our code may break down 79 ! for large order and argument's auxiliary functions 80 integer, parameter :: max_order_aux = 24 81 ! number of terms of Taylor expansion for the tabulated auxiliary functions 82 integer, parameter :: num_taylor_boys = 6, num_taylor_gfun = 6 83 ! NB: Please do not change the following variables, otherwise, you have to 84 ! update gen1int_aux_boys & gen1int_aux_gfun !! 85 ! number of sampling points for the tabulated auxiliary functions 86 integer, parameter :: num_sample_boys = 120, num_sample_gfun = 120 87 ! we now use sampling points from 0 to 12, interval is 0.1 88 real(xp), parameter :: sample0_boys = zero, interval_boys = tenth, & 89 sample0_gfun = zero, interval_gfun = tenth 90 ! the last sampling points 91 real(xp), parameter :: sample1_boys = sample0_boys+interval_boys*num_sample_boys, & 92 sample1_gfun = sample0_gfun+interval_gfun*num_sample_gfun 93 ! tabulated auxiliary functions 94 real(xp), allocatable :: g1_boys_tab(:,:), g1_gfun_tab(:,:) 95 ! maximum argument for power series expansions 96 ! 97 ! NOTE!!!! 98 ! 1. For Boys function F_n(T), according to our tests, the current asymptotic 99 ! series expansions (gen1int_boys_asymp) can achieve a 10^-12 relative 100 ! accuracy after T>25 (BUT for quite smaller n if T is not large enough), 101 ! we may need to provide a new asymptotic series expansions later ...; 102 ! the modified asymptotic series expansions do not have a stable upward 103 ! recursion relations when n>=2T (n>=T/2 for 45.0<=T<=100.0, but the 104 ! value of Boys function is quite small) according to our tests (negtive 105 ! values may occur); but for the time being, it seems to be enough; 106 ! last, for the power series expansions, we have used downward recursion 107 ! relation, which seems to be very stable according to our tests (0<=T<=400, 108 ! n=10, 1000), even the value of largest order is not correct (for T<24.0 109 ! it accurate enough for n=1000; and it is accurate enough for T<28.0, n=100; 110 ! according to our tests, for larger T, it may become reasonable after n<2T) 111 ! 2. For G function G_n(T), since we have used downward recursion relation after 112 ! the power series expansions of the largest order N, the downward recursion 113 ! relation is not stable for large T; according to our tests, it is stable for 114 ! even 1000th order with T<20.0 (indeed, it becomes better for larger N), but 115 ! it is unstable when T>=20.0, even the value of the largest order is still 116 ! correct (with an appropriate "max_pterms_gfun", it is accurate enough even 117 ! for n=1000, T=400.0); the asymptotic series expansions seems to only work 118 ! for T>30.0, but it, as well as the rational Chebyshev approximation use 119 ! upward recursion relations, which is not stable for large n, according to 120 ! our tests, it seems that n must satisfy n<2T+a (a=10 for T<=40.0; a may become 121 ! smaller for larger T); for the time being, it seems to be enough, since we use 122 ! rational Chebyshev approximation after T>=12, that means we can arrive at least 123 ! 24th order for T==12; however, if higher orders are required, a non-efficient 124 ! way may use power series expansions for each order (oh my god ;-) ), well, 125 ! we may need to improve it larter ... 126 real(xp), parameter :: max_arg_pow_boys = 30.0D0 127 real(xp), parameter :: max_arg_pow_gfun = 30.0D0 128 ! maximum terms for the power series expansions 129!fixme: how to determine max_pterms, which depends on the order & argument 130!fixme: 200 terms are enough for arguments smaller than, like 20 131 integer, parameter :: max_pterms_boys = 200 132 integer, parameter :: max_pterms_gfun = 200 133 ! cut of for auxiliary function 134 real(xp), parameter :: g1_cut_aux = 1.0D-15 135 136 ! At present, this module contains the following public subroutines, 137 ! 1. gen1int_aux_tabulate: tabulate the auxiliary functions 138 ! 2. gen1int_aux_test: test the accuracy of calcualtions of auxiliary functions 139 ! 3. gen1int_aux_order: calculate the maximum order of auxiliary functions 140 ! from the given property, and/or from *calcualted* pre-defined properties 141 ! 4. gen1int_aux_boys: calculate Boys function F_N(T) using (1) Taylor expansions 142 ! and downward recursion relation for small T, and (2) modified asymptotic series 143 ! and upward recursion relation for large T 144 ! 5. gen1int_aux_gfun: calculate G function G_N(T) for operator 1/r^2, using 145 ! (1) Taylor expansions and downward recursion relation for small T, and 146 ! (2) rational Chebyshev approximation and upward recursion relation for 147 ! large T 148 ! 149 ! and the following private subroutines, 150 ! 1. gen1int_boys_power: power series expansions of Boys function 151 ! 2. gen1int_boys_asymp: asymptotic series expansions of Boys funciton 152 ! 3. gen1int_gfun_power: power series expansions of G function 153 ! 4. gen1int_gfun_asymp: asymptotic series expansions of G function 154 public :: gen1int_aux_tabulate, & 155 gen1int_aux_test, & 156 gen1int_aux_order, & 157 gen1int_aux_boys, & 158 gen1int_aux_gfun 159 private :: gen1int_boys_power, & 160 gen1int_boys_asymp, & 161 gen1int_gfun_power, & 162 gen1int_gfun_asymp 163 164 contains 165 166 ! tabulate the auxiliary functions 167 subroutine gen1int_aux_tabulate 168 implicit none 169 ! push current subroutine into the stack 170 call xsub_enter('gen1int_aux_tabulate') 171 ! check if we got max_angular_num 172 if ( max_angular_num.lt.0 ) then 173 write(x_lupri,1000)'Maximum angular number: ',max_angular_num 174 call xstop('Wrong maximum angular number!') 175 end if 176 ! determine the maximum orders of the auxiliary functions from 177 ! *calculated* property integrals, if they have not been done 178 if ( (g1_boys_used .and. op1_order_boys.lt.0) .or. & 179 (g1_gfun_used .and. op1_order_gfun.lt.0) ) & 180 call gen1int_aux_order( check_prop = .true. ) 181 ! now we get the maximum order of the auxiliary functions 182!fixme: maximum orders depend on op1_orders, max_angular_num, and the initial accuracy of F_N(T); 183!fixme: it is used to increase the accuracy (should be less than g1_cut_aux) after downward recursions 184!fixme: so, if max(2,2*max_angular_num) ok?? 185 if (g1_boys_used) then 186 max_order_boys = op1_order_boys + max(2,2*max_angular_num) 187 ! check order 188 if ( max_order_boys.gt.max_order_aux ) then 189 write(x_lupri,1000)'Maximum order of auxiliary functions: ',max_order_aux 190 write(x_lupri,1000)'Order of Boys functions: ',max_order_boys 191 write(x_lupri,1000)'Be careful! The results of high order may be unstable!!' 192 call xstop('Increase max_order_aux in gen1int_aux.F90') 193 end if 194 end if 195 if (g1_gfun_used) then 196 max_order_gfun = op1_order_gfun + max(2,2*max_angular_num) 197 ! check order 198 if ( max_order_gfun.gt.max_order_aux ) then 199 write(x_lupri,1000)'Maximum order of auxiliary functions: ',max_order_aux 200 write(x_lupri,1000)'Order of G functions: ',max_order_gfun 201 write(x_lupri,1000)'Be careful! The results of high order may be unstable!!' 202 call xstop('Increase max_order_aux in gen1int_aux.F90') 203 end if 204 end if 205 ! dump 206 if (xlib_lprint.ge.10) then 207 if (g1_boys_used) write(x_lupri,1000)'Maximum order of Boys function: ',max_order_boys 208 if (g1_gfun_used) write(x_lupri,1000)'Maximum order of G function: ',max_order_gfun 209 write(x_lupri,'()') 210 end if 211 ! tabulate Boys function 212 if (g1_boys_used) then 213 allocate( g1_boys_tab(0:max_order_boys+num_taylor_boys, & 214 0:num_sample_boys), stat=x_ierr) 215 if ( x_ierr.ne.0 ) call xstop('allocate g1_boys_tab!') 216 call gen1int_boys_power( 0, max_order_boys+num_taylor_boys, & 217 sample0_boys, interval_boys, num_sample_boys, g1_boys_tab ) 218 end if 219 ! tabulate G function 220 if (g1_gfun_used) then 221 allocate( g1_gfun_tab(0:max_order_gfun+num_taylor_gfun, & 222 0:num_sample_gfun), stat=x_ierr) 223 if ( x_ierr.ne.0 ) call xstop('allocate g1_gfun_tab!') 224 call gen1int_gfun_power( 0, max_order_gfun+num_taylor_gfun, & 225 sample0_gfun, interval_gfun, num_sample_gfun, g1_gfun_tab ) 226 end if 227 ! pop the stack 228 call xsub_exit 229 return 2301000 format('ATAB ',A,I6) 231 end subroutine gen1int_aux_tabulate 232 233 ! test the accuracy of calcualtions of auxiliary functions 234 subroutine gen1int_aux_test 235!fixme: think a much better referenced results, some software could provide infinite accuracy 236 implicit none 237 !//////////////////////////////////////////////////////////////// 238 ! local variables 239 !//////////////////////////////////////////////////////////////// 240 ! check orders 241 integer, parameter :: my_order = 50 242 ! number of arguments 243 integer, parameter :: my_num_arg = 1000 244 ! incremental recorder for orders 245 integer iord 246 ! incremental recorder for arguments 247 integer iarg 248 ! arguments 249 real(xp) my_arg 250 ! results from power series expansions or asymptotic series expansions 251 real(xp), allocatable :: my_refs(:,:) 252 ! results from subroutines used in the calculations 253 real(xp), allocatable :: my_chks(:) 254 ! error 255 real(xp) my_err 256 ! push current subroutine into the stack 257 call xsub_enter('gen1int_aux_test') 258 write(x_lupri,1010)'Threshold of test: ',g1_cut_aux 259 ! results from power series expansions or asymptotic series expansions 260 allocate(my_refs(0:my_order,0:my_num_arg), stat=x_ierr) 261 if ( x_ierr.ne.0 ) call xstop('allocate my_refs!') 262 ! results from subroutines used in the calculations 263 allocate(my_chks(0:my_order), stat=x_ierr) 264 if ( x_ierr.ne.0 ) call xstop('allocate my_chks!') 265 ! tabulate Boys function 266 if ( allocated(g1_boys_tab) ) deallocate(g1_boys_tab) 267 allocate(g1_boys_tab(0:my_order+num_taylor_boys,0:num_sample_boys), stat=x_ierr) 268 if ( x_ierr.ne.0 ) call xstop('allocate g1_boys_tab!') 269 call gen1int_boys_power( 0, my_order+num_taylor_boys, & 270 sample0_boys, interval_boys, num_sample_boys, g1_boys_tab ) 271 ! test Boys functions 272 write(x_lupri,1000) 273 write(x_lupri,1000)'Test the accuracy of Boys functions ...' 274 write(x_lupri,1000)'Order Argument Result Reference log10(relative error)' 275 ! we use power series expansions and asymptotic series expansions as the referenced 276 ! results, but we should know that they may be not accurate enough sometimes!! 277 call gen1int_boys_power( 0, my_order, zero, one, my_num_arg, my_refs ) 278 ! initialize argument 279 my_arg = zero 280 ! loop over arguments 281 do iarg = 0, my_num_arg 282 ! results from subroutines used in the calculations 283 call gen1int_aux_boys( my_order, my_arg, my_chks ) 284 ! compare them 285 do iord = 0, my_order 286 my_err = abs( (my_chks(iord)-my_refs(iord,iarg))/my_refs(iord,iarg) ) 287 ! Boys function is the monotonically decreasing function of T 288 ! F_n(0) = 1/(2n+1) <= 1 289 if ( my_chks(iord).lt.zero .or. my_refs(iord,iarg).lt.zero .or. & 290 my_chks(iord).gt.one .or. my_refs(iord,iarg).gt.one .or. & 291 my_err.gt.g1_cut_aux ) then 292 write(x_lupri,1020) iord,my_arg,my_chks(iord),my_refs(iord,iarg),log10(my_err) 293 end if 294 end do 295 ! increase argument 296 my_arg = my_arg + one 297 end do 298 ! clean tabulated functions 299 deallocate(g1_boys_tab) 300 ! tabulate G function 301 if ( allocated(g1_gfun_tab) ) deallocate(g1_gfun_tab) 302 allocate(g1_gfun_tab(0:my_order+num_taylor_gfun,0:num_sample_gfun), stat=x_ierr) 303 if ( x_ierr.ne.0 ) call xstop('allocate g1_gfun_tab!') 304 call gen1int_gfun_power( 0, my_order+num_taylor_gfun, & 305 sample0_gfun, interval_gfun, num_sample_gfun, g1_gfun_tab ) 306 ! test G functions 307 write(x_lupri,1000) 308 write(x_lupri,1000)'Test the accuracy of G functions ...' 309 write(x_lupri,1000)'Order Argument Result Reference log10(relative error)' 310 ! we use power series expansions and asymptotic series expansions as the referenced 311 ! results, but we should know that they may be not accurate enough sometimes!! 312 call gen1int_gfun_power( 0, my_order, zero, one, my_num_arg, my_refs ) 313 ! initialize argument 314 my_arg = zero 315 ! loop over arguments 316 do iarg = 0, my_num_arg 317 ! results from subroutines used in the calculations 318 call gen1int_aux_gfun( my_order, my_arg, my_chks ) 319 ! compare them 320 do iord = 0, my_order 321 my_err = abs( (my_chks(iord)-my_refs(iord,iarg))/my_refs(iord,iarg) ) 322 ! G function is the monotonically decreasing function of T 323 ! G_n(0) = 1/(2n+1) <= 1 324 if ( my_chks(iord).lt.zero .or. my_refs(iord,iarg).lt.zero .or. & 325 my_chks(iord).gt.one .or. my_refs(iord,iarg).gt.one .or. & 326 my_err.gt.g1_cut_aux ) then 327 write(x_lupri,1020) iord,my_arg,my_chks(iord),my_refs(iord,iarg),log10(my_err) 328 end if 329 end do 330 ! increase argument 331 my_arg = my_arg + one 332 end do 333 ! clean test results 334 deallocate(my_refs) 335 deallocate(my_chks) 336 ! clean tabulated functions 337 deallocate(g1_gfun_tab) 338 write(x_lupri,'()') 339 ! tabulate auxiliary functions again 340 call gen1int_aux_tabulate 341 ! pop the stack 342 call xsub_exit 343 return 3441000 format('AXTT ',A) 3451010 format('AXTT ',A,E20.12) 3461020 format('AXTT ',I4,F10.3,2E20.12,F10.3) 347 end subroutine gen1int_aux_test 348 349 ! calculate the maximum order of auxiliary functions from the given property, 350 ! and/or from *calcualted* pre-defined properties 351 subroutine gen1int_aux_order( aux_prop, check_prop ) 352 ! property integrals 353 use gen1int_property 354 ! the maximum order of auxiliary function should be the sum of 355 ! (1) powers of multipole moments, 356 ! (2) orders of derivatives of 1/r_C^m (m=1,2) function, 357 ! (3) orders of derivatives with respect to electrons, 358 ! (4) orders of geometric derivatives 359 implicit none 360 ! property 361 type(g1prop), optional, intent(in) :: aux_prop 362 ! calcualte the order of auxiliary functions from *calcualted* pre-defined properties 363 logical, optional, intent(in) :: check_prop 364 !//////////////////////////////////////////////////////////////// 365 ! local variables 366 !//////////////////////////////////////////////////////////////// 367 ! temporary order of auxiliary functions 368 integer tmp_order 369 ! incremental recorder over integrals 370 integer int1 371 ! incremental recorder over operators and *calculated* pre-defined properties 372 integer iop1 373 ! incremental recorder 374 integer ior1 375 ! temporary stuff for the type of r_C function 376 character chk_prop_rfun 377 ! push current subroutine into the stack 378 call xsub_enter('gen1int_aux_order') 379 ! calcualte the order from given operators 380 if ( present(aux_prop) ) then 381 if ( allocated(aux_prop%range_g1int) ) then 382 ! loop over integrals 383 do int1 = aux_prop%range_g1int(1), aux_prop%range_g1int(2) 384 ! loop over operators 385 do iop1 = 1, aux_prop%g1prop_int1(int1)%num_g1opt 386!fixme: consider (1/r_K)*(1/r_L) 387 select case ( aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%r_func(1) ) 388 ! 1/r_C 389 case ('1') 390 g1_boys_used = .true. 391 ! orders of geometric derivatives 392 if ( allocated(aux_prop%order_gdcent) ) then 393 tmp_order = max( sum(aux_prop%order_gdcent), & 394 aux_prop%order_gderv ) 395 else 396 tmp_order = aux_prop%order_gderv 397 end if 398 ! orders of derivatives with respect to electrons 399 tmp_order = tmp_order & 400 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%ederv_order(4) 401 do ior1 = 1, 3 402 ! powers of multipole moments 403 tmp_order = tmp_order & 404 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%mult_power(ior1) 405 ! orders of derivatives of 1/r_C^m (m=1,2) function 406!fixme: consider (1/r_K)*(1/r_L) 407 tmp_order = tmp_order & 408 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%rderv_order(ior1,1) 409 end do 410 ! update order of auxiliary funtions 411 if ( tmp_order.gt.op1_order_boys ) op1_order_boys = tmp_order 412 ! 1/r_C^2 413 case ('2') 414 g1_gfun_used = .true. 415 ! orders of geometric derivatives 416 if ( allocated(aux_prop%order_gdcent) ) then 417 tmp_order = max( sum(aux_prop%order_gdcent), & 418 aux_prop%order_gderv ) 419 else 420 tmp_order = aux_prop%order_gderv 421 end if 422 ! orders of derivatives with respect to electrons 423 tmp_order = tmp_order & 424 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%ederv_order(4) 425 do ior1 = 1, 3 426 ! powers of multipole moments 427 tmp_order = tmp_order & 428 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%mult_power(ior1) 429 ! orders of derivatives of 1/r_C^m (m=1,2) function 430!fixme: consider (1/r_K)*(1/r_L) 431 tmp_order = tmp_order & 432 + aux_prop%g1prop_int1(int1)%g1int_opt(iop1)%rderv_order(ior1,1) 433 end do 434 ! update order of auxiliary funtions 435 if ( tmp_order.gt.op1_order_gfun ) op1_order_gfun = tmp_order 436 end select 437 end do ! loop over operators 438 end do ! loop over integrals 439 end if 440 end if 441 ! calcualte the order from *calcualted* pre-defined properties 442 if ( present(check_prop) ) then 443 if ( check_prop ) then 444 do iop1 = 1, num_calc_prop 445 ! get the order from the operators of this *calcualted* pre-defined property 446 call gen1int_proper_create( prop_name = def_prop(iop1)%prop_name, & 447 prop_rfun = chk_prop_rfun, & 448 prop_aux_order = tmp_order ) 449 ! type of r_C function 450 select case (chk_prop_rfun) 451 case ('1') 452 g1_boys_used = .true. 453 ! consider the order of geometric derivatives 454 if ( allocated(def_prop(iop1)%order_gdcent) ) then 455 tmp_order = tmp_order & 456 + max( sum(def_prop(iop1)%order_gdcent), def_prop(iop1)%order_gderv ) 457 else 458 tmp_order = tmp_order + def_prop(iop1)%order_gderv 459 end if 460 ! update order of auxiliary funtions 461 if ( tmp_order.gt.op1_order_boys ) op1_order_boys = tmp_order 462 case ('2') 463 g1_gfun_used = .true. 464 ! consider the order of geometric derivatives 465 if ( allocated(def_prop(iop1)%order_gdcent) ) then 466 tmp_order = tmp_order & 467 + max( sum(def_prop(iop1)%order_gdcent), def_prop(iop1)%order_gderv ) 468 else 469 tmp_order = tmp_order + def_prop(iop1)%order_gderv 470 end if 471 ! update order of auxiliary funtions 472 if ( tmp_order.gt.op1_order_gfun ) op1_order_gfun = tmp_order 473 end select 474 end do 475 end if 476 end if 477 ! pop the stack 478 call xsub_exit 479 return 4801000 format('AXOR ',A) 481 end subroutine gen1int_aux_order 482 483 ! calculate Boys function F_N(T) using (1) Taylor expansions and downward 484 ! recursion relation for small T, and (2) modified asymptotic series and 485 ! upward recursion relation for large T 486 subroutine gen1int_aux_boys(order_boys,arg_boys,val_boys) 487 implicit none 488 ! order of Boys function 489 integer, intent(in) :: order_boys 490 ! argument of Boys function 491 real(xp), intent(in) :: arg_boys 492 ! value of Boys function 493 real(xp), intent(out) :: val_boys(0:order_boys) 494 !//////////////////////////////////////////////////////////////// 495 ! local variables 496 !//////////////////////////////////////////////////////////////// 497 ! recorder for the argument 498 integer iarg 499 ! recorder for the orders 500 integer iord 501 ! (T*-T) 502 real(xp) my_dif_arg 503 ! twice of the argument 504 real(xp) my_arg2 505 ! F_J(T) 506 real(xp) val_boys_max 507 ! exp(-T) 508 real(xp) my_pre_arg 509 ! 1/T 510 real(xp) my_iarg 511 ! 1/(2T) 512 real(xp) my_i2arg 513 ! divisor in the downward recursion relations 514 real(xp) div_down 515 ! four term modified asymptotic series 516 real(xp), parameter :: g_asym3(0:3) = (/0.4999489092D0,-0.2473631686D0, & 517 0.321180909D0,-0.3811559346D0/) 518 ! check if we got the tabulated auxiliary functions 519 if ( g1_boys_used .and. .not. allocated(g1_boys_tab) ) & 520 call xstop('Boys function is not tabulated!') 521 ! 0 <= T < 12: using Taylor expansion and downward recursion relation 522 if ( arg_boys.lt.sample1_boys ) then 523 ! the nearest T*, NB: since sample0_boys == 0! 524 iarg = nint(arg_boys/interval_boys) 525 ! T*-T 526 my_dif_arg = iarg*interval_boys - arg_boys 527 ! (T*-T)^k/k!, k==0 528 my_pre_arg = one 529 ! the 1st Taylor's term 530!N val_boys_max = g1_boys_tab(max_order_boys,iarg) 531 val_boys(order_boys) = g1_boys_tab(order_boys,iarg) 532 ! the left Taylor's terms 533 do iord = 1, num_taylor_boys 534 my_pre_arg = my_pre_arg*my_dif_arg/xfloat(iord) 535!N val_boys_max = val_boys_max + my_pre_arg*g1_boys_tab(iord+max_order_boys,iarg) 536 val_boys(order_boys) = val_boys(order_boys) + my_pre_arg*g1_boys_tab(iord+order_boys,iarg) 537 end do 538 ! exp(-T) 539 my_pre_arg = exp(-arg_boys) 540 ! 2*T 541 my_arg2 = arg_boys + arg_boys 542!N ! use downward recursion relations to get F_j(T) from F_J(T) (val_boys_max) 543!N div_down = xfloat(max_order_boys+max_order_boys+1) 544 div_down = xfloat(order_boys+order_boys+1) 545!N do iord = max_order_boys-1, order_boys, -1 546!N div_down = div_down - two 547!N val_boys(order_boys) = (my_arg2*val_boys_max+my_pre_arg)/div_down 548!N val_boys_max = val_boys(order_boys) 549!N end do 550 ! use downward recursion relations for other F_j(T) 551 do iord = order_boys-1, 0, -1 552 div_down = div_down - two 553 val_boys(iord) = (my_arg2*val_boys(iord+1)+my_pre_arg)/div_down 554 end do ! loop over orders 555 ! T >= 12: using modified asymptotic series expansions and upward recursion relation 556 ! for T <= 2J+36, we use four asymptotic terms, and exact upward recursion relation 557!fixme: this modified asymptotic series expansion seems to be bad for larger order, like 50th, 558!fixme: with smaller arguments, like 12.0, 13.0, ... 559 else if ( arg_boys.le.xfloat(max_order_boys+max_order_boys+36) ) then 560 if ( order_boys.ge.floor(arg_boys+arg_boys) ) then 561 write(x_lupri,1000)'The maximum order of Boys functions: ',order_boys 562 write(x_lupri,1010)'The argument of Boys functions: ',arg_boys 563 write(x_lupri,1000)'The upward recursion relation is unstable for them!' 564#ifndef TEST_AUX 565 call xstop('Write to author asking them to improve it!') 566#endif 567 end if 568 ! 1/T 569 my_iarg = one/arg_boys 570 ! exp(-T) 571 my_pre_arg = exp(-arg_boys) 572 ! g 573 val_boys_max = g_asym3(0) + my_iarg*( g_asym3(1) + my_iarg*( g_asym3(2) + my_iarg*g_asym3(3) ) ) 574 ! F_0(T) = sqrt(pi/T)/2 - exp(-T)*g/T 575 val_boys(0) = half*sqrt(pi*my_iarg) - my_pre_arg*val_boys_max*my_iarg 576 ! 1/(2T) 577 my_i2arg = half*my_iarg 578 ! exp(-T)/(2T) 579 my_pre_arg = my_i2arg*my_pre_arg 580 ! upward recursion relation for other F_j(T) 581 do iord = 1, order_boys 582 val_boys(iord) = my_i2arg*val_boys(iord-1) - my_pre_arg 583 my_i2arg = my_iarg + my_i2arg 584 end do 585 ! for T > 2J+36, we use F_0(T) = sqrt(pi/T)/2, and approximate upward recursion relation 586 ! F_{j+1}(T) = (2T)^{-1}(2j+1)F_j(T) 587 else 588 if ( order_boys.ge.floor(half*arg_boys) ) then 589 write(x_lupri,1000)'The maximum order of Boys functions: ',order_boys 590 write(x_lupri,1010)'The argument of Boys functions: ',arg_boys 591 write(x_lupri,1000)'The upward recursion relation is unstable for them!' 592#ifndef TEST_AUX 593 call xstop('Write to author asking them to improve it!') 594#endif 595 end if 596 ! 1/T 597 my_iarg = one/arg_boys 598 val_boys(0) = half*sqrt(pi*my_iarg) 599 ! 1/(2T) 600 my_i2arg = half*my_iarg 601 ! upward approximate recursion relation for other F_j(T) 602 do iord = 1, order_boys 603 val_boys(iord) = my_i2arg*val_boys(iord-1) 604 my_i2arg = my_iarg + my_i2arg 605 end do 606 end if 607 return 6081000 format('BOYS ',A,I6) 6091010 format('BOYS ',A,F14.8) 610 end subroutine gen1int_aux_boys 611 612 ! calculate G function G_N(T) for operator 1/r^2, using (1) Taylor expansions 613 ! and downward recursion relation for small T, and (2) rational Chebyshev 614 ! approximation and upward recursion relation for large T 615 subroutine gen1int_aux_gfun( order_gfun, arg_gfun, val_gfun ) 616 implicit none 617 ! order of G function 618 integer, intent(in) :: order_gfun 619 ! argument of G function 620 real(xp), intent(in) :: arg_gfun 621 ! values of G function 622 real(xp), intent(out) :: val_gfun(0:order_gfun) 623 !//////////////////////////////////////////////////////////////// 624 ! local variables 625 !//////////////////////////////////////////////////////////////// 626 ! recorder for the argument 627 integer iarg 628 ! recorder for the orders 629 integer iord 630 ! (T*-T) 631 real(xp) my_dif_arg 632 ! twice of the argument 633 real(xp) my_arg2 634 ! 1/T 635 real(xp) my_iarg 636 ! 1/(2T) 637 real(xp) my_i2arg 638 ! k/T 639 real(xp) my_pre_arg 640 ! divisor in the downward recursion relations 641 real(xp) div_down 642 ! coefficients of rational Chebyshev approximation of 0th order G function 643 real(xp), parameter :: cheb0p1(0:5) = (/ 0.1828591449264953D0, & 644 -0.0208111908526560D0, & 645 0.1333867391625493D0, & 646 -0.0069563805854603D0, & 647 -0.0006399406464297D0, & 648 0.0002262208851365D0/) 649 real(xp), parameter :: cheb0q1(1:6) = (/-0.5001558049033900D0, & 650 -0.0422332795713542D0, & 651 0.2634740589566963D0, & 652 -0.0130363366131583D0, & 653 -0.0015170455558345D0, & 654 0.0004525493560858D0/) 655 real(xp), parameter :: cheb0p2(0:4) = (/-0.0237209201842862D0, & 656 0.0369293462259353D0, & 657 -0.1782265603862903D0, & 658 0.0305190498240735D0, & 659 -0.0019473305442354D0/) 660 real(xp), parameter :: cheb0q2(1:5) = (/-0.0929707538344341D0, & 661 0.2362992390268955D0, & 662 -0.3854264829228933D0, & 663 0.0629939171730423D0, & 664 -0.0038947347910983D0/) 665 real(xp), parameter :: cheb0p3(0:2) = (/-0.7973613445962394D0, & 666 -0.1763813954302433D0, & 667 0.0562314564328704D0/) 668 real(xp), parameter :: cheb0q3(1:4) = (/-1.4886748949591277D0, & 669 -0.4084032309951818D0, & 670 0.1124504768032380D0, & 671 0.0000001053836247D0/) 672 real(xp), parameter :: cheb0p4(0:2) = (/-1.6901358291748010D0, & 673 0.9095465039538543D0, & 674 -0.0821332073172411D0/) 675 real(xp), parameter :: cheb0q4(1:3) = (/-4.2082066572381995D0, & 676 1.9012356348172637D0, & 677 -0.1642664818611426D0/) 678 real(xp), parameter :: cheb0p5(0:2) = (/-1.7168354719492460D0, & 679 0.9510390010759011D0, & 680 -0.0888365521252258D0/) 681 real(xp), parameter :: cheb0q5(1:3) = (/-4.2963007179996406D0, & 682 1.9909219459253031D0, & 683 -0.1776731550443674D0/) 684 real(xp), parameter :: cheb0p6(0:2) = (/-1.8429493331531948D0, & 685 1.1598159456555421D0, & 686 -0.1246514775701103D0/) 687 real(xp), parameter :: cheb0q6(1:3) = (/-4.7212032569199138D0, & 688 2.4442852443824989D0, & 689 -0.2493029652188083D0/) 690 real(xp), parameter :: cheb0p7(0:2) = (/-1.9744562443420302D0, & 691 1.3986614326037992D0, & 692 -0.1691251258024091D0/) 693 real(xp), parameter :: cheb0q7(1:3) = (/-5.1784724951222785D0, & 694 2.9664481966513980D0, & 695 -0.3382502523273182D0/) 696 real(xp), parameter :: cheb0p8(0:2) = (/-2.1099021081579972D0, & 697 1.6701190095543932D0, & 698 -0.2240542865988810D0/) 699 real(xp), parameter :: cheb0q8(1:3) = (/-5.6658694074754621D0, & 700 3.5642923069732153D0, & 701 -0.4481085731988854D0/) 702 real(xp), parameter :: cheb0p9(0:1) = (/-1.6616791513685247D0, & 703 0.6616886439852467D0/) 704 real(xp), parameter :: cheb0q9(1:2) = (/-3.9850469515266829D0, & 705 1.3233772879713688D0/) 706 ! check if we got the tabulated auxiliary functions 707 if ( g1_gfun_used .and. .not. allocated(g1_gfun_tab) ) & 708 call xstop('G function is not tabulated!') 709 ! 0 <= T < 12: using Taylor expansions for the largest order and downward 710 ! recursion relations for others 711 if ( arg_gfun.lt.sample1_gfun ) then 712 ! the nearest T*, NB: since sample0_gfun == 0! 713 iarg = nint(arg_gfun/interval_gfun) 714 ! T*-T 715 my_dif_arg = iarg*interval_gfun - arg_gfun 716 ! (T*-T)^k/k!, k==0 717 my_iarg = one 718 ! the 1st Taylor's term 719 val_gfun(order_gfun) = g1_gfun_tab(order_gfun,iarg) 720 ! the left Taylor's terms 721 do iord = 1, num_taylor_gfun 722 my_iarg = my_iarg*my_dif_arg/xfloat(iord) 723 val_gfun(order_gfun) = val_gfun(order_gfun) + my_iarg*g1_gfun_tab(iord+order_gfun,iarg) 724 end do 725 ! left orders using downward recursions 726 ! 2*T 727 my_arg2 = arg_gfun + arg_gfun 728 ! use downward recursion relations for other orders 729 div_down = xfloat(order_gfun+order_gfun+1) 730 do iord = order_gfun-1, 0, -1 731 div_down = div_down - two 732 val_gfun(iord) = (one-my_arg2*val_gfun(iord+1))/div_down 733 end do ! loop over orders 734 ! T >= 12: using rational Chebyshev approximation for the 0th order 735 ! and upward recursion relations for others 736!fixme: the upward recursion relation is unstable for larger order (seems to be >=2T+10) 737!fixme: well, it may not be exactly 2T+10, but depends on T (larger T smaller n) 738 else 739 if ( order_gfun.ge.floor(arg_gfun+arg_gfun)+10 ) then 740 write(x_lupri,1000)'The maximum order of G functions: ',order_gfun 741 write(x_lupri,1010)'The argument of G functions: ',arg_gfun 742 write(x_lupri,1000)'The upward recursion relation is unstable for them!' 743#ifndef TEST_AUX 744 call xstop('Write to author asking them to improve it!') 745#endif 746 end if 747!fixme: we may use asymptotic series expansions for very large T, like > 10,000 748 !----------- 749 ! 0th order 750 !----------- 751 ! 12 <= T <= 15.8 (fitted from 12 <= T <= 16) 752 if ( arg_gfun.le. 15.8D0 ) then 753 my_iarg = cheb0p1(0) + arg_gfun*( cheb0p1(1) + arg_gfun*( cheb0p1(2) & 754 + arg_gfun*( cheb0p1(3) + arg_gfun*( cheb0p1(4) + arg_gfun*cheb0p1(5) ) ) ) ) 755 my_i2arg = one + arg_gfun*( cheb0q1(1) + arg_gfun*( cheb0q1(2) & 756 + arg_gfun*( cheb0q1(3) + arg_gfun*( cheb0q1(4) + arg_gfun*( cheb0q1(5) & 757 + arg_gfun*cheb0q1(6) ) ) ) ) ) 758 ! 15.8 < T <= 18.9 (fitted from 16 <= T <= 19) 759 else if ( arg_gfun.le.18.9D0 ) then 760 my_iarg = cheb0p2(0) + arg_gfun*( cheb0p2(1) + arg_gfun*( cheb0p2(2) & 761 + arg_gfun*( cheb0p2(3) + arg_gfun*cheb0p2(4) ) ) ) 762 my_i2arg = one + arg_gfun*( cheb0q2(1) + arg_gfun*( cheb0q2(2) & 763 + arg_gfun*( cheb0q2(3) + arg_gfun*( cheb0q2(4) + arg_gfun*cheb0q2(5) ) ) ) ) 764 ! 18.9 < T <= 21.9 (fitted from 19 <= T <= 22) 765 else if ( arg_gfun.le.21.9D0 ) then 766 my_iarg = cheb0p3(0) + arg_gfun*( cheb0p3(1) + arg_gfun*cheb0p3(2) ) 767 my_i2arg = one + arg_gfun*( cheb0q3(1) + arg_gfun*( cheb0q3(2) & 768 + arg_gfun*( cheb0q3(3) + arg_gfun*cheb0q3(4) ) ) ) 769 ! 21.9 < T <= 25 (fitted from 22 <= T <= 25) 770 else if ( arg_gfun.le.25.0D0) then 771 my_iarg = cheb0p4(0) + arg_gfun*( cheb0p4(1) + arg_gfun*cheb0p4(2) ) 772 my_i2arg = one + arg_gfun*( cheb0q4(1) + arg_gfun*( cheb0q4(2) + arg_gfun*cheb0q4(3) ) ) 773 ! 25 < T <= 30 (fitted from 25 <= T <= 30) 774 else if ( arg_gfun.le.30.0D0 ) then 775 my_iarg = cheb0p5(0) + arg_gfun*( cheb0p5(1) + arg_gfun*cheb0p5(2) ) 776 my_i2arg = one + arg_gfun*( cheb0q5(1) + arg_gfun*( cheb0q5(2) + arg_gfun*cheb0q5(3) ) ) 777 ! 30 < T <= 39 (fitted from 30 <= T <= 39) 778 else if ( arg_gfun.le.39.0D0 ) then 779 my_iarg = cheb0p6(0) + arg_gfun*( cheb0p6(1) + arg_gfun*cheb0p6(2) ) 780 my_i2arg = one + arg_gfun*( cheb0q6(1) + arg_gfun*( cheb0q6(2) + arg_gfun*cheb0q6(3) ) ) 781 ! 39 < T <= 65.3 (fitted from 39 <= T <= 65) 782 else if ( arg_gfun.le.65.3D0 ) then 783 my_iarg = cheb0p7(0) + arg_gfun*( cheb0p7(1) + arg_gfun*cheb0p7(2) ) 784 my_i2arg = one + arg_gfun*( cheb0q7(1) + arg_gfun*( cheb0q7(2) + arg_gfun*cheb0q7(3) ) ) 785 ! 65.3 < T <= 614.9 (fitted from 65 <= T <= 580) 786 else if ( arg_gfun.le.614.9D0 ) then 787 my_iarg = cheb0p8(0) + arg_gfun*( cheb0p8(1) + arg_gfun*cheb0p8(2) ) 788 my_i2arg = one + arg_gfun*( cheb0q8(1) + arg_gfun*( cheb0q8(2) + arg_gfun*cheb0q8(3) ) ) 789 ! T > 614.9 (fitted from 600 <= T <= 10000) 790 else 791 my_iarg = cheb0p9(0) + arg_gfun*cheb0p9(1) 792 my_i2arg = one + arg_gfun*( cheb0q9(1) + arg_gfun*cheb0q9(2) ) 793 end if 794 val_gfun(0) = my_iarg/my_i2arg 795 ! 1/T 796 my_i2arg = half/arg_gfun 797 my_iarg = my_i2arg + my_i2arg 798 ! recorder of (2n+1)/2T 799 div_down = -my_i2arg 800 ! loop over other orders using upward recursion relations 801 do iord = 0, order_gfun-1 802 ! plus 1/T 803 div_down = div_down + my_iarg 804 ! perform recursion relation 805 val_gfun(iord+1) = my_i2arg - div_down*val_gfun(iord) 806 end do ! loop over other orders 807 end if 808 return 8091000 format('GFUN ',A,I6) 8101010 format('GFUN ',A,F14.8) 811 end subroutine gen1int_aux_gfun 812 813 ! power series expansions of Boys function (for small argument!!) 814 subroutine gen1int_boys_power( strt_order, end_order, strt_arg, step_arg, num_arg, val_boys ) 815 ! the power series expansions can be found, for example, in: 816 ! V. R. Saunders. An introduction to molecular integral evaluation. 817 ! In G.H.F. Diercksen, B.T. Sutcliffe, and A. Veillard, editors, 818 ! Computational Techniques in Quantum Chemistry and Molecular Physics, 819 ! Eq. (39), page 347, 1975. 820 implicit none 821 ! start order of Boys function 822 integer, intent(in) :: strt_order 823 ! end order of Boys function 824 integer, intent(in) :: end_order 825 ! start argument of Boys function 826 real(xp), intent(in) :: strt_arg 827 ! step of the argument 828 real(xp), intent(in) :: step_arg 829 ! number of arguments 830 integer, intent(in) :: num_arg 831 ! values of Boys function 832 real(xp), intent(out) :: val_boys(strt_order:end_order,0:num_arg) 833 !//////////////////////////////////////////////////////////////// 834 ! local variables 835 !//////////////////////////////////////////////////////////////// 836 ! argument of Boys function at each step 837 real(xp) my_arg 838 ! twice of the argument 839 real(xp) my_arg2 840 ! incremental recorder for the orders 841 integer iord 842 ! incremental recorder for the arguments 843 integer iarg 844 ! value of power series expansions' term 845 real(xp) val_power 846 ! divisor in the term of power series expansions 847 real(xp) div_power 848 ! incremental recorder for the summation 849 integer ipower 850 ! exp(-T) 851 real(xp) my_pre_arg 852 ! push current subroutine into the stack 853 call xsub_enter('gen1int_boys_power') 854 ! the minimum argument 855 my_arg = min( strt_arg+step_arg*xfloat(num_arg), strt_arg ) 856 ! the maximum argument 857 my_arg2 = max( strt_arg+step_arg*xfloat(num_arg), strt_arg ) 858 ! check if the value of largest order is reasonable 859 if ( my_arg2.gt.twenty3 .and. end_order.ge.nint(my_arg+my_arg) ) then 860 write(x_lupri,1000)'Your minimum argument of Boys function: ',my_arg 861 write(x_lupri,1010)'Your maximum order of Boys function: ',end_order 862 write(x_lupri,1000)'Warning! Some values of the largest order may not be accurate or correct!!' 863 end if 864 ! loop over number of arguments 865 do iarg = 0, num_arg 866 ! for accuracy, this is much better than add "step_arg" to "my_arg" at each time 867 my_arg = strt_arg + step_arg*xfloat(iarg) 868 ! if the argument exceeds "max_arg_pow_boys", we change to asymptotic series expansions 869 if ( my_arg.gt.max_arg_pow_boys ) then 870 if (xlib_lprint.ge.200) then 871 write(x_lupri,1000)'Your argument of Boys function: ',my_arg 872 write(x_lupri,1000)'It is too large, we change to asymptotic series expansions' 873 end if 874 call gen1int_boys_asymp( strt_order, end_order, my_arg, step_arg, 0, val_boys(:,iarg) ) 875 else 876 ! 2*T 877 my_arg2 = my_arg + my_arg 878 ! we use power series expansions for the largest order, and use downward recursion 879 ! relations for others 880 div_power = xfloat(end_order+end_order+1) 881 val_power = one/div_power 882 ! fisrt term 883 val_boys(end_order,iarg) = val_power 884 ! loop over power series expansions 885 do ipower = 1, max_pterms_boys 886 div_power = div_power + two 887 val_power = val_power*my_arg2/div_power 888 val_boys(end_order,iarg) = val_boys(end_order,iarg) + val_power 889 if (val_power.le.g1_cut_aux) exit 890 end do 891 ! exp(-T) 892 my_pre_arg = exp(-my_arg) 893 val_boys(end_order,iarg) = my_pre_arg*val_boys(end_order,iarg) 894 ! use downward recursion relations for other orders 895 div_power = xfloat(end_order+end_order+1) 896 do iord = end_order-1, strt_order, -1 897 div_power = div_power - two 898 val_boys(iord,iarg) = (my_arg2*val_boys(iord+1,iarg)+my_pre_arg)/div_power 899 end do ! loop over orders 900 end if 901 end do ! loop over number of arguments 902 ! pop the stack 903 call xsub_exit 904 return 9051000 format('BPOW ',A,F14.8) 9061010 format('BPOW ',A,I6) 907 end subroutine gen1int_boys_power 908 909 ! asymptotic series expansions of Boys funciton (for large argument!!) 910!fixme: we may think about another asymptotic series expansions 911 subroutine gen1int_boys_asymp( strt_order, end_order, strt_arg, step_arg, num_arg, val_boys ) 912 ! we use F_n(T) = (2n-1)!!/2^(n+1)*sqrt(pi/T^(2n+1)), see, for example, 913 ! Trygve Helgaker, Poul Jorgensen, Jeppe Olsen, Molecular Electronic Structure Theory, 914 ! Eq. (9.8.9), page 365 915 implicit none 916 ! start order of Boys function 917 integer, intent(in) :: strt_order 918 ! end order of Boys function 919 integer, intent(in) :: end_order 920 ! start argument of Boys function 921 real(xp), intent(in) :: strt_arg 922 ! step of the argument 923 real(xp), intent(in) :: step_arg 924 ! number of arguments 925 integer, intent(in) :: num_arg 926 ! values of Boys function 927 real(xp), intent(out) :: val_boys(strt_order:end_order,0:num_arg) 928 !//////////////////////////////////////////////////////////////// 929 ! local variables 930 !//////////////////////////////////////////////////////////////// 931 ! argument of Boys function at each step 932 real(xp) my_arg 933 ! 1/T 934 real(xp) my_iarg 935 ! incremental recorder for the orders 936 integer iord 937 ! incremental recorder for the arguments 938 integer iarg 939 ! dividend in the asymptotic term 940 real(xp) div_asymp 941 ! incremental recorder for the asymptotic series expansions 942 integer iasymp 943 ! sqrt(pi/T)/2 944 real(xp) my_pre_arg 945 ! push current subroutine into the stack 946 call xsub_enter('gen1int_boys_asymp') 947 ! the minimum argument 948 my_arg = min( strt_arg+step_arg*xfloat(num_arg), strt_arg ) 949 ! check if the minimum argument is large enough 950 if ( my_arg.le.max_arg_pow_boys ) then 951 write(x_lupri,1000)'Your minimum argument of Boys function: ',my_arg 952 call xstop('Too small argument for asymptotic series expansions!') 953 end if 954 ! loop over number of arguments 955 do iarg = 0, num_arg 956 my_arg = strt_arg + step_arg*xfloat(iarg) 957 ! 1/T 958 my_iarg = one/my_arg 959 ! -1/(2T) 960 div_asymp = -half*my_iarg 961 ! sqrt(pi/T)/2 962 val_boys(strt_order,iarg) = half*sqrt(pi/my_arg) 963 ! the smallest order 964 do iasymp = 1, strt_order 965 ! plus 1/T 966 div_asymp = div_asymp + my_iarg 967 ! times (2n-1)/(2T) 968 val_boys(strt_order,iarg) = val_boys(strt_order,iarg)*div_asymp 969 end do 970 ! the left orders 971 do iord = strt_order, end_order-1 972 ! plus 1/T 973 div_asymp = div_asymp + my_iarg 974 val_boys(iord+1,iarg) = val_boys(iord,iarg)*div_asymp 975 end do ! loop over orders 976 end do ! loop over number of arguments 977 ! pop the stack 978 call xsub_exit 979 return 9801000 format('BASY ',A,F14.8) 981 end subroutine gen1int_boys_asymp 982 983 ! power series expansions of G function (for small argument!!) 984 ! NB: the downward recursion relation is very bad for large arguments!! 985 subroutine gen1int_gfun_power( strt_order, end_order, strt_arg, step_arg, num_arg, val_gfun ) 986!fixme: give the reference 987 implicit none 988 ! start order of G function 989 integer, intent(in) :: strt_order 990 ! end order of G function 991 integer, intent(in) :: end_order 992 ! start argument of G function 993 real(xp), intent(in) :: strt_arg 994 ! step of the argument 995 real(xp), intent(in) :: step_arg 996 ! number of arguments 997 integer, intent(in) :: num_arg 998 ! values of G function 999 real(xp), intent(out) :: val_gfun(strt_order:end_order,0:num_arg) 1000 !//////////////////////////////////////////////////////////////// 1001 ! local variables 1002 !//////////////////////////////////////////////////////////////// 1003 ! argument of G function at each step 1004 real(xp) my_arg 1005 ! twice of the argument 1006 real(xp) my_arg2 1007 ! incremental recorder for the orders 1008 integer iord 1009 ! incremental recorder for the arguments 1010 integer iarg 1011 ! value of power series expansions' term 1012 real(xp) val_power 1013 ! divisor in the term of power series expansions 1014 real(xp) div_power 1015 ! incremental recorder for the summation 1016 integer ipower 1017 ! exp(-T) 1018 real(xp) my_pre_arg 1019 ! push current subroutine into the stack 1020 call xsub_enter('gen1int_gfun_power') 1021 ! the maximum argument 1022 my_arg = max( strt_arg+step_arg*xfloat(num_arg), strt_arg ) 1023 ! check if the value of largest order is reasonable 1024 if ( my_arg.ge.twenty ) then 1025 write(x_lupri,1000)'Your maximum argument of G function: ',my_arg 1026 write(x_lupri,1010)'Your maximum order of G function: ',end_order 1027 write(x_lupri,1000)'Warning! The downward recursion relation may be unstable!!' 1028 end if 1029 ! loop over number of arguments 1030 do iarg = 0, num_arg 1031 ! for accuracy, this is much better than add "step_arg" to "my_arg" at each time 1032 my_arg = strt_arg + step_arg*xfloat(iarg) 1033 ! if the argument exceeds "max_arg_pow_gfun", we change to asymptotic series expansions 1034 if ( my_arg.gt.max_arg_pow_gfun ) then 1035 if (xlib_lprint.ge.200) then 1036 write(x_lupri,1000)'Your argument of G function: ',my_arg 1037 write(x_lupri,1000)'It is too large, we change to asymptotic series expansions' 1038 end if 1039 call gen1int_gfun_asymp( strt_order, end_order, my_arg, step_arg, 0, val_gfun(:,iarg) ) 1040 else 1041 ! we use power series expansions for the largest order, and use downward recursion 1042 ! relations for others 1043 ! recorder of 2n+2k+1 1044 div_power = xfloat(end_order+end_order+1) 1045 ! recorder for T^k/k! 1046 my_arg2 = one 1047 ! fisrt term 1/(2n+1) 1048 val_gfun(end_order,iarg) = my_arg2/div_power 1049 ! loop over power series expansions 1050 do ipower = 1, max_pterms_gfun 1051 div_power = div_power + two 1052 ! time T/k 1053 my_arg2 = my_arg2*my_arg/xfloat(ipower) 1054 ! divided by 2n+2k+1 1055 val_power = my_arg2/div_power 1056 val_gfun(end_order,iarg) = val_gfun(end_order,iarg) + val_power 1057 if (val_power.le.g1_cut_aux) exit 1058 end do 1059 ! exp(-T) 1060 my_pre_arg = exp(-my_arg) 1061 val_gfun(end_order,iarg) = val_gfun(end_order,iarg)*my_pre_arg 1062 ! 2*T 1063 my_arg2 = my_arg + my_arg 1064 ! use downward recursion relations for other orders 1065 div_power = xfloat(end_order+end_order+1) 1066 do iord = end_order-1, strt_order, -1 1067 div_power = div_power - two 1068 val_gfun(iord,iarg) = (one-my_arg2*val_gfun(iord+1,iarg))/div_power 1069 end do ! loop over orders 1070 end if 1071 end do ! loop over number of arguments 1072 ! pop the stack 1073 call xsub_exit 1074 return 10751000 format('GPOW ',A,F14.8) 10761010 format('GPOW ',A,I6) 1077 end subroutine gen1int_gfun_power 1078 1079 ! asymptotic series expansions of G function (for large argument!! > 30) 1080 subroutine gen1int_gfun_asymp( strt_order, end_order, strt_arg, step_arg, num_arg, val_gfun ) 1081 ! we use G_n(T) = 0.5*sum_{j=0}^{J} (0.5-n)_{j} (1/T)^{j+1} = sum_{j=0}^{J} V_{j} 1082 ! V_{0} = 1/(2T), V_{j} = (0.5-n+j-1)*V_{j-1}/T 1083!fixme: give the reference 1084 implicit none 1085 ! start order of G function 1086 integer, intent(in) :: strt_order 1087 ! end order of G function 1088 integer, intent(in) :: end_order 1089 ! start argument of G function 1090 real(xp), intent(in) :: strt_arg 1091 ! step of the argument 1092 real(xp), intent(in) :: step_arg 1093 ! number of arguments 1094 integer, intent(in) :: num_arg 1095 ! values of G function 1096 real(xp), intent(out) :: val_gfun(strt_order:end_order,0:num_arg) 1097 !//////////////////////////////////////////////////////////////// 1098 ! local variables 1099 !//////////////////////////////////////////////////////////////// 1100 ! argument of G function at each step 1101 real(xp) my_arg 1102 ! 1/T 1103 real(xp) my_iarg 1104 ! 1/(2T) 1105 real(xp) my_i2arg 1106 ! dividend in the asymptotic term 1107 real(xp) div_asymp 1108 ! incremental recorder for the orders 1109 integer iord 1110 ! incremental recorder for the arguments 1111 integer iarg 1112 ! value of asymptotic term 1113 real(xp) val_asymp 1114 ! maximum terms for asymptotic series expansions, for safety 1115!fixme: how to determine it, depends on the order and argument 1116 integer, parameter :: max_aterms_gfun = 100 1117 ! incremental recorder for the asymptotic series expansions 1118 integer iasymp 1119 ! push current subroutine into the stack 1120 call xsub_enter('gen1int_gfun_asymp') 1121 ! the minimum argument 1122 my_arg = min( strt_arg+step_arg*xfloat(num_arg), strt_arg ) 1123 ! check if the minimum argument is large enough 1124 if ( my_arg.le.max_arg_pow_gfun ) then 1125 write(x_lupri,1000)'Your minimum argument of G function: ',my_arg 1126 call xstop('Too small argument for asymptotic series expansions!') 1127 end if 1128 ! check if the minimum argument and maximum order are good for a stable upward recursion relation 1129 if ( end_order.ge.floor(my_arg+my_arg)+10 ) then 1130 write(x_lupri,1010)'The maximum order of G functions: ',end_order 1131 write(x_lupri,1000)'The minimum argument of G functions: ',my_arg 1132 write(x_lupri,1000)'The upward recursion relation is unstable for them!' 1133#ifndef TEST_AUX 1134 call xstop('Write to author asking them to improve it!') 1135#endif 1136 end if 1137 ! loop over number of arguments 1138 do iarg = 0, num_arg 1139 my_arg = strt_arg + step_arg*xfloat(iarg) 1140 ! we use asymptotic series expansions for the smallest order, and upward recursion relations 1141 ! for others (if there are) 1142 ! 1/(2T) 1143 my_i2arg = half/my_arg 1144 ! 1/T 1145 my_iarg = my_i2arg + my_i2arg 1146 ! (0.5-n+j-1)/T == -(2n+1)/(2T), j == 0 1147 div_asymp = -xfloat(strt_order+strt_order+1)*my_i2arg 1148 ! first term in asymptotic series expansions: 1/T 1149 val_asymp = my_i2arg 1150 val_gfun(strt_order,iarg) = val_asymp 1151 ! left terms 1152 do iasymp = 1, max_aterms_gfun 1153 ! plus 1/T 1154 div_asymp = div_asymp + my_iarg 1155 ! times (0.5-n+j-1)/T 1156 val_asymp = val_asymp*div_asymp 1157 val_gfun(strt_order,iarg) = val_gfun(strt_order,iarg) + val_asymp 1158 if ( abs(val_asymp).lt.g1_cut_aux ) exit 1159 end do 1160 ! recorder of (2n+1)/2T 1161 div_asymp = xfloat(strt_order+strt_order-1)*my_i2arg 1162 ! loop over other orders using upward recursion relations 1163 do iord = strt_order, end_order-1 1164 ! plus 1/T 1165 div_asymp = div_asymp + my_iarg 1166 ! perform recursion relation 1167 val_gfun(iord+1,iarg) = my_i2arg - div_asymp*val_gfun(iord,iarg) 1168 end do ! loop over other orders 1169 end do ! loop over number of arguments 1170 ! pop the stack 1171 call xsub_exit 1172 return 11731000 format('GASY ',A,F14.8) 11741010 format('GASY ',A,I6) 1175 end subroutine gen1int_gfun_asymp 1176 1177end module gen1int_rfunction 1178 1179