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!! Tests subroutines in delta_hket.F90. 18!! 19!! 2012-03-16, Bin Gao 20!! * first version 21 22 !> \brief tests subroutines in delta_hket.F90 23 !> \author Bin Gao 24 !> \date 2012-03-16 25 !> \param io_log is the IO unit of log file 26 !> \return test_failed indicates if the test is failed 27 subroutine test_delta_hket(io_log, test_failed) 28 use xkind 29 ! module of HTML test log routines 30 use html_log 31 implicit none 32 integer, intent(in) :: io_log 33 logical, intent(out) :: test_failed 34! parameters of testing primitive integrals 35#include "test_f90/prim_gto.h" 36 integer, parameter :: ORDER_ELEC = 2 !order of electronic derivatives 37 integer orders_hgto_ket(2) !orders of HGTOs on ket center 38 integer orders_geo_pot(2) !orders of geometric derivatives 39 ! assumes the test will pass 40 test_failed = .false. 41 ! dumps information of tests 42 call html_log_real_array(log_text="Coordinates of bra center:", & 43 log_real=COORD_BRA, fmt_real="F16.8", io_log=io_log) 44 call html_log_real_number("Orbital exponent of bra center:", & 45 EXPONENT_BRA, "F16.8", io_log) 46 call html_log_real_array(log_text="Coordinates of ket center:", & 47 log_real=COORD_KET, fmt_real="F16.8", io_log=io_log) 48 call html_log_real_number("Orbital exponent of ket center:", & 49 EXPONENT_KET, "F16.8", io_log) 50 call html_log_int_number("Order of electronic derivatives:", & 51 ORDER_ELEC, "I3", io_log) 52 call html_log_real_number("Scaling constant on Dirac delta function:", & 53 SCAL_CONST, "F16.8", io_log) 54 ! the following tests are for individual cases in \fn(delta_hket) 55 orders_hgto_ket = (/0,0/) 56 orders_geo_pot = (/0,3/) 57 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 58 io_log, test_failed) 59 orders_hgto_ket = (/0,1/) 60 orders_geo_pot = (/0,0/) 61 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 62 io_log, test_failed) 63 orders_hgto_ket = (/0,2/) 64 orders_geo_pot = (/0,0/) 65 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 66 io_log, test_failed) 67 orders_hgto_ket = (/0,3/) 68 orders_geo_pot = (/0,0/) 69 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 70 io_log, test_failed) 71 orders_hgto_ket = (/0,5/) 72 orders_geo_pot = (/0,0/) 73 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 74 io_log, test_failed) 75 orders_hgto_ket = (/0,1/) 76 orders_geo_pot = (/4,4/) 77 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 78 io_log, test_failed) 79 orders_hgto_ket = (/1,1/) 80 orders_geo_pot = (/4,4/) 81 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 82 io_log, test_failed) 83 orders_hgto_ket = (/0,2/) 84 orders_geo_pot = (/4,4/) 85 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 86 io_log, test_failed) 87 orders_hgto_ket = (/0,3/) 88 orders_geo_pot = (/4,4/) 89 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 90 io_log, test_failed) 91 orders_hgto_ket = (/0,4/) 92 orders_geo_pot = (/4,4/) 93 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 94 io_log, test_failed) 95 orders_hgto_ket = (/4,4/) 96 orders_geo_pot = (/4,4/) 97 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 98 io_log, test_failed) 99 orders_hgto_ket = (/0,5/) 100 orders_geo_pot = (/4,4/) 101 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 102 io_log, test_failed) 103 orders_hgto_ket = (/4,5/) 104 orders_geo_pot = (/4,4/) 105 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 106 io_log, test_failed) 107 orders_hgto_ket = (/5,5/) 108 orders_geo_pot = (/4,4/) 109 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 110 io_log, test_failed) 111 orders_hgto_ket = (/5,6/) 112 orders_geo_pot = (/4,4/) 113 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 114 io_log, test_failed) 115 orders_hgto_ket = (/0,1/) 116 orders_geo_pot = (/0,4/) 117 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 118 io_log, test_failed) 119 orders_hgto_ket = (/1,1/) 120 orders_geo_pot = (/1,4/) 121 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 122 io_log, test_failed) 123 orders_hgto_ket = (/0,2/) 124 orders_geo_pot = (/0,4/) 125 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 126 io_log, test_failed) 127 orders_hgto_ket = (/0,3/) 128 orders_geo_pot = (/1,4/) 129 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 130 io_log, test_failed) 131 orders_hgto_ket = (/0,4/) 132 orders_geo_pot = (/0,4/) 133 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 134 io_log, test_failed) 135 orders_hgto_ket = (/4,4/) 136 orders_geo_pot = (/1,4/) 137 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 138 io_log, test_failed) 139 orders_hgto_ket = (/0,5/) 140 orders_geo_pot = (/0,4/) 141 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 142 io_log, test_failed) 143 orders_hgto_ket = (/4,5/) 144 orders_geo_pot = (/1,4/) 145 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 146 io_log, test_failed) 147 orders_hgto_ket = (/5,5/) 148 orders_geo_pot = (/0,4/) 149 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 150 io_log, test_failed) 151 orders_hgto_ket = (/5,6/) 152 orders_geo_pot = (/1,4/) 153 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 154 io_log, test_failed) 155 orders_hgto_ket = (/5,6/) 156 orders_geo_pot = (/1,2/) 157 call sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, ORDER_ELEC, & 158 io_log, test_failed) 159 return 160 end subroutine test_delta_hket 161 162 !> \brief tests subroutines in delta_hket.F90 with specific orders of HGTOs on 163 !> ket center and geometric derivatives on Dirac delta function 164 !> \author Bin Gao 165 !> \date 2012-03-16 166 !> \param orders_hgto_ket contains the orders of HGTOs on ket center 167 !> \param orders_geo_pot contains the orders of geometric derivatives on Dirac delta function 168 !> \param order_elec is the order of electronic derivatives 169 !> \param io_log is the IO unit of log file 170 !> \return test_failed indicates if the test is failed 171 subroutine sub_test_delta_hket(orders_hgto_ket, orders_geo_pot, order_elec, & 172 io_log, test_failed) 173 use xkind 174 ! module of HTML test log routines 175 use html_log 176 ! recursive functions of Dirac delta function integrals 177 use recur_delta 178 implicit none 179 integer, intent(in) :: orders_hgto_ket(2) 180 integer, intent(in) :: orders_geo_pot(2) 181 integer, intent(in) :: order_elec 182 integer, intent(in) :: io_log 183 logical, intent(inout) :: test_failed 184! parameters of testing primitive integrals 185#include "test_f90/prim_gto.h" 186 integer orders_hgto_bra(2) !orders of HGTOs on bra center 187 integer order_hbra(3) !orders of xyz components of HGTOs on bra center 188 integer max_dim_geo !maximum dimension of geometric derivatives on ... 189 !Dirac delta function 190 integer min_delta_geom !minimum order of geometric derivatives on Dirac ... 191 !delta function from \fn(delta_geom) 192 integer dim_geo_pot !dimension of integrals from \fn(delta_geom) 193 real(REALK), allocatable :: geo_pot_pints(:) !integrals from \fn(delta_geom) 194 integer min_geo_hbra !minimum order of geometric derivatives on Dirac ... 195 !delta function after recovering HGTOs on bra center 196 integer dim_hgto_bra !dimension of HGTOs on bra center 197 integer dim_geo_hbra !dimension of geometric derivatives on Dirac delta ... 198 !function after recovering HGTOs on bra center 199 real(REALK), allocatable :: hbra_pints(:,:) !integrals after recovering HGTOs on bra center 200 integer dim_hgto_ket !dimension of HGTOs on ket center 201 integer dim_geo_hket !dimension of geometric derivatives on Dirac delta ... 202 !function after recovering HGTOs on ket center 203 real(REALK), allocatable :: hket_pints(:,:,:) !integrals after recovering HGTOs on ket center 204 integer order_geo !order of geometric derivatives on Dirac delta function 205 integer x_geo, y_geo, z_geo !orders of xyz components of geometric derivatives 206 integer order_hket !order of HGTOs on ket center 207 integer x_hket, y_hket, z_hket !orders of xyz components of HGTOs on ket center 208 real(REALK) recur_pint !results from recursive function \fn(rec_delta_hket) 209 integer addr_hket, addr_geo !addresses of \var(hket_pints) 210 integer ierr !error information 211 real(REALK) begin_time !begin of CPU time 212 real(REALK) end_time !end of CPU time 213 logical different !if result from \fn(delta_hket) is different from reference 214 ! dumps information of tests 215 call html_log_int_array(log_text="Orders of HGTOs on ket center:", & 216 log_int=orders_hgto_ket, fmt_int="I3", io_log=io_log, & 217 separation="→") 218 call html_log_int_array(log_text="Orders of geometric derivatives:", & 219 log_int=orders_geo_pot, fmt_int="I3", io_log=io_log, & 220 separation="→") 221 ! sets the orders of HGTOs on bra center 222 orders_hgto_bra = (/0,0/) 223 order_hbra = (/0,0,0/) 224 ! gets the begin time 225 call xtimer_set(begin_time) 226 ! sets the minimum order of geometric derivatives on Dirac delta function 227 ! after recovering HGTOs on bra center 228 min_geo_hbra = max(0,orders_geo_pot(1)-orders_hgto_ket(2)) 229 ! sets the minimum order of geometric derivatives on Dirac delta function from \fn(delta_geom) 230 min_delta_geom = max(0,min_geo_hbra-orders_hgto_bra(2)) 231 ! sets the maximum dimension of geometric derivatives on Dirac delta function 232 max_dim_geo = (orders_geo_pot(2)+1)*(orders_geo_pot(2)+2)*(orders_geo_pot(2)+3)/6 233 ! allocates memory for integrals from \fn(delta_geom) 234 if (min_delta_geom==0) then 235 dim_geo_pot = max_dim_geo 236 else 237 dim_geo_pot = max_dim_geo & 238 - min_delta_geom*(min_delta_geom+1)*(min_delta_geom+2)/6 239 end if 240 allocate(geo_pot_pints(dim_geo_pot), stat=ierr) 241 if (ierr/=0) then 242 call html_log_int_number(log_text="Failed to allocate geo_pot_pints:", & 243 log_int=dim_geo_pot, fmt_int="I12", & 244 io_log=io_log, font_color="red") 245 test_failed = .true. 246 return 247 end if 248 ! recovers the geometric derivatives on Dirac delta function using \fn(delta_geom) 249 call delta_geom(COORD_BRA, EXPONENT_BRA, COORD_KET, EXPONENT_KET, & 250 DELORG, SCAL_CONST, (/min_delta_geom,orders_geo_pot(2)/), & 251 order_elec, dim_geo_pot, geo_pot_pints) 252 ! allocates memory for integrals from \fn(delta_hket) on bra center 253 if (orders_hgto_bra(1)==0) then 254 dim_hgto_bra = (orders_hgto_bra(2)+1)*(orders_hgto_bra(2)+2) & 255 * (orders_hgto_bra(2)+3)/6 256 else 257 dim_hgto_bra = ((orders_hgto_bra(2)+1)*(orders_hgto_bra(2)+2) & 258 * (orders_hgto_bra(2)+3) & 259 - orders_hgto_bra(1)*(orders_hgto_bra(1)+1) & 260 * (orders_hgto_bra(1)+2))/6 261 end if 262 if (min_geo_hbra==0) then 263 dim_geo_hbra = max_dim_geo 264 else 265 dim_geo_hbra = max_dim_geo-min_geo_hbra*(min_geo_hbra+1)*(min_geo_hbra+2)/6 266 end if 267 allocate(hbra_pints(dim_hgto_bra,dim_geo_hbra), stat=ierr) 268 if (ierr/=0) then 269 call html_log_int_number(log_text="Failed to allocate hbra_pints:", & 270 log_int=dim_hgto_bra*dim_geo_hbra, & 271 fmt_int="I12", io_log=io_log, font_color="red") 272 deallocate(geo_pot_pints) 273 test_failed = .true. 274 return 275 end if 276 ! recovers the HGTOs on bra center using \fn(delta_hket) 277 call delta_hket(orders_hgto_bra, (/min_geo_hbra,orders_geo_pot(2)/), & 278 COORD_BRA, EXPONENT_BRA, DELORG, 1, dim_geo_pot, & 279 geo_pot_pints, dim_hgto_bra, dim_geo_hbra, hbra_pints) 280 deallocate(geo_pot_pints) 281 ! allocates memory for integrals from \fn(delta_hket) on ket center 282 if (orders_hgto_ket(1)==0) then 283 dim_hgto_ket = (orders_hgto_ket(2)+1)*(orders_hgto_ket(2)+2) & 284 * (orders_hgto_ket(2)+3)/6 285 else 286 dim_hgto_ket = ((orders_hgto_ket(2)+1)*(orders_hgto_ket(2)+2) & 287 * (orders_hgto_ket(2)+3) & 288 - orders_hgto_ket(1)*(orders_hgto_ket(1)+1) & 289 * (orders_hgto_ket(1)+2))/6 290 end if 291 if (orders_geo_pot(1)==0) then 292 dim_geo_hket = max_dim_geo 293 else 294 dim_geo_hket = max_dim_geo & 295 - orders_geo_pot(1)*(orders_geo_pot(1)+1)*(orders_geo_pot(1)+2)/6 296 end if 297 allocate(hket_pints(dim_hgto_bra,dim_hgto_ket,dim_geo_hket), stat=ierr) 298 if (ierr/=0) then 299 call html_log_int_number(log_text="Failed to allocate hket_pints:", & 300 log_int=dim_hgto_bra*dim_hgto_ket*dim_geo_hket, & 301 fmt_int="I12", io_log=io_log, font_color="red") 302 deallocate(hbra_pints) 303 test_failed = .true. 304 return 305 end if 306 ! recovers the HGTOs on ket center using \fn(delta_hket) 307 call delta_hket(orders_hgto_ket, orders_geo_pot, COORD_KET, EXPONENT_KET, & 308 DELORG, dim_hgto_bra, dim_geo_hbra, hbra_pints, & 309 dim_hgto_ket, dim_geo_hket, hket_pints) 310 deallocate(hbra_pints) 311 ! gets the end time 312 call xtimer_set(end_time) 313 ! prints the CPU elapsed time 314 call html_log_real_number(log_text="Time (s) used for delta_geom and delta_hket:", & 315 log_real=end_time-begin_time, fmt_real="F10.4", & 316 io_log=io_log, font_color="blue") 317 ! resets the begin time 318 call xtimer_set(begin_time) 319 ! sets the data used for the module \fn(recur_delta) 320 call recur_delta_create(COORD_BRA, EXPONENT_BRA, COORD_KET, EXPONENT_KET, & 321 DELORG, DIPORG, SCAL_CONST, order_elec) 322 ! initializes the address of integrals from \fn(delta_hket) 323 addr_geo = 0 324 ! loops over the orders of geometric derivatives 325 do order_geo = orders_geo_pot(1), orders_geo_pot(2) 326 do z_geo = 0, order_geo 327 do y_geo = 0, order_geo-z_geo 328 x_geo = order_geo-(z_geo+y_geo) 329 addr_geo = addr_geo+1 330 ! initializes the address of integrals from \fn(delta_hket) 331 addr_hket = 0 332 do order_hket = orders_hgto_ket(1), orders_hgto_ket(2) 333 do z_hket = 0, order_hket 334 do y_hket = 0, order_hket-z_hket 335 x_hket = order_hket-(z_hket+y_hket) 336 ! integral from \fn(rec_delta_hket) 337 recur_pint = recur_delta_hket(order_hbra, (/x_hket,y_hket,z_hket/), & 338 (/x_geo,y_geo,z_geo/)) 339 addr_hket = addr_hket+1 340 call check_difference(recur_pint, hket_pints(1,addr_hket,addr_geo), & 341 different) 342 if (different) then 343 call html_log_int_array(log_text="Orders of HGTOs on ket center:", & 344 log_int=(/x_hket,y_hket,z_hket/), & 345 fmt_int="I3", io_log=io_log) 346 call html_log_int_array(log_text="Orders of geometric derivatives:", & 347 log_int=(/x_geo,y_geo,z_geo/), fmt_int="I3", & 348 io_log=io_log) 349 call html_log_real_number( & 350 log_text="Reference from recur_delta_hket:", & 351 log_real=recur_pint, fmt_real="Es20.12", & 352 io_log=io_log, font_color="blue") 353 call html_log_real_number(log_text="Result from delta_hket:", & 354 log_real=hket_pints(1,addr_hket,addr_geo), & 355 fmt_real="Es20.12", io_log=io_log, & 356 font_color="red") 357 test_failed = .true. 358 end if 359 end do 360 end do 361 end do 362 end do 363 end do 364 end do 365 ! gets the end time 366 call xtimer_set(end_time) 367 ! prints the CPU elapsed time 368 call html_log_real_number(log_text="Time (s) used for recur_delta_hket:", & 369 log_real=end_time-begin_time, fmt_real="F10.4", & 370 io_log=io_log, font_color="blue") 371 ! cleans 372 deallocate(hket_pints) 373 return 374 end subroutine sub_test_delta_hket 375