Update dependencies

cmake/modules/Finddftd4.cmake

@@ -0,0 +1,33 @@
+# This file is part of xtb.
+# SPDX-Identifier: LGPL-3.0-or-later
+#
+# xtb is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# xtb is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with xtb.  If not, see <https://www.gnu.org/licenses/>.
+
+set(_lib "dftd4")
+set(_pkg "DFTD4")
+set(_url "https://github.com/dftd4/dftd4")
+set(_rev "v4.0.1")
+
+if(NOT DEFINED "${_pkg}_FIND_METHOD")
+   set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch")
+endif()
+
+include("${CMAKE_CURRENT_LIST_DIR}/xtb-utils.cmake")
+
+xtb_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}" "${_rev}")
+
+unset(_lib)
+unset(_pkg)
+unset(_url)
+unset(_rev)

cmake/modules/Findmctc-lib.cmake

@@ -17,7 +17,7 @@
 set(_lib "mctc-lib")
 set(_pkg "MCTCLIB")
 set(_url "https://github.com/grimme-lab/mctc-lib")
-set(_rev "v0.5.0")
+set(_rev "v0.5.1")
 
 if(NOT DEFINED "${_pkg}_FIND_METHOD")
    set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch")

cmake/modules/Findmulticharge.cmake

@@ -0,0 +1,33 @@
+# This file is part of xtb.
+# SPDX-Identifier: LGPL-3.0-or-later
+#
+# xtb is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# xtb is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with xtb.  If not, see <https://www.gnu.org/licenses/>.
+
+set(_lib "multicharge")
+set(_pkg "MULTICHARGE")
+set(_url "https://github.com/grimme-lab/multicharge")
+set(_rev "v0.5.0")
+
+if(NOT DEFINED "${_pkg}_FIND_METHOD")
+   set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch")
+endif()
+
+include("${CMAKE_CURRENT_LIST_DIR}/xtb-utils.cmake")
+
+xtb_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}" "${_rev}")
+
+unset(_lib)
+unset(_pkg)
+unset(_url)
+unset(_rev)

cmake/modules/Findtblite.cmake

@@ -17,7 +17,7 @@
 set(_lib "tblite")
 set(_pkg "TBLITE")
 set(_url "https://github.com/tblite/tblite")
-set(_rev "v0.4.0")
+set(_rev "HEAD")
 
 if(NOT DEFINED "${_pkg}_FIND_METHOD")
    set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch")

src/dipro.F90

@@ -105,7 +105,7 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
    type(wavefunction_type), allocatable :: wfx(:)
 
    !> Molecular gradient, strain derivatives
-   real(wp), allocatable :: gradient(:, :), sigma(:,:), nel(:)
+   real(wp), allocatable :: gradient(:, :), sigma(:,:), nel(:), homo(:,:)
    real(wp), allocatable :: overlap(:, :), trans(:, :), wbo(:, :), chrg(:), p2mat(:,:), coeff2(:,:),loc(:,:)
    real(wp), allocatable :: orbital(:, :, :), scmat(:, :), fdim(:, :), scratch(:), efrag(:),y(:,:),Edim(:,:)
    integer, allocatable :: spinfrag(:), start_index(:),end_index(:),orbprint(:)
@@ -206,7 +206,7 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
    nao = size(wfn%emo)
    allocate(orbital(nao, nfrag,nao), efrag(nfrag), scmat(nao, nao), fdim(nao, nao), &
       & scratch(nao), mfrag(nfrag), wfx(nfrag), chrg(nfrag), spinfrag(nfrag), &
-      & start_index(nfrag),end_index(nfrag),orbprint(nfrag),nel(nfrag))
+      & start_index(nfrag), end_index(nfrag), orbprint(nfrag), nel(nfrag), homo(nfrag, 2))
 
 !==================================external files CHRG & UHF read-in====================================
 
@@ -270,7 +270,7 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
       call new_wavefunction(wfx(ifr), mfrag(ifr)%nat, fcalc(ifr)%bas%nsh, fcalc(ifr)%bas%nao, &
          & 1, set%etemp * ktoau)
 
