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!> Code to calculate forces for several different types of calculation (non-scc, scc, sDFTB etc)
11module dftbp_forces
12  use dftbp_assert
13  use dftbp_accuracy
14  use dftbp_nonscc, only : NonSccDiff
15  use dftbp_scc
16  use dftbp_commontypes
17  use dftbp_slakocont
18  use dftbp_schedule
19  use dftbp_environment
20  implicit none
21
22  private
23
24  public :: derivative_shift
25
26
27  !> forces with shift vectors present
28  interface derivative_shift
29
30    !> derivatives without any shift
31    module procedure derivative_nonSCC
32
33    !> derivatives with shift
34    module procedure derivative_block
35
36    !> derivatives with complex shift
37    module procedure derivative_iBlock
38
39  end interface derivative_shift
40
41contains
42
43
44  !> The non-SCC electronic force contribution for all atoms from the matrix derivatives and the
45  !> density and energy-density matrices
46  subroutine derivative_nonSCC(env, deriv, derivator, DM, EDM, skHamCont, skOverCont, coords,&
47      & species, iNeighbour, nNeighbourSK, img2CentCell, iPair, orb)
48
49    !> Computational environment settings
50    type(TEnvironment), intent(in) :: env
51
52    !> x,y,z derivatives for each real atom in the system
53    real(dp), intent(out) :: deriv(:,:)
54
55    !> Differentiatior for the non-scc components
56    class(NonSccDiff), intent(in) :: derivator
57
58    !> density matrix in packed format
59    real(dp), intent(in) :: DM(:)
60
61    !> energy-weighted density matrix in packed format
62    real(dp), intent(in) :: EDM(:)
63
64    !> Container for SK Hamiltonian integrals
65    type(OSlakoCont), intent(in) :: skHamCont
66
67    !> Container for SK overlap integrals
68    type(OSlakoCont), intent(in) :: skOverCont
69
70    !> list of all atomic coordinates
71    real(dp), intent(in) :: coords(:,:)
72
73    !> list of all atomic species
74    integer, intent(in) :: species(:)
75
76    !> neighbour list for atoms
77    integer, intent(in) :: iNeighbour(0:,:)
78
79    !> number of neighbours of each atom
80    integer, intent(in) :: nNeighbourSK(:)
81
82    !> indexing array for periodic image atoms
83    integer, intent(in) :: img2CentCell(:)
84
85    !> indexing array for the Hamiltonian
86    integer, intent(in) :: iPair(0:,:)
87
88    !> Information about the shells and orbitals in the system.
89    type(TOrbitals), intent(in) :: orb
90
91    integer :: iOrig, ii
92    integer :: nAtom, iNeigh, iAtom1, iAtom2, iAtom2f
93    integer :: nOrb1, nOrb2
94    real(dp) :: sqrDMTmp(orb%mOrb,orb%mOrb), sqrEDMTmp(orb%mOrb,orb%mOrb)
95    real(dp) :: hPrimeTmp(orb%mOrb,orb%mOrb,3), sPrimeTmp(orb%mOrb,orb%mOrb,3)
96    integer :: iAtFirst, iAtLast
97
98    @:ASSERT(size(deriv,dim=1) == 3)
99
100    nAtom = size(orb%nOrbAtom)
101    deriv(:,:) = 0.0_dp
102
103    call distributeRangeInChunks(env, 1, nAtom, iAtFirst, iAtLast)
104
105    !$OMP PARALLEL DO PRIVATE(iAtom1,nOrb1,iNeigh,iAtom2,iAtom2f,nOrb2,iOrig,sqrDMTmp,sqrEDMTmp, &
106    !$OMP& hPrimeTmp,sPrimeTmp,ii) DEFAULT(SHARED) SCHEDULE(RUNTIME) REDUCTION(+:deriv)
107    do iAtom1 = iAtFirst, iAtLast
108      nOrb1 = orb%nOrbAtom(iAtom1)
109      !! loop from 1 as no contribution from the atom itself
110      do iNeigh = 1, nNeighbourSK(iAtom1)
111        iAtom2 = iNeighbour(iNeigh, iAtom1)
112        iAtom2f = img2CentCell(iAtom2)
113        if (iAtom1 /= iAtom2f) then
114          nOrb2 = orb%nOrbAtom(iAtom2f)
115          iOrig = iPair(iNeigh,iAtom1)
116          sqrDMTmp(:,:) = 0.0_dp
117          sqrEDMTmp(:,:) = 0.0_dp
118          hPrimeTmp(:,:,:) = 0.0_dp
119          sPrimeTmp(:,:,:) = 0.0_dp
120          sqrDMTmp(1:nOrb2,1:nOrb1) = reshape(DM(iOrig+1:iOrig+nOrb1*nOrb2), (/nOrb2,nOrb1/))
121          sqrEDMTmp(1:nOrb2,1:nOrb1) = reshape(EDM(iOrig+1:iOrig+nOrb1*nOrb2), (/nOrb2,nOrb1/))
122          call derivator%getFirstDeriv(hPrimeTmp, skHamCont, coords, species, iAtom1, iAtom2, orb)
123          call derivator%getFirstDeriv(sPrimeTmp, skOverCont, coords, species, iAtom1, iAtom2, orb)
124          ! note factor of 2 for implicit summation over lower triangle of density matrix:
125          do ii = 1, 3
126            deriv(ii,iAtom1) = deriv(ii,iAtom1)&
127                & + sum(sqrDMTmp(1:nOrb2,1:nOrb1) * 2.0_dp*hPrimeTmp(1:nOrb2,1:nOrb1,ii))&
128                & - sum(sqrEDMTmp(1:nOrb2,1:nOrb1) * 2.0_dp*sPrimeTmp(1:nOrb2,1:nOrb1,ii))
129          end do
130          ! Add contribution to the force from atom 1 onto atom 2f using the symmetry in the blocks,
131          ! and note that the skew symmetry in the derivatives is being used
132          do ii = 1, 3
133            deriv(ii,iAtom2f) = deriv(ii,iAtom2f)&
134                & - sum(sqrDMTmp(1:nOrb2,1:nOrb1) * 2.0_dp*hPrimeTmp(1:nOrb2,1:nOrb1,ii))&
135                & + sum(sqrEDMTmp(1:nOrb2,1:nOrb1) * 2.0_dp*sPrimeTmp(1:nOrb2,1:nOrb1,ii))
136          end do
137        end if
138      end do
139    end do
140    !$OMP END PARALLEL DO
141
142    call assembleChunks(env, deriv)
143
144  end subroutine derivative_nonSCC
145
146
147  !> The SCC and spin electronic force contribution for all atoms from the matrix derivatives, self
148  !> consistent potential and the density and energy-density matrices
149  subroutine derivative_block(env, deriv, derivator, DM, EDM, skHamCont, skOverCont, coords,&
150      & species, iNeighbour, nNeighbourSK, img2CentCell, iPair, orb, shift)
151
152    !> Computational environment settings
153    type(TEnvironment), intent(in) :: env
154
155    !> x,y,z derivatives for each real atom in the system
156    real(dp), intent(out) :: deriv(:,:)
157
158    !> Differentiatior for the non-scc components
159    class(NonSccDiff), intent(in) :: derivator
160
161    !> density matrix in packed format
162    real(dp), intent(in) :: DM(:,:)
163
164    !> energy-weighted density matrix in packed format
165    real(dp), intent(in) :: EDM(:)
166
167    !> Container for SK Hamiltonian integrals
168    type(OSlakoCont) :: skHamCont
169
170    !> Container for SK overlap integrals
171    type(OSlakoCont) :: skOverCont
172
173    !> list of all atomic coordinates
174    real(dp), intent(in) :: coords(:,:)
175
176    !> list of all atomic species
177    integer, intent(in) :: species(:)
178
179    !> neighbour list for atoms
180    integer, intent(in) :: iNeighbour(0:,:)
181
182    !> number of neighbours of each atom
183    integer, intent(in) :: nNeighbourSK(:)
184
185    !> indexing array for periodic image atoms
186    integer, intent(in) :: img2CentCell(:)
187
188    !> indexing array for the Hamiltonian
189    integer, intent(in) :: iPair(0:,:)
190
191    !> Information about the shells and orbitals in the system.
192    type(TOrbitals), intent(in) :: orb
193
194    !> block shift from the potential
195    real(dp), intent(in) :: shift(:,:,:,:)
196
197    integer :: iOrig, iSpin, ii, nSpin, nAtom
198    integer :: iNeigh, iAtom1, iAtom2, iAtom2f, iSp1, iSp2
199    integer :: nOrb1, nOrb2
200
201    real(dp) :: sqrDMTmp(orb%mOrb,orb%mOrb), sqrEDMTmp(orb%mOrb,orb%mOrb)
202    real(dp) :: shiftSprime(orb%mOrb,orb%mOrb)
203    real(dp) :: hPrimeTmp(orb%mOrb,orb%mOrb,3), sPrimeTmp(orb%mOrb,orb%mOrb,3)
204    real(dp) :: derivTmp(3)
205    integer :: iAtFirst, iAtLast
206
207    nAtom = size(orb%nOrbAtom)
208    nSpin = size(shift,dim=4)
209    @:ASSERT(nSpin == 1 .or. nSpin == 2 .or. nSpin ==4)
210    @:ASSERT(size(deriv,dim=1) == 3)
211    @:ASSERT(size(deriv,dim=2)==nAtom)
212    @:ASSERT(size(DM,dim=1)==size(EDM,dim=1))
213    @:ASSERT(size(shift,dim=1)==orb%mOrb)
214    @:ASSERT(size(shift,dim=2)==orb%mOrb)
215    @:ASSERT(size(shift,dim=3)==nAtom)
216    @:ASSERT(size(DM,dim=2)==nSpin)
217
218    deriv(:,:) = 0.0_dp
219
220    call distributeRangeInChunks(env, 1, nAtom, iAtFirst, iAtLast)
221
222    !$OMP PARALLEL DO PRIVATE(iAtom1,iSp1,nOrb1,iNeigh,iAtom2,iAtom2f,iSp2,nOrb2,iOrig,sqrDMTmp, &
223    !$OMP& sqrEDMTmp,hPrimeTmp,sPrimeTmp,derivTmp,shiftSprime,iSpin,ii) DEFAULT(SHARED) &
224    !$OMP& SCHEDULE(RUNTIME) REDUCTION(+:deriv)
225    do iAtom1 = iAtFirst, iAtLast
226      iSp1 = species(iAtom1)
227      nOrb1 = orb%nOrbSpecies(iSp1)
228      do iNeigh = 1, nNeighbourSK(iAtom1)
229        iAtom2 = iNeighbour(iNeigh, iAtom1)
230        iAtom2f = img2CentCell(iAtom2)
231        iSp2 = species(iAtom2f)
232        if (iAtom1 /= iAtom2f) then
233          nOrb2 = orb%nOrbSpecies(iSp2)
234          iOrig = iPair(iNeigh,iAtom1) + 1
235          sqrDMTmp(1:nOrb2,1:nOrb1) = reshape(DM(iOrig:iOrig+nOrb1*nOrb2-1,1),(/nOrb2,nOrb1/))
236          sqrEDMTmp(1:nOrb2,1:nOrb1) = reshape(EDM(iOrig:iOrig+nOrb1*nOrb2-1),(/nOrb2,nOrb1/))
237          call derivator%getFirstDeriv(hPrimeTmp, skHamCont, coords, species, iAtom1, iAtom2, orb)
238          call derivator%getFirstDeriv(sPrimeTmp, skOverCont, coords, species, iAtom1, iAtom2, orb)
239
240          derivTmp(:) = 0.0_dp
241          ! note factor of 2 for implicit summation over lower triangle of density matrix:
242          do ii = 1, 3
243            derivTmp(ii) = 2.0_dp * (&
244                & sum(sqrDMTmp(1:nOrb2,1:nOrb1)*hPrimeTmp(1:nOrb2,1:nOrb1,ii))&
245                & - sum(sqrEDMTmp(1:nOrb2,1:nOrb1)*sPrimeTmp(1:nOrb2,1:nOrb1,ii)))
246          end do
247
248          do iSpin = 1, nSpin
249            do ii = 1, 3
250              shiftSprime(1:nOrb2,1:nOrb1) = 0.5_dp * (&
251                  & matmul(sPrimeTmp(1:nOrb2,1:nOrb1,ii), shift(1:nOrb1,1:nOrb1,iAtom1,iSpin) )&
252                  & + matmul(shift(1:nOrb2,1:nOrb2,iAtom2f,iSpin), sPrimeTmp(1:nOrb2,1:nOrb1,ii)) )
253              ! again factor of 2 from lower triangle, cf published force expressions for SCC:
254              derivTmp(ii) = derivTmp(ii) + 2.0_dp * ( sum(shiftSprime(1:nOrb2,1:nOrb1) *&
255                  & reshape(DM(iOrig:iOrig+nOrb1*nOrb2-1,iSpin),(/nOrb2,nOrb1/)) ) )
256            end do
257          end do
258
259          ! forces from atom 1 on atom 2f and 2f onto 1
260          deriv(:,iAtom1) = deriv(:,iAtom1) + derivTmp(:)
261          deriv(:,iAtom2f) = deriv(:,iAtom2f) - derivTmp(:)
262
263        end if
264      enddo
265    enddo
266    !$OMP END PARALLEL DO
267
268    call assembleChunks(env, deriv)
269
270  end subroutine derivative_block
271
272
273  !> The SCC and spin electronic force contribution for all atoms, including complex contributions,
274  !> for example from spin-orbit
275  subroutine derivative_iBlock(env, deriv, derivator, DM, iDM, EDM, skHamCont, skOverCont, coords,&
276      & species, iNeighbour, nNeighbourSK, img2CentCell, iPair, orb, shift, iShift)
277
278    !> Computational environment settings
279    type(TEnvironment), intent(in) :: env
280
281    !> x,y,z derivatives for each real atom in the system
282    real(dp), intent(out) :: deriv(:,:)
283
284    !> Differentiatior for the non-scc components
285    class(NonSccDiff), intent(in) :: derivator
286
287    !> density matrix in packed format
288    real(dp), intent(in) :: DM(:,:)
289
290    !> imaginary part of the density matrix in packed format
291    real(dp), intent(in) :: iDM(:,:)
292
293    !> energy-weighted density matrix in packed format
294    real(dp), intent(in) :: EDM(:)
295
296    !> Container for SK Hamiltonian integrals
297    type(OSlakoCont) :: skHamCont
298
299    !> Container for SK overlap integrals
300    type(OSlakoCont) :: skOverCont
301
302    !> list of all atomic coordinates
303    real(dp), intent(in) :: coords(:,:)
304
305    !> list of all atomic species
306    integer, intent(in) :: species(:)
307
308    !> neighbour list for atoms
309    integer, intent(in) :: iNeighbour(0:,:)
310
311    !> number of neighbours of each atom
312    integer, intent(in) :: nNeighbourSK(:)
313
314    !> indexing array for periodic image atoms
315    integer, intent(in) :: img2CentCell(:)
316
317    !> indexing array for the Hamiltonian
318    integer, intent(in) :: iPair(0:,:)
319
320    !> Information about the shells and orbitals in the system.
321    type(TOrbitals), intent(in) :: orb
322
323    !> block shift from the potential
324    real(dp), intent(in) :: shift(:,:,:,:)
325
326    !> imaginary block shift from the potential
327    real(dp), intent(in) :: iShift(:,:,:,:)
328
329    integer :: iOrig, iSpin, ii, nSpin, nAtom
330    integer :: iNeigh, iAtom1, iAtom2, iAtom2f, iSp1, iSp2
331    integer :: nOrb1, nOrb2
332
333    real(dp) :: sqrDMTmp(orb%mOrb,orb%mOrb)
334    real(dp) :: sqrEDMTmp(orb%mOrb,orb%mOrb)
335    complex(dp) :: shiftSprime(orb%mOrb,orb%mOrb)
336    real(dp) :: hPrimeTmp(orb%mOrb,orb%mOrb,3),sPrimeTmp(orb%mOrb,orb%mOrb,3)
337    real(dp) :: derivTmp(3)
338    complex(dp), parameter :: i = (0.0_dp,1.0_dp)
339    integer :: iAtFirst, iAtLast
340
341    nAtom = size(orb%nOrbAtom)
342    nSpin = size(shift,dim=4)
343    @:ASSERT(nSpin == 1 .or. nSpin == 2 .or. nSpin ==4)
344    @:ASSERT(size(deriv,dim=1) == 3)
345    @:ASSERT(size(deriv,dim=2)==nAtom)
346    @:ASSERT(size(DM,dim=1)==size(EDM,dim=1))
347    @:ASSERT(size(DM,dim=2)==nSpin)
348    @:ASSERT(all(shape(iDM)==shape(DM)))
349    @:ASSERT(size(shift,dim=1)==orb%mOrb)
350    @:ASSERT(size(shift,dim=2)==orb%mOrb)
351    @:ASSERT(size(shift,dim=3)==nAtom)
352    @:ASSERT(all(shape(iShift)==shape(shift)))
353
354    deriv(:,:) = 0.0_dp
355
356    call distributeRangeInChunks(env, 1, nAtom, iAtFirst, iAtLast)
357
358    !$OMP PARALLEL DO PRIVATE(iAtom1,iSp1,nOrb1,iNeigh,iAtom2,iAtom2f,iSp2,nOrb2,iOrig,sqrDMTmp, &
359    !$OMP& sqrEDMTmp,hPrimeTmp,sPrimeTmp,derivTmp,shiftSprime,iSpin,ii) DEFAULT(SHARED) &
360    !$OMP& SCHEDULE(RUNTIME) REDUCTION(+:deriv)
361    do iAtom1 = iAtFirst, iAtLast
362      iSp1 = species(iAtom1)
363      nOrb1 = orb%nOrbSpecies(iSp1)
364      do iNeigh = 1, nNeighbourSK(iAtom1)
365        iAtom2 = iNeighbour(iNeigh, iAtom1)
366        iAtom2f = img2CentCell(iAtom2)
367        iSp2 = species(iAtom2f)
368        if (iAtom1 /= iAtom2f) then
369          nOrb2 = orb%nOrbSpecies(iSp2)
370          iOrig = iPair(iNeigh,iAtom1) + 1
371          sqrDMTmp(1:nOrb2,1:nOrb1) = reshape(DM(iOrig:iOrig+nOrb1*nOrb2-1,1),(/nOrb2,nOrb1/))
372          sqrEDMTmp(1:nOrb2,1:nOrb1) = reshape(EDM(iOrig:iOrig+nOrb1*nOrb2-1),(/nOrb2,nOrb1/))
373          call derivator%getFirstDeriv(hPrimeTmp, skHamCont, coords, species, iAtom1, iAtom2, orb)
374          call derivator%getFirstDeriv(sPrimeTmp, skOverCont, coords, species, iAtom1, iAtom2, orb)
375
376          derivTmp(:) = 0.0_dp
377          ! note factor of 2 for implicit summation over lower triangle of density matrix:
378          do ii = 1, 3
379            derivTmp(ii) = 2.0_dp * (&
380                & sum(sqrDMTmp(1:nOrb2,1:nOrb1)*hPrimeTmp(1:nOrb2,1:nOrb1,ii))&
381                & - sum(sqrEDMTmp(1:nOrb2,1:nOrb1)*sPrimeTmp(1:nOrb2,1:nOrb1,ii)))
382          end do
383
384          do iSpin = 1, nSpin
385            do ii = 1, 3
386              shiftSprime(1:nOrb2,1:nOrb1) = 0.5_dp * (&
387                  & matmul(sPrimeTmp(1:nOrb2,1:nOrb1,ii), shift(1:nOrb1,1:nOrb1,iAtom1,iSpin) )&
388                  & + matmul(shift(1:nOrb2,1:nOrb2,iAtom2f,iSpin), sPrimeTmp(1:nOrb2,1:nOrb1,ii)) )
389              ! again factor of 2 from lower triangle sum of DM
390              derivTmp(ii) = derivTmp(ii)&
391                  & + 2.0_dp* ( real(sum(shiftSprime(1:nOrb2,1:nOrb1)&
392                  & * reshape(DM(iOrig:iOrig+nOrb1*nOrb2-1,iSpin), (/nOrb2,nOrb1/)))) )
393            end do
394          end do
395
396          do iSpin = 1, nSpin
397            do ii = 1, 3
398              shiftSprime(1:nOrb2,1:nOrb1) = 0.5_dp *  (&
399                  & matmul(sPrimeTmp(1:nOrb2,1:nOrb1,ii), ishift(1:nOrb1,1:nOrb1,iAtom1,iSpin) )&
400                  & + matmul(ishift(1:nOrb2,1:nOrb2,iAtom2f,iSpin), sPrimeTmp(1:nOrb2,1:nOrb1,ii)) )
401              ! again factor of 2 from lower triangle sum of DM
402              derivTmp(ii) = derivTmp(ii)&
403                  & + 2.0_dp * real(sum(shiftSprime(1:nOrb2,1:nOrb1) *&
404                  & reshape(iDM(iOrig:iOrig+nOrb1*nOrb2-1,iSpin), (/nOrb2,nOrb1/))))
405            end do
406          end do
407
408          deriv(:,iAtom1) = deriv(:,iAtom1) + derivTmp(:)
409          deriv(:,iAtom2f) = deriv(:,iAtom2f) - derivTmp(:)
410
411        end if
412      enddo
413    enddo
414    !$OMP END PARALLEL DO
415
416    call assembleChunks(env, deriv)
417
418  end subroutine derivative_iBlock
419
420end module dftbp_forces
421