From 96e6a2547fb50b30859e52416505e096a037d209 Mon Sep 17 00:00:00 2001 From: albert <92109627+Albkat@users.noreply.github.com> Date: Tue, 14 May 2024 13:25:45 +0200 Subject: [PATCH] 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 | 27 +------ src/tblite/CMakeLists.txt | 1 + src/tblite/calculator.F90 | 20 +++-- src/tblite/mapping.F90 | 161 ++++++++++++++++++++++++++++++++++++++ src/tblite/meson.build | 1 + src/type/data.f90 | 28 +++++++ src/type/wavefunction.f90 | 138 ++++++++++++++++++++++++++++---- 7 files changed, 325 insertions(+), 51 deletions(-) create mode 100644 src/tblite/mapping.F90 diff --git a/src/ptb/calculator.F90 b/src/ptb/calculator.F90 index c6538ebcb..1c3bd9a0e 100644 --- a/src/ptb/calculator.F90 +++ b/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 diff --git a/src/tblite/CMakeLists.txt b/src/tblite/CMakeLists.txt index 7353809b6..8b2db7176 100644 --- a/src/tblite/CMakeLists.txt +++ b/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) diff --git a/src/tblite/calculator.F90 b/src/tblite/calculator.F90 index 7ad8542fb..653e6b285 100644 --- a/src/tblite/calculator.F90 +++ b/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 diff --git a/src/tblite/mapping.F90 b/src/tblite/mapping.F90 new file mode 100644 index 000000000..9431dd2cf --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/tblite/meson.build b/src/tblite/meson.build index d9aeabe4c..56529ca06 100644 --- a/src/tblite/meson.build +++ b/src/tblite/meson.build @@ -16,5 +16,6 @@ srcs += files( 'calculator.F90', 'restart.F90', + 'mapping.F90', ) diff --git a/src/type/data.f90 b/src/type/data.f90 index 2da7963d2..59bd6a196 100644 --- a/src/type/data.f90 +++ b/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 diff --git a/src/type/wavefunction.f90 b/src/type/wavefunction.f90 index 71712d450..6b6788247 100644 --- a/src/type/wavefunction.f90 +++ b/src/type/wavefunction.f90 @@ -24,45 +24,65 @@ module xtb_type_wavefunction private type :: TWavefunction + + !> number of atoms integer :: n = 0 - !! Number of atoms + + !> number of electrons integer :: nel = 0 - !! Number of elctrons + + !> number of unpaired electrons integer :: nopen = 0 - !! Number of unpaired electrons + + !> number of atomic orbitals integer :: nao = 0 - !! Number of atomic orbitals + + !> number of shells integer :: nshell = 0 - !! Number of shells + + !> density matrix real(wp),allocatable :: P(:,:) - !! Density matrix + + !> partial charges real(wp),allocatable :: q(:) - !! Partial charges + + !> shell charges real(wp),allocatable :: qsh(:) - !! Shell charges + + !> dipole moments real(wp),allocatable :: dipm(:,:) - !! Dipole moments + + !> quadrupole moments real(wp),allocatable :: qp(:,:) - !! Quadrupole moments + + !> wiberg bond orders real(wp),allocatable :: wbo(:,:) - !! Wiberg bond orders + + !> HOMO position integer :: ihomo = 0,ihomoa = 0,ihomob = 0 - !! HOMO position + + !> fermi real(wp) :: efa = 0.0_wp, efb = 0.0_wp + + !> fractional occupation real(wp),allocatable :: focc(:) - !! Fractional occupation + + !> alpha space real(wp),allocatable :: focca(:) - !! For alpha space + + !> beta space real(wp),allocatable :: foccb(:) - !! For beta space + + !> orbital energies real(wp),allocatable :: emo(:) - !! Orbital energies + + !> molecular orbitals real(wp),allocatable :: C(:,:) - !! Molecular orbitals contains procedure :: allocate => allocate_wavefunction procedure :: deallocate => deallocate_wavefunction + procedure :: print => print_wavefunction end type TWavefunction contains @@ -105,4 +125,88 @@ subroutine deallocate_wavefunction(self) if(allocated(self%C)) deallocate(self%C) end subroutine deallocate_wavefunction +!> print content of wavefunction (for debug) +subroutine print_wavefunction(self, level, unit) + use iso_fortran_env, only : output_unit + class(TWavefunction),intent(in) :: self + integer,intent(in) :: level + integer,intent(in),optional :: unit + + character(len=*), parameter :: fmt0 = '(3x,a,3x,i0)' + character(len=30):: fmt1_n + character(len=30):: fmt1_nao + character(len=30):: fmt1_nshell + integer :: out, i, ndim + + if (self%nao>20) then + ndim=20 + else + ndim=self%nao + end if + + write(fmt1_n,'(a,i0,a)') '(3x,a,3x,',self%n,'(f6.2,1x))' + write(fmt1_nao,'(a,i0,a)') '(3x,a,3x,',ndim,'(f6.2,1x))' + write(fmt1_nshell,'(a,i0,a)') '(3x,a,3x,',self%nshell,'(f6.2,1x))' + + if (present(unit)) then + out = unit + else + out = output_unit + end if + + + if (level >= 0 .and. level <= 2) then + write(out, '(3x,a,/,2x,12("="))') 'Wavefunction' + + ! scalar values ! + write(out,fmt0) 'n :', self%n + write(out,fmt0) 'nel :', self%nel + write(out,fmt0) 'nopen :', self%nopen + write(out,fmt0) 'nao :', self%nao + write(out,fmt0) 'nshell :', self%nshell + write(out,fmt0) 'ihomo :', self%ihomo + write(out,fmt0) 'ihomoa :', self%ihomoa + write(out,fmt0) 'ihomob :', self%ihomob + + ! vector values ! + if (level>0) then + write(out,fmt1_n) 'q :', self%q + write(out,fmt1_nshell) 'qsh :', self%qsh + write(out,fmt1_nao) 'focc :', self%focc(:ndim) + write(out,fmt1_nao) 'emo :', self%emo(:ndim) + + ! matrix values ! + if (level>1) then + + write(out,'(/)') + + do i=1,3 + write(out,fmt1_n) 'dipm :', self%dipm(i,:) + end do + + write(out,'(/)') + + do i=1,6 + write(out,fmt1_n) 'qp :', self%qp(i,:) + end do + + write(out,'(/)') + + do i=1,ndim + write(out,fmt1_nao) 'C :', self%C(:ndim,i) + end do + + write(out,'(/)') + + do i=1,ndim + write(out,fmt1_nao) 'P :', self%P(:ndim,i) + end do + + write(out,'(/)') + + endif + endif + endif +end subroutine print_wavefunction + end module xtb_type_wavefunction