allow tblite parameter readin. parameter files in toml format located in assets/tblite/ (#303)

Signed-off-by: Philipp Pracht <pp555@cam.ac.uk>

Author: pprcht

assets/tblite/gfn1-si.toml

@@ -77,21 +77,22 @@ subroutine crest_singlepoint(env,tim)
 
   allocate (grad(3,mol%nat),source=0.0_wp)
   calc = env%calc
+  calc%calcs(:)%prstdout = .true.
 
 !>--- print some info about the calculation
   call calc%info(stdout)
 
 !>--- and then start it
   write (stdout,'(a)') repeat('-',80)
-  write (stdout,'(a)',advance='no') '> Performing singlepoint calculations ... '
+  write (stdout,'(a)',advance='yes') '> Performing singlepoint calculations ... '
   flush (stdout)
 !>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<!
 !>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<!
   call engrad(mol,calc,energy,grad,io)
 !>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<!
 !>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<!
   call tim%stop(14)
-  write (stdout,*) 'done.'
+  write (stdout,'(a)') '> done.'
   write (atmp,'(a)') '> Total wall time for calculations'
   call tim%write_timing(stdout,14,trim(atmp),.true.)
   write (stdout,'(a)') repeat('-',80)

assets/tblite/gfn1.toml

@@ -77,9 +77,14 @@ subroutine tblite_engrad(mol,calc,energy,grad,iostatus)
 !>--- tblite printout handling
     call api_handle_output(calc,'tblite.out',mol,pr)
     if (pr) then
-      !> tblite uses its context (ctx)( type, rather than calc%prch
+      !> tblite uses its context (ctx) type, rather than calc%prch
       calc%tblite%ctx%unit = calc%prch
       calc%tblite%ctx%verbosity = 1
+      if(calc%prstdout)then
+        !> special case, fwd to stdout (be carefule with this!)
+        calc%tblite%ctx%unit = stdout
+        calc%tblite%ctx%verbosity = 2
+      endif
     else
       calc%tblite%ctx%verbosity = 0
     end if
@@ -100,7 +105,7 @@ subroutine tblite_engrad(mol,calc,energy,grad,iostatus)
     call tblite_singlepoint(mol,calc%chrg,calc%uhf,calc%tblite, &
     &                       energy,grad,iostatus)
     if (iostatus /= 0) return
-    call api_print_e_grd(pr,calc%tblite%ctx%unit,mol,energy,grad)
+    call api_print_e_grd(pr,calc%prch,mol,energy,grad)
 
 !>--- postprocessing, getting other data
     !$omp critical

assets/tblite/gfn2.toml

@@ -149,6 +149,9 @@ subroutine tblite_init(calc,loadnew)
     loadnew = .false.
     if (.not.allocated(calc%tblite)) then
       allocate (calc%tblite)
+      if(allocated(calc%tbliteparam))then
+        calc%tblite%paramfile = calc%tbliteparam
+      endif
       loadnew = .true.
     end if
     if (calc%apiclean) loadnew = .true.

assets/tblite/ipea1.toml

@@ -76,6 +76,7 @@ module calc_type
     integer :: prch = stdout  !> printout channel
     logical :: pr = .false.   !> allow the calculation to produce printout? Results in a lot I/O
     logical :: prappend = .false. !> append printout
+    logical :: prstdout = .false. !> special case, fwd some printout to stdout
     integer :: refine_lvl = 0 !> to allow defining different refinement levels
 
     integer :: chrg = 0          !> molecular charge
@@ -139,6 +140,7 @@ module calc_type
 
 !>--- tblite data
     type(tblite_data),allocatable :: tblite
+    character(len=:),allocatable :: tbliteparam
 
 !>--- GFN0-xTB data
     type(gfn0_data),allocatable          :: g0calc
@@ -685,6 +687,8 @@ subroutine calculation_settings_shortflag(self)
         self%shortflag = 'CEH'
       case (xtblvl%eeq)
         self%shortflag = 'EEQ(D4)'
+      case (xtblvl%param)
+        self%shortflag = 'parameter file: '//trim(self%tbliteparam)
       end select
     case( jobtype%gfn0 )
       self%shortflag =  'GFN0-xTB' 

src/algos/singlepoint.f90

@@ -32,6 +32,8 @@ module tblite_api
   use tblite_wavefunction_type,only:wavefunction_type,new_wavefunction
   use tblite_wavefunction,only:sad_guess,eeq_guess
   use tblite_xtb,xtb_calculator => xtb_calculator
+  use tblite_xtb_calculator, only: new_xtb_calculator
+  use tblite_param, only : param_record
   use tblite_results,only:tblite_resultstype => results_type
   use tblite_wavefunction_mulliken,only:get_molecular_dipole_moment
   use tblite_ceh_singlepoint,only:ceh_guess
@@ -65,6 +67,7 @@ module tblite_api
   type :: tblite_data
     integer  :: lvl = 0
     real(wp) :: accuracy = 1.0_wp
+    character(len=:),allocatable :: paramfile
     type(wavefunction_type)     :: wfn
     type(xtb_calculator)        :: calc
     type(tblite_ctx)            :: ctx
@@ -81,17 +84,13 @@ module tblite_api
     !> the guesses can be used for charges, but NOT for e+grd!
     integer :: eeq = 4
     integer :: ceh = 5
+    integer :: param = 6
   end type enum_tblite_method
   type(enum_tblite_method),parameter,public :: xtblvl = enum_tblite_method()
 
   !> Conversion factor from Kelvin to Hartree
   real(wp),parameter :: ktoau = 3.166808578545117e-06_wp
 
-  integer :: verbosity = 0
-  !> IMPORTANT: tblite is not entirely thread-safe
-  !> if verbosity is >0. We'll have to turn it off.
-  !> At least for statically compiled binaries
-
   public :: wavefunction_type,xtb_calculator
   public :: tblite_ctx,tblite_resultstype
   public :: tblite_setup,tblite_singlepoint,tblite_addsettings
@@ -121,10 +120,12 @@ subroutine tblite_setup(mol,chrg,uhf,lvl,etemp,tblite)
 #ifdef WITH_TBLITE
     type(structure_type) :: mctcmol
     type(error_type),allocatable :: error
+    type(param_record) :: param
 
     real(wp) :: etemp_au,energy
     real(wp),allocatable :: grad(:,:)
     logical :: pr
+    integer :: io
 
     pr = (tblite%ctx%verbosity > 0)
 
@@ -149,6 +150,19 @@ subroutine tblite_setup(mol,chrg,uhf,lvl,etemp,tblite)
     case (xtblvl%eeq)
       if (pr) call tblite%ctx%message("tblite> setting up D4 EEQ charges calculation")
       call new_ceh_calculator(tblite%calc,mctcmol) !> doesn't matter but needs initialization
+    case (xtblvl%param)
+      if (pr) call tblite%ctx%message("tblite> setting up xtb calculator from parameter file")
+      if(allocated(tblite%paramfile))then
+        call tblite_read_param_record(tblite%paramfile,param,io)
+        call new_xtb_calculator(tblite%calc, mctcmol, param, error)
+        if(allocated(error))then
+          write(stdout,*) 'Could not read tblite parameter file '//tblite%paramfile
+          error stop
+        endif
+      else
+        if (pr) call tblite%ctx%message("tblite> parameter file does not exist, defaulting to GFN2-xTB")
+        call new_gfn2_calculator(tblite%calc,mctcmol)
+      endif
     case default
       call tblite%ctx%message("Error: Unknown method in tblite!")
       error stop
@@ -202,7 +216,7 @@ subroutine tblite_add_solv(mol,chrg,uhf,tblite,smodel,solvent)
     pr = (tblite%ctx%verbosity > 0)
 
 !>--- some tblite calculators have nothing to do with implicit solvation
-    if (tblite%lvl > 3) then
+    if (tblite%lvl > 3 .and. tblite%lvl.ne.xtblvl%param) then
       if (pr) call tblite%ctx%message("tblite> skipping implicit solvation setup for this potential")
       return
     end if
@@ -282,11 +296,17 @@ subroutine tblite_singlepoint(mol,chrg,uhf,tblite,energy,gradient,iostatus)
     type(error_type),allocatable :: error
     real(wp) :: sigma(3,3)
     logical :: pr
+    integer :: verbosity
 
     iostatus = 0
     energy = 0.0_wp
     gradient(:,:) = 0.0_wp
     pr = (tblite%ctx%verbosity > 0)
+    if(tblite%ctx%verbosity>1)then
+      verbosity = tblite%ctx%verbosity
+    else
+      verbosity = 0
+    endif
 
 !>--- make an mctcmol object from mol
     call tblite_mol2mol(mol,chrg,uhf,mctcmol)
@@ -454,6 +474,42 @@ subroutine tblite_getdipole(mol,chrg,uhf,tblite,dipole)
 #endif
   end subroutine tblite_getdipole
 
+!========================================================================================!
+
+#ifdef WITH_TBLITE
+subroutine tblite_read_param_record(paramfile,param,io)
+   use tomlf
+   implicit none
+   character(len=*),intent(in) :: paramfile
+   type(param_record),intent(out) :: param
+   integer,intent(out) :: io
+   type(error_type),allocatable :: error
+   type(toml_table),allocatable :: table
+   type(toml_error),allocatable :: terror
+   type(toml_context) :: tcontext
+   logical,parameter :: color = .true.
+
+   io=1
+
+   call toml_load(table,paramfile,error=terror,context=tcontext, &
+   & config=toml_parser_config(color=color))
+   if(allocated(terror))then
+     io=1
+     return
+   endif
+
+   call param%load_from_toml(table,error)
+
+   if(allocated(error))then
+     io=1
+   else
+     io=0
+   endif
+   if(allocated(table))deallocate(table)
+
+end subroutine tblite_read_param_record
+#endif
+
 !========================================================================================!
 !========================================================================================!
 end module tblite_api

src/calculator/api_engrad.f90

@@ -343,6 +343,9 @@ subroutine parse_setting_c(env,job,key,val)
       case default
         job%tblitelvl = xtblvl%unknown
       end select
+    case('tblite_param')
+      job%tbliteparam = val
+      job%tblitelvl = xtblvl%param
 
     case ('orca_cmd')
       job%id = jobtype%orca