diff --git a/.github/workflows/fortran-build.yml b/.github/workflows/fortran-build.yml
index d8c9c6045..55300bc1d 100644
--- a/.github/workflows/fortran-build.yml
+++ b/.github/workflows/fortran-build.yml
@@ -224,7 +224,7 @@ jobs:
             DESTDIR: ${{ env.PWD }}/xtb-bleed
          - name: Upload binary
            if: github.event_name == 'push'
-           uses: actions/upload-artifact@v2
+           uses: actions/upload-artifact@v4
            with:
             name: xtb-bleed.tar.xz
             path: xtb-bleed.tar.xz
diff --git a/src/hessian.F90 b/src/hessian.F90
index abf6d181b..c95d2bbb7 100644
--- a/src/hessian.F90
+++ b/src/hessian.F90
@@ -25,7 +25,7 @@ module xtb_hessian
 
    public :: numhess
    public :: trproj, rdhess, g98fake2, distort, write_tm_vibspectrum
-
+   public :: rescale_freq
 contains
 
 subroutine numhess( &
@@ -102,7 +102,6 @@ subroutine numhess( &
    real(wp),allocatable :: fc_tb(:)
    real(wp),allocatable :: fc_bias(:)
    real(wp),allocatable :: v(:)
-   real(wp),allocatable :: fc_tmp(:)
    real(wp),allocatable :: freq_scal(:)
    real(wp),allocatable :: aux (:)
    real(wp),allocatable :: isqm(:)
@@ -129,7 +128,7 @@ subroutine numhess( &
    allocate(hss(n3*(n3+1)/2),hsb(n3*(n3+1)/2),h(n3,n3),htb(n3,n3),hbias(n3,n3), &
       & gl(3,mol%n),isqm(n3),xyzsave(3,mol%n),dipd(3,n3), &
       & pold(n3),nb(20,mol%n),indx(mol%n),molvec(mol%n),bond(mol%n,mol%n), &
-      & v(n3),fc_tmp(n3),freq_scal(n3),fc_tb(n3),fc_bias(n3),amass(n3), h_dummy(n3,n3))
+      & freq_scal(n3),fc_tb(n3),fc_bias(n3),amass(n3), h_dummy(n3,n3))
 
    if (set%elprop == p_elprop_alpha) then
       allocate(dalphadr(6,n3), source = 0.0_wp)
@@ -355,6 +354,14 @@ subroutine numhess( &
    write(env%unit,'("writing file <",a,">.")') hname
    call wrhess(n3,hss,hname)
 
+   ! non mass weigthed biased Hessian in hsb
+   if (set%runtyp .eq. p_run_bhess) then
+      hname = 'hessian_sph'
+      write (env%unit, '(a)')
+      write (env%unit, '("writing file <",a,">.")') hname
+      call wrhess(n3, hsb, hname)
+   end if
+
    ! include masses
    k=0
    do i=1,n3
@@ -386,26 +393,13 @@ subroutine numhess( &
       return
    end if
 
-   ! calculate fc_tb and fc_bias
-   alp1=1.27_wp
-   alp2=1.5d-4
+   
+
    if (set%runtyp.eq.p_run_bhess) then
-      do j=1,n3
-         v(1:n3) = res%hess(1:n3,j) ! modes
-         call mctc_gemv(htb,v,fc_tmp)
-         fc_tb(j) = mctc_dot(v,fc_tmp)
-         call mctc_gemv(hbias,v,fc_tmp)
-         fc_bias(j) = mctc_dot(v,fc_tmp)
-         if (abs(res%freq(j)).gt.1.0d-6) then
-            freq_scal(j) = sqrt( (fc_tb(j)+alp2) / ( (fc_tb(j)+alp2) +  alp1*fc_bias(j) ) )
-            if (fc_tb(j).lt.0.and.fc_bias(j).ne.0) then
-               freq_scal(j) = -sqrt( (abs(fc_tb(j))+alp2) / ( (abs(fc_tb(j))+alp2) + alp1*fc_bias(j) ) )
-            end if
-         else
-            freq_scal(j) = 1.0_wp
-         end if
-      end do
-   end if
+      call rescale_freq(n3,htb,res%hess,hbias,res%freq,fc_tb,fc_bias,freq_scal)
+   else
+      freq_scal(1:n3) = 1.0_wp
+   end if 
 
    write(env%unit,'(a)')
    if(res%linear)then
@@ -721,6 +715,43 @@ subroutine distort(mol,freq,u)
 
 end subroutine distort
 
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+subroutine rescale_freq(n3,htb,hess,hbias,freq,fc_tb,fc_bias,freq_scal)
+   use xtb_mctc_blas
+   implicit none
+   real(wp),intent(in) :: htb (n3,n3)
+   real(wp),intent(in) :: hess(n3,n3)
+   real(wp),intent(in) :: hbias(n3,n3)
+   real(wp),intent(in) :: freq(n3)
+   real(wp), intent(out) :: fc_tb(n3),fc_bias(n3)
+   real(wp), intent(out) :: freq_scal(n3) 
+   real(wp),allocatable :: v(:)
+   real(wp),allocatable :: fc_tmp(:)
+   real(wp), parameter :: alp1=1.27_wp, alp2=1.5e-4_wp
+   integer, intent(in) :: n3
+   integer :: j
+   
+   allocate(fc_tmp(n3),v(n3))
+   ! calculate fc_tb and fc_bias
+   do j=1,n3
+      v(1:n3) = hess(1:n3,j) ! modes
+      call mctc_gemv(htb,v,fc_tmp)
+      fc_tb(j) = mctc_dot(v,fc_tmp)
+      call mctc_gemv(hbias,v,fc_tmp)
+      fc_bias(j) = mctc_dot(v,fc_tmp)
+      if (abs(freq(j)) .gt. 1.0e-6_wp) then
+         freq_scal(j) = sqrt( (fc_tb(j)+alp2) / ( (fc_tb(j)+alp2) +  alp1*fc_bias(j) ) )
+         if (fc_tb(j) .lt. 0.0_wp .and. fc_bias(j) .ne. 0.0_wp) then
+            freq_scal(j) = -sqrt( (abs(fc_tb(j))+alp2) / ( (abs(fc_tb(j))+alp2) + alp1*fc_bias(j) ) )
+         end if
+      else
+         freq_scal(j) = 1.0_wp   
+      end if
+   end do
+end subroutine rescale_freq
+   
+
 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/src/prog/thermo.f90 b/src/prog/thermo.f90
index 43c6d935d..348acd63b 100644
--- a/src/prog/thermo.f90
+++ b/src/prog/thermo.f90
@@ -17,312 +17,375 @@
 
 !> Thermodynamic functions and hessian post processing
 module xtb_prog_thermo
-   use xtb_mctc_accuracy, only : wp
-   use xtb_mctc_convert, only : autorcm, autoamu
-   use xtb_mctc_filetypes, only : getFileType, hessType
+   use xtb_mctc_blas, only: mctc_gemv, mctc_dot
+   use xtb_mctc_accuracy, only: wp
+   use xtb_mctc_convert, only: autorcm, autoamu, amutoau
+   use xtb_mctc_filetypes, only: getFileType, hessType
    use xtb_mctc_timings
-   use xtb_mctc_version, only : version, author, date
-   use xtb_io_reader, only : readMolecule, readHessian
-   use xtb_type_environment, only : TEnvironment
-   use xtb_type_molecule, only : TMolecule
-   use xtb_type_reader, only : TReader
-   use xtb_prog_argparser, only : TArgParser
-   use xtb_freq_project, only : projectHessian
-   use xtb_freq_utils, only : massWeightHessian, diagHessian
-   use xtb_propertyoutput, only : print_thermo
-   use xtb_splitparam, only : atmass
-   use xtb_setmod, only : set_thermo, set_symmetry
+   use xtb_mctc_version, only: version, author, date
+   use xtb_io_reader, only: readMolecule, readHessian
+   use xtb_type_environment, only: TEnvironment
+   use xtb_type_molecule, only: TMolecule
+   use xtb_type_reader, only: TReader
+   use xtb_prog_argparser, only: TArgParser
+   use xtb_freq_project, only: projectHessian
+   use xtb_freq_utils, only: massWeightHessian, diagHessian
+   use xtb_propertyoutput, only: print_thermo
+   use xtb_splitparam, only: atmass
+   use xtb_setmod, only: set_thermo, set_symmetry
+   use xtb_freq_io, only: rdhess
    implicit none
    private
 
    public :: xtbThermo
 
-
 contains
 
-
 !> Main program of the thermo submodule
-subroutine xtbThermo(env, argParser)
-
-   !> Source of errors in the main program unit
-   character(len=*), parameter :: source = "prog_thermo"
-
-   !> Calculation environment
-   type(TEnvironment), intent(inout) :: env
-
-   !> Commandline argument parser
-   type(TArgParser), intent(inout) :: argParser
+   subroutine xtbThermo(env, argParser)
 
-   !> Molecular structure data
-   type(TMolecule) :: mol
+      !> Source of errors in the main program unit
+      character(len=*), parameter :: source = "prog_thermo"
 
-   type(TReader) :: reader
-   character(len=:), allocatable :: file
-   real(wp), allocatable :: freq(:), hessian(:, :)
-   real(wp) :: htot, gtot, zp
-   integer :: nFiles, ftype, hFormat, nimag
-   logical :: massWeighted
+      !> Calculation environment
+      type(TEnvironment), intent(inout) :: env
 
-   call parseArguments(env, argParser, hFormat, massWeighted)
+      !> Commandline argument parser
+      type(TArgParser), intent(inout) :: argParser
 
-   nFiles = argParser%countFiles()
-   if (nFiles == 0) then
-      call env%error("No input file given, so there is nothing to do", source)
-      call thermoHelp(env%unit)
-   end if
-   if (nFiles > 2) then
-      call env%error("Multiple input files present, aborting", source)
-   end if
+      !> Molecular structure data
+      type(TMolecule) :: mol
 
-   call env%checkpoint("Command line argument parsing failed")
+      type(TReader) :: reader
+      character(len=:), allocatable :: file
+      real(wp), allocatable :: freq(:), hessian(:, :)
+      real(wp) :: htot, gtot, zp
+      integer :: nFiles, ftype, hFormat, nimag
+      logical :: massWeighted, bhess
 
-   call argParser%nextFile(file)
+      call parseArguments(env, argParser, hFormat, massWeighted, bhess)
 
-   !> Determine the file type for processing in the reader
-   ftype = getFileType(file)
+      nFiles = argParser%countFiles()
+      if (nFiles == 0) then
+         call env%error("No input file given, so there is nothing to do", source)
+         call thermoHelp(env%unit)
+      end if
+      if (nFiles > 2 .and. .not. bhess) then
+         call env%error("Multiple input files present, aborting", source)
+      end if
 
-   !> Generate the reader by connecting the file to the instance
-   call reader%open(file)
-   !> Process the file to obtain the molecular structure data
-   call readMolecule(env, mol, reader%unit, ftype)
-   !> Close the file, we are done with it here
-   call reader%close
+      if (nFiles > 3) then
+         call env%error("Multiple input files present, aborting", source)
+      end if
 
-   call env%checkpoint("Could not read geometry from '"//file//"'")
+      call env%checkpoint("Command line argument parsing failed")
 
-   !> Print an informative banner
-   call thermoHeader(env%unit)
-   !> print current time
-   call prdate('S')
+      call argParser%nextFile(file)
 
-   !> Also store a copy of the atomic masses in amu in global storage
-   atmass = mol%atmass * autoamu
+      !> Determine the file type for processing in the reader
+      ftype = getFileType(file)
 
-   call argParser%nextFile(file)
-   allocate(freq(3*mol%n))
-   if (allocated(file)) then
-      allocate(hessian(3*mol%n, 3*mol%n))
+      !> Generate the reader by connecting the file to the instance
       call reader%open(file)
-      call readHessian(env, mol, hessian, reader, format=hFormat)
+      !> Process the file to obtain the molecular structure data
+      call readMolecule(env, mol, reader%unit, ftype)
+      !> Close the file, we are done with it here
       call reader%close
 
-      if (.not.massWeighted) then
-         call massWeightHessian(hessian, mol%atmass)
-      end if
+      call env%checkpoint("Could not read geometry from '"//file//"'")
 
-      if (hFormat == hessType%dftbplus .or. hFormat == hessType%orca) then
-         call projectHessian(hessian, mol, .true., .true.)
-      end if
-
-      call diagHessian(env, hessian, freq)
-
-      call preigf(env%unit, freq, autorcm)
+      !> Print an informative banner
+      call thermoHeader(env%unit)
+      !> print current time
+      call prdate('S')
 
-      call env%checkpoint("Could not read hessian from '"//file//"'")
-   else
-      freq(:) = -1.0_wp
-   end if
-   freq(:) = freq * autorcm
+      !> Also store a copy of the atomic masses in amu in global storage
+      atmass = mol%atmass * autoamu
 
-   call print_thermo(env%unit,mol%n, size(freq), mol%at, mol%xyz, freq, 0.0_wp, &
-      & htot, gtot, nimag, .true., zp)
+      call argParser%nextFile(file)
+      allocate (freq(3 * mol%n))
+      if (allocated(file)) then
+         allocate (hessian(3 * mol%n, 3 * mol%n))
+         call reader%open(file)
+         call readHessian(env, mol, hessian, reader, format=hFormat)
+         call reader%close
 
-   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 (bhess) then
+            call remove_sphshift(env, argParser, mol, hessian, freq) ! frequencies have to be rescaled
+         else
+            if (.not. massWeighted) then
+               call massWeightHessian(hessian, mol%atmass)
+            end if
 
-   write(env%unit,'(a)')
-   call terminate(0)
+            if (hFormat == hessType%dftbplus .or. hFormat == hessType%orca) then
+               call projectHessian(hessian, mol, .true., .true.)
+            end if
 
-end subroutine xtbThermo
+            call diagHessian(env, hessian, freq)
 
+            call preigf(env%unit, freq, autorcm)
 
-subroutine preigf(unit, freq, scale)
-   real(wp), intent(in) :: freq(:)
-   real(wp), intent(in) :: scale
-   integer, intent(in) :: unit
-   integer, parameter :: nRows = 6
-   integer :: ntimes, nrest, i, j, n2, k
-   ntimes = size(freq)/nRows
-   nrest = mod(size(freq), nRows)
-   if (ntimes.eq.0) nrest = size(freq)
-   j = 1
-   n2 = nRows
-   do k = 1, ntimes
-      write(unit,'("eigval :",3x,10f9.2)') (freq(i)*scale,i=j,n2)
-      j = j + nRows
-      n2 = n2 + nRows
-   end do
-   if (nrest.gt.0.or.ntimes.eq.0) then
-      write(unit,'("eigval :",3x,10f9.2)') (freq(i)*scale,i=j,j+nrest-1)
-   end if
+            call env%checkpoint("Could not read hessian from '"//file//"'")
+         end if
+      else
+         freq(:) = -1.0_wp
+      end if
+      freq(:) = freq * autorcm
+
+      call print_thermo(env%unit, mol%n, size(freq), mol%at, mol%xyz, freq, 0.0_wp, &
+         & htot, gtot, nimag, .true., zp)
+
+      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')
+
+      write (env%unit, '(a)')
+      call terminate(0)
+
+   end subroutine xtbThermo
+
+   subroutine preigf(unit, freq, scale)
+      real(wp), intent(in) :: freq(:)
+      real(wp), intent(in) :: scale
+      integer, intent(in) :: unit
+      integer, parameter :: nRows = 6
+      integer :: ntimes, nrest, i, j, n2, k
+      ntimes = size(freq) / nRows
+      nrest = mod(size(freq), nRows)
+      if (ntimes == 0) nrest = size(freq)
+      j = 1
+      n2 = nRows
+      do k = 1, ntimes
+         write (unit, '("eigval :",3x,10f9.2)') (freq(i) * scale, i=j, n2)
+         j = j + nRows
+         n2 = n2 + nRows
+      end do
+      if (nrest > 0 .or. ntimes == 0) then
+         write (unit, '("eigval :",3x,10f9.2)') (freq(i) * scale, i=j, j + nrest - 1)
+      end if
 
-end subroutine preigf
+   end subroutine preigf
+
+! Remove the frequency shift caused by the bias RMSD employed in the SPH calculation:
+   subroutine remove_sphshift(env, argParser, mol, hessian, freq)
+      use xtb_hessian
+      real(wp), intent(inout) :: freq(:), hessian(:, :)
+      real(wp), allocatable :: freq_sph(:), hessian_sph(:, :) ! biased hessian from SPH
+      real(wp), allocatable :: htb(:, :), v(:), fc_tmp(:), fc_tb(:), fc_bias(:), freq_scal(:)
+      real(wp) :: alp1, alp2 ! factors for SPH scaling
+      integer :: i, j, n3
+      type(TEnvironment), intent(inout) :: env
+      !> Commandline argument parser
+      type(TArgParser), intent(inout) :: argParser
+      type(TReader) :: reader
+      type(TMolecule) :: mol
+      character(len=:), allocatable :: file
+      integer :: hFormat
+      logical :: ex
+
+      n3 = 3 * mol%n
+      hFormat = 1 ! hessian_sph must stem from xtb SPH calculation!
+      call argParser%nextFile(file)
+      write (env%unit, '(a)') 'Reading hessian_sph file: '//file
+      write (env%unit, '(a)') 'Removing the frequency shift caused by the bias RMSD:'
+      allocate (freq_sph(n3))
+      if (allocated(file)) then
+         allocate (hessian_sph(n3, n3))
+         call reader%open(file)
+         call readHessian(env, mol, hessian_sph, reader, format=hFormat)
+         call reader%close
+         call massWeightHessian(hessian_sph, atmass) ! keep units similar to hessian.f90 routines in main program
+         call massWeightHessian(hessian, atmass) ! keep units similar to hessian.f90 routines in main program
+         htb = hessian - hessian_sph
+
+         call diagHessian(env, hessian, freq)
+
+         allocate (v(n3), fc_tmp(n3), freq_scal(n3), fc_tb(n3), fc_bias(n3))
+
+         call rescale_freq(n3, htb, hessian, hessian_sph, freq, fc_tb, fc_bias, freq_scal)
+
+         ! scale frequencies and convert now to atomic units
+         do j = 1, n3
+            freq(j) = freq(j) * freq_scal(j) / sqrt(amutoau)
+         end do
+
+         call preigf(env%unit, freq, autorcm)
+      else
+         print *, "no hessian_sph file found, --bhess mode can not be used"
+         call terminate(1)
+      end if
 
+   end subroutine remove_sphshift
 
 !> Parse command line arguments
-subroutine parseArguments(env, args, htype, massWeighted)
+   subroutine parseArguments(env, args, htype, massWeighted, bhess)
 
-   !> Name of error producer
-   character(len=*), parameter :: source = "prog_thermo_parseArguments"
+      !> Name of error producer
+      character(len=*), parameter :: source = "prog_thermo_parseArguments"
 
-   !> Calculation environment
-   type(TEnvironment) :: env
+      !> Calculation environment
+      type(TEnvironment) :: env
 
-   !> Command line argument parser
-   type(TArgParser) :: args
+      !> Command line argument parser
+      type(TArgParser) :: args
 
-   !> File type of the hessian
-   integer, intent(out) :: htype
+      !> File type of the hessian
+      integer, intent(out) :: htype
 
-   !> Hessian is already mass weighted
-   logical, intent(out) :: massWeighted
+      !> Hessian is already mass weighted
+      logical, intent(out) :: massWeighted
 
-   integer :: nFlags
-   character(len=:), allocatable :: flag, sec
+      !> Hessian stems from a biased SPH calculation
+      logical, intent(out) :: bhess
 
-   htype = hessType%tmol
-   massWeighted = .false.
+      integer :: nFlags
+      character(len=:), allocatable :: flag, sec
 
-   nFlags = args%countFlags()
-   call args%nextFlag(flag)
-   do while(allocated(flag))
-      if (len(flag) > 2 .and. flag(1:1) == '-' .and. flag(1:2) /= '--') then
-         call env%warning("the use of '"//flag//"' is discouraged, "// &
-            & "please use '-"//flag//"' next time", source)
-         flag = '-'//flag
-      end if
-      select case(flag)
-      case default
-         call env%warning("Unknown option '"//flag//"' provided", source)
+      htype = hessType%tmol
+      massWeighted = .false.
+      bhess = .false.
 
-      case('--help', '-h')
-         call thermoHelp(env%unit)
-         call terminate(0)
-
-      case('--dftb+','--dftbplus')
-         htype = hessType%dftbplus
-
-      case('--turbomole')
-         htype = hessType%tmol
-         massWeighted = .true.
-     
-      case('--orca')
-         htype = hessType%orca
-
-      case('--scale')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_thermo(env, 'scale', sec)
-         else
-            call env%error("Freq. scaling factor missing", source)
-         end if
-
-      case('--sthr')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_thermo(env, 'sthr', sec)
-         else
-            call env%error("Rotor cutoff is missing", source)
-         end if
-     
-      case('--ithr')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_thermo(env, 'imagthr','-'//sec)
-         else
-            call env%error("Imaginary cutoff is missing", source)
-         end if
-      
-      case('--desy')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_symmetry(env,'desy',sec)
-         else
-            call env%error("Threshold for symmetry recognition missing", source)
-         end if
-
-      case('--temp')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_thermo(env, 'temp', sec)
-         else
-            call env%error("Temperature in K is missing", source)
-         end if
-
-      end select
+      nFlags = args%countFlags()
       call args%nextFlag(flag)
-   end do
-
-end subroutine parseArguments
-
+      do while (allocated(flag))
+         if (len(flag) > 2 .and. flag(1:1) == '-' .and. flag(1:2) /= '--') then
+            call env%warning("the use of '"//flag//"' is discouraged, "// &
+               & "please use '-"//flag//"' next time", source)
+            flag = '-'//flag
+         end if
+         select case (flag)
+         case default
+            call env%warning("Unknown option '"//flag//"' provided", source)
+
+         case ('--help', '-h')
+            call thermoHelp(env%unit)
+            call terminate(0)
+
+         case ('--dftb+', '--dftbplus')
+            htype = hessType%dftbplus
+
+         case ('--turbomole')
+            htype = hessType%tmol
+            massWeighted = .true.
+
+         case ('--orca')
+            htype = hessType%orca
+
+         case ('--scale')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               call set_thermo(env, 'scale', sec)
+            else
+               call env%error("Freq. scaling factor missing", source)
+            end if
+
+         case ('--sthr')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               call set_thermo(env, 'sthr', sec)
+            else
+               call env%error("Rotor cutoff is missing", source)
+            end if
+
+         case ('--ithr')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               call set_thermo(env, 'imagthr', '-'//sec)
+            else
+               call env%error("Imaginary cutoff is missing", source)
+            end if
+
+         case ('--desy')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               call set_symmetry(env, 'desy', sec)
+            else
+               call env%error("Threshold for symmetry recognition missing", source)
+            end if
+
+         case ('--temp')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               call set_thermo(env, 'temp', sec)
+            else
+               call env%error("Temperature in K is missing", source)
+            end if
+         case ('--bhess')
+            bhess = .true.
+
+         end select
+         call args%nextFlag(flag)
+      end do
+
+   end subroutine parseArguments
 
 !> Header for this submodule
-subroutine thermoHeader(unit)
-
-   !> IO unit
-   integer, intent(in) :: unit
-
-   write(unit,'(a)') &
-      "      -----------------------------------------------------------",&
-      "     |                   =====================                   |",&
-      "     |                        T H E R M O                        |",&
-      "     |                   =====================                   |",&
-      "     |                         S. Ehlert                         |",&
-      "     |          Mulliken Center for Theoretical Chemistry        |",&
-      "     |                    University of Bonn                     |",&
-      "      -----------------------------------------------------------",""
-   write(unit,'(3x,"*",*(1x,a))') &
-      & "xtb-thermo version", version, "compiled by", author, "on", date
-   write(unit,'(a)')
-
-end subroutine thermoHeader
-
-
-subroutine thermoHelp(unit)
-
-   !> IO unit
-   integer, intent(in) :: unit
-
-   write(unit, '(a)') &
-   "Usage: xtb thermo [options] <geometry> [hessian]", &
-   "",&
-   "<geometry> may be provided as any valid input to xtb", &
-   "the [hessian] file is read and processed depending on the selected options,", &
-   "the default format is a non-massweighted Turbomole hessian file", &
-   "",&
-   "Options",&
-   "",&
-   "   --sthr REAL         Rotor cutoff for RRHO partition function in rcm",&
-   "",&
-   "   --ithr REAL         Imag. freq. cutoff for RRHO in rcm",&
-   "                       enter a positive value, the sign will be inverted",&
-   "",&
-   "   --scale REAL        Frequency scaling factor for RRHO",&
-   "",&
-   "   --desy REAL         Threshold for symmetry recognition",&
-   "",&
-   "   --temp REAL[,...]   Temperature for thermodynamic functions in K,",&
-   "                       takes a comma separated list of temperatures",&
-   "",&
-   "   --dftbplus             Read a DFTB+ hessian.out file, implies projection",&
-   "",&
-   "   --turbomole         Read a Turbomole Hessian file",&
-   "                       use this only when $nomw is not present in control",&
-   "",&
-   "   --orca         Read an Orca Hessian file",&
-   "",&
-   "Output Conventions:",&
-   "",&
-   "total energies are given in atomic units (Eh), frequencies in reciprocal cm",&
-   "",&
-   "More information can be obtained at https://xtb-docs.readthedocs.io/",&
-   ""
-end subroutine thermoHelp
-
+   subroutine thermoHeader(unit)
+
+      !> IO unit
+      integer, intent(in) :: unit
+
+      write (unit, '(a)') &
+         "      -----------------------------------------------------------", &
+         "     |                   =====================                   |", &
+         "     |                        T H E R M O                        |", &
+         "     |                   =====================                   |", &
+         "     |                         S. Ehlert                         |", &
+         "     |          Mulliken Center for Theoretical Chemistry        |", &
+         "     |                    University of Bonn                     |", &
+         "      -----------------------------------------------------------", ""
+      write (unit, '(3x,"*",*(1x,a))') &
+         & "xtb-thermo version", version, "compiled by", author, "on", date
+      write (unit, '(a)')
+
+   end subroutine thermoHeader
+
+   subroutine thermoHelp(unit)
+
+      !> IO unit
+      integer, intent(in) :: unit
+
+      write (unit, '(a)') &
+         "Usage: xtb thermo [options] <geometry> [hessian]", &
+         "", &
+         "<geometry> may be provided as any valid input to xtb", &
+         "the [hessian] file is read and processed depending on the selected options,", &
+         "the default format is a non-massweighted Turbomole hessian file", &
+         "", &
+         "Options", &
+         "", &
+         "   --sthr REAL         Rotor cutoff for RRHO partition function in rcm", &
+         "", &
+         "   --ithr REAL         Imag. freq. cutoff for RRHO in rcm", &
+         "                       enter a positive value, the sign will be inverted", &
+         "", &
+         "   --scale REAL        Frequency scaling factor for RRHO", &
+         "", &
+         "   --desy REAL         Threshold for symmetry recognition", &
+         "", &
+         "   --temp REAL[,...]   Temperature for thermodynamic functions in K,", &
+         "                       takes a comma separated list of temperatures", &
+         "", &
+         "   --dftbplus             Read a DFTB+ hessian.out file, implies projection", &
+         "", &
+         "   --turbomole         Read a Turbomole Hessian file", &
+         "                       use this only when $nomw is not present in control", &
+         "", &
+         "   --orca         Read an Orca Hessian file", &
+         "", &
+         "   --bhess         Read in an <hessian_sph> file from a previous xTB SPH calculation ", &
+         "                       and scale frequencies accordingly", &
+         "", &
+         "Output Conventions:", &
+         "", &
+         "total energies are given in atomic units (Eh), frequencies in reciprocal cm", &
+         "", &
+         "More information can be obtained at https://xtb-docs.readthedocs.io/", &
+         ""
+   end subroutine thermoHelp
 
 end module xtb_prog_thermo