diff --git a/.gitmodules b/.gitmodules
index 30a29d29..0e699866 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -19,5 +19,5 @@
 	url = https://github.com/fortran-lang/test-drive.git
 [submodule "subprojects/tblite"]
 	path = subprojects/tblite
-	url = https://github.com/tblite/tblite.git
-    branch = main
+	url = https://github.com/pprcht/tblite.git
+    branch = xtb_solvation
diff --git a/src/calculator/calculator.F90 b/src/calculator/calculator.F90
index 8f038d12..66157a14 100644
--- a/src/calculator/calculator.F90
+++ b/src/calculator/calculator.F90
@@ -141,10 +141,10 @@ subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
     iostatus = 0
     dum1 = 1.0_wp
     dum2 = 1.0_wp
-    if(n>0)then
-    calc%etmp = 0.0_wp
-    calc%grdtmp = 0.0_wp
-    endif
+    if (n > 0) then
+      calc%etmp = 0.0_wp
+      calc%grdtmp = 0.0_wp
+    end if
     !$omp end critical
 
 !**************************************
@@ -260,13 +260,13 @@ subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
 
       !!$omp parallel &
       !!$omp shared(calc%cons,mol,n,energy,gradient) &
-      !!$omp private(efix,i,calc%grdfix) 
+      !!$omp private(efix,i,calc%grdfix)
       do i = 1,calc%nconstraints
         efix = 0.0_wp
         calc%grdfix = 0.0_wp
-        if(calc%nfreeze > 0)then
-           call calc%cons(i)%addfreeze(calc%freezelist)
-        endif
+        if (calc%nfreeze > 0) then
+          call calc%cons(i)%addfreeze(calc%freezelist)
+        end if
         if (calc%cons(i)%type >= 0) then
           !>--- structural constraints
           call calc_constraint(mol%nat,mol%xyz,calc%cons(i),efix,calc%grdfix)
@@ -363,13 +363,13 @@ subroutine potential_core(molptr,calc,id,iostatus)
       !> ORCA-style SPs
       call orca_engrad(molptr,calc%calcs(id),calc%etmp(id),calc%grdtmp(:,1:pnat,id),iostatus)
 
-    case (jobtype%lj) 
+    case (jobtype%lj)
       !> Lennard-Jones potential calculation
       if (allocated(calc%calcs(id)%other)) then
         read (calc%calcs(id)%other,*) dum1,dum2
       else
         dum2 = 3.4_wp !> Argon-LJ distance (σ) in Angstroem
-        dum1 = 0.238129_wp/627.50947428_wp !> Argon-LJ energy (ε) in atomic units  
+        dum1 = 0.238129_wp/627.50947428_wp !> Argon-LJ energy (ε) in atomic units
       end if
       call lj_engrad(molptr%nat,molptr%xyz*autoaa,dum1,dum2, &
       &              calc%etmp(id),calc%grdtmp(:,1:pnat,id))
@@ -386,52 +386,62 @@ end subroutine potential_core
 !========================================================================================!
 !========================================================================================!
 
-  subroutine numgrad(mol,calc,angrad)
+  subroutine numgrad(mol,calc,angrad,numgradient)
 !*******************************************************
 !* subroutine numgrad
 !* routine to perform a numerical gradient calculation
 !*******************************************************
     implicit none
-
-    type(coord) :: mol
-    type(calcdata) :: calc
-    real(wp) :: angrad(3,mol%nat)
+    type(coord),intent(inout) :: mol
+    type(calcdata),intent(inout) :: calc
+    real(wp),intent(in),optional :: angrad(3,mol%nat)
+    real(wp),allocatable,intent(out),optional :: numgradient(:,:)
     integer :: i,j,k,l,ich,och,io
 
     real(wp) :: energy,el,er
     real(wp),allocatable :: grad(:,:)
-    real(wp),allocatable :: numgrd(:,:)
-    real(wp),parameter :: step = 0.00001_wp,step2 = 0.5_wp/step
+    real(wp),allocatable :: ngrd(:,:)
+    real(wp),parameter :: step = 0.0005_wp
+    real(wp),parameter :: step2 = 0.5_wp/step
 
     allocate (grad(3,mol%nat),source=0.0_wp)
