1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Automatic generation of auxiliary basis sets of different kind 8!> \author JGH 9!> 10!> <b>Modification history:</b> 11!> - 11.2017 creation [JGH] 12! ************************************************************************************************** 13MODULE auto_basis 14 USE aux_basis_set, ONLY: create_aux_basis 15 USE basis_set_types, ONLY: get_gto_basis_set,& 16 gto_basis_set_type 17 USE bibliography, ONLY: Stoychev2016,& 18 cite_reference 19 USE kinds, ONLY: default_string_length,& 20 dp 21 USE mathconstants, ONLY: dfac,& 22 fac,& 23 gamma1,& 24 pi,& 25 rootpi 26 USE orbital_pointers, ONLY: init_orbital_pointers 27 USE periodic_table, ONLY: get_ptable_info 28 USE powell, ONLY: opt_state_type,& 29 powell_optimize 30 USE qs_kind_types, ONLY: get_qs_kind,& 31 qs_kind_type 32#include "./base/base_uses.f90" 33 34 IMPLICIT NONE 35 36 PRIVATE 37 38 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'auto_basis' 39 40 PUBLIC :: create_ri_aux_basis_set, create_aux_fit_basis_set, create_lri_aux_basis_set 41 42CONTAINS 43 44! ************************************************************************************************** 45!> \brief Create a RI_AUX basis set using some heuristics 46!> \param ri_aux_basis_set ... 47!> \param qs_kind ... 48!> \param basis_cntrl ... 49!> \date 01.11.2017 50!> \author JGH 51! ************************************************************************************************** 52 SUBROUTINE create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, basis_cntrl) 53 TYPE(gto_basis_set_type), POINTER :: ri_aux_basis_set 54 TYPE(qs_kind_type), INTENT(IN) :: qs_kind 55 INTEGER, INTENT(IN) :: basis_cntrl 56 57 CHARACTER(len=*), PARAMETER :: routineN = 'create_ri_aux_basis_set', & 58 routineP = moduleN//':'//routineN 59 60 CHARACTER(LEN=2) :: element_symbol 61 CHARACTER(LEN=default_string_length) :: bsname 62 INTEGER :: i, j, jj, l, laux, linc, lmax, lval, lx, & 63 nsets, nx, z 64 INTEGER, DIMENSION(0:18) :: nval 65 INTEGER, DIMENSION(0:9, 1:20) :: nl 66 INTEGER, DIMENSION(1:3) :: ls1, ls2, npgf 67 INTEGER, DIMENSION(:), POINTER :: econf 68 REAL(KIND=dp) :: xv, zval 69 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zet 70 REAL(KIND=dp), DIMENSION(0:18) :: bv, bval, fv, peff, pend, pmax, pmin 71 REAL(KIND=dp), DIMENSION(0:9) :: zeff, zmax, zmin 72 REAL(KIND=dp), DIMENSION(3) :: amax, amin, bmin 73 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 74 75 ! 76 CALL cite_reference(Stoychev2016) 77 ! 78 bv(0:18) = (/1.8_dp, 2.0_dp, 2.2_dp, 2.2_dp, 2.3_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, & 79 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp, 3.0_dp/) 80 fv(0:18) = (/20.0_dp, 4.0_dp, 4.0_dp, 3.5_dp, 2.5_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, & 81 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 2.0_dp/) 82 ! 83 CPASSERT(.NOT. ASSOCIATED(ri_aux_basis_set)) 84 NULLIFY (orb_basis_set) 85 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB") 86 IF (ASSOCIATED(orb_basis_set)) THEN 87 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff) 88 !Note: RI basis coud require lmax up to 2*orb_lmax. This ensures that all orbital pointers 89 ! are properly initialized before building the basis 90 CALL init_orbital_pointers(2*lmax) 91 CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff) 92 CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol) 93 CALL get_ptable_info(element_symbol, ielement=z) 94 lval = 0 95 DO l = 0, MAXVAL(UBOUND(econf)) 96 IF (econf(l) > 0) lval = l 97 END DO 98 IF (SUM(econf) /= NINT(zval)) THEN 99 CPWARN("Valence charge and electron configuration not consistent") 100 END IF 101 pend = 0.0_dp 102 linc = 1 103 IF (z > 18) linc = 2 104 SELECT CASE (basis_cntrl) 105 CASE (0) 106 laux = MAX(2*lval, lmax + linc) 107 CASE (1) 108 laux = MAX(2*lval, lmax + linc) 109 CASE (2) 110 laux = MAX(2*lval, lmax + linc + 1) 111 CASE (3) 112 laux = MAX(2*lmax, lmax + linc + 2) 113 CASE DEFAULT 114 CPABORT("Invalid value of control variable") 115 END SELECT 116 ! 117 DO l = 2*lmax + 1, laux 118 xv = peff(2*lmax) 119 pmin(l) = xv 120 pmax(l) = xv 121 peff(l) = xv 122 pend(l) = xv 123 END DO 124 ! 125 DO l = 0, laux 126 IF (l <= 2*lval) THEN 127 pend(l) = MIN(fv(l)*peff(l), pmax(l)) 128 bval(l) = 1.8_dp 129 ELSE 130 pend(l) = peff(l) 131 bval(l) = bv(l) 132 END IF 133 xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp 134 nval(l) = MAX(CEILING(xv), 0) 135 END DO 136 ! first set include valence only 137 nsets = 1 138 ls1(1) = 0 139 ls2(1) = lval 140 DO l = lval + 1, laux 141 IF (nval(l) < nval(lval) - 1) EXIT 142 ls2(1) = l 143 END DO 144 ! second set up to 2*lval 145 IF (laux > ls2(1)) THEN 146 IF (lval == 0 .OR. 2*lval <= ls2(1) + 1) THEN 147 nsets = 2 148 ls1(2) = ls2(1) + 1 149 ls2(2) = laux 150 ELSE 151 nsets = 2 152 ls1(2) = ls2(1) + 1 153 ls2(2) = MIN(2*lval, laux) 154 lx = ls2(2) 155 DO l = lx + 1, laux 156 IF (nval(l) < nval(lx) - 1) EXIT 157 ls2(2) = l 158 END DO 159 IF (laux > ls2(2)) THEN 160 nsets = 3 161 ls1(3) = ls2(2) + 1 162 ls2(3) = laux 163 END IF 164 END IF 165 END IF 166 ! 167 amax = 0.0 168 amin = HUGE(0.0_dp) 169 bmin = HUGE(0.0_dp) 170 DO i = 1, nsets 171 DO j = ls1(i), ls2(i) 172 amax(i) = MAX(amax(i), pend(j)) 173 amin(i) = MIN(amin(i), pmin(j)) 174 bmin(i) = MIN(bmin(i), bval(j)) 175 END DO 176 xv = LOG(amax(i)/amin(i))/LOG(bmin(i)) + 1.e-10_dp 177 npgf(i) = MAX(CEILING(xv), 0) 178 END DO 179 nx = MAXVAL(npgf(1:nsets)) 180 ALLOCATE (zet(nx, nsets)) 181 zet = 0.0_dp 182 nl = 0 183 DO i = 1, nsets 184 DO j = 1, npgf(i) 185 jj = npgf(i) - j + 1 186 zet(jj, i) = amin(i)*bmin(i)**(j - 1) 187 END DO 188 DO l = ls1(i), ls2(i) 189 nl(l, i) = nval(l) 190 END DO 191 END DO 192 bsname = TRIM(element_symbol)//"-RI-AUX-"//TRIM(orb_basis_set%name) 193 ! 194 CALL create_aux_basis(ri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet) 195 ! 196 DEALLOCATE (zet) 197 END IF 198 199 END SUBROUTINE create_ri_aux_basis_set 200! ************************************************************************************************** 201!> \brief Create a LRI_AUX basis set using some heuristics 202!> \param lri_aux_basis_set ... 203!> \param qs_kind ... 204!> \param basis_cntrl ... 205!> \param exact_1c_terms ... 206!> \date 01.11.2017 207!> \author JGH 208! ************************************************************************************************** 209 SUBROUTINE create_lri_aux_basis_set(lri_aux_basis_set, qs_kind, basis_cntrl, exact_1c_terms) 210 TYPE(gto_basis_set_type), POINTER :: lri_aux_basis_set 211 TYPE(qs_kind_type), INTENT(IN) :: qs_kind 212 INTEGER, INTENT(IN) :: basis_cntrl 213 LOGICAL, INTENT(IN) :: exact_1c_terms 214 215 CHARACTER(len=*), PARAMETER :: routineN = 'create_lri_aux_basis_set', & 216 routineP = moduleN//':'//routineN 217 218 CHARACTER(LEN=2) :: element_symbol 219 CHARACTER(LEN=default_string_length) :: bsname 220 INTEGER :: i, j, l, laux, linc, lm, lmax, lval, n1, & 221 n2, nsets, z 222 INTEGER, DIMENSION(0:18) :: nval 223 INTEGER, DIMENSION(0:9, 1:50) :: nl 224 INTEGER, DIMENSION(1:50) :: ls1, ls2, npgf 225 INTEGER, DIMENSION(:), POINTER :: econf 226 REAL(KIND=dp) :: xv, zval 227 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zet 228 REAL(KIND=dp), DIMENSION(0:18) :: bval, peff, pend, pmax, pmin 229 REAL(KIND=dp), DIMENSION(0:9) :: zeff, zmax, zmin 230 REAL(KIND=dp), DIMENSION(4) :: bv, bx 231 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 232 233 ! 234 bv(1:4) = (/2.00_dp, 1.90_dp, 1.80_dp, 1.80_dp/) 235 bx(1:4) = (/2.60_dp, 2.40_dp, 2.20_dp, 2.20_dp/) 236 ! 237 CPASSERT(.NOT. ASSOCIATED(lri_aux_basis_set)) 238 NULLIFY (orb_basis_set) 239 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB") 240 IF (ASSOCIATED(orb_basis_set)) THEN 241 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff) 242 CALL get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff) 243 CALL get_qs_kind(qs_kind, zeff=zval, elec_conf=econf, element_symbol=element_symbol) 244 CALL get_ptable_info(element_symbol, ielement=z) 245 lval = 0 246 DO l = 0, MAXVAL(UBOUND(econf)) 247 IF (econf(l) > 0) lval = l 248 END DO 249 IF (SUM(econf) /= NINT(zval)) THEN 250 CPWARN("Valence charge and electron configuration not consistent") 251 END IF 252 ! 253 pend = 0.0_dp 254 linc = 1 255 IF (z > 18) linc = 2 256 SELECT CASE (basis_cntrl) 257 CASE (0) 258 laux = MAX(2*lval, lmax + linc) 259 laux = MIN(laux, 2 + linc) 260 CASE (1) 261 laux = MAX(2*lval, lmax + linc) 262 laux = MIN(laux, 3 + linc) 263 CASE (2) 264 laux = MAX(2*lval, lmax + linc + 1) 265 laux = MIN(laux, 4 + linc) 266 CASE (3) 267 laux = MAX(2*lval, lmax + linc + 1) 268 laux = MIN(laux, 4 + linc) 269 CASE DEFAULT 270 CPABORT("Invalid value of control variable") 271 END SELECT 272 ! 273 DO l = 2*lmax + 1, laux 274 pmin(l) = pmin(2*lmax) 275 pmax(l) = pmax(2*lmax) 276 peff(l) = peff(2*lmax) 277 END DO 278 ! 279 IF (exact_1c_terms) THEN 280 DO l = 0, laux 281 IF (l <= lval + 1) THEN 282 pend(l) = zmax(l) + 1.0_dp 283 bval(l) = bv(basis_cntrl + 1) 284 ELSE 285 pend(l) = 2.0_dp*peff(l) 286 bval(l) = bx(basis_cntrl + 1) 287 END IF 288 pmin(l) = zmin(l) 289 xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp 290 nval(l) = MAX(CEILING(xv), 0) 291 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l)) 292 END DO 293 ELSE 294 DO l = 0, laux 295 IF (l <= lval + 1) THEN 296 pend(l) = pmax(l) 297 bval(l) = bv(basis_cntrl + 1) 298 pmin(l) = zmin(l) 299 ELSE 300 pend(l) = 4.0_dp*peff(l) 301 bval(l) = bx(basis_cntrl + 1) 302 END IF 303 xv = LOG(pend(l)/pmin(l))/LOG(bval(l)) + 1.e-10_dp 304 nval(l) = MAX(CEILING(xv), 0) 305 bval(l) = (pend(l)/pmin(l))**(1._dp/nval(l)) 306 END DO 307 END IF 308 ! 309 lm = MIN(2*lval, 3) 310 n1 = MAXVAL(nval(0:lm)) 311 n2 = MAXVAL(nval(lm + 1:laux)) 312 nsets = n1 + n2 313 ALLOCATE (zet(1, nsets)) 314 zet = 0.0_dp 315 nl = 0 316 j = MAXVAL(MAXLOC(nval(0:lm))) 317 DO i = 1, n1 318 ls1(i) = 0 319 ls2(i) = lm 320 npgf(i) = 1 321 zet(1, i) = pmin(j)*bval(j)**(i - 1) 322 DO l = 0, lm 323 nl(l, i) = 1 324 END DO 325 END DO 326 j = lm + 1 327 DO i = n1 + 1, nsets 328 ls1(i) = lm + 1 329 ls2(i) = laux 330 npgf(i) = 1 331 zet(1, i) = pmin(j)*bval(j)**(i - n1 - 1) 332 DO l = lm + 1, laux 333 nl(l, i) = 1 334 END DO 335 END DO 336 ! 337 bsname = TRIM(element_symbol)//"-LRI-AUX-"//TRIM(orb_basis_set%name) 338 ! 339 CALL create_aux_basis(lri_aux_basis_set, bsname, nsets, ls1, ls2, nl, npgf, zet) 340 ! 341 DEALLOCATE (zet) 342 END IF 343 344 END SUBROUTINE create_lri_aux_basis_set 345! ************************************************************************************************** 346!> \brief Create a AUX_FIT basis set using some heuristics 347!> \param aux_fit_basis ... 348!> \param qs_kind ... 349!> \param basis_cntrl ... 350!> \date 01.11.2017 351!> \author JGH 352! ************************************************************************************************** 353 SUBROUTINE create_aux_fit_basis_set(aux_fit_basis, qs_kind, basis_cntrl) 354 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis 355 TYPE(qs_kind_type), INTENT(IN) :: qs_kind 356 INTEGER, INTENT(IN) :: basis_cntrl 357 358 CHARACTER(len=*), PARAMETER :: routineN = 'create_aux_fit_basis_set', & 359 routineP = moduleN//':'//routineN 360 361 CHARACTER(LEN=2) :: element_symbol 362 CHARACTER(LEN=default_string_length) :: bsname 363 INTEGER :: i, iset, l, laux, lmax, lval, maxpgf, & 364 mx, nf, np, nsets, nx, z 365 INTEGER, DIMENSION(0:9) :: nfun, nval 366 INTEGER, DIMENSION(0:9, 1:20) :: nl 367 INTEGER, DIMENSION(1:20) :: lset, npgf 368 INTEGER, DIMENSION(:), POINTER :: econf 369 REAL(KIND=dp) :: amet, z1, z2, zvel 370 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zval 371 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gcval, zet 372 REAL(KIND=dp), DIMENSION(0:9) :: zeff, zmax, zmin 373 REAL(KIND=dp), DIMENSION(25) :: afit, xval 374 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 375 TYPE(opt_state_type) :: ostate 376 377 ! 378 CPABORT("Automatic basis set generation not activated") 379 ! 380 CPASSERT(.NOT. ASSOCIATED(aux_fit_basis)) 381 NULLIFY (orb_basis_set) 382 CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, basis_type="ORB") 383 IF (ASSOCIATED(orb_basis_set)) THEN 384 CALL get_basis_keyfigures(orb_basis_set, lmax, zmin, zmax, zeff) 385 CALL get_qs_kind(qs_kind, zeff=zvel, elec_conf=econf, element_symbol=element_symbol) 386 CALL get_ptable_info(element_symbol, ielement=z) 387 lval = 0 388 DO l = 0, MAXVAL(UBOUND(econf)) 389 IF (econf(l) > 0) lval = l 390 END DO 391 IF (SUM(econf) /= NINT(zvel)) THEN 392 CPWARN("Valence charge and electron configuration not consistent") 393 END IF 394 nval = 0 395 DO l = 0, lval 396 nx = econf(l) 397 DO 398 IF (nx > 0) THEN 399 nval(l) = nval(l) + 1 400 nx = nx - 2*(2*l + 1) 401 ELSE 402 EXIT 403 END IF 404 END DO 405 END DO 406 nfun = nval 407 SELECT CASE (basis_cntrl) 408 CASE (0) 409 laux = lval 410 DO l = 0, lval 411 nfun(l) = nfun(l) + 1 412 END DO 413 CASE (1) 414 laux = MIN(lval + 1, lmax) 415 DO l = 0, lval 416 nfun(l) = nfun(l) + 1 417 END DO 418 IF (laux > lval) nfun(laux) = 1 419 CASE (2) 420 laux = MIN(lval + 1, lmax) 421 DO l = 0, lval 422 nfun(l) = nfun(l) + 2 423 END DO 424 IF (laux > lval) nfun(laux) = 1 425 CASE (3) 426 laux = MIN(lval + 2, lmax) 427 DO l = 0, lval 428 nfun(l) = nfun(l) + 3 429 END DO 430 IF (laux > lval) nfun(lval + 1) = 2 431 IF (laux > lval + 1) nfun(laux) = 1 432 CASE DEFAULT 433 CPABORT("Invalid value of control variable") 434 END SELECT 435 ! 436 nsets = 0 437 maxpgf = 0 438 DO l = 0, lval 439 z1 = MAX(zmin(l), 0.10_dp + l*0.025_dp) 440 z2 = zmax(l) 441 mx = CEILING(LOG(z2/z1)/LOG(5.0_dp)) 442 IF (nval(l) > 1) THEN 443 nsets = nsets + 2 444 maxpgf = MAX(maxpgf, nval(l), mx, 3) 445 ELSEIF (nval(l) == 1) THEN 446 nsets = nsets + 1 447 maxpgf = MAX(maxpgf, mx, 3) 448 END IF 449 DO i = nval(l) + 1, nfun(l) 450 maxpgf = MAX(maxpgf, mx, 1) 451 nsets = nsets + 1 452 END DO 453 END DO 454 DO l = lval + 1, laux 455 maxpgf = MAX(maxpgf, 1) 456 nsets = nsets + nfun(l) 457 END DO 458 ! 459 ALLOCATE (zet(maxpgf, nsets)) 460 zet = 0.0_dp 461 nl = 0 462 iset = 0 463 DO l = 0, laux 464 amet = 2.50_dp 465 ! optimize exponensts 466 z1 = MAX(zmin(l), 0.20_dp + l*0.025_dp) 467 z2 = zmax(l) 468 mx = CEILING(LOG(z2/z1)/LOG(4.0_dp)) 469 IF (nval(l) > 1) THEN 470 nx = MAX(nfun(l), mx, 3) 471 ELSE IF (nval(l) == 1) THEN 472 nx = MAX(nfun(l), mx, 2) 473 ELSE 474 nx = nfun(l) 475 END IF 476 CPASSERT(nx <= 25) 477 ! Powell optimizer 478 CALL get_basis_functions(orb_basis_set, l, np, nf, zval, gcval) 479 ostate%nf = 0 480 ostate%nvar = nx 481 xval(:) = 5.00_dp 482 xval(1) = z1 483 ostate%rhoend = 1.e-8_dp 484 ostate%rhobeg = 5.e-1_dp 485 ostate%maxfun = 1000 486 ostate%iprint = 0 487 ostate%unit = 0 488 ostate%state = 0 489 DO 490 IF (ostate%state == 2) THEN 491 afit(1) = xval(1) 492 DO i = 2, nx 493 afit(i) = afit(i - 1)*xval(i) 494 END DO 495 CALL overlap_maximum(l, np, nf, zval, gcval, nx, afit, amet, ostate%f) 496 CALL neb_potential(xval, ostate%nvar, ostate%f) 497 END IF 498 IF (ostate%state == -1) EXIT 499 CALL powell_optimize(ostate%nvar, xval, ostate) 500 END DO 501 ostate%state = 8 502 CALL powell_optimize(ostate%nvar, xval, ostate) 503 afit(1) = xval(1) 504 DO i = 2, nx 505 afit(i) = afit(i - 1)*xval(i) 506 END DO 507 DEALLOCATE (zval, gcval) 508 ! 509 IF (nval(l) > 1) THEN 510 ! split set 511 iset = iset + 1 512 lset(iset) = l 513 npgf(iset) = nx - 1 514 nl(l, iset) = nval(l) - 1 515 zet(1:nx - 1, iset) = afit(1:nx - 1) 516 ! new set 517 iset = iset + 1 518 lset(iset) = l 519 npgf(iset) = nx - 1 520 nl(l, iset) = 1 521 zet(1:nx - 1, iset) = afit(2:nx) 522 DO i = 1, nfun(l) - 2 523 iset = iset + 1 524 lset(iset) = l 525 npgf(iset) = 1 526 zet(1, iset) = afit(nx - i + 1) 527 nl(l, iset) = 1 528 END DO 529 ELSEIF (nval(l) == 1) THEN 530 iset = iset + 1 531 lset(iset) = l 532 npgf(iset) = nx 533 zet(1:nx, iset) = afit(1:nx) 534 nl(l, iset) = 1 535 DO i = 1, nfun(l) - 1 536 iset = iset + 1 537 lset(iset) = l 538 npgf(iset) = 1 539 zet(1, iset) = afit(nx - i + 1) 540 nl(l, iset) = 1 541 END DO 542 ELSE 543 DO i = 1, nfun(l) 544 iset = iset + 1 545 lset(iset) = l 546 npgf(iset) = 1 547 zet(1, iset) = afit(i) 548 nl(l, iset) = 1 549 END DO 550 END IF 551 END DO 552 ! 553 bsname = TRIM(element_symbol)//"-AUX-FIT-"//TRIM(orb_basis_set%name) 554 ! 555 CALL create_aux_basis(aux_fit_basis, bsname, nsets, lset, lset, nl, npgf, zet) 556 ! 557 DEALLOCATE (zet) 558 ! 559 END IF 560 561 END SUBROUTINE create_aux_fit_basis_set 562! ************************************************************************************************** 563!> \brief ... 564!> \param basis_set ... 565!> \param lmax ... 566!> \param zmin ... 567!> \param zmax ... 568!> \param zeff ... 569! ************************************************************************************************** 570 SUBROUTINE get_basis_keyfigures(basis_set, lmax, zmin, zmax, zeff) 571 TYPE(gto_basis_set_type), POINTER :: basis_set 572 INTEGER, INTENT(OUT) :: lmax 573 REAL(KIND=dp), DIMENSION(0:9), INTENT(OUT) :: zmin, zmax, zeff 574 575 CHARACTER(len=*), PARAMETER :: routineN = 'get_basis_keyfigures', & 576 routineP = moduleN//':'//routineN 577 578 INTEGER :: i, ipgf, iset, ishell, j, l, nset 579 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell 580 INTEGER, DIMENSION(:, :), POINTER :: lshell 581 REAL(KIND=dp) :: aeff, gcca, gccb, kval, rexp, rint, rno, & 582 zeta 583 REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet 584 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc 585 586 CALL get_gto_basis_set(gto_basis_set=basis_set, & 587 nset=nset, & 588 nshell=nshell, & 589 npgf=npgf, & 590 l=lshell, & 591 lmax=lm, & 592 zet=zet, & 593 gcc=gcc) 594 595 lmax = MAXVAL(lm) 596 CPASSERT(lmax <= 9) 597 598 zmax = 0.0_dp 599 zmin = HUGE(0.0_dp) 600 zeff = 0.0_dp 601 602 DO iset = 1, nset 603 ! zmin zmax 604 DO ipgf = 1, npgf(iset) 605 DO ishell = 1, nshell(iset) 606 l = lshell(ishell, iset) 607 zeta = zet(ipgf, iset) 608 zmax(l) = MAX(zmax(l), zeta) 609 zmin(l) = MIN(zmin(l), zeta) 610 END DO 611 END DO 612 ! zeff 613 DO ishell = 1, nshell(iset) 614 l = lshell(ishell, iset) 615 kval = fac(l + 1)**2*2._dp**(2*l + 1)/fac(2*l + 2) 616 rexp = 0.0_dp 617 rno = 0.0_dp 618 DO i = 1, npgf(iset) 619 gcca = gcc(i, ishell, iset) 620 DO j = 1, npgf(iset) 621 zeta = zet(i, iset) + zet(j, iset) 622 gccb = gcc(j, ishell, iset) 623 rint = 0.5_dp*fac(l + 1)/zeta**(l + 2) 624 rexp = rexp + gcca*gccb*rint 625 rint = rootpi*0.5_dp**(l + 2)*dfac(2*l + 1)/zeta**(l + 1.5_dp) 626 rno = rno + gcca*gccb*rint 627 END DO 628 END DO 629 rexp = rexp/rno 630 aeff = (fac(l + 1)/dfac(2*l + 1))**2*2._dp**(2*l + 1)/(pi*rexp**2) 631 zeff(l) = MAX(zeff(l), aeff) 632 END DO 633 END DO 634 635 END SUBROUTINE get_basis_keyfigures 636 637! ************************************************************************************************** 638!> \brief ... 639!> \param lmax ... 640!> \param zmin ... 641!> \param zmax ... 642!> \param zeff ... 643!> \param pmin ... 644!> \param pmax ... 645!> \param peff ... 646! ************************************************************************************************** 647 SUBROUTINE get_basis_products(lmax, zmin, zmax, zeff, pmin, pmax, peff) 648 INTEGER, INTENT(IN) :: lmax 649 REAL(KIND=dp), DIMENSION(0:9), INTENT(IN) :: zmin, zmax, zeff 650 REAL(KIND=dp), DIMENSION(0:18), INTENT(OUT) :: pmin, pmax, peff 651 652 CHARACTER(len=*), PARAMETER :: routineN = 'get_basis_products', & 653 routineP = moduleN//':'//routineN 654 655 INTEGER :: l1, l2, la 656 657 pmin = HUGE(0.0_dp) 658 pmax = 0.0_dp 659 peff = 0.0_dp 660 661 DO l1 = 0, lmax 662 DO l2 = l1, lmax 663 DO la = l2 - l1, l2 + l1 664 pmax(la) = MAX(pmax(la), zmax(l1) + zmax(l2)) 665 pmin(la) = MIN(pmin(la), zmin(l1) + zmin(l2)) 666 peff(la) = MAX(peff(la), zeff(l1) + zeff(l2)) 667 END DO 668 END DO 669 END DO 670 671 END SUBROUTINE get_basis_products 672! ************************************************************************************************** 673!> \brief ... 674!> \param lm ... 675!> \param npgf ... 676!> \param nfun ... 677!> \param zet ... 678!> \param gcc ... 679!> \param nfit ... 680!> \param afit ... 681!> \param amet ... 682!> \param eval ... 683! ************************************************************************************************** 684 SUBROUTINE overlap_maximum(lm, npgf, nfun, zet, gcc, nfit, afit, amet, eval) 685 INTEGER, INTENT(IN) :: lm, npgf, nfun 686 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zet 687 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: gcc 688 INTEGER, INTENT(IN) :: nfit 689 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: afit 690 REAL(KIND=dp), INTENT(IN) :: amet 691 REAL(KIND=dp), INTENT(OUT) :: eval 692 693 CHARACTER(len=*), PARAMETER :: routineN = 'overlap_maximum', & 694 routineP = moduleN//':'//routineN 695 696 INTEGER :: i, ia, ib, info 697 REAL(KIND=dp) :: fij, fxij, intab, p, xij 698 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fx, tx, x2, xx 699 700 ! SUM_i(fi M fi) 701 fij = 0.0_dp 702 DO ia = 1, npgf 703 DO ib = 1, npgf 704 p = zet(ia) + zet(ib) + amet 705 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1) 706 DO i = 1, nfun 707 fij = fij + gcc(ia, i)*gcc(ib, i)*intab 708 END DO 709 END DO 710 END DO 711 712 !Integrals (fi M xj) 713 ALLOCATE (fx(nfit, nfun), tx(nfit, nfun)) 714 fx = 0.0_dp 715 DO ia = 1, npgf 716 DO ib = 1, nfit 717 p = zet(ia) + afit(ib) + amet 718 intab = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1) 719 DO i = 1, nfun 720 fx(ib, i) = fx(ib, i) + gcc(ia, i)*intab 721 END DO 722 END DO 723 END DO 724 725 !Integrals (xi M xj) 726 ALLOCATE (xx(nfit, nfit), x2(nfit, nfit)) 727 DO ia = 1, nfit 728 DO ib = 1, nfit 729 p = afit(ia) + afit(ib) + amet 730 xx(ia, ib) = 0.5_dp/p**(lm + 1.5_dp)*gamma1(lm + 1) 731 END DO 732 END DO 733 734 !Solve for tab 735 tx(1:nfit, 1:nfun) = fx(1:nfit, 1:nfun) 736 x2(1:nfit, 1:nfit) = xx(1:nfit, 1:nfit) 737 CALL DPOSV("U", nfit, nfun, x2, nfit, tx, nfit, info) 738 IF (info == 0) THEN 739 ! value t*xx*t 740 xij = 0.0_dp 741 DO i = 1, nfun 742 xij = xij + DOT_PRODUCT(tx(:, i), MATMUL(xx, tx(:, i))) 743 END DO 744 ! value t*fx 745 fxij = 0.0_dp 746 DO i = 1, nfun 747 fxij = fxij + DOT_PRODUCT(tx(:, i), fx(:, i)) 748 END DO 749 ! 750 eval = fij - 2.0_dp*fxij + xij 751 ELSE 752 ! error in solving for max overlap 753 eval = 1.0e10_dp 754 END IF 755 756 DEALLOCATE (fx, xx, x2, tx) 757 758 END SUBROUTINE overlap_maximum 759! ************************************************************************************************** 760!> \brief ... 761!> \param x ... 762!> \param n ... 763!> \param eval ... 764! ************************************************************************************************** 765 SUBROUTINE neb_potential(x, n, eval) 766 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x 767 INTEGER, INTENT(IN) :: n 768 REAL(KIND=dp), INTENT(INOUT) :: eval 769 770 CHARACTER(len=*), PARAMETER :: routineN = 'neb_potential', routineP = moduleN//':'//routineN 771 772 INTEGER :: i 773 774 DO i = 2, n 775 IF (x(i) < 1.5_dp) THEN 776 eval = eval + 10.0_dp*(1.5_dp - x(i))**2 777 END IF 778 END DO 779 780 END SUBROUTINE neb_potential 781! ************************************************************************************************** 782!> \brief ... 783!> \param basis_set ... 784!> \param lin ... 785!> \param np ... 786!> \param nf ... 787!> \param zval ... 788!> \param gcval ... 789! ************************************************************************************************** 790 SUBROUTINE get_basis_functions(basis_set, lin, np, nf, zval, gcval) 791 TYPE(gto_basis_set_type), POINTER :: basis_set 792 INTEGER, INTENT(IN) :: lin 793 INTEGER, INTENT(OUT) :: np, nf 794 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zval 795 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gcval 796 797 CHARACTER(len=*), PARAMETER :: routineN = 'get_basis_functions', & 798 routineP = moduleN//':'//routineN 799 800 INTEGER :: iset, ishell, j1, j2, jf, jp, l, nset 801 INTEGER, DIMENSION(:), POINTER :: lm, npgf, nshell 802 INTEGER, DIMENSION(:, :), POINTER :: lshell 803 LOGICAL :: toadd 804 REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet 805 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc 806 807 CALL get_gto_basis_set(gto_basis_set=basis_set, & 808 nset=nset, & 809 nshell=nshell, & 810 npgf=npgf, & 811 l=lshell, & 812 lmax=lm, & 813 zet=zet, & 814 gcc=gcc) 815 816 np = 0 817 nf = 0 818 DO iset = 1, nset 819 toadd = .TRUE. 820 DO ishell = 1, nshell(iset) 821 l = lshell(ishell, iset) 822 IF (l == lin) THEN 823 nf = nf + 1 824 IF (toadd) THEN 825 np = np + npgf(iset) 826 toadd = .FALSE. 827 END IF 828 END IF 829 END DO 830 END DO 831 ALLOCATE (zval(np), gcval(np, nf)) 832 zval = 0.0_dp 833 gcval = 0.0_dp 834 ! 835 jp = 0 836 jf = 0 837 DO iset = 1, nset 838 toadd = .TRUE. 839 DO ishell = 1, nshell(iset) 840 l = lshell(ishell, iset) 841 IF (l == lin) THEN 842 jf = jf + 1 843 IF (toadd) THEN 844 j1 = jp + 1 845 j2 = jp + npgf(iset) 846 zval(j1:j2) = zet(1:npgf(iset), iset) 847 jp = jp + npgf(iset) 848 toadd = .FALSE. 849 END IF 850 gcval(j1:j2, jf) = gcc(1:npgf(iset), ishell, iset) 851 END IF 852 END DO 853 END DO 854 855 END SUBROUTINE get_basis_functions 856 857END MODULE auto_basis 858