Skip to content

Commit

Permalink
CEH gradients (#1153)
Browse files Browse the repository at this point in the history 
* add finalize_xtb for dynamic termination

Signed-off-by: albert <[email protected]>

* add num_grad_chrg for CEH gradients

Signed-off-by: albert <[email protected]>

* Test CI

* Renamed file ceh.grad to ceh.charges.numgrad

* add CLI description

Signed-off-by: albert <[email protected]>

* Modified order of ceh numgrad written to file

---------

Signed-off-by: albert <[email protected]>
Co-authored-by: lukaswittmann <[email protected]>
Co-authored-by: Lukas Wittmann <[email protected]>
  • Loading branch information
@Albkat @lukaswittmann
3 people authored Feb 5, 2025
1 parent 9435410 commit e8cac25
Show file tree
Hide file tree
Showing 6 changed files with 182 additions and 61 deletions.
6 changes: 4 additions & 2 deletions man/xtb.1.adoc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
125 changes: 82 additions & 43 deletions src/prog/main.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1230,44 +1240,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

Expand Down Expand Up @@ -1468,7 +1441,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
Expand Down Expand Up @@ -2035,4 +2026,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
4 changes: 3 additions & 1 deletion src/set_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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'
Expand Down
1 change: 1 addition & 0 deletions src/setparam.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
101 changes: 88 additions & 13 deletions src/tblite/calculator.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand All @@ -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
Expand Down
6 changes: 4 additions & 2 deletions src/xhelp.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",&
Expand Down

0 comments on commit e8cac25

Please sign in to comment.