-    allocate (numgrd(3,mol%nat),source=0.0_wp)
+    allocate (ngrd(3,mol%nat),source=0.0_wp)
 
     do i = 1,mol%nat
       do j = 1,3
+        call calc%dealloc_params()
         mol%xyz(j,i) = mol%xyz(j,i)+step
         call engrad(mol,calc,er,grad,io)
 
+        call calc%dealloc_params()
         mol%xyz(j,i) = mol%xyz(j,i)-2*step
         call engrad(mol,calc,el,grad,io)
 
         mol%xyz(j,i) = mol%xyz(j,i)+step
-        numgrd(j,i) = step2*(er-el)
+        ngrd(j,i) = step2*(er-el)
       end do
     end do
 
-    write (stdout,*)
-    write (stdout,*) 'Numerical Gradient:'
-    do i = 1,mol%nat
-      write (*,'(3f18.8)') numgrd(1:3,i)
-    end do
+    if (present(angrad)) then
+      write (stdout,*)
+      write (stdout,*) 'Numerical Gradient:'
+      do i = 1,mol%nat
+        write (*,'(3f18.8)') ngrd(1:3,i)
+      end do
 
-    write (stdout,*)
-    write (stdout,*) 'Gradient Difference:'
-    do i = 1,mol%nat
-      write (stdout,'(3f18.8)') numgrd(1:3,i)-angrad(1:3,i)
-    end do
+      write (stdout,*)
+      write (stdout,*) 'Gradient Difference:'
+      do i = 1,mol%nat
+        write (stdout,'(3f18.8)') ngrd(1:3,i)-angrad(1:3,i)
+      end do
+    end if
 
-    deallocate (numgrd,grad)
+    if (present(numgradient)) then
+      call move_alloc(ngrd,numgradient)
+    else
+      deallocate (ngrd)
+    end if
+    deallocate (grad)
 
     return
   end subroutine numgrad
@@ -537,7 +547,7 @@ subroutine numhess2(nat,at,xyz,calc,hess,io)
     mol%nat = nat
     mol%at = at
     mol%xyz = xyz
-    
+
     allocate (gradr(3,mol%nat),source=0.0_wp) !dummy
     allocate (gradl(3,mol%nat),source=0.0_wp) !dummy
 
@@ -546,15 +556,15 @@ subroutine numhess2(nat,at,xyz,calc,hess,io)
 
 !>--- prepare potential dipole gradient readout
     rddipoles = (any(calc%calcs(:)%rddip))
-    if(rddipoles)then
-      do m=1,calc%ncalculations
-        if(calc%calcs(m)%rddip)then
-          calc%calcs(m)%rddipgrad = .true. 
-          if(allocated(calc%calcs(m)%dipgrad)) deallocate(calc%calcs(m)%dipgrad)
-          allocate(calc%calcs(m)%dipgrad(3,nat*3))  
-        endif
-      enddo
-    endif
+    if (rddipoles) then
+      do m = 1,calc%ncalculations
+        if (calc%calcs(m)%rddip) then
+          calc%calcs(m)%rddipgrad = .true.
+          if (allocated(calc%calcs(m)%dipgrad)) deallocate (calc%calcs(m)%dipgrad)
+          allocate (calc%calcs(m)%dipgrad(3,nat*3))
+        end if
+      end do
+    end if
 
 !>--- actual numerical Hessian loop
     do i = 1,mol%nat
@@ -564,13 +574,13 @@ subroutine numhess2(nat,at,xyz,calc,hess,io)
         call engrad(mol,calc,er,gradr,io)
 
         gradr_tmp = calc%grdtmp
-        if(rddipoles) call get_dipoles(calc,dipr)
+        if (rddipoles) call get_dipoles(calc,dipr)
 
         mol%xyz(j,i) = mol%xyz(j,i)-2.0_wp*step
         call engrad(mol,calc,el,gradl,io)
 
         gradl_tmp = calc%grdtmp
-        if(rddipoles) call get_dipoles(calc,dipl)
+        if (rddipoles) call get_dipoles(calc,dipl)
 
         mol%xyz(j,i) = mol%xyz(j,i)+step
 
