Add mDFTB

src/dftbp/dftb/CMakeLists.txt

@@ -34,6 +34,7 @@ set(sources-fpp
   ${curdir}/halogenx.F90
   ${curdir}/hamiltonian.F90
   ${curdir}/hybridxc.F90
+  ${curdir}/multiexpan.F90
   ${curdir}/nonscc.F90
   ${curdir}/onscorrection.F90
   ${curdir}/orbitalequiv.F90

src/dftbp/dftb/energytypes.F90

@@ -35,6 +35,14 @@ module dftbp_dftb_energytypes
     !> Range-separation energy
     real(dp) :: Efock = 0.0_dp
 
+    !> MultiPole energy
+    real(dp) :: EDftbMultiExpan = 0.0_dp
+    real(dp) :: EDftbMultiExpanMD = 0.0_dp
+    real(dp) :: EDftbMultiExpanDD = 0.0_dp
+    real(dp) :: EDftbMultiExpanMQ = 0.0_dp
+    real(dp) :: EDftbMultiExpanDQ = 0.0_dp
+    real(dp) :: EDftbMultiExpanQQ = 0.0_dp
+
     !> Spin orbit energy
     real(dp) :: ELS = 0.0_dp
 
@@ -121,6 +129,9 @@ module dftbp_dftb_energytypes
     !> Atom resolved spin
     real(dp), allocatable :: atomSpin(:)
 
+    !> Atom resolved Multipole
+    real(dp), allocatable :: atomDftbMultiExpan(:)
+
     !> Atom resolved spin orbit
     real(dp), allocatable :: atomLS(:)
 
@@ -192,6 +203,7 @@ subroutine TEnergies_init(this, nAtom, nSpin)
     allocate(this%atomOnSite(nAtom))
     allocate(this%atomHalogenX(nAtom))
     allocate(this%atom3rd(nAtom))
+    allocate(this%atomDftbMultiExpan(nAtom))
     allocate(this%atomSolv(nAtom))
     allocate(this%atomTotal(nAtom))
     this%atomRep(:) = 0.0_dp
@@ -214,6 +226,12 @@ subroutine TEnergies_init(this, nAtom, nSpin)
     this%ESCC = 0.0_dp
     this%Espin = 0.0_dp
     this%Efock = 0.0_dp
+    this%EDftbMultiExpan = 0.0_dp
+    this%EDftbMultiExpanMD = 0.0_dp
+    this%EDftbMultiExpanDD = 0.0_dp
+    this%EDftbMultiExpanMQ = 0.0_dp
+    this%EDftbMultiExpanDQ = 0.0_dp
+    this%EDftbMultiExpanQQ = 0.0_dp
     this%ELS = 0.0_dp
     this%Edftbu = 0.0_dp
     this%Eext = 0.0_dp

src/dftbp/dftb/getenergies.F90

@@ -18,6 +18,7 @@ module dftbp_dftb_getenergies
   use dftbp_dftb_dftbplusu, only : TDftbU
   use dftbp_dftb_dispiface, only : TDispersionIface
   use dftbp_dftb_energytypes, only : TEnergies
+  use dftbp_dftb_multiexpan, only : TDftbMultiExpan
   use dftbp_dftb_onsitecorrection, only : getEons
   use dftbp_dftb_periodic, only : TNeighbourList
   use dftbp_dftb_populations, only : mulliken
@@ -46,8 +47,8 @@ module dftbp_dftb_getenergies
 
 
   !> Calculates various energy contribution that can potentially update for the same geometry
-  subroutine calcEnergies(env, sccCalc, tblite, qOrb, q0, chargePerShell, multipole, species,&
-      & isExtField, isXlbomd, dftbU, tDualSpinOrbit, rhoPrim, H0, orb, neighbourList,&
+  subroutine calcEnergies(env, sccCalc, tblite, qOrb, q0, chargePerShell, multipole, dftbMultiExpan,&
+      & species, isExtField, isXlbomd, dftbU, tDualSpinOrbit, rhoPrim, H0, orb, neighbourList,&
       & nNeighbourSK, img2CentCell, iSparseStart, cellVol, extPressure, TS, potential,&
       & energy, thirdOrd, solvation, hybridXc, reks, qDepExtPot, qBlock, qiBlock, xi,&
       & iAtInCentralRegion, tFixEf, Ef, tRealHS, onSiteElements, errStatus, qNetAtom,&
@@ -74,6 +75,10 @@ subroutine calcEnergies(env, sccCalc, tblite, qOrb, q0, chargePerShell, multipol
     !> Multipole moments
     type(TMultipole), intent(in) :: multipole
 
+    !> DFTB multipole moments
+    type(TDftbMultiExpan), intent(inout), allocatable :: dftbMultiExpan
+
+    !> chemical species
     !> Chemical species
     integer, intent(in) :: species(:)
 
@@ -279,6 +284,13 @@ subroutine calcEnergies(env, sccCalc, tblite, qOrb, q0, chargePerShell, multipol
       energy%ELS = sum(energy%atomLS(iAtInCentralRegion))
     end if
 
+    ! Add contribution for DFTB multipole calculations
+    if (allocated(dftbMultiExpan)) then
+      call dftbMultiExpan%addMultiExpanEnergy(energy%atomDftbMultiExpan, energy%EDftbMultiExpanMD,&
+          & energy%EDftbMultiExpanDD, energy%EDftbMultiExpanMQ, energy%EDftbMultiExpanDQ,&
+          & energy%EDftbMultiExpanQQ, energy%EDftbMultiExpan)
+    end if
+
     ! Add exchange contribution for range separated calculations
     if (allocated(hybridXc) .and. .not. allocated(reks)) then
       if (tRealHS) then
@@ -351,11 +363,12 @@ subroutine sumEnergies(energy)
     type(TEnergies), intent(inout) :: energy
 
     energy%Eelec = energy%EnonSCC + energy%ESCC + energy%Espin + energy%ELS + energy%Edftbu&
-        & + energy%Eext + energy%e3rd + energy%eOnSite + energy%ESolv + energy%Efock
+        & + energy%Eext + energy%e3rd + energy%eOnSite + energy%ESolv + energy%Efock&
+        & + energy%EDftbMultiExpan
 
     energy%atomElec(:) = energy%atomNonSCC + energy%atomSCC + energy%atomSpin + energy%atomDftbu&
-        & + energy%atomLS + energy%atomExt + energy%atom3rd + energy%atomOnSite &
-        & + energy%atomSolv
+        & + energy%atomLS + energy%atomExt + energy%atom3rd + energy%atomOnSite&
+        & + energy%atomSolv + energy%atomDftbMultiExpan
     energy%atomTotal(:) = energy%atomElec + energy%atomRep + energy%atomDisp + energy%atomHalogenX
     energy%Etotal = energy%Eelec + energy%Erep + energy%eDisp + energy%eHalogenX
     energy%EMermin = energy%Etotal - sum(energy%TS)

src/dftbp/dftb/hamiltonian.F90

@@ -17,6 +17,7 @@ module dftbp_dftb_hamiltonian
   use dftbp_dftb_dftbplusu, only : TDftbU
   use dftbp_dftb_dispersions, only : TDispersionIface
   use dftbp_dftb_extfields, only : TEField
+  use dftbp_dftb_multiexpan, only : TDftbMultiExpan
   use dftbp_dftb_periodic, only : TNeighbourList
   use dftbp_dftb_potentials, only : TPotentials
   use dftbp_dftb_scc, only : TScc
@@ -148,7 +149,7 @@ end subroutine resetInternalPotentials
   !> spin, and where relevant dispersion
   subroutine addChargePotentials(env, sccCalc, tblite, updateScc, qInput, q0, chargePerShell,&
       & orb, multipole, species, neighbourList, img2CentCell, spinW, solvation, thirdOrd,&
-      & dispersion, potential, errStatus)
+      & dftbMultiExpan, dispersion, potential, errStatus)
 
     !> Environment settings
     type(TEnvironment), intent(inout) :: env
@@ -195,6 +196,9 @@ subroutine addChargePotentials(env, sccCalc, tblite, updateScc, qInput, q0, char
     !> Third order SCC interactions
     type(TThirdOrder), allocatable, intent(inout) :: thirdOrd
 
+    !> Multipole expansion
+    type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan
+
     !> Potentials acting
     type(TPotentials), intent(inout) :: potential
 
@@ -207,6 +211,7 @@ subroutine addChargePotentials(env, sccCalc, tblite, updateScc, qInput, q0, char
     ! local variables
     real(dp), allocatable :: atomPot(:,:)
     real(dp), allocatable :: shellPot(:,:,:)
+    real(dp), allocatable :: deltaMAtom(:)
     real(dp), allocatable :: dipPot(:,:), quadPot(:,:)
     integer, pointer :: pSpecies0(:)
     integer :: nAtom, nSpin
@@ -264,6 +269,13 @@ subroutine addChargePotentials(env, sccCalc, tblite, updateScc, qInput, q0, char
       end if
     end if
 
+    if (allocated(dftbMultiExpan)) then
+      allocate(deltaMAtom(nAtom))
+      call sccCalc%getDeltaQAtom(deltaMAtom)
+      call dftbMultiExpan%updateDQPotentials(deltaMAtom)
+      deallocate(deltaMAtom)
+    end if
+
     if (allocated(thirdOrd)) then
       call thirdOrd%updateCharges(pSpecies0, neighbourList, qInput, q0, img2CentCell, orb)
       call thirdOrd%getShifts(atomPot(:,1), shellPot(:,:,1))
@@ -331,7 +343,7 @@ end subroutine addBlockChargePotentials
 
   !> Returns the Hamiltonian for the given scc iteration
   subroutine getSccHamiltonian(env, H0, ints, nNeighbourSK, neighbourList, species, orb,&
-      & iSparseStart, img2CentCell, potential, isREKS, ham, iHam)
+      & iSparseStart, img2CentCell, potential, dftbMultiExpan, isREKS, ham, iHam)
 
     !> Environment settings
     type(TEnvironment), intent(in) :: env
@@ -363,6 +375,9 @@ subroutine getSccHamiltonian(env, H0, ints, nNeighbourSK, neighbourList, species
     !> Potential acting on system
     type(TPotentials), intent(in) :: potential
 
+    !> DFTB multipole expansion
+    type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan
+
     !> Is this DFTB/SSR formalism
     logical, intent(in) :: isREKS
 
@@ -413,6 +428,11 @@ subroutine getSccHamiltonian(env, H0, ints, nNeighbourSK, neighbourList, species
           & potential%quadrupoleAtom)
     end if
 
+    if (allocated(dftbMultiExpan)) then
+      call dftbMultiExpan%addMultiExpanHamiltonian(ham, ints%overlap, nNeighbourSK,&
+          & neighbourList%iNeighbour, species, orb, iSparseStart, nAtom, img2CentCell)
+    end if
+
     if (allocated(iHam)) then
       iHam(:,:) = 0.0_dp
       call addShift(env, iHam, ints%overlap, nNeighbourSK, neighbourList%iNeighbour, species, orb,&