mapping between tblite and xtb (#1026)

* add wfn && results print for debug

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* map tblite wfn && results to xtb

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* dipole moments tblite vs ptb

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

---------

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

src/ptb/calculator.F90

@@ -36,6 +36,8 @@ module xtb_ptb_calculator
    use mctc_io_convert, only: autoev
 
 #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 tblite_basis_type, only: basis_type
@@ -290,29 +292,8 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
             return
          end if
       end if
-
-      call chk%wfn%allocate(mol%n, self%bas%nsh, self%bas%nao)
-      chk%wfn%n = mctcmol%nat
-      chk%wfn%nel = nint(chk%tblite%nocc)
-      chk%wfn%nopen = mctcmol%uhf
-      chk%wfn%nshell = self%bas%nsh
-      chk%wfn%nao = self%bas%nao
-      chk%wfn%P = chk%tblite%density(:, :, 1)
-      chk%wfn%q = chk%tblite%qat(:, 1)
-      chk%wfn%qsh = chk%tblite%qsh(:, 1)
-      chk%wfn%focca = chk%tblite%focc(:, 1)
-      chk%wfn%foccb = 0.0_wp
-      chk%wfn%focc(:) = chk%tblite%focc(:, 1)
-      chk%wfn%emo = chk%tblite%emo(:, 1) * autoev
-      chk%wfn%C = chk%tblite%coeff(:, :, 1)
-      chk%wfn%ihomo = chk%tblite%homo(1)
-      chk%wfn%ihomoa = chk%tblite%homo(1)
-      chk%wfn%ihomob = chk%tblite%homo(2)
-      chk%wfn%wbo = wbo(:, :, 1)
-      chk%wfn%dipm = chk%tblite%dpat(:, :, 1)
-      chk%wfn%qp = chk%tblite%qpat(:, :, 1)
-
-      results%hl_gap = (chk%tblite%emo(chk%tblite%homo(1) + 1, 1) - chk%tblite%emo(chk%tblite%homo(1), 1)) * autoev
+      call convert_tblite_to_wfn(env, self%bas, mol, chk, wbo=wbo)
+      call convert_tblite_to_results(results,mol,chk,energy,.true.)
       hlgap = results%hl_gap
 
       !> polarizability by simple perturbative treatment

src/tblite/CMakeLists.txt

@@ -18,6 +18,7 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 list(APPEND srcs
   "${dir}/calculator.F90"
   "${dir}/restart.F90"
+  "${dir}/mapping.F90"
 )
 
 set(srcs ${srcs} PARENT_SCOPE)

src/tblite/calculator.F90

@@ -40,8 +40,10 @@ module xtb_tblite_calculator
    use tblite_xtb_gfn1, only : new_gfn1_calculator, export_gfn1_param
    use tblite_xtb_ipea1, only : new_ipea1_calculator, export_ipea1_param
    use tblite_xtb_singlepoint, only : xtb_singlepoint
-   use tblite_data_spin, only : get_spin_constant
+   use tblite_data_spin, only : get_spin_constant 
+   use xtb_tblite_mapping, only : convert_tblite_to_wfn
 #endif
+   use xtb_tblite_mapping, only : convert_tblite_to_results
    use xtb_mctc_accuracy, only : wp
    use xtb_type_calculator, only : TCalculator
    use xtb_type_data
@@ -108,7 +110,7 @@ module xtb_tblite_calculator
 contains
 
 
-!> Create a new instance of an tblite based xTB calculator
+!> Create a new instance of tblite based xTB calculator
 subroutine newTBLiteCalculator(env, mol, calc, input)
    !> Source of the generated errors
    character(len=*), parameter :: source = 'tblite_calculator_newTBLiteCalculator'
@@ -319,17 +321,13 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
       return
    end if
 
-   results%e_total = energy
-   results%converged = .true.
-   results%dipole = sum(chk%tblite%dpat(:, :, 1), 2) + matmul(struc%xyz, chk%tblite%qat(: ,1))
-   results%gnorm = norm2(gradient)
+   ! convert tblite results into xtb data !
+   call convert_tblite_to_wfn(env, self%tblite%bas, mol, chk)
+   call convert_tblite_to_results(results,mol,chk,energy,.true.,gradient=gradient)
+   hlgap = results%hl_gap
 
-!   if (printlevel > 2) then
-!      call ascii_levels(ctx%unit, printlevel, chk%tblite%homo, chk%tblite%emo, &
-!         & chk%tblite%focc, 7)
-!   end if
 #else
-    call feature_not_implemented(env)
+   call feature_not_implemented(env)
 #endif
 end subroutine singlepoint
 

src/tblite/mapping.F90

@@ -0,0 +1,161 @@
+#ifndef WITH_TBLITE
+#define WITH_TBLITE 0
+#endif
+
+module xtb_tblite_mapping
+
+   use xtb_type_environment, only: TEnvironment
+   use xtb_type_restart, only: TRestart
+   use xtb_type_data, only: scc_results
+   use xtb_type_molecule, only: TMolecule
+   use xtb_mctc_accuracy, only: wp   
+   use mctc_io_convert, only: autoev
+   use xtb_type_wavefunction, only: TWavefunction
+
+#if WITH_TBLITE
+   use tblite_basis_type, only: basis_type
+#endif
+   implicit none
+   private
+   public :: convert_tblite_to_results
+
+#if WITH_TBLITE
+   public :: convert_tblite_to_wfn
+   interface assignment(=)
+      module procedure :: from_tblite_wfn
+      module procedure :: from_tblite_basis
+      module procedure :: from_struc
+   end interface assignment(=)
+#endif
+contains
+
+#if WITH_TBLITE
+!> convert tblite wavefunction to xtb wavefunction
+subroutine convert_tblite_to_wfn(env, bas, mol, chk, wbo)
+
+   !> computational environment
+   type(TEnvironment), intent(inout) :: env
+
+   !> basis set
+   type(basis_type), intent(in) :: bas
+
+   !> molecular structure
+   type(TMolecule), intent(in) :: mol
+
+   !> restart data
+   type(TRestart), intent(inout) :: chk
+
+   !> wiberg bond orders
+   real(wp), optional, intent(in) :: wbo(:,:,:)
+
+   call chk%wfn%allocate(mol%n, bas%nsh, bas%nao)
+
+   ! assignment interfaces !
+   chk%wfn = chk%tblite
+   chk%wfn = bas
+   chk%wfn = mol
+
+   if (present(wbo)) chk%wfn%wbo = wbo(:,:,1) ! only 
+
+end subroutine convert_tblite_to_wfn
+#endif
+
+!> convert tblite results to xtb results
+subroutine convert_tblite_to_results(results, mol, chk, energy, converged, gradient) 
+
+   !> scc results
+   type(scc_results), intent(inout) :: results
+
+   !> molecular structure
+   type(TMolecule), intent(in) :: mol
+
+   !> restart data
+   type(TRestart), intent(in) :: chk
+
+   !> SCC energy
+   real(wp), intent(in) :: energy
+
+   !> convergence flag
+   logical, intent(in) :: converged
+
+   !> (analytical) gradients
+   real(wp), optional, intent(in) :: gradient(:,:)
+
+
+   results%e_total = energy
+   results%converged = converged
+
+#if WITH_TBLITE
+
+   ! do not overwrite dipole moments, if calculated (PTB case) ! 
+   if (all(results%dipole == 0.0_wp)) then
+      results%dipole = sum(chk%tblite%dpat(:, :, 1), 2) + matmul(mol%xyz, chk%tblite%qat(: ,1))
+   endif
+   results%hl_gap = (chk%tblite%emo(chk%tblite%homo(1) + 1, 1) - chk%tblite%emo(chk%tblite%homo(1), 1)) * autoev
+
+#endif
+
+   if (present(gradient)) results%gnorm = norm2(gradient)
+
+end subroutine convert_tblite_to_results
+
+
+#if WITH_TBLITE
+
+subroutine from_tblite_wfn(wfn, tblite)
+
+   use tblite_wavefunction, only : wavefunction_type
+
+   type(TWavefunction), intent(inout) :: wfn 
+   type(wavefunction_type), intent(in) :: tblite
+
+   wfn%dipm = tblite%dpat(:, :, 1)
+   wfn%nel = nint(tblite%nocc)
+   wfn%P = tblite%density(:, :, 1)
+   wfn%q = tblite%qat(:, 1)
+   wfn%qsh = tblite%qsh(:, 1)
+   wfn%focca = tblite%focc(:, 1)
+   wfn%foccb = 0.0_wp
+   wfn%focc(:) = tblite%focc(:, 1)
+   wfn%emo = tblite%emo(:, 1) * autoev
+   wfn%C = tblite%coeff(:, :, 1)
+   wfn%ihomo = tblite%homo(1)
+   wfn%ihomoa = tblite%homo(1)
+   wfn%ihomob = tblite%homo(2)
+   wfn%qp = tblite%qpat(:, :, 1)
+
+end subroutine from_tblite_wfn
+
+subroutine from_tblite_basis(wfn, bas)
+
+   use tblite_basis_type, only: basis_type
+
+   type(TWavefunction), intent(inout) :: wfn
+   type(basis_type), intent(in) :: bas
+
+   wfn%nshell = bas%nsh
+   wfn%nao = bas%nao
+
+end subroutine from_tblite_basis
+
+subroutine from_struc(wfn, mol)
+
+   type(TWavefunction), intent(inout) :: wfn
+   type(TMolecule), intent(in) :: mol
+
+   wfn%n = mol%n
+   wfn%nopen = mol%uhf
+
+end subroutine from_struc
+#endif
+
+#if ! WITH_TBLITE
+subroutine feature_not_implemented(env)
+   !> Computational environment
+   type(TEnvironment), intent(inout) :: env
+
+   call env%error("Compiled without support for tblite library")
+end subroutine feature_not_implemented
+#endif
+
+end module xtb_tblite_mapping

src/tblite/meson.build

@@ -16,5 +16,6 @@
 srcs += files(
   'calculator.F90',
   'restart.F90',
+  'mapping.F90', 
 )
 

src/type/data.f90

@@ -88,6 +88,8 @@ module xtb_type_data
       real(wp) :: e_batm = 0.0_wp
       real(wp) :: e_ext = 0.0_wp
       type(TIFFResults), allocatable :: iff_results
+   contains 
+      procedure :: print => print_scc_results
    end type scc_results
 
    type freq_results
@@ -172,4 +174,30 @@ subroutine deallocate_freq_results(self)
    if (allocated( self%pg   )) deallocate( self%pg   )
 end subroutine deallocate_freq_results
 
+!> print SCC results (for debug)
+subroutine print_scc_results(self, unit)
+   use iso_fortran_env, only : output_unit
+   class(scc_results) :: self
+   integer, optional :: unit
+
+   integer :: out ! output unit holder
+   character(len=*), parameter :: dfmt = '(3x,a,2x,f12.6)' ! format for double precision
+   integer :: i
+
+   if (present(unit)) then
+      out = unit
+   else
+      out = output_unit
+   endif
+
+   write(out, '(3x,a,/,2x,11("="))') 'SCC Results'
+
+   write(out, dfmt) 'e_total   :', self%e_total
+   write(out, dfmt) 'hl_gap    :', self%hl_gap
+   do i=1,3
+      write(out, dfmt) 'dipole    :', self%dipole(i)
+   enddo
+
+end subroutine print_scc_results
+
 end module xtb_type_data