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