diff --git a/man/xtb.1.adoc b/man/xtb.1.adoc
index 4f417b48a..0f0af18db 100644
--- a/man/xtb.1.adoc
+++ b/man/xtb.1.adoc
@@ -81,8 +81,10 @@ 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
+*--ceh <grad> [REAL]* ::
+      calculate CEH (Charge-Extended Hückel model) charges and write them to the ceh.charges file,
+      optionally, calculate numerical gradients and write them to the ceh.charges.numgrad file
+      with an adjustable step size for the numerical gradients.
 
 *--ptb* ::
     performs single-point calculation with the density tight-binding method PTB.
diff --git a/src/prog/main.F90 b/src/prog/main.F90
index 746278bea..d043b9cae 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -94,7 +94,7 @@ module xtb_prog_main
    use xtb_oniom, only: oniom_input, TOniomCalculator, calculateCharge
    use xtb_vertical, only: vfukui
    use xtb_tblite_calculator, only: TTBLiteCalculator, TTBLiteInput, &
-         & newTBLiteWavefunction, get_ceh
+         & newTBLiteWavefunction, get_ceh, num_grad_chrg
    use xtb_ptb_calculator, only: TPTBCalculator
    use xtb_solv_cpx, only: TCpcmx
    use xtb_dipro, only: get_jab, jab_input
@@ -544,7 +544,7 @@ subroutine xtbMain(env, argParser)
          select case (set%runtyp)
          case default
             call env%terminate('This is an internal error, please define your runtypes!')
