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