@@ -585,13 +595,13 @@ subroutine numhess2(nat,at,xyz,calc,hess,io)
         end do
 
 !>--- create dipole gradients
-        if(rddipoles)then
-          do m=1,calc%ncalculations
-           if(calc%calcs(m)%rddip)then
-             calc%calcs(m)%dipgrad(1:3, ii) = (dipr(1:3,m) - dipl(1:3,m)) * step2           
-           endif
-          enddo
-        endif
+        if (rddipoles) then
+          do m = 1,calc%ncalculations
+            if (calc%calcs(m)%rddip) then
+              calc%calcs(m)%dipgrad(1:3,ii) = (dipr(1:3,m)-dipl(1:3,m))*step2
+            end if
+          end do
+        end if
 
       end do
     end do
@@ -609,8 +619,8 @@ subroutine numhess2(nat,at,xyz,calc,hess,io)
 
     call engrad(mol,calc,el,gradl,io) !>- to get the gradient of the non-displaced s
 
-    if(allocated(dipl)) deallocate(dipl)
-    if(allocated(dipr)) deallocate(dipr)
+    if (allocated(dipl)) deallocate (dipl)
+    if (allocated(dipr)) deallocate (dipr)
     deallocate (gradl_tmp,gradr_tmp)
     deallocate (gradl,gradr)
     call mol%deallocate()
@@ -652,18 +662,18 @@ subroutine constrhess(nat,at,xyz,calc,phess)
     dummycalc%id = -1000
     dummycalc%ncalculations = 0
     dummycalc%pr_energies = .false.
-    
+
     !> copy the constraints
-    dummycalc%nfreeze = calc%nfreeze 
+    dummycalc%nfreeze = calc%nfreeze
     dummycalc%nconstraints = ncons
     !$omp critical
-    allocate(dummycalc%cons( ncons ))
-    do i=1,ncons
-      dummycalc%cons(i) = calc%cons(i)   
-    enddo
-    if(calc%nfreeze > 0)then
-     dummycalc%freezelist = calc%freezelist
-    endif
+    allocate (dummycalc%cons(ncons))
+    do i = 1,ncons
+      dummycalc%cons(i) = calc%cons(i)
+    end do
+    if (calc%nfreeze > 0) then
+      dummycalc%freezelist = calc%freezelist
+    end if
     n3 = nat*3
     allocate (hess(n3,n3),source=0.0_wp)
     !$omp end critical
@@ -677,7 +687,7 @@ subroutine constrhess(nat,at,xyz,calc,phess)
         phess(k) = phess(k)+0.5_wp*(hess(j,i)+hess(i,j))
       end do
     end do
-    
+
     deallocate (hess)
     return
   end subroutine constrhess
@@ -706,7 +716,7 @@ subroutine calc_print_energies(calc,energy,energies,gnorm,chnl)
       write (atmp,'(f20.12)') energies(i)
       btmp = btmp//atmp
     end do
-    write(atmp,'(e20.5)') gnorm
+    write (atmp,'(e20.5)') gnorm
     btmp = btmp//atmp
     if (present(chnl)) then
       write (chnl,'(a)') btmp
diff --git a/src/calculator/tblite_api.F90 b/src/calculator/tblite_api.F90
index 3e7a4edf..d16e5c56 100644
--- a/src/calculator/tblite_api.F90
+++ b/src/calculator/tblite_api.F90
@@ -32,8 +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_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_singlepoint
@@ -152,17 +152,17 @@ subroutine tblite_setup(mol,chrg,uhf,lvl,etemp,tblite)
       call new_ceh_calculator(tblite%calc,mctcmol,error) !> 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
+      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
+        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
+        end if
       else
         if (pr) call tblite%ctx%message("tblite> parameter file does not exist, defaulting to GFN2-xTB")
         call new_gfn2_calculator(tblite%calc,mctcmol,error)
-      endif
+      end if
     case default
       call tblite%ctx%message("Error: Unknown method in tblite!")
       error stop
@@ -189,8 +189,9 @@ subroutine tblite_add_solv(mol,chrg,uhf,tblite,smodel,solvent)
 #ifdef WITH_TBLITE
     use tblite_container,only:container_type
     use tblite_solvation,only:new_solvation,tblite_solvation_type => solvation_type, &