-         case (p_run_scc, p_run_grad, p_run_opt, p_run_hess, p_run_ohess, p_run_bhess, &
+         case (p_run_prescc,p_run_scc, p_run_grad, p_run_opt, p_run_hess, p_run_ohess, p_run_bhess, &
                p_run_md, p_run_omd, p_run_path, p_run_screen, &
                p_run_modef, p_run_mdopt, p_run_metaopt)
             if (set%mode_extrun == p_ext_gfnff) then
@@ -629,8 +629,18 @@ subroutine xtbMain(env, argParser)
       end select
 
       ! get CEH charges !
-      if (tblite%ceh) &
-         call get_ceh(env,mol,tblite)
+      if (allocated(tblite%ceh)) then
+         ! if numercal charges are requested !
+         if (tblite%ceh%grad) then
+            call num_grad_chrg(env,mol,tblite)
+         ! only charges !
+         else
+            call get_ceh(env,mol,tblite)
+         endif 
+
+         ! stop calculation !
+         call finalize_xtb(env)
+      end if
 
       ! ------------------------------------------------------------------------
       !> printout a header for the exttyp
@@ -1233,44 +1243,7 @@ subroutine xtbMain(env, argParser)
 
       ! ------------------------------------------------------------------------
       !  make some post processing afterward, show some timings and stuff
-      write (env%unit, '(a)')
-      write (env%unit, '(72("-"))')
-      call stop_timing_run
-      call stop_timing(1)
-      call prdate('E')
-      write (env%unit, '(72("-"))')
-      call prtiming(1, 'total')
-      call prtiming(2, 'SCF')
-      if ((set%runtyp == p_run_opt) .or. (set%runtyp == p_run_ohess) .or. &
-         &   (set%runtyp == p_run_omd) .or. (set%runtyp == p_run_metaopt)) then
-         call prtiming(3, 'ANC optimizer')
-      end if
-      if (set%runtyp == p_run_path) then
-         call prtiming(4, 'path finder')
-      end if
-      if (((set%runtyp == p_run_hess) .or. (set%runtyp == p_run_ohess) .or. (set%runtyp == p_run_bhess))) then
-         if (set%mode_extrun /= p_ext_turbomole) then
-            call prtiming(5, 'analytical hessian')
-         else
-            call prtiming(5, 'numerical hessian')
-         end if
-      end if
-      if ((set%runtyp == p_run_md) .or. (set%runtyp == p_run_omd) .or. &
-          (set%runtyp == p_run_metaopt)) then
-         call prtiming(6, 'MD')
-      end if
-      if (set%runtyp == p_run_screen) then
-         call prtiming(8, 'screen')
-      end if
-      if (set%runtyp == p_run_modef) then
-         call prtiming(9, 'mode following')
-      end if
-      if (set%runtyp == p_run_mdopt) then
-         call prtiming(10, 'MD opt.')
-      end if
-
-      write (env%unit, '(a)')
-      call terminate(0)
+      call finalize_xtb(env) 
 
    end subroutine xtbMain
 
@@ -1471,7 +1444,25 @@ subroutine parseArguments(env, args, inputFile, paramFile, lgrad, &
          
          case('--ceh')
             if (get_xtb_feature('tblite')) then
-               tblite%ceh = .true.
+               allocate(tblite%ceh)
+               call set_runtyp('prescc')
+               ! check if numerical gradients should be calculated !
+               call args%nextArg(sec)
+               if (allocated(sec)) then
+                  if (sec == 'grad') then
+                     tblite%ceh%grad = .true.
+                     ! check if step size is provided !
+                     call args%nextArg(sec)
+                     if (allocated(sec)) then
+                        if (getValue(env, sec, ddum)) then
+                           tblite%ceh%step = ddum
+                        end if
+                     end if
+
+                  else
+                     call env%warning("Unknown CEH option '"//sec//"' provided", source)
+                  end if
+               end if 
             else
                call env%error("CEH charges are only available through tblite library", source)
                return
@@ -2038,4 +2029,52 @@ subroutine ptb_feature_not_implemented(env)
       call env%error("Please recompile without '-Dtblite=disabled' option or change meson setup.")
    end subroutine ptb_feature_not_implemented
 
+   !>  make some post processing afterward, show some timings and stuff
+   subroutine finalize_xtb(env)
+    
+      !> Calculation environment
+      type(TEnvironment), intent(in) :: env
+
+      write (env%unit, '(a)')
+      write (env%unit, '(72("-"))')
+      call stop_timing_run
+      call stop_timing(1)
+      call prdate('E')
+      write (env%unit, '(72("-"))')
+      call prtiming(1, 'total')
+      if (.not. set%runtyp == p_run_prescc) &
+         call prtiming(2, 'SCF')
+      if ((set%runtyp == p_run_opt) .or. (set%runtyp == p_run_ohess) .or. &
+         &   (set%runtyp == p_run_omd) .or. (set%runtyp == p_run_metaopt)) then
+         call prtiming(3, 'ANC optimizer')
+      end if
+      if (set%runtyp == p_run_path) then
+         call prtiming(4, 'path finder')
+      end if
+      if (((set%runtyp == p_run_hess) .or. (set%runtyp == p_run_ohess) .or. (set%runtyp == p_run_bhess))) then
+         if (set%mode_extrun /= p_ext_turbomole) then
+            call prtiming(5, 'analytical hessian')
+         else
+            call prtiming(5, 'numerical hessian')
+         end if
+      end if
+      if ((set%runtyp == p_run_md) .or. (set%runtyp == p_run_omd) .or. &
+          (set%runtyp == p_run_metaopt)) then
+         call prtiming(6, 'MD')
+      end if
+      if (set%runtyp == p_run_screen) then
+         call prtiming(8, 'screen')
+      end if
+      if (set%runtyp == p_run_modef) then
+         call prtiming(9, 'mode following')
+      end if
+      if (set%runtyp == p_run_mdopt) then
+         call prtiming(10, 'MD opt.')
+      end if
+
+      write (env%unit, '(a)')
+      call terminate(0)
+
+   end subroutine finalize_xtb
+
 end module xtb_prog_main
diff --git a/src/set_module.f90 b/src/set_module.f90
index dd5f93ca2..ad5f81ade 100644
--- a/src/set_module.f90
+++ b/src/set_module.f90
@@ -1031,7 +1031,8 @@ subroutine set_runtyp(typ)
    select case(typ)
    case default ! do nothing !
       call raise('E',typ//' is no valid runtyp (internal error)')
-
+   case ('prescc')
+      set%runtyp = p_run_prescc
    case('scc')
       set%runtyp = p_run_scc
 
@@ -1198,6 +1199,7 @@ subroutine set_chrg(env,val)
 
 end subroutine set_chrg
 
+
 subroutine set_spin(env,val)
    implicit none
    character(len=*), parameter :: source = 'set_spin'
diff --git a/src/setparam.f90 b/src/setparam.f90
index 20f4e0baf..811baf0b4 100644
--- a/src/setparam.f90
+++ b/src/setparam.f90
@@ -188,6 +188,7 @@ module xtb_setparam
    integer, parameter :: p_ext_ptb       = 17
    integer, parameter :: p_ext_mcgfnff   = 18
 
+   integer, parameter :: p_run_prescc  =   1
    integer, parameter :: p_run_scc    =   2
    integer, parameter :: p_run_grad   =   3
    integer, parameter :: p_run_opt    =   4
diff --git a/src/tblite/calculator.F90 b/src/tblite/calculator.F90
index f05087414..270e2ed8f 100644
--- a/src/tblite/calculator.F90
+++ b/src/tblite/calculator.F90
@@ -60,7 +60,14 @@ module xtb_tblite_calculator
    private
 
    public :: TTBLiteCalculator, TTBLiteInput, newTBLiteCalculator, newTBLiteWavefunction
-   public :: get_ceh
+   public :: get_ceh, num_grad_chrg
+
+   type, private :: ceh
+      !> numerical gradients
+      logical :: grad
+      !> step size for numerical gradients
+      real(wp) :: step = 0.00001_wp
+   endtype
 
    !> Input for tblite library
    type :: TTBLiteInput
@@ -77,7 +84,7 @@ module xtb_tblite_calculator
       !> Colorful output
       logical :: color = .false.
       !> CEH charges
-      logical :: ceh = .false.
+      type(ceh), allocatable :: ceh
    end type TTBLiteInput
 
    !> Calculator interface for xTB based methods
@@ -250,10 +257,11 @@ subroutine newTBLiteWavefunction(env, mol, calc, chk)
             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)
+            ! temporary turn off colorful output
+            ! ctx%terminal = context_terminal(calc%color)
+            ! write (env%unit, '(0x,a)') escape(ctx%terminal%cyan) // "Calculation of CEH charges" // &
+            !    & escape(ctx%terminal%reset)
 
             call ceh_singlepoint(ctx, calc%tblite, struc, wfn, calc%accuracy, 1)
          end block
@@ -387,7 +395,7 @@ end subroutine get_spin_constants
 #endif
 
 !> get CEH charges via tblite
-subroutine get_ceh(env,mol,tblite)
+subroutine get_ceh(env,mol,tblite, ceh_chrg)
 
    use xtb_propertyoutput, only : print_charges
 
@@ -400,6 +408,9 @@ subroutine get_ceh(env,mol,tblite)
    !> tblite input
    type(TTBLiteInput), intent(in) :: tblite
 
+   !> CEH charges, if present assumed the numerical gradients are requested
+   real(wp), allocatable, optional, intent(out) :: ceh_chrg(:)
+
    !> initialize temporary calculator for CEH
    type(TTBLiteCalculator) :: calc_ceh
             
@@ -414,23 +425,87 @@ subroutine get_ceh(env,mol,tblite)
    tblite_ceh = tblite        ! copy the tblite input
    tblite_ceh%method = "ceh" 
 #if WITH_TBLITE
-   write(env%unit, '(a)') repeat('-', 36)
+
+   if (.not. present(ceh_chrg)) &
+      write(env%unit, '(1x, a, /, a)') "Calculation of CEH charges",repeat('-', 36)
    
    call newTBLiteCalculator(env, mol, calc_ceh, tblite_ceh)
    call newTBLiteWavefunction(env, mol, calc_ceh, chk_ceh)
+   
+   if (present(ceh_chrg)) then
+      allocate(ceh_chrg(mol%n))
+      ceh_chrg = chk_ceh%tblite%qat(:,1)  
+   else
+      ! 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)') "CEH charges written to ceh.charges"
+   endif
 
-   ! 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
 
+!> get numerical gradients for charges
+subroutine num_grad_chrg(env, mol, tblite)
+   !> Calculation environment
+   type(TEnvironment), intent(inout) :: env
+   
+   !> Molecular structure data
+   type(TMolecule), intent(inout) :: mol
+   
+   !> step size for numerical gradients
+   type(TTBLiteInput), intent(in) :: tblite
+   
+   !> numerical gradients
+   real(wp) :: numgrad(3, mol%n, mol%n)
+   
+   real(wp), allocatable, dimension(:) :: cehr, cehl
+   !
+   integer :: i,j,k, ich
+
+   real(wp) :: step, step2 ! for numerical gradient
+   
+   numgrad=0.0_wp
+   step = tblite%ceh%step
+   step2 = 0.5_wp / step
+   call get_ceh(env,mol,tblite)
+   
+   !$omp parallel do private(j,cehr,cehl) shared(env, numgrad, mol, tblite, step, step2) 
+   do i = 1, mol%n
+      do j = 1, 3
+         mol%xyz(j,i) = mol%xyz(j,i) + step
+         call get_ceh(env, mol, tblite, cehr)
+
+         mol%xyz(j,i) = mol%xyz(j,i) - 2*step
+         call get_ceh(env, mol, tblite, cehl)
+         
+         numgrad(j,i,:) = step2 * (cehr - cehl) ! numerical gradient
+         mol%xyz(j,i) = mol%xyz(j,i) + step ! reset the coordinates
+         
+      enddo
+   enddo
+   !$omp end parallel do
+
+   ! write the numerical gradient to the ceh.charges.numgrad file
+   call open_file(ich, 'ceh.charges.numgrad', 'w')
+   do i = 1, mol%n
+      do j = 1, mol%n
+         do k = 1, 3
+            write(ich, '(3f12.6)') numgrad(k,j,i)
+         enddo
+      enddo
+   enddo
+   call close_file(ich)
+   write(env%unit, '(1x, a)') "CEH gradients written to ceh.charges.numgrad"
+
+ end subroutine num_grad_chrg
+ 
+
+
 #if ! WITH_TBLITE
 subroutine feature_not_implemented(env)
    !> Computational environment
diff --git a/src/xhelp.f90 b/src/xhelp.f90
index 98fdd44d0..84f3db4dd 100644
--- a/src/xhelp.f90
+++ b/src/xhelp.f90
@@ -110,8 +110,10 @@ 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",&
+   "--ceh <grad> [REAL]",&
+   "    calculate CEH (Charge-Extended Hückel model) charges and write them to the ceh.charges file",&
+   "    optionally, calculate numerical gradients and write them to the ceh.charges.numgrad file",&
+   "    with an adjustable step size for the numerical gradients.",&
    "",&
    "-u, --uhf INT",&
    "    specify number of unpaired electrons as INT, overrides .UHF file and xcontrol option",&