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