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!> Contains routines to manipulate orbital equivalency relations.
11!>
12!> An orbital equivalency relation is a mapping, which maps the orbitals of the atoms onto a one
13!> dimensional vector, where equivalent orbitals are mapped on the same element in the vector. Two
14!> orbitals are equivalent, if charge can be transferred between the orbitals, without changing the
15!> resulting Hamiltonian or the total energy. The mapping is an (mmAng, nAtom, nSpin) shaped integer
16!> array, where the integer for (iOrb, iAtom, iSpin) specifies the position in the 1D vector for
17!> orbital iOrb on atom iAtom for spin iSpin. Values must be positive integers and continuous. Zeros
18!> in the mapping vector stand for non-existent orbitals.
19module dftbp_orbitalequiv
20  use dftbp_assert
21  use dftbp_accuracy
22  use dftbp_commontypes
23  implicit none
24  private
25
26  public :: OrbitalEquiv_merge, OrbitalEquiv_reduce, OrbitalEquiv_expand
27
28contains
29
30
31  !> Merges two equivalency arrays by finding the intersection in the equivalences.
32  !> Note: Equivalences marked below 0 are mapped to 0 in the final equivalence matrix
33  subroutine OrbitalEquiv_merge(equiv1, equiv2, orb, equivNew)
34
35    !> First equivalency array.
36    integer, intent(in) :: equiv1(:,:,:)
37
38    !> Second equivalency array
39    integer, intent(in) :: equiv2(:,:,:)
40
41    !> Contains information about the atomic orbitals in the system.
42    type(TOrbitals), intent(in) :: orb
43
44    !> New equivalency array on return.
45    integer, intent(out) :: equivNew(:,:,:)
46
47    integer :: nAtom, nSpin
48    integer :: newInd, iS, iAt, iOrb
49    logical, allocatable :: mask(:,:,:), tmpMask(:,:,:)
50
51    nAtom = size(equiv1, dim=2)
52    nSpin = size(equiv1, dim=3)
53
54    @:ASSERT(all(shape(equiv1) == shape(equiv2)))
55    @:ASSERT(all(shape(equiv1) == shape(equivNew)))
56    @:ASSERT(size(equiv1, dim=1) >= orb%mOrb)
57    @:ASSERT(all((equiv1 /= 0) .eqv. (equiv2 /=0)))
58
59    allocate(mask(size(equiv1, dim=1), size(equiv1, dim=2), size(equiv1, dim=3)))
60    allocate(tmpMask(size(equiv1, dim=1), size(equiv1, dim=2), size(equiv1, dim=3)))
61
62    ! True for the elements to be processed
63    mask(:,:,:) = (equiv1 > 0 .and. equiv2 > 0)
64    equivNew(:,:,:) = 0
65    ! Position in the reduced vector
66    newInd = 1
67    do iS = 1, nSpin
68      do iAt = 1, nAtom
69        do iOrb = 1, orb%nOrbAtom(iAt)
70          ! If element had been already processed, skip
71          if (.not. mask(iOrb, iAt, iS)) then
72            cycle
73          end if
74          ! Get mask for orbitals equivalent to (iOrb, iAt, iS) in both arrays
75          tmpMask = (equiv1 == equiv1(iOrb, iAt, iS) &
76              &.and. equiv2 == equiv2(iOrb, iAt, iS))
77          ! If equivalent orbitals exist, map them to reduced vector
78          if (any(tmpMask)) then
79            where (tmpMask)
80              equivNew = newInd
81            end where
82            newInd = newInd + 1
83            mask = mask .and. .not. tmpMask
84          end if
85        end do
86      end do
87    end do
88
89  end subroutine OrbitalEquiv_merge
90
91
92  !> Reduces passed orbital property by summing up on equivalent orbitals.
93  !> Note: equivalences of 0 or below are not reduced
94  subroutine OrbitalEquiv_reduce(input, equiv, orb, output)
95
96    !> The vector containing the values for all orbitals
97    real(dp), intent(in) :: input(:,:,:)
98
99    !> Vector containing equivalency relations between the orbitals.
100    integer, intent(in) :: equiv(:,:,:)
101
102    !> Contains information about the atomic orbitals in the system
103    type(TOrbitals), intent(in) :: orb
104
105    !> Contains the equivalent orbital sums at exit as vector.
106    real(dp), intent(out) :: output(:)
107
108    integer :: nAtom, nSpin
109    integer :: iS, iOrb, iAt
110
111    nAtom = size(input, dim=2)
112    nSpin = size(input, dim=3)
113
114    @:ASSERT(size(input, dim=1) == orb%mOrb)
115    @:ASSERT(all(shape(equiv) == (/ orb%mOrb, nAtom, nSpin /)))
116    @:ASSERT(size(output) == maxval(equiv))
117
118    output(:) = 0.0_dp
119    do iS = 1, nSpin
120      do iAt = 1, nAtom
121        do iOrb = 1, orb%nOrbAtom(iAt)
122          if (equiv(iOrb, iAt, iS) > 0) then
123            output(equiv(iOrb, iAt, iS)) = output(equiv(iOrb, iAt, iS)) &
124                & + input(iOrb, iAt, iS)
125          end if
126        end do
127      end do
128    end do
129
130  end subroutine OrbitalEquiv_reduce
131
132
133  !> Expands a reduced vector by putting every element of it into the first corresponding orbital
134  !> and putting zero for all other equivalent orbitals.
135  !> Note: equivalences of 0 or below are not expanded
136  subroutine OrbitalEquiv_expand(input, equiv, orb, output)
137
138    !> Reduced vector.
139    real(dp), intent(in) :: input(:)
140
141    !> Equivalency relation.
142    integer, intent(in) :: equiv(:,:,:)
143
144    !> Contains information about the atomic orbitals in the system
145    type(TOrbitals), intent(in) :: orb
146
147    !> Expanded array.
148    real(dp), intent(out) :: output(:,:,:)
149
150    integer :: nSpin, nAtom
151    integer ::iS, iAt, iOrb
152    logical, allocatable :: mask(:)
153
154    nSpin = size(output, dim=3)
155    nAtom = size(output, dim=2)
156
157    @:ASSERT(all(shape(equiv) == shape(output)))
158    @:ASSERT(maxval(equiv) == size(input))
159
160    allocate(mask(0:size(input)))
161
162    mask(:) = .true.
163    output(:,:,:) = 0.0_dp
164    do iS = 1, nSpin
165      do iAt = 1, nAtom
166        do iOrb = 1, orb%nOrbAtom(iAt)
167          if (mask(equiv(iOrb, iAt, iS))) then
168            if (equiv(iOrb, iAt, iS) > 0) then
169              output(iOrb, iAt, iS) = input(equiv(iOrb, iAt, iS))
170              mask(equiv(iOrb, iAt, iS)) = .false.
171            end if
172          end if
173        end do
174      end do
175    end do
176
177  end subroutine OrbitalEquiv_expand
178
179end module dftbp_orbitalequiv
180