From a7e0c9f2a6d9072dea487ee4a1072a3c2f73ca6c Mon Sep 17 00:00:00 2001
From: albert <92109627+Albkat@users.noreply.github.com>
Date: Thu, 11 Apr 2024 19:11:16 +0200
Subject: [PATCH] add --ceh flag

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>
---
 man/xtb.1.adoc            |  6 ++++
 src/prog/main.F90         | 17 +++++++--
 src/tblite/calculator.F90 | 73 ++++++++++++++++++++++++++++++++++++++-
 src/xhelp.f90             |  3 ++
 4 files changed, 96 insertions(+), 3 deletions(-)

diff --git a/man/xtb.1.adoc b/man/xtb.1.adoc
index 7d88206cb..b829f3291 100644
--- a/man/xtb.1.adoc
+++ b/man/xtb.1.adoc
@@ -80,6 +80,9 @@ OPTIONS
 
 *--tblite* ::
     use tblite library as implementation for xTB
+    
+--ceh* ::
+    calculate CEH (Charge-Extended Hückel model) charges and write them to ceh.charges file
 
 *--ptb* ::
     performs single-point calculation with the density tight-binding method PTB.
@@ -346,6 +349,9 @@ OUTPUT
 *charges*::
    contains Mulliken partial charges calculated in SCC
 
+*ceh.charges*::
+   contains CEH (Charge-Extended Hückel) charges
+
 *wbo*::
    contains Wiberg bond order calculated in SCC, *--namespace* compatible.
 
diff --git a/src/prog/main.F90 b/src/prog/main.F90
index 96465feb3..4e5821882 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -93,7 +93,8 @@ module xtb_prog_main
    use xtb_iff_data, only: TIFFData
    use xtb_oniom, only: oniom_input, TOniomCalculator, calculateCharge
    use xtb_vertical, only: vfukui
-   use xtb_tblite_calculator, only: TTBLiteCalculator, TTBLiteInput, newTBLiteWavefunction
+   use xtb_tblite_calculator, only: TTBLiteCalculator, TTBLiteInput, &
+         & newTBLiteWavefunction, get_ceh
    use xtb_ptb_calculator, only: TPTBCalculator
    use xtb_solv_cpx, only: TCpcmx
    use xtb_dipro, only: get_jab, jab_input
@@ -165,7 +166,6 @@ subroutine xtbMain(env, argParser)
       real(wp), allocatable :: q(:)
       real(wp), allocatable :: ql(:)
       real(wp), allocatable :: qr(:)
-
 !! ------------------------------------------------------------------------
       integer, external :: ncore
 
@@ -617,6 +617,10 @@ subroutine xtbMain(env, argParser)
          call newTBLiteWavefunction(env, mol, calc, chk)
       end select
 
+      ! get CEH charges !
+      if (tblite%ceh) &
+         call get_ceh(env,mol,tblite)
+
       ! ------------------------------------------------------------------------
       !> printout a header for the exttyp
       call calc%writeInfo(env%unit, mol)