-     !> mol%type (dimer) == mfrag%type (fragments), wfn (dimer) == wfx (fragments), calc (dimer)==fcalc(fragments)
+      !> mol%type (dimer) == mfrag%type (fragments), wfn (dimer) == wfx (fragments), calc (dimer)==fcalc(fragments)
       wfx%nspin=1
       call xtb_singlepoint(ctx, mfrag(ifr), fcalc(ifr), wfx(ifr), tblite%accuracy, energy)
       if (ctx%failed()) then
@@ -280,6 +280,14 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
 
       nel(ifr)=wfx(ifr)%nel(1)+wfx(ifr)%nel(2)
 
+      ! Find HOMO index
+      homo(ifr,1) = merge(wfx(ifr)%nel(1)+1, wfx(ifr)%nel(1), mod(wfx(ifr)%nel(1), 1.0_wp) > 0.5_wp)
+      if (size(wfx(ifr)%nel, 1) == 2) then
+         homo(ifr,2) = merge(wfx(ifr)%nel(2)+1, wfx(ifr)%nel(2), mod(wfx(ifr)%nel(2), 1.0_wp) > 0.5_wp)
+      else
+         homo(ifr,2) = homo(ifr,1)
+      end if
+
 !==================================DIPRO==================================================
 
       do j = 1, fcalc(ifr)%bas%nao
@@ -313,8 +321,8 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
 
    do ifr=1,nfrag
       do j = 1, fcalc(ifr)%bas%nao
-         if (wfx(ifr)%emo(j,1) .ge. (wfx(ifr)%emo(wfx(ifr)%homo(2),1) - dipro%othr/autoev) .and.&
-           & wfx(ifr)%emo(j,1) .le. (wfx(ifr)%emo(wfx(ifr)%homo(2)+1,1) + dipro%othr/autoev)) then
+         if (wfx(ifr)%emo(j,1) .ge. (wfx(ifr)%emo(homo(ifr,2),1) - dipro%othr/autoev) .and.&
+           & wfx(ifr)%emo(j,1) .le. (wfx(ifr)%emo(homo(ifr,2)+1,1) + dipro%othr/autoev)) then
            if (start_index(ifr).eq.-1) then 
                    start_index(ifr) = j
               end if
@@ -343,7 +351,7 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
    !> scmat=S_dim*C_dim
    call gemm(overlap, coeff2, scmat)
    do j = start_index(1), end_index(1)
-      orbprint(1)=wfx(1)%homo(max(2,1))-j
+      orbprint(1)=homo(1,2)-j
 
       y(:,1)=0
       !> gemv(amat, xvec,yvec,a1,a2,transa): X=a1*Amat*xvec+a2*yvec
@@ -355,7 +363,7 @@ subroutine get_jab(env, tblite, mol, fragment, dipro)
       efrag(1)=dot( y(:,1), scratch)
 
       do k = start_index(2), end_index(2)
-         orbprint(2)=wfx(2)%homo(max(2,1))-k
+         orbprint(2)=homo(2,2)-k
 
          y(:,2)=0
          !> y_mon2(ifr)=C_mon2(k)*S_dim*C_dim

src/ptb/CMakeLists.txt

@@ -35,6 +35,7 @@ list(APPEND srcs
   "${dir}/plusu.F90"
   "${dir}/property.F90"
   "${dir}/response.F90"
+  "${dir}/solver.F90"
 )
 
 set(srcs ${srcs} PARENT_SCOPE)

src/ptb/calculator.F90

@@ -38,7 +38,7 @@ module xtb_ptb_calculator
 #if WITH_TBLITE
    use xtb_tblite_mapping, only : convert_tblite_to_wfn, convert_tblite_to_results
 
-   use multicharge_model, only: new_mchrg_model, mchrg_model_type
+   use multicharge_model, only: new_eeq_model, mchrg_model_type, eeq_model
 
    use tblite_basis_type, only: basis_type
    use tblite_context, only: context_type
