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