-    &                            solvent_data,get_solvent_data,solvation_input,  &
-    &                            cpcm_input,alpb_input,alpb_solvation
+    &                         solvent_data,get_solvent_data,solvation_input,  &
+    &                         cpcm_input,alpb_input,alpb_solvation, &
+    &                         cds_input,new_solvation_cds,shift_input,new_solvation_shift
 #endif
     implicit none
     type(coord),intent(in)  :: mol
@@ -207,7 +208,10 @@ subroutine tblite_add_solv(mol,chrg,uhf,tblite,smodel,solvent)
     class(tblite_solvation_type),allocatable :: solv
     type(solvation_input),allocatable :: solv_inp
     type(solvent_data) :: solv_data
-    character(len=:),allocatable :: str,solvdum
+    type(alpb_input)  :: alpb_tmp
+    type(cds_input)   :: cds_tmp
+    type(shift_input) :: shift_tmp
+    character(len=:),allocatable :: str,solvdum,method
     logical :: pr
 
     if (.not.allocated(smodel).or..not.allocated(solvent)) then
@@ -216,10 +220,16 @@ 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 .and. tblite%lvl.ne.xtblvl%param) 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
+    select case (tblite%lvl)
+    case (xtblvl%gfn1)
+      method = 'gfn1'
+    case (xtblvl%gfn2)
+      method = 'gfn2'
+    end select
 
 !>--- make an mctcmol object from mol
     call tblite_mol2mol(mol,chrg,uhf,mctcmol)