@@ -1453,6 +1457,15 @@ subroutine parseArguments(env, args, inputFile, paramFile, lgrad, &
             else
                call env%error("Molecular charge is not provided", source)
             end if
+         
+         case('--ceh')
+            if (get_xtb_feature('tblite')) then
+               tblite%ceh = .true.
+            else
+               call env%error("CEH charges are only available through tblite library", source)
+               return
+            end if
+
 
          case ('-u', '--uhf')
             call args%nextArg(sec)
diff --git a/src/tblite/calculator.F90 b/src/tblite/calculator.F90
index ff45b96a3..7ad8542fb 100644
--- a/src/tblite/calculator.F90
+++ b/src/tblite/calculator.F90
@@ -26,6 +26,7 @@ module xtb_tblite_calculator
    use mctc_io, only : structure_type, read_structure, filetype
 #if WITH_TBLITE
    use tblite_basis_type, only : basis_type
+   use tblite_ceh_ceh, only : new_ceh_calculator
    use tblite_container, only : container_type
    use tblite_context, only : context_type, context_terminal, escape
    use tblite_external_field, only : electric_field
@@ -57,6 +58,7 @@ module xtb_tblite_calculator
    private
 
    public :: TTBLiteCalculator, TTBLiteInput, newTBLiteCalculator, newTBLiteWavefunction
+   public :: get_ceh
 
    !> Input for tblite library
    type :: TTBLiteInput
@@ -72,6 +74,8 @@ module xtb_tblite_calculator
       character(len=:), allocatable :: method
       !> Colorful output
       logical :: color = .false.
+      !> CEH charges
+      logical :: ceh = .false.
    end type TTBLiteInput
 
    !> Calculator interface for xTB based methods
@@ -149,6 +153,11 @@ subroutine newTBLiteCalculator(env, mol, calc, input)
          call new_gfn1_calculator(calc%tblite, struc)
       case("ipea1")
          call new_ipea1_calculator(calc%tblite, struc)
+      case("ceh")
+         calc%guess = method
+         calc%nspin = 1
+         calc%etemp = 5000.0_wp * kt
+         call new_ceh_calculator(calc%tblite, struc)
       end select
    end if
    if (allocated(error)) then
@@ -208,7 +217,7 @@ subroutine newTBLiteWavefunction(env, mol, calc, chk)
    !> Molecular structure data
    type(TMolecule), intent(in) :: mol
    !> Instance of the new calculator
-   type(TTBLiteCalculator), intent(in) :: calc
+   type(TTBLiteCalculator), intent(inout) :: calc
    !> Wavefunction data
    type(TRestart), intent(inout) :: chk
 
@@ -229,6 +238,23 @@ subroutine newTBLiteWavefunction(env, mol, calc, chk)
          call sad_guess(struc, calc%tblite, wfn)
       case("eeq")
          call eeq_guess(struc, calc%tblite, wfn)
+      case("ceh")
+         block 
+            use tblite_context, only : context_type, context_terminal
+            use tblite_context_terminal, only : escape
+            use tblite_ceh_singlepoint, only : ceh_guess
+            use tblite_lapack_solver, only : lapack_solver 
+            use tblite_lapack_solver, only : lapack_algorithm
+            type(context_type) :: ctx
+            
+            ctx%solver = lapack_solver(lapack_algorithm%gvd)
+            ctx%terminal = context_terminal(calc%color)
+
+            write (env%unit, '(1x,a)') escape(ctx%terminal%cyan) // "Calculation of CEH charges" // &
+               & escape(ctx%terminal%reset)
+            
+            call ceh_guess(ctx, calc%tblite, struc, error, wfn, calc%accuracy, 1)
+         end block
       end select
    end associate
    if (allocated(error)) then
@@ -362,6 +388,51 @@ subroutine get_spin_constants(wll, mol, bas)
 end subroutine get_spin_constants
 #endif
 
+!> get CEH charges via tblite
+subroutine get_ceh(env,mol,tblite)
+
+   use xtb_propertyoutput, only : print_charges
+
+   !> computational environment
+   type(TEnvironment), intent(inout) :: env
+
+   !> molecular structure data
+   type(TMolecule), intent(in) :: mol
+
+   !> tblite input
+   type(TTBLiteInput), intent(in) :: tblite
+
+   !> initialize temporary calculator for CEH
+   type(TTBLiteCalculator) :: calc_ceh
+            
+   !> initialize temporary input for CEH
+   type(TTBLiteInput) :: tblite_ceh
+   
+   !> initialize temporary wfn for CEH
+   type(TRestart) :: chk_ceh
+
+   integer :: ich
+
+   tblite_ceh = tblite        ! copy the tblite input
+   tblite_ceh%method = "ceh" 
+#if WITH_TBLITE
+   write(env%unit, '(a)') repeat('-', 36)
+   
+   call newTBLiteCalculator(env, mol, calc_ceh, tblite_ceh)
+   call newTBLiteWavefunction(env, mol, calc_ceh, chk_ceh)
+
+   ! create ceh.charges file !
+   call open_file(ich, 'ceh.charges', 'w') 
+   call print_charges(ich, mol%n, chk_ceh%tblite%qat(:,1))
+   call close_file(ich)
+
+   write(env%unit, '(1x, a, /, a, /)') "CEH charges written to ceh.charges", repeat('-', 36)
+#else
+    call feature_not_implemented(env)
+#endif
+
+end subroutine get_ceh
+
 #if ! WITH_TBLITE
 subroutine feature_not_implemented(env)
    !> Computational environment
diff --git a/src/xhelp.f90 b/src/xhelp.f90
index e752a57f1..65ea587e6 100644
--- a/src/xhelp.f90
+++ b/src/xhelp.f90
@@ -110,6 +110,9 @@ subroutine help(iunit)
    "-c, --chrg INT",&
    "    specify molecular charge as INT, overrides .CHRG file and xcontrol option",&
    "",&
+   "--ceh",&
+   "    calculate CEH (Charge-Extended Hückel model) charges and write them to ceh.charges file",&
+   "",&
    "-u, --uhf INT",&
    "    specify number of unpaired electrons as INT, overrides .UHF file and xcontrol option",&
    "",&