Skip to content

Commit

Permalink
Full ALPB/GBSA with empirical parameters (tblite#142)
Browse files Browse the repository at this point in the history 
Signed-off-by: haneug <[email protected]>
Signed-off-by: Philipp Pracht <[email protected]>
Co-authored-by: Philipp Pracht <[email protected]>
Co-authored-by: thfroitzheim <[email protected]>
Co-authored-by: Pit Steinbach <[email protected]>
Co-authored-by: Pit Steinbach <[email protected]>
Co-authored-by: Sebastian Ehlert <[email protected]>
  • Loading branch information
@haneug @pprcht
@thfroitzheim @pitsteinbach @awvwgk
6 people authored Dec 11, 2024
1 parent b12d609 commit 9ca469a
Show file tree
Hide file tree
Showing 240 changed files with 15,049 additions and 2,666 deletions.
138 changes: 120 additions & 18 deletions app/cli.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ module tblite_cli
use tblite_features, only : get_tblite_feature
use tblite_lapack_solver, only : lapack_algorithm
use tblite_solvation, only : solvation_input, cpcm_input, alpb_input, &
& solvent_data, get_solvent_data
& cds_input, shift_input, solvent_data, get_solvent_data, solution_state, born_kernel
use tblite_version, only : get_tblite_version
implicit none
private
Expand Down Expand Up @@ -248,8 +248,10 @@ subroutine get_run_arguments(config, list, start, error)
integer :: iarg, narg
logical :: getopts
character(len=:), allocatable :: arg
logical :: alpb
type(solvent_data) :: solvent
logical :: solvent_not_found, parametrized_solvation
logical, allocatable :: alpb
integer, allocatable :: kernel, sol_state
type(solvent_data), allocatable :: solvent

iarg = start
getopts = .true.
Expand Down Expand Up @@ -335,46 +337,107 @@ subroutine get_run_arguments(config, list, start, error)
iarg = iarg + 1
call list%get(iarg, config%guess)

case("--solv-state")
iarg = iarg + 1
call list%get(iarg, arg)
if (.not.allocated(arg)) then
call fatal_error(error, "Missing solution state")
exit
end if

select case(arg)
case("gsolv")
sol_state = solution_state%gsolv
case("bar1mol")
sol_state = solution_state%bar1mol
case("reference")
sol_state = solution_state%reference
case default
call fatal_error(error, "Unknown solution state '"//arg//"'")
exit
end select

case("--born-kernel")
iarg = iarg + 1
call list%get(iarg, arg)
if (.not.allocated(arg)) then
call fatal_error(error, "Missing Born kernel")
exit
end if

select case(arg)
case("p16")
kernel = born_kernel%p16
case("still")
kernel = born_kernel%still
case default
call fatal_error(error, "Unknown Born kernel '"//arg//"'")
exit
end select

case("--cpcm")
if (allocated(config%solvation)) then
call fatal_error(error, "Cannot use CPCM if ALPB/GBSA is enabled")
if (allocated(solvent)) then
call fatal_error(error, "Cannot use multiple solvation models")
exit
end if
parametrized_solvation = .false.

! Check for solvent information
iarg = iarg + 1
call list%get(iarg, arg)
if (.not.allocated(arg)) then
call fatal_error(error, "Missing argument for CPCM")
exit
end if

solvent_not_found = .false.
allocate(solvent)
solvent = get_solvent_data(arg)
if (solvent%eps <= 0.0_wp) then
solvent_not_found = .true.
call get_argument_as_real(arg, solvent%eps, error)
end if
if (allocated(error)) exit
allocate(config%solvation)
config%solvation%cpcm = cpcm_input(solvent%eps)

case("--alpb", "--gbsa")
if (allocated(config%solvation)) then
call fatal_error(error, "Cannot use ALPB/GBSA if CPCM is enabled")
case("--gb", "--gbe")
if (allocated(solvent)) then
call fatal_error(error, "Cannot use multiple solvation models")
exit
end if
alpb = arg == "--alpb"
end if
parametrized_solvation = .false.
alpb = arg == "--gbe"

! Check for solvent information
iarg = iarg + 1
call list%get(iarg, arg)
if (.not.allocated(arg)) then
call fatal_error(error, "Missing argument for ALPB/GBSA")
exit
end if

allocate(solvent)
solvent = get_solvent_data(arg)
if (solvent%eps <= 0.0_wp) then
call get_argument_as_real(arg, solvent%eps, error)
end if
if (allocated(error)) exit
allocate(config%solvation)
config%solvation%alpb = alpb_input(solvent%eps, alpb=alpb)

case("--alpb", "--gbsa")
if (allocated(solvent)) then
call fatal_error(error, "Cannot use multiple solvation models")
exit
end if
parametrized_solvation = .true.
alpb = arg == "--alpb"

! Check for solvent information
iarg = iarg + 1
call list%get(iarg, arg)
if (.not.allocated(arg)) then
call fatal_error(error, "Missing argument for ALPB/GBSA")
exit
end if
allocate(solvent)
solvent = get_solvent_data(arg)
if (allocated(error)) exit

case("--param")
if (allocated(config%param)) then
Expand Down Expand Up @@ -448,10 +511,12 @@ subroutine get_run_arguments(config, list, start, error)

case("--grad")
config%grad = .true.
config%grad_output = "tblite.txt"
iarg = iarg + 1
call list%get(iarg, arg)
if (allocated(arg)) then
if (arg(1:1) == "-") then
if (arg(1:1) == "-" .or. &
& iarg == narg .and. .not.allocated(config%input)) then
iarg = iarg - 1
cycle
end if
Expand All @@ -464,7 +529,8 @@ subroutine get_run_arguments(config, list, start, error)
iarg = iarg + 1
call list%get(iarg, arg)
if (allocated(arg)) then
if (arg(1:1) == "-") then
if (arg(1:1) == "-" .or. &
& iarg == narg .and. .not.allocated(config%input)) then
iarg = iarg - 1
cycle
end if
Expand All @@ -480,6 +546,42 @@ subroutine get_run_arguments(config, list, start, error)
end if
end if

if (allocated(solvent)) then
if (.not.allocated(sol_state)) then
sol_state = solution_state%gsolv
end if
if (allocated(alpb)) then
! ALPB/GBSA solvation model
if (.not.allocated(kernel)) then
kernel = merge(born_kernel%still, born_kernel%p16, alpb)
end if

if (.not.parametrized_solvation .and. sol_state /= solution_state%gsolv) then
call fatal_error(error, "Solution state shift is only supported for named solvents")
return
end if

allocate(config%solvation)
if (parametrized_solvation) then
config%solvation%alpb = alpb_input(solvent%eps, solvent=solvent%solvent, &
& kernel=kernel, alpb=alpb)
config%solvation%cds = cds_input(alpb=alpb, solvent=solvent%solvent)
config%solvation%shift = shift_input(alpb=alpb, solvent=solvent%solvent, &
& state=sol_state)
else
config%solvation%alpb = alpb_input(solvent%eps, kernel=kernel, alpb=alpb)
end if
else
! CPCM solvation model
if (sol_state /= solution_state%gsolv) then
call fatal_error(error, "Solution state shift not supported for CPCM")
return
end if
allocate(config%solvation)
config%solvation%cpcm = cpcm_input(solvent%eps)
end if
end if

if (.not.allocated(config%guess)) config%guess = "sad"

if (config%grad .and. .not.allocated(config%grad_output) .and. .not.config%json) then
Expand Down
34 changes: 23 additions & 11 deletions app/cli_help.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,7 @@ module tblite_cli_help
" --method <name> Parametrization of the xTB Hamiltonian to use"//nl//&
" Available methods: gfn1, gfn2, ipea1 (Default: gfn2)"//nl//&
" --param <file> Parametrization file to use for calculation"//nl//&
" --acc <real> Convergence criterion for SCF (Default: 1.0)"//nl//&
" --etemp <real> Electronic temperature for calculation (Default: 300K)"//nl//&
" --guess <name> Guess for the initial populations, possible options:"//nl//&
" sad (default), eeq, ceh (Charge-Extended Hückel method)"//nl//&
Expand All @@ -121,16 +122,27 @@ module tblite_cli_help
" gvd (default), and gvr"//nl//&
" --efield <real>,<real>,<real>"//nl//&
" Homogeneous electric field in V/Å."//nl//&
"--alpb <real> Use analytical linearized Poisson-Boltzmann solvation."//nl//&
" Solvent is specified by dielectric constant."//nl//&
"--gbsa <real> Use generalized Born solvation model."//nl//&
" Solvent is specified by dielectric constant."//nl//&
"--cpcm <real> Use polarizable continuum solvation model."//nl//&
" Solvent is specified by dielectric constant."//nl//&
" --alpb <name> Use analytical linearized Poisson-Boltzmann (ALPB) solvation model."//nl//&
" Solvent is specified by the solvent name."//nl//&
" Uses parametrized ALPB with CDS and empirical shift."//nl//&
" --gbsa <name> Use generalized Born solvation model (GBSA)."//nl//&
" Solvent is specified by the solvent name."//nl//&
" Uses parametrized GB with CDS and empirical shift."//nl//&
" --gbe <real>/<name> Use generalized Born for finite epsilion (GBE) solvation model."//nl//&
" Solvent is specified by dielectric constant or the solvent name."//nl//&
" --gb <real>/<name> Use generalized Born solvation model (GB)."//nl//&
" Solvent is specified by dielectric constant or solvent name."//nl//&
" --cpcm Use polarizable continuum solvation model (CPCM)."//nl//&
" Solvent is specified by dielectric constant or solvent name."//nl//&
" --born-kernel <name> Specify Born kernel to use with ALPB, GBSA or GB solvation model."//nl//&
" Possible options are p16 (default for ALPB) and still (default for GB/GBSA)."//nl//&
" --solv-state <name> Solution state correction: gsolv (default), bar1mol, reference."//nl//&
" --spin-polarized Use spin-polarized xTB Hamiltonian"//nl//&
"--post-processing <file> Add post processing methods to the calculation"//nl//&
" --post-processing <file>"//nl//&
" Add post processing methods to the calculation"//nl//&
" by using a toml file as input."//nl//&
"--post-processing <name> Add post processing methods to the calculation,"//nl//&
" --post-processing <name>"//nl//&
" Add post processing methods to the calculation,"//nl//&
" Mayer-Wiberg bond orders are computed by default."//nl//&
" Options: molmom (molecular multipole moments)"//nl//&
" --grad [file] Evaluate molecular gradient and virial"//nl//&
Expand Down Expand Up @@ -158,14 +170,14 @@ module tblite_cli_help
" -c, --charge <real> Set charge to molecule, overwrites .CHRG file"//nl//&
" --spin <int> Set number of unpaired electrons, overwrites .UHF file"//nl//&
" --method <name> Guess for the initial populations, possible options:"//nl//&
" sad (default), eeq, ceh (Charge-Extended Hückel method)"//nl//&
" --etemp-guess <real> Electronic temperature for ceh-guess (Default: 5000K)"//nl//&
" sad, eeq, ceh (Charge-Extended Hückel method, default)"//nl//&
" --etemp-guess <real> Electronic temperature for ceh-guess (Default: 4000K)"//nl//&
" --solver <name> Electronic solver for SCF, possible options:"//nl//&
" gvd (default), and gvr"//nl//&
" --efield <real>,<real>,<real>"//nl//&
" Homogeneous electric field in V/Å."//nl//&
" --grad Evaluate analytic gradient of charges."//nl//&
! " --json [file] Dump results as JSON output (default: tblite.json)"//nl//&
" --json [file] Dump results as JSON output (default: tblite.json)"//nl//&
" -i, --input <format> Hint for the format of the input file"//nl//&
" -v, --verbose Increase verbosity of printout"//nl//&
" -s, --silent Reduce verbosity of printout"//nl//&
Expand Down
44 changes: 31 additions & 13 deletions app/driver_run.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ module tblite_driver_run
use tblite_param, only : param_record
use tblite_results, only : results_type
use tblite_spin, only : spin_polarization, new_spin_polarization
use tblite_solvation, only : new_solvation, solvation_type
use tblite_solvation, only : new_solvation, new_solvation_cds, new_solvation_shift, solvation_type
use tblite_wavefunction, only : wavefunction_type, new_wavefunction, &
& sad_guess, eeq_guess, shell_partition
use tblite_xtb_calculator, only : xtb_calculator, new_xtb_calculator
Expand Down Expand Up @@ -159,10 +159,6 @@ subroutine run_main(config, error)
call new_ceh_calculator(calc_ceh, mol, error)
if (allocated(error)) return
call new_wavefunction(wfn_ceh, mol%nat, calc_ceh%bas%nsh, calc_ceh%bas%nao, 1, config%etemp_guess * kt)
if (config%grad) then
call ctx%message("WARNING: CEH gradient not yet implemented. Stopping.")
return
end if
end if

if (allocated(config%efield)) then
Expand All @@ -171,13 +167,13 @@ subroutine run_main(config, error)
cont = electric_field(config%efield*vatoau)
call calc%push_back(cont)
end block
if (config%guess == "ceh") then
block
class(container_type), allocatable :: cont
cont = electric_field(config%efield*vatoau)
call calc_ceh%push_back(cont)
end block
end if
if (config%guess == "ceh") then
block
class(container_type), allocatable :: cont
cont = electric_field(config%efield*vatoau)
call calc_ceh%push_back(cont)
end block
end if
end if

select case(config%guess)
Expand Down Expand Up @@ -224,14 +220,36 @@ subroutine run_main(config, error)
end if

if (allocated(config%solvation)) then
method = "gfn2"
if (allocated(config%method)) method = config%method
block
class(container_type), allocatable :: cont
class(solvation_type), allocatable :: solv
call new_solvation(solv, mol, config%solvation, error)
call new_solvation(solv, mol, config%solvation, error, method)
if (allocated(error)) return
call move_alloc(solv, cont)
call calc%push_back(cont)
end block
if (allocated(config%solvation%cds)) then
block
class(container_type), allocatable :: cont
class(solvation_type), allocatable :: cds
call new_solvation_cds(cds, mol, config%solvation, error, method)
if (allocated(error)) return
call move_alloc(cds, cont)
call calc%push_back(cont)
end block
end if
if (allocated(config%solvation%shift)) then
block
class(container_type), allocatable :: cont
class(solvation_type), allocatable :: shift
call new_solvation_shift(shift, config%solvation, error, method)
if (allocated(error)) return
call move_alloc(shift, cont)
call calc%push_back(cont)
end block
end if
end if

wbo_label = "bond-orders"
Expand Down
Loading

0 comments on commit 9ca469a

Please sign in to comment.