1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include 'common.fypp' 9 10!> Routines implementing the full 3rd order DFTB. 11module dftbp_thirdorder_module 12 use dftbp_assert 13 use dftbp_accuracy 14 use dftbp_commontypes, only : TOrbitals 15 use dftbp_shortgamma, only : expGammaCutoff 16 use dftbp_periodic, only : TNeighbourList, getNrOfNeighbours 17 use dftbp_charges 18 implicit none 19 private 20 21 public :: ThirdOrderInp, ThirdOrder, ThirdOrder_init 22 23 24 !> Input for the 3rd order module 25 type ThirdOrderInp 26 27 !> Orbital information 28 type(TOrbitals), pointer :: orb 29 30 !> Hubbard U values. Shape: [nShell, nSpecies] 31 real(dp), allocatable :: hubbUs(:,:) 32 33 !> Hubbard U derivatives. Shape: [nShell, nSpecies] 34 real(dp), allocatable :: hubbUDerivs(:,:) 35 36 !> Whether interaction should be damped when given atom is involved. 37 logical, allocatable :: damped(:) 38 39 !> Exponention of the damping 40 real(dp) :: dampExp 41 42 !> Whether third order should be considered shell resolved. If not, only the first shell of each 43 !> atom in hubbUs and hubbUDerivs is used. 44 logical :: shellResolved 45 46 end type ThirdOrderInp 47 48 49 !> Internal status of third order. 50 type ThirdOrder 51 integer :: nSpecies, nAtoms, mShells, mShellsReal 52 logical :: shellResolved 53 integer, allocatable :: nShells(:) 54 real(dp), allocatable :: UU(:,:) 55 real(dp), allocatable :: dUdQ(:,:) 56 real(dp), allocatable :: shift1(:,:), shift2(:,:), shift3(:) 57 real(dp), allocatable :: chargesPerAtom(:) 58 real(dp), allocatable :: chargesPerShell(:,:) 59 real(dp), allocatable :: cutoffs(:,:) 60 integer, allocatable :: nNeigh(:,:), nNeighMax(:) 61 real(dp), allocatable :: gamma3ab(:,:,:,:), gamma3ba(:,:,:,:) 62 logical, allocatable :: damped(:) 63 real(dp) :: dampExp 64 real(dp) :: maxCutoff 65 contains 66 procedure :: getCutoff 67 procedure :: updateCoords 68 procedure :: updateCharges 69 procedure :: getShifts 70 procedure :: getEnergyPerAtom 71 procedure :: getEnergyPerAtomXlbomd 72 procedure :: addGradientDc 73 procedure :: addGradientDcXlbomd 74 procedure :: addStressDc 75 end type ThirdOrder 76 77contains 78 79 80 !> Initializes instance. 81 subroutine ThirdOrder_init(this, inp) 82 83 !> Instance. 84 type(ThirdOrder), intent(out) :: this 85 86 !> Input data. 87 type(ThirdOrderInp), intent(in) :: inp 88 89 this%nAtoms = size(inp%orb%nOrbAtom) 90 this%mShellsReal = inp%orb%mShell 91 this%nSpecies = size(inp%hubbUs, dim=2) 92 this%shellResolved = inp%shellResolved 93 if (this%shellResolved) then 94 this%mShells = this%mShellsReal 95 this%nShells = inp%orb%nShell 96 this%UU = inp%hubbUs 97 this%dUdQ = inp%hubbUDerivs 98 else 99 this%mShells = 1 100 allocate(this%nShells(this%nSpecies)) 101 this%nShells(:) = 1 102 this%UU = inp%hubbUs(1:1, :) 103 this%dUdQ = inp%hubbUDerivs(1:1, :) 104 end if 105 106 allocate(this%cutoffs(this%nSpecies, this%nSpecies)) 107 call calcCutoffs(this%UU, this%nShells, this%cutoffs) 108 this%maxCutoff = maxval(this%cutoffs) 109 110 allocate(this%nNeigh(this%nSpecies, this%nAtoms)) 111 allocate(this%nNeighMax(this%nAtoms)) 112 allocate(this%chargesPerAtom(this%nAtoms)) 113 allocate(this%chargesPerShell(this%mShells, this%nAtoms)) 114 allocate(this%shift1(this%mShells, this%nAtoms)) 115 allocate(this%shift2(this%mShells, this%nAtoms)) 116 allocate(this%shift3(this%nAtoms)) 117 allocate(this%gamma3ab(this%mShells, this%mShells, 0, this%nAtoms)) 118 allocate(this%gamma3ba(this%mShells, this%mShells, 0, this%nAtoms)) 119 allocate(this%damped(this%nSpecies)) 120 this%damped = inp%damped 121 this%dampExp = inp%dampExp 122 123 end subroutine ThirdOrder_init 124 125 126 !> Returns real space cutoff. 127 function getCutoff(this) result(cutoff) 128 129 !> Instance 130 class(ThirdOrder), intent(inout) :: this 131 132 !> Cutoff 133 real(dp) :: cutoff 134 135 cutoff = this%maxCutoff 136 137 end function getCutoff 138 139 140 !> Updates data structures if there are changed coordinates for the instance. 141 subroutine updateCoords(this, neighList, species) 142 143 !> Instance. 144 class(ThirdOrder), intent(inout) :: this 145 146 !> Neighbour list. 147 type(TNeighbourList), intent(in) :: neighList 148 149 !> Species for all atoms, shape: [nAllAtom]. 150 integer, intent(in) :: species(:) 151 152 integer :: iNeigh, iAt1, iAt2, iSp1, iSp2, iSh1, iSh2 153 logical :: damping 154 real(dp) :: rr 155 156 this%nNeigh(:,:) = 0 157 do iAt1 = 1, this%nAtoms 158 iSp1 = species(iAt1) 159 do iSp2 = 1, this%nSpecies 160 this%nNeigh(iSp2, iAt1) = getNrOfNeighbours(neighList, this%cutoffs(iSp2, iSp1), iAt1) 161 end do 162 end do 163 this%nNeighMax = maxval(this%nNeigh, dim=1) 164 165 if (size(this%gamma3ab, dim=3) < maxval(this%nNeighMax) + 1) then 166 deallocate(this%gamma3ab) 167 deallocate(this%gamma3ba) 168 allocate(this%gamma3ab(this%mShells, this%mShells, 0:maxval(this%nNeighMax), this%nAtoms)) 169 allocate(this%gamma3ba(this%mShells, this%mShells, 0:maxval(this%nNeighMax), this%nAtoms)) 170 end if 171 this%gamma3ab(:,:,:,:) = 0.0_dp 172 this%gamma3ba(:,:,:,:) = 0.0_dp 173 do iAt1 = 1, this%nAtoms 174 iSp1 = species(iAt1) 175 do iNeigh = 0, this%nNeighMax(iAt1) 176 iAt2 = neighList%iNeighbour(iNeigh, iAt1) 177 iSp2 = species(iAt2) 178 if (iNeigh <= this%nNeigh(iSp2, iAt1)) then 179 rr = sqrt(neighList%neighDist2(iNeigh, iAt1)) 180 damping = this%damped(iSp1) .or. this%damped(iSp2) 181 do iSh1 = 1, this%nShells(iSp1) 182 do iSh2 = 1, this%nShells(iSp2) 183 this%gamma3ab(iSh2, iSh1, iNeigh, iAt1) = gamma3(this%UU(iSh1, iSp1),& 184 & this%UU(iSh2, iSp2), this%dUdQ(iSh1, iSp1), rr, damping, this%dampExp) 185 this%gamma3ba(iSh2, iSh1, iNeigh, iAt1) = gamma3(this%UU(iSh2, iSp2),& 186 & this%UU(iSh1, iSp1), this%dUdQ(iSh2, iSp2), rr, damping, this%dampExp) 187 end do 188 end do 189 end if 190 end do 191 end do 192 193 end subroutine updateCoords 194 195 196 !> Updates with changed charges for the instance. 197 subroutine updateCharges(this, species, neighList, qq, q0, img2CentCell, orb) 198 199 !> Instance 200 class(ThirdOrder), intent(inout) :: this 201 202 !> Species, shape: [nAtom] 203 integer, intent(in) :: species(:) 204 205 !> Neighbour list. 206 type(TNeighbourList), intent(in) :: neighList 207 208 !> Orbital charges. 209 real(dp), intent(in) :: qq(:,:,:) 210 211 !> Reference orbital charges. 212 real(dp), intent(in) :: q0(:,:,:) 213 214 !> Mapping on atoms in central cell. 215 integer, intent(in) :: img2CentCell(:) 216 217 !> Orbital information 218 type(TOrbitals), intent(in) :: orb 219 220 real(dp), allocatable :: chargesPerShell(:,:) 221 integer :: iAt1, iAt2f, iSp1, iSp2, iSh1, iSh2, iNeigh 222 223 @:ASSERT(size(species) == this%nAtoms) 224 @:ASSERT(size(qq, dim=2) == this%nAtoms) 225 @:ASSERT(size(q0, dim=2) == this%nAtoms) 226 227 if (this%shellResolved) then 228 call getSummedCharges(species, orb, qq, q0, dQAtom=this%chargesPerAtom,& 229 & dQShell=this%chargesPerShell) 230 else 231 ! First (only) component of this%chargesPerShell contains atomic charge 232 allocate(chargesPerShell(this%mShellsReal, this%nAtoms)) 233 call getSummedCharges(species, orb, qq, q0, dQAtom=this%chargesPerAtom,& 234 & dQShell=chargesPerShell) 235 this%chargesPerShell(1,:) = sum(chargesPerShell, dim=1) 236 end if 237 238 this%shift1(:,:) = 0.0_dp 239 this%shift2(:,:) = 0.0_dp 240 this%shift3(:) = 0.0_dp 241 242 do iAt1 = 1, this%nAtoms 243 iSp1 = species(iAt1) 244 do iNeigh = 0, this%nNeighMax(iAt1) 245 iAt2f = img2CentCell(neighList%iNeighbour(iNeigh, iAt1)) 246 iSp2 = species(iAt2f) 247 do iSh1 = 1, this%nShells(iSp1) 248 do iSh2 = 1, this%nShells(iSp2) 249 this%shift1(iSh1, iAt1) = this%shift1(iSh1, iAt1)& 250 & + this%gamma3ab(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 251 & * this%chargesPerAtom(iAt1) 252 this%shift2(iSh1, iAt1) = this%shift2(iSh1, iAt1)& 253 & + this%gamma3ba(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 254 & * this%chargesPerAtom(iAt2f) 255 this%shift3(iAt1) = this%shift3(iAt1)& 256 & + this%gamma3ab(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 257 & * this%chargesPerShell(iSh1, iAt1) 258 if (iAt2f /= iAt1) then 259 this%shift1(iSh2, iAt2f) = this%shift1(iSh2, iAt2f)& 260 & + this%gamma3ba(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh1, iAt1)& 261 & * this%chargesPerAtom(iAt2f) 262 this%shift2(iSh2, iAt2f) = this%shift2(iSh2, iAt2f)& 263 & + this%gamma3ab(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh1, iAt1)& 264 & * this%chargesPerAtom(iAt1) 265 this%shift3(iAt2f) = this%shift3(iAt2f)& 266 & + this%gamma3ba(iSh2, iSh1, iNeigh, iAt1) * this%chargesPerShell(iSh1, iAt1)& 267 & * this%chargesPerShell(iSh2, iAt2f) 268 end if 269 end do 270 end do 271 end do 272 end do 273 this%shift1(:,:) = 1.0_dp / 3.0_dp * this%shift1 274 this%shift2(:,:) = 1.0_dp / 3.0_dp * this%shift2 275 this%shift3(:) = 1.0_dp / 3.0_dp * this%shift3 276 277 end subroutine updateCharges 278 279 280 !> Returns shifts per atom. 281 subroutine getShifts(this, shiftPerAtom, shiftPerShell) 282 283 !> Instance. 284 class(ThirdOrder), intent(inout) :: this 285 286 !> Shift per atom. 287 real(dp), intent(out) :: shiftPerAtom(:) 288 289 !> Shift per shell. 290 real(dp), intent(out) :: shiftPerShell(:,:) 291 292 @:ASSERT(size(shiftPerAtom) == this%nAtoms) 293 @:ASSERT(size(shiftPerShell, dim=1) == this%mShellsReal) 294 295 if (this%shellResolved) then 296 shiftPerAtom(:) = this%shift3 297 shiftPerShell(:,:) = this%shift1 + this%shift2 298 else 299 shiftPerAtom(:) = this%shift1(1,:) + this%shift2(1,:) + this%shift3 300 shiftPerShell(:,:) = 0.0_dp 301 end if 302 303 end subroutine getShifts 304 305 306 !> Returns energy per atom. 307 subroutine getEnergyPerAtom(this, energyPerAtom) 308 class(ThirdOrder), intent(inout) :: this 309 real(dp), intent(out) :: energyPerAtom(:) 310 311 @:ASSERT(size(energyPerAtom) == this%nAtoms) 312 313 energyPerAtom(:) = (1.0_dp / 3.0_dp) * (& 314 & sum((this%shift1 + this%shift2) * this%chargesPerShell, dim=1)& 315 & + this%shift3 * this%chargesPerAtom) 316 317 end subroutine getEnergyPerAtom 318 319 320 !> Returns the energy per atom for linearized 3rd order Hamiltonian. 321 !> Note: When using the linear XLBOMD form, charges should not be updated via updateCharges after 322 !> the diagonalization, so that the shift vectors remain the ones built with the input 323 !> charges. However, since for calculating energy/forces, the output charges are needed, they must 324 !> be passed explicitely here. 325 subroutine getEnergyPerAtomXlbomd(this, qOut, q0, species, orb, energyPerAtom) 326 327 !> Instance. 328 class(ThirdOrder), intent(inout) :: this 329 330 !> Output populations determined after the diagonalization. 331 real(dp), intent(in) :: qOut(:,:,:) 332 333 !> Reference populations. 334 real(dp), intent(in) :: q0(:,:,:) 335 336 !> Species of each atom. 337 integer, intent(in) :: species(:) 338 339 !> Orbital information 340 type(TOrbitals), intent(in) :: orb 341 342 !> Energy per atom for linearized case. 343 real(dp), intent(out) :: energyPerAtom(:) 344 345 real(dp), allocatable :: qOutAtom(:), qOutShell(:,:), qOutShellTmp(:,:) 346 347 allocate(qOutAtom(this%nAtoms)) 348 allocate(qOutShell(this%mShells, this%nAtoms)) 349 if (this%shellResolved) then 350 call getSummedCharges(species, orb, qOut, q0, dQAtom=qOutAtom, dQShell=qOutShell) 351 else 352 ! First (only) component of qOutShell contains atomic charge 353 allocate(qOutShellTmp(this%mShellsReal, this%nAtoms)) 354 call getSummedCharges(species, orb, qOut, q0, dQAtom=qOutAtom, dQShell=qOutShellTmp) 355 qOutShell(1,:) = sum(qOutShellTmp, dim=1) 356 end if 357 energyPerAtom(:) = sum(this%shift1 * qOutShell, dim=1)& 358 & + sum(this%shift2 * (qOutShell - this%chargesPerShell), dim=1)& 359 & + this%shift3 * (qOutAtom - this%chargesPerAtom) 360 361 end subroutine getEnergyPerAtomXlbomd 362 363 364 !> Add gradient component resulting from the derivative of the potential. 365 subroutine addGradientDc(this, neighList, species, coords, img2CentCell, derivs) 366 367 !> Instance. 368 class(ThirdOrder), intent(inout) :: this 369 370 !> Neighbour list. 371 type(TNeighbourList), intent(in) :: neighList 372 373 !> Specie for each atom. 374 integer, intent(in) :: species(:) 375 376 !> Coordinate of each atom. 377 real(dp), intent(in) :: coords(:,:) 378 379 !> Mapping of atoms to cetnral cell. 380 integer, intent(in) :: img2CentCell(:) 381 382 !> Gradient on exit. 383 real(dp), intent(inout) :: derivs(:,:) 384 385 integer :: iAt1, iAt2, iAt2f, iSp1, iSp2, iSh1, iSh2, iNeigh 386 real(dp) :: rab, tmp, tmp3(3) 387 logical :: damping 388 389 do iAt1 = 1, this%nAtoms 390 iSp1 = species(iAt1) 391 do iNeigh = 1, this%nNeighMax(iAt1) 392 iAt2 = neighList%iNeighbour(iNeigh, iAt1) 393 iAt2f = img2CentCell(iAt2) 394 iSp2 = species(iAt2f) 395 if (iAt1 == iAt2f .or. iNeigh > this%nNeigh(iSp2, iAt1)) then 396 cycle 397 end if 398 rab = sqrt(neighList%neighDist2(iNeigh, iAt1)) 399 damping = this%damped(iSp1) .or. this%damped(iSp2) 400 tmp = 0.0_dp 401 do iSh1 = 1, this%nShells(iSp1) 402 do iSh2 = 1, this%nShells(iSp2) 403 tmp = tmp + this%chargesPerShell(iSh1, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 404 & * (this%chargesPerAtom(iAt1)& 405 & * gamma3pR(this%UU(iSh1, iSp1), this%UU(iSh2, iSp2),& 406 & this%dUdQ(iSh1, iSp1), rab, damping, this%dampExp)& 407 & + this%chargesPerAtom(iAt2f)& 408 & * gamma3pR(this%UU(iSh2, iSp2), this%UU(iSh1, iSp1),& 409 & this%dUdQ(iSh2, iSp2), rab, damping, this%dampExp)) 410 end do 411 end do 412 tmp3(:) = tmp / (3.0_dp * rab) * (coords(:, iAt1) - coords(:, iAt2)) 413 derivs(:, iAt1) = derivs(:, iAt1) + tmp3 414 derivs(:, iAt2f) = derivs(:, iAt2f) - tmp3 415 end do 416 end do 417 418 end subroutine addGradientDc 419 420 421 !> Add gradient component resulting from the derivative of the potential for the linearized 422 !> (XLBOMD) case. 423 subroutine addGradientDcXlbomd(this, neighList, species, coords, img2CentCell, qOut, q0, orb,& 424 & derivs) 425 426 !> Instance. 427 class(ThirdOrder), intent(inout) :: this 428 429 !> Neighbour list. 430 type(TNeighbourList), intent(in) :: neighList 431 432 !> Specie for each atom. 433 integer, intent(in) :: species(:) 434 435 !> Coordinate of each atom. 436 real(dp), intent(in) :: coords(:,:) 437 438 !> Mapping of atoms to cetnral cell. 439 integer, intent(in) :: img2CentCell(:) 440 441 !> Output populations determined after the diagonalization. 442 real(dp), intent(in) :: qOut(:,:,:) 443 444 !> Reference populations. 445 real(dp), intent(in) :: q0(:,:,:) 446 447 !> Orbital information 448 type(TOrbitals), intent(in) :: orb 449 450 !> Modified gradient on exit. 451 real(dp), intent(inout) :: derivs(:,:) 452 453 real(dp), allocatable :: qOutAtom(:), qOutShell(:,:), qOutShellTmp(:,:) 454 real(dp), allocatable :: qDiffShell(:,:), qDiffAtom(:) 455 integer :: iAt1, iAt2, iAt2f, iSp1, iSp2, iSh1, iSh2, iNeigh 456 real(dp) :: gammaDeriv1, gammaderiv2, rab, tmp 457 real(dp) :: tmp3(3) 458 logical :: damping 459 460 allocate(qOutAtom(this%nAtoms)) 461 allocate(qOutShell(this%mShells, this%nAtoms)) 462 allocate(qDiffAtom(this%nAtoms)) 463 allocate(qDiffShell(this%mShells, this%nAtoms)) 464 if (this%shellResolved) then 465 call getSummedCharges(species, orb, qOut, q0, dQAtom=qOutAtom, dQShell=qOutShell) 466 else 467 ! First (only) component of qOutShell contains atomic charge 468 allocate(qOutShellTmp(this%mShellsReal, this%nAtoms)) 469 call getSummedCharges(species, orb, qOut, q0, dQAtom=qOutAtom, dQShell=qOutShellTmp) 470 qOutShell(1,:) = sum(qOutShellTmp, dim=1) 471 end if 472 473 qDiffAtom(:) = qOutAtom - this%chargesPerAtom 474 qDiffShell(:,:) = qOutShell - this%chargesPerShell 475 do iAt1 = 1, this%nAtoms 476 iSp1 = species(iAt1) 477 do iNeigh = 1, this%nNeighMax(iAt1) 478 iAt2 = neighList%iNeighbour(iNeigh, iAt1) 479 iAt2f = img2CentCell(iAt2) 480 iSp2 = species(iAt2f) 481 if (iAt1 == iAt2f .or. iNeigh > this%nNeigh(iSp2, iAt1)) then 482 cycle 483 end if 484 rab = sqrt(neighList%neighDist2(iNeigh, iAt1)) 485 damping = this%damped(iSp1) .or. this%damped(iSp2) 486 tmp = 0.0_dp 487 do iSh1 = 1, this%nShells(iSp1) 488 do iSh2 = 1, this%nShells(iSp2) 489 gammaDeriv1 =& 490 & gamma3pR(this%UU(iSh1, iSp1), this%UU(iSh2, iSp2),& 491 & this%dUdQ(iSh1, iSp1), rab, damping, this%dampExp) 492 gammaDeriv2 =& 493 & gamma3pR(this%UU(iSh2, iSp2), this%UU(iSh1, iSp1),& 494 & this%dUdQ(iSh2, iSp2), rab, damping, this%dampExp) 495 tmp = tmp + gammaDeriv1& 496 & * (this%chargesPerAtom(iAt1) * this%chargesPerShell(iSh2, iAt2f)& 497 & * qOutShell(iSh1, iAt1)& 498 & + this%chargesPerAtom(iAt1) * this%chargesPerShell(iSh1, iAt1)& 499 & * qDiffShell(iSh2, iAt2f)& 500 & + this%chargesPerShell(iSh1, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 501 & * qDiffAtom(iAt1)) 502 tmp = tmp + gammaDeriv2& 503 & * (this%chargesPerAtom(iAt2f) * this%chargesPerShell(iSh1, iAt1)& 504 & * qOutShell(iSh2, iAt2f)& 505 & + this%chargesPerAtom(iAt2f) * this%chargesPerShell(iSh2, iAt2f)& 506 & * qDiffShell(iSh1, iAt1)& 507 & + this%chargesPerShell(iSh2, iAt2f) * this%chargesPerShell(iSh1, iAt1)& 508 & * qDiffAtom(iAt2f)) 509 end do 510 end do 511 tmp3 = tmp / (3.0_dp * rab) * (coords(:, iAt1) - coords(:, iAt2)) 512 derivs(:, iAt1) = derivs(:, iAt1) + tmp3 513 derivs(:, iAt2f) = derivs(:, iAt2f) - tmp3 514 end do 515 end do 516 517 end subroutine addGradientDcXlbomd 518 519 520 !> Add stress component resulting from the derivative of the potential. 521 subroutine addStressDc(this, neighList, species, coords, img2CentCell, cellVol, st) 522 523 !> Instance. 524 class(ThirdOrder), intent(inout) :: this 525 526 !> Neighbour list. 527 type(TNeighbourList), intent(in) :: neighList 528 529 !> Specie for each atom. 530 integer, intent(in) :: species(:) 531 532 !> Coordinate of each atom. 533 real(dp), intent(in) :: coords(:,:) 534 535 !> Mapping of atoms to cetnral cell. 536 integer, intent(in) :: img2CentCell(:) 537 538 !> cell volume. 539 real(dp), intent(in) :: cellVol 540 541 !> Gradient on exit. 542 real(dp), intent(inout) :: st(:,:) 543 544 integer :: iAt1, iAt2, iAt2f, iSp1, iSp2, iSh1, iSh2, iNeigh, ii 545 real(dp) :: rab, tmp, tmp3(3), stTmp(3,3), prefac, vect(3) 546 logical :: damping 547 548 stTmp(:,:) = 0.0_dp 549 do iAt1 = 1, this%nAtoms 550 iSp1 = species(iAt1) 551 do iNeigh = 1, this%nNeighMax(iAt1) 552 iAt2 = neighList%iNeighbour(iNeigh, iAt1) 553 iAt2f = img2CentCell(iAt2) 554 iSp2 = species(iAt2f) 555 if (iNeigh > this%nNeigh(iSp2, iAt1)) then 556 cycle 557 end if 558 vect(:) = coords(:,iAt1) - coords(:,iAt2) 559 if (iAt1 == iAt2f) then 560 prefac = 0.5_dp 561 else 562 prefac = 1.0_dp 563 end if 564 rab = sqrt(neighList%neighDist2(iNeigh, iAt1)) 565 damping = this%damped(iSp1) .or. this%damped(iSp2) 566 tmp = 0.0_dp 567 do iSh1 = 1, this%nShells(iSp1) 568 do iSh2 = 1, this%nShells(iSp2) 569 tmp = tmp + this%chargesPerShell(iSh1, iAt1) * this%chargesPerShell(iSh2, iAt2f)& 570 & * (this%chargesPerAtom(iAt1)& 571 & * gamma3pR(this%UU(iSh1, iSp1), this%UU(iSh2, iSp2),& 572 & this%dUdQ(iSh1, iSp1), rab, damping, this%dampExp)& 573 & + this%chargesPerAtom(iAt2f)& 574 & * gamma3pR(this%UU(iSh2, iSp2), this%UU(iSh1, iSp1),& 575 & this%dUdQ(iSh2, iSp2), rab, damping, this%dampExp)) 576 end do 577 end do 578 tmp3(:) = tmp / (3.0_dp * rab) * (coords(:, iAt1) - coords(:, iAt2)) 579 do ii = 1, 3 580 stTmp(:, ii) = stTmp(:, ii) + prefac * tmp3(:) * vect(ii) 581 end do 582 end do 583 end do 584 585 st = st - stTmp / cellVol 586 587 end subroutine addStressDc 588 589 590! Private routines 591 592 593 !> calculate short range cut off distance 594 subroutine calcCutoffs(hubbUs, nShells, cutoffs) 595 596 !> Hubard U values 597 real(dp), intent(in) :: hubbUs(:,:) 598 599 !> Shells on atoms 600 integer, intent(in) :: nShells(:) 601 602 !> resulting cutoff distances 603 real(dp), intent(out) :: cutoffs(:,:) 604 605 integer :: nSpecies 606 real(dp) :: cutoff 607 real(dp), allocatable :: minUs(:) 608 integer :: iSp1, iSp2 609 610 nSpecies = size(cutoffs, dim=1) 611 allocate(minUs(nSpecies)) 612 do iSp1 = 1, nSpecies 613 minUs(iSp1) = minval(hubbUs(1:nShells(iSp1), iSp1)) 614 end do 615 do iSp1 = 1, nSpecies 616 do iSp2 = iSp1, nSpecies 617 cutoff = expGammaCutoff(minUs(iSp2), minUs(iSp1)) 618 cutoffs(iSp2, iSp1) = cutoff 619 cutoffs(iSp1, iSp2) = cutoff 620 end do 621 end do 622 623 end subroutine calcCutoffs 624 625 626 !> Gamma_AB = dgamma_AB/dUa * (dUa/dQa) 627 function gamma3(Ua, Ub, dUa, rab, damping, xi) result(res) 628 real(dp), intent(in) :: Ua 629 real(dp), intent(in) :: Ub 630 real(dp), intent(in) :: dUa 631 real(dp), intent(in) :: rab 632 logical, intent(in) :: damping 633 real(dp), intent(in) :: xi 634 real(dp) :: res 635 636 res = gamma2pU(Ua, Ub, rab, damping, xi) * dUa 637 638 end function gamma3 639 640 641 !> dGamma_AB/dr 642 function gamma3pR(Ua, Ub, dUa, rab, damping, xi) result(res) 643 real(dp), intent(in) :: Ua, Ub, dUa, rab 644 logical, intent(in) :: damping 645 real(dp), intent(in) :: xi 646 real(dp) :: res 647 648 res = gamma2pUpR(Ua, Ub, rab, damping, xi) * dUa 649 650 end function gamma3pR 651 652 653 !> dgamma_AB/dUa 654 !> Sign convention: routine delivers dgamma_AB/dUa with the right sign. 655 !> Energy contributions must be therefore summed with *positive* sign. 656 function gamma2pU(Ua, Ub, rab, damping, xi) result(res) 657 real(dp), intent(in) :: Ua, Ub, rab 658 logical, intent(in) :: damping 659 real(dp), intent(in) :: xi 660 real(dp) :: res 661 662 real(dp) :: tauA, tauB, tau, uu 663 664 tauA = 3.2_dp * Ua ! 16/5 * Ua 665 tauB = 3.2_dp * Ub 666 667 if (rab < tolSameDist) then 668 if (abs(Ua - Ub) < minHubDiff) then 669 ! Limiting case for dG/dU with Ua=Ub and Rab = 0 670 res = 0.5_dp 671 else 672 res = dGdUr0(tauA, tauB) 673 end if 674 else if (abs(Ua - Ub) < minHubDiff) then 675 tau = 0.5_dp * (tauA + tauB) 676 res = -3.2_dp * shortpT_2(tau, rab) 677 if (damping) then 678 uu = 0.5_dp * (Ua + Ub) 679 res = res * hh(uu, uu, rab, xi) - short_2(tau, rab) * hpU(uu, uu, rab, xi) 680 end if 681 else 682 res = -3.2_dp * shortpT_1(tauA, tauB, rab) 683 if (damping) then 684 res = res * hh(Ua, Ub, rab, xi) - short_1(tauA, tauB, rab) * hpU(Ua, Ub, rab, xi) 685 end if 686 end if 687 688 end function gamma2pU 689 690 691 !> d^2gamma_AB/dUa*dr 692 !> Sign convention: routine delivers d^2gamma_AB/dUa*dr with the right sign. 693 !> Gradient contributions must be therefore summed with *positive* sign. 694 function gamma2pUpR(Ua, Ub, rab, damping, xi) result(res) 695 real(dp), intent(in) :: Ua, Ub, rab 696 logical, intent(in) :: damping 697 real(dp), intent(in) :: xi 698 real(dp) :: res 699 700 real(dp) :: tauA, tauB, tau, uu 701 702 tauA = 3.2_dp * Ua 703 tauB = 3.2_dp * Ub 704 705 if (rab < tolSameDist) then 706 res = 0.0_dp 707 else if (abs(Ua - Ub) < minHubDiff) then 708 tau = 0.5_dp * (tauA + tauB) 709 res = -3.2_dp * shortpTpR_2(tau, rab) 710 if (damping) then 711 uu = 0.5_dp * (Ua + Ub) 712 res = res * hh(uu, uu, rab, xi)& 713 & - 3.2_dp * shortpT_2(tau, rab) * hpR(uu, uu, rab, xi)& 714 & - shortpR_2(tau, rab) * hpU(uu, uu, rab, xi)& 715 & - short_2(tau, rab) * hpUpR(uu, uu, rab, xi) 716 end if 717 else 718 res = -3.2_dp *shortpTpR_1(tauA, tauB, rab) 719 if (damping) then 720 res = res * hh(Ua, Ub, rab, xi)& 721 & - 3.2_dp * shortpT_1(tauA, tauB, rab) * hpR(Ua, Ub, rab, xi)& 722 & - shortpR_1(tauA, tauB, rab) * hpU(Ua, Ub, rab, xi)& 723 & - short_1(tauA, tauB, rab) * hpUpR(Ua, Ub, rab, xi) 724 end if 725 end if 726 727 end function gamma2pUpR 728 729 730 !> \frac{d\gamma}{dU_{l_a}} for r = 0 731 !> Eq S7 in Gaus et al. (2015) JCTC 11:4205-4219, DOI: 10.1021/acs.jctc.5b00600 732 function dGdUr0(tauA, tauB) result(res) 733 real(dp), intent(in) :: tauA, tauB 734 real(dp) :: res 735 736 real(dp) :: invTauSum 737 738 invTauSum = 1.0_dp / (tauA + tauB) 739 740 res = 1.6_dp * invTauSum * (tauB& 741 & + invTauSum * (-tauA * tauB& 742 & + invTauSum * (2.0_dp * tauA * tauB**2& 743 & + invTauSum * (-3.0_dp * tauA**2 * tauB**2)))) 744 745 end function dGdUr0 746 747 748 !> S1(tauA,tauB,r): Short range SCC when tauA <> tauB and r <> 0 749 function short_1(tauA, tauB, rab) result(res) 750 real(dp), intent(in) :: tauA, tauB, rab 751 real(dp) :: res 752 753 res = exp(-tauA * rab) * ff(tauA, tauB, rab) + exp(-tauB * rab) * ff(tauB, tauA, rab) 754 755 end function short_1 756 757 758 !> S2(tau,r), short range SCC when tauA = tauB = tau and r <> 0. 759 function short_2(tau, rab) result(res) 760 real(dp), intent(in) :: tau, rab 761 real(dp) :: res 762 763 res = exp(-tau * rab) * gg(tau, rab) 764 765 end function short_2 766 767 768 !> dS1(tauA,tauB,r)/dtauA 769 function shortpT_1(tauA, tauB, rab) result(res) 770 real(dp), intent(in) :: tauA, tauB, rab 771 real(dp) :: res 772 773 res = exp(-tauA * rab) * (fpT1(tauA, tauB, rab) - rab * ff(tauA, tauB, rab))& 774 & + exp(-tauB * rab) * fpT2(tauB, tauA, rab) 775 776 end function shortpT_1 777 778 779 !> dS2(tauA,tauB,r)/dtauA 780 function shortpT_2(tau, rab) result(res) 781 real(dp), intent(in) :: tau, rab 782 real(dp) :: res 783 784 res = exp(-tau * rab) * (gpT(tau, rab) - rab * gg(tau, rab)) 785 786 end function shortpT_2 787 788 789 !> dS1(tauA,tauB,r)/dr 790 function shortpR_1(tauA, tauB, rab) result(res) 791 real(dp), intent(in) :: tauA, tauB, rab 792 real(dp) :: res 793 794 res = exp(-tauA * rab) * (fpR(tauA, tauB, rab) - tauA * ff(tauA, tauB, rab))& 795 & + exp(-tauB * rab) * (fpR(tauB, tauA, rab) - tauB * ff(tauB, tauA, rab)) 796 797 end function shortpR_1 798 799 800 !> dS2(tauA,tauB,r)/dr 801 function shortpR_2(tau, rab) result(res) 802 real(dp), intent(in) :: tau, rab 803 real(dp) :: res 804 805 res = exp(-tau * rab) * (gpR(tau, rab) - tau * gg(tau, rab)) 806 807 end function shortpR_2 808 809 810 !> d^2S1(tauA,tauB,r)/dtauA*dr 811 function shortpTpR_1(tauA, tauB, rab) result(res) 812 real(dp), intent(in) :: tauA, tauB, rab 813 real(dp) :: res 814 815 res = exp(-tauA * rab) * (ff(tauA, tauB, rab) * (tauA * rab - 1.0_dp)& 816 & - tauA * fpT1(tauA, tauB, rab) + fpT1pR(tauA, tauB, rab) - rab * fpR(tauA, tauB, rab))& 817 & + exp(-tauB * rab) * (fpT2pR(tauB, tauA, rab) - tauB * fpT2(tauB, tauA, rab)) 818 819 end function shortpTpR_1 820 821 822 !> d^2S2(tau,r)/dtau*dr 823 function shortpTpR_2(tau, rab) result(res) 824 real(dp), intent(in) :: tau, rab 825 real(dp) :: res 826 827 res = exp(-tau * rab) * ((tau * rab - 1.0_dp) * gg(tau, rab) - tau * gpT(tau, rab)& 828 & + gpTpR(tau, rab) - rab * gpR(tau, rab)) 829 830 end function shortpTpR_2 831 832 833 !> f(tauA,tauB,r) 834 function ff(tauA, tauB, rab) result(res) 835 real(dp), intent(in) :: tauA, tauB, rab 836 real(dp) :: res 837 838 res = 0.5_dp * tauA * tauB**4 / (tauA**2 - tauB**2)**2& 839 & - (tauB**6 - 3.0_dp * tauA**2 * tauB**4) / ((tauA**2 - tauB**2)**3 * rab) 840 841 end function ff 842 843 844 !> df(tauA,tauB,r)/dtauA 845 function fpT1(tauA, tauB, rab) result(res) 846 real(dp), intent(in) :: tauA, tauB, rab 847 real(dp) :: res 848 849 res = -0.5_dp * (tauB**6 + 3.0_dp * tauA**2 * tauB**4) / (tauA**2 - tauB**2)**3& 850 & - 12.0_dp * tauA**3 * tauB**4 / ((tauA**2 - tauB**2)**4 * rab) 851 852 end function fpT1 853 854 855 !> df(tauB,tauA,rab)/dtauA 856 function fpT2(tauB, tauA , rab) result(res) 857 real(dp), intent(in) :: tauA, tauB, rab 858 real(dp) :: res 859 860 res = 2.0_dp * tauB**3 * tauA**3 / (tauB**2 - tauA**2)**3& 861 & + 12.0_dp * tauB**4 * tauA**3 / ((tauB**2 - tauA**2)**4 * rab) 862 863 end function fpT2 864 865 866 !> df(tauA, tauB,r)/dr 867 function fpR(tauA, tauB, rab) result(res) 868 real(dp), intent(in) :: tauA, tauB, rab 869 real(dp) :: res 870 871 res = (tauB**6 - 3.0_dp * tauB**4 * tauA**2) / (rab**2 * (tauA**2 - tauB**2)**3) 872 873 end function fpR 874 875 876 !> d^2f(tauA,tauB,r)/dtauA*dr 877 function fpT1pR(tauA, tauB, rab) result(res) 878 real(dp), intent(in) :: tauA, tauB, rab 879 real(dp) :: res 880 881 res = 12.0_dp * tauA**3 * tauB**4 / (rab**2 * (tauA**2 - tauB**2)**4) 882 883 end function fpT1pR 884 885 886 !> d^2f(tauB,tauA,r)/dtauA*dr 887 function fpT2pR(tauB, tauA, rab) result(res) 888 real(dp), intent(in) :: tauB, tauA, rab 889 real(dp) :: res 890 891 res = -12.0_dp * tauA**3 * tauB**4 / (rab**2 * (tauA**2 - tauB**2)**4) 892 893 end function fpT2pR 894 895 896 !> g(tau,r) 897 function gg(tau, rab) result(res) 898 real(dp), intent(in) :: tau, rab 899 real(dp) :: res 900 901 res = 1.0_dp / (48.0_dp * rab) * (48.0_dp + 33.0_dp * tau * rab + 9.0_dp * tau**2 * rab**2& 902 & + tau**3 * rab**3) 903 904 end function gg 905 906 907 !> dg(tau,rab)/dtau 908 function gpT(tau, rab) result(res) 909 real(dp), intent(in) :: tau, rab 910 real(dp) :: res 911 912 res = 1.0_dp / 48.0_dp * (33.0_dp + 18.0_dp * tau * rab + 3.0_dp * tau**2 * rab**2) 913 914 end function gpT 915 916 917 !> dg(tau,r)/dr 918 function gpR(tau, rab) result(res) 919 real(dp), intent(in) :: tau, rab 920 real(dp) :: res 921 922 res = (-48.0_dp + 9.0_dp * (tau * rab)**2 + 2.0_dp * (tau * rab)**3) / (48.0_dp * rab**2) 923 924 end function gpR 925 926 927 !> d^2g(tau,r)/dtau*dr 928 function gpTpR(tau, rab) result(res) 929 real(dp), intent(in) :: tau, rab 930 real(dp) :: res 931 932 res = (3.0_dp * tau + tau**2 * rab) / 8.0_dp 933 934 end function gpTpR 935 936 937 !> Damping: h(Ua,Ub) 938 function hh(Ua, Ub, rab, xi) result(res) 939 real(dp), intent(in) :: Ua, Ub, rab, xi 940 real(dp) :: res 941 942 res = exp(-(0.5_dp * (Ua + Ub))**xi * rab**2) 943 944 end function hh 945 946 947 !> dh(Ua,Ub)/dUa 948 function hpU(Ua, Ub, rab, xi) result(res) 949 real(dp), intent(in) :: Ua, Ub, rab, xi 950 real(dp) :: res 951 952 res = -0.5_dp * xi * rab**2 * (0.5_dp * (Ua + Ub))**(xi - 1.0_dp) * hh(Ua, Ub, rab, xi) 953 954 end function hpU 955 956 957 !> dh(Ua,Ub)/dr 958 function hpR(Ua, Ub, rab, xi) result(res) 959 real(dp), intent(in) :: Ua, Ub, rab, xi 960 real(dp) :: res 961 962 res = -2.0_dp * rab * (0.5_dp * (Ua + Ub))**xi * hh(Ua, Ub, rab, xi) 963 964 end function hpR 965 966 967 !> dh(Ua,Ub)/dUa*dr 968 function hpUpR(Ua, Ub, rab, xi) result(res) 969 real(dp), intent(in) :: Ua, Ub, rab, xi 970 real(dp) :: res 971 972 res = xi * rab * (0.5_dp * (Ua + Ub))**(xi - 1.0_dp)& 973 & * (rab**2 * (0.5_dp * (Ua + Ub))**xi - 1.0_dp) * hh(Ua, Ub, rab, xi) 974 975 end function hpUpR 976 977end module dftbp_thirdorder_module 978