@@ -240,32 +250,77 @@ subroutine tblite_add_solv(mol,chrg,uhf,tblite,smodel,solvent)
     select case (trim(smodel))
     case ('gbsa')
       if (pr) call tblite%ctx%message("tblite> using GBSA/"//solvdum)
-      allocate (solv_inp%alpb)
-      solv_inp%alpb = alpb_input(solv_data%eps,alpb=.false.)
+      alpb_tmp%dielectric_const = solv_data%eps
+      alpb_tmp%alpb=.false.
+      alpb_tmp%method=method
+      alpb_tmp%solvent=solv_data%solvent
+      alpb_tmp%xtb=.true.
+      allocate (solv_inp%alpb, source=alpb_tmp)
+      cds_tmp%alpb=.false.
+      cds_tmp%solvent=solv_data%solvent
+      cds_tmp%method=method 
+      allocate (solv_inp%cds, source=cds_tmp)
+      shift_tmp%alpb=.false.
+      shift_tmp%solvent=solv_data%solvent
+      shift_tmp%method=method
+      allocate (solv_inp%shift, source=shift_tmp)
     case ('cpcm')
       if (pr) call tblite%ctx%message("tblite> using CPCM/"//solvdum)
       allocate (solv_inp%cpcm)
       solv_inp%cpcm = cpcm_input(solv_data%eps)
     case ('alpb')
       if (pr) call tblite%ctx%message("tblite> using ALPB/"//solvdum)
-      allocate (solv_inp%alpb)
-      solv_inp%alpb = alpb_input(solv_data%eps,alpb=.true.)
+      alpb_tmp%dielectric_const = solv_data%eps
+      alpb_tmp%alpb=.true.
+      alpb_tmp%method=method
+      alpb_tmp%solvent=solv_data%solvent
+      alpb_tmp%xtb=.true.
+      allocate (solv_inp%alpb, source=alpb_tmp)
+      cds_tmp%alpb=.true.
+      cds_tmp%solvent=solv_data%solvent
+      cds_tmp%method=method 
+      allocate (solv_inp%cds, source=cds_tmp)
+      shift_tmp%alpb=.true.
+      shift_tmp%solvent=solv_data%solvent
+      shift_tmp%method=method
+      allocate (solv_inp%shift, source=shift_tmp)
     case default
       if (pr) call tblite%ctx%message("tblite> Unknown tblite implicit solvation model!")
       return
     end select
+
     str = 'tblite> WARNING: implicit solvation energies are not entirely '// &
     &'consistent with the xtb implementation.'
     if (pr) call tblite%ctx%message(str)
 
-!>--- add to calculator
-    call new_solvation(solv,mctcmol,solv_inp,error)
+!>--- add electrostatic (Born part) to calculator
+    call new_solvation(solv,mctcmol,solv_inp,error,method)
     if (allocated(error)) then
       if (pr) call tblite%ctx%message("tblite> failed to set up tblite implicit solvation!")
       return
     end if
     call move_alloc(solv,cont)
     call tblite%calc%push_back(cont)
+!>--- add hbond and dispersion pert to calculator
+    if (allocated(solv_inp%cds)) then
+      block
+        class(tblite_solvation_type),allocatable :: cds
+        call new_solvation_cds(cds,mctcmol,solv_inp,error,method)
+        if (allocated(error)) return
+        call move_alloc(cds,cont)
+        call tblite%calc%push_back(cont)
+      end block
+    end if
+!>--- add gsolv shift to calculator
+    if (allocated(solv_inp%shift)) then
+      block
+        class(tblite_solvation_type),allocatable :: shift
+        call new_solvation_shift(shift,solv_inp,error,method)
+        if (allocated(error)) return
+        call move_alloc(shift,cont)
+        call tblite%calc%push_back(cont)
+      end block
+    end if
 
     deallocate (solv_inp)
 
@@ -302,11 +357,11 @@ subroutine tblite_singlepoint(mol,chrg,uhf,tblite,energy,gradient,iostatus)
     energy = 0.0_wp
     gradient(:,:) = 0.0_wp
     pr = (tblite%ctx%verbosity > 0)
-    if(tblite%ctx%verbosity>1)then
+    if (tblite%ctx%verbosity > 1) then
       verbosity = tblite%ctx%verbosity
     else
       verbosity = 0
-    endif
+    end if
 
 !>--- make an mctcmol object from mol
     call tblite_mol2mol(mol,chrg,uhf,mctcmol)
@@ -320,8 +375,8 @@ subroutine tblite_singlepoint(mol,chrg,uhf,tblite,energy,gradient,iostatus)
     case (xtblvl%ceh)
       call ceh_singlepoint(tblite%ctx,tblite%calc,mctcmol,tblite%wfn, &
       &              tblite%accuracy,verbosity)
-    case(xtblvl%eeq)
-      call eeq_guess(mctcmol, tblite%calc, tblite%wfn)
+    case (xtblvl%eeq)
+      call eeq_guess(mctcmol,tblite%calc,tblite%wfn)
     end select
 
     if (tblite%ctx%failed()) then
@@ -409,25 +464,25 @@ subroutine tblite_getwbos(tblite,nat,wbo)
     case default
 
       nao = tblite%calc%bas%nao
-      allocate(Pa(nao,nao),Pb(nao,nao))
-      call split_foccab(nao,tblite%wfn%focc, tblite%wfn%nel(1), tblite%wfn%nel(2), &
-      & focca, foccb)
+      allocate (Pa(nao,nao),Pb(nao,nao))
+      call split_foccab(nao,tblite%wfn%focc,tblite%wfn%nel(1),tblite%wfn%nel(2), &
+      & focca,foccb)
       call density_matrix(nao,focca,tblite%wfn%coeff(:,:,1),Pa)
       call density_matrix(nao,foccb,tblite%wfn%coeff(:,:,1),Pb)
-      call get_wbo(nat, nao, Pa,Pb, tblite%res%overlap, tblite%calc%bas%ao2at, wbo)
+      call get_wbo(nat,nao,Pa,Pb,tblite%res%overlap,tblite%calc%bas%ao2at,wbo)
 
     case (xtblvl%ceh)
-    !> no external access to the overlap in CEH, hence use the Wiberg BO with S=I
+      !> no external access to the overlap in CEH, hence use the Wiberg BO with S=I
       nao = tblite%calc%bas%nao
-      allocate(S(nao,nao), source=0.0_wp) 
-      do i=1,nao
+      allocate (S(nao,nao),source=0.0_wp)
+      do i = 1,nao
         S(i,i) = 1.0_wp
-      enddo
+      end do
       call get_wbo_rhf(nat,tblite%calc%bas%nao,tblite%wfn%density, &
       &                S,tblite%calc%bas%ao2at,wbo)
       wbo = wbo*2.0_wp !> somehow this is much better
 
-    case( xtblvl%eeq )
+    case (xtblvl%eeq)
       wbo = 0.0_wp
     end select
 #endif
@@ -477,37 +532,37 @@ 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
+  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
+    end if
+
+    call param%load_from_toml(table,error)
+
+    if (allocated(error)) then
+      io = 1
+    else
+      io = 0
+    end if
+    if (allocated(table)) deallocate (table)
+
+  end subroutine tblite_read_param_record
 #endif
 
 !========================================================================================!
diff --git a/src/ompmklset.F90 b/src/ompmklset.F90
index 75298733..6365cc38 100644
--- a/src/ompmklset.F90
+++ b/src/ompmklset.F90
@@ -105,9 +105,8 @@ subroutine new_ompautoset(env,modus,maxjobs,parallel_jobs,cores_per_job)
       call mkl_set_dynamic(0)
 #endif
     end if
-#ifdef WITH_OPENBLAS
-      call openblas_set_num_threads(1)
-#endif
+    call openblasset(1)
+
   case ('max')
     !> Both intern and environment variable threads to max
     parallel_jobs = T
@@ -137,9 +136,10 @@ subroutine new_ompautoset(env,modus,maxjobs,parallel_jobs,cores_per_job)
   call ompmklset(parallel_jobs)
   call ompenvset(cores_per_job)
 #ifdef WITH_OPENBLAS
-  if(modus.ne.'auto'.and.modus.ne.'auto_nested')then
-      call openblas_set_num_threads(cores_per_job)
-  endif
+!  if(modus.ne.'auto'.and.modus.ne.'auto_nested')then
+!      call openblasset(cores_per_job)
+!  endif
+   call openblasset(1)
 #endif
 end subroutine new_ompautoset
 
diff --git a/subprojects/tblite b/subprojects/tblite
index 1665c595..4556e28f 160000
--- a/subprojects/tblite
+++ b/subprojects/tblite
@@ -1 +1 @@
-Subproject commit 1665c5959a588a3b2388b420baca6e9198773cbd
+Subproject commit 4556e28f391573b3c80c94beb7c56313005d5269
diff --git a/subprojects/tblite.wrap b/subprojects/tblite.wrap
index 261920a9..bb401a77 100644
--- a/subprojects/tblite.wrap
+++ b/subprojects/tblite.wrap
@@ -1,4 +1,4 @@
 [wrap-git]
 directory = tblite
-url = https://github.com/tblite/tblite
-revision = main
+url = https://github.com/ppracht/tblite
+revision = xtb_solvation
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 4b06e776..f1e0acc6 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -11,6 +11,7 @@ set(
 set(
   test-srcs
   "testmol.f90"
+  "helpers.f90"
   "main.f90"
 )
 foreach(t IN LISTS tests)
diff --git a/test/helpers.f90 b/test/helpers.f90
new file mode 100644
index 00000000..dbaa4b8a
--- /dev/null
+++ b/test/helpers.f90
@@ -0,0 +1,87 @@
+module test_helpers
+  use testdrive,only:error_type,check,test_failed
+  use crest_parameters
+  use crest_testmol
+  use crest_calculator
+  use strucrd
+  implicit none
+  private
+
+  public :: test_e,test_g
+
+!========================================================================================!
+!========================================================================================!
+contains  !> Helper routines for testing
+!========================================================================================!
+!========================================================================================!
+
+  subroutine test_e(mol,calc,error,ref)
+!***********************************************
+!* calculate energy and cross-check with
+!* a given reference
+!***********************************************
+    !> Error handling
+    type(error_type),allocatable,intent(out) :: error
+    !> calculation data
+    type(calcdata),intent(inout) :: calc
+    !> molecule
+    type(coord),intent(inout) :: mol
+    !> reference energy
+    real(wp),intent(in) :: ref
+    !> LOCAL
+    real(wp),parameter :: thr = sqrt(epsilon(1.0_wp))
+    real(wp),allocatable :: gradient(:,:)
+    integer :: io
+    real(wp) :: energy
+
+    energy=0.0_wp
+
+    !> calculation
+    allocate (gradient(3,mol%nat),source=0.0_wp)
+    call engrad(mol,calc,energy,gradient,io)
+
+    if (abs(energy-ref) > thr) then
+      !call test_failed(error,"energy does not match reference")
+      call calc%info(stdout)
+      call check(error,energy,ref,thr=thr)
+    end if
+
+  end subroutine test_e
+
+  subroutine test_g(mol,calc,error)
+!***********************************************
+!* calculate numerical gradient and cross-check
+!* with the one returned by calculator
+!***********************************************
+    !> Error handling
+    type(error_type),allocatable,intent(out) :: error
+    !> calculation data
+    type(calcdata),intent(inout) :: calc
+    !> molecule
+    type(coord),intent(inout) :: mol
+    !> LOCAL
+    real(wp),parameter :: thr = 1.0E-6_wp
+    real(wp),allocatable :: gradient(:,:),numg(:,:)
+    integer :: io
+    real(wp) :: energy
+
+    !> calculation
+    allocate (gradient(3,mol%nat),source=0.0_wp)
+    call engrad(mol,calc,energy,gradient,io)
+
+    call numgrad(mol,calc,numgradient=numg)
+
+    if (any(abs(gradient-numg) > thr)) then
+      call test_failed(error,"Gradient does not match")
+      call calc%info(stdout)
+      print '(3es20.13)',gradient
+      print '(a)',"---"
+      print '(3es20.13)',numg
+      print '(a)',"---"
+      print '(3es20.13)',gradient-numg
+    end if
+  end subroutine test_g
+
+!========================================================================================!
+!========================================================================================!
+end module test_helpers
diff --git a/test/test_tblite.F90 b/test/test_tblite.F90
index ec6225d0..fecefd74 100644
--- a/test/test_tblite.F90
+++ b/test/test_tblite.F90
@@ -4,6 +4,7 @@ module test_tblite
   use crest_calculator
   use strucrd
   use crest_testmol, only: get_testmol
+  use test_helpers, only: test_e, test_g
   implicit none
   private
 
@@ -323,64 +324,23 @@ subroutine test_gfn2_sp_alpb(error)
     real(wp) :: energy
     real(wp),allocatable :: grad(:,:)
     integer :: io
-!&<
-    real(wp),parameter :: e_ref = -42.150818113966622_wp
-    real(wp),parameter :: g_ref(3,24) = reshape([&
-    &   -0.004041852805369_wp, -0.001459116937018_wp,  0.000004212952013_wp, & 
-    &    0.005431782185647_wp, -0.026769409144275_wp, -0.000053933849249_wp, &
-    &    0.005232547908947_wp,  0.019316152153876_wp,  0.000051828956087_wp, &
-    &    0.003984299183311_wp,  0.004654018323578_wp, -0.000031672641741_wp, &
-    &   -0.034791239964625_wp, -0.009664440508457_wp, -0.000039268100546_wp, &
-    &    0.006792385075909_wp, -0.003127178836761_wp,  0.000032791144117_wp, &
-    &    0.015306198183764_wp,  0.005195424722641_wp,  0.000035383591075_wp, &
-    &   -0.011644038936970_wp,  0.009526167376125_wp, -0.000103656945404_wp, &
-    &   -0.001073718470205_wp, -0.010467260905049_wp,  0.000071170633309_wp, &
-    &   -0.012208257598543_wp,  0.007792225857756_wp, -0.000026835529946_wp, &
-    &    0.024297093142032_wp, -0.020057317289964_wp, -0.000036190038701_wp, &
-    &    0.002288207714506_wp,  0.015784323612706_wp, -0.000015109378387_wp, &
-    &   -0.001839503827186_wp, -0.003557956433963_wp,  0.000003971003294_wp, &
-    &   -0.006475695156935_wp,  0.007869132396390_wp, -0.000061014145496_wp, &
-    &   -0.001323378271276_wp,  0.002797437673167_wp,  0.000004384755727_wp, &
-    &    0.002794743156382_wp,  0.000752555767144_wp,  0.003757718239733_wp, &
-    &    0.002796550406646_wp,  0.000753562108628_wp, -0.003760695667456_wp, &
-    &   -0.003950802127381_wp,  0.008247313874474_wp,  0.000028580331243_wp, &
-    &    0.007699255283409_wp,  0.001924889526582_wp,  0.000000878165662_wp, &
-    &    0.000911988480912_wp, -0.000170582007313_wp,  0.003208806647171_wp, &
-    &    0.000936315963916_wp, -0.000147196304551_wp, -0.003212846696159_wp, &
-    &   -0.002208060498652_wp, -0.009259192346010_wp,  0.000102993956974_wp, &
-    &    0.000531577240056_wp,  0.000047927691044_wp,  0.002360617394228_wp, &
-    &    0.000553603731705_wp,  0.000018519629250_wp, -0.002322114777547_wp &
-    & ], shape(g_ref))
-!&>
+    real(wp),parameter :: e_ref = -23.983234299793384_wp
 
     !> setup
     call sett%create('gfn2')
     sett%solvmodel = 'alpb'
     sett%solvent = 'water'
     call calc%add(sett)
-    call get_testmol('caffeine',mol)
-    allocate (grad(3,mol%nat))
+    call get_testmol('cytosine',mol)
 
-    !> calculation
-    call engrad(mol,calc,energy,grad,io)
-    !write(*,'(F25.15)') energy
-    !write(*,'(3F25.15)') grad
-    call check(error,io,0)
+    !> energy
+    call test_e(mol,calc,error,e_ref)
     if (allocated(error)) return
 
-    call check(error,energy,e_ref,thr=1e-7_wp)
+    !> gradient
+    call test_g(mol,calc,error)
     if (allocated(error)) return
-
-    if (any(abs(grad-g_ref) > thr2)) then
-      call test_failed(error,"Gradient of energy does not match")
-      print'(3es21.14)',grad
-      print'("---")'
-      print'(3es21.14)',g_ref
-      print'("---")'
-      print'(3es21.14)',grad-g_ref
-    end if
-
-    deallocate (grad)
+ 
   end subroutine test_gfn2_sp_alpb
 
 !========================================================================================!
diff --git a/test/testmol.f90 b/test/testmol.f90
index e380085c..452fef0c 100644
--- a/test/testmol.f90
+++ b/test/testmol.f90
@@ -52,6 +52,28 @@ module crest_testmol
     &  shape(methane_xyz))
 !&>
 
+
+!&<
+  !> distorted methane
+  integer,parameter ::  cytosine_nat = 13
+  integer, parameter :: cytosine_at(cytosine_nat) = [8,6,7,6,6,6,7,7,1,1,1,1,1]
+  real(wp),parameter :: cytosine_xyz(3,cytosine_nat) = reshape(&
+    & [-4.142195348814,  2.008151058305, -0.018151734373, & 
+    &  -2.160173249661,  0.802767457051, -0.000058193762, &  
+    &  -2.268542712629, -1.794135272435,  0.003855688769, &
+    &  -0.191488773092, -3.309463007997, -0.005095763012, &
+    &   2.109332653259, -2.284418067990, -0.007062235327, &
+    &   2.185078861327,  0.516677179272,  0.028019148471, &
+    &   4.467560431809,  1.727180498582,  0.071380211323, &
+    &   0.177657703530,  1.926043596048,  0.012000916144, &
+    &  -4.036103295159, -2.521391305153,  0.014335396651, &
+    &  -0.507864921184, -5.333664751031, -0.008697813741, &
+    &   3.837497023245, -3.367210871956, -0.007065809214, &
+    &   5.961039176371,  0.856469245914, -0.742491063662, &
+    &   4.318834327087,  3.589804055892, -0.368535366892],&
+    &  shape(cytosine_xyz))
+!&>
+
 contains
 
   subroutine get_testmol(name,mol)
@@ -66,6 +88,14 @@ subroutine get_testmol(name,mol)
       mol%xyz(:,:) = methane_xyz(:,:)
       mol%chrg = 0
       mol%uhf = 0
+    case ('cytosine')
+      mol%nat = cytosine_nat
+      allocate (mol%at(mol%nat))
+      mol%at(:) = cytosine_at(:)
+      allocate (mol%xyz(3,mol%nat))
+      mol%xyz(:,:) = cytosine_xyz(:,:)
+      mol%chrg = 0
+      mol%uhf = 0
     case default
       !> caffeine
       mol%nat = caffeine_nat