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!> Helper routines for transition charges between levels.
11module dftbp_transcharges
12  use dftbp_assert
13  use dftbp_accuracy
14  implicit none
15  private
16
17  public :: TTransCharges, TTransCharges_init
18  public :: transq
19
20  !> Internal data of the transition charges
21  type :: TTransCharges
22    private
23
24    !> should transition charges be cached in memory or evaluated when needed?
25    logical :: tCacheCharges
26
27    !> storage if caching the occupied -> virtual transition charges
28    real(dp), allocatable :: qCacheOccVirt(:,:)
29
30    !> Number of transitions in the cache
31    integer :: nTransitions
32
33    !> Number of atoms in the system
34    integer :: nAtom
35
36    !> Number of transitions within the spin-up block
37    integer :: nMatUp
38
39  contains
40
41    procedure :: qTransIJ =>TTransCharges_qTransIJ
42    procedure :: qMatVec => TTransCharges_qMatVec
43    procedure :: qVecMat => TTransCharges_qVecMat
44
45  end type TTransCharges
46
47
48contains
49
50  !> initialise the cache/on-the fly transition charge evaluator
51  subroutine TTransCharges_init(this, iAtomStart, sTimesGrndEigVecs, grndEigVecs, nTrans, nMatUp,&
52      & getij, win, tStore)
53
54    !> Instance
55    type(TTransCharges), intent(out) :: this
56
57    !> Starting position of each atom in the list of orbitals
58    integer, intent(in) :: iAtomStart(:)
59
60    !> Overlap times eigenvector: sum_m Smn cmi (nOrb, nOrb)
61    real(dp), intent(in) :: sTimesGrndEigVecs(:,:,:)
62
63    !> Eigenvectors (nOrb, nOrb)
64    real(dp), intent(in) :: grndEigVecs(:,:,:)
65
66    !> number of transitions in the system
67    integer, intent(in) :: nTrans
68
69    !> number of up-up excitations
70    integer, intent(in) :: nMatUp
71
72    !> should transitions be stored?
73    logical, intent(in) :: tStore
74
75    !> index array for for single particle excitations
76    integer, intent(in) :: getij(:,:)
77
78    !> index array for single particle excitions that are included
79    integer, intent(in) :: win(:)
80
81    integer :: ij, ii, jj, kk
82    logical :: updwn
83
84    this%nTransitions = nTrans
85    this%nAtom = size(iAtomStart) - 1
86    this%nMatUp = nMatUp
87
88    if (tStore) then
89
90      @:ASSERT(.not.allocated(this%qCacheOccVirt))
91      allocate(this%qCacheOccVirt(this%nAtom, nTrans))
92
93      !!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ij,ii,jj,kk,updwn) SCHEDULE(RUNTIME)
94      do ij = 1, nTrans
95        kk = win(ij)
96        ii = getij(kk,1)
97        jj = getij(kk,2)
98        !call indxov(win, ij, getij, ii, jj)
99        updwn = (kk <= this%nMatUp)
100        this%qCacheOccVirt(:,ij) = transq(ii, jj, iAtomStart, updwn,  sTimesGrndEigVecs,&
101            & grndEigVecs)
102      end do
103      !!$OMP  END PARALLEL DO
104
105      this%tCacheCharges = .true.
106
107    else
108
109      this%tCacheCharges = .false.
110
111    end if
112
113  end subroutine TTransCharges_init
114
115
116  !> returns transtion charges between single particle levels
117  pure function TTransCharges_qTransIJ(this, ij, iAtomStart, sTimesGrndEigVecs, grndEigVecs, getij,&
118      & win) result(q)
119
120    !> instance of the transition charge object
121    class(TTransCharges), intent(in) :: this
122
123    !> Index of transition
124    integer, intent(in) :: ij
125
126    !> Starting position of each atom in the list of orbitals
127    integer, intent(in) :: iAtomStart(:)
128
129    !> Overlap times eigenvector: sum_m Smn cmi (nOrb, nOrb)
130    real(dp), intent(in) :: sTimesGrndEigVecs(:,:,:)
131
132    !> Eigenvectors (nOrb, nOrb)
133    real(dp), intent(in) :: grndEigVecs(:,:,:)
134
135    !> index array for for single particle excitations
136    integer, intent(in) :: getij(:,:)
137
138    !> index array for single particle excitions that are included
139    integer, intent(in) :: win(:)
140
141    !> Transition charge on exit. (nAtom)
142    real(dp), dimension(size(iAtomStart)-1) :: q
143
144    logical :: updwn
145    integer :: ii, jj, kk
146
147    if (allocated(this%qCacheOccVirt)) then
148      q(:) = this%qCacheOccVirt(:, ij)
149    else
150      kk = win(ij)
151      ii = getij(kk,1)
152      jj = getij(kk,2)
153      updwn = (kk <= this%nMatUp)
154      q(:) = transq(ii, jj, iAtomStart, updwn, sTimesgrndEigVecs, grndEigVecs)
155    end if
156
157  end function TTransCharges_qTransIJ
158
159
160  !> Transition charges left producted with a vector Q * v
161  pure subroutine TTransCharges_qMatVec(this, iAtomStart, sTimesGrndEigVecs, grndEigVecs, getij,&
162      & win, vector, qProduct)
163
164    !> instance of the transition charge object
165    class(TTransCharges), intent(in) :: this
166
167    !> Starting position of each atom in the list of orbitals
168    integer, intent(in) :: iAtomStart(:)
169
170    !> Overlap times eigenvector: sum_m Smn cmi (nOrb, nOrb)
171    real(dp), intent(in) :: sTimesGrndEigVecs(:,:,:)
172
173    !> Eigenvectors (nOrb, nOrb)
174    real(dp), intent(in) :: grndEigVecs(:,:,:)
175
176    !> index array for for single particle excitations
177    integer, intent(in) :: getij(:,:)
178
179    !> index array for single particle excitions that are included
180    integer, intent(in) :: win(:)
181
182    !> vector to product with the transition charges
183    real(dp), intent(in) :: vector(:)
184
185    !> Product on exit
186    real(dp), intent(inout) :: qProduct(:)
187
188    real(dp), allocatable :: qij(:)
189    integer :: ii, jj, ij, kk
190    logical :: updwn
191
192    if (this%tCacheCharges) then
193
194      qProduct(:) = qProduct + matmul(this%qCacheOccVirt, vector)
195
196    else
197
198      allocate(qij(this%nAtom))
199
200      !!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ij,ii,jj,kk,updwn,qij)&
201      !!$OMP& SCHEDULE(RUNTIME) REDUCTION(+:qProduct)
202      do ij = 1, this%nTransitions
203        kk = win(ij)
204        ii = getij(kk,1)
205        jj = getij(kk,2)
206        updwn = (kk <= this%nMatUp)
207        qij(:) = transq(ii, jj, iAtomStart, updwn, sTimesGrndEigVecs, grndEigVecs)
208        qProduct(:) = qProduct + qij * vector(ij)
209      end do
210      !!$OMP  END PARALLEL DO
211
212      deallocate(qij)
213
214    end if
215
216  end subroutine TTransCharges_qMatVec
217
218
219  !> Transition charges right producted with a vector v * Q
220  pure subroutine TTransCharges_qVecMat(this, iAtomStart, sTimesGrndEigVecs, grndEigVecs, getij,&
221      & win, vector, qProduct)
222
223    !> instance of the transition charge object
224    class(TTransCharges), intent(in) :: this
225
226    !> Starting position of each atom in the list of orbitals
227    integer, intent(in) :: iAtomStart(:)
228
229    !> Overlap times eigenvector: sum_m Smn cmi (nOrb, nOrb)
230    real(dp), intent(in) :: sTimesGrndEigVecs(:,:,:)
231
232    !> Eigenvectors (nOrb, nOrb)
233    real(dp), intent(in) :: grndEigVecs(:,:,:)
234
235    !> index array for for single particle excitations
236    integer, intent(in) :: getij(:,:)
237
238    !> index array for single particle excitions that are included
239    integer, intent(in) :: win(:)
240
241    !> vector to product with the transition charges
242    real(dp), intent(in) :: vector(:)
243
244    !> Product on exit
245    real(dp), intent(inout) :: qProduct(:)
246
247    real(dp), allocatable :: qij(:)
248    integer :: ii, jj, ij, kk
249    logical :: updwn
250
251    if (this%tCacheCharges) then
252
253      qProduct(:) = qProduct + matmul(vector, this%qCacheOccVirt)
254
255    else
256
257      allocate(qij(this%nAtom))
258
259      !!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ij,ii,jj,kk,updwn,qij)&
260      !!$OMP& SCHEDULE(RUNTIME)
261      do ij = 1, this%nTransitions
262        kk = win(ij)
263        ii = getij(kk,1)
264        jj = getij(kk,2)
265        updwn = (kk <= this%nMatUp)
266        qij(:) = transq(ii, jj, iAtomStart, updwn, sTimesGrndEigVecs, grndEigVecs)
267        qProduct(ij) = qProduct(ij) + dot_product(qij, vector)
268      end do
269      !!$OMP  END PARALLEL DO
270
271      deallocate(qij)
272
273    end if
274
275  end subroutine TTransCharges_qVecMat
276
277
278  !> Calculates atomic transition charges for a specified excitation.
279  !> Calculates qij = 0.5 * (c_i S c_j + c_j S c_i) where c_i and c_j are selected eigenvectors, and
280  !> S the overlap matrix.
281  !> Since qij is atomic quantity (so far) the corresponding values for the atom are summed up.
282  !> Note: the parameters 'updwn' were added for spin alpha and beta channels.
283  pure function transq(ii, jj, iAtomStart, updwn, sTimesGrndEigVecs, grndEigVecs) result(qij)
284
285    !> Index of inital state.
286    integer, intent(in) :: ii
287
288    !> Index of final state.
289    integer, intent(in) :: jj
290
291    !> Starting position of each atom in the list of orbitals.
292    integer, intent(in) :: iAtomStart(:)
293
294    !> up spin channel (T) or down spin channel (F)
295    logical, intent(in) :: updwn
296
297    !> Overlap times eigenvector: sum_m Smn cmi (nOrb, nOrb).
298    real(dp), intent(in) :: sTimesGrndEigVecs(:,:,:)
299
300    !> Eigenvectors (nOrb, nOrb)
301    real(dp), intent(in) :: grndEigVecs(:,:,:)
302
303    !> Transition charge on exit. (nAtom)
304    real(dp) :: qij(size(iAtomStart)-1)
305
306    integer :: kk, aa, bb, ss
307    real(dp) :: qTmp(size(grndEigVecs,dim=1))
308
309    ss = 1
310    if (.not. updwn) then
311      ss = 2
312    end if
313
314    qTmp(:) =  grndEigVecs(:,ii,ss) * sTimesGrndEigVecs(:,jj,ss)&
315        & + grndEigVecs(:,jj,ss) * sTimesGrndEigVecs(:,ii,ss)
316    do kk = 1, size(qij)
317      aa = iAtomStart(kk)
318      bb = iAtomStart(kk + 1) -1
319      qij(kk) = 0.5_dp * sum(qTmp(aa:bb))
320    end do
321
322  end function transq
323
324
325end module dftbp_transcharges
326