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 provides the Fortran 90 module of geometric derivatives. 18!! 19!! 2012-03-20, Bin Gao: 20!! * first version 21 22#include "stdout.h" 23#include "err_info.h" 24 25!> \brief Fortran 90 module of geometric derivatives 26!> \author Bin Gao 27!> \date 2012-03-20 28module gen1int_geom 29 30 use xkind 31 use london_ao 32 implicit none 33 34 ! information of N-ary tree for geometric derivatives 35 type, public :: nary_tree_t 36 private 37 ! number of atoms 38 integer :: num_atoms = 0 39 ! indices of atoms to be used as differentiated centers 40 integer, allocatable :: idx_atoms(:) 41 ! order of geometric derivatives 42 integer :: order_geo = 0 43 ! maximum number of differentiated centers 44 integer :: max_num_cent = 0 45 ! total number of different allowed paths in N-ary tree 46 integer :: num_paths = 0 47 ! number of unique geometric derivatives in the N-ary tree 48 integer :: num_unique_geo = 0 49 ! index of current path 50 integer :: idx_path = 0 51 ! height of atom to visit 52 integer visit_height 53 ! indices of the selected atom nodes 54 integer, allocatable :: idx_node(:) 55 ! weights of the selected atom nodes 56 integer, allocatable :: wt_node(:) 57 ! indices of the generated differentiated centers 58 integer, allocatable :: idx_cent(:) 59 ! orders of derivatives of the differentiated centers 60 integer, allocatable :: order_cent(:) 61 ! offset of unique derivatives of current path 62 integer :: path_offset = 0 63 ! number of unique derivatives of current path 64 integer :: path_num_unique = 0 65 ! number of redundant derivatives of current path 66 integer :: path_num_redunt = 0 67 end type nary_tree_t 68 69 public :: NaryTreeCreate 70 public :: NaryTreeSetAtoms 71 public :: NaryTreeGetNumAtoms 72 public :: NaryTreeGetOrder 73 public :: NaryTreeGetMaxNumCenters 74 public :: NaryTreeGetNumPaths 75 public :: NaryTreeGetNumGeo 76 public :: NaryTreeSearch 77 public :: NaryTreePathGetIndex 78 public :: NaryTreePathGetNumCenters 79 public :: NaryTreePathGetIdxCent 80 public :: NaryTreePathGetOrderCent 81 public :: NaryTreePathGetOffset 82 public :: NaryTreePathGetNumUnique 83 public :: NaryTreePathGetNumRedunt 84 public :: NaryTreePathGetReduntList 85 public :: NaryTreePathSetReduntExpt 86 !-public :: NaryTreePathReorderInts 87 public :: NaryTreeView 88 public :: NaryTreeDestroy 89 90 public :: IntGetCARMOM 91 public :: IntGetNUCPOT 92 public :: IntGetGAUPOT 93 public :: IntGetODST 94 95 contains 96 97 !> \brief initializes the information of N-ary tree for geometric derivatives, and 98 !> returns the total number of different paths, and generates the first path 99 !> \author Bin Gao 100 !> \date 2011-12-12 101 !> \param num_atoms is the number of atoms 102 !> \param order_geo is the order of geometric derivatives 103 !> \param max_num_cent is the maximum number of differentiated centers 104 !> \return nary_tree contains the information of N-ary tree for geometric derivatives 105 !> \return info_geom (==ERR_INFO) indicates the N-ary tree is not successfully created 106 !> \return num_paths is the total number of different paths 107 !> \return path_num_unique is the number of unique derivatives of the first path 108 !> \return path_num_redunt is the number of redundant geometric derivatives of the first path 109 subroutine NaryTreeCreate(num_atoms, order_geo, max_num_cent, nary_tree, info_geom, & 110 num_paths, path_num_unique, path_num_redunt) 111 integer, intent(in) :: num_atoms 112 integer, intent(in) :: order_geo 113 integer, intent(in) :: max_num_cent 114 type(nary_tree_t), intent(inout) :: nary_tree 115 integer, intent(out) :: info_geom 116 integer, optional, intent(out) :: num_paths 117 integer, optional, intent(out) :: path_num_unique 118 integer, optional, intent(out) :: path_num_redunt 119 integer num_cent !number of differentiated centers 120#if defined(XTIME) 121 real(REALK) curr_time !current CPU time 122 ! sets current CPU time 123 call xtimer_set(curr_time) 124#endif 125 if (num_atoms<1 .or. order_geo<0) then 126 info_geom = ERR_INFO 127 else 128 info_geom = 0 129 nary_tree%num_atoms = num_atoms 130 nary_tree%order_geo = order_geo 131 ! no geometric derivatives 132 if (nary_tree%order_geo==0) then 133 nary_tree%max_num_cent = 0 134 allocate(nary_tree%idx_node(0:0), stat=info_geom) 135 if (info_geom/=0) return 136 allocate(nary_tree%wt_node(0:0), stat=info_geom) 137 if (info_geom/=0) then 138 deallocate(nary_tree%idx_node) 139 return 140 end if 141 allocate(nary_tree%idx_cent(1), stat=info_geom) 142 if (info_geom/=0) then 143 deallocate(nary_tree%idx_node) 144 deallocate(nary_tree%wt_node) 145 return 146 end if 147 allocate(nary_tree%order_cent(1), stat=info_geom) 148 if (info_geom/=0) then 149 deallocate(nary_tree%idx_node) 150 deallocate(nary_tree%wt_node) 151 deallocate(nary_tree%idx_cent) 152 return 153 end if 154 nary_tree%num_paths = 1 155 if (present(num_paths)) num_paths = nary_tree%num_paths 156 nary_tree%idx_path = 1 157 nary_tree%visit_height = 0 158 nary_tree%idx_node = 0 159 nary_tree%wt_node = 0 160 nary_tree%idx_cent = 0 161 nary_tree%order_cent = 0 162 nary_tree%path_offset = 0 163 nary_tree%path_num_unique = 1 164 nary_tree%path_num_redunt = 1 165 nary_tree%num_unique_geo = 1 166 if (present(path_num_unique)) path_num_unique = nary_tree%path_num_unique 167 if (present(path_num_redunt)) path_num_redunt = nary_tree%path_num_redunt 168 ! having geometric derivatives 169 else 170 nary_tree%max_num_cent = max_num_cent 171 ! if the number of atoms is less than the number of differentiated centers 172 nary_tree%max_num_cent = min(nary_tree%num_atoms,nary_tree%max_num_cent) 173 if (nary_tree%max_num_cent<1) then 174 info_geom = ERR_INFO 175 return 176 end if 177 allocate(nary_tree%idx_node(order_geo), stat=info_geom) 178 if (info_geom/=0) return 179 allocate(nary_tree%wt_node(order_geo), stat=info_geom) 180 if (info_geom/=0) then 181 deallocate(nary_tree%idx_node) 182 return 183 end if 184 allocate(nary_tree%idx_cent(nary_tree%max_num_cent), stat=info_geom) 185 if (info_geom/=0) then 186 deallocate(nary_tree%idx_node) 187 deallocate(nary_tree%wt_node) 188 return 189 end if 190 allocate(nary_tree%order_cent(nary_tree%max_num_cent), stat=info_geom) 191 if (info_geom/=0) then 192 deallocate(nary_tree%idx_node) 193 deallocate(nary_tree%wt_node) 194 deallocate(nary_tree%idx_cent) 195 return 196 end if 197 ! returns the total number of different paths, and generates the first path 198 call geom_total_tree_init(num_atoms, order_geo, & 199 nary_tree%max_num_cent, & 200 nary_tree%num_paths, & 201 nary_tree%visit_height, & 202 nary_tree%idx_node, & 203 nary_tree%wt_node, & 204 nary_tree%idx_cent, & 205 nary_tree%order_cent, & 206 nary_tree%path_num_unique) 207 ! gets the number of redundant derivatives of current path 208 num_cent = nary_tree%wt_node(order_geo) 209 call geom_total_num_redunt(num_cent, nary_tree%order_cent(1:num_cent), & 210 nary_tree%path_num_redunt) 211 if (present(path_num_unique)) path_num_unique = nary_tree%path_num_unique 212 if (present(path_num_redunt)) path_num_redunt = nary_tree%path_num_redunt 213 if (present(num_paths)) num_paths = nary_tree%num_paths 214 nary_tree%idx_path = 1 215 ! gets the number of unique geometric derivatives in the N-ary tree 216 call geom_total_num_derv(order_geo, nary_tree%max_num_cent, num_atoms, & 217 nary_tree%num_unique_geo) 218 end if 219 end if 220#if defined(XTIME) 221 ! prints the CPU elapsed time 222 call xtimer_view(curr_time, "NaryTreeCreate", STDOUT) 223#endif 224 end subroutine NaryTreeCreate 225 226 !> \brief sets the indices of atoms to be used as differentiated centers 227 !> \author Bin Gao 228 !> \date 2012-05-09 229 !> \param num_atoms is the number of atoms 230 !> \param idx_atoms contains the indices of atoms to be used as differentiated centers 231 !> \return nary_tree contains the information of N-ary tree for geometric derivatives 232 !> \return info_geom (==ERR_INFO) indicates the N-ary tree is not successfully created 233 subroutine NaryTreeSetAtoms(num_atoms, idx_atoms, nary_tree, info_geom) 234 integer, intent(in) :: num_atoms 235 integer, intent(in) :: idx_atoms(num_atoms) 236 type(nary_tree_t), intent(inout) :: nary_tree 237 integer, intent(out) :: info_geom 238#if defined(XTIME) 239 real(REALK) curr_time !current CPU time 240 ! sets current CPU time 241 call xtimer_set(curr_time) 242#endif 243 if (num_atoms==nary_tree%num_atoms) then 244 allocate(nary_tree%idx_atoms(num_atoms), stat=info_geom) 245 if (info_geom==0) then 246 nary_tree%idx_atoms = idx_atoms 247 end if 248 else 249 info_geom = num_atoms 250 end if 251#if defined(XTIME) 252 ! prints the CPU elapsed time 253 call xtimer_view(curr_time, "NaryTreeSetAtoms", STDOUT) 254#endif 255 end subroutine NaryTreeSetAtoms 256 257 !> \brief gets the number of atoms for a given N-ary tree 258 !> \author Bin Gao 259 !> \date 2012-05-09 260 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 261 !> \return num_atoms is the number of atoms 262 subroutine NaryTreeGetNumAtoms(nary_tree, num_atoms) 263 type(nary_tree_t), intent(in) :: nary_tree 264 integer, intent(out) :: num_atoms 265#if defined(XTIME) 266 real(REALK) curr_time !current CPU time 267 ! sets current CPU time 268 call xtimer_set(curr_time) 269#endif 270 num_atoms = nary_tree%num_atoms 271#if defined(XTIME) 272 ! prints the CPU elapsed time 273 call xtimer_view(curr_time, "NaryTreeGetNumAtoms", STDOUT) 274#endif 275 end subroutine NaryTreeGetNumAtoms 276 277 !> \brief gets the order of geometric derivatives for a given N-ary tree 278 !> \author Bin Gao 279 !> \date 2012-05-09 280 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 281 !> \return order_geo is the order of geometric derivatives 282 subroutine NaryTreeGetOrder(nary_tree, order_geo) 283 type(nary_tree_t), intent(in) :: nary_tree 284 integer, intent(out) :: order_geo 285#if defined(XTIME) 286 real(REALK) curr_time !current CPU time 287 ! sets current CPU time 288 call xtimer_set(curr_time) 289#endif 290 order_geo = nary_tree%order_geo 291#if defined(XTIME) 292 ! prints the CPU elapsed time 293 call xtimer_view(curr_time, "NaryTreeGetOrder", STDOUT) 294#endif 295 end subroutine NaryTreeGetOrder 296 297 !> \brief gets the maximum number of differentiated centers for a given N-ary tree 298 !> \author Bin Gao 299 !> \date 2012-05-09 300 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 301 !> \return max_num_cent is the maximum number of differentiated centers 302 subroutine NaryTreeGetMaxNumCenters(nary_tree, max_num_cent) 303 type(nary_tree_t), intent(in) :: nary_tree 304 integer, intent(out) :: max_num_cent 305#if defined(XTIME) 306 real(REALK) curr_time !current CPU time 307 ! sets current CPU time 308 call xtimer_set(curr_time) 309#endif 310 max_num_cent = nary_tree%max_num_cent 311#if defined(XTIME) 312 ! prints the CPU elapsed time 313 call xtimer_view(curr_time, "NaryTreeGetMaxNumCenters", STDOUT) 314#endif 315 end subroutine NaryTreeGetMaxNumCenters 316 317 !> \brief gets the total number of different allowed paths in N-ary tree 318 !> \author Bin Gao 319 !> \date 2012-03-09 320 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 321 !> \return num_paths is the total number of different allowed paths in N-ary tree 322 subroutine NaryTreeGetNumPaths(nary_tree, num_paths) 323 type(nary_tree_t), intent(in) :: nary_tree 324 integer, intent(out) :: num_paths 325#if defined(XTIME) 326 real(REALK) curr_time !current CPU time 327 ! sets current CPU time 328 call xtimer_set(curr_time) 329#endif 330 ! checks if the N-ary tree is created 331 if (nary_tree%idx_path==0) then 332 call error_stop("NaryTreeGetNumPaths", "N-ary tree is not created", 0) 333 else 334 num_paths = nary_tree%num_paths 335 end if 336#if defined(XTIME) 337 ! prints the CPU elapsed time 338 call xtimer_view(curr_time, "NaryTreeGetNumPaths", STDOUT) 339#endif 340 end subroutine NaryTreeGetNumPaths 341 342 !> \brief gets the number of unique geometric derivatives in the N-ary tree 343 !> \author Bin Gao 344 !> \date 2012-01-12 345 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 346 !> \return num_unique_geo is the number of unique geometric derivatives in the N-ary tree 347 subroutine NaryTreeGetNumGeo(nary_tree, num_unique_geo) 348 type(nary_tree_t), intent(in) :: nary_tree 349 integer, intent(out) :: num_unique_geo 350#if defined(XTIME) 351 real(REALK) curr_time !current CPU time 352 ! sets current CPU time 353 call xtimer_set(curr_time) 354#endif 355 ! checks if the N-ary tree is created 356 if (nary_tree%idx_path==0) then 357 call error_stop("NaryTreeGetNumGeo", "N-ary tree is not created", 0) 358 else 359 num_unique_geo = nary_tree%num_unique_geo 360 end if 361#if defined(XTIME) 362 ! prints the CPU elapsed time 363 call xtimer_view(curr_time, "NaryTreeGetNumGeo", STDOUT) 364#endif 365 end subroutine NaryTreeGetNumGeo 366 367 !> \brief searches for the next satisfied path from a given path, could be called recursively 368 !> \author Bin Gao 369 !> \date 2011-12-12 370 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 371 !> \return path_num_unique is the number of unique derivatives of current path 372 !> \return path_num_redunt is the number of redundant derivatives of current path 373 subroutine NaryTreeSearch(nary_tree, path_num_unique, path_num_redunt) 374 type(nary_tree_t), intent(inout) :: nary_tree 375 integer, optional, intent(out) :: path_num_unique 376 integer, optional, intent(out) :: path_num_redunt 377 integer num_cent !number of differentiated centers 378#if defined(XTIME) 379 real(REALK) curr_time !current CPU time 380 ! sets current CPU time 381 call xtimer_set(curr_time) 382#endif 383 ! checks if the N-ary tree is created 384 if (nary_tree%idx_path==0) then 385 call error_stop("NaryTreeSearch", "N-ary tree is not created", 0) 386 else if (nary_tree%order_geo==0) then 387 if (present(path_num_unique)) path_num_unique = nary_tree%path_num_unique 388 if (present(path_num_redunt)) path_num_redunt = nary_tree%path_num_redunt 389 else 390 ! updates the offset of unique derivatives of current path 391 if (nary_tree%idx_path==nary_tree%num_paths) then 392 nary_tree%path_offset = 0 393 nary_tree%idx_path = 0 394 else 395 nary_tree%path_offset = nary_tree%path_offset+nary_tree%path_num_unique 396 end if 397 ! searches for the next satisfied path from the given path 398 call geom_total_tree_search(nary_tree%num_atoms, nary_tree%order_geo, & 399 nary_tree%max_num_cent, nary_tree%visit_height, & 400 nary_tree%idx_node, nary_tree%wt_node, & 401 nary_tree%idx_cent, nary_tree%order_cent, & 402 nary_tree%path_num_unique) 403 ! gets the number of redundant derivatives of current path 404 num_cent = nary_tree%wt_node(nary_tree%order_geo) 405 call geom_total_num_redunt(num_cent, nary_tree%order_cent(1:num_cent), & 406 nary_tree%path_num_redunt) 407 if (present(path_num_unique)) path_num_unique = nary_tree%path_num_unique 408 if (present(path_num_redunt)) path_num_redunt = nary_tree%path_num_redunt 409 nary_tree%idx_path = nary_tree%idx_path+1 410 end if 411#if defined(XTIME) 412 ! prints the CPU elapsed time 413 call xtimer_view(curr_time, "NaryTreeSearch", STDOUT) 414#endif 415 end subroutine NaryTreeSearch 416 417 !> \brief gets the index of current path 418 !> \author Bin Gao 419 !> \date 2012-03-09 420 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 421 !> \return idx_path is the index of current path 422 subroutine NaryTreePathGetIndex(nary_tree, idx_path) 423 type(nary_tree_t), intent(in) :: nary_tree 424 integer, intent(out) :: idx_path 425#if defined(XTIME) 426 real(REALK) curr_time !current CPU time 427 ! sets current CPU time 428 call xtimer_set(curr_time) 429#endif 430 ! checks if the N-ary tree is created 431 if (nary_tree%idx_path==0) then 432 call error_stop("NaryTreePathGetIndex", "N-ary tree is not created", 0) 433 else 434 idx_path = nary_tree%idx_path 435 end if 436#if defined(XTIME) 437 ! prints the CPU elapsed time 438 call xtimer_view(curr_time, "NaryTreePathGetIndex", STDOUT) 439#endif 440 end subroutine NaryTreePathGetIndex 441 442 !> \brief returns the number of differentiated centers of current path for given N-ary tree 443 !> \author Bin Gao 444 !> \date 2012-03-20 445 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 446 !> \return num_centers is the number of differentiated centers of current path 447 subroutine NaryTreePathGetNumCenters(nary_tree, num_centers) 448 type(nary_tree_t), intent(in) :: nary_tree 449 integer, intent(out) :: num_centers 450#if defined(XTIME) 451 real(REALK) curr_time !current CPU time 452 ! sets current CPU time 453 call xtimer_set(curr_time) 454#endif 455 ! checks if the N-ary tree is created 456 if (nary_tree%idx_path==0) then 457 call error_stop("NaryTreePathGetNumCenters", "N-ary tree is not created", 0) 458 else 459 num_centers = nary_tree%wt_node(nary_tree%order_geo) 460 end if 461#if defined(XTIME) 462 ! prints the CPU elapsed time 463 call xtimer_view(curr_time, "NaryTreePathGetNumCenters", STDOUT) 464#endif 465 end subroutine NaryTreePathGetNumCenters 466 467 !> \brief returns the indices of differentiated centers of current path for given N-ary tree 468 !> \author Bin Gao 469 !> \date 2013-05-07 470 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 471 !> \return idx_cent contains the indices of differentiated centers of current path 472 subroutine NaryTreePathGetIdxCent(nary_tree, idx_cent) 473 type(nary_tree_t), intent(in) :: nary_tree 474 integer, intent(out) :: idx_cent(:) 475 integer num_centers !number of differentiated centers of current path 476#if defined(XTIME) 477 real(REALK) curr_time !current CPU time 478 ! sets current CPU time 479 call xtimer_set(curr_time) 480#endif 481 ! checks if the N-ary tree is created 482 if (nary_tree%idx_path==0) then 483 call error_stop("NaryTreePathGetIdxCent", "N-ary tree is not created", 0) 484 else 485 num_centers = nary_tree%wt_node(nary_tree%order_geo) 486 if (size(idx_cent)<num_centers) then 487 call error_stop("NaryTreePathGetIdxCent", & 488 "size of idx_cent not enough", & 489 num_centers) 490 else 491 idx_cent(1:num_centers) = nary_tree%idx_cent(1:num_centers) 492 end if 493 end if 494#if defined(XTIME) 495 ! prints the CPU elapsed time 496 call xtimer_view(curr_time, "NaryTreePathGetIdxCent", STDOUT) 497#endif 498 end subroutine NaryTreePathGetIdxCent 499 500 !> \brief returns the orders of differentiated centers of current path for given N-ary tree 501 !> \author Bin Gao 502 !> \date 2013-05-07 503 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 504 !> \return order_cent contains the orders of differentiated centers of current path 505 subroutine NaryTreePathGetOrderCent(nary_tree, order_cent) 506 type(nary_tree_t), intent(in) :: nary_tree 507 integer, intent(out) :: order_cent(:) 508 integer num_centers !number of differentiated centers of current path 509#if defined(XTIME) 510 real(REALK) curr_time !current CPU time 511 ! sets current CPU time 512 call xtimer_set(curr_time) 513#endif 514 ! checks if the N-ary tree is created 515 if (nary_tree%idx_path==0) then 516 call error_stop("NaryTreePathGetOrderCent", "N-ary tree is not created", 0) 517 else 518 num_centers = nary_tree%wt_node(nary_tree%order_geo) 519 if (size(order_cent)<num_centers) then 520 call error_stop("NaryTreePathGetOrderCent", & 521 "size of order_cent not enough", & 522 num_centers) 523 else 524 order_cent(1:num_centers) = nary_tree%order_cent(1:num_centers) 525 end if 526 end if 527#if defined(XTIME) 528 ! prints the CPU elapsed time 529 call xtimer_view(curr_time, "NaryTreePathGetOrderCent", STDOUT) 530#endif 531 end subroutine NaryTreePathGetOrderCent 532 533 !> \brief returns the offset of unique derivatives of current path 534 !> \author Bin Gao 535 !> \date 2012-01-12 536 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 537 !> \return path_offset is the offset of unique derivatives of current path 538 subroutine NaryTreePathGetOffset(nary_tree, path_offset) 539 type(nary_tree_t), intent(in) :: nary_tree 540 integer, intent(out) :: path_offset 541#if defined(XTIME) 542 real(REALK) curr_time !current CPU time 543 ! sets current CPU time 544 call xtimer_set(curr_time) 545#endif 546 ! checks if the N-ary tree is created 547 if (nary_tree%idx_path==0) then 548 call error_stop("NaryTreePathGetOffset", "N-ary tree is not created", 0) 549 else 550 path_offset = nary_tree%path_offset 551 end if 552#if defined(XTIME) 553 ! prints the CPU elapsed time 554 call xtimer_view(curr_time, "NaryTreePathGetOffset", STDOUT) 555#endif 556 end subroutine NaryTreePathGetOffset 557 558 !> \brief returns the number of unique derivatives of current path 559 !> \author Bin Gao 560 !> \date 2012-01-12 561 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 562 !> \return path_num_unique is the number of unique derivatives of current path 563 subroutine NaryTreePathGetNumUnique(nary_tree, path_num_unique) 564 type(nary_tree_t), intent(in) :: nary_tree 565 integer, intent(out) :: path_num_unique 566#if defined(XTIME) 567 real(REALK) curr_time !current CPU time 568 ! sets current CPU time 569 call xtimer_set(curr_time) 570#endif 571 ! checks if the N-ary tree is created 572 if (nary_tree%idx_path==0) then 573 call error_stop("NaryTreePathGetNumUnique", "N-ary tree is not created", 0) 574 else 575 path_num_unique = nary_tree%path_num_unique 576 end if 577#if defined(XTIME) 578 ! prints the CPU elapsed time 579 call xtimer_view(curr_time, "NaryTreePathGetNumUnique", STDOUT) 580#endif 581 end subroutine NaryTreePathGetNumUnique 582 583 !> \brief returns the number of redundant derivatives of current path 584 !> \author Bin Gao 585 !> \date 2012-01-12 586 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 587 !> \return path_num_redunt is the number of redundant derivatives of current path 588 subroutine NaryTreePathGetNumRedunt(nary_tree, path_num_redunt) 589 type(nary_tree_t), intent(in) :: nary_tree 590 integer, intent(out) :: path_num_redunt 591#if defined(XTIME) 592 real(REALK) curr_time !current CPU time 593 ! sets current CPU time 594 call xtimer_set(curr_time) 595#endif 596 ! checks if the N-ary tree is created 597 if (nary_tree%idx_path==0) then 598 call error_stop("NaryTreePathGetNumRedunt", "N-ary tree is not created", 0) 599 else 600 path_num_redunt = nary_tree%path_num_redunt 601 end if 602#if defined(XTIME) 603 ! prints the CPU elapsed time 604 call xtimer_view(curr_time, "NaryTreePathGetNumRedunt", STDOUT) 605#endif 606 end subroutine NaryTreePathGetNumRedunt 607 608 !> \brief gets the list addresses of redundant derivatives of current path 609 !> \author Bin Gao 610 !> \date 2011-12-12 611 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 612 !> \return redunt_list contains the list addresses of redundant derivatives of current path 613 subroutine NaryTreePathGetReduntList(nary_tree, redunt_list) 614 type(nary_tree_t), intent(in) :: nary_tree 615 integer, intent(out) :: redunt_list(:,:) 616 integer path_num_redunt !number of redundant derivatives of current path 617 integer num_cent !number of differentiated centers 618#if defined(XTIME) 619 real(REALK) curr_time !current CPU time 620 ! sets current CPU time 621 call xtimer_set(curr_time) 622#endif 623 ! checks if the N-ary tree is created 624 if (nary_tree%idx_path==0) then 625 call error_stop("NaryTreePathGetReduntList", "N-ary tree is not created", 0) 626 else if (nary_tree%order_geo==0) then 627 redunt_list = 1 628 else 629 ! checks the sizes of the list addresses 630 if (size(redunt_list,1)/=2) & 631 call error_stop("NaryTreePathGetReduntList", & 632 "incorrect size of variable redunt_list", & 633 size(redunt_list,1)) 634 path_num_redunt = size(redunt_list,2) 635 if (path_num_redunt/=nary_tree%path_num_redunt) & 636 call error_stop("NaryTreePathGetReduntList", & 637 "incorrect number of redundant derivatives", & 638 path_num_redunt) 639 ! gets the list addresses of redundant derivatives of current path 640 num_cent = nary_tree%wt_node(nary_tree%order_geo) 641 call geom_total_redunt_list(nary_tree%num_atoms, num_cent, & 642 nary_tree%idx_cent(1:num_cent), & 643 nary_tree%order_cent(1:num_cent), & 644 .true., path_num_redunt, redunt_list) 645 end if 646#if defined(XTIME) 647 ! prints the CPU elapsed time 648 call xtimer_view(curr_time, "NaryTreePathGetReduntList", STDOUT) 649#endif 650 end subroutine NaryTreePathGetReduntList 651 652 !> \brief returns the expectation values with redundant derivatives of current path 653 !> \author Bin Gao 654 !> \date 2011-01-14 655 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 656 !> \param num_opt is the number of operators 657 !> \param path_num_unique is the number of unique derivatives of current path 658 !> \param num_dens is the number of AO density matrices 659 !> \param unique_expt contains the expectation values of unique geometric derivatives 660 !> \param num_redunt_geo is the total number of redundant derivatives, equals to 661 !> \var(3*num_atoms)^sum(\var(order_cent)) 662 !> \return redunt_expt contains the updated expectation values of redundant derivatives on exit 663 subroutine NaryTreePathSetReduntExpt(nary_tree, num_opt, path_num_unique, num_dens, & 664 unique_expt, num_redunt_geo, redunt_expt) 665 type(nary_tree_t), intent(in) :: nary_tree 666 integer, intent(in) :: num_opt 667 integer, intent(in) :: path_num_unique 668 integer, intent(in) :: num_dens 669 real(REALK), intent(in) :: unique_expt(num_opt,path_num_unique,num_dens) 670 integer, intent(in) :: num_redunt_geo 671 real(REALK), intent(inout) :: redunt_expt(num_opt,num_redunt_geo,num_dens) 672 integer num_cent !number of differentiated centers 673#if defined(XTIME) 674 real(REALK) curr_time !current CPU time 675 ! sets current CPU time 676 call xtimer_set(curr_time) 677#endif 678 ! checks if the N-ary tree is created 679 if (nary_tree%idx_path==0) then 680 call error_stop("NaryTreePathSetReduntExpt", "N-ary tree is not created", 0) 681 else if (nary_tree%order_geo==0) then 682 if (num_redunt_geo/=path_num_unique) & 683 call error_stop("NaryTreePathSetReduntExpt", & 684 "incorrect size of expectation values", & 685 path_num_unique) 686 redunt_expt = redunt_expt+unique_expt 687 else 688 ! checks the sizes of the expectation values 689 if (path_num_unique/=nary_tree%path_num_unique) & 690 call error_stop("NaryTreePathSetReduntExpt", & 691 "incorrect variable path_num_unique", & 692 path_num_unique) 693 if (num_redunt_geo/=(3*nary_tree%num_atoms)**nary_tree%order_geo) & 694 call error_stop("NaryTreePathSetReduntExpt", & 695 "incorrect variable num_redunt_geo", & 696 num_redunt_geo) 697 ! gets the expectation values with redundant geometric derivatives for current path 698 num_cent = nary_tree%wt_node(nary_tree%order_geo) 699 call geom_total_redunt_expectation(nary_tree%num_atoms, num_cent, & 700 nary_tree%idx_cent(1:num_cent), & 701 nary_tree%order_cent(1:num_cent), & 702 num_opt, path_num_unique, num_dens, & 703 unique_expt, num_redunt_geo, redunt_expt) 704 end if 705#if defined(XTIME) 706 ! prints the CPU elapsed time 707 call xtimer_view(curr_time, "NaryTreePathSetReduntExpt", STDOUT) 708#endif 709 end subroutine NaryTreePathSetReduntExpt 710 711!- !> \brief reorders geometric derivatives 712!- !> \author Bin Gao 713!- !> \date 2012-08-30 714!- !> \param nary_tree contains the information of N-ary tree for geometric derivatives 715!- !> \param dim_ints is the dimension of integrals 716!- !> \param path_num_unique is the number of unique derivatives of current path 717!- !> \param dim_opt is the dimension of operators 718!- !> \return contr_ints contains the contracted integrals 719!- subroutine NaryTreePathReorderInts(nary_tree, dim_ints, path_num_unique, dim_opt, contr_ints) 720!- type(nary_tree_t), intent(in) :: nary_tree 721!- integer, intent(in) :: dim_ints 722!- integer, intent(in) :: path_num_unique 723!- integer, intent(in) :: dim_opt 724!- real(REALK), intent(inout) :: contr_ints(dim_ints,path_num_unique,dim_opt) 725!- integer icent !incremental recorder over differentiated centers 726!- integer num_ints !number of integrals 727!- integer num_geo !number of geometric derivatives 728!- integer num_opt !number of operators 729!-#if defined(XTIME) 730!- real(REALK) curr_time !current CPU time 731!- ! sets current CPU time 732!- call xtimer_set(curr_time) 733!-#endif 734!- ! checks if the N-ary tree is created 735!- if (nary_tree%idx_path==0) then 736!- call error_stop("NaryTreePathReorderInts", "N-ary tree is not created", 0) 737!- else if (nary_tree%order_geo>1) then 738!- ! checks the sizes of the expectation values 739!- if (path_num_unique/=nary_tree%path_num_unique) & 740!- call error_stop("NaryTreePathReorderInts", & 741!- "incorrect variable path_num_unique", & 742!- path_num_unique) 743!- num_ints = dim_ints 744!- num_opt = path_num_unique*dim_opt 745!- do icent = 1, nary_tree%wt_node(nary_tree%order_geo) 746!- num_geo = (nary_tree%order_cent(icent)+1)*(nary_tree%order_cent(icent)+2)/2 747!- num_opt = num_opt/num_geo 748!- call geom_part_reorder(nary_tree%order_cent(icent), & 749!- num_ints, num_opt, contr_ints) 750!- num_ints = num_ints*num_geo 751!- end do 752!- end if 753!-#if defined(XTIME) 754!- ! prints the CPU elapsed time 755!- call xtimer_view(curr_time, "NaryTreePathReorderInts", STDOUT) 756!-#endif 757!- end subroutine NaryTreePathReorderInts 758 759 !> \brief visualizes the information of N-ary tree for geometric derivatives 760 !> \author Bin Gao 761 !> \date 2011-12-12 762 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 763 !> \param io_viewer is the logical unit number of the viewer 764 subroutine NaryTreeView(nary_tree, io_viewer) 765 type(nary_tree_t), intent(in) :: nary_tree 766 integer, intent(in) :: io_viewer 767 integer icent !incremental recorder over centers 768#if defined(XTIME) 769 real(REALK) curr_time !current CPU time 770 ! sets current CPU time 771 call xtimer_set(curr_time) 772#endif 773 ! checks if the N-ary tree is created 774 if (nary_tree%idx_path==0) then 775 write(io_viewer,100) "N-ary tree is not created!" 776 return 777 else 778 write(io_viewer,100) "number of atoms", nary_tree%num_atoms 779 if (allocated(nary_tree%idx_atoms)) then 780 do icent = 1, nary_tree%num_atoms 781 write(io_viewer,100) "indices of atoms to be used as differentiated centers", & 782 nary_tree%idx_atoms(icent) 783 end do 784 end if 785 write(io_viewer,100) "order of geometric derivatives", nary_tree%order_geo 786 write(io_viewer,100) "maximum number of differentiated centers", nary_tree%max_num_cent 787 write(io_viewer,100) "total number of different allowed paths in N-ary tree", & 788 nary_tree%num_paths 789 write(io_viewer,100) "number of unique geometric derivatives in the N-ary tree", & 790 nary_tree%num_unique_geo 791 write(io_viewer,100) "index of current path", nary_tree%idx_path 792 if (nary_tree%order_geo>0) then 793 write(io_viewer,110) (nary_tree%idx_node(icent),"(",nary_tree%wt_node(icent), & 794 ")",icent=1,nary_tree%order_geo) 795 write(io_viewer,120) (nary_tree%idx_cent(icent),"(",nary_tree%order_cent(icent), & 796 ")",icent=1,nary_tree%wt_node(nary_tree%order_geo)) 797 end if 798 write(io_viewer,100) "offset of unique derivatives of current path", & 799 nary_tree%path_offset 800 write(io_viewer,100) "number of unique derivatives of current path", & 801 nary_tree%path_num_unique 802 write(io_viewer,100) "number of redundant derivatives of current path", & 803 nary_tree%path_num_redunt 804 end if 805#if defined(XTIME) 806 ! prints the CPU elapsed time 807 call xtimer_view(curr_time, "NaryTreeView", STDOUT) 808#endif 809100 format("NaryTreeView>> ",A,I8) 810110 format("NaryTreeView>> indices and weights of atom nodes",12(I4,A,I3,A)) 811120 format("NaryTreeView>> generated differentiated centers",16(I3,A,I2,A)) 812 end subroutine NaryTreeView 813 814 !> \brief frees space taken by a N-ary tree for geometric derivatives 815 !> \author Bin Gao 816 !> \date 2011-12-12 817 !> \param nary_tree contains the information of N-ary tree for geometric derivatives 818 subroutine NaryTreeDestroy(nary_tree) 819 type(nary_tree_t), intent(inout) :: nary_tree 820#if defined(XTIME) 821 real(REALK) curr_time !current CPU time 822 ! sets current CPU time 823 call xtimer_set(curr_time) 824#endif 825 if (allocated(nary_tree%idx_atoms)) deallocate(nary_tree%idx_atoms) 826 if (allocated(nary_tree%idx_node)) deallocate(nary_tree%idx_node) 827 if (allocated(nary_tree%wt_node)) deallocate(nary_tree%wt_node) 828 if (allocated(nary_tree%idx_cent)) deallocate(nary_tree%idx_cent) 829 if (allocated(nary_tree%order_cent)) deallocate(nary_tree%order_cent) 830 nary_tree%num_atoms = 0 831 nary_tree%order_geo = 0 832 nary_tree%max_num_cent = 0 833 nary_tree%num_paths = 0 834 nary_tree%idx_path = 0 835#if defined(XTIME) 836 ! prints the CPU elapsed time 837 call xtimer_view(curr_time, "NaryTreeDestroy", STDOUT) 838#endif 839 end subroutine NaryTreeDestroy 840 841 !> \brief calculates the Cartesian multipole moment integrals using contracted 842 !> Gaussian type orbitals (GTOs) 843 !> \author Bin Gao 844 !> \date 2011-12-12 845 !> \param idx_bra is the atomic index of bra center 846 !> \param coord_bra contains the coordinates of bra center 847 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 848 !> \param num_prim_bra is the number of primitive Gaussians of bra center 849 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 850 !> \param num_contr_bra is the number of contractions of bra center 851 !> \param contr_coef_bra contains the contraction coefficients of bra center 852 !> \param idx_ket is the atomic index of ket center 853 !> \param coord_ket contains the coordinates of ket center 854 !> \param angular_ket is the angular number of ket center 855 !> \param num_prim_ket is the number of primitive Gaussians of ket center 856 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 857 !> \param num_contr_ket is the number of contractions of ket center 858 !> \param contr_coef_ket contains the contraction coefficients of ket center 859 !> \param spher_gto indicates if using spherical GTOs, otherwise Cartesian GTOs 860 !> \param info_LAO contains the information of London atomic orbital 861 !> \param order_elec is the order of electronic derivatives 862 !> \param idx_diporg is the atomic center of dipole origin (<1 for non-atomic center) 863 !> \param dipole_origin contains the coordinates of dipole origin 864 !> \param scal_const is the scale constant for Cartesian multipole moments 865 !> \param order_mom is the order of Cartesian multipole moments 866 !> \param order_mag_bra is the order of magnetic derivatives on bra center 867 !> \param order_mag_ket is the order of magnetic derivatives on ket center 868 !> \param order_mag_total is the order of total magnetic derivatives 869 !> \param order_ram_bra is the order of derivatives w.r.t. total rotational angular momentum on bra center 870 !> \param order_ram_ket is the order of derivatives w.r.t. total rotational angular momentum on ket center 871 !> \param order_ram_total is the order of total derivatives w.r.t. total rotational angular momentum 872 !> \param order_geo_bra is the order of geometric derivatives with respect to bra center 873 !> \param order_geo_ket is the order of geometric derivatives with respect to ket center 874 !> \param order_geo_mom is the order of geometric derivatives on dipole origin 875 !> \param nary_tree_total contains the information of N-ary tree for total geometric derivatives 876 !> \param num_gto_bra is the number of spherical/Cartesian GTOs on bra center 877 !> \param num_gto_ket is the number of spherical/Cartesian GTOs on ket center 878 !> \param num_opt is the number of operators including derivatives 879 !> \param mag_num_bra contains the magnetic numbers of spherical GTOs on bra center 880 !> \param mag_num_ket contains the magnetic numbers of spherical GTOs on ket center 881 !> \param powers_bra contains the Cartesian powers of Cartesian GTOs on bra center 882 !> \param powers_ket contains the Cartesian powers of Cartesian GTOs on ket center 883 !> \return contr_ints contains the calculated contracted integrals 884 subroutine IntGetCARMOM(idx_bra, coord_bra, angular_bra, num_prim_bra, & 885 exponent_bra, num_contr_bra, contr_coef_bra, & 886 idx_ket, coord_ket, angular_ket, num_prim_ket, & 887 exponent_ket, num_contr_ket, contr_coef_ket, & 888 spher_gto, info_LAO, order_elec, idx_diporg, & 889 dipole_origin, scal_const, order_mom, & 890 order_mag_bra, order_mag_ket, order_mag_total, & 891 order_ram_bra, order_ram_ket, order_ram_total, & 892 order_geo_bra, order_geo_ket, order_geo_mom, & 893 nary_tree_total, num_gto_bra, num_gto_ket, & 894 num_opt, contr_ints, mag_num_bra, mag_num_ket, & 895 powers_bra, powers_ket) 896 integer, intent(in) :: idx_bra 897 real(REALK), intent(in) :: coord_bra(3) 898 integer, intent(in) :: angular_bra 899 integer, intent(in) :: num_prim_bra 900 real(REALK), intent(in) :: exponent_bra(num_prim_bra) 901 integer, intent(in) :: num_contr_bra 902 real(REALK), intent(in) :: contr_coef_bra(num_contr_bra,num_prim_bra) 903 integer, intent(in) :: idx_ket 904 real(REALK), intent(in) :: coord_ket(3) 905 integer, intent(in) :: angular_ket 906 integer, intent(in) :: num_prim_ket 907 real(REALK), intent(in) :: exponent_ket(num_prim_ket) 908 integer, intent(in) :: num_contr_ket 909 real(REALK), intent(in) :: contr_coef_ket(num_contr_ket,num_prim_ket) 910 logical, optional, intent(in) :: spher_gto 911 type(london_ao_t), intent(in) :: info_LAO 912 integer, intent(in) :: order_elec 913 integer, intent(in) :: idx_diporg 914 real(REALK), intent(in) :: dipole_origin(3) 915 real(REALK), intent(in) :: scal_const 916 integer, intent(in) :: order_mom 917 integer, intent(in) :: order_mag_bra 918 integer, intent(in) :: order_mag_ket 919 integer, intent(in) :: order_mag_total 920 integer, intent(in) :: order_ram_bra 921 integer, intent(in) :: order_ram_ket 922 integer, intent(in) :: order_ram_total 923 integer, optional, intent(in) :: order_geo_bra 924 integer, optional, intent(in) :: order_geo_ket 925 integer, optional, intent(in) :: order_geo_mom 926 type(nary_tree_t), optional, intent(in) :: nary_tree_total 927 integer, intent(in) :: num_gto_bra 928 integer, intent(in) :: num_gto_ket 929 integer, intent(in) :: num_opt 930 real(REALK), intent(out) :: contr_ints(num_gto_bra,num_contr_bra, & 931 num_gto_ket,num_contr_ket,num_opt) 932 integer, optional, intent(in) :: mag_num_bra(num_gto_bra) 933 integer, optional, intent(in) :: mag_num_ket(num_gto_ket) 934 integer, optional, intent(in) :: powers_bra(3,num_gto_bra) 935 integer, optional, intent(in) :: powers_ket(3,num_gto_ket) 936#include "max_idx_non.h" 937 real(REALK) gauge_origin(3) !gauge origin of the magnetic vector potential 938 real(REALK) origin_London_PF(3) !origin of the London phase factor 939 integer iopt !incremental recorder over operators 940 logical p_spher_gto !arguments for Gen1Int (local) 941 integer p_order_geo_bra 942 integer p_order_geo_ket 943 integer p_order_geo_mom 944 integer p_num_cents 945 integer, allocatable :: p_idx_cent(:) 946 integer, allocatable :: p_order_cent(:) 947 integer icent 948 real(REALK), allocatable :: tmp_ints(:,:,:,:,:) !contracted integrals from Gen1Int 949 integer ierr !error information 950 ! sets the arguments for Gen1Int (local) 951 if (present(spher_gto)) then 952 p_spher_gto = spher_gto 953 else 954 p_spher_gto = .true. 955 end if 956 if (present(order_geo_bra)) then 957 p_order_geo_bra = order_geo_bra 958 else 959 p_order_geo_bra = 0 960 end if 961 if (present(order_geo_ket)) then 962 p_order_geo_ket = order_geo_ket 963 else 964 p_order_geo_ket = 0 965 end if 966 if (present(order_geo_mom)) then 967 p_order_geo_mom = order_geo_mom 968 else 969 p_order_geo_mom = 0 970 end if 971 if (present(nary_tree_total)) then 972 p_num_cents = nary_tree_total%wt_node(nary_tree_total%order_geo) 973 else 974 p_num_cents = 0 975 end if 976 if (p_num_cents>0) then 977 allocate(p_idx_cent(p_num_cents), stat=ierr) 978 if (ierr/=0) & 979 call error_stop("IntGetCARMOM", "failed to allocate p_idx_cent", p_num_cents) 980 allocate(p_order_cent(p_num_cents), stat=ierr) 981 if (ierr/=0) & 982 call error_stop("IntGetCARMOM", "failed to allocate p_order_cent", p_num_cents) 983 if (allocated(nary_tree_total%idx_atoms)) then 984 do icent = 1, p_num_cents 985 p_idx_cent(icent) = nary_tree_total%idx_atoms(nary_tree_total%idx_cent(icent)) 986 end do 987 else 988 p_idx_cent = nary_tree_total%idx_cent(1:p_num_cents) 989 end if 990 p_order_cent = nary_tree_total%order_cent(1:p_num_cents) 991 else 992 allocate(p_idx_cent(1), stat=ierr) 993 if (ierr/=0) call error_stop("IntGetCARMOM", "failed to allocate p_idx_cent", 1) 994 allocate(p_order_cent(1), stat=ierr) 995 if (ierr/=0) call error_stop("IntGetCARMOM", "failed to allocate p_order_cent", 1) 996 p_idx_cent = 0 997 p_order_cent = 0 998 end if 999 ! magnetic derivatives or derivatives with respect to rotational angular momentum 1000 if (order_mag_total>0 .or. order_mag_bra>0 .or. order_mag_ket>0 .or. & 1001 order_ram_total>0 .or. order_ram_bra>0 .or. order_ram_ket>0) then 1002 ! London atomic orbitals 1003 if (LondonAOUsed(info_LAO)) then 1004 if ((order_mag_total==0 .and. order_mom==0 .and. & 1005 ((order_mag_bra==1 .and. order_mag_ket==0) .or. & 1006 (order_mag_bra==0 .and. order_mag_ket==1)) .and. & 1007 order_elec==0 .and. p_order_geo_bra==0 .and. & 1008 p_order_geo_ket==0) .or. & 1009 (order_mag_total==1 .and. order_mom==0 .and. & 1010 order_mag_bra==0 .and. order_mag_ket==0 .and. & 1011 order_elec==0 .and. p_order_geo_bra==0 .and. & 1012 p_order_geo_ket==0)) then 1013 call LondonAOGetLPFOrigin(info_LAO, origin_London_PF) 1014 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1015 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1016 if (ierr/=0) & 1017 call error_stop("IntGetCARMOM", "failed to allocate tmp_ints", & 1018 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1019 ! SGTOs 1020 if (p_spher_gto) then 1021 if ((angular_bra==0 .and. angular_ket==0) .or. & 1022 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1023 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1024 exponent_bra, num_contr_bra, contr_coef_bra, & 1025 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1026 exponent_ket, num_contr_ket, contr_coef_ket, & 1027 order_elec, MAX_IDX_NON, origin_London_PF, & 1028 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1029 p_order_geo_ket, p_order_geo_mom, & 1030 p_num_cents, p_idx_cent, p_order_cent, & 1031 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1032 ! reorders integrals if required 1033 else 1034 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1035 exponent_bra, num_contr_bra, contr_coef_bra, & 1036 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1037 exponent_ket, num_contr_ket, contr_coef_ket, & 1038 order_elec, MAX_IDX_NON, origin_London_PF, & 1039 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1040 p_order_geo_ket, p_order_geo_mom, & 1041 p_num_cents, p_idx_cent, p_order_cent, & 1042 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1043 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1044 angular_ket>0 .and. present(mag_num_ket)) then 1045 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1046 angular_ket, num_gto_ket, mag_num_ket, & 1047 num_contr_bra, num_contr_ket, num_opt, & 1048 contr_ints, tmp_ints) 1049 else if (angular_bra>0 .and. present(mag_num_bra)) then 1050 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1051 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 1052 else if (angular_ket>0 .and. present(mag_num_ket)) then 1053 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1054 num_gto_bra*num_contr_bra, num_contr_ket, & 1055 num_opt, contr_ints, tmp_ints) 1056 end if 1057 end if 1058 ! CGTOs 1059 else 1060 if ((angular_bra==0 .and. angular_ket==0) .or. & 1061 .not.(present(powers_bra) .or. present(powers_ket))) then 1062 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1063 exponent_bra, num_contr_bra, contr_coef_bra, & 1064 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1065 exponent_ket, num_contr_ket, contr_coef_ket, & 1066 order_elec, MAX_IDX_NON, origin_London_PF, & 1067 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1068 p_order_geo_ket, p_order_geo_mom, & 1069 p_num_cents, p_idx_cent, p_order_cent, & 1070 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1071 ! reorders integrals if required 1072 else 1073 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1074 exponent_bra, num_contr_bra, contr_coef_bra, & 1075 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1076 exponent_ket, num_contr_ket, contr_coef_ket, & 1077 order_elec, MAX_IDX_NON, origin_London_PF, & 1078 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1079 p_order_geo_ket, p_order_geo_mom, & 1080 p_num_cents, p_idx_cent, p_order_cent, & 1081 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1082 if (angular_bra>0 .and. present(powers_bra) .and. & 1083 angular_ket>0 .and. present(powers_ket)) then 1084 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1085 angular_ket, num_gto_ket, powers_ket, & 1086 num_contr_bra, num_contr_ket, num_opt, & 1087 contr_ints, tmp_ints) 1088 else if (angular_bra>0 .and. present(powers_bra)) then 1089 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1090 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 1091 else if (angular_ket>0 .and. present(powers_ket)) then 1092 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1093 num_gto_bra*num_contr_bra, num_contr_ket, & 1094 num_opt, contr_ints, tmp_ints) 1095 end if 1096 end if 1097 end if 1098 ! sets the first order partial magnetic derivatives 1099 call LondonAOGetGaugeOrigin(info_LAO, gauge_origin) 1100 if (order_mag_bra==1 .and. order_mag_ket==0) then 1101 do iopt = 1, num_opt-2, 3 1102 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) & 1103 - (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) 1104 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) & 1105 - (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) 1106 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) & 1107 - (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) 1108 end do 1109 else if (order_mag_bra==0 .and. order_mag_ket==1) then 1110 do iopt = 1, num_opt-2, 3 1111 contr_ints(:,:,:,:,iopt) = (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) & 1112 - (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) 1113 contr_ints(:,:,:,:,iopt+1) = (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) & 1114 - (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) 1115 contr_ints(:,:,:,:,iopt+2) = (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) & 1116 - (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) 1117 end do 1118 else 1119 do iopt = 1, num_opt-2, 3 1120 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt+2) & 1121 - (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt+1) 1122 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt) & 1123 - (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+2) 1124 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+1) & 1125 - (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt) 1126 end do 1127 end if 1128 deallocate(tmp_ints) 1129 !FIXME: quick implementation of CM-1 integrals 1130 else if ((order_mag_total==1 .and. order_mom==1 .and. & 1131 order_mag_bra==0 .and. order_mag_ket==0 .and. & 1132 order_elec==0 .and. p_order_geo_bra==0 .and. & 1133 p_order_geo_ket==0)) then 1134 call LondonAOGetLPFOrigin(info_LAO, origin_London_PF) 1135 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1136 num_gto_ket,num_contr_ket,9), stat=ierr) 1137 if (ierr/=0) & 1138 call error_stop("IntGetCARMOM", "failed to allocate tmp_ints", & 1139 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*9) 1140 ! SGTOs 1141 if (p_spher_gto) then 1142 if ((angular_bra==0 .and. angular_ket==0) .or. & 1143 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1144 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1145 exponent_bra, num_contr_bra, contr_coef_bra, & 1146 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1147 exponent_ket, num_contr_ket, contr_coef_ket, & 1148 order_elec, MAX_IDX_NON, origin_London_PF, & 1149 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1150 p_order_geo_ket, p_order_geo_mom, & 1151 p_num_cents, p_idx_cent, p_order_cent, & 1152 num_gto_bra, num_gto_ket, 3, tmp_ints) 1153 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1154 exponent_bra, num_contr_bra, contr_coef_bra, & 1155 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1156 exponent_ket, num_contr_ket, contr_coef_ket, & 1157 order_elec, MAX_IDX_NON, origin_London_PF, & 1158 scal_const*0.5_REALK, 2, p_order_geo_bra, & 1159 p_order_geo_ket, p_order_geo_mom, & 1160 p_num_cents, p_idx_cent, p_order_cent, & 1161 num_gto_bra, num_gto_ket, 6, tmp_ints(:,:,:,:,4:9)) 1162 ! reorders integrals if required 1163 else 1164 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1165 exponent_bra, num_contr_bra, contr_coef_bra, & 1166 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1167 exponent_ket, num_contr_ket, contr_coef_ket, & 1168 order_elec, MAX_IDX_NON, origin_London_PF, & 1169 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1170 p_order_geo_ket, p_order_geo_mom, & 1171 p_num_cents, p_idx_cent, p_order_cent, & 1172 num_gto_bra, num_gto_ket, 3, contr_ints) 1173 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1174 exponent_bra, num_contr_bra, contr_coef_bra, & 1175 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1176 exponent_ket, num_contr_ket, contr_coef_ket, & 1177 order_elec, MAX_IDX_NON, origin_London_PF, & 1178 scal_const*0.5_REALK, 2, p_order_geo_bra, & 1179 p_order_geo_ket, p_order_geo_mom, & 1180 p_num_cents, p_idx_cent, p_order_cent, & 1181 num_gto_bra, num_gto_ket, 6, contr_ints(:,:,:,:,4:9)) 1182 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1183 angular_ket>0 .and. present(mag_num_ket)) then 1184 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1185 angular_ket, num_gto_ket, mag_num_ket, & 1186 num_contr_bra, num_contr_ket, 9, & 1187 contr_ints, tmp_ints) 1188 else if (angular_bra>0 .and. present(mag_num_bra)) then 1189 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1190 num_gto_ket*num_contr_ket*9, contr_ints, tmp_ints) 1191 else if (angular_ket>0 .and. present(mag_num_ket)) then 1192 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1193 num_gto_bra*num_contr_bra, num_contr_ket, & 1194 9, contr_ints, tmp_ints) 1195 end if 1196 end if 1197 ! CGTOs 1198 else 1199 if ((angular_bra==0 .and. angular_ket==0) .or. & 1200 .not.(present(powers_bra) .or. present(powers_ket))) then 1201 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1202 exponent_bra, num_contr_bra, contr_coef_bra, & 1203 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1204 exponent_ket, num_contr_ket, contr_coef_ket, & 1205 order_elec, MAX_IDX_NON, origin_London_PF, & 1206 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1207 p_order_geo_ket, p_order_geo_mom, & 1208 p_num_cents, p_idx_cent, p_order_cent, & 1209 num_gto_bra, num_gto_ket, 3, tmp_ints) 1210 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1211 exponent_bra, num_contr_bra, contr_coef_bra, & 1212 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1213 exponent_ket, num_contr_ket, contr_coef_ket, & 1214 order_elec, MAX_IDX_NON, origin_London_PF, & 1215 scal_const*0.5_REALK, 2, p_order_geo_bra, & 1216 p_order_geo_ket, p_order_geo_mom, & 1217 p_num_cents, p_idx_cent, p_order_cent, & 1218 num_gto_bra, num_gto_ket, 6, tmp_ints(:,:,:,:,4:9)) 1219 ! reorders integrals if required 1220 else 1221 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1222 exponent_bra, num_contr_bra, contr_coef_bra, & 1223 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1224 exponent_ket, num_contr_ket, contr_coef_ket, & 1225 order_elec, MAX_IDX_NON, origin_London_PF, & 1226 scal_const*0.5_REALK, 1, p_order_geo_bra, & 1227 p_order_geo_ket, p_order_geo_mom, & 1228 p_num_cents, p_idx_cent, p_order_cent, & 1229 num_gto_bra, num_gto_ket, 3, contr_ints) 1230 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1231 exponent_bra, num_contr_bra, contr_coef_bra, & 1232 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1233 exponent_ket, num_contr_ket, contr_coef_ket, & 1234 order_elec, MAX_IDX_NON, origin_London_PF, & 1235 scal_const*0.5_REALK, 2, p_order_geo_bra, & 1236 p_order_geo_ket, p_order_geo_mom, & 1237 p_num_cents, p_idx_cent, p_order_cent, & 1238 num_gto_bra, num_gto_ket, 6, contr_ints(:,:,:,:,4:9)) 1239 if (angular_bra>0 .and. present(powers_bra) .and. & 1240 angular_ket>0 .and. present(powers_ket)) then 1241 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1242 angular_ket, num_gto_ket, powers_ket, & 1243 num_contr_bra, num_contr_ket, 9, & 1244 contr_ints, tmp_ints) 1245 else if (angular_bra>0 .and. present(powers_bra)) then 1246 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1247 num_gto_ket*num_contr_ket*9, contr_ints, tmp_ints) 1248 else if (angular_ket>0 .and. present(powers_ket)) then 1249 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1250 num_gto_bra*num_contr_bra, num_contr_ket, & 1251 9, contr_ints, tmp_ints) 1252 end if 1253 end if 1254 end if 1255 ! sets the first order magnetic derivatives 1256 origin_London_PF(1) = origin_London_PF(1)-dipole_origin(1) 1257 origin_London_PF(2) = origin_London_PF(2)-dipole_origin(2) 1258 origin_London_PF(3) = origin_London_PF(3)-dipole_origin(3) 1259 ! x,Bx 1260 contr_ints(:,:,:,:,1) = (coord_bra(2)-coord_ket(2)) & 1261 * (tmp_ints(:,:,:,:,7)+origin_London_PF(1)*tmp_ints(:,:,:,:,3)) & 1262 - (coord_bra(3)-coord_ket(3)) & 1263 * (tmp_ints(:,:,:,:,5)+origin_London_PF(1)*tmp_ints(:,:,:,:,2)) 1264 ! y,Bx 1265 contr_ints(:,:,:,:,2) = (coord_bra(2)-coord_ket(2)) & 1266 * (tmp_ints(:,:,:,:,8)+origin_London_PF(2)*tmp_ints(:,:,:,:,3)) & 1267 - (coord_bra(3)-coord_ket(3)) & 1268 * (tmp_ints(:,:,:,:,6)+origin_London_PF(2)*tmp_ints(:,:,:,:,2)) 1269 ! z,Bx 1270 contr_ints(:,:,:,:,3) = (coord_bra(2)-coord_ket(2)) & 1271 * (tmp_ints(:,:,:,:,9)+origin_London_PF(3)*tmp_ints(:,:,:,:,3)) & 1272 - (coord_bra(3)-coord_ket(3)) & 1273 * (tmp_ints(:,:,:,:,8)+origin_London_PF(3)*tmp_ints(:,:,:,:,2)) 1274 ! x,By 1275 contr_ints(:,:,:,:,4) = (coord_bra(3)-coord_ket(3)) & 1276 * (tmp_ints(:,:,:,:,4)+origin_London_PF(1)*tmp_ints(:,:,:,:,1)) & 1277 - (coord_bra(1)-coord_ket(1)) & 1278 * (tmp_ints(:,:,:,:,7)+origin_London_PF(1)*tmp_ints(:,:,:,:,3)) 1279 ! y,By 1280 contr_ints(:,:,:,:,5) = (coord_bra(3)-coord_ket(3)) & 1281 * (tmp_ints(:,:,:,:,5)+origin_London_PF(2)*tmp_ints(:,:,:,:,1)) & 1282 - (coord_bra(1)-coord_ket(1)) & 1283 * (tmp_ints(:,:,:,:,8)+origin_London_PF(2)*tmp_ints(:,:,:,:,3)) 1284 ! z,By 1285 contr_ints(:,:,:,:,6) = (coord_bra(3)-coord_ket(3)) & 1286 * (tmp_ints(:,:,:,:,7)+origin_London_PF(3)*tmp_ints(:,:,:,:,1)) & 1287 - (coord_bra(1)-coord_ket(1)) & 1288 * (tmp_ints(:,:,:,:,9)+origin_London_PF(3)*tmp_ints(:,:,:,:,3)) 1289 ! x,Bz 1290 contr_ints(:,:,:,:,7) = (coord_bra(1)-coord_ket(1)) & 1291 * (tmp_ints(:,:,:,:,5)+origin_London_PF(1)*tmp_ints(:,:,:,:,2)) & 1292 - (coord_bra(2)-coord_ket(2)) & 1293 * (tmp_ints(:,:,:,:,4)+origin_London_PF(1)*tmp_ints(:,:,:,:,1)) 1294 ! y,Bz 1295 contr_ints(:,:,:,:,8) = (coord_bra(1)-coord_ket(1)) & 1296 * (tmp_ints(:,:,:,:,6)+origin_London_PF(2)*tmp_ints(:,:,:,:,2)) & 1297 - (coord_bra(2)-coord_ket(2)) & 1298 * (tmp_ints(:,:,:,:,5)+origin_London_PF(2)*tmp_ints(:,:,:,:,1)) 1299 ! z,Bz 1300 contr_ints(:,:,:,:,9) = (coord_bra(1)-coord_ket(1)) & 1301 * (tmp_ints(:,:,:,:,8)+origin_London_PF(3)*tmp_ints(:,:,:,:,2)) & 1302 - (coord_bra(2)-coord_ket(2)) & 1303 * (tmp_ints(:,:,:,:,7)+origin_London_PF(3)*tmp_ints(:,:,:,:,1)) 1304 deallocate(tmp_ints) 1305 else 1306 call error_stop("IntGetCARMOM", "LAO is not implemented", -1) 1307 end if 1308 ! non-LAOs 1309 else 1310 contr_ints = 0.0 1311 end if 1312 else 1313 ! SGTOs 1314 if (p_spher_gto) then 1315 if ((angular_bra==0 .and. angular_ket==0) .or. & 1316 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1317 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1318 exponent_bra, num_contr_bra, contr_coef_bra, & 1319 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1320 exponent_ket, num_contr_ket, contr_coef_ket, & 1321 order_elec, idx_diporg, dipole_origin, & 1322 scal_const, order_mom, p_order_geo_bra, & 1323 p_order_geo_ket, p_order_geo_mom, & 1324 p_num_cents, p_idx_cent, p_order_cent, & 1325 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1326 ! reorders integrals if required 1327 else 1328 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1329 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1330 if (ierr/=0) & 1331 call error_stop("IntGetCARMOM", "failed to allocate tmp_ints", & 1332 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1333 call contr_sgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1334 exponent_bra, num_contr_bra, contr_coef_bra, & 1335 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1336 exponent_ket, num_contr_ket, contr_coef_ket, & 1337 order_elec, idx_diporg, dipole_origin, & 1338 scal_const, order_mom, p_order_geo_bra, & 1339 p_order_geo_ket, p_order_geo_mom, & 1340 p_num_cents, p_idx_cent, p_order_cent, & 1341 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1342 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1343 angular_ket>0 .and. present(mag_num_ket)) then 1344 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1345 angular_ket, num_gto_ket, mag_num_ket, & 1346 num_contr_bra, num_contr_ket, num_opt, & 1347 tmp_ints, contr_ints) 1348 else if (angular_bra>0 .and. present(mag_num_bra)) then 1349 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1350 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 1351 else if (angular_ket>0 .and. present(mag_num_ket)) then 1352 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1353 num_gto_bra*num_contr_bra, num_contr_ket, & 1354 num_opt, tmp_ints, contr_ints) 1355 end if 1356 deallocate(tmp_ints) 1357 end if 1358 ! CGTOs 1359 else 1360 if ((angular_bra==0 .and. angular_ket==0) .or. & 1361 .not.(present(powers_bra) .or. present(powers_ket))) then 1362 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1363 exponent_bra, num_contr_bra, contr_coef_bra, & 1364 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1365 exponent_ket, num_contr_ket, contr_coef_ket, & 1366 order_elec, idx_diporg, dipole_origin, & 1367 scal_const, order_mom, p_order_geo_bra, & 1368 p_order_geo_ket, p_order_geo_mom, & 1369 p_num_cents, p_idx_cent, p_order_cent, & 1370 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1371 ! reorders integrals if required 1372 else 1373 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1374 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1375 if (ierr/=0) & 1376 call error_stop("IntGetCARMOM", "failed to allocate tmp_ints", & 1377 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1378 call contr_cgto_carmom(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1379 exponent_bra, num_contr_bra, contr_coef_bra, & 1380 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1381 exponent_ket, num_contr_ket, contr_coef_ket, & 1382 order_elec, idx_diporg, dipole_origin, & 1383 scal_const, order_mom, p_order_geo_bra, & 1384 p_order_geo_ket, p_order_geo_mom, & 1385 p_num_cents, p_idx_cent, p_order_cent, & 1386 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1387 if (angular_bra>0 .and. present(powers_bra) .and. & 1388 angular_ket>0 .and. present(powers_ket)) then 1389 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1390 angular_ket, num_gto_ket, powers_ket, & 1391 num_contr_bra, num_contr_ket, num_opt, & 1392 tmp_ints, contr_ints) 1393 else if (angular_bra>0 .and. present(powers_bra)) then 1394 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1395 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 1396 else if (angular_ket>0 .and. present(powers_ket)) then 1397 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1398 num_gto_bra*num_contr_bra, num_contr_ket, & 1399 num_opt, tmp_ints, contr_ints) 1400 end if 1401 deallocate(tmp_ints) 1402 end if 1403 end if 1404 end if 1405 deallocate(p_idx_cent) 1406 deallocate(p_order_cent) 1407 end subroutine IntGetCARMOM 1408 1409 !> \brief calculates the nuclear attraction potential integrals using contracted 1410 !> Gaussian type orbitals (GTOs) 1411 !> \author Bin Gao 1412 !> \date 2011-12-12 1413 !> \param idx_bra is the atomic index of bra center 1414 !> \param coord_bra contains the coordinates of bra center 1415 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 1416 !> \param num_prim_bra is the number of primitive Gaussians of bra center 1417 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 1418 !> \param num_contr_bra is the number of contractions of bra center 1419 !> \param contr_coef_bra contains the contraction coefficients of bra center 1420 !> \param idx_ket is the atomic index of ket center 1421 !> \param coord_ket contains the coordinates of ket center 1422 !> \param angular_ket is the angular number of ket center 1423 !> \param num_prim_ket is the number of primitive Gaussians of ket center 1424 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 1425 !> \param num_contr_ket is the number of contractions of ket center 1426 !> \param contr_coef_ket contains the contraction coefficients of ket center 1427 !> \param spher_gto indicates if using spherical GTOs, otherwise Cartesian GTOs 1428 !> \param info_LAO contains the information of London atomic orbital 1429 !> \param order_elec is the order of electronic derivatives 1430 !> \param idx_nucorg is the atomic center of nuclear potential origin (<1 for non-atomic center) 1431 !> \param nucpot_origin contains the coordinates of nuclear potential origin 1432 !> \param idx_diporg is the atomic center of dipole origin (<1 for non-atomic center) 1433 !> \param dipole_origin contains the coordinates of dipole origin 1434 !> \param scal_const is the scale constant for nuclear attraction potential operators 1435 !> \param order_mom is the order of Cartesian multipole moments 1436 !> \param order_mag_bra is the order of magnetic derivatives on bra center 1437 !> \param order_mag_ket is the order of magnetic derivatives on ket center 1438 !> \param order_mag_total is the order of total magnetic derivatives 1439 !> \param order_ram_bra is the order of derivatives w.r.t. total rotational angular momentum on bra center 1440 !> \param order_ram_ket is the order of derivatives w.r.t. total rotational angular momentum on ket center 1441 !> \param order_ram_total is the order of total derivatives w.r.t. total rotational angular momentum 1442 !> \param order_geo_bra is the order of geometric derivatives with respect to bra center 1443 !> \param order_geo_ket is the order of geometric derivatives with respect to ket center 1444 !> \param order_geo_pot is the order of geometric derivatives on nuclear attraction potential origin 1445 !> \param order_geo_mom is the order of geometric derivatives on dipole origin 1446 !> \param nary_tree_total contains the information of N-ary tree for total geometric derivatives 1447 !> \param num_gto_bra is the number of spherical/Cartesian GTOs on bra center 1448 !> \param num_gto_ket is the number of spherical/Cartesian GTOs on ket center 1449 !> \param num_opt is the number of operators including derivatives 1450 !> \param mag_num_bra contains the magnetic numbers of spherical GTOs on bra center 1451 !> \param mag_num_ket contains the magnetic numbers of spherical GTOs on ket center 1452 !> \param powers_bra contains the Cartesian powers of Cartesian GTOs on bra center 1453 !> \param powers_ket contains the Cartesian powers of Cartesian GTOs on ket center 1454 !> \return contr_ints contains the calculated contracted integrals 1455 subroutine IntGetNUCPOT(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1456 exponent_bra, num_contr_bra, contr_coef_bra, & 1457 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1458 exponent_ket, num_contr_ket, contr_coef_ket, & 1459 spher_gto, info_LAO, order_elec, idx_nucorg, & 1460 nucpot_origin, idx_diporg, dipole_origin, & 1461 scal_const, order_mom, & 1462 order_mag_bra, order_mag_ket, order_mag_total, & 1463 order_ram_bra, order_ram_ket, order_ram_total, & 1464 order_geo_bra, order_geo_ket, & 1465 order_geo_pot, order_geo_mom, & 1466 nary_tree_total, num_gto_bra, num_gto_ket, & 1467 num_opt, contr_ints, mag_num_bra, mag_num_ket, & 1468 powers_bra, powers_ket) 1469 integer, intent(in) :: idx_bra 1470 real(REALK), intent(in) :: coord_bra(3) 1471 integer, intent(in) :: angular_bra 1472 integer, intent(in) :: num_prim_bra 1473 real(REALK), intent(in) :: exponent_bra(num_prim_bra) 1474 integer, intent(in) :: num_contr_bra 1475 real(REALK), intent(in) :: contr_coef_bra(num_contr_bra,num_prim_bra) 1476 integer, intent(in) :: idx_ket 1477 real(REALK), intent(in) :: coord_ket(3) 1478 integer, intent(in) :: angular_ket 1479 integer, intent(in) :: num_prim_ket 1480 real(REALK), intent(in) :: exponent_ket(num_prim_ket) 1481 integer, intent(in) :: num_contr_ket 1482 real(REALK), intent(in) :: contr_coef_ket(num_contr_ket,num_prim_ket) 1483 logical, optional, intent(in) :: spher_gto 1484 type(london_ao_t), intent(in) :: info_LAO 1485 integer, intent(in) :: order_elec 1486 integer, intent(in) :: idx_nucorg 1487 real(REALK), intent(in) :: nucpot_origin(3) 1488 integer, intent(in) :: idx_diporg 1489 real(REALK), intent(in) :: dipole_origin(3) 1490 real(REALK), intent(in) :: scal_const 1491 integer, intent(in) :: order_mom 1492 integer, intent(in) :: order_mag_bra 1493 integer, intent(in) :: order_mag_ket 1494 integer, intent(in) :: order_mag_total 1495 integer, intent(in) :: order_ram_bra 1496 integer, intent(in) :: order_ram_ket 1497 integer, intent(in) :: order_ram_total 1498 integer, optional, intent(in) :: order_geo_bra 1499 integer, optional, intent(in) :: order_geo_ket 1500 integer, optional, intent(in) :: order_geo_pot 1501 integer, optional, intent(in) :: order_geo_mom 1502 type(nary_tree_t), optional, intent(in) :: nary_tree_total 1503 integer, intent(in) :: num_gto_bra 1504 integer, intent(in) :: num_gto_ket 1505 integer, intent(in) :: num_opt 1506 real(REALK), intent(out) :: contr_ints(num_gto_bra,num_contr_bra, & 1507 num_gto_ket,num_contr_ket,num_opt) 1508 integer, optional, intent(in) :: mag_num_bra(num_gto_bra) 1509 integer, optional, intent(in) :: mag_num_ket(num_gto_ket) 1510 integer, optional, intent(in) :: powers_bra(3,num_gto_bra) 1511 integer, optional, intent(in) :: powers_ket(3,num_gto_ket) 1512#include "max_idx_non.h" 1513 real(REALK) gauge_origin(3) !gauge origin of the magnetic vector potential 1514 real(REALK) origin_London_PF(3) !origin of the London phase factor 1515 integer iopt !incremental recorder over operators 1516 logical p_spher_gto !arguments for Gen1Int (local) 1517 integer p_order_geo_bra 1518 integer p_order_geo_ket 1519 integer p_order_geo_mom 1520 integer p_order_geo_pot 1521 integer p_num_cents 1522 integer, allocatable :: p_idx_cent(:) 1523 integer, allocatable :: p_order_cent(:) 1524 integer icent 1525 real(REALK), allocatable :: tmp_ints(:,:,:,:,:) !contracted integrals from Gen1Int 1526 real(REALK), allocatable :: ro_ints(:,:,:,:,:) 1527 integer ierr !error information 1528 ! sets the arguments for Gen1Int (local) 1529 if (present(spher_gto)) then 1530 p_spher_gto = spher_gto 1531 else 1532 p_spher_gto = .true. 1533 end if 1534 if (present(order_geo_bra)) then 1535 p_order_geo_bra = order_geo_bra 1536 else 1537 p_order_geo_bra = 0 1538 end if 1539 if (present(order_geo_ket)) then 1540 p_order_geo_ket = order_geo_ket 1541 else 1542 p_order_geo_ket = 0 1543 end if 1544 if (present(order_geo_pot)) then 1545 p_order_geo_pot = order_geo_pot 1546 else 1547 p_order_geo_pot = 0 1548 end if 1549 if (present(order_geo_mom)) then 1550 p_order_geo_mom = order_geo_mom 1551 else 1552 p_order_geo_mom = 0 1553 end if 1554 if (present(nary_tree_total)) then 1555 p_num_cents = nary_tree_total%wt_node(nary_tree_total%order_geo) 1556 else 1557 p_num_cents = 0 1558 end if 1559 if (p_num_cents>0) then 1560 allocate(p_idx_cent(p_num_cents), stat=ierr) 1561 if (ierr/=0) & 1562 call error_stop("IntGetNUCPOT", "failed to allocate p_idx_cent", p_num_cents) 1563 allocate(p_order_cent(p_num_cents), stat=ierr) 1564 if (ierr/=0) & 1565 call error_stop("IntGetNUCPOT", "failed to allocate p_order_cent", p_num_cents) 1566 if (allocated(nary_tree_total%idx_atoms)) then 1567 do icent = 1, p_num_cents 1568 p_idx_cent(icent) = nary_tree_total%idx_atoms(nary_tree_total%idx_cent(icent)) 1569 end do 1570 else 1571 p_idx_cent = nary_tree_total%idx_cent(1:p_num_cents) 1572 end if 1573 p_order_cent = nary_tree_total%order_cent(1:p_num_cents) 1574 else 1575 allocate(p_idx_cent(1), stat=ierr) 1576 if (ierr/=0) call error_stop("IntGetNUCPOT", "failed to allocate p_idx_cent", 1) 1577 allocate(p_order_cent(1), stat=ierr) 1578 if (ierr/=0) call error_stop("IntGetNUCPOT", "failed to allocate p_order_cent", 1) 1579 p_idx_cent = 0 1580 p_order_cent = 0 1581 end if 1582 ! magnetic derivatives or derivatives with respect to rotational angular momentum 1583 if (order_mag_total>0 .or. order_mag_bra>0 .or. order_mag_ket>0 .or. & 1584 order_ram_total>0 .or. order_ram_bra>0 .or. order_ram_ket>0) then 1585 ! London atomic orbitals 1586 if (LondonAOUsed(info_LAO)) then 1587 if ((order_mag_total==0 .and. order_mom==0 .and. & 1588 ((order_mag_bra==1 .and. order_mag_ket==0) .or. & 1589 (order_mag_bra==0 .and. order_mag_ket==1)) .and. & 1590 order_elec==0 .and. p_order_geo_bra==0 .and. & 1591 p_order_geo_ket==0 .and. p_num_cents<2) .or. & 1592 (order_mag_total==1 .and. order_mom==0 .and. & 1593 order_mag_bra==0 .and. order_mag_ket==0 .and. & 1594 order_elec==0 .and. p_order_geo_bra==0 .and. & 1595 p_order_geo_ket==0 .and. p_num_cents<2)) then 1596 if ((order_mag_total==1 .and. idx_bra/=idx_ket) .or. order_mag_total==0) then 1597 call LondonAOGetLPFOrigin(info_LAO, origin_London_PF) 1598 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1599 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1600 if (ierr/=0) & 1601 call error_stop("IntGetNUCPOT", "failed to allocate tmp_ints", & 1602 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1603 ! SGTOs 1604 if (p_spher_gto) then 1605 if ((angular_bra==0 .and. angular_ket==0) .or. & 1606 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1607 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1608 exponent_bra, num_contr_bra, contr_coef_bra, & 1609 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1610 exponent_ket, num_contr_ket, contr_coef_ket, & 1611 order_elec, idx_nucorg, nucpot_origin, & 1612 MAX_IDX_NON, origin_London_PF, & 1613 scal_const*0.5_REALK, 1, & 1614 p_order_geo_bra, p_order_geo_ket, & 1615 p_order_geo_pot, p_order_geo_mom, & 1616 p_num_cents, p_idx_cent, p_order_cent, & 1617 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1618 ! reorders integrals if required 1619 else 1620 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1621 exponent_bra, num_contr_bra, contr_coef_bra, & 1622 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1623 exponent_ket, num_contr_ket, contr_coef_ket, & 1624 order_elec, idx_nucorg, nucpot_origin, & 1625 MAX_IDX_NON, origin_London_PF, & 1626 scal_const*0.5_REALK, 1, & 1627 p_order_geo_bra, p_order_geo_ket, & 1628 p_order_geo_pot, p_order_geo_mom, & 1629 p_num_cents, p_idx_cent, p_order_cent, & 1630 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1631 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1632 angular_ket>0 .and. present(mag_num_ket)) then 1633 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1634 angular_ket, num_gto_ket, mag_num_ket, & 1635 num_contr_bra, num_contr_ket, num_opt, & 1636 contr_ints, tmp_ints) 1637 else if (angular_bra>0 .and. present(mag_num_bra)) then 1638 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1639 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 1640 else if (angular_ket>0 .and. present(mag_num_ket)) then 1641 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1642 num_gto_bra*num_contr_bra, num_contr_ket, & 1643 num_opt, contr_ints, tmp_ints) 1644 end if 1645 end if 1646 ! CGTOs 1647 else 1648 if ((angular_bra==0 .and. angular_ket==0) .or. & 1649 .not.(present(powers_bra) .or. present(powers_ket))) then 1650 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1651 exponent_bra, num_contr_bra, contr_coef_bra, & 1652 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1653 exponent_ket, num_contr_ket, contr_coef_ket, & 1654 order_elec, idx_nucorg, nucpot_origin, & 1655 MAX_IDX_NON, origin_London_PF, & 1656 scal_const*0.5_REALK, 1, & 1657 p_order_geo_bra, p_order_geo_ket, & 1658 p_order_geo_pot, p_order_geo_mom, & 1659 p_num_cents, p_idx_cent, p_order_cent, & 1660 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1661 ! reorders integrals if required 1662 else 1663 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1664 exponent_bra, num_contr_bra, contr_coef_bra, & 1665 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1666 exponent_ket, num_contr_ket, contr_coef_ket, & 1667 order_elec, idx_nucorg, nucpot_origin, & 1668 MAX_IDX_NON, origin_London_PF, & 1669 scal_const*0.5_REALK, 1, & 1670 p_order_geo_bra, p_order_geo_ket, & 1671 p_order_geo_pot, p_order_geo_mom, & 1672 p_num_cents, p_idx_cent, p_order_cent, & 1673 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1674 if (angular_bra>0 .and. present(powers_bra) .and. & 1675 angular_ket>0 .and. present(powers_ket)) then 1676 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1677 angular_ket, num_gto_ket, powers_ket, & 1678 num_contr_bra, num_contr_ket, num_opt, & 1679 contr_ints, tmp_ints) 1680 else if (angular_bra>0 .and. present(powers_bra)) then 1681 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1682 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 1683 else if (angular_ket>0 .and. present(powers_ket)) then 1684 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1685 num_gto_bra*num_contr_bra, num_contr_ket, & 1686 num_opt, contr_ints, tmp_ints) 1687 end if 1688 end if 1689 end if 1690 ! sets the first order partial magnetic derivatives 1691 call LondonAOGetGaugeOrigin(info_LAO, gauge_origin) 1692 if (order_mag_bra==1) then 1693 do iopt = 1, num_opt-2, 3 1694 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) & 1695 - (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) 1696 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) & 1697 - (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) 1698 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) & 1699 - (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) 1700 end do 1701 else if (order_mag_ket==1) then 1702 do iopt = 1, num_opt-2, 3 1703 contr_ints(:,:,:,:,iopt) = (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) & 1704 - (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) 1705 contr_ints(:,:,:,:,iopt+1) = (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) & 1706 - (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) 1707 contr_ints(:,:,:,:,iopt+2) = (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) & 1708 - (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) 1709 end do 1710 else 1711 do iopt = 1, num_opt-2, 3 1712 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt+2) & 1713 - (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt+1) 1714 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt) & 1715 - (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+2) 1716 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+1) & 1717 - (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt) 1718 end do 1719 end if 1720 deallocate(tmp_ints) 1721 ! mixed total geometric derivatives and magnetic derivatives 1722 if (p_num_cents==1) then 1723 if (p_order_cent(1)==1) then 1724 if ((p_idx_cent(1)==idx_bra .and. (order_mag_bra==1 .or. order_mag_total==1)) .or. & 1725 (p_idx_cent(1)==idx_ket .and. (order_mag_ket==1 .or. order_mag_total==1))) then 1726 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1727 num_gto_ket,num_contr_ket,num_opt/3), stat=ierr) 1728 if (ierr/=0) & 1729 call error_stop("IntGetNUCPOT", "failed to allocate tmp_ints", & 1730 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt/3) 1731 allocate(ro_ints(num_gto_bra,num_contr_bra, & 1732 num_gto_ket,num_contr_ket,num_opt/3), stat=ierr) 1733 if (ierr/=0) & 1734 call error_stop("IntGetNUCPOT", "failed to allocate ro_ints", & 1735 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt/3) 1736 ! SGTOs 1737 if (p_spher_gto) then 1738 if ((angular_bra==0 .and. angular_ket==0) .or. & 1739 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1740 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1741 exponent_bra, num_contr_bra, contr_coef_bra, & 1742 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1743 exponent_ket, num_contr_ket, contr_coef_ket, & 1744 order_elec, idx_nucorg, nucpot_origin, & 1745 MAX_IDX_NON, origin_London_PF, & 1746 scal_const*0.5_REALK, 1, & 1747 p_order_geo_bra, p_order_geo_ket, & 1748 p_order_geo_pot, p_order_geo_mom, & 1749 0, (/0/), (/0/), & 1750 num_gto_bra, num_gto_ket, num_opt/3, ro_ints) 1751 ! reorders integrals if required 1752 else 1753 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1754 exponent_bra, num_contr_bra, contr_coef_bra, & 1755 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1756 exponent_ket, num_contr_ket, contr_coef_ket, & 1757 order_elec, idx_nucorg, nucpot_origin, & 1758 MAX_IDX_NON, origin_London_PF, & 1759 scal_const*0.5_REALK, 1, & 1760 p_order_geo_bra, p_order_geo_ket, & 1761 p_order_geo_pot, p_order_geo_mom, & 1762 0, (/0/), (/0/), & 1763 num_gto_bra, num_gto_ket, num_opt/3, tmp_ints) 1764 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1765 angular_ket>0 .and. present(mag_num_ket)) then 1766 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1767 angular_ket, num_gto_ket, mag_num_ket, & 1768 num_contr_bra, num_contr_ket, num_opt/3, & 1769 tmp_ints, ro_ints) 1770 else if (angular_bra>0 .and. present(mag_num_bra)) then 1771 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1772 num_gto_ket*num_contr_ket*num_opt/3, tmp_ints, ro_ints) 1773 else if (angular_ket>0 .and. present(mag_num_ket)) then 1774 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1775 num_gto_bra*num_contr_bra, num_contr_ket, & 1776 num_opt/3, tmp_ints, ro_ints) 1777 end if 1778 end if 1779 ! CGTOs 1780 else 1781 if ((angular_bra==0 .and. angular_ket==0) .or. & 1782 .not.(present(powers_bra) .or. present(powers_ket))) then 1783 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1784 exponent_bra, num_contr_bra, contr_coef_bra, & 1785 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1786 exponent_ket, num_contr_ket, contr_coef_ket, & 1787 order_elec, idx_nucorg, nucpot_origin, & 1788 MAX_IDX_NON, origin_London_PF, & 1789 scal_const*0.5_REALK, 1, & 1790 p_order_geo_bra, p_order_geo_ket, & 1791 p_order_geo_pot, p_order_geo_mom, & 1792 0, (/0/), (/0/), & 1793 num_gto_bra, num_gto_ket, num_opt/3, ro_ints) 1794 ! reorders integrals if required 1795 else 1796 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1797 exponent_bra, num_contr_bra, contr_coef_bra, & 1798 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1799 exponent_ket, num_contr_ket, contr_coef_ket, & 1800 order_elec, idx_nucorg, nucpot_origin, & 1801 MAX_IDX_NON, origin_London_PF, & 1802 scal_const*0.5_REALK, 1, & 1803 p_order_geo_bra, p_order_geo_ket, & 1804 p_order_geo_pot, p_order_geo_mom, & 1805 0, (/0/), (/0/), & 1806 num_gto_bra, num_gto_ket, num_opt/3, tmp_ints) 1807 if (angular_bra>0 .and. present(powers_bra) .and. & 1808 angular_ket>0 .and. present(powers_ket)) then 1809 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1810 angular_ket, num_gto_ket, powers_ket, & 1811 num_contr_bra, num_contr_ket, num_opt/3, & 1812 tmp_ints, ro_ints) 1813 else if (angular_bra>0 .and. present(powers_bra)) then 1814 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1815 num_gto_ket*num_contr_ket*num_opt/3, tmp_ints, ro_ints) 1816 else if (angular_ket>0 .and. present(powers_ket)) then 1817 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1818 num_gto_bra*num_contr_bra, num_contr_ket, & 1819 num_opt/3, tmp_ints, ro_ints) 1820 end if 1821 end if 1822 end if 1823 deallocate(tmp_ints) 1824 ! first order partial magnetic derivatives on the bra center or total magnetic derivatives 1825 if (p_idx_cent(1)==idx_bra) then 1826 ! Bx, Gx 1827 ! By, Gx 1828 contr_ints(:,:,:,:,2) = contr_ints(:,:,:,:,2)-ro_ints(:,:,:,:,3) 1829 ! Bz, Gx 1830 contr_ints(:,:,:,:,3) = contr_ints(:,:,:,:,3)+ro_ints(:,:,:,:,2) 1831 ! Bx, Gy 1832 contr_ints(:,:,:,:,4) = contr_ints(:,:,:,:,4)+ro_ints(:,:,:,:,3) 1833 ! By, Gy 1834 ! Bz, Gy 1835 contr_ints(:,:,:,:,6) = contr_ints(:,:,:,:,6)-ro_ints(:,:,:,:,1) 1836 ! Bx, Gz 1837 contr_ints(:,:,:,:,7) = contr_ints(:,:,:,:,7)-ro_ints(:,:,:,:,2) 1838 ! By, Gz 1839 contr_ints(:,:,:,:,8) = contr_ints(:,:,:,:,8)+ro_ints(:,:,:,:,1) 1840 ! Bz, Gz 1841 ! first order partial magnetic derivatives on the ket center or total magnetic derivatives 1842 else 1843 ! Bx, Gx 1844 ! By, Gx 1845 contr_ints(:,:,:,:,2) = contr_ints(:,:,:,:,2)+ro_ints(:,:,:,:,3) 1846 ! Bz, Gx 1847 contr_ints(:,:,:,:,3) = contr_ints(:,:,:,:,3)-ro_ints(:,:,:,:,2) 1848 ! Bx, Gy 1849 contr_ints(:,:,:,:,4) = contr_ints(:,:,:,:,4)-ro_ints(:,:,:,:,3) 1850 ! By, Gy 1851 ! Bz, Gy 1852 contr_ints(:,:,:,:,6) = contr_ints(:,:,:,:,6)+ro_ints(:,:,:,:,1) 1853 ! Bx, Gz 1854 contr_ints(:,:,:,:,7) = contr_ints(:,:,:,:,7)+ro_ints(:,:,:,:,2) 1855 ! By, Gz 1856 contr_ints(:,:,:,:,8) = contr_ints(:,:,:,:,8)-ro_ints(:,:,:,:,1) 1857 ! Bz, Gz 1858 end if 1859 deallocate(ro_ints) 1860 end if 1861 else 1862 call error_stop("IntGetNUCPOT", "LAO is not implemented", -1) 1863 end if 1864 end if 1865 else 1866 contr_ints = 0.0 1867 end if 1868 else 1869 call error_stop("IntGetNUCPOT", "LAO is not implemented", -1) 1870 end if 1871 ! non-LAOs 1872 else 1873 contr_ints = 0.0 1874 end if 1875 else 1876 ! SGTOs 1877 if (p_spher_gto) then 1878 if ((angular_bra==0 .and. angular_ket==0) .or. & 1879 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 1880 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1881 exponent_bra, num_contr_bra, contr_coef_bra, & 1882 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1883 exponent_ket, num_contr_ket, contr_coef_ket, & 1884 order_elec, idx_nucorg, nucpot_origin, & 1885 idx_diporg, dipole_origin, & 1886 scal_const, order_mom, & 1887 p_order_geo_bra, p_order_geo_ket, & 1888 p_order_geo_pot, p_order_geo_mom, & 1889 p_num_cents, p_idx_cent, p_order_cent, & 1890 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1891 ! reorders integrals if required 1892 else 1893 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1894 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1895 if (ierr/=0) & 1896 call error_stop("IntGetNUCPOT", "failed to allocate tmp_ints", & 1897 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1898 call contr_sgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1899 exponent_bra, num_contr_bra, contr_coef_bra, & 1900 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1901 exponent_ket, num_contr_ket, contr_coef_ket, & 1902 order_elec, idx_nucorg, nucpot_origin, & 1903 idx_diporg, dipole_origin, & 1904 scal_const, order_mom, & 1905 p_order_geo_bra, p_order_geo_ket, & 1906 p_order_geo_pot, p_order_geo_mom, & 1907 p_num_cents, p_idx_cent, p_order_cent, & 1908 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1909 if (angular_bra>0 .and. present(mag_num_bra) .and. & 1910 angular_ket>0 .and. present(mag_num_ket)) then 1911 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 1912 angular_ket, num_gto_ket, mag_num_ket, & 1913 num_contr_bra, num_contr_ket, num_opt, & 1914 tmp_ints, contr_ints) 1915 else if (angular_bra>0 .and. present(mag_num_bra)) then 1916 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 1917 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 1918 else if (angular_ket>0 .and. present(mag_num_ket)) then 1919 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 1920 num_gto_bra*num_contr_bra, num_contr_ket, & 1921 num_opt, tmp_ints, contr_ints) 1922 end if 1923 deallocate(tmp_ints) 1924 end if 1925 ! CGTOs 1926 else 1927 if ((angular_bra==0 .and. angular_ket==0) .or. & 1928 .not.(present(powers_bra) .or. present(powers_ket))) then 1929 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1930 exponent_bra, num_contr_bra, contr_coef_bra, & 1931 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1932 exponent_ket, num_contr_ket, contr_coef_ket, & 1933 order_elec, idx_nucorg, nucpot_origin, & 1934 idx_diporg, dipole_origin, & 1935 scal_const, order_mom, & 1936 p_order_geo_bra, p_order_geo_ket, & 1937 p_order_geo_pot, p_order_geo_mom, & 1938 p_num_cents, p_idx_cent, p_order_cent, & 1939 num_gto_bra, num_gto_ket, num_opt, contr_ints) 1940 ! reorders integrals if required 1941 else 1942 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 1943 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 1944 if (ierr/=0) & 1945 call error_stop("IntGetNUCPOT", "failed to allocate tmp_ints", & 1946 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 1947 call contr_cgto_nucpot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 1948 exponent_bra, num_contr_bra, contr_coef_bra, & 1949 idx_ket, coord_ket, angular_ket, num_prim_ket, & 1950 exponent_ket, num_contr_ket, contr_coef_ket, & 1951 order_elec, idx_nucorg, nucpot_origin, & 1952 idx_diporg, dipole_origin, & 1953 scal_const, order_mom, & 1954 p_order_geo_bra, p_order_geo_ket, & 1955 p_order_geo_pot, p_order_geo_mom, & 1956 p_num_cents, p_idx_cent, p_order_cent, & 1957 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 1958 if (angular_bra>0 .and. present(powers_bra) .and. & 1959 angular_ket>0 .and. present(powers_ket)) then 1960 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 1961 angular_ket, num_gto_ket, powers_ket, & 1962 num_contr_bra, num_contr_ket, num_opt, & 1963 tmp_ints, contr_ints) 1964 else if (angular_bra>0 .and. present(powers_bra)) then 1965 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 1966 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 1967 else if (angular_ket>0 .and. present(powers_ket)) then 1968 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 1969 num_gto_bra*num_contr_bra, num_contr_ket, & 1970 num_opt, tmp_ints, contr_ints) 1971 end if 1972 deallocate(tmp_ints) 1973 end if 1974 end if 1975 end if 1976 deallocate(p_idx_cent) 1977 deallocate(p_order_cent) 1978 end subroutine IntGetNUCPOT 1979 1980 !> \brief calculates the Gaussian charge potential integrals using contracted 1981 !> Gaussian type orbitals (GTOs) 1982 !> \author Bin Gao 1983 !> \date 2011-12-12 1984 !> \param idx_bra is the atomic index of bra center 1985 !> \param coord_bra contains the coordinates of bra center 1986 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 1987 !> \param num_prim_bra is the number of primitive Gaussians of bra center 1988 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 1989 !> \param num_contr_bra is the number of contractions of bra center 1990 !> \param contr_coef_bra contains the contraction coefficients of bra center 1991 !> \param idx_ket is the atomic index of ket center 1992 !> \param coord_ket contains the coordinates of ket center 1993 !> \param angular_ket is the angular number of ket center 1994 !> \param num_prim_ket is the number of primitive Gaussians of ket center 1995 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 1996 !> \param num_contr_ket is the number of contractions of ket center 1997 !> \param contr_coef_ket contains the contraction coefficients of ket center 1998 !> \param spher_gto indicates if using spherical GTOs, otherwise Cartesian GTOs 1999 !> \param info_LAO contains the information of London atomic orbital 2000 !> \param order_elec is the order of electronic derivatives 2001 !> \param idx_gauorg is the atomic center of Gaussian charge potential origin (<1 for non-atomic center) 2002 !> \param gaupot_origin contains the coordinates of Gaussian charge potential origin 2003 !> \param gaupot_expt is the exponent used in the Gaussian broadening function of the charge 2004 !> \param idx_diporg is the atomic center of dipole origin (<1 for non-atomic center) 2005 !> \param dipole_origin contains the coordinates of dipole origin 2006 !> \param scal_const is the scale constant for potential operators 2007 !> \param order_mom is the order of Cartesian multipole moments 2008 !> \param order_mag_bra is the order of magnetic derivatives on bra center 2009 !> \param order_mag_ket is the order of magnetic derivatives on ket center 2010 !> \param order_mag_total is the order of total magnetic derivatives 2011 !> \param order_ram_bra is the order of derivatives w.r.t. total rotational angular momentum on bra center 2012 !> \param order_ram_ket is the order of derivatives w.r.t. total rotational angular momentum on ket center 2013 !> \param order_ram_total is the order of total derivatives w.r.t. total rotational angular momentum 2014 !> \param order_geo_bra is the order of geometric derivatives with respect to bra center 2015 !> \param order_geo_ket is the order of geometric derivatives with respect to ket center 2016 !> \param order_geo_pot is the order of geometric derivatives on potential origin 2017 !> \param order_geo_mom is the order of geometric derivatives on dipole origin 2018 !> \param nary_tree_total contains the information of N-ary tree for total geometric derivatives 2019 !> \param num_gto_bra is the number of spherical/Cartesian GTOs on bra center 2020 !> \param num_gto_ket is the number of spherical/Cartesian GTOs on ket center 2021 !> \param num_opt is the number of operators including derivatives 2022 !> \param mag_num_bra contains the magnetic numbers of spherical GTOs on bra center 2023 !> \param mag_num_ket contains the magnetic numbers of spherical GTOs on ket center 2024 !> \param powers_bra contains the Cartesian powers of Cartesian GTOs on bra center 2025 !> \param powers_ket contains the Cartesian powers of Cartesian GTOs on ket center 2026 !> \return contr_ints contains the calculated contracted integrals 2027 subroutine IntGetGAUPOT(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2028 exponent_bra, num_contr_bra, contr_coef_bra, & 2029 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2030 exponent_ket, num_contr_ket, contr_coef_ket, & 2031 spher_gto, info_LAO, order_elec, idx_gauorg, & 2032 gaupot_origin, gaupot_expt, idx_diporg, & 2033 dipole_origin, scal_const, order_mom, & 2034 order_mag_bra, order_mag_ket, order_mag_total, & 2035 order_ram_bra, order_ram_ket, order_ram_total, & 2036 order_geo_bra, order_geo_ket, & 2037 order_geo_pot, order_geo_mom, & 2038 nary_tree_total, num_gto_bra, num_gto_ket, & 2039 num_opt, contr_ints, mag_num_bra, mag_num_ket, & 2040 powers_bra, powers_ket) 2041 integer, intent(in) :: idx_bra 2042 real(REALK), intent(in) :: coord_bra(3) 2043 integer, intent(in) :: angular_bra 2044 integer, intent(in) :: num_prim_bra 2045 real(REALK), intent(in) :: exponent_bra(num_prim_bra) 2046 integer, intent(in) :: num_contr_bra 2047 real(REALK), intent(in) :: contr_coef_bra(num_contr_bra,num_prim_bra) 2048 integer, intent(in) :: idx_ket 2049 real(REALK), intent(in) :: coord_ket(3) 2050 integer, intent(in) :: angular_ket 2051 integer, intent(in) :: num_prim_ket 2052 real(REALK), intent(in) :: exponent_ket(num_prim_ket) 2053 integer, intent(in) :: num_contr_ket 2054 real(REALK), intent(in) :: contr_coef_ket(num_contr_ket,num_prim_ket) 2055 logical, optional, intent(in) :: spher_gto 2056 type(london_ao_t), intent(in) :: info_LAO 2057 integer, intent(in) :: order_elec 2058 integer, intent(in) :: idx_gauorg 2059 real(REALK), intent(in) :: gaupot_origin(3) 2060 real(REALK), intent(in) :: gaupot_expt 2061 integer, intent(in) :: idx_diporg 2062 real(REALK), intent(in) :: dipole_origin(3) 2063 real(REALK), intent(in) :: scal_const 2064 integer, intent(in) :: order_mom 2065 integer, intent(in) :: order_mag_bra 2066 integer, intent(in) :: order_mag_ket 2067 integer, intent(in) :: order_mag_total 2068 integer, intent(in) :: order_ram_bra 2069 integer, intent(in) :: order_ram_ket 2070 integer, intent(in) :: order_ram_total 2071 integer, optional, intent(in) :: order_geo_bra 2072 integer, optional, intent(in) :: order_geo_ket 2073 integer, optional, intent(in) :: order_geo_pot 2074 integer, optional, intent(in) :: order_geo_mom 2075 type(nary_tree_t), optional, intent(in) :: nary_tree_total 2076 integer, intent(in) :: num_gto_bra 2077 integer, intent(in) :: num_gto_ket 2078 integer, intent(in) :: num_opt 2079 real(REALK), intent(out) :: contr_ints(num_gto_bra,num_contr_bra, & 2080 num_gto_ket,num_contr_ket,num_opt) 2081 integer, optional, intent(in) :: mag_num_bra(num_gto_bra) 2082 integer, optional, intent(in) :: mag_num_ket(num_gto_ket) 2083 integer, optional, intent(in) :: powers_bra(3,num_gto_bra) 2084 integer, optional, intent(in) :: powers_ket(3,num_gto_ket) 2085#include "max_idx_non.h" 2086 real(REALK) gauge_origin(3) !gauge origin of the magnetic vector potential 2087 real(REALK) origin_London_PF(3) !origin of the London phase factor 2088 integer iopt !incremental recorder over operators 2089 logical p_spher_gto !arguments for Gen1Int (local) 2090 integer p_order_geo_bra 2091 integer p_order_geo_ket 2092 integer p_order_geo_mom 2093 integer p_order_geo_pot 2094 integer p_num_cents 2095 integer, allocatable :: p_idx_cent(:) 2096 integer, allocatable :: p_order_cent(:) 2097 integer icent 2098 real(REALK), allocatable :: tmp_ints(:,:,:,:,:) !contracted integrals from Gen1Int 2099 integer ierr !error information 2100 ! sets the arguments for Gen1Int (local) 2101 if (present(spher_gto)) then 2102 p_spher_gto = spher_gto 2103 else 2104 p_spher_gto = .true. 2105 end if 2106 if (present(order_geo_bra)) then 2107 p_order_geo_bra = order_geo_bra 2108 else 2109 p_order_geo_bra = 0 2110 end if 2111 if (present(order_geo_ket)) then 2112 p_order_geo_ket = order_geo_ket 2113 else 2114 p_order_geo_ket = 0 2115 end if 2116 if (present(order_geo_pot)) then 2117 p_order_geo_pot = order_geo_pot 2118 else 2119 p_order_geo_pot = 0 2120 end if 2121 if (present(order_geo_mom)) then 2122 p_order_geo_mom = order_geo_mom 2123 else 2124 p_order_geo_mom = 0 2125 end if 2126 if (present(nary_tree_total)) then 2127 p_num_cents = nary_tree_total%wt_node(nary_tree_total%order_geo) 2128 else 2129 p_num_cents = 0 2130 end if 2131 if (p_num_cents>0) then 2132 allocate(p_idx_cent(p_num_cents), stat=ierr) 2133 if (ierr/=0) & 2134 call error_stop("IntGetGAUPOT", "failed to allocate p_idx_cent", p_num_cents) 2135 allocate(p_order_cent(p_num_cents), stat=ierr) 2136 if (ierr/=0) & 2137 call error_stop("IntGetGAUPOT", "failed to allocate p_order_cent", p_num_cents) 2138 if (allocated(nary_tree_total%idx_atoms)) then 2139 do icent = 1, p_num_cents 2140 p_idx_cent(icent) = nary_tree_total%idx_atoms(nary_tree_total%idx_cent(icent)) 2141 end do 2142 else 2143 p_idx_cent = nary_tree_total%idx_cent(1:p_num_cents) 2144 end if 2145 p_order_cent = nary_tree_total%order_cent(1:p_num_cents) 2146 else 2147 allocate(p_idx_cent(1), stat=ierr) 2148 if (ierr/=0) call error_stop("IntGetGAUPOT", "failed to allocate p_idx_cent", 1) 2149 allocate(p_order_cent(1), stat=ierr) 2150 if (ierr/=0) call error_stop("IntGetGAUPOT", "failed to allocate p_order_cent", 1) 2151 p_idx_cent = 0 2152 p_order_cent = 0 2153 end if 2154 ! magnetic derivatives or derivatives with respect to rotational angular momentum 2155 if (order_mag_total>0 .or. order_mag_bra>0 .or. order_mag_ket>0 .or. & 2156 order_ram_total>0 .or. order_ram_bra>0 .or. order_ram_ket>0) then 2157 ! London atomic orbitals 2158 if (LondonAOUsed(info_LAO)) then 2159 if ((order_mag_total==0 .and. order_mom==0 .and. & 2160 ((order_mag_bra==1 .and. order_mag_ket==0) .or. & 2161 (order_mag_bra==0 .and. order_mag_ket==1)) .and. & 2162 order_elec==0 .and. p_order_geo_bra==0 .and. & 2163 p_order_geo_ket==0 .and. p_num_cents==0) .or. & 2164 (order_mag_total==1 .and. order_mom==0 .and. & 2165 order_mag_bra==0 .and. order_mag_ket==0 .and. & 2166 order_elec==0 .and. p_order_geo_bra==0 .and. & 2167 p_order_geo_ket==0 .and. p_num_cents==0)) then 2168 call LondonAOGetLPFOrigin(info_LAO, origin_London_PF) 2169 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 2170 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 2171 if (ierr/=0) & 2172 call error_stop("IntGetGAUPOT", "failed to allocate tmp_ints", & 2173 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 2174 ! SGTOs 2175 if (p_spher_gto) then 2176 if ((angular_bra==0 .and. angular_ket==0) .or. & 2177 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 2178 call contr_sgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2179 exponent_bra, num_contr_bra, contr_coef_bra, & 2180 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2181 exponent_ket, num_contr_ket, contr_coef_ket, & 2182 order_elec, idx_gauorg, gaupot_origin, & 2183 MAX_IDX_NON, origin_London_PF, & 2184 scal_const*0.5_REALK, 1, & 2185 p_order_geo_bra, p_order_geo_ket, & 2186 p_order_geo_pot, p_order_geo_mom, & 2187 p_num_cents, p_idx_cent, p_order_cent, & 2188 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 2189 ! reorders integrals if required 2190 else 2191 call contr_sgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2192 exponent_bra, num_contr_bra, contr_coef_bra, & 2193 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2194 exponent_ket, num_contr_ket, contr_coef_ket, & 2195 order_elec, idx_gauorg, gaupot_origin, & 2196 MAX_IDX_NON, origin_London_PF, & 2197 scal_const*0.5_REALK, 1, & 2198 p_order_geo_bra, p_order_geo_ket, & 2199 p_order_geo_pot, p_order_geo_mom, & 2200 p_num_cents, p_idx_cent, p_order_cent, & 2201 num_gto_bra, num_gto_ket, num_opt, contr_ints) 2202 if (angular_bra>0 .and. present(mag_num_bra) .and. & 2203 angular_ket>0 .and. present(mag_num_ket)) then 2204 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 2205 angular_ket, num_gto_ket, mag_num_ket, & 2206 num_contr_bra, num_contr_ket, num_opt, & 2207 contr_ints, tmp_ints) 2208 else if (angular_bra>0 .and. present(mag_num_bra)) then 2209 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 2210 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 2211 else if (angular_ket>0 .and. present(mag_num_ket)) then 2212 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 2213 num_gto_bra*num_contr_bra, num_contr_ket, & 2214 num_opt, contr_ints, tmp_ints) 2215 end if 2216 end if 2217 ! CGTOs 2218 else 2219 if ((angular_bra==0 .and. angular_ket==0) .or. & 2220 .not.(present(powers_bra) .or. present(powers_ket))) then 2221 call contr_cgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2222 exponent_bra, num_contr_bra, contr_coef_bra, & 2223 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2224 exponent_ket, num_contr_ket, contr_coef_ket, & 2225 order_elec, idx_gauorg, gaupot_origin, & 2226 MAX_IDX_NON, origin_London_PF, & 2227 scal_const*0.5_REALK, 1, & 2228 p_order_geo_bra, p_order_geo_ket, & 2229 p_order_geo_pot, p_order_geo_mom, & 2230 p_num_cents, p_idx_cent, p_order_cent, & 2231 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 2232 ! reorders integrals if required 2233 else 2234 call contr_cgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2235 exponent_bra, num_contr_bra, contr_coef_bra, & 2236 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2237 exponent_ket, num_contr_ket, contr_coef_ket, & 2238 order_elec, idx_gauorg, gaupot_origin, & 2239 MAX_IDX_NON, origin_London_PF, & 2240 scal_const*0.5_REALK, 1, & 2241 p_order_geo_bra, p_order_geo_ket, & 2242 p_order_geo_pot, p_order_geo_mom, & 2243 p_num_cents, p_idx_cent, p_order_cent, & 2244 num_gto_bra, num_gto_ket, num_opt, contr_ints) 2245 if (angular_bra>0 .and. present(powers_bra) .and. & 2246 angular_ket>0 .and. present(powers_ket)) then 2247 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 2248 angular_ket, num_gto_ket, powers_ket, & 2249 num_contr_bra, num_contr_ket, num_opt, & 2250 contr_ints, tmp_ints) 2251 else if (angular_bra>0 .and. present(powers_bra)) then 2252 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 2253 num_gto_ket*num_contr_ket*num_opt, contr_ints, tmp_ints) 2254 else if (angular_ket>0 .and. present(powers_ket)) then 2255 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 2256 num_gto_bra*num_contr_bra, num_contr_ket, & 2257 num_opt, contr_ints, tmp_ints) 2258 end if 2259 end if 2260 end if 2261 ! sets the first order partial magnetic derivatives 2262 call LondonAOGetGaugeOrigin(info_LAO, gauge_origin) 2263 if (order_mag_bra==1 .and. order_mag_ket==0) then 2264 do iopt = 1, num_opt-2, 3 2265 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) & 2266 - (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) 2267 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) & 2268 - (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) 2269 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) & 2270 - (coord_bra(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) 2271 end do 2272 else if (order_mag_bra==0 .and. order_mag_ket==1) then 2273 do iopt = 1, num_opt-2, 3 2274 contr_ints(:,:,:,:,iopt) = (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt+1) & 2275 - (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt+2) 2276 contr_ints(:,:,:,:,iopt+1) = (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+2) & 2277 - (coord_ket(3)-gauge_origin(3))*tmp_ints(:,:,:,:,iopt) 2278 contr_ints(:,:,:,:,iopt+2) = (coord_ket(2)-gauge_origin(2))*tmp_ints(:,:,:,:,iopt) & 2279 - (coord_ket(1)-gauge_origin(1))*tmp_ints(:,:,:,:,iopt+1) 2280 end do 2281 else 2282 do iopt = 1, num_opt-2, 3 2283 contr_ints(:,:,:,:,iopt) = (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt+2) & 2284 - (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt+1) 2285 contr_ints(:,:,:,:,iopt+1) = (coord_bra(3)-coord_ket(3))*tmp_ints(:,:,:,:,iopt) & 2286 - (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+2) 2287 contr_ints(:,:,:,:,iopt+2) = (coord_bra(1)-coord_ket(1))*tmp_ints(:,:,:,:,iopt+1) & 2288 - (coord_bra(2)-coord_ket(2))*tmp_ints(:,:,:,:,iopt) 2289 end do 2290 end if 2291 deallocate(tmp_ints) 2292 else 2293 call error_stop("IntGetGAUPOT", "LAO is not implemented", -1) 2294 end if 2295 ! non-LAOs 2296 else 2297 contr_ints = 0.0 2298 end if 2299 else 2300 ! SGTOs 2301 if (p_spher_gto) then 2302 if ((angular_bra==0 .and. angular_ket==0) .or. & 2303 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 2304 call contr_sgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2305 exponent_bra, num_contr_bra, contr_coef_bra, & 2306 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2307 exponent_ket, num_contr_ket, contr_coef_ket, & 2308 order_elec, idx_gauorg, gaupot_origin, & 2309 gaupot_expt, idx_diporg, dipole_origin, & 2310 scal_const, order_mom, & 2311 p_order_geo_bra, p_order_geo_ket, & 2312 p_order_geo_pot, p_order_geo_mom, & 2313 p_num_cents, p_idx_cent, p_order_cent, & 2314 num_gto_bra, num_gto_ket, num_opt, contr_ints) 2315 ! reorders integrals if required 2316 else 2317 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 2318 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 2319 if (ierr/=0) & 2320 call error_stop("IntGetGAUPOT", "failed to allocate tmp_ints", & 2321 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 2322 call contr_sgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2323 exponent_bra, num_contr_bra, contr_coef_bra, & 2324 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2325 exponent_ket, num_contr_ket, contr_coef_ket, & 2326 order_elec, idx_gauorg, gaupot_origin, & 2327 gaupot_expt, idx_diporg, dipole_origin, & 2328 scal_const, order_mom, & 2329 p_order_geo_bra, p_order_geo_ket, & 2330 p_order_geo_pot, p_order_geo_mom, & 2331 p_num_cents, p_idx_cent, p_order_cent, & 2332 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 2333 if (angular_bra>0 .and. present(mag_num_bra) .and. & 2334 angular_ket>0 .and. present(mag_num_ket)) then 2335 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 2336 angular_ket, num_gto_ket, mag_num_ket, & 2337 num_contr_bra, num_contr_ket, num_opt, & 2338 tmp_ints, contr_ints) 2339 else if (angular_bra>0 .and. present(mag_num_bra)) then 2340 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 2341 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 2342 else if (angular_ket>0 .and. present(mag_num_ket)) then 2343 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 2344 num_gto_bra*num_contr_bra, num_contr_ket, & 2345 num_opt, tmp_ints, contr_ints) 2346 end if 2347 deallocate(tmp_ints) 2348 end if 2349 ! CGTOs 2350 else 2351 if ((angular_bra==0 .and. angular_ket==0) .or. & 2352 .not.(present(powers_bra) .or. present(powers_ket))) then 2353 call contr_cgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2354 exponent_bra, num_contr_bra, contr_coef_bra, & 2355 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2356 exponent_ket, num_contr_ket, contr_coef_ket, & 2357 order_elec, idx_gauorg, gaupot_origin, & 2358 gaupot_expt, idx_diporg, dipole_origin, & 2359 scal_const, order_mom, & 2360 p_order_geo_bra, p_order_geo_ket, & 2361 p_order_geo_pot, p_order_geo_mom, & 2362 p_num_cents, p_idx_cent, p_order_cent, & 2363 num_gto_bra, num_gto_ket, num_opt, contr_ints) 2364 ! reorders integrals if required 2365 else 2366 allocate(tmp_ints(num_gto_bra,num_contr_bra, & 2367 num_gto_ket,num_contr_ket,num_opt), stat=ierr) 2368 if (ierr/=0) & 2369 call error_stop("IntGetGAUPOT", "failed to allocate tmp_ints", & 2370 num_gto_bra*num_contr_bra*num_gto_ket*num_contr_ket*num_opt) 2371 call contr_cgto_gaupot(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2372 exponent_bra, num_contr_bra, contr_coef_bra, & 2373 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2374 exponent_ket, num_contr_ket, contr_coef_ket, & 2375 order_elec, idx_gauorg, gaupot_origin, & 2376 gaupot_expt, idx_diporg, dipole_origin, & 2377 scal_const, order_mom, & 2378 p_order_geo_bra, p_order_geo_ket, & 2379 p_order_geo_pot, p_order_geo_mom, & 2380 p_num_cents, p_idx_cent, p_order_cent, & 2381 num_gto_bra, num_gto_ket, num_opt, tmp_ints) 2382 if (angular_bra>0 .and. present(powers_bra) .and. & 2383 angular_ket>0 .and. present(powers_ket)) then 2384 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 2385 angular_ket, num_gto_ket, powers_ket, & 2386 num_contr_bra, num_contr_ket, num_opt, & 2387 tmp_ints, contr_ints) 2388 else if (angular_bra>0 .and. present(powers_bra)) then 2389 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 2390 num_gto_ket*num_contr_ket*num_opt, tmp_ints, contr_ints) 2391 else if (angular_ket>0 .and. present(powers_ket)) then 2392 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 2393 num_gto_bra*num_contr_bra, num_contr_ket, & 2394 num_opt, tmp_ints, contr_ints) 2395 end if 2396 deallocate(tmp_ints) 2397 end if 2398 end if 2399 end if 2400 deallocate(p_idx_cent) 2401 deallocate(p_order_cent) 2402 end subroutine IntGetGAUPOT 2403 2404 !> \brief calculates the overlap distribution using contracted Gaussian type orbitals (GTOs) 2405 !> \author Bin Gao 2406 !> \date 2012-02-10 2407 !> \param idx_bra is the atomic index of bra center 2408 !> \param coord_bra contains the coordinates of bra center 2409 !> \param angular_bra is the angular number of bra center (s=0, p=1, d=2, ...) 2410 !> \param num_prim_bra is the number of primitive Gaussians of bra center 2411 !> \param exponent_bra contains the exponents of primitive Gaussians of bra center 2412 !> \param num_contr_bra is the number of contractions of bra center 2413 !> \param contr_coef_bra contains the contraction coefficients of bra center 2414 !> \param idx_ket is the atomic index of ket center 2415 !> \param coord_ket contains the coordinates of ket center 2416 !> \param angular_ket is the angular number of ket center 2417 !> \param num_prim_ket is the number of primitive Gaussians of ket center 2418 !> \param exponent_ket contains the exponents of primitive Gaussians of ket center 2419 !> \param num_contr_ket is the number of contractions of ket center 2420 !> \param contr_coef_ket contains the contraction coefficients of ket center 2421 !> \param spher_gto indicates if using spherical GTOs, otherwise Cartesian GTOs 2422 !> \param info_LAO contains the information of London atomic orbital 2423 !> \param order_mag_bra is the order of magnetic derivatives on bra center 2424 !> \param order_mag_ket is the order of magnetic derivatives on ket center 2425 !> \param order_mag_total is the order of total magnetic derivatives 2426 !> \param order_ram_bra is the order of derivatives w.r.t. total rotational angular momentum on bra center 2427 !> \param order_ram_ket is the order of derivatives w.r.t. total rotational angular momentum on ket center 2428 !> \param order_ram_total is the order of total derivatives w.r.t. total rotational angular momentum 2429 !> \param order_geo_bra is the order of geometric derivatives with respect to bra center 2430 !> \param order_geo_ket is the order of geometric derivatives with respect to ket center 2431 !> \param nary_tree_total contains the information of N-ary tree for total geometric derivatives 2432 !> \param num_points is the number of grid points 2433 !> \param grid_points contains the coordinates of grid points 2434 !> \param num_gto_bra is the number of spherical/Cartesian GTOs on bra center 2435 !> \param num_gto_ket is the number of spherical/Cartesian GTOs on ket center 2436 !> \param num_derv is the number of derivatives 2437 !> \param mag_num_bra contains the magnetic numbers of spherical GTOs on bra center 2438 !> \param mag_num_ket contains the magnetic numbers of spherical GTOs on ket center 2439 !> \param powers_bra contains the Cartesian powers of Cartesian GTOs on bra center 2440 !> \param powers_ket contains the Cartesian powers of Cartesian GTOs on ket center 2441 !> \return contr_odist contains the calculated contracted overlap distribution 2442 subroutine IntGetODST(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2443 exponent_bra, num_contr_bra, contr_coef_bra, & 2444 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2445 exponent_ket, num_contr_ket, contr_coef_ket, & 2446 spher_gto, info_LAO, & 2447 order_mag_bra, order_mag_ket, order_mag_total, & 2448 order_ram_bra, order_ram_ket, order_ram_total, & 2449 order_geo_bra, order_geo_ket, nary_tree_total, & 2450 num_points, grid_points, num_gto_bra, & 2451 num_gto_ket, num_derv, contr_odist, & 2452 mag_num_bra, mag_num_ket, powers_bra, powers_ket) 2453 integer, intent(in) :: idx_bra 2454 real(REALK), intent(in) :: coord_bra(3) 2455 integer, intent(in) :: angular_bra 2456 integer, intent(in) :: num_prim_bra 2457 real(REALK), intent(in) :: exponent_bra(num_prim_bra) 2458 integer, intent(in) :: num_contr_bra 2459 real(REALK), intent(in) :: contr_coef_bra(num_contr_bra,num_prim_bra) 2460 integer, intent(in) :: idx_ket 2461 real(REALK), intent(in) :: coord_ket(3) 2462 integer, intent(in) :: angular_ket 2463 integer, intent(in) :: num_prim_ket 2464 real(REALK), intent(in) :: exponent_ket(num_prim_ket) 2465 integer, intent(in) :: num_contr_ket 2466 real(REALK), intent(in) :: contr_coef_ket(num_contr_ket,num_prim_ket) 2467 logical, optional, intent(in) :: spher_gto 2468 type(london_ao_t), intent(in) :: info_LAO 2469 integer, intent(in) :: order_mag_bra 2470 integer, intent(in) :: order_mag_ket 2471 integer, intent(in) :: order_mag_total 2472 integer, intent(in) :: order_ram_bra 2473 integer, intent(in) :: order_ram_ket 2474 integer, intent(in) :: order_ram_total 2475 integer, optional, intent(in) :: order_geo_bra 2476 integer, optional, intent(in) :: order_geo_ket 2477 type(nary_tree_t), optional, intent(in) :: nary_tree_total 2478 integer, intent(in) :: num_points 2479 real(REALK), intent(in) :: grid_points(3,num_points) 2480 integer, intent(in) :: num_gto_bra 2481 integer, intent(in) :: num_gto_ket 2482 integer, intent(in) :: num_derv 2483 real(REALK), intent(out) :: contr_odist(num_gto_bra,num_contr_bra, & 2484 num_gto_ket,num_contr_ket, & 2485 num_points,num_derv) 2486 integer, optional, intent(in) :: mag_num_bra(num_gto_bra) 2487 integer, optional, intent(in) :: mag_num_ket(num_gto_ket) 2488 integer, optional, intent(in) :: powers_bra(3,num_gto_bra) 2489 integer, optional, intent(in) :: powers_ket(3,num_gto_ket) 2490 logical p_spher_gto !arguments for Gen1Int (local) 2491 integer p_order_geo_bra 2492 integer p_order_geo_ket 2493 integer p_num_cents 2494 integer, allocatable :: p_idx_cent(:) 2495 integer, allocatable :: p_order_cent(:) 2496 integer icent 2497 real(REALK), allocatable :: tmp_odist(:,:,:,:,:,:) !contracted overlap distribution from Gen1Int 2498 integer ierr !error information 2499 ! sets the arguments for Gen1Int (local) 2500 if (present(spher_gto)) then 2501 p_spher_gto = spher_gto 2502 else 2503 p_spher_gto = .true. 2504 end if 2505 if (present(order_geo_bra)) then 2506 p_order_geo_bra = order_geo_bra 2507 else 2508 p_order_geo_bra = 0 2509 end if 2510 if (present(order_geo_ket)) then 2511 p_order_geo_ket = order_geo_ket 2512 else 2513 p_order_geo_ket = 0 2514 end if 2515 if (present(nary_tree_total)) then 2516 p_num_cents = nary_tree_total%wt_node(nary_tree_total%order_geo) 2517 else 2518 p_num_cents = 0 2519 end if 2520 if (p_num_cents>0) then 2521 allocate(p_idx_cent(p_num_cents), stat=ierr) 2522 if (ierr/=0) & 2523 call error_stop("IntGetODST", "failed to allocate p_idx_cent", p_num_cents) 2524 allocate(p_order_cent(p_num_cents), stat=ierr) 2525 if (ierr/=0) & 2526 call error_stop("IntGetODST", "failed to allocate p_order_cent", p_num_cents) 2527 if (allocated(nary_tree_total%idx_atoms)) then 2528 do icent = 1, p_num_cents 2529 p_idx_cent(icent) = nary_tree_total%idx_atoms(nary_tree_total%idx_cent(icent)) 2530 end do 2531 else 2532 p_idx_cent = nary_tree_total%idx_cent(1:p_num_cents) 2533 end if 2534 p_order_cent = nary_tree_total%order_cent(1:p_num_cents) 2535 else 2536 allocate(p_idx_cent(1), stat=ierr) 2537 if (ierr/=0) call error_stop("IntGetODST", "failed to allocate p_idx_cent", 1) 2538 allocate(p_order_cent(1), stat=ierr) 2539 if (ierr/=0) call error_stop("IntGetODST", "failed to allocate p_order_cent", 1) 2540 p_idx_cent = 0 2541 p_order_cent = 0 2542 end if 2543 ! magnetic derivatives or derivatives with respect to rotational angular momentum 2544 if (order_mag_total>0 .or. order_mag_bra>0 .or. order_mag_ket>0 .or. & 2545 order_ram_total>0 .or. order_ram_bra>0 .or. order_ram_ket>0) then 2546 ! London atomic orbitals 2547 if (LondonAOUsed(info_LAO)) then 2548 call error_stop("IntGetODST", "LAO is not implemented", -1) 2549 ! non-LAOs 2550 else 2551 contr_odist = 0.0 2552 end if 2553 else 2554 ! SGTOs 2555 if (p_spher_gto) then 2556 if ((angular_bra==0 .and. angular_ket==0) .or. & 2557 .not.(present(mag_num_bra) .or. present(mag_num_ket))) then 2558 call contr_sgto_odist(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2559 exponent_bra, num_contr_bra, contr_coef_bra, & 2560 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2561 exponent_ket, num_contr_ket, contr_coef_ket, & 2562 p_order_geo_bra, p_order_geo_ket, & 2563 p_num_cents, p_idx_cent, p_order_cent, & 2564 num_points, grid_points, & 2565 num_gto_bra, num_gto_ket, num_derv, contr_odist) 2566 ! reorders integrals if required 2567 else 2568 allocate(tmp_odist(num_gto_bra,num_contr_bra, & 2569 num_gto_ket,num_contr_ket, & 2570 num_points,num_derv), stat=ierr) 2571 if (ierr/=0) & 2572 call error_stop("IntGetODST", "failed to allocate tmp_odist", & 2573 num_gto_bra*num_contr_bra*num_gto_ket & 2574 *num_contr_ket*num_points*num_derv) 2575 call contr_sgto_odist(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2576 exponent_bra, num_contr_bra, contr_coef_bra, & 2577 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2578 exponent_ket, num_contr_ket, contr_coef_ket, & 2579 p_order_geo_bra, p_order_geo_ket, & 2580 p_num_cents, p_idx_cent, p_order_cent, & 2581 num_points, grid_points, & 2582 num_gto_bra, num_gto_ket, num_derv, tmp_odist) 2583 if (angular_bra>0 .and. present(mag_num_bra) .and. & 2584 angular_ket>0 .and. present(mag_num_ket)) then 2585 call reorder_sgto_ints(angular_bra, num_gto_bra, mag_num_bra, & 2586 angular_ket, num_gto_ket, mag_num_ket, & 2587 num_contr_bra, num_contr_ket, & 2588 num_points*num_derv, tmp_odist, contr_odist) 2589 else if (angular_bra>0 .and. present(mag_num_bra)) then 2590 call reorder_sgtos(angular_bra, num_gto_bra, mag_num_bra, 1, num_contr_bra, & 2591 num_gto_ket*num_contr_ket*num_points*num_derv, & 2592 tmp_odist, contr_odist) 2593 else if (angular_ket>0 .and. present(mag_num_ket)) then 2594 call reorder_sgtos(angular_ket, num_gto_ket, mag_num_ket, & 2595 num_gto_bra*num_contr_bra, num_contr_ket, & 2596 num_points*num_derv, tmp_odist, contr_odist) 2597 end if 2598 deallocate(tmp_odist) 2599 end if 2600 ! CGTOs 2601 else 2602 if ((angular_bra==0 .and. angular_ket==0) .or. & 2603 .not.(present(powers_bra) .or. present(powers_ket))) then 2604 call contr_cgto_odist(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2605 exponent_bra, num_contr_bra, contr_coef_bra, & 2606 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2607 exponent_ket, num_contr_ket, contr_coef_ket, & 2608 p_order_geo_bra, p_order_geo_ket, & 2609 p_num_cents, p_idx_cent, p_order_cent, & 2610 num_points, grid_points, & 2611 num_gto_bra, num_gto_ket, num_derv, contr_odist) 2612 ! reorders integrals if required 2613 else 2614 allocate(tmp_odist(num_gto_bra,num_contr_bra, & 2615 num_gto_ket,num_contr_ket, & 2616 num_points,num_derv), stat=ierr) 2617 if (ierr/=0) & 2618 call error_stop("IntGetODST", "failed to allocate tmp_odist", & 2619 num_gto_bra*num_contr_bra*num_gto_ket & 2620 *num_contr_ket*num_points*num_derv) 2621 call contr_cgto_odist(idx_bra, coord_bra, angular_bra, num_prim_bra, & 2622 exponent_bra, num_contr_bra, contr_coef_bra, & 2623 idx_ket, coord_ket, angular_ket, num_prim_ket, & 2624 exponent_ket, num_contr_ket, contr_coef_ket, & 2625 p_order_geo_bra, p_order_geo_ket, & 2626 p_num_cents, p_idx_cent, p_order_cent, & 2627 num_points, grid_points, & 2628 num_gto_bra, num_gto_ket, num_derv, tmp_odist) 2629 if (angular_bra>0 .and. present(powers_bra) .and. & 2630 angular_ket>0 .and. present(powers_ket)) then 2631 call reorder_cgto_ints(angular_bra, num_gto_bra, powers_bra, & 2632 angular_ket, num_gto_ket, powers_ket, & 2633 num_contr_bra, num_contr_ket, & 2634 num_points*num_derv, tmp_odist, contr_odist) 2635 else if (angular_bra>0 .and. present(powers_bra)) then 2636 call reorder_cgtos(angular_bra, num_gto_bra, powers_bra, 1, num_contr_bra, & 2637 num_gto_ket*num_contr_ket*num_points*num_derv, tmp_odist, & 2638 contr_odist) 2639 else if (angular_ket>0 .and. present(powers_ket)) then 2640 call reorder_cgtos(angular_ket, num_gto_ket, powers_ket, & 2641 num_gto_bra*num_contr_bra, num_contr_ket, & 2642 num_points*num_derv, tmp_odist, contr_odist) 2643 end if 2644 deallocate(tmp_odist) 2645 end if 2646 end if 2647 end if 2648 deallocate(p_idx_cent) 2649 deallocate(p_order_cent) 2650 end subroutine IntGetODST 2651 2652end module gen1int_geom 2653