@@ -71,7 +71,7 @@ module xtb_ptb_calculator
       type(basis_type) :: bas, cbas
 
       !> EEQ Model
-      type(mchrg_model_type) :: eeqmodel
+      class(mchrg_model_type), allocatable :: eeqmodel
 #endif
 
       !> Parametrisation data base
@@ -127,6 +127,8 @@ subroutine newPTBCalculator(env, struc, calc, accuracy)
 #if WITH_TBLITE
       !> mctc-io structure type
       type(structure_type) :: mol
+      type(eeq_model), allocatable :: tmp_eeqmodel
+      type(error_type), allocatable :: error
 
       mol = struc
 
@@ -164,8 +166,14 @@ subroutine newPTBCalculator(env, struc, calc, accuracy)
       call add_core_basis(mol, calc%ptbData%corepotential, calc%cbas)
 
       !> set up the EEQ model
-      call new_mchrg_model(calc%eeqmodel, chi=calc%ptbData%eeq%chi, &
-      & rad=calc%ptbData%eeq%alp, eta=calc%ptbData%eeq%gam, kcn=calc%ptbData%eeq%cnf)
+      allocate(tmp_eeqmodel)
+      call new_eeq_model(tmp_eeqmodel, mol, error, chi=calc%ptbData%eeq%chi, &
+      & rad=calc%ptbData%eeq%alp, eta=calc%ptbData%eeq%gam, kcnchi=calc%ptbData%eeq%cnf)
+      if (allocated(error)) then
+         call env%error('Could not initialize EEQ model', source)
+         return
+      end if
+      call move_alloc(tmp_eeqmodel, calc%eeqmodel)
 
       !> check for external point charge field
       !> not implemented yet

src/ptb/meson.build

@@ -33,4 +33,5 @@ srcs += files(
   'plusu.F90',
   'property.F90',
   'response.F90',
+  'solver.F90',
 )

src/ptb/ncoord.F90

@@ -26,7 +26,7 @@ module xtb_ptb_ncoord
    use mctc_env, only: wp
    use mctc_io, only: structure_type
 
-   use tblite_data_covrad, only: get_covalent_rad
+   use mctc_data, only: get_covalent_rad
 
    implicit none
 

src/ptb/response.F90

@@ -45,13 +45,13 @@ module xtb_ptb_response
    use xtb_ptb_data, only: TPTBData
    use xtb_ptb_integral_types, only: aux_integral_type
    use xtb_ptb_param, only: ptbGlobals
-   use xtb_ptb_scf, only: get_density
    use xtb_ptb_mmlpopanalysis, only: get_mml_shell_charges
    use xtb_ptb_hamiltonian, only: get_hamiltonian
    use xtb_ptb_guess, only: get_psh_from_qsh
    use xtb_ptb_paulixc, only: calc_Vxc_pauli
    use xtb_ptb_coulomb, only: coulomb_potential
    use xtb_ptb_plusu, only: plusu_potential_type
+   use xtb_ptb_solver, only : ptb_solver_type, new_ptb_solver
 
    implicit none
    private
