From 8b02aa38bfb94e3e13e9b3562158bd019c0e14da Mon Sep 17 00:00:00 2001 From: Van-Quan Vuong Date: Mon, 13 Jan 2025 02:57:40 +0100 Subject: [PATCH] Address code review feedback --- src/dftbp/dftb/CMakeLists.txt | 2 +- src/dftbp/dftb/energytypes.F90 | 59 +- src/dftbp/dftb/getenergies.F90 | 4 +- src/dftbp/dftb/hamiltonian.F90 | 4 +- src/dftbp/dftb/{multiexpan.F90 => mdftb.F90} | 656 +++++++----------- src/dftbp/dftbplus/initprogram.F90 | 99 ++- src/dftbp/dftbplus/inputdata.F90 | 26 +- src/dftbp/dftbplus/main.F90 | 12 +- src/dftbp/dftbplus/mainio.F90 | 69 +- src/dftbp/dftbplus/parser.F90 | 248 ++++--- src/dftbp/timedep/timeprop.F90 | 6 +- test/app/dftb+/mdftb/2H2O/dftb_in.hsd | 12 +- test/app/dftb+/mdftb/32H2O/dftb_in.hsd | 12 +- test/app/dftb+/mdftb/Glycine/dftb_in.hsd | 12 +- .../dftb+/mdftb/Glycine_Ordering/dftb_in.hsd | 12 +- .../dftb+/mdftb/Glycine_RotTran1/dftb_in.hsd | 12 +- .../dftb+/mdftb/Glycine_RotTran2/dftb_in.hsd | 12 +- .../dftb+/mdftb/Glycine_RotTran3/dftb_in.hsd | 12 +- .../dftb+/mdftb/Glycine_mDFTBoff1/dftb_in.hsd | 15 +- .../dftb+/mdftb/Glycine_mDFTBoff2/dftb_in.hsd | 15 +- .../dftb+/mdftb/Glycine_mDFTBoff3/dftb_in.hsd | 15 +- test/app/dftb+/mdftb/NH4+H2O/dftb_in.hsd | 12 +- .../dftb+/mdftb/NH4+H2O_RotTran1/dftb_in.hsd | 12 +- .../dftb+/mdftb/NH4+H2O_RotTran2/dftb_in.hsd | 12 +- .../dftb+/mdftb/NH4+H2O_RotTran3/dftb_in.hsd | 12 +- test/app/dftb+/mdftb/OH-H2O/dftb_in.hsd | 12 +- .../dftb+/mdftb/OH-H2O_RotTran1/dftb_in.hsd | 12 +- .../dftb+/mdftb/OH-H2O_RotTran2/dftb_in.hsd | 12 +- .../dftb+/mdftb/OH-H2O_RotTran3/dftb_in.hsd | 12 +- test/app/dftb+/tests | 4 +- 30 files changed, 618 insertions(+), 796 deletions(-) rename src/dftbp/dftb/{multiexpan.F90 => mdftb.F90} (75%) diff --git a/src/dftbp/dftb/CMakeLists.txt b/src/dftbp/dftb/CMakeLists.txt index b519a76cb4..bdcd85a4e7 100644 --- a/src/dftbp/dftb/CMakeLists.txt +++ b/src/dftbp/dftb/CMakeLists.txt @@ -34,7 +34,7 @@ set(sources-fpp ${curdir}/halogenx.F90 ${curdir}/hamiltonian.F90 ${curdir}/hybridxc.F90 - ${curdir}/multiexpan.F90 + ${curdir}/mdftb.F90 ${curdir}/nonscc.F90 ${curdir}/onscorrection.F90 ${curdir}/orbitalequiv.F90 diff --git a/src/dftbp/dftb/energytypes.F90 b/src/dftbp/dftb/energytypes.F90 index 05132aa254..15236628ff 100644 --- a/src/dftbp/dftb/energytypes.F90 +++ b/src/dftbp/dftb/energytypes.F90 @@ -35,12 +35,22 @@ module dftbp_dftb_energytypes !> Range-separation energy real(dp) :: Efock = 0.0_dp - !> MultiPole energy + !> Total energy due to all multipolar interactions except the monopole-monopole. real(dp) :: EDftbMultiExpan = 0.0_dp + + !> Energy contribution from monopole-dipole interactions. real(dp) :: EDftbMultiExpanMD = 0.0_dp + + !> Energy contribution from dipole-dipole interactions. real(dp) :: EDftbMultiExpanDD = 0.0_dp + + !> Energy contribution from monopole-quadrupole interactions. real(dp) :: EDftbMultiExpanMQ = 0.0_dp + + !> Energy contribution from dipole-quadrupole interactions. real(dp) :: EDftbMultiExpanDQ = 0.0_dp + + !> Energy contribution from quadrupole-quadrupole interactions. real(dp) :: EDftbMultiExpanQQ = 0.0_dp !> Spin orbit energy @@ -129,7 +139,7 @@ module dftbp_dftb_energytypes !> Atom resolved spin real(dp), allocatable :: atomSpin(:) - !> Atom resolved Multipole + !> Atom resolved total multipolar energy real(dp), allocatable :: atomDftbMultiExpan(:) !> Atom resolved spin orbit @@ -191,36 +201,21 @@ subroutine TEnergies_init(this, nAtom, nSpin) this%E0(:) = 0.0_dp this%EBand(:) = 0.0_dp - allocate(this%atomRep(nAtom)) - allocate(this%atomNonSCC(nAtom)) - allocate(this%atomSCC(nAtom)) - allocate(this%atomSpin(nAtom)) - allocate(this%atomLS(nAtom)) - allocate(this%atomDftbu(nAtom)) - allocate(this%atomExt(nAtom)) - allocate(this%atomElec(nAtom)) - allocate(this%atomDisp(nAtom)) - 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 - this%atomNonSCC(:) = 0.0_dp - this%atomSCC(:) = 0.0_dp - this%atomSpin(:) = 0.0_dp - this%atomLS(:) = 0.0_dp - this%atomDftbu(:) = 0.0_dp - this%atomExt(:) = 0.0_dp - this%atomElec(:) = 0.0_dp - this%atomDisp(:) = 0.0_dp - this%atomOnSite(:) = 0.0_dp - this%atomHalogenX(:) = 0.0_dp - this%atom3rd(:) = 0.0_dp - this%atomDftbMultiExpan(:) = 0.0_dp - this%atomSolv(:) = 0.0_dp - this%atomTotal(:) = 0.0_dp + allocate(this%atomRep(nAtom), source=0.0_dp) + allocate(this%atomNonSCC(nAtom), source=0.0_dp) + allocate(this%atomSCC(nAtom), source=0.0_dp) + allocate(this%atomSpin(nAtom), source=0.0_dp) + allocate(this%atomLS(nAtom), source=0.0_dp) + allocate(this%atomDftbu(nAtom), source=0.0_dp) + allocate(this%atomExt(nAtom), source=0.0_dp) + allocate(this%atomElec(nAtom), source=0.0_dp) + allocate(this%atomDisp(nAtom), source=0.0_dp) + allocate(this%atomOnSite(nAtom), source=0.0_dp) + allocate(this%atomHalogenX(nAtom), source=0.0_dp) + allocate(this%atom3rd(nAtom), source=0.0_dp) + allocate(this%atomDftbMultiExpan(nAtom), source=0.0_dp) + allocate(this%atomSolv(nAtom), source=0.0_dp) + allocate(this%atomTotal(nAtom), source=0.0_dp) this%Erep = 0.0_dp this%EnonSCC = 0.0_dp diff --git a/src/dftbp/dftb/getenergies.F90 b/src/dftbp/dftb/getenergies.F90 index f11d9a43f9..5e85f4c556 100644 --- a/src/dftbp/dftb/getenergies.F90 +++ b/src/dftbp/dftb/getenergies.F90 @@ -18,7 +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_multiexpan, only : TMdftb use dftbp_dftb_onsitecorrection, only : getEons use dftbp_dftb_periodic, only : TNeighbourList use dftbp_dftb_populations, only : mulliken @@ -76,7 +76,7 @@ subroutine calcEnergies(env, sccCalc, tblite, qOrb, q0, chargePerShell, multipol type(TMultipole), intent(in) :: multipole !> DFTB multipole moments - type(TDftbMultiExpan), intent(inout), allocatable :: dftbMultiExpan + type(TMdftb), intent(inout), allocatable :: dftbMultiExpan !> chemical species !> Chemical species diff --git a/src/dftbp/dftb/hamiltonian.F90 b/src/dftbp/dftb/hamiltonian.F90 index 814465008f..aafc0671aa 100644 --- a/src/dftbp/dftb/hamiltonian.F90 +++ b/src/dftbp/dftb/hamiltonian.F90 @@ -17,7 +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_multiexpan, only : TMdftb use dftbp_dftb_periodic, only : TNeighbourList use dftbp_dftb_potentials, only : TPotentials use dftbp_dftb_scc, only : TScc @@ -365,7 +365,7 @@ subroutine getSccHamiltonian(env, H0, ints, nNeighbourSK, neighbourList, species type(TPotentials), intent(in) :: potential !> DFTB multipole expansion - type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan + type(TMdftb), allocatable, intent(inout) :: dftbMultiExpan !> Is this DFTB/SSR formalism logical, intent(in) :: isREKS diff --git a/src/dftbp/dftb/multiexpan.F90 b/src/dftbp/dftb/mdftb.F90 similarity index 75% rename from src/dftbp/dftb/multiexpan.F90 rename to src/dftbp/dftb/mdftb.F90 index 768f032921..b01c1440b2 100644 --- a/src/dftbp/dftb/multiexpan.F90 +++ b/src/dftbp/dftb/mdftb.F90 @@ -24,10 +24,36 @@ module dftbp_dftb_multiexpan implicit none private - public :: TDftbMultiExpanInp, TDftbMultiExpan, dftbMultiExpan_init + public :: TMdftbAtomicIntegrals, TMdftbInp, TMdftb, TMdftb_init + + !> Store one-center atomic integrals required by mdftb + type TMdftbAtomicIntegrals + real(dp), allocatable :: DScaling(:) + real(dp), allocatable :: QScaling(:) + real(dp), allocatable :: SXPx(:) + real(dp), allocatable :: PxXDxxyy(:) + real(dp), allocatable :: PxXDzz(:) + real(dp), allocatable :: PyYDxxyy(:) + real(dp), allocatable :: PzZDzz(:) + real(dp), allocatable :: SXXS(:) + real(dp), allocatable :: PxXXPx(:) + real(dp), allocatable :: PyXXPy(:) + real(dp), allocatable :: SXXDxxyy(:) + real(dp), allocatable :: SXXDzz(:) + real(dp), allocatable :: SYYDxxyy(:) + real(dp), allocatable :: SZZDzz(:) + real(dp), allocatable :: DxyXXDxy(:) + real(dp), allocatable :: DyzXXDyz(:) + real(dp), allocatable :: DxxyyXXDzz(:) + real(dp), allocatable :: DzzXXDzz(:) + real(dp), allocatable :: DxxyyYYDzz(:) + real(dp), allocatable :: DzzZZDzz(:) + real(dp), allocatable :: DxzXZDzz(:) + real(dp), allocatable :: DyzYZDxxyy(:) + end type TMdftbAtomicIntegrals !> Input for the MultiExpan module - type TDftbMultiExpanInp + type TMdftbInp !> Orbital information integer :: nOrb, nSpin, nSpecies @@ -39,10 +65,9 @@ module dftbp_dftb_multiexpan !> Species of atoms integer, allocatable :: species(:) - !> One-center atomic intergral data + !> One-center atomic integral data real(dp), allocatable :: atomicDIntgrlScaling(:) real(dp), allocatable :: atomicQIntgrlScaling(:) - real(dp), allocatable :: atomicOnsiteScaling(:) real(dp), allocatable :: atomicSXPxIntgrl(:) real(dp), allocatable :: atomicPxXDxxyyIntgrl(:) real(dp), allocatable :: atomicPxXDzzIntgrl(:) @@ -64,12 +89,11 @@ module dftbp_dftb_multiexpan real(dp), allocatable :: atomicDxzXZDzzIntgrl(:) real(dp), allocatable :: atomicDyzYZDxxyyIntgrl(:) - end type TDftbMultiExpanInp + end type TMdftbInp !> Internal status of TMultiExpan. - type TDftbMultiExpan + type TMdftb integer :: nAtoms, nSpecies, mShells, mShellsReal, nOrb, nSpin - integer, allocatable :: nShells(:) !> Hubbard U values for atoms real(dp), allocatable :: hubbu(:) @@ -84,8 +108,6 @@ module dftbp_dftb_multiexpan real(dp), allocatable :: atomicDIntgrl(:,:,:,:) real(dp), allocatable :: atomicQIntgrl(:,:,:,:,:) - real(dp), allocatable :: atomicOnsiteScaling(:) - !> Mulliken charge per atom real(dp), allocatable :: deltaMAtom(:) !> Dipole charge per atom @@ -131,7 +153,7 @@ module dftbp_dftb_multiexpan procedure :: addAtomicQuadrupoleMoment procedure :: getMultiExpanInfo procedure :: getOrbitalEquiv - end type TDftbMultiExpan + end type TMdftb !> Number of dipole components used in tblite library (x, y, z) integer, parameter :: dimDipole = 3 @@ -145,7 +167,7 @@ module dftbp_dftb_multiexpan subroutine getMultiExpanInfo(this, nDipole, nQuadrupole) !> Data structure - class(TDftbMultiExpan), intent(in) :: this + class(TMdftb), intent(in) :: this !> Number of dipole moment components integer, intent(out) :: nDipole @@ -158,6 +180,7 @@ subroutine getMultiExpanInfo(this, nDipole, nQuadrupole) nDipole = dimDipole nQuadrupole = dimQuadrupole + end subroutine getMultiExpanInfo @@ -165,7 +188,7 @@ end subroutine getMultiExpanInfo subroutine getOrbitalEquiv(this, equivDip, equivQuad) !> Data structure - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> The equivalence vector for cumulative atomic dipole populations integer, intent(out) :: equivDip(:,:) @@ -203,64 +226,51 @@ end subroutine getOrbitalEquiv !> Initializes instance. - subroutine dftbMultiExpan_init(this, inp) + subroutine TMdftb_init(this, inp) !> Instance. - type(TDftbMultiExpan), intent(out) :: this + type(TMdftb), intent(out) :: this !> Input data. - type(TDftbMultiExpanInp), intent(in) :: inp + type(TMdftbInp), intent(in) :: inp - integer, parameter :: xyzLen = 3 integer, parameter :: icx = 1, icy = 2, icz = 3 integer, parameter :: ios = 1, iopy = 2, iopz = 3, iopx = 4 integer, parameter :: iodxy = 5, iodyz = 6, iodzz = 7, iodxz = 8, iodxxyy = 9 real(dp), parameter :: tolZero = 1.0e-15_dp - integer :: nAtoms, maxNOrb, iAt1, iAt2, nOrb1, nOrb2, ii, jj, iSp1, iSp2 + integer :: nAtoms, mOrb, iAt1, iAt2, nOrb1, nOrb2, ii, jj, iSp1, iSp2 integer :: mu, nu, mm, nn real(dp) tmpIntgrl, tmpTrace this%nAtoms = size(inp%orb%nOrbAtom) this%nSpecies = inp%nSpecies - allocate(this%nOrbSpecies(this%nSpecies)) - this%nOrbSpecies(:) = inp%orb%nOrbSpecies(:) + ! allocate(this%nOrbSpecies(this%nSpecies)) + this%nOrbSpecies = inp%orb%nOrbSpecies(:) this%mShells = 1 this%mShellsReal = inp%orb%mShell this%nOrb = inp%nOrb - maxNOrb = maxval(inp%orb%nOrbSpecies(:)) + mOrb = inp%orb%mOrb this%nSpin = inp%nSpin - allocate(this%nShells(this%nSpecies)) - this%nShells(:) = 1 - allocate(this%hubbu(size(inp%hubbu))) - this%hubbu(:) = inp%hubbu + ! allocate(this%hubbu(size(inp%hubbu))) + this%hubbu = inp%hubbu - allocate(this%species(this%nAtoms)) - this%species(:) = inp%species + ! allocate(this%species(this%nAtoms)) + this%species = inp%species - allocate(this%onDQOCharges(3, this%nSpecies)) - this%onDQOCharges(1,:) = 1 - this%onDQOCharges(2,:) = 1 + allocate(this%onDQOCharges(3, this%nSpecies), source=1) this%onDQOCharges(3,:) = 0 - allocate(this%coords(xyzLen, this%nAtoms)) - this%coords(:,:) = 0.0_dp - - allocate(this%atomicDIntgrl(xyzLen, maxNOrb, maxNOrb, this%nSpecies)) - allocate(this%atomicQIntgrl(xyzLen, xyzLen, maxNOrb, maxNOrb, this%nSpecies)) - this%atomicDIntgrl(:,:,:,:) = 0.0_dp - this%atomicQIntgrl(:,:,:,:,:) = 0.0_dp - - allocate(this%atomicOnsiteScaling(this%nSpecies)) - this%atomicOnsiteScaling(:) = 1.0_dp - + allocate(this%coords(3, this%nAtoms), source=0.0_dp) + allocate(this%atomicDIntgrl(3, mOrb, mOrb, this%nSpecies), source=0.0_dp) + allocate(this%atomicQIntgrl(3,3, mOrb, mOrb, this%nSpecies), source=0.0_dp) do iSp1 = 1, this%nSpecies nOrb1 = this%nOrbSpecies(iSp1) - !Dipole + ! Dipole if ( nOrb1 >= 4 ) then - !assign + ! Assign tmpIntgrl = inp%atomicSXPxIntgrl(iSp1) this%atomicDIntgrl(icx,ios,iopx,iSp1) = tmpIntgrl this%atomicDIntgrl(icy,ios,iopy,iSp1) = tmpIntgrl @@ -268,7 +278,7 @@ subroutine dftbMultiExpan_init(this, inp) end if if ( nOrb1 >= 9 ) then - !assign + ! Assign tmpIntgrl = inp%atomicPxXDxxyyIntgrl(iSp1) this%atomicDIntgrl(icx,iopx,iodxxyy,iSp1) = tmpIntgrl this%atomicDIntgrl(icx,iopy,iodxy,iSp1) = tmpIntgrl @@ -278,23 +288,23 @@ subroutine dftbMultiExpan_init(this, inp) this%atomicDIntgrl(icz,iopx,iodxz,iSp1) = tmpIntgrl this%atomicDIntgrl(icz,iopy,iodyz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicPxXDzzIntgrl(iSp1) this%atomicDIntgrl(icx,iopx,iodzz,iSp1) = tmpIntgrl this%atomicDIntgrl(icy,iopy,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicPyYDxxyyIntgrl(iSp1) this%atomicDIntgrl(icy,iopy,iodxxyy,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicPzZDzzIntgrl(iSp1) this%atomicDIntgrl(icz,iopz,iodzz,iSp1) = tmpIntgrl end if - !Quadrupole + ! Quadrupole if ( nOrb1 >= 1 ) then - !assign + ! Assign tmpIntgrl = inp%atomicSXXSIntgrl(iSp1) this%atomicQIntgrl(icx,icx,ios,ios,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icy,ios,ios,iSp1) = tmpIntgrl @@ -302,13 +312,13 @@ subroutine dftbMultiExpan_init(this, inp) end if if ( nOrb1 >= 4 ) then - !assign + ! Assign tmpIntgrl = inp%atomicPxXXPxIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iopx,iopx,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icy,iopy,iopy,iSp1) = tmpIntgrl this%atomicQIntgrl(icz,icz,iopz,iopz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicPyXXPyIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iopy,iopy,iSp1) = tmpIntgrl this%atomicQIntgrl(icx,icx,iopz,iopz,iSp1) = tmpIntgrl @@ -322,27 +332,27 @@ subroutine dftbMultiExpan_init(this, inp) end if if ( nOrb1 >= 9 ) then - !assign + ! Assign tmpIntgrl = inp%atomicSXXDxxyyIntgrl(iSp1) this%atomicQIntgrl(icx,icx,ios,iodxxyy,iSp1) = tmpIntgrl this%atomicQIntgrl(icx,icy,ios,iodxy,iSp1) = tmpIntgrl this%atomicQIntgrl(icx,icz,ios,iodxz,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icz,ios,iodyz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicSXXDzzIntgrl(iSp1) this%atomicQIntgrl(icx,icx,ios,iodzz,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icy,ios,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicSYYDxxyyIntgrl(iSp1) this%atomicQIntgrl(icy,icy,ios,iodxxyy,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicSZZDzzIntgrl(iSp1) this%atomicQIntgrl(icz,icz,ios,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDxyXXDxyIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iodxy,iodxy,iSp1) = tmpIntgrl this%atomicQIntgrl(icx,icx,iodxz,iodxz,iSp1) = tmpIntgrl @@ -353,7 +363,7 @@ subroutine dftbMultiExpan_init(this, inp) this%atomicQIntgrl(icz,icz,iodxz,iodxz,iSp1) = tmpIntgrl this%atomicQIntgrl(icz,icz,iodyz,iodyz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDyzXXDyzIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iodyz,iodyz,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icy,iodxz,iodxz,iSp1) = tmpIntgrl @@ -364,30 +374,30 @@ subroutine dftbMultiExpan_init(this, inp) this%atomicQIntgrl(icx,icz,iodxz,iodxxyy,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icz,iodxy,iodxz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDxxyyXXDzzIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iodxxyy,iodzz,iSp1) = tmpIntgrl this%atomicQIntgrl(icx,icy,iodxy,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDzzXXDzzIntgrl(iSp1) this%atomicQIntgrl(icx,icx,iodzz,iodzz,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icy,iodzz,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDxxyyYYDzzIntgrl(iSp1) this%atomicQIntgrl(icy,icy,iodxxyy,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDzzZZDzzIntgrl(iSp1) this%atomicQIntgrl(icz,icz,iodzz,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDxzXZDzzIntgrl(iSp1) this%atomicQIntgrl(icx,icz,iodxz,iodzz,iSp1) = tmpIntgrl this%atomicQIntgrl(icy,icz,iodyz,iodzz,iSp1) = tmpIntgrl - !assign + ! Assign tmpIntgrl = inp%atomicDyzYZDxxyyIntgrl(iSp1) this%atomicQIntgrl(icy,icz,iodyz,iodxxyy,iSp1) = tmpIntgrl end if @@ -397,11 +407,11 @@ subroutine dftbMultiExpan_init(this, inp) nOrb1 = this%nOrbSpecies(iSp1) do mm = 1, nOrb1 do nn = 1, nOrb1 - do ii = 1, xyzLen + do ii = 1, 3 if (abs(this%atomicDIntgrl(ii,mm,nn,iSp1)) > tolZero ) then this%atomicDIntgrl(ii,nn,mm,iSp1) = this%atomicDIntgrl(ii,mm,nn,iSp1) end if - do jj = 1, xyzLen + do jj = 1, 3 if (abs(this%atomicQIntgrl(ii,jj,mm,nn,iSp1)) > tolZero ) then this%atomicQIntgrl(ii,jj,nn,mm,iSp1) = this%atomicQIntgrl(ii,jj,mm,nn,iSp1) this%atomicQIntgrl(jj,ii,mm,nn,iSp1) = this%atomicQIntgrl(ii,jj,mm,nn,iSp1) @@ -413,7 +423,7 @@ subroutine dftbMultiExpan_init(this, inp) end do end do - !> scaling atomic Intgrls + !> Scaling atomic Intgrls do iSp1 = 1, this%nSpecies this%atomicDIntgrl(:,:,:,iSp1) = inp%atomicDIntgrlScaling(iSp1)& & * this%atomicDIntgrl(:,:,:,iSp1) @@ -428,7 +438,7 @@ subroutine dftbMultiExpan_init(this, inp) do nn = 1, nOrb1 tmpTrace = (this%atomicQIntgrl(1,1,mm,nn,iSp1) + this%atomicQIntgrl(2,2,mm,nn,iSp1)& & + this%atomicQIntgrl(3,3,mm,nn,iSp1)) / 3.0_dp - do ii = 1, xyzLen + do ii = 1, 3 this%atomicQIntgrl(ii,ii,mm,nn,iSp1) = this%atomicQIntgrl(ii,ii,mm,nn,iSp1) - tmpTrace end do end do @@ -444,89 +454,55 @@ subroutine dftbMultiExpan_init(this, inp) end if end do - this%atomicOnsiteScaling(:) = inp%atomicOnsiteScaling(:) - - - allocate(this%deltaMAtom(this%nAtoms)) - allocate(this%deltaDAtom(xyzLen, this%nAtoms)) - allocate(this%deltaQAtom(xyzLen, xyzLen, this%nAtoms)) - allocate(this%f10AB(xyzLen, this%nAtoms, this%nAtoms)) - allocate(this%f20AB(xyzLen, xyzLen, this%nAtoms, this%nAtoms)) - allocate(this%f30AB(xyzLen, xyzLen, xyzLen, this%nAtoms, this%nAtoms)) - allocate(this%f40AB(xyzLen, xyzLen, xyzLen, xyzLen, this%nAtoms, this%nAtoms)) - allocate(this%f50AB(xyzLen, xyzLen, xyzLen, xyzLen, xyzLen, this%nAtoms, this%nAtoms)) - - allocate(this%pot10x1Atom(this%nAtoms)) - allocate(this%pot20x2Atom(this%nAtoms)) - allocate(this%pot10x0Atom(xyzLen, this%nAtoms)) - allocate(this%pot11x1Atom(xyzLen, this%nAtoms)) - allocate(this%pot21x2Atom(xyzLen, this%nAtoms)) - allocate(this%pot20x0Atom(xyzLen, xyzLen, this%nAtoms)) - allocate(this%pot21x1Atom(xyzLen, xyzLen, this%nAtoms)) - allocate(this%pot22x2Atom(xyzLen, xyzLen, this%nAtoms)) - allocate(this%pot30x0Atom(xyzLen, xyzLen, xyzLen, this%nAtoms)) - allocate(this%pot31x1Atom(xyzLen, xyzLen, xyzLen, this%nAtoms)) - allocate(this%pot32x2Atom(xyzLen, xyzLen, xyzLen, this%nAtoms)) - - this%deltaMAtom(:) = 0.0_dp - this%deltaDAtom(:,:) = 0.0_dp - this%deltaQAtom(:,:,:) = 0.0_dp - this%f10AB(:,:,:) = 0.0_dp - this%f20AB(:,:,:,:) = 0.0_dp - this%f30AB(:,:,:,:,:) = 0.0_dp - this%f40AB(:,:,:,:,:,:) = 0.0_dp - this%f50AB(:,:,:,:,:,:,:) = 0.0_dp - - this%pot10x1Atom(:) = 0.0_dp - this%pot20x2Atom(:) = 0.0_dp - this%pot10x0Atom(:,:) = 0.0_dp - this%pot11x1Atom(:,:) = 0.0_dp - this%pot21x2Atom(:,:) = 0.0_dp - this%pot20x0Atom(:,:,:) = 0.0_dp - this%pot21x1Atom(:,:,:) = 0.0_dp - this%pot22x2Atom(:,:,:) = 0.0_dp - this%pot30x0Atom(:,:,:,:) = 0.0_dp - this%pot31x1Atom(:,:,:,:) = 0.0_dp - this%pot32x2Atom(:,:,:,:) = 0.0_dp - - end subroutine dftbMultiExpan_init + allocate(this%deltaMAtom(this%nAtoms), source=0.0_dp) + allocate(this%deltaDAtom(3, this%nAtoms), source=0.0_dp) + allocate(this%deltaQAtom(3,3, this%nAtoms), source=0.0_dp) + allocate(this%f10AB(3, this%nAtoms, this%nAtoms), source=0.0_dp) + allocate(this%f20AB(3,3, this%nAtoms, this%nAtoms), source=0.0_dp) + allocate(this%f30AB(3,3,3, this%nAtoms, this%nAtoms), source=0.0_dp) + allocate(this%f40AB(3,3,3,3, this%nAtoms, this%nAtoms), source=0.0_dp) + allocate(this%f50AB(3,3,3,3,3, this%nAtoms, this%nAtoms), source=0.0_dp) + + allocate(this%pot10x1Atom(this%nAtoms), source=0.0_dp) + allocate(this%pot20x2Atom(this%nAtoms), source=0.0_dp) + allocate(this%pot10x0Atom(3, this%nAtoms), source=0.0_dp) + allocate(this%pot11x1Atom(3, this%nAtoms), source=0.0_dp) + allocate(this%pot21x2Atom(3, this%nAtoms), source=0.0_dp) + allocate(this%pot20x0Atom(3,3, this%nAtoms), source=0.0_dp) + allocate(this%pot21x1Atom(3,3, this%nAtoms), source=0.0_dp) + allocate(this%pot22x2Atom(3,3, this%nAtoms), source=0.0_dp) + allocate(this%pot30x0Atom(3,3,3, this%nAtoms), source=0.0_dp) + allocate(this%pot31x1Atom(3,3,3, this%nAtoms), source=0.0_dp) + allocate(this%pot32x2Atom(3,3,3, this%nAtoms), source=0.0_dp) + + end subroutine TMdftb_init !> Updates data structures if there are changed coordinates for the instance. subroutine updateCoords(this, coords) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> list of atomic coordinates real(dp), intent(in) :: coords(:,:) - integer, parameter :: xyzLen = 3 - integer :: nAtoms, iAt1, iAt2, iSp1, iSp2, nOrb1, nOrb2, ii, jj, ll, mm, nn, kk real(dp) :: rab, u1, u2, coeffTerm1, coeffTerm2, coeffTerm3 real(dp) :: gammaPrime, gammaDoublePrime, gammaTriplePrime, gammaQuadruplePrime,& & gammaQuintuplePrime - real(dp) :: vRabx3(xyzLen), workVx3(xyzLen) - real(dp) :: mI3x3(xyzLen, xyzLen), mRRab3x3(xyzLen, xyzLen), workM3x3(xyzLen, xyzLen),& - & workM3x3x3(xyzLen, xyzLen, xyzLen) - real(dp) :: mRoIab3x3x3(xyzLen, xyzLen, xyzLen), mIotRab3x3x3(xyzLen, xyzLen, xyzLen),& - & mIoRab3x3x3(xyzLen, xyzLen, xyzLen) - real(dp) :: mRIotoRab3x3x3(xyzLen, xyzLen, xyzLen), mRRRab3x3x3(xyzLen, xyzLen, xyzLen) - - real(dp) :: mI3x3otI3x3(xyzLen, xyzLen, xyzLen, xyzLen), mI3x3ohI3x3(xyzLen, xyzLen, xyzLen,& - & xyzLen) - real(dp) :: mI3x3oI3x3(xyzLen, xyzLen, xyzLen, xyzLen), mI3x3otohoI3x3(xyzLen, xyzLen, xyzLen,& - & xyzLen) - real(dp) :: workM3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen), f40Term1M3x3x3x3(xyzLen, xyzLen,& - & xyzLen, xyzLen) - real(dp) :: f40Term3M3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen), f40Term2M3x3x3x3(xyzLen, xyzLen,& - & xyzLen, xyzLen) - real(dp) :: workM3x3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen, xyzLen),& - & f50Term1M3x3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen, xyzLen) - real(dp) :: f50Term3M3x3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen, xyzLen),& - & f50Term2M3x3x3x3x3(xyzLen, xyzLen, xyzLen, xyzLen, xyzLen) + real(dp) :: vRabx3(3), workVx3(3) + real(dp) :: mI3x3(3,3), mRRab3x3(3,3), workM3x3(3,3), workM3x3x3(3,3,3) + real(dp) :: mRoIab3x3x3(3,3,3), mIotRab3x3x3(3,3,3), mIoRab3x3x3(3,3,3) + real(dp) :: mRIotoRab3x3x3(3,3,3), mRRRab3x3x3(3,3,3) + + real(dp) :: mI3x3otI3x3(3,3,3,3), mI3x3ohI3x3(3,3,3,3) + real(dp) :: mI3x3oI3x3(3,3,3,3), mI3x3otohoI3x3(3,3,3,3) + real(dp) :: workM3x3x3x3(3,3,3,3), f40Term1M3x3x3x3(3,3,3,3) + real(dp) :: f40Term3M3x3x3x3(3,3,3,3), f40Term2M3x3x3x3(3,3,3,3) + real(dp) :: workM3x3x3x3x3(3,3,3,3,3), f50Term1M3x3x3x3x3(3,3,3,3,3) + real(dp) :: f50Term3M3x3x3x3x3(3,3,3,3,3), f50Term2M3x3x3x3x3(3,3,3,3,3) this%coords(:,:) = coords this%f10AB(:,:,:) = 0.0_dp @@ -536,7 +512,7 @@ subroutine updateCoords(this, coords) this%f50AB(:,:,:,:,:,:,:) = 0.0_dp mI3x3(:,:) = 0.0_dp - do ii = 1, xyzLen + do ii = 1, 3 mI3x3(ii, ii) = 1.0_dp end do call makeA3x3oB3x3(mI3x3oI3x3, mI3x3, mI3x3) @@ -593,9 +569,9 @@ subroutine updateCoords(this, coords) coeffTerm1 = coeffTerm1 / 6.0_dp coeffTerm2 = coeffTerm2 / 6.0_dp - do ll = 1, xyzLen - do jj = 1, xyzLen - do ii = 1, xyzLen + do ll = 1, 3 + do jj = 1, 3 + do ii = 1, 3 mRoIab3x3x3(ii,jj,ll) = vRabx3(ii) * mI3x3(jj,ll) mIotRab3x3x3(ii,jj,ll) = mI3x3(ii,ll) * vRabx3(jj) mIoRab3x3x3(ii,jj,ll) = mI3x3(ii,jj) * vRabx3(ll) @@ -604,8 +580,7 @@ subroutine updateCoords(this, coords) end do end do mRIotoRab3x3x3(:,:,:) = mRoIab3x3x3(:,:,:) + mIotRab3x3x3(:,:,:) + mIoRab3x3x3(:,:,:) - workM3x3x3(:,:,:) = coeffTerm1 * mRIotoRab3x3x3(:,:,:) & - &+ coeffTerm2 * mRRRab3x3x3(:,:,:) + workM3x3x3(:,:,:) = coeffTerm1 * mRIotoRab3x3x3(:,:,:) + coeffTerm2 * mRRRab3x3x3(:,:,:) this%f30AB(:,:,:,iAt2,iAt1) = workM3x3x3(:,:,:) this%f30AB(:,:,:,iAt1,iAt2) = -workM3x3x3(:,:,:) @@ -630,10 +605,10 @@ subroutine updateCoords(this, coords) f40Term1M3x3x3x3(:,:,:,:) = mI3x3otohoI3x3(:,:,:,:) f40Term2M3x3x3x3(:,:,:,:) = 0.0_dp f40Term3M3x3x3x3(:,:,:,:) = 0.0_dp - do mm = 1, xyzLen - do ll = 1, xyzLen - do jj = 1, xyzLen - do ii = 1, xyzLen + do mm = 1, 3 + do ll = 1, 3 + do jj = 1, 3 + do ii = 1, 3 f40Term2M3x3x3x3(ii,jj,ll,mm) = f40Term2M3x3x3x3(ii,jj,ll,mm) & & + mRRab3x3(ii,jj) * mI3x3(ll,mm) & & + mRRab3x3(ii,ll) * mI3x3(jj,mm) & @@ -669,11 +644,11 @@ subroutine updateCoords(this, coords) f50Term1M3x3x3x3x3(:,:,:,:,:) = 0.0_dp f50Term2M3x3x3x3x3(:,:,:,:,:) = 0.0_dp f50Term3M3x3x3x3x3(:,:,:,:,:) = 0.0_dp - do nn = 1, xyzLen - do mm = 1, xyzLen - do ll = 1, xyzLen - do jj = 1, xyzLen - do ii = 1, xyzLen + do nn = 1, 3 + do mm = 1, 3 + do ll = 1, 3 + do jj = 1, 3 + do ii = 1, 3 f50Term1M3x3x3x3x3(ii,jj,ll,mm,nn) = f50Term1M3x3x3x3x3(ii,jj,ll,mm,nn) & & + mRoIab3x3x3(ii,jj,ll) * mI3x3(mm,nn) & & + mRoIab3x3x3(ii,jj,mm) * mI3x3(ll,nn) & @@ -708,10 +683,10 @@ subroutine updateCoords(this, coords) this%f50AB(:,:,:,:,:,iAt1,iAt2) = -workM3x3x3x3x3(:,:,:,:,:) end do - coeffTerm1 = -this%atomicOnsiteScaling(iSp1) * (3.2_dp * u1)**3 / 96.0_dp + coeffTerm1 = -(3.2_dp * u1)**3 / 96.0_dp this%f20AB(:,:,iAt1,iAt1) = coeffTerm1 * mI3x3(:,:) - coeffTerm1 = this%atomicOnsiteScaling(iSp1) * (3.2_dp * u1)**5 / 5760.0_dp + coeffTerm1 = (3.2_dp * u1)**5 / 5760.0_dp this%f40AB(:,:,:,:,iAt1,iAt1) = coeffTerm1 * mI3x3otohoI3x3(:,:,:,:) end do @@ -724,7 +699,7 @@ subroutine updateDeltaDQAtom(this, over, rho, orb, iNeighbour, nNeighbourSK, img & q0) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> Overlap matrix in packed format real(dp), intent(in) :: over(:) @@ -750,21 +725,21 @@ subroutine updateDeltaDQAtom(this, over, rho, orb, iNeighbour, nNeighbourSK, img !> reference charges real(dp), intent(in) :: q0(:,:,:) - real(dp), allocatable :: mqPerOrbital(:,:) - real(dp), allocatable :: dqPerOrbital(:,:,:) - real(dp), allocatable :: qqPerOrbital(:,:,:,:) - real(dp), allocatable :: mq(:) - real(dp), allocatable :: dq(:,:) - real(dp), allocatable :: qq(:,:,:) + ! real(dp), allocatable :: mqPerOrbital(:,:) + ! real(dp), allocatable :: mq(:) + ! real(dp), allocatable :: dqPerOrbital(:,:,:) + ! real(dp), allocatable :: dq(:,:) + ! real(dp), allocatable :: qqPerOrbital(:,:,:,:) + ! real(dp), allocatable :: qq(:,:,:) real(dp), allocatable :: tmpOvr(:,:), tmpDRho(:,:) - integer, parameter :: xyzLen = 3, iStart = 1, iEnd = 2, iNOrb = 3 + integer, parameter :: iStart = 1, iEnd = 2, iNOrb = 3 - integer :: matrixSize, nAtoms, maxNOrb + integer :: nAtoms integer :: iAt1, iAt2, iSp1, iSp2 integer :: ii, jj, mu, nu, mm, nn, kk, iStartOrb, iEndOrb - real(dp) :: tmpVx3(xyzLen), tmpM3x3(xyzLen, xyzLen) + real(dp) :: tmpVx3(3), tmpM3x3(3,3) real(dp) :: dPSSdP, tmpTrace integer :: iOrig @@ -781,14 +756,8 @@ subroutine updateDeltaDQAtom(this, over, rho, orb, iNeighbour, nNeighbourSK, img nAtom = this%nAtoms @:ASSERT(size(over) == size(rho)) - ! matrixSize = size(overlap, dim = 1) - ! maxNOrb = maxval(orb%nOrbSpecies(:)) - ! maxNOrb = orb%mOrb - - allocate(mqPerOrbital(orb%mOrb, nAtom)) - allocate(mq(nAtom)) - mqPerOrbital(:,:) = 0.0_dp - mq(:) = 0.0_dp + ! allocate(mqPerOrbital(orb%mOrb, nAtom), source=0.0_dp) + ! allocate(mq(nAtom), source=0.0_dp) this%deltaDAtom(:,:) = 0.0_dp this%deltaQAtom(:,:,:) = 0.0_dp @@ -805,11 +774,11 @@ subroutine updateDeltaDQAtom(this, over, rho, orb, iNeighbour, nNeighbourSK, img mulTmp(:) = 0.0_dp mulTmp(1:nOrb1*nOrb2) = over(iOrig:iOrig+nOrb1*nOrb2-1) * rho(iOrig:iOrig+nOrb1*nOrb2-1) sqrTmp(1:nOrb2,1:nOrb1) = reshape(mulTmp(1:nOrb1*nOrb2), (/nOrb2,nOrb1/)) - mqPerOrbital(1:nOrb1,iAtom1) = mqPerOrbital(1:nOrb1,iAtom1) + sum(sqrTmp(1:nOrb2,1:nOrb1), dim=1) - ! Add contribution to the other triangle sum, using the symmetry, but only when off diagonal - if (iAtom1 /= iAtom2f) then - mqPerOrbital(1:nOrb2,iAtom2f) = mqPerOrbital(1:nOrb2,iAtom2f) + sum(sqrTmp(1:nOrb2,1:nOrb1), dim=2) - end if + ! mqPerOrbital(1:nOrb1,iAtom1) = mqPerOrbital(1:nOrb1,iAtom1) + sum(sqrTmp(1:nOrb2,1:nOrb1), dim=1) + ! ! Add contribution to the other triangle sum, using the symmetry, but only when off diagonal + ! if (iAtom1 /= iAtom2f) then + ! mqPerOrbital(1:nOrb2,iAtom2f) = mqPerOrbital(1:nOrb2,iAtom2f) + sum(sqrTmp(1:nOrb2,1:nOrb1), dim=2) + ! end if end do iSp1 = this%species(iAtom1) @@ -856,7 +825,7 @@ subroutine updateDeltaDQAtom(this, over, rho, orb, iNeighbour, nNeighbourSK, img end do end do - mq(:) = mq + sum(mqPerOrbital, dim=1) + ! mq(:) = mq + sum(mqPerOrbital, dim=1) ! do iAtom1 = 1, nAtom ! nOrb1 = orb%nOrbAtom(iAtom1) @@ -876,40 +845,34 @@ end subroutine updateDeltaDQAtom subroutine pullDeltaDQAtom(this, multiExpanData) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> Multipole moments to pull type(TMultipole), intent(inout) :: multiExpanData - integer, parameter :: xyzLen = 3 integer :: nAtoms integer :: iAt1, iSp1, iSpin - real(dp) :: tmpVx3(xyzLen), tmpM3x3(xyzLen, xyzLen) + real(dp) :: tmpVx3(3), tmpM3x3(3,3) iSpin = 1 nAtoms = this%nAtoms do iAt1 = 1, nAtoms iSp1 = this%species(iAt1) - if ((this%onDQOCharges(1, iSp1) == 0) .and. (this%onDQOCharges(2, iSp1) == 0)) then - cycle - end if + ! if ((this%onDQOCharges(1, iSp1) == 0) .and. (this%onDQOCharges(2, iSp1) == 0)) then + ! cycle + ! end if this%deltaDAtom(:,iAt1) = multiExpanData%dipoleAtom(:,iAt1,iSpin) - ! this%deltaDAtom(1,iAt1) = multiExpanData%dipoleAtom(1,iAt1,iSpin) - ! this%deltaDAtom(2,iAt1) = multiExpanData%dipoleAtom(2,iAt1,iSpin) - ! this%deltaDAtom(3,iAt1) = multiExpanData%dipoleAtom(3,iAt1,iSpin) tmpM3x3(:,:) = 0.0_dp + !> Quadrupole components used (xx, xy, yy, xz, yz, zz) tmpM3x3(1,1) = multiExpanData%quadrupoleAtom(1,iAt1,iSpin) - tmpM3x3(1,2) = multiExpanData%quadrupoleAtom(2,iAt1,iSpin) - ! tmpM3x3(2,1) = tmpM3x3(1,2) + tmpM3x3(2,1) = multiExpanData%quadrupoleAtom(2,iAt1,iSpin) tmpM3x3(2,2) = multiExpanData%quadrupoleAtom(3,iAt1,iSpin) - tmpM3x3(1,3) = multiExpanData%quadrupoleAtom(4,iAt1,iSpin) - ! tmpM3x3(3,1) = tmpM3x3(1,3) - tmpM3x3(2,3) = multiExpanData%quadrupoleAtom(5,iAt1,iSpin) - ! tmpM3x3(3,2) = tmpM3x3(2,3) + tmpM3x3(3,1) = multiExpanData%quadrupoleAtom(4,iAt1,iSpin) + tmpM3x3(3,2) = multiExpanData%quadrupoleAtom(5,iAt1,iSpin) tmpM3x3(3,3) = multiExpanData%quadrupoleAtom(6,iAt1,iSpin) - call symmetrizeSquareMatrix(tmpM3x3) + call adjointLowerTriangle(tmpM3x3) this%deltaQAtom(:,:,iAt1) = tmpM3x3(:,:) end do @@ -918,35 +881,31 @@ end subroutine pullDeltaDQAtom subroutine pushDeltaDQAtom(this, multiExpanData) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> Multipole moments push type(TMultipole), intent(inout) :: multiExpanData - integer, parameter :: xyzLen = 3 integer :: nAtoms integer :: iAt1, iSp1, iSpin - real(dp) :: tmpVx3(xyzLen), tmpM3x3(xyzLen, xyzLen) + real(dp) :: tmpVx3(3), tmpM3x3(3,3) iSpin = 1 nAtoms = this%nAtoms do iAt1 = 1, nAtoms iSp1 = this%species(iAt1) - if ((this%onDQOCharges(1, iSp1) == 0) .and. (this%onDQOCharges(2, iSp1) == 0)) then - cycle - end if + ! if ((this%onDQOCharges(1, iSp1) == 0) .and. (this%onDQOCharges(2, iSp1) == 0)) then + ! cycle + ! end if multiExpanData%dipoleAtom(:,iAt1,iSpin) = this%deltaDAtom(:,iAt1) - ! multiExpanData%dipoleAtom(1,iAt1,iSpin) = this%deltaDAtom(1,iAt1) - ! multiExpanData%dipoleAtom(2,iAt1,iSpin) = this%deltaDAtom(2,iAt1) - ! multiExpanData%dipoleAtom(3,iAt1,iSpin) = this%deltaDAtom(3,iAt1) tmpM3x3(:,:) = this%deltaQAtom(:,:,iAt1) multiExpanData%quadrupoleAtom(1,iAt1,iSpin) = tmpM3x3(1,1) - multiExpanData%quadrupoleAtom(2,iAt1,iSpin) = tmpM3x3(1,2) + multiExpanData%quadrupoleAtom(2,iAt1,iSpin) = tmpM3x3(2,1) multiExpanData%quadrupoleAtom(3,iAt1,iSpin) = tmpM3x3(2,2) - multiExpanData%quadrupoleAtom(4,iAt1,iSpin) = tmpM3x3(1,3) - multiExpanData%quadrupoleAtom(5,iAt1,iSpin) = tmpM3x3(2,3) + multiExpanData%quadrupoleAtom(4,iAt1,iSpin) = tmpM3x3(3,1) + multiExpanData%quadrupoleAtom(5,iAt1,iSpin) = tmpM3x3(3,2) multiExpanData%quadrupoleAtom(6,iAt1,iSpin) = tmpM3x3(3,3) end do @@ -955,14 +914,12 @@ end subroutine pushDeltaDQAtom subroutine updateDQPotentials(this, deltaMAtom) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> Delta charge per atom real(dp), intent(in) :: deltaMAtom(:) - integer, parameter :: xyzLen = 3 integer :: nAtoms, iAt1, iAt2, ii, jj, iSp1, iSp2 - real(dp) :: vRabx3(xyzLen), tmpVx3(xyzLen) this%deltaMAtom(:) = deltaMAtom nAtoms = this%nAtoms @@ -1000,12 +957,12 @@ subroutine updateDQPotentials(this, deltaMAtom) & + sum(this%f10AB(:,iAt2,iAt1) * this%deltaDAtom(:,iAt2)) this%pot20x2Atom(iAt1) = this%pot20x2Atom(iAt1)& & + sum(this%f20AB(:,:,iAt2,iAt1) * this%deltaQAtom(:,:,iAt2)) - do ii = 1, xyzLen + do ii = 1, 3 this%pot11x1Atom(ii,iAt1) = this%pot11x1Atom(ii,iAt1)& & - 2.0_dp * sum(this%f20AB(:,ii,iAt2,iAt1) * this%deltaDAtom(:,iAt2)) this%pot21x2Atom(ii,iAt1) = this%pot21x2Atom(ii,iAt1)& & - 3.0_dp * sum(this%f30AB(:,:,ii,iAt2,iAt1) * this%deltaQAtom(:,:,iAt2)) - do jj = 1, xyzLen + do jj = 1, 3 this%pot21x1Atom(jj,ii,iAt1) = this%pot21x1Atom(jj,ii,iAt1)& & - 3.0_dp * sum(this%f30AB(:,jj,ii,iAt2,iAt1) * this%deltaDAtom(:,iAt2)) end do @@ -1014,7 +971,7 @@ subroutine updateDQPotentials(this, deltaMAtom) end do !$OMP END PARALLEL DO - !> Calculate pot21x2 + ! Calculate pot21x2 this%pot22x2Atom(:,:,:) = 0.0_dp !$OMP PARALLEL DO PRIVATE(iAt1, iSp1, iAt2, iSp2, ii, jj) DEFAULT(SHARED) SCHEDULE(RUNTIME) @@ -1028,8 +985,8 @@ subroutine updateDQPotentials(this, deltaMAtom) if ((this%onDQOCharges(1, iSp2) == 0) .and. (this%onDQOCharges(2, iSp2) == 0)) then cycle end if - do ii = 1, xyzLen - do jj = 1, xyzLen + do ii = 1, 3 + do jj = 1, 3 this%pot22x2Atom(jj,ii,iAt1) = this%pot22x2Atom(jj,ii,iAt1)& & + 6.0_dp * sum(this%f40AB(:,:,jj,ii,iAt2,iAt1) * this%deltaQAtom(:,:,iAt2)) end do @@ -1045,7 +1002,7 @@ subroutine addMultiExpanHamiltonian(this, ham, over, nNeighbour, iNeighbour, spe & nAtom, img2CentCell) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> The resulting Hamiltonian contribution. real(dp), intent(inout) :: ham(:,:) @@ -1074,23 +1031,16 @@ subroutine addMultiExpanHamiltonian(this, ham, over, nNeighbour, iNeighbour, spe !> Index mapping atoms onto the central cell atoms. integer, intent(in) :: img2CentCell(:) - integer, parameter :: xyzLen = 3 - integer :: iAt1, iAt2, img, ind, nBlk, iBlk, iSp1, iSp2, iOrb1, iOrb2, iNeigh - integer :: iOrig + integer :: mu, nu, kappa, mm, nn, kk integer :: iAtom1, iAtom2, iAtom2f + integer :: iSp1, iSp2, iNeigh, iOrig integer :: nOrb1, nOrb2 - integer :: mu, nu, kappa, mm, nn, kk - real(dp) :: sqrTmp(orb%mOrb,orb%mOrb) - real(dp) :: sqrTmpRho(orb%mOrb,orb%mOrb) real(dp) :: sqrTmpOver(orb%mOrb,orb%mOrb) real(dp) :: sqrTmpOverT(orb%mOrb,orb%mOrb) real(dp) :: sqrTmpHam(orb%mOrb,orb%mOrb) - real(dp) :: mulTmp(orb%mOrb**2) - real(dp) :: tmpVx3(xyzLen), tmpM3x3(xyzLen, xyzLen) - real(dp) :: tmpadS(xyzLen), tmpSad(xyzLen) - real(dp) :: tmpaQS(xyzLen, xyzLen), tmpSaQ(xyzLen, xyzLen) - real(dp) :: tmpCoAtMu1(xyzLen, xyzLen), tmpCoAtNu1(xyzLen, xyzLen) - real(dp) :: tmpCoAtMu2(xyzLen, xyzLen), tmpCoAtNu2(xyzLen, xyzLen) + real(dp) :: tmpVx3(3), tmpM3x3(3,3) + real(dp) :: tmpadS(3), tmpSad(3) + real(dp) :: tmpaQS(3,3), tmpSaQ(3,3) @:ASSERT(size(nNeighbour)==nAtom) @:ASSERT(size(iNeighbour,dim=2)==nAtom) @@ -1132,30 +1082,29 @@ subroutine addMultiExpanHamiltonian(this, ham, over, nNeighbour, iNeighbour, spe tmpSaQ(:,:) = tmpSaQ(:,:) + this%atomicQIntgrl(:,:,kk,nu,iSp2) * sqrTmpOver(kk,mu) end do - ! add Monopole-Dipole contribution + ! Add Monopole-Dipole contribution sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) - (this%pot10x1Atom(iAtom1)& & + this%pot10x1Atom(iAtom2f)) * sqrTmpOver(nu,mu) - sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + sum(this%pot10x0Atom(:,iAtom1) * tmpadS(:)) & & + sum(this%pot10x0Atom(:,iAtom2f) * tmpSad(:)) - !>add Dipole-Dipole contribution + ! Add Dipole-Dipole contribution sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + sum(this%pot11x1Atom(:,iAtom1) * tmpadS(:)) & & + sum(this%pot11x1Atom(:,iAtom2f) * tmpSad(:)) - !>add Monopole-Quadrupole contribution + ! Add Monopole-Quadrupole contribution sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + (this%pot20x2Atom(iAtom1)& & + this%pot20x2Atom(iAtom2f)) * sqrTmpOver(nu,mu) sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + sum(this%pot20x0Atom(:,:,iAtom1) * tmpaQS(:,:)) & & + sum(this%pot20x0Atom(:,:,iAtom2f) * tmpSaQ(:,:)) - !>add Dipole-Quadrupole contribution + ! Add Dipole-Quadrupole contribution sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) - sum(this%pot21x2Atom(:,iAtom1) * tmpadS(:)) & & - sum(this%pot21x2Atom(:,iAtom2f) * tmpSad(:)) sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + sum(this%pot21x1Atom(:,:,iAtom1) * tmpaQS(:,:)) & & + sum(this%pot21x1Atom(:,:,iAtom2f) * tmpSaQ(:,:)) - !>add Quadrupole-Quadrupole contribution + ! Add Quadrupole-Quadrupole contribution sqrTmpHam(nu,mu) = sqrTmpHam(nu,mu) + sum(this%pot22x2Atom(:,:,iAtom1) * tmpaQS(:,:)) & & + sum(this%pot22x2Atom(:,:,iAtom2f) * tmpSaQ(:,:)) end do @@ -1173,7 +1122,7 @@ subroutine addMultiExpanEnergy(this, energyPerAtom, energyMD, energyDD, energyMQ & energyQQ, energyTT) !> RangeSep class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this real(dp), intent(out) :: energyPerAtom(:) @@ -1195,40 +1144,35 @@ end subroutine addMultiExpanEnergy !> Returns energy per atom. subroutine getEnergyPerAtom(this, energyPerAtom) - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this real(dp), intent(out) :: energyPerAtom(:) real(dp), allocatable :: EnergyMDAtom(:), EnergyDDAtom(:), EnergyMQAtom(:) real(dp), allocatable :: EnergyDQAtom(:), EnergyQQAtom(:) - integer, parameter :: xyzLen = 3, iStart = 1, iEnd = 2, iNOrb = 3 integer :: nAtoms integer :: iAt1 - @:ASSERT(size(energyPerAtom) == this%nAtoms) nAtoms = this%nAtoms - allocate(EnergyMDAtom(nAtoms)) - allocate(EnergyDDAtom(nAtoms)) - allocate(EnergyMQAtom(nAtoms)) - allocate(EnergyDQAtom(nAtoms)) - allocate(EnergyQQAtom(nAtoms)) - EnergyMDAtom(:) = 0.0_dp - EnergyDDAtom(:) = 0.0_dp - EnergyMQAtom(:) = 0.0_dp - EnergyDQAtom(:) = 0.0_dp - EnergyQQAtom(:) = 0.0_dp + @:ASSERT(size(energyPerAtom) == nAtoms) + + allocate(EnergyMDAtom(nAtoms), source=0.0_dp) + allocate(EnergyDDAtom(nAtoms), source=0.0_dp) + allocate(EnergyMQAtom(nAtoms), source=0.0_dp) + allocate(EnergyDQAtom(nAtoms), source=0.0_dp) + allocate(EnergyQQAtom(nAtoms), source=0.0_dp) do iAt1 = 1, nAtoms - !>Monopole-Dipole contribution + ! Monopole-Dipole contribution EnergyMDAtom(iAt1) = sum(this%pot10x0Atom(:,iAt1) * this%deltaDAtom(:,iAt1)) - !>Dipole-Dipole contribution + ! Dipole-Dipole contribution EnergyDDAtom(iAt1) = 0.5_dp * sum(this%pot11x1Atom(:,iAt1) * this%deltaDAtom(:,iAt1)) - !>Monopole-Quadrupole contribution + ! Monopole-Quadrupole contribution EnergyMQAtom(iAt1) = sum(this%pot20x0Atom(:,:,iAt1) * this%deltaQAtom(:,:,iAt1)) - !>Dipole-Quadrupole contribution + ! Dipole-Quadrupole contribution EnergyDQAtom(iAt1) = sum(this%pot21x1Atom(:,:,iAt1) * this%deltaQAtom(:,:,iAt1)) - !>Quadrupole-Quadrupole contribution + ! Quadrupole-Quadrupole contribution EnergyQQAtom(iAt1) = 0.5_dp * sum(this%pot22x2Atom(:,:,iAt1) * this%deltaQAtom(:,:,iAt1)) end do @@ -1240,20 +1184,16 @@ subroutine getEnergyPerAtom(this, energyPerAtom) this%energyQQ = sum(EnergyQQAtom) this%energyTT = this%energyMD + this%energyDD + this%energyMQ + this%energyDQ + this%energyQQ - deallocate(EnergyMDAtom) - deallocate(EnergyDDAtom) - deallocate(EnergyMQAtom) - deallocate(EnergyDQAtom) - deallocate(EnergyQQAtom) - + deallocate(EnergyMDAtom, EnergyDDAtom, EnergyMQAtom, EnergyDQAtom, EnergyQQAtom) end subroutine getEnergyPerAtom + subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, species,& & iNeighbour, nNeighbourSK, img2CentCell, iPair, orb, coords) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> Gradient on exit. real(dp), intent(inout) :: derivs(:,:) @@ -1288,29 +1228,18 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec !> Coordinate of each atom. real(dp), intent(in) :: coords(:,:) - integer, parameter :: xyzLen = 3 + integer :: ii, jj, ll, mu, nu, kappa, mm, nn, kk + integer :: nOrb1, nOrb2 integer :: iAt1, iAt2, iAt2f integer :: iOrig, iSpin, nSpin, nAtom, iNeigh, iSp1, iSp2 - integer :: nOrb1, nOrb2 - integer :: ii, jj, ll, mu, nu, kappa, mm, nn, kk real(dp) :: sqrDMTmp(orb%mOrb,orb%mOrb), sqrDMTmpT(orb%mOrb,orb%mOrb) real(dp) :: sPrimeTmp(orb%mOrb,orb%mOrb,3) real(dp) :: derivTmp(3) - real(dp) :: rab, vRabx3(3), tmp, tmp3(3) - - real(dp) :: sPrimeTmp12(orb%mOrb,orb%mOrb,3) - real(dp) :: tmpDerivMuNu, tmpDeriv(xyzLen), tmpVx3(xyzLen), tmpM3x3(xyzLen, xyzLen) - + real(dp) :: tmpDerivMuNu - ! real(dp) :: sqrTmp(orb%mOrb,orb%mOrb) - ! real(dp) :: sqrTmpRho(orb%mOrb,orb%mOrb) - ! real(dp) :: sqrTmpOver(orb%mOrb,orb%mOrb) - ! real(dp) :: sqrTmpOverT(orb%mOrb,orb%mOrb) - ! real(dp) :: sqrTmpHam(orb%mOrb,orb%mOrb) - ! real(dp) :: mulTmp(orb%mOrb**2) - real(dp) :: tmpPad(xyzLen), tmpadP(xyzLen) - real(dp) :: tmpPaQ(xyzLen, xyzLen), tmpaQP(xyzLen, xyzLen) + real(dp) :: tmpPad(3), tmpadP(3) + real(dp) :: tmpPaQ(3,3), tmpaQP(3,3) !> Calculate pot30x0, pot31x1, and pot32x2 this%pot30x0Atom(:,:,:,:) = 0.0_dp @@ -1328,9 +1257,9 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec & .or. (this%onDQOCharges(1, iSp2) == 0 .and. this%onDQOCharges(2, iSp2) == 0)) then cycle end if - do ii = 1, xyzLen - do jj = 1, xyzLen - do ll = 1, xyzLen + do ii = 1, 3 + do jj = 1, 3 + do ll = 1, 3 this%pot31x1Atom(ll,jj,ii,iAt1) = this%pot31x1Atom(ll,jj,ii,iAt1)& & - 4.0_dp * sum(this%f40AB(:,ll,jj,ii,iAt2,iAt1) * this%deltaDAtom(:,iAt2)) this%pot32x2Atom(ll,jj,ii,iAt1) = this%pot32x2Atom(ll,jj,ii,iAt1)& @@ -1343,27 +1272,27 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec !$OMP END PARALLEL DO do iAt1 = 1, this%nAtoms - !> M-D contribution + ! M-D contribution derivs(:,iAt1) = derivs(:,iAt1) + this%pot11x1Atom(:,iAt1) * this%deltaMAtom(iAt1) - !> M-Q contribution + ! M-Q contribution derivs(:,iAt1) = derivs(:,iAt1) - this%pot21x2Atom(:,iAt1) * this%deltaMAtom(iAt1) - do ii = 1, xyzLen - !> M-D contribution + do ii = 1, 3 + ! D-M contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 2.0_dp * sum(this%pot20x0Atom(:,ii,iAt1) * this%deltaDAtom(:,iAt1)) - !> D-D contribution + ! D-D contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 2.0_dp * sum(this%pot21x1Atom(:,ii,iAt1) * this%deltaDAtom(:,iAt1)) - !> D-Q contribution + ! D-Q contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 2.0_dp * sum(this%pot22x2Atom(:,ii,iAt1) * this%deltaDAtom(:,iAt1)) - !> M-Q contribution + ! Q-M contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 3.0_dp * sum(this%pot30x0Atom(:,:,ii,iAt1) * this%deltaQAtom(:,:,iAt1)) - !> D-Q contribution + ! Q-D contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 3.0_dp * sum(this%pot31x1Atom(:,:,ii,iAt1) * this%deltaQAtom(:,:,iAt1)) - !> Q-Q contribution + ! Q-Q contribution derivs(ii,iAt1) = derivs(ii,iAt1)& & + 3.0_dp * sum(this%pot32x2Atom(:,:,ii,iAt1) * this%deltaQAtom(:,:,iAt1)) end do @@ -1383,26 +1312,6 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec sqrDMTmpT(1:nOrb1,1:nOrb2) = transpose(sqrDMTmp(1:nOrb2,1:nOrb1)) call derivator%getFirstDeriv(sPrimeTmp, skOverCont, coords, species, iAt1, iAt2, orb) - tmpDeriv(:) = 0.0_dp - derivTmp(:) = 0.0_dp - ! note factor of 2 for implicit summation over lower triangle of density matrix: - ! do ii = 1, 3 - ! derivTmp(ii) = 2.0_dp * (& - ! & - sum(sqrEDMTmp(1:nOrb2,1:nOrb1)*sPrimeTmp(1:nOrb2,1:nOrb1,ii))) - ! end do - - ! iSpin = 1 - ! do ii = 1, 3 - ! shiftSprime(1:nOrb2,1:nOrb1) = 0.5_dp * (& - ! & matmul(sPrimeTmp(1:nOrb2,1:nOrb1,ii), shift(1:nOrb1,1:nOrb1,iAt1,iSpin) )& - ! & + matmul(shift(1:nOrb2,1:nOrb2,iAt2f,iSpin), sPrimeTmp(1:nOrb2,1:nOrb1,ii))) - ! ! again factor of 2 from lower triangle, cf published force expressions for SCC: - ! derivTmp(ii) = derivTmp(ii) + 2.0_dp * ( sum(shiftSprime(1:nOrb2,1:nOrb1) *& - ! & reshape(rho(iOrig:iOrig+nOrb1*nOrb2-1), [nOrb2,nOrb1]) ) ) - ! end do - - ! vRabx3(:) = coords(:,iAt1) - coords(:,iAt2) - ! rab = norm2(vRabx3(:)) derivTmp(:) = 0.0_dp do mu = 1, nOrb1 do nu = 1, nOrb2 @@ -1412,50 +1321,49 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec tmpPaQ(:,:) = 0.0_dp tmpaQP(:,:) = 0.0_dp do kk = 1, nOrb1 - tmpPad(:) = tmpPad(:) + this%atomicDIntgrl(:,kk,mu,iSp1) * sqrDMTmpT(kk,nu) - tmpPaQ(:,:) = tmpPaQ(:,:) + this%atomicQIntgrl(:,:,kk,mu,iSp1) * sqrDMTmpT(kk,nu) + tmpPad(:) = tmpPad + this%atomicDIntgrl(:,kk,mu,iSp1) * sqrDMTmpT(kk,nu) + tmpPaQ(:,:) = tmpPaQ + this%atomicQIntgrl(:,:,kk,mu,iSp1) * sqrDMTmpT(kk,nu) end do do kk = 1, nOrb2 - tmpadP(:) = tmpadP(:) + this%atomicDIntgrl(:,kk,nu,iSp2) * sqrDMTmp(kk,mu) - tmpaQP(:,:) = tmpaQP(:,:) + this%atomicQIntgrl(:,:,kk,nu,iSp2) * sqrDMTmp(kk,mu) + tmpadP(:) = tmpadP + this%atomicDIntgrl(:,kk,nu,iSp2) * sqrDMTmp(kk,mu) + tmpaQP(:,:) = tmpaQP + this%atomicQIntgrl(:,:,kk,nu,iSp2) * sqrDMTmp(kk,mu) end do tmpDerivMuNu = 0.0_dp - !> M-D contribution + ! M-D contribution tmpDerivMuNu = tmpDerivMuNu - (this%pot10x1Atom(iAt1) + this%pot10x1Atom(iAt2))& & * sqrDMTmp(nu,mu) - !> M-Q contribution + ! M-Q contribution tmpDerivMuNu = tmpDerivMuNu + (this%pot20x2Atom(iAt1) + this%pot20x2Atom(iAt2))& & * sqrDMTmp(nu,mu) - !> M-D contribution + ! D-M contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot10x0Atom(:,iAt1) * tmpPad(:)) - !> D-D contribution + ! D-D contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot11x1Atom(:,iAt1) * tmpPad(:)) - !> D-Q contribution + ! D-Q contribution tmpDerivMuNu = tmpDerivMuNu - sum(this%pot21x2Atom(:,iAt1) * tmpPad(:)) - !> M-Q contribution + ! M-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot20x0Atom(:,:,iAt1) * tmpPaQ(:,:)) - !> D-Q contribution + ! D-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot21x1Atom(:,:,iAt1) * tmpPaQ(:,:)) - !> Q-Q contribution + ! Q-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot22x2Atom(:,:,iAt1) * tmpPaQ(:,:)) - !> M-D contribution + ! M-D contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot10x0Atom(:,iAt2) * tmpadP(:)) - !> D-D contribution + ! D-D contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot11x1Atom(:,iAt2) * tmpadP(:)) - !> D-Q contribution + ! D-Q contribution tmpDerivMuNu = tmpDerivMuNu - sum(this%pot21x2Atom(:,iAt2) * tmpadP(:)) - !> M-Q contribution + ! M-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot20x0Atom(:,:,iAt2) * tmpaQP(:,:)) - !> D-Q contribution + ! D-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot21x1Atom(:,:,iAt2) * tmpaQP(:,:)) - !> Q-Q contribution + ! Q-Q contribution tmpDerivMuNu = tmpDerivMuNu + sum(this%pot22x2Atom(:,:,iAt2) * tmpaQP(:,:)) - ! tmpDeriv(:) = tmpDeriv(:) + tmpDerivMuNu * sPrimeTmp12(nn,mm,:) - derivTmp(:) = derivTmp(:) + tmpDerivMuNu * sPrimeTmp(nu,mu,:) + derivTmp(:) = derivTmp + tmpDerivMuNu * sPrimeTmp(nu,mu,:) end do end do @@ -1464,24 +1372,19 @@ subroutine addMultiExpanGradients(this, derivs, derivator, skOverCont, rho, spec derivs(:,iAt2f) = derivs(:,iAt2f) - derivTmp end do end do + end subroutine addMultiExpanGradients subroutine addAtomicDipoleMoment(this, dipoleMoment) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> energy gradients real(dp), intent(inout) :: dipoleMoment(:) - integer :: nAtoms - integer :: iAt1 - - nAtoms = this%nAtoms - do iAt1 = 1, nAtoms - dipoleMoment(:) = dipoleMoment(:) - this%deltaDAtom(:,iAt1) - end do + dipoleMoment(:) = dipoleMoment - sum(this%deltaDAtom(:,:), dim=2) end subroutine addAtomicDipoleMoment @@ -1489,18 +1392,12 @@ end subroutine addAtomicDipoleMoment subroutine addAtomicQuadrupoleMoment(this, quadrupoleMoment) !> class instance - class(TDftbMultiExpan), intent(inout) :: this + class(TMdftb), intent(inout) :: this !> energy gradients real(dp), intent(inout) :: quadrupoleMoment(:,:) - integer :: nAtoms - integer :: iAt1 - - nAtoms = this%nAtoms - do iAt1 = 1, nAtoms - quadrupoleMoment(:,:) = quadrupoleMoment(:,:) - this%deltaQAtom(:,:,iAt1) - end do + quadrupoleMoment(:,:) = quadrupoleMoment - sum(this%deltaQAtom(:,:,:), dim=3) end subroutine addAtomicQuadrupoleMoment @@ -1516,8 +1413,8 @@ subroutine makeA3oB3(M3x3, A3, B3) @:ASSERT(size(B3) == dimSize) @:ASSERT(size(M3x3, dim=1) == dimSize) @:ASSERT(size(M3x3, dim=2) == dimSize) - do ii = 1, dimSize - do jj = 1, dimSize + do jj = 1, dimSize + do ii = 1, dimSize M3x3(ii, jj) = A3(ii) * B3(jj) end do end do @@ -1538,9 +1435,9 @@ subroutine makeA3oB3oC3(M3x3x3, A3, B3, C3) @:ASSERT(size(M3x3x3, dim=1) == dimSize) @:ASSERT(size(M3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3, dim=3) == dimSize) - do ii = 1, dimSize + do ll = 1, dimSize do jj = 1, dimSize - do ll = 1, dimSize + do ii = 1, dimSize M3x3x3(ii,jj,ll) = A3(ii) * B3(jj) * C3(ll) end do end do @@ -1562,9 +1459,9 @@ subroutine makeA3oB3x3(M3x3x3, A3, B3x3) @:ASSERT(size(M3x3x3, dim=1) == dimSize) @:ASSERT(size(M3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3, dim=3) == dimSize) - do ii = 1, dimSize + do ll = 1, dimSize do jj = 1, dimSize - do ll = 1, dimSize + do ii = 1, dimSize M3x3x3(ii,jj,ll) = A3(ii) * B3x3(jj,ll) end do end do @@ -1586,9 +1483,9 @@ subroutine makeA3x3oB3(M3x3x3, A3x3, B3) @:ASSERT(size(M3x3x3, dim=1) == dimSize) @:ASSERT(size(M3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3, dim=3) == dimSize) - do ii = 1, dimSize + do ll = 1, dimSize do jj = 1, dimSize - do ll = 1, dimSize + do ii = 1, dimSize M3x3x3(ii,jj,ll) = A3x3(ii,jj) * B3(ll) end do end do @@ -1610,9 +1507,9 @@ subroutine makeA3x3otB3(M3x3x3, A3x3, B3) @:ASSERT(size(M3x3x3, dim=1) == dimSize) @:ASSERT(size(M3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3, dim=3) == dimSize) - do ii = 1, dimSize + do ll = 1, dimSize do jj = 1, dimSize - do ll = 1, dimSize + do ii = 1, dimSize M3x3x3(ii,jj,ll) = A3x3(ii,ll) * B3(jj) end do end do @@ -1636,10 +1533,10 @@ subroutine makeA3x3oB3x3(M3x3x3x3, A3x3, B3x3) @:ASSERT(size(M3x3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3x3, dim=3) == dimSize) @:ASSERT(size(M3x3x3x3, dim=4) == dimSize) - do ii = 1, dimSize - do jj = 1, dimSize - do ll = 1, dimSize - do mm = 1, dimSize + do mm = 1, dimSize + do ll = 1, dimSize + do jj = 1, dimSize + do ii = 1, dimSize M3x3x3x3(ii,jj,ll,mm) = A3x3(ii,jj) * B3x3(ll,mm) end do end do @@ -1664,10 +1561,10 @@ subroutine makeA3x3ohB3x3(M3x3x3x3, A3x3, B3x3) @:ASSERT(size(M3x3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3x3, dim=3) == dimSize) @:ASSERT(size(M3x3x3x3, dim=4) == dimSize) - do ii = 1, dimSize - do jj = 1, dimSize - do ll = 1, dimSize - do mm = 1, dimSize + do mm = 1, dimSize + do ll = 1, dimSize + do jj = 1, dimSize + do ii = 1, dimSize M3x3x3x3(ii,jj,ll,mm) = A3x3(ii,ll) * B3x3(jj,mm) end do end do @@ -1692,10 +1589,10 @@ subroutine makeA3x3otB3x3(M3x3x3x3, A3x3, B3x3) @:ASSERT(size(M3x3x3x3, dim=2) == dimSize) @:ASSERT(size(M3x3x3x3, dim=3) == dimSize) @:ASSERT(size(M3x3x3x3, dim=4) == dimSize) - do ii = 1, dimSize - do jj = 1, dimSize - do ll = 1, dimSize - do mm = 1, dimSize + do mm = 1, dimSize + do ll = 1, dimSize + do jj = 1, dimSize + do ii = 1, dimSize M3x3x3x3(ii,jj,ll,mm) = A3x3(ii,mm) * B3x3(jj,ll) end do end do @@ -1704,39 +1601,4 @@ subroutine makeA3x3otB3x3(M3x3x3x3, A3x3, B3x3) end subroutine makeA3x3otB3x3 - !> copy lower triangle to upper for a square matrix - subroutine symmetrizeSquareMatrix(matrix) - - !> matrix to symmetrize - real(dp), intent(inout) :: matrix(:,:) - integer :: ii, matSize - - @:ASSERT(size(matrix, dim=1) == size(matrix, dim=2)) - matSize = size(matrix, dim = 1) - - !!$OMP PARALLEL DO PRIVATE(ii) DEFAULT(SHARED) SCHEDULE(RUNTIME) - do ii = 1, matSize - 1 - matrix(ii,ii+1:matSize) = matrix(ii + 1 : matSize, ii) - end do - !!$OMP END PARALLEL DO - - end subroutine symmetrizeSquareMatrix - - - !> location of relevant atomic block indices in a dense matrix - pure function getDescriptor(iAt, iSquare) result(desc) - - !> relevant atom - integer, intent(in) :: iAt - - !> indexing array for start of atom orbitals - integer, intent(in) :: iSquare(:) - - !> resulting location ranges - integer :: desc(3) - - desc(:) = [iSquare(iAt), iSquare(iAt + 1) - 1, iSquare(iAt + 1) - iSquare(iAt)] - - end function getDescriptor - end module dftbp_dftb_multiexpan diff --git a/src/dftbp/dftbplus/initprogram.F90 b/src/dftbp/dftbplus/initprogram.F90 index 23977bd14f..33a4798f99 100644 --- a/src/dftbp/dftbplus/initprogram.F90 +++ b/src/dftbp/dftbplus/initprogram.F90 @@ -41,7 +41,7 @@ module dftbp_dftbplus_initprogram use dftbp_dftb_h5correction, only : TH5CorrectionInput use dftbp_dftb_halogenx, only : THalogenX, THalogenX_init use dftbp_dftb_hamiltonian, only : TRefExtPot - use dftbp_dftb_multiexpan, only : TDftbMultiExpan, TDftbMultiExpanInp, dftbMultiExpan_init + use dftbp_dftb_multiexpan, only : TMdftb, TMdftbAtomicIntegrals, TMdftbInp, TMdftb_init use dftbp_dftb_nonscc, only : TNonSccDiff, NonSccDiff_init, diffTypes use dftbp_dftb_onsitecorrection, only : Ons_getOrbitalEquiv, Ons_blockIndx use dftbp_dftb_orbitalequiv, only : OrbitalEquiv_merge, OrbitalEquiv_reduce @@ -805,10 +805,10 @@ module dftbp_dftbplus_initprogram logical :: isDftbMultiExpan !> input data structure for multipole - type(TDftbMultiExpanInp) :: dftbMultiExpanInp + type(TMdftbInp) :: dftbMultiExpanInp !> data structure for multipole - type(TDftbMultiExpan), allocatable :: dftbMultiExpan + type(TMdftb), allocatable :: dftbMultiExpan !> If initial charges/dens mtx. from external file. logical :: tReadChrg @@ -1773,62 +1773,12 @@ subroutine initProgramVariables(this, input, env) this%dftbMultiExpanInp%nOrb = this%nOrb this%dftbMultiExpanInp%nSpin = this%nSpin this%dftbMultiExpanInp%nSpecies = this%nType - - allocate(this%dftbMultiExpanInp%atomicDIntgrlScaling(this%nType)) - allocate(this%dftbMultiExpanInp%atomicQIntgrlScaling(this%nType)) - allocate(this%dftbMultiExpanInp%atomicOnsiteScaling(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSXPxIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPxXDxxyyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPxXDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPyYDxxyyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPzZDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSXXSIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPxXXPxIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicPyXXPyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSXXDxxyyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSXXDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSYYDxxyyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicSZZDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDxyXXDxyIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDyzXXDyzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDxxyyXXDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDzzXXDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDxxyyYYDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDzzZZDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDxzXZDzzIntgrl(this%nType)) - allocate(this%dftbMultiExpanInp%atomicDyzYZDxxyyIntgrl(this%nType)) - - this%dftbMultiExpanInp%atomicDIntgrlScaling(:) = input%ctrl%atomicDIntgrlScaling - this%dftbMultiExpanInp%atomicQIntgrlScaling(:) = input%ctrl%atomicQIntgrlScaling - this%dftbMultiExpanInp%atomicOnsiteScaling(:) = input%ctrl%atomicOnsiteScaling - this%dftbMultiExpanInp%atomicSXPxIntgrl(:) = input%ctrl%atomicSXPxIntgrl - this%dftbMultiExpanInp%atomicPxXDxxyyIntgrl(:) = input%ctrl%atomicPxXDxxyyIntgrl - this%dftbMultiExpanInp%atomicPxXDzzIntgrl(:) = input%ctrl%atomicPxXDzzIntgrl - this%dftbMultiExpanInp%atomicPyYDxxyyIntgrl(:) = input%ctrl%atomicPyYDxxyyIntgrl - this%dftbMultiExpanInp%atomicPzZDzzIntgrl(:) = input%ctrl%atomicPzZDzzIntgrl - this%dftbMultiExpanInp%atomicSXXSIntgrl(:) = input%ctrl%atomicSXXSIntgrl - this%dftbMultiExpanInp%atomicPxXXPxIntgrl(:) = input%ctrl%atomicPxXXPxIntgrl - this%dftbMultiExpanInp%atomicPyXXPyIntgrl(:) = input%ctrl%atomicPyXXPyIntgrl - this%dftbMultiExpanInp%atomicSXXDxxyyIntgrl(:) = input%ctrl%atomicSXXDxxyyIntgrl - this%dftbMultiExpanInp%atomicSXXDzzIntgrl(:) = input%ctrl%atomicSXXDzzIntgrl - this%dftbMultiExpanInp%atomicSYYDxxyyIntgrl(:) = input%ctrl%atomicSYYDxxyyIntgrl - this%dftbMultiExpanInp%atomicSZZDzzIntgrl(:) = input%ctrl%atomicSZZDzzIntgrl - this%dftbMultiExpanInp%atomicDxyXXDxyIntgrl(:) = input%ctrl%atomicDxyXXDxyIntgrl - this%dftbMultiExpanInp%atomicDyzXXDyzIntgrl(:) = input%ctrl%atomicDyzXXDyzIntgrl - this%dftbMultiExpanInp%atomicDxxyyXXDzzIntgrl(:) = input%ctrl%atomicDxxyyXXDzzIntgrl - this%dftbMultiExpanInp%atomicDzzXXDzzIntgrl(:) = input%ctrl%atomicDzzXXDzzIntgrl - this%dftbMultiExpanInp%atomicDxxyyYYDzzIntgrl(:) = input%ctrl%atomicDxxyyYYDzzIntgrl - this%dftbMultiExpanInp%atomicDzzZZDzzIntgrl(:) = input%ctrl%atomicDzzZZDzzIntgrl - this%dftbMultiExpanInp%atomicDxzXZDzzIntgrl(:) = input%ctrl%atomicDxzXZDzzIntgrl - this%dftbMultiExpanInp%atomicDyzYZDxxyyIntgrl(:) = input%ctrl%atomicDyzYZDxxyyIntgrl - - allocate(this%dftbMultiExpanInp%hubbU(size(hubbU, dim=2))) this%dftbMultiExpanInp%hubbU = hubbU(1,:) - allocate(this%dftbMultiExpanInp%species(this%nAtom)) - this%dftbMultiExpanInp%species(:) = input%geom%species + this%dftbMultiExpanInp%species = input%geom%species(:) + call initMdftbAtomicIntegrals(this%dftbMultiExpanInp, input%ctrl) allocate(this%dftbMultiExpan) - allocate(this%quadrupoleMoment(3,3)) - call dftbMultiExpan_init(this%dftbMultiExpan, this%dftbMultiExpanInp) + allocate(this%quadrupoleMoment(3,3), source=0.0_dp) + call TMdftb_init(this%dftbMultiExpan, this%dftbMultiExpanInp) end if end if @@ -3420,7 +3370,7 @@ subroutine initProgramVariables(this, input, env) ! write(stdout, "(A,':',T30,E14.6)") "Ewald alpha parameter", this%scc%getEwaldPar() !end if if (this%isDftbMultiExpan) then - write(stdOut, "(A,':',T30,A)") "MultiPole expansion", "Yes" + write(stdOut, "(A,':',T30,A)") "Multipole expansion", "Yes" end if if (input%ctrl%tShellResolved) then write(stdOut, "(A,':',T30,A)") "Shell resolved Hubbard", "Yes" @@ -7034,6 +6984,39 @@ subroutine initCoulombInput_(env, ewaldAlpha, tolEwald, boundaryCond, coulombInp end subroutine initCoulombInput_ + subroutine initMdftbAtomicIntegrals(dftbMultiExpanInp, ctrl) + + !> input data structure for multipole + type(TMdftbInp), intent(inout) :: dftbMultiExpanInp + + !> Data control structure + type(TControl), intent(in) :: ctrl + + dftbMultiExpanInp%atomicDIntgrlScaling = ctrl%mdftbAtomicIntegrals%DScaling(:) + dftbMultiExpanInp%atomicQIntgrlScaling = ctrl%mdftbAtomicIntegrals%QScaling(:) + dftbMultiExpanInp%atomicSXPxIntgrl = ctrl%mdftbAtomicIntegrals%SXPx(:) + dftbMultiExpanInp%atomicPxXDxxyyIntgrl = ctrl%mdftbAtomicIntegrals%PxXDxxyy(:) + dftbMultiExpanInp%atomicPxXDzzIntgrl = ctrl%mdftbAtomicIntegrals%PxXDzz(:) + dftbMultiExpanInp%atomicPyYDxxyyIntgrl = ctrl%mdftbAtomicIntegrals%PyYDxxyy(:) + dftbMultiExpanInp%atomicPzZDzzIntgrl = ctrl%mdftbAtomicIntegrals%PzZDzz(:) + dftbMultiExpanInp%atomicSXXSIntgrl = ctrl%mdftbAtomicIntegrals%SXXS(:) + dftbMultiExpanInp%atomicPxXXPxIntgrl = ctrl%mdftbAtomicIntegrals%PxXXPx(:) + dftbMultiExpanInp%atomicPyXXPyIntgrl = ctrl%mdftbAtomicIntegrals%PyXXPy(:) + dftbMultiExpanInp%atomicSXXDxxyyIntgrl = ctrl%mdftbAtomicIntegrals%SXXDxxyy(:) + dftbMultiExpanInp%atomicSXXDzzIntgrl = ctrl%mdftbAtomicIntegrals%SXXDzz(:) + dftbMultiExpanInp%atomicSYYDxxyyIntgrl = ctrl%mdftbAtomicIntegrals%SYYDxxyy(:) + dftbMultiExpanInp%atomicSZZDzzIntgrl = ctrl%mdftbAtomicIntegrals%SZZDzz(:) + dftbMultiExpanInp%atomicDxyXXDxyIntgrl = ctrl%mdftbAtomicIntegrals%DxyXXDxy(:) + dftbMultiExpanInp%atomicDyzXXDyzIntgrl = ctrl%mdftbAtomicIntegrals%DyzXXDyz(:) + dftbMultiExpanInp%atomicDxxyyXXDzzIntgrl = ctrl%mdftbAtomicIntegrals%DxxyyXXDzz(:) + dftbMultiExpanInp%atomicDzzXXDzzIntgrl = ctrl%mdftbAtomicIntegrals%DzzXXDzz(:) + dftbMultiExpanInp%atomicDxxyyYYDzzIntgrl = ctrl%mdftbAtomicIntegrals%DxxyyYYDzz(:) + dftbMultiExpanInp%atomicDzzZZDzzIntgrl = ctrl%mdftbAtomicIntegrals%DzzZZDzz(:) + dftbMultiExpanInp%atomicDxzXZDzzIntgrl = ctrl%mdftbAtomicIntegrals%DxzXZDzz(:) + dftbMultiExpanInp%atomicDyzYZDxxyyIntgrl = ctrl%mdftbAtomicIntegrals%DyzYZDxxyy(:) + + end subroutine initMdftbAtomicIntegrals + #:if WITH_POISSON diff --git a/src/dftbp/dftbplus/inputdata.F90 b/src/dftbp/dftbplus/inputdata.F90 index da36a2d179..098e5ef3d8 100644 --- a/src/dftbp/dftbplus/inputdata.F90 +++ b/src/dftbp/dftbplus/inputdata.F90 @@ -24,6 +24,7 @@ module dftbp_dftbplus_inputdata use dftbp_dftb_pmlocalisation, only : TPipekMezeyInp use dftbp_dftb_potentials, only : TAtomExtPotInput use dftbp_dftb_slakocont, only : TSlakoCont + use dftbp_dftb_multiexpan, only : TMdftbAtomicIntegrals use dftbp_dftbplus_input_geoopt, only : TGeoOptInput use dftbp_elecsolvers_elecsolvers, only : TElectronicSolverInp use dftbp_extlibs_poisson, only : TPoissonInfo @@ -552,32 +553,9 @@ module dftbp_dftbplus_inputdata !> Hybrid xc-functional input type(THybridXcInp), allocatable :: hybridXcInp - !> Multipole expansion logical :: isDftbMultiExpan = .false. - real(dp), allocatable :: atomicDIntgrlScaling(:) - real(dp), allocatable :: atomicQIntgrlScaling(:) - real(dp), allocatable :: atomicOnsiteScaling(:) - real(dp), allocatable :: atomicSXPxIntgrl(:) - real(dp), allocatable :: atomicPxXDxxyyIntgrl(:) - real(dp), allocatable :: atomicPxXDzzIntgrl(:) - real(dp), allocatable :: atomicPyYDxxyyIntgrl(:) - real(dp), allocatable :: atomicPzZDzzIntgrl(:) - real(dp), allocatable :: atomicSXXSIntgrl(:) - real(dp), allocatable :: atomicPxXXPxIntgrl(:) - real(dp), allocatable :: atomicPyXXPyIntgrl(:) - real(dp), allocatable :: atomicSXXDxxyyIntgrl(:) - real(dp), allocatable :: atomicSXXDzzIntgrl(:) - real(dp), allocatable :: atomicSYYDxxyyIntgrl(:) - real(dp), allocatable :: atomicSZZDzzIntgrl(:) - real(dp), allocatable :: atomicDxyXXDxyIntgrl(:) - real(dp), allocatable :: atomicDyzXXDyzIntgrl(:) - real(dp), allocatable :: atomicDxxyyXXDzzIntgrl(:) - real(dp), allocatable :: atomicDzzXXDzzIntgrl(:) - real(dp), allocatable :: atomicDxxyyYYDzzIntgrl(:) - real(dp), allocatable :: atomicDzzZZDzzIntgrl(:) - real(dp), allocatable :: atomicDxzXZDzzIntgrl(:) - real(dp), allocatable :: atomicDyzYZDxxyyIntgrl(:) + type(TMdftbAtomicIntegrals), allocatable :: mdftbAtomicIntegrals #:if WITH_SOCKETS !> Socket communication diff --git a/src/dftbp/dftbplus/main.F90 b/src/dftbp/dftbplus/main.F90 index f5d66a5ed8..8e06d2baf5 100644 --- a/src/dftbp/dftbplus/main.F90 +++ b/src/dftbp/dftbplus/main.F90 @@ -34,7 +34,7 @@ module dftbp_dftbplus_main use dftbp_dftb_hamiltonian, only : resetInternalPotentials, addChargePotentials,& & getSccHamiltonian, mergeExternalPotentials, resetExternalPotentials,& & addBlockChargePotentials, constrainSccHamiltonian - use dftbp_dftb_multiexpan, only : TDftbMultiExpan + use dftbp_dftb_multiexpan, only : TMdftb use dftbp_dftb_hybridxc, only : THybridXcFunc, hybridXcAlgo use dftbp_dftb_nonscc, only : TNonSccDiff, buildS, buildH0 use dftbp_dftb_onsitecorrection, only : Onsblock_expand, onsBlock_reduce, addOnsShift @@ -2265,7 +2265,7 @@ subroutine handleCoordinateChange(env, boundaryCond, coord0, latVec, invLatVec, type(TReksCalc), allocatable, intent(inout) :: reks !> multipole contributions - type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan + type(TMdftb), allocatable, intent(inout) :: dftbMultiExpan !> Image atoms to their equivalent in the central cell integer, allocatable, intent(inout) :: img2CentCell(:) @@ -5608,8 +5608,8 @@ subroutine getQuadrupoleMoment(qOutput, q0, coord, quadrupoleMoment, iAtInCentra do ii = 1, size(iAtInCentralRegion) iAtom = iAtInCentralRegion(ii) dqAtom = sum(q0(:, iAtom, 1) - qOutput(:, iAtom, 1)) - do jj = 1, 3 - do ll = 1, 3 + do ll = 1, 3 + do jj = 1, 3 quadrupoleMoment(jj,ll) = quadrupoleMoment(jj,ll)& & + dqAtom * coord(jj,iAtom) * coord(ll,iAtom) end do @@ -6580,7 +6580,7 @@ subroutine getGradients(env, parallelKS, boundaryConds, sccCalc, tblite, isExtFi class(THybridXcFunc), intent(inout), allocatable :: hybridXc !> Multipole contributions - type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan + type(TMdftb), allocatable, intent(inout) :: dftbMultiExpan !> Dense overlap matrix, required for hybridXc real(dp), intent(inout), allocatable :: SSqrReal(:,:) @@ -8208,7 +8208,7 @@ subroutine getHamiltonianLandEnergyL(env, denseDesc, sccCalc, tblite, orb, speci type(TMultipole), intent(inout) :: multipole !> DFTB multipole expansion - type(TDftbMultiExpan), allocatable, intent(inout) :: dftbMultiExpan + type(TMdftb), allocatable, intent(inout) :: dftbMultiExpan !> Data type for REKS type(TReksCalc), allocatable, intent(inout) :: reks diff --git a/src/dftbp/dftbplus/mainio.F90 b/src/dftbp/dftbplus/mainio.F90 index d06c182b5f..ec23429ca1 100644 --- a/src/dftbp/dftbplus/mainio.F90 +++ b/src/dftbp/dftbplus/mainio.F90 @@ -3299,7 +3299,7 @@ subroutine writeDetailedOut3(fd, qInput, qOutput, energy, species, tDFTBU, tPrin !> Is this a solvation model used? logical, intent(in) :: tSolv - !> Are there energy contributions up to quadrupole in the DFTB model? + !> Are there energy contributions from mdftb? logical, intent(in) :: isDftbMultiExpan real(dp), allocatable :: qInputUpDown(:,:,:), qOutputUpDown(:,:,:) @@ -3809,21 +3809,7 @@ subroutine writeDetailedOut7(fd, tGeoOpt, tGeomEnd, tMd, tDerivs, eField, dipole end if if (allocated(quadrupoleMoment)) then - write(fd, "(A)") ' Traceless Quadrupole moment in au' - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XX', quadrupoleMoment(1,1), ' YY',& - & quadrupoleMoment(2,2), ' ZZ', quadrupoleMoment(3,3) - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XY', quadrupoleMoment(1,2), ' XZ',& - & quadrupoleMoment(1,3), ' YZ', quadrupoleMoment(2,3) - write(fd, "(A)") ' Traceless Quadrupole moment in Debye*Ang' - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XX',& - & quadrupoleMoment(1,1) * au__Debye * Bohr__AA, ' YY',& - & quadrupoleMoment(2,2) * au__Debye * Bohr__AA, ' ZZ',& - & quadrupoleMoment(3,3) * au__Debye * Bohr__AA - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XY',& - & quadrupoleMoment(1,2) * au__Debye * Bohr__AA, ' XZ',& - & quadrupoleMoment(1,3) * au__Debye * Bohr__AA, ' YZ',& - & quadrupoleMoment(2,3) * au__Debye * Bohr__AA - write(fd, *) + call printQuadrupoleMoment(quadrupoleMoment, fd) end if if (allocated(eField)) then @@ -4173,21 +4159,7 @@ subroutine writeMdOut2(fd, isPeriodic, printForces, hasStress, withBarostat, isL end if if (allocated(quadrupoleMoment)) then - write(fd, "(A)") ' Traceless Quadrupole moment in au' - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XX', quadrupoleMoment(1,1), ' YY',& - & quadrupoleMoment(2,2), ' ZZ', quadrupoleMoment(3,3) - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XY', quadrupoleMoment(1,2), ' XZ',& - & quadrupoleMoment(1,3), ' YZ', quadrupoleMoment(2,3) - write(fd, "(A)") ' Traceless Quadrupole moment in Debye*Ang' - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XX',& - & quadrupoleMoment(1,1) * au__Debye * Bohr__AA, ' YY',& - & quadrupoleMoment(2,2) * au__Debye * Bohr__AA, ' ZZ',& - & quadrupoleMoment(3,3) * au__Debye * Bohr__AA - write(fd, "(A, F14.8, A, F14.8, A, F14.8)") ' XY',& - & quadrupoleMoment(1,2) * au__Debye * Bohr__AA, ' XZ',& - & quadrupoleMoment(1,3) * au__Debye * Bohr__AA, ' YZ',& - & quadrupoleMoment(2,3) * au__Debye * Bohr__AA - write(fd, *) + call printQuadrupoleMoment(quadrupoleMoment, fd) end if end subroutine writeMdOut2 @@ -5978,4 +5950,39 @@ subroutine writeCosmoFile(solvation, species0, speciesNames, coords0, energy) end subroutine writeCosmoFile + subroutine printQuadrupoleMoment(quadrupoleMoment, outUnit) + + !> quadrupole moment + real(dp), intent(in), allocatable :: quadrupoleMoment(:,:) + + !> Optional unit to print out the quadrupoleMoment + integer, intent(in), optional :: outUnit + + integer :: iUnit + + if (present(outUnit)) then + iUnit = outUnit + else + iUnit = stdOut + end if + + write(iUnit, "(A)") ' Traceless Quadrupole moment in au' + write(iUnit, "(A, F14.8, A, F14.8, A, F14.8)") ' XX', quadrupoleMoment(1,1), ' YY',& + & quadrupoleMoment(2,2), ' ZZ', quadrupoleMoment(3,3) + write(iUnit, "(A, F14.8, A, F14.8, A, F14.8)") ' XY', quadrupoleMoment(1,2), ' XZ',& + & quadrupoleMoment(1,3), ' YZ', quadrupoleMoment(2,3) + + write(iUnit, "(A)") ' Traceless Quadrupole moment in Buckingham or Debye*Ang' + write(iUnit, "(A, F14.8, A, F14.8, A, F14.8)") ' XX',& + & quadrupoleMoment(1,1) * au__Debye * Bohr__AA, ' YY',& + & quadrupoleMoment(2,2) * au__Debye * Bohr__AA, ' ZZ',& + & quadrupoleMoment(3,3) * au__Debye * Bohr__AA + write(iUnit, "(A, F14.8, A, F14.8, A, F14.8)") ' XY',& + & quadrupoleMoment(1,2) * au__Debye * Bohr__AA, ' XZ',& + & quadrupoleMoment(1,3) * au__Debye * Bohr__AA, ' YZ',& + & quadrupoleMoment(2,3) * au__Debye * Bohr__AA + write(iUnit, *) + + end subroutine printQuadrupoleMoment + end module dftbp_dftbplus_mainio diff --git a/src/dftbp/dftbplus/parser.F90 b/src/dftbp/dftbplus/parser.F90 index 4f355432f9..0337e478a9 100644 --- a/src/dftbp/dftbplus/parser.F90 +++ b/src/dftbp/dftbplus/parser.F90 @@ -39,6 +39,7 @@ module dftbp_dftbplus_parser use dftbp_dftb_repulsive_splinerep, only : TSplineRepInp, TSplineRep use dftbp_dftb_slakocont, only : init, addTable use dftbp_dftb_slakoeqgrid, only : skEqGridNew, skEqGridOld, TSlakoEqGrid, init + use dftbp_dftb_multiexpan, only : TMdftbAtomicIntegrals use dftbp_dftbplus_forcetypes, only : forceTypes use dftbp_dftbplus_input_fileaccess, only : readBinaryAccessTypes use dftbp_dftbplus_inputconversion, only : transformpdosregioninfo @@ -1748,136 +1749,7 @@ subroutine readDFTBHam(node, ctrl, geo, slako, poisson, errStatus) ! Multipole expansion ctrl%isDftbMultiExpan = .false. - if (ctrl%tSCC) then - !call getChildValue(node, "Multipole", ctrl%isDftbMultiExpan, .false.) - call getChildValue(node, "Multipole", value1, "", child=child, allowEmptyValue=.true.,& - & dummyValue=.true.) - if (associated(value1)) then - call getNodeName(value1, buffer) - select case(char(buffer)) - case("onecenterapproximation") - ctrl%isDftbMultiExpan = .true. - allocate(ctrl%atomicDIntgrlScaling(geo%nSpecies)) - allocate(ctrl%atomicQIntgrlScaling(geo%nSpecies)) - allocate(ctrl%atomicOnsiteScaling(geo%nSpecies)) - allocate(ctrl%atomicSXPxIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPxXDxxyyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPxXDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPyYDxxyyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPzZDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicSXXSIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPxXXPxIntgrl(geo%nSpecies)) - allocate(ctrl%atomicPyXXPyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicSXXDxxyyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicSXXDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicSYYDxxyyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicSZZDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDxyXXDxyIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDyzXXDyzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDxxyyXXDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDzzXXDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDxxyyYYDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDzzZZDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDxzXZDzzIntgrl(geo%nSpecies)) - allocate(ctrl%atomicDyzYZDxxyyIntgrl(geo%nSpecies)) - ctrl%atomicDIntgrlScaling(:) = 1.0_dp - ctrl%atomicQIntgrlScaling(:) = 1.0_dp - ctrl%atomicOnsiteScaling(:) = 1.0_dp - ctrl%atomicSXPxIntgrl(:) = 0.0_dp - ctrl%atomicPxXDxxyyIntgrl(:) = 0.0_dp - ctrl%atomicPxXDzzIntgrl(:) = 0.0_dp - ctrl%atomicPyYDxxyyIntgrl(:) = 0.0_dp - ctrl%atomicPzZDzzIntgrl(:) = 0.0_dp - ctrl%atomicSXXSIntgrl(:) = 0.0_dp - ctrl%atomicPxXXPxIntgrl(:) = 0.0_dp - ctrl%atomicPyXXPyIntgrl(:) = 0.0_dp - ctrl%atomicSXXDxxyyIntgrl(:) = 0.0_dp - ctrl%atomicSXXDzzIntgrl(:) = 0.0_dp - ctrl%atomicSYYDxxyyIntgrl(:) = 0.0_dp - ctrl%atomicSZZDzzIntgrl(:) = 0.0_dp - ctrl%atomicDxyXXDxyIntgrl(:) = 0.0_dp - ctrl%atomicDyzXXDyzIntgrl(:) = 0.0_dp - ctrl%atomicDxxyyXXDzzIntgrl(:) = 0.0_dp - ctrl%atomicDzzXXDzzIntgrl(:) = 0.0_dp - ctrl%atomicDxxyyYYDzzIntgrl(:) = 0.0_dp - ctrl%atomicDzzZZDzzIntgrl(:) = 0.0_dp - ctrl%atomicDxzXZDzzIntgrl(:) = 0.0_dp - ctrl%atomicDyzYZDxxyyIntgrl(:) = 0.0_dp - - call getChild(value1, 'AtomDIntergralScalings', child2, requested=.false.) - if (associated(child2)) then - do iSp1 = 1, geo%nSpecies - call getChildValue(child2, trim(geo%speciesNames(iSp1)),& - & ctrl%atomicDIntgrlScaling(iSp1), 1.0_dp) - end do - end if - - call getChild(value1, 'AtomQIntergralScalings', child2, requested=.false.) - if (associated(child2)) then - do iSp1 = 1, geo%nSpecies - call getChildValue(child2, trim(geo%speciesNames(iSp1)),& - & ctrl%atomicQIntgrlScaling(iSp1), 1.0_dp) - end do - end if - - call getChild(value1, 'AtomOnsiteScalings', child2, requested=.false.) - if (associated(child2)) then - do iSp1 = 1, geo%nSpecies - call getChildValue(child2, trim(geo%speciesNames(iSp1)),& - & ctrl%atomicOnsiteScaling(iSp1), 1.0_dp) - end do - end if - - call getChild(value1, 'OneCenterAtomIntergrals', child2, requested=.true.) - do iSp1 = 1, geo%nSpecies - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_X_Px",& - & ctrl%atomicSXPxIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_X_Dxx-yy",& - & ctrl%atomicPxXDxxyyIntgrl (iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_X_Dzz",& - & ctrl%atomicPxXDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Py_Y_Dxx-yy",& - & ctrl%atomicPyYDxxyyIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Pz_Z_Dzz",& - & ctrl%atomicPzZDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_S",& - & ctrl%atomicSXXSIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_XX_Px",& - & ctrl%atomicPxXXPxIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Py_XX_Py",& - & ctrl%atomicPyXXPyIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_Dxx-yy",& - & ctrl%atomicSXXDxxyyIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_Dzz",& - & ctrl%atomicSXXDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_YY_Dxx-yy",& - & ctrl%atomicSYYDxxyyIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_ZZ_Dzz",& - & ctrl%atomicSZZDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxy_XX_Dxy",& - & ctrl%atomicDxyXXDxyIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dyz_XX_Dyz",& - & ctrl%atomicDyzXXDyzIntgrl(iSp1), 0.0_dp) - !call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxx-yy_XX_Dzz",& - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_XX_Dxx-yy",& - & ctrl%atomicDxxyyXXDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_XX_Dzz",& - & ctrl%atomicDzzXXDzzIntgrl(iSp1), 0.0_dp) - !call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxx-yy_YY_Dzz",& - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_YY_Dxx-yy",& - & ctrl%atomicDxxyyYYDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_ZZ_Dzz",& - & ctrl%atomicDzzZZDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxz_XZ_Dzz",& - & ctrl%atomicDxzXZDzzIntgrl(iSp1), 0.0_dp) - call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dyz_YZ_Dxx-yy",& - & ctrl%atomicDyzYZDxxyyIntgrl(iSp1), 0.0_dp) - end do - case default - call detailedError(child,"Unknown functions :"// char(buffer)) - end select - end if - end if + call readMdftb(node, ctrl, geo) ! Third order stuff ctrl%t3rd = .false. @@ -2412,6 +2284,122 @@ subroutine readElectrostatics(node, ctrl, geo, poisson) end subroutine readElectrostatics + subroutine readMdftb(node, ctrl, geo) + + !> Node to get the information from + type(fnode), pointer :: node + + !> Control structure to be filled + type(TControl), intent(inout) :: ctrl + + !> Geometry structure to be filled + type(TGeometry), intent(in) :: geo + + type(fnode), pointer :: value1, child, child2 + type(string) :: buffer + integer :: iSp1 + + ctrl%isDftbMultiExpan = .false. + if (ctrl%tSCC) then + call getChildValue(node, "Multipole", value1, "", child=child, allowEmptyValue=.true.,& + & dummyValue=.false.) + if (associated(value1)) then + call getNodeName(value1, buffer) + select case(char(buffer)) + case("onecenterapproximation") + ctrl%isDftbMultiExpan = .true. + allocate(ctrl%mdftbAtomicIntegrals) + allocate(ctrl%mdftbAtomicIntegrals%DScaling(geo%nSpecies), source=1.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%QScaling(geo%nSpecies), source=1.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SXPx(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PxXDxxyy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PxXDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PyYDxxyy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PzZDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SXXS(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PxXXPx(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%PyXXPy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SXXDxxyy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SXXDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SYYDxxyy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%SZZDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DxyXXDxy(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DyzXXDyz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DxxyyXXDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DzzXXDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DxxyyYYDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DzzZZDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DxzXZDzz(geo%nSpecies), source=0.0_dp) + allocate(ctrl%mdftbAtomicIntegrals%DyzYZDxxyy(geo%nSpecies), source=0.0_dp) + + call getChild(value1, 'AtomDIntegralScalings', child2, requested=.false.) + if (associated(child2)) then + do iSp1 = 1, geo%nSpecies + call getChildValue(child2, trim(geo%speciesNames(iSp1)),& + & ctrl%mdftbAtomicIntegrals%DScaling(iSp1), 1.0_dp) + end do + end if + + call getChild(value1, 'AtomQIntegralScalings', child2, requested=.false.) + if (associated(child2)) then + do iSp1 = 1, geo%nSpecies + call getChildValue(child2, trim(geo%speciesNames(iSp1)),& + & ctrl%mdftbAtomicIntegrals%QScaling(iSp1), 1.0_dp) + end do + end if + + call getChild(value1, 'OneCenterAtomIntegrals', child2, requested=.true.) + do iSp1 = 1, geo%nSpecies + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_X_Px",& + & ctrl%mdftbAtomicIntegrals%SXPx(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_X_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%PxXDxxyy (iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_X_Dzz",& + & ctrl%mdftbAtomicIntegrals%PxXDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Py_Y_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%PyYDxxyy(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Pz_Z_Dzz",& + & ctrl%mdftbAtomicIntegrals%PzZDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_S",& + & ctrl%mdftbAtomicIntegrals%SXXS(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Px_XX_Px",& + & ctrl%mdftbAtomicIntegrals%PxXXPx(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Py_XX_Py",& + & ctrl%mdftbAtomicIntegrals%PyXXPy(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%SXXDxxyy(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_XX_Dzz",& + & ctrl%mdftbAtomicIntegrals%SXXDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_YY_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%SYYDxxyy(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":S_ZZ_Dzz",& + & ctrl%mdftbAtomicIntegrals%SZZDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxy_XX_Dxy",& + & ctrl%mdftbAtomicIntegrals%DxyXXDxy(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dyz_XX_Dyz",& + & ctrl%mdftbAtomicIntegrals%DyzXXDyz(iSp1), 0.0_dp) + !call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxx-yy_XX_Dzz",& + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_XX_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%DxxyyXXDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_XX_Dzz",& + & ctrl%mdftbAtomicIntegrals%DzzXXDzz(iSp1), 0.0_dp) + !call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxx-yy_YY_Dzz",& + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_YY_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%DxxyyYYDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dzz_ZZ_Dzz",& + & ctrl%mdftbAtomicIntegrals%DzzZZDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dxz_XZ_Dzz",& + & ctrl%mdftbAtomicIntegrals%DxzXZDzz(iSp1), 0.0_dp) + call getChildValue(child2, trim(geo%speciesNames(iSp1))//":Dyz_YZ_Dxx-yy",& + & ctrl%mdftbAtomicIntegrals%DyzYZDxxyy(iSp1), 0.0_dp) + end do + case default + call detailedError(child,"Unknown functions :"// char(buffer)) + end select + end if + end if + + end subroutine readMdftb !> Spin calculation #:if WITH_TRANSPORT diff --git a/src/dftbp/timedep/timeprop.F90 b/src/dftbp/timedep/timeprop.F90 index 29ba65102e..140b630761 100644 --- a/src/dftbp/timedep/timeprop.F90 +++ b/src/dftbp/timedep/timeprop.F90 @@ -32,7 +32,7 @@ module dftbp_timedep_timeprop use dftbp_dftb_getenergies, only : calcEnergies, calcDispersionEnergy, sumEnergies use dftbp_dftb_hamiltonian, only : TRefExtPot, resetExternalPotentials, resetInternalPotentials,& & addBlockChargePotentials, addChargePotentials, getSccHamiltonian - use dftbp_dftb_multiexpan, only : TDftbMultiExpan + use dftbp_dftb_multiexpan, only : TMdftb use dftbp_dftb_hybridxc, only : THybridXcFunc use dftbp_dftb_nonscc, only : TNonSccDiff, buildH0, buildS use dftbp_dftb_onsitecorrection, only : addOnsShift @@ -1472,7 +1472,7 @@ subroutine updateH(this, H1, ints, H0, speciesAll, qq, q0, coord, orb, potential integer :: iAtom, iEatom, iSpin, iKS, iK logical :: tImHam !! Multipole expansion - type(TDftbMultiExpan), allocatable :: dftbMultiExpan + type(TMdftb), allocatable :: dftbMultiExpan allocate(T2(this%nOrbs,this%nOrbs)) @@ -2015,7 +2015,7 @@ subroutine getTDEnergy(this, env, energy, rhoPrim, rho, neighbourList, nNeighbou type(TReksCalc), allocatable :: reks ! never allocated !! Multipole expansion - type(TDftbMultiExpan), allocatable :: dftbMultiExpan + type(TMdftb), allocatable :: dftbMultiExpan ! if Forces are calculated, rhoPrim has already been calculated ! check allways that calcEnergy is called AFTER getForces diff --git a/test/app/dftb+/mdftb/2H2O/dftb_in.hsd b/test/app/dftb+/mdftb/2H2O/dftb_in.hsd index 5b82dfe77e..f07a358708 100644 --- a/test/app/dftb+/mdftb/2H2O/dftb_in.hsd +++ b/test/app/dftb+/mdftb/2H2O/dftb_in.hsd @@ -30,15 +30,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -55,12 +55,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/32H2O/dftb_in.hsd b/test/app/dftb+/mdftb/32H2O/dftb_in.hsd index ff3cbe31f2..e8e414e592 100644 --- a/test/app/dftb+/mdftb/32H2O/dftb_in.hsd +++ b/test/app/dftb+/mdftb/32H2O/dftb_in.hsd @@ -120,15 +120,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -145,12 +145,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine/dftb_in.hsd index 26d2d15f10..203f2642f2 100644 --- a/test/app/dftb+/mdftb/Glycine/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine/dftb_in.hsd @@ -38,19 +38,19 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +75,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_Ordering/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_Ordering/dftb_in.hsd index 05c7ab6fdb..732eedbdad 100644 --- a/test/app/dftb+/mdftb/Glycine_Ordering/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_Ordering/dftb_in.hsd @@ -38,19 +38,19 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +75,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_RotTran1/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_RotTran1/dftb_in.hsd index 582bf197cf..33db191ca3 100644 --- a/test/app/dftb+/mdftb/Glycine_RotTran1/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_RotTran1/dftb_in.hsd @@ -38,19 +38,19 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +75,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_RotTran2/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_RotTran2/dftb_in.hsd index b6a34aab06..446fa58e78 100644 --- a/test/app/dftb+/mdftb/Glycine_RotTran2/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_RotTran2/dftb_in.hsd @@ -38,19 +38,19 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +75,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_RotTran3/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_RotTran3/dftb_in.hsd index 2cecb0ab4d..30c77cf973 100644 --- a/test/app/dftb+/mdftb/Glycine_RotTran3/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_RotTran3/dftb_in.hsd @@ -38,19 +38,19 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +75,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_mDFTBoff1/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_mDFTBoff1/dftb_in.hsd index a32b6fd629..a5cb1dd66a 100644 --- a/test/app/dftb+/mdftb/Glycine_mDFTBoff1/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_mDFTBoff1/dftb_in.hsd @@ -38,19 +38,22 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + # Disable all atomic dipole and quadrupole moments + # to ensure that mDFTB does not interfere with other parts of the code. + # The results should be identical to those of DFTB2. + AtomDIntegralScalings = { H = 0.0 C = 0.0 N = 0.0 O = 0.0 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 0.0 N = 0.0 O = 0.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 C:S_X_Px = 0.79621 C:S_XX_S = 0.76689 @@ -75,12 +78,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_mDFTBoff2/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_mDFTBoff2/dftb_in.hsd index 139f089cc6..9f2e3040d5 100644 --- a/test/app/dftb+/mdftb/Glycine_mDFTBoff2/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_mDFTBoff2/dftb_in.hsd @@ -38,19 +38,22 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + # Disable all atomic dipole and quadrupole moments + # to ensure that mDFTB does not interfere with other parts of the code. + # The results should be identical to those of DFTB2. + AtomDIntegralScalings = { H = 0.0 C = 0.8 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 2.6 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.0 C:S_X_Px = 0.0 C:S_XX_S = 0.0 @@ -75,12 +78,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/Glycine_mDFTBoff3/dftb_in.hsd b/test/app/dftb+/mdftb/Glycine_mDFTBoff3/dftb_in.hsd index 25fd9475db..b9a3d680ea 100644 --- a/test/app/dftb+/mdftb/Glycine_mDFTBoff3/dftb_in.hsd +++ b/test/app/dftb+/mdftb/Glycine_mDFTBoff3/dftb_in.hsd @@ -38,19 +38,22 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + # Disable all atomic dipole and quadrupole moments + # to ensure that mDFTB does not interfere with other parts of the code. + # The results should be identical to those of DFTB2. + AtomDIntegralScalings = { H = 0.0 C = 0.0 N = 0.0 O = 0.0 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 C = 0.0 N = 0.0 O = 0.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.0 C:S_X_Px = 0.0 C:S_XX_S = 0.0 @@ -75,12 +78,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/NH4+H2O/dftb_in.hsd b/test/app/dftb+/mdftb/NH4+H2O/dftb_in.hsd index dc94c52b3a..9a667c239d 100644 --- a/test/app/dftb+/mdftb/NH4+H2O/dftb_in.hsd +++ b/test/app/dftb+/mdftb/NH4+H2O/dftb_in.hsd @@ -34,17 +34,17 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 N:S_X_Px = 0.68857 N:S_XX_S = 0.57638 @@ -65,12 +65,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/NH4+H2O_RotTran1/dftb_in.hsd b/test/app/dftb+/mdftb/NH4+H2O_RotTran1/dftb_in.hsd index f56db4433f..6125414d65 100644 --- a/test/app/dftb+/mdftb/NH4+H2O_RotTran1/dftb_in.hsd +++ b/test/app/dftb+/mdftb/NH4+H2O_RotTran1/dftb_in.hsd @@ -34,17 +34,17 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 N:S_X_Px = 0.68857 N:S_XX_S = 0.57638 @@ -65,12 +65,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/NH4+H2O_RotTran2/dftb_in.hsd b/test/app/dftb+/mdftb/NH4+H2O_RotTran2/dftb_in.hsd index 7194c28091..d87ebe8ecd 100644 --- a/test/app/dftb+/mdftb/NH4+H2O_RotTran2/dftb_in.hsd +++ b/test/app/dftb+/mdftb/NH4+H2O_RotTran2/dftb_in.hsd @@ -34,17 +34,17 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 N:S_X_Px = 0.68857 N:S_XX_S = 0.57638 @@ -65,12 +65,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/NH4+H2O_RotTran3/dftb_in.hsd b/test/app/dftb+/mdftb/NH4+H2O_RotTran3/dftb_in.hsd index 82ad873c17..8a45e9c1f4 100644 --- a/test/app/dftb+/mdftb/NH4+H2O_RotTran3/dftb_in.hsd +++ b/test/app/dftb+/mdftb/NH4+H2O_RotTran3/dftb_in.hsd @@ -34,17 +34,17 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 N = 0.6 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 N = 3.4 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 N:S_X_Px = 0.68857 N:S_XX_S = 0.57638 @@ -65,12 +65,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/OH-H2O/dftb_in.hsd b/test/app/dftb+/mdftb/OH-H2O/dftb_in.hsd index 91d40e53fa..3546024e2e 100644 --- a/test/app/dftb+/mdftb/OH-H2O/dftb_in.hsd +++ b/test/app/dftb+/mdftb/OH-H2O/dftb_in.hsd @@ -29,15 +29,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -54,12 +54,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/OH-H2O_RotTran1/dftb_in.hsd b/test/app/dftb+/mdftb/OH-H2O_RotTran1/dftb_in.hsd index 2f5e8aea64..afc2486346 100644 --- a/test/app/dftb+/mdftb/OH-H2O_RotTran1/dftb_in.hsd +++ b/test/app/dftb+/mdftb/OH-H2O_RotTran1/dftb_in.hsd @@ -29,15 +29,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -54,12 +54,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/OH-H2O_RotTran2/dftb_in.hsd b/test/app/dftb+/mdftb/OH-H2O_RotTran2/dftb_in.hsd index 7670d1d08a..a021b3cce3 100644 --- a/test/app/dftb+/mdftb/OH-H2O_RotTran2/dftb_in.hsd +++ b/test/app/dftb+/mdftb/OH-H2O_RotTran2/dftb_in.hsd @@ -29,15 +29,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -54,12 +54,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/mdftb/OH-H2O_RotTran3/dftb_in.hsd b/test/app/dftb+/mdftb/OH-H2O_RotTran3/dftb_in.hsd index 550b01e8a2..654dc145fb 100644 --- a/test/app/dftb+/mdftb/OH-H2O_RotTran3/dftb_in.hsd +++ b/test/app/dftb+/mdftb/OH-H2O_RotTran3/dftb_in.hsd @@ -29,15 +29,15 @@ Hamiltonian = DFTB { } Multipole = OneCenterApproximation { - AtomDIntergralScalings = { + AtomDIntegralScalings = { H = 0.0 O = 0.2 } - AtomQIntergralScalings = { + AtomQIntegralScalings = { H = 0.0 O = 3.0 } - OneCenterAtomIntergrals = { + OneCenterAtomIntegrals = { H:S_XX_S = 0.55976 O:S_X_Px = 0.60894 O:S_XX_S = 0.45215 @@ -54,12 +54,12 @@ Hamiltonian = DFTB { } ParserOptions = { - ParserVersion = 8 - IgnoreUnprocessedNodes = Yes + ParserVersion = 14 + IgnoreUnprocessedNodes = No } Analysis = { - CalculateForces = Yes + PrintForces = Yes } Options = { diff --git a/test/app/dftb+/tests b/test/app/dftb+/tests index c089d26ca4..750a45f633 100644 --- a/test/app/dftb+/tests +++ b/test/app/dftb+/tests @@ -200,8 +200,8 @@ reks/Thymine_single #? not WITH_MPI reks/Thymine_T1 #? not WITH_MPI reks/Thymine_T1_rangesep #? not WITH_MPI -mdftb/2H2O #? not WITH_MPI -mdftb/32H2O #? not WITH_MPI +mdftb/2H2O #? OMP_THREADS <= 2 and not WITH_MPI +mdftb/32H2O #? MPI_PROCS <= 2 mdftb/Glycine #? not WITH_MPI mdftb/Glycine_mDFTBoff1 #? not WITH_MPI mdftb/Glycine_mDFTBoff2 #? not WITH_MPI