1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Working with the DFTB parameter types. 8!> \author JGH (24.02.2007) 9! ************************************************************************************************** 10MODULE qs_dftb_utils 11 12 USE cp_log_handling, ONLY: cp_get_default_logger,& 13 cp_logger_type 14 USE cp_output_handling, ONLY: cp_p_file,& 15 cp_print_key_finished_output,& 16 cp_print_key_should_output,& 17 cp_print_key_unit_nr 18 USE input_section_types, ONLY: section_vals_type 19 USE kinds, ONLY: default_string_length,& 20 dp 21 USE qs_dftb_types, ONLY: qs_dftb_atom_type 22#include "./base/base_uses.f90" 23 24 IMPLICIT NONE 25 26 PRIVATE 27 28 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_utils' 29 30 ! Maximum number of points used for interpolation 31 INTEGER, PARAMETER :: max_inter = 5 32 ! Maximum number of points used for extrapolation 33 INTEGER, PARAMETER :: max_extra = 9 34 ! see also qs_dftb_parameters 35 REAL(dp), PARAMETER :: slako_d0 = 1._dp 36 ! pointer to skab 37 INTEGER, DIMENSION(0:3, 0:3, 0:3, 0:3, 0:3):: iptr 38 ! small real number 39 REAL(dp), PARAMETER :: rtiny = 1.e-10_dp 40 ! eta(0) for mm atoms and non-scc qm atoms 41 REAL(dp), PARAMETER :: eta_mm = 0.47_dp 42 ! step size for qmmm finite difference 43 REAL(dp), PARAMETER :: ddrmm = 0.0001_dp 44 45 PUBLIC :: allocate_dftb_atom_param, & 46 deallocate_dftb_atom_param, & 47 get_dftb_atom_param, & 48 set_dftb_atom_param, & 49 write_dftb_atom_param 50 PUBLIC :: compute_block_sk, & 51 urep_egr, iptr 52 53CONTAINS 54 55! ************************************************************************************************** 56!> \brief ... 57!> \param dftb_parameter ... 58! ************************************************************************************************** 59 SUBROUTINE allocate_dftb_atom_param(dftb_parameter) 60 61 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 62 63 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_dftb_atom_param', & 64 routineP = moduleN//':'//routineN 65 66 IF (ASSOCIATED(dftb_parameter)) & 67 CALL deallocate_dftb_atom_param(dftb_parameter) 68 69 ALLOCATE (dftb_parameter) 70 71 dftb_parameter%defined = .FALSE. 72 dftb_parameter%name = "" 73 dftb_parameter%typ = "NONE" 74 dftb_parameter%z = -1 75 dftb_parameter%zeff = -1.0_dp 76 dftb_parameter%natorb = 0 77 dftb_parameter%lmax = -1 78 dftb_parameter%skself = 0.0_dp 79 dftb_parameter%occupation = 0.0_dp 80 dftb_parameter%eta = 0.0_dp 81 dftb_parameter%energy = 0.0_dp 82 dftb_parameter%xi = 0.0_dp 83 dftb_parameter%di = 0.0_dp 84 dftb_parameter%rcdisp = 0.0_dp 85 dftb_parameter%dudq = 0.0_dp 86 87 END SUBROUTINE allocate_dftb_atom_param 88 89! ************************************************************************************************** 90!> \brief ... 91!> \param dftb_parameter ... 92! ************************************************************************************************** 93 SUBROUTINE deallocate_dftb_atom_param(dftb_parameter) 94 95 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 96 97 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_dftb_atom_param', & 98 routineP = moduleN//':'//routineN 99 100 CPASSERT(ASSOCIATED(dftb_parameter)) 101 DEALLOCATE (dftb_parameter) 102 103 END SUBROUTINE deallocate_dftb_atom_param 104 105! ************************************************************************************************** 106!> \brief ... 107!> \param dftb_parameter ... 108!> \param name ... 109!> \param typ ... 110!> \param defined ... 111!> \param z ... 112!> \param zeff ... 113!> \param natorb ... 114!> \param lmax ... 115!> \param skself ... 116!> \param occupation ... 117!> \param eta ... 118!> \param energy ... 119!> \param cutoff ... 120!> \param xi ... 121!> \param di ... 122!> \param rcdisp ... 123!> \param dudq ... 124! ************************************************************************************************** 125 SUBROUTINE get_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, & 126 lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq) 127 128 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 129 CHARACTER(LEN=default_string_length), & 130 INTENT(OUT), OPTIONAL :: name, typ 131 LOGICAL, INTENT(OUT), OPTIONAL :: defined 132 INTEGER, INTENT(OUT), OPTIONAL :: z 133 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: zeff 134 INTEGER, INTENT(OUT), OPTIONAL :: natorb, lmax 135 REAL(KIND=dp), DIMENSION(0:3), OPTIONAL :: skself, occupation, eta 136 REAL(KIND=dp), OPTIONAL :: energy, cutoff, xi, di, rcdisp, dudq 137 138 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_dftb_atom_param', & 139 routineP = moduleN//':'//routineN 140 141 CPASSERT(ASSOCIATED(dftb_parameter)) 142 143 IF (PRESENT(name)) name = dftb_parameter%name 144 IF (PRESENT(typ)) typ = dftb_parameter%typ 145 IF (PRESENT(defined)) defined = dftb_parameter%defined 146 IF (PRESENT(z)) z = dftb_parameter%z 147 IF (PRESENT(zeff)) zeff = dftb_parameter%zeff 148 IF (PRESENT(natorb)) natorb = dftb_parameter%natorb 149 IF (PRESENT(lmax)) lmax = dftb_parameter%lmax 150 IF (PRESENT(skself)) skself = dftb_parameter%skself 151 IF (PRESENT(eta)) eta = dftb_parameter%eta 152 IF (PRESENT(energy)) energy = dftb_parameter%energy 153 IF (PRESENT(cutoff)) cutoff = dftb_parameter%cutoff 154 IF (PRESENT(occupation)) occupation = dftb_parameter%occupation 155 IF (PRESENT(xi)) xi = dftb_parameter%xi 156 IF (PRESENT(di)) di = dftb_parameter%di 157 IF (PRESENT(rcdisp)) rcdisp = dftb_parameter%rcdisp 158 IF (PRESENT(dudq)) dudq = dftb_parameter%dudq 159 160 END SUBROUTINE get_dftb_atom_param 161 162! ************************************************************************************************** 163!> \brief ... 164!> \param dftb_parameter ... 165!> \param name ... 166!> \param typ ... 167!> \param defined ... 168!> \param z ... 169!> \param zeff ... 170!> \param natorb ... 171!> \param lmax ... 172!> \param skself ... 173!> \param occupation ... 174!> \param eta ... 175!> \param energy ... 176!> \param cutoff ... 177!> \param xi ... 178!> \param di ... 179!> \param rcdisp ... 180!> \param dudq ... 181! ************************************************************************************************** 182 SUBROUTINE set_dftb_atom_param(dftb_parameter, name, typ, defined, z, zeff, natorb, & 183 lmax, skself, occupation, eta, energy, cutoff, xi, di, rcdisp, dudq) 184 185 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 186 CHARACTER(LEN=default_string_length), INTENT(IN), & 187 OPTIONAL :: name, typ 188 LOGICAL, INTENT(IN), OPTIONAL :: defined 189 INTEGER, INTENT(IN), OPTIONAL :: z 190 REAL(KIND=dp), INTENT(IN), OPTIONAL :: zeff 191 INTEGER, INTENT(IN), OPTIONAL :: natorb, lmax 192 REAL(KIND=dp), DIMENSION(0:3), OPTIONAL :: skself, occupation, eta 193 REAL(KIND=dp), OPTIONAL :: energy, cutoff, xi, di, rcdisp, dudq 194 195 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_dftb_atom_param', & 196 routineP = moduleN//':'//routineN 197 198 CPASSERT(ASSOCIATED(dftb_parameter)) 199 200 IF (PRESENT(name)) dftb_parameter%name = name 201 IF (PRESENT(typ)) dftb_parameter%typ = typ 202 IF (PRESENT(defined)) dftb_parameter%defined = defined 203 IF (PRESENT(z)) dftb_parameter%z = z 204 IF (PRESENT(zeff)) dftb_parameter%zeff = zeff 205 IF (PRESENT(natorb)) dftb_parameter%natorb = natorb 206 IF (PRESENT(lmax)) dftb_parameter%lmax = lmax 207 IF (PRESENT(skself)) dftb_parameter%skself = skself 208 IF (PRESENT(eta)) dftb_parameter%eta = eta 209 IF (PRESENT(occupation)) dftb_parameter%occupation = occupation 210 IF (PRESENT(energy)) dftb_parameter%energy = energy 211 IF (PRESENT(cutoff)) dftb_parameter%cutoff = cutoff 212 IF (PRESENT(xi)) dftb_parameter%xi = xi 213 IF (PRESENT(di)) dftb_parameter%di = di 214 IF (PRESENT(rcdisp)) dftb_parameter%rcdisp = rcdisp 215 IF (PRESENT(dudq)) dftb_parameter%dudq = dudq 216 217 END SUBROUTINE set_dftb_atom_param 218 219! ************************************************************************************************** 220!> \brief ... 221!> \param dftb_parameter ... 222!> \param subsys_section ... 223! ************************************************************************************************** 224 SUBROUTINE write_dftb_atom_param(dftb_parameter, subsys_section) 225 226 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 227 TYPE(section_vals_type), POINTER :: subsys_section 228 229 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_dftb_atom_param', & 230 routineP = moduleN//':'//routineN 231 232 CHARACTER(LEN=default_string_length) :: name, typ 233 INTEGER :: lmax, natorb, output_unit, z 234 LOGICAL :: defined 235 REAL(dp) :: zeff 236 TYPE(cp_logger_type), POINTER :: logger 237 238 NULLIFY (logger) 239 logger => cp_get_default_logger() 240 IF (ASSOCIATED(dftb_parameter) .AND. & 241 BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, & 242 "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN 243 244 output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", & 245 extension=".Log") 246 247 IF (output_unit > 0) THEN 248 CALL get_dftb_atom_param(dftb_parameter, name=name, typ=typ, defined=defined, & 249 z=z, zeff=zeff, natorb=natorb, lmax=lmax) 250 251 WRITE (UNIT=output_unit, FMT="(/,A,T67,A14)") & 252 " DFTB parameters: ", TRIM(name) 253 IF (defined) THEN 254 WRITE (UNIT=output_unit, FMT="(T16,A,T71,F10.2)") & 255 "Effective core charge:", zeff 256 WRITE (UNIT=output_unit, FMT="(T16,A,T71,I10)") & 257 "Number of orbitals:", natorb 258 ELSE 259 WRITE (UNIT=output_unit, FMT="(T55,A)") & 260 "Parameters are not defined" 261 END IF 262 END IF 263 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 264 "PRINT%KINDS") 265 END IF 266 267 END SUBROUTINE write_dftb_atom_param 268 269! ************************************************************************************************** 270!> \brief ... 271!> \param block ... 272!> \param smatij ... 273!> \param smatji ... 274!> \param rij ... 275!> \param ngrd ... 276!> \param ngrdcut ... 277!> \param dgrd ... 278!> \param llm ... 279!> \param lmaxi ... 280!> \param lmaxj ... 281!> \param irow ... 282!> \param iatom ... 283! ************************************************************************************************** 284 SUBROUTINE compute_block_sk(block, smatij, smatji, rij, ngrd, ngrdcut, dgrd, & 285 llm, lmaxi, lmaxj, irow, iatom) 286 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, smatij, smatji 287 REAL(KIND=dp), DIMENSION(3) :: rij 288 INTEGER :: ngrd, ngrdcut 289 REAL(KIND=dp) :: dgrd 290 INTEGER :: llm, lmaxi, lmaxj, irow, iatom 291 292 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_block_sk', & 293 routineP = moduleN//':'//routineN 294 295 REAL(KIND=dp) :: dr 296 REAL(KIND=dp), DIMENSION(20) :: skabij, skabji 297 298 dr = SQRT(SUM(rij(:)**2)) 299 CALL getskz(smatij, skabij, dr, ngrd, ngrdcut, dgrd, llm) 300 CALL getskz(smatji, skabji, dr, ngrd, ngrdcut, dgrd, llm) 301 IF (irow == iatom) THEN 302 CALL turnsk(block, skabji, skabij, rij, dr, lmaxi, lmaxj) 303 ELSE 304 CALL turnsk(block, skabij, skabji, -rij, dr, lmaxj, lmaxi) 305 END IF 306 307 END SUBROUTINE compute_block_sk 308 309! ************************************************************************************************** 310!> \brief Gets matrix elements on z axis, as they are stored in the tables 311!> \param slakotab ... 312!> \param skpar ... 313!> \param dx ... 314!> \param ngrd ... 315!> \param ngrdcut ... 316!> \param dgrd ... 317!> \param llm ... 318!> \author 07. Feb. 2004, TH 319! ************************************************************************************************** 320 SUBROUTINE getskz(slakotab, skpar, dx, ngrd, ngrdcut, dgrd, llm) 321 REAL(dp), INTENT(in) :: slakotab(:, :), dx 322 INTEGER, INTENT(in) :: ngrd, ngrdcut 323 REAL(dp), INTENT(in) :: dgrd 324 INTEGER, INTENT(in) :: llm 325 REAL(dp), INTENT(out) :: skpar(llm) 326 327 INTEGER :: clgp 328 329 skpar = 0._dp 330 ! 331 ! Determine closest grid point 332 ! 333 clgp = NINT(dx/dgrd) 334 ! 335 ! Screen elements which are too far away 336 ! 337 IF (clgp > ngrdcut) RETURN 338 ! 339 ! The grid point is either contained in the table --> matrix element 340 ! can be interpolated, or it is outside the table --> matrix element 341 ! needs to be extrapolated. 342 ! 343 IF (clgp > ngrd) THEN 344 ! 345 ! Extrapolate external matrix elements if table does not finish with zero 346 ! 347 CALL extrapol(slakotab, skpar, dx, ngrd, dgrd, llm) 348 ELSE 349 ! 350 ! Interpolate tabulated matrix elements 351 ! 352 CALL interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp) 353 END IF 354 END SUBROUTINE getskz 355 356! ************************************************************************************************** 357!> \brief ... 358!> \param slakotab ... 359!> \param skpar ... 360!> \param dx ... 361!> \param ngrd ... 362!> \param dgrd ... 363!> \param llm ... 364!> \param clgp ... 365! ************************************************************************************************** 366 SUBROUTINE interpol(slakotab, skpar, dx, ngrd, dgrd, llm, clgp) 367 REAL(dp), INTENT(in) :: slakotab(:, :), dx 368 INTEGER, INTENT(in) :: ngrd 369 REAL(dp), INTENT(in) :: dgrd 370 INTEGER, INTENT(in) :: llm 371 REAL(dp), INTENT(out) :: skpar(llm) 372 INTEGER, INTENT(in) :: clgp 373 374 INTEGER :: fgpm, k, l, lgpm 375 REAL(dp) :: error, xa(max_inter), ya(max_inter) 376 377 lgpm = MIN(clgp + INT(max_inter/2.0), ngrd) 378 fgpm = lgpm - max_inter + 1 379 DO k = 0, max_inter - 1 380 xa(k + 1) = (fgpm + k)*dgrd 381 END DO 382 ! 383 ! Interpolate matrix elements for all orbitals 384 ! 385 DO l = 1, llm 386 ! 387 ! Read SK parameters from table 388 ! 389 ya(1:max_inter) = slakotab(fgpm:lgpm, l) 390 CALL polint(xa, ya, max_inter, dx, skpar(l), error) 391 END DO 392 END SUBROUTINE interpol 393 394! ************************************************************************************************** 395!> \brief ... 396!> \param slakotab ... 397!> \param skpar ... 398!> \param dx ... 399!> \param ngrd ... 400!> \param dgrd ... 401!> \param llm ... 402! ************************************************************************************************** 403 SUBROUTINE extrapol(slakotab, skpar, dx, ngrd, dgrd, llm) 404 REAL(dp), INTENT(in) :: slakotab(:, :), dx 405 INTEGER, INTENT(in) :: ngrd 406 REAL(dp), INTENT(in) :: dgrd 407 INTEGER, INTENT(in) :: llm 408 REAL(dp), INTENT(out) :: skpar(llm) 409 410 INTEGER :: fgp, k, l, lgp, ntable, nzero 411 REAL(dp) :: error, xa(max_extra), ya(max_extra) 412 413 nzero = max_extra/3 414 ntable = max_extra - nzero 415 ! 416 ! Get the three last distances from the table 417 ! 418 DO k = 1, ntable 419 xa(k) = (ngrd - (max_extra - 3) + k)*dgrd 420 END DO 421 DO k = 1, nzero 422 xa(ntable + k) = (ngrd + k - 1)*dgrd + slako_d0 423 ya(ntable + k) = 0.0 424 END DO 425 ! 426 ! Extrapolate matrix elements for all orbitals 427 ! 428 DO l = 1, llm 429 ! 430 ! Read SK parameters from table 431 ! 432 fgp = ngrd + 1 - (max_extra - 3) 433 lgp = ngrd 434 ya(1:max_extra - 3) = slakotab(fgp:lgp, l) 435 CALL polint(xa, ya, max_extra, dx, skpar(l), error) 436 END DO 437 END SUBROUTINE extrapol 438 439! ************************************************************************************************** 440!> \brief Turn matrix element from z-axis to orientation of dxv 441!> \param mat ... 442!> \param skab1 ... 443!> \param skab2 ... 444!> \param dxv ... 445!> \param dx ... 446!> \param lmaxa ... 447!> \param lmaxb ... 448!> \date 13. Jan 2004 449!> \par Notes 450!> These routines are taken from an old TB code (unknown to TH). 451!> They are highly optimised and taken because they are time critical. 452!> They are explicit, so not recursive, and work up to d functions. 453!> 454!> Set variables necessary for rotation of matrix elements 455!> 456!> r_i^2/r, replicated in rr2(4:6) for index convenience later 457!> r_i/r, direction vector, rr(4:6) are replicated from 1:3 458!> lmax of A and B 459!> \author TH 460!> \version 1.0 461! ************************************************************************************************** 462 SUBROUTINE turnsk(mat, skab1, skab2, dxv, dx, lmaxa, lmaxb) 463 REAL(dp), INTENT(inout) :: mat(:, :) 464 REAL(dp), INTENT(in) :: skab1(:), skab2(:), dxv(3), dx 465 INTEGER, INTENT(in) :: lmaxa, lmaxb 466 467 CHARACTER(len=*), PARAMETER :: routineN = 'turnsk', & 468 routineP = moduleN//':'//routineN 469 470 INTEGER :: lmaxab, minlmaxab 471 REAL(dp) :: rinv, rr(6), rr2(6) 472 473 lmaxab = MAX(lmaxa, lmaxb) 474 ! Determine l quantum limits. 475 IF (lmaxab .GT. 2) CPABORT('lmax=2') 476 minlmaxab = MIN(lmaxa, lmaxb) 477 ! 478 ! s-s interaction 479 ! 480 CALL skss(skab1, mat) 481 ! 482 IF (lmaxab .LE. 0) RETURN 483 ! 484 rr2(1:3) = dxv(1:3)**2 485 rr(1:3) = dxv(1:3) 486 rinv = 1.0_dp/dx 487 ! 488 rr(1:3) = rr(1:3)*rinv 489 rr(4:6) = rr(1:3) 490 rr2(1:3) = rr2(1:3)*rinv**2 491 rr2(4:6) = rr2(1:3) 492 ! 493 ! s-p, p-s and p-p interaction 494 ! 495 IF (minlmaxab .GE. 1) THEN 496 CALL skpp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb)) 497 CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 498 CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 499 ELSE 500 IF (lmaxb .GE. 1) THEN 501 CALL sksp(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 502 ELSE 503 CALL sksp(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 504 END IF 505 END IF 506 ! 507 ! If there is only s-p interaction we have finished 508 ! 509 IF (lmaxab .LE. 1) RETURN 510 ! 511 ! at least one atom has d functions 512 ! 513 IF (minlmaxab .EQ. 2) THEN 514 ! 515 ! in case both atoms have d functions 516 ! 517 CALL skdd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb)) 518 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 519 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 520 CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 521 CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 522 ELSE 523 ! 524 ! One atom has d functions, the other has s or s and p functions 525 ! 526 IF (lmaxa .EQ. 0) THEN 527 ! 528 ! atom b has d, the atom a only s functions 529 ! 530 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 531 ELSE IF (lmaxa .EQ. 1) THEN 532 ! 533 ! atom b has d, the atom a s and p functions 534 ! 535 CALL sksd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 536 CALL skpd(skab2, mat, iptr(:, :, :, lmaxa, lmaxb), .TRUE.) 537 ELSE 538 ! 539 ! atom a has d functions 540 ! 541 IF (lmaxb .EQ. 0) THEN 542 ! 543 ! atom a has d, atom b has only s functions 544 ! 545 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 546 ELSE 547 ! 548 ! atom a has d, atom b has s and p functions 549 ! 550 CALL sksd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 551 CALL skpd(skab1, mat, iptr(:, :, :, lmaxa, lmaxb), .FALSE.) 552 END IF 553 END IF 554 END IF 555 ! 556 CONTAINS 557 ! 558 ! The subroutines to turn the matrix elements are taken as internal subroutines 559 ! as it is beneficial to inline them. 560 ! 561 ! They are both turning the matrix elements and placing them appropriately 562 ! into the matrix block 563 ! 564! ************************************************************************************************** 565!> \brief s-s interaction (no rotation necessary) 566!> \param skpar ... 567!> \param mat ... 568!> \version 1.0 569! ************************************************************************************************** 570 SUBROUTINE skss(skpar, mat) 571 REAL(dp), INTENT(in) :: skpar(:) 572 REAL(dp), INTENT(inout) :: mat(:, :) 573 574 mat(1, 1) = mat(1, 1) + skpar(1) 575 ! 576 END SUBROUTINE skss 577 578! ************************************************************************************************** 579!> \brief s-p interaction (simple rotation) 580!> \param skpar ... 581!> \param mat ... 582!> \param ind ... 583!> \param transposed ... 584!> \version 1.0 585! ************************************************************************************************** 586 SUBROUTINE sksp(skpar, mat, ind, transposed) 587 REAL(dp), INTENT(in) :: skpar(:) 588 REAL(dp), INTENT(inout) :: mat(:, :) 589 INTEGER, INTENT(in) :: ind(0:, 0:, 0:) 590 LOGICAL, INTENT(in) :: transposed 591 592 INTEGER :: l 593 REAL(dp) :: skp 594 595 skp = skpar(ind(1, 0, 0)) 596 IF (transposed) THEN 597 DO l = 1, 3 598 mat(1, l + 1) = mat(1, l + 1) + rr(l)*skp 599 END DO 600 ELSE 601 DO l = 1, 3 602 mat(l + 1, 1) = mat(l + 1, 1) - rr(l)*skp 603 END DO 604 END IF 605 ! 606 END SUBROUTINE sksp 607 608! ************************************************************************************************** 609!> \brief ... 610!> \param skpar ... 611!> \param mat ... 612!> \param ind ... 613! ************************************************************************************************** 614 SUBROUTINE skpp(skpar, mat, ind) 615 REAL(dp), INTENT(in) :: skpar(:) 616 REAL(dp), INTENT(inout) :: mat(:, :) 617 INTEGER, INTENT(in) :: ind(0:, 0:, 0:) 618 619 INTEGER :: ii, ir, is, k, l 620 REAL(dp) :: epp(6), matel(6), skppp, skpps 621 622 epp(1:3) = rr2(1:3) 623 DO l = 1, 3 624 epp(l + 3) = rr(l)*rr(l + 1) 625 END DO 626 skppp = skpar(ind(1, 1, 1)) 627 skpps = skpar(ind(1, 1, 0)) 628 ! 629 DO l = 1, 3 630 matel(l) = epp(l)*skpps + (1._dp - epp(l))*skppp 631 END DO 632 DO l = 4, 6 633 matel(l) = epp(l)*(skpps - skppp) 634 END DO 635 ! 636 DO ir = 1, 3 637 DO is = 1, ir - 1 638 ii = ir - is 639 k = 3*ii - (ii*(ii - 1))/2 + is 640 mat(is + 1, ir + 1) = mat(is + 1, ir + 1) + matel(k) 641 mat(ir + 1, is + 1) = mat(ir + 1, is + 1) + matel(k) 642 END DO 643 mat(ir + 1, ir + 1) = mat(ir + 1, ir + 1) + matel(ir) 644 END DO 645 END SUBROUTINE skpp 646 647! ************************************************************************************************** 648!> \brief ... 649!> \param skpar ... 650!> \param mat ... 651!> \param ind ... 652!> \param transposed ... 653! ************************************************************************************************** 654 SUBROUTINE sksd(skpar, mat, ind, transposed) 655 REAL(dp), INTENT(in) :: skpar(:) 656 REAL(dp), INTENT(inout) :: mat(:, :) 657 INTEGER, INTENT(in) :: ind(0:, 0:, 0:) 658 LOGICAL, INTENT(in) :: transposed 659 660 INTEGER :: l 661 REAL(dp) :: d4, d5, es(5), r3, sksds 662 663 sksds = skpar(ind(2, 0, 0)) 664 r3 = SQRT(3._dp) 665 d4 = rr2(3) - 0.5_dp*(rr2(1) + rr2(2)) 666 d5 = rr2(1) - rr2(2) 667 ! 668 DO l = 1, 3 669 es(l) = r3*rr(l)*rr(l + 1) 670 END DO 671 es(4) = 0.5_dp*r3*d5 672 es(5) = d4 673 ! 674 IF (transposed) THEN 675 DO l = 1, 5 676 mat(1, l + 4) = mat(1, l + 4) + es(l)*sksds 677 END DO 678 ELSE 679 DO l = 1, 5 680 mat(l + 4, 1) = mat(l + 4, 1) + es(l)*sksds 681 END DO 682 END IF 683 END SUBROUTINE sksd 684 685! ************************************************************************************************** 686!> \brief ... 687!> \param skpar ... 688!> \param mat ... 689!> \param ind ... 690!> \param transposed ... 691! ************************************************************************************************** 692 SUBROUTINE skpd(skpar, mat, ind, transposed) 693 REAL(dp), INTENT(in) :: skpar(:) 694 REAL(dp), INTENT(inout) :: mat(:, :) 695 INTEGER, INTENT(in) :: ind(0:, 0:, 0:) 696 LOGICAL, INTENT(in) :: transposed 697 698 INTEGER :: ir, is, k, l, m 699 REAL(dp) :: d3, d4, d5, d6, dm(15), epd(13, 2), r3, & 700 sktmp 701 702 r3 = SQRT(3.0_dp) 703 d3 = rr2(1) + rr2(2) 704 d4 = rr2(3) - 0.5_dp*d3 705 d5 = rr2(1) - rr2(2) 706 d6 = rr(1)*rr(2)*rr(3) 707 DO l = 1, 3 708 epd(l, 1) = r3*rr2(l)*rr(l + 1) 709 epd(l, 2) = rr(l + 1)*(1.0_dp - 2._dp*rr2(l)) 710 epd(l + 4, 1) = r3*rr2(l)*rr(l + 2) 711 epd(l + 4, 2) = rr(l + 2)*(1.0_dp - 2*rr2(l)) 712 epd(l + 7, 1) = 0.5_dp*r3*rr(l)*d5 713 epd(l + 10, 1) = rr(l)*d4 714 END DO 715 ! 716 epd(4, 1) = r3*d6 717 epd(4, 2) = -2._dp*d6 718 epd(8, 2) = rr(1)*(1.0_dp - d5) 719 epd(9, 2) = -rr(2)*(1.0_dp + d5) 720 epd(10, 2) = -rr(3)*d5 721 epd(11, 2) = -r3*rr(1)*rr2(3) 722 epd(12, 2) = -r3*rr(2)*rr2(3) 723 epd(13, 2) = r3*rr(3)*d3 724 ! 725 dm(1:15) = 0.0_dp 726 ! 727 DO m = 1, 2 728 sktmp = skpar(ind(2, 1, m - 1)) 729 dm(1) = dm(1) + epd(1, m)*sktmp 730 dm(2) = dm(2) + epd(6, m)*sktmp 731 dm(3) = dm(3) + epd(4, m)*sktmp 732 dm(5) = dm(5) + epd(2, m)*sktmp 733 dm(6) = dm(6) + epd(7, m)*sktmp 734 dm(7) = dm(7) + epd(5, m)*sktmp 735 dm(9) = dm(9) + epd(3, m)*sktmp 736 DO l = 8, 13 737 dm(l + 2) = dm(l + 2) + epd(l, m)*sktmp 738 END DO 739 END DO 740 ! 741 dm(4) = dm(3) 742 dm(8) = dm(3) 743 ! 744 IF (transposed) THEN 745 DO ir = 1, 5 746 DO is = 1, 3 747 k = 3*(ir - 1) + is 748 mat(is + 1, ir + 4) = mat(is + 1, ir + 4) + dm(k) 749 END DO 750 END DO 751 ELSE 752 DO ir = 1, 5 753 DO is = 1, 3 754 k = 3*(ir - 1) + is 755 mat(ir + 4, is + 1) = mat(ir + 4, is + 1) - dm(k) 756 END DO 757 END DO 758 END IF 759 ! 760 END SUBROUTINE skpd 761 762! ************************************************************************************************** 763!> \brief ... 764!> \param skpar ... 765!> \param mat ... 766!> \param ind ... 767! ************************************************************************************************** 768 SUBROUTINE skdd(skpar, mat, ind) 769 REAL(dp), INTENT(in) :: skpar(:) 770 REAL(dp), INTENT(inout) :: mat(:, :) 771 INTEGER, INTENT(in) :: ind(0:, 0:, 0:) 772 773 INTEGER :: ii, ir, is, k, l, m 774 REAL(dp) :: d3, d4, d5, dd(3), dm(15), e(15, 3), r3 775 776 r3 = SQRT(3._dp) 777 d3 = rr2(1) + rr2(2) 778 d4 = rr2(3) - 0.5_dp*d3 779 d5 = rr2(1) - rr2(2) 780 DO l = 1, 3 781 e(l, 1) = rr2(l)*rr2(l + 1) 782 e(l, 2) = rr2(l) + rr2(l + 1) - 4._dp*e(l, 1) 783 e(l, 3) = rr2(l + 2) + e(l, 1) 784 e(l, 1) = 3._dp*e(l, 1) 785 END DO 786 e(4, 1) = d5**2 787 e(4, 2) = d3 - e(4, 1) 788 e(4, 3) = rr2(3) + 0.25_dp*e(4, 1) 789 e(4, 1) = 0.75_dp*e(4, 1) 790 e(5, 1) = d4**2 791 e(5, 2) = 3._dp*rr2(3)*d3 792 e(5, 3) = 0.75_dp*d3**2 793 dd(1) = rr(1)*rr(3) 794 dd(2) = rr(2)*rr(1) 795 dd(3) = rr(3)*rr(2) 796 DO l = 1, 2 797 e(l + 5, 1) = 3._dp*rr2(l + 1)*dd(l) 798 e(l + 5, 2) = dd(l)*(1._dp - 4._dp*rr2(l + 1)) 799 e(l + 5, 3) = dd(l)*(rr2(l + 1) - 1._dp) 800 END DO 801 e(8, 1) = dd(1)*d5*1.5_dp 802 e(8, 2) = dd(1)*(1.0_dp - 2.0_dp*d5) 803 e(8, 3) = dd(1)*(0.5_dp*d5 - 1.0_dp) 804 e(9, 1) = d5*0.5_dp*d4*r3 805 e(9, 2) = -d5*rr2(3)*r3 806 e(9, 3) = d5*0.25_dp*(1.0_dp + rr2(3))*r3 807 e(10, 1) = rr2(1)*dd(3)*3.0_dp 808 e(10, 2) = (0.25_dp - rr2(1))*dd(3)*4.0_dp 809 e(10, 3) = dd(3)*(rr2(1) - 1.0_dp) 810 e(11, 1) = 1.5_dp*dd(3)*d5 811 e(11, 2) = -dd(3)*(1.0_dp + 2.0_dp*d5) 812 e(11, 3) = dd(3)*(1.0_dp + 0.5_dp*d5) 813 e(13, 3) = 0.5_dp*d5*dd(2) 814 e(13, 2) = -2.0_dp*dd(2)*d5 815 e(13, 1) = e(13, 3)*3.0_dp 816 e(12, 1) = d4*dd(1)*r3 817 e(14, 1) = d4*dd(3)*r3 818 e(15, 1) = d4*dd(2)*r3 819 e(15, 2) = -2.0_dp*r3*dd(2)*rr2(3) 820 e(15, 3) = 0.5_dp*r3*(1.0_dp + rr2(3))*dd(2) 821 e(14, 2) = r3*dd(3)*(d3 - rr2(3)) 822 e(14, 3) = -r3*0.5_dp*dd(3)*d3 823 e(12, 2) = r3*dd(1)*(d3 - rr2(3)) 824 e(12, 3) = -r3*0.5_dp*dd(1)*d3 825 ! 826 dm(1:15) = 0._dp 827 DO l = 1, 15 828 DO m = 1, 3 829 dm(l) = dm(l) + e(l, m)*skpar(ind(2, 2, m - 1)) 830 END DO 831 END DO 832 ! 833 DO ir = 1, 5 834 DO is = 1, ir - 1 835 ii = ir - is 836 k = 5*ii - (ii*(ii - 1))/2 + is 837 mat(ir + 4, is + 4) = mat(ir + 4, is + 4) + dm(k) 838 mat(is + 4, ir + 4) = mat(is + 4, ir + 4) + dm(k) 839 END DO 840 mat(ir + 4, ir + 4) = mat(ir + 4, ir + 4) + dm(ir) 841 END DO 842 END SUBROUTINE skdd 843 ! 844 END SUBROUTINE turnsk 845 846! ************************************************************************************************** 847!> \brief ... 848!> \param xa ... 849!> \param ya ... 850!> \param n ... 851!> \param x ... 852!> \param y ... 853!> \param dy ... 854! ************************************************************************************************** 855 SUBROUTINE polint(xa, ya, n, x, y, dy) 856 INTEGER, INTENT(in) :: n 857 REAL(dp), INTENT(in) :: ya(n), xa(n), x 858 REAL(dp), INTENT(out) :: y, dy 859 860 CHARACTER(len=*), PARAMETER :: routineN = 'polint', routineP = moduleN//':'//routineN 861 862 INTEGER :: i, m, ns 863 REAL(dp) :: c(n), d(n), den, dif, dift, ho, hp, w 864 865! 866! 867 868 ns = 1 869 870 dif = ABS(x - xa(1)) 871 DO i = 1, n 872 dift = ABS(x - xa(i)) 873 IF (dift .LT. dif) THEN 874 ns = i 875 dif = dift 876 ENDIF 877 c(i) = ya(i) 878 d(i) = ya(i) 879 END DO 880 ! 881 y = ya(ns) 882 ns = ns - 1 883 DO m = 1, n - 1 884 DO i = 1, n - m 885 ho = xa(i) - x 886 hp = xa(i + m) - x 887 w = c(i + 1) - d(i) 888 den = ho - hp 889 CPASSERT(den /= 0.0_dp) 890 den = w/den 891 d(i) = hp*den 892 c(i) = ho*den 893 END DO 894 IF (2*ns .LT. n - m) THEN 895 dy = c(ns + 1) 896 ELSE 897 dy = d(ns) 898 ns = ns - 1 899 ENDIF 900 y = y + dy 901 END DO 902! 903 RETURN 904 END SUBROUTINE polint 905 906! ************************************************************************************************** 907!> \brief ... 908!> \param rv ... 909!> \param r ... 910!> \param erep ... 911!> \param derep ... 912!> \param n_urpoly ... 913!> \param urep ... 914!> \param spdim ... 915!> \param s_cut ... 916!> \param srep ... 917!> \param spxr ... 918!> \param scoeff ... 919!> \param surr ... 920!> \param dograd ... 921! ************************************************************************************************** 922 SUBROUTINE urep_egr(rv, r, erep, derep, & 923 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, dograd) 924 925 REAL(dp), INTENT(in) :: rv(3), r 926 REAL(dp), INTENT(inout) :: erep, derep(3) 927 INTEGER, INTENT(in) :: n_urpoly 928 REAL(dp), INTENT(in) :: urep(:) 929 INTEGER, INTENT(in) :: spdim 930 REAL(dp), INTENT(in) :: s_cut, srep(3) 931 REAL(dp), POINTER :: spxr(:, :), scoeff(:, :) 932 REAL(dp), INTENT(in) :: surr(2) 933 LOGICAL, INTENT(in) :: dograd 934 935 INTEGER :: ic, isp, jsp, nsp 936 REAL(dp) :: de_z, rz 937 938 derep = 0._dp 939 de_z = 0._dp 940 IF (n_urpoly > 0) THEN 941 ! 942 ! polynomial part 943 ! 944 rz = urep(1) - r 945 IF (rz <= rtiny) RETURN 946 DO ic = 2, n_urpoly 947 erep = erep + urep(ic)*rz**(ic) 948 END DO 949 IF (dograd) THEN 950 DO ic = 2, n_urpoly 951 de_z = de_z - ic*urep(ic)*rz**(ic - 1) 952 END DO 953 END IF 954 ELSE IF (spdim > 0) THEN 955 ! 956 ! spline part 957 ! 958 ! This part is kind of proprietary Paderborn code and I won't reverse-engineer 959 ! everything in detail. What is obvious is documented. 960 ! 961 ! This part has 4 regions: 962 ! a) very long range is screened 963 ! b) short-range is extrapolated with e-functions 964 ! ca) normal range is approximated with a spline 965 ! cb) longer range is extrapolated with an higher degree spline 966 ! 967 IF (r > s_cut) RETURN ! screening (condition a) 968 ! 969 IF (r < spxr(1, 1)) THEN 970 ! a) short range 971 erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3) 972 IF (dograd) de_z = de_z - srep(1)*EXP(-srep(1)*r + srep(2)) 973 ELSE 974 ! 975 ! condition c). First determine between which places the spline is located: 976 ! 977 ispg: DO isp = 1, spdim ! condition ca) 978 IF (r < spxr(isp, 1)) CYCLE ispg ! distance is smaller than this spline range 979 IF (r >= spxr(isp, 2)) CYCLE ispg ! distance is larger than this spline range 980 ! at this point we have found the correct spline interval 981 rz = r - spxr(isp, 1) 982 IF (isp /= spdim) THEN 983 nsp = 3 ! condition ca 984 DO jsp = 0, nsp 985 erep = erep + scoeff(isp, jsp + 1)*rz**(jsp) 986 END DO 987 IF (dograd) THEN 988 DO jsp = 1, nsp 989 de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1) 990 END DO 991 END IF 992 ELSE 993 nsp = 5 ! condition cb 994 DO jsp = 0, nsp 995 IF (jsp <= 3) THEN 996 erep = erep + scoeff(isp, jsp + 1)*rz**(jsp) 997 ELSE 998 erep = erep + surr(jsp - 3)*rz**(jsp) 999 ENDIF 1000 END DO 1001 IF (dograd) THEN 1002 DO jsp = 1, nsp 1003 IF (jsp <= 3) THEN 1004 de_z = de_z + jsp*scoeff(isp, jsp + 1)*rz**(jsp - 1) 1005 ELSE 1006 de_z = de_z + jsp*surr(jsp - 3)*rz**(jsp - 1) 1007 ENDIF 1008 END DO 1009 END IF 1010 END IF 1011 EXIT ispg 1012 END DO ispg 1013 END IF 1014 END IF 1015 ! 1016 IF (dograd) THEN 1017 IF (r > 1.e-12_dp) derep(1:3) = (de_z/r)*rv(1:3) 1018 END IF 1019 1020 END SUBROUTINE urep_egr 1021 1022END MODULE qs_dftb_utils 1023 1024