@@ -106,10 +106,8 @@ subroutine numgrad_polarizability(ctx, data, mol, bas, wfn, ints, auxints, &
       type(container_cache) :: icache
       !> Electric field object
       type(electric_field) :: efield_object
-      !> Electronic solver
-      class(solver_type), allocatable :: solver
-      !> Electronic entropy
-      real(wp) :: ts
+      !> PTB electronic solver
+      class(ptb_solver_type), allocatable :: ptbsolver
       !> Molecular dipole moment
       real(wp) :: dip_plus(3), dip_minus(3)
       real(wp) :: eff_ef(3), tmp_ef(3)
@@ -121,8 +119,8 @@ subroutine numgrad_polarizability(ctx, data, mol, bas, wfn, ints, auxints, &
       logical, parameter :: debug(4) = &
                [ .false., .false., .false., .false. ]
 
-      !> Solver for the effective Hamiltonian
-      call ctx%new_solver(solver, bas%nao)
+      !> Solver for the effective Hamiltonian (specified for each diagonalization)
+      allocate(ptbsolver)
 
       alpha = 0.0_wp
       if (present(efield)) then
@@ -158,7 +156,8 @@ subroutine numgrad_polarizability(ctx, data, mol, bas, wfn, ints, auxints, &
          endif
 
          !> Solve effective Hamiltonian including the electric field
-         call get_density(wfn_tmp, solver, ints, ts, error, ptbGlobals%geps, ptbGlobals%geps0)
+         call new_ptb_solver(ptbsolver, wfn%nel, wfn%kt, ptbGlobals%geps, ptbGlobals%geps0)
+         call ptbsolver%get_density(wfn_tmp%coeff, ints%overlap, wfn_tmp%emo, wfn_tmp%focc, wfn_tmp%density, error)
          if (allocated(error)) then
             call ctx%set_error(error)
             return
@@ -209,7 +208,8 @@ subroutine numgrad_polarizability(ctx, data, mol, bas, wfn, ints, auxints, &
          endif
 
          !> Solve effective Hamiltonian including the electric field
-         call get_density(wfn_tmp, solver, ints, ts, error, ptbGlobals%geps, ptbGlobals%geps0)
+         call new_ptb_solver(ptbsolver, wfn%nel, wfn%kt, ptbGlobals%geps, ptbGlobals%geps0)
+         call ptbsolver%get_density(wfn_tmp%coeff, ints%overlap, wfn_tmp%emo, wfn_tmp%focc, wfn_tmp%density, error)
          if (allocated(error)) then
             call ctx%set_error(error)
             return
@@ -235,7 +235,8 @@ subroutine numgrad_polarizability(ctx, data, mol, bas, wfn, ints, auxints, &
             & ves_twostepscf, CN_plusU, efield_object, dip_plus)
          if (ctx%failed()) return
 
-         alpha(k, 1:3) = -(dip_minus - dip_plus) / (2.0_wp * delta)                               ! numerical diff. dmu/dfield
+         ! numerical diff. dmu/dfield
+         alpha(k, 1:3) = -(dip_minus - dip_plus) / (2.0_wp * delta)
       end do cart_coord_loop
 
       !> Symmetrization of polarizability tensor
@@ -278,10 +279,8 @@ subroutine onestepscf(ctx, data, mol, bas, wfn, ints, auxints, vecp, list, level
       type(electric_field), intent(in) :: efield
       !> Dipole moment
       real(wp), intent(out) :: dipole(3)
-      !> Electronic solver
-      class(solver_type), allocatable :: solver
-      !> Electronic entropy
-      real(wp) :: ts
+      !> PTB electronic solver
+      class(ptb_solver_type), allocatable :: ptbsolver
       !> Error type
       type(error_type), allocatable :: error
       !> Restart data for interaction containers
@@ -303,8 +302,8 @@ subroutine onestepscf(ctx, data, mol, bas, wfn, ints, auxints, vecp, list, level
       !> debug mode
       logical, parameter :: debug = .false.
 
-      !> Solver for the effective Hamiltonian
-      call ctx%new_solver(solver, bas%nao)
+      !> Solver for the effective Hamiltonian (specified for each diagonalization)
+      allocate(ptbsolver)
 
       !> Reset H0 matrix
       ints%hamiltonian = 0.0_wp
@@ -362,7 +361,8 @@ subroutine onestepscf(ctx, data, mol, bas, wfn, ints, auxints, vecp, list, level
          end do
       endif
 
-      call get_density(wfn, solver, ints, ts, error, ptbGlobals%geps, ptbGlobals%geps0)
+      call new_ptb_solver(ptbsolver, wfn%nel, wfn%kt, ptbGlobals%geps, ptbGlobals%geps0)
+      call ptbsolver%get_density(wfn%coeff, ints%overlap, wfn%emo, wfn%focc, wfn%density, error)
       if (allocated(error)) then
          call ctx%set_error(error)
          return