diff --git a/src/hessian.F90 b/src/hessian.F90
index 0e274c801..78b306e5b 100644
--- a/src/hessian.F90
+++ b/src/hessian.F90
@@ -504,6 +504,12 @@ subroutine numhess( &
       write(env%unit, '(A)') "DFTB+ style hessian.out written"
    end if
 
+   ! If we are currently computing a Hessian for Gaussian to use, we return it
+   ! right here and now and skip all the post-processing.
+   if (geometry_inputfile .eq. p_geo_gaussian) then
+      return
+   end if
+
    ! prepare all for diag
    ! copy
    k=0
diff --git a/src/io/reader/gaussian.f90 b/src/io/reader/gaussian.f90
index 8ad3ede05..3cbeb7489 100644
--- a/src/io/reader/gaussian.f90
+++ b/src/io/reader/gaussian.f90
@@ -29,6 +29,7 @@ module xtb_io_reader_gaussian
 
 subroutine readMoleculeGaussianExternal(mol, unit, status, iomsg)
    use xtb_mctc_accuracy, only : wp
+   use xtb_setmod, only : set_geopref
    logical, parameter :: debug = .false.
    class(TMolecule), intent(out) :: mol
    integer, intent(in) :: unit
@@ -44,6 +45,8 @@ subroutine readMoleculeGaussianExternal(mol, unit, status, iomsg)
 
    status = .false.
 
+   call set_geopref('gaussian')
+
    read(unit, '(4i10)', iostat=err) n, mode, chrg, spin
    if (err.ne.0) then
       iomsg = "Could not read number of atoms, check format!"
diff --git a/src/prog/main.F90 b/src/prog/main.F90
index 0238828d9..aeb9327c4 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -265,7 +265,7 @@ subroutine xtbMain(env, argParser)
          call close_file(ich)
       end if
    endif
-
+   
    !> efield read: gfnff only
    call open_file(ich,'.EFIELD','r')
    if (ich.ne.-1) then
@@ -321,7 +321,7 @@ subroutine xtbMain(env, argParser)
       fname = 'caffeine'
       call get_coffee(mol)
       call generateFileMetaInfo(fname, directory, basename, extension)
-   else
+   else                                                              
       call generateFileMetaInfo(fname, directory, basename, extension)
       ftype = getFileType(basename, extension)
       call open_file(ich, fname, 'r')
@@ -584,7 +584,7 @@ subroutine xtbMain(env, argParser)
    end select
    call env%checkpoint("Single point calculation terminated")
 
-   !> write 2d => 3d converted structure
+   !> write 2d => 3d converted structure   
    if (struc_conversion_done) then
       call generateFileName(tmpname, 'gfnff_convert', extension, mol%ftype)
       write(env%unit,'(10x,a,1x,a,/)') &
@@ -593,7 +593,7 @@ subroutine xtbMain(env, argParser)
       call writeMolecule(mol, ich, energy=res%e_total, gnorm=res%gnorm)
       call close_file(ich)
    end if
-
+   
    ! ========================================================================
    !> determine kopt for bhess including final biased geometry optimization
    if (runtyp.eq.p_run_bhess) then
@@ -857,7 +857,8 @@ subroutine xtbMain(env, argParser)
       call main_geometry(iprop,mol)
    endif
 
-   if ((runtyp.eq.p_run_hess).or.(runtyp.eq.p_run_ohess).or.(runtyp.eq.p_run_bhess)) then
+   if (((runtyp.eq.p_run_hess).or.(runtyp.eq.p_run_ohess).or.(runtyp.eq.p_run_bhess)) &
+      .and.(geometry_inputfile.ne.p_geo_gaussian)) then
       call generic_header(iprop,'Frequency Printout',49,10)
       call main_freq(iprop,mol,chk%wfn,fres)
    endif
@@ -888,7 +889,7 @@ subroutine xtbMain(env, argParser)
    type is(TGFFCalculator)
       call write_energy_gff(env%unit,res,fres, &
         & (runtyp.eq.p_run_hess).or.(runtyp.eq.p_run_ohess).or.(runtyp.eq.p_run_bhess))
-   end select
+   end select  
 
 
    ! ------------------------------------------------------------------------
@@ -1278,10 +1279,10 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
          call set_exttyp('eht')
          call env%warning("The use of '"//flag//"' is discouraged, " //&
             & "please use '--gfn 0' next time", source)
-
+      
       case('--gfnff')
          call set_exttyp('ff')
-
+      
       case('--gff')
          call set_exttyp('ff')
 
@@ -1315,7 +1316,7 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
 
       case('--json')
          call set_write(env,'json','true')
-
+       
       case('--ceasefiles')
          restart = .false.
          verbose=.false.
@@ -1326,7 +1327,7 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
          call set_opt(env, 'logfile', 'NUL')
 #else
          call set_opt(env, 'logfile', '/dev/null')
-#endif
+#endif         
 
       case('--orca')
          call set_exttyp('orca')
@@ -1461,7 +1462,7 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
          if (allocated(sec)) then
             call set_opt(env,'optlevel',sec)
          endif
-
+      
       case('--bhess')
          call set_runtyp('bhess')
          call args%nextArg(sec)
diff --git a/src/set_module.f90 b/src/set_module.f90
index 357847e3c..a688758cf 100644
--- a/src/set_module.f90
+++ b/src/set_module.f90
@@ -947,6 +947,9 @@ subroutine set_geopref(typ)
    case('vasp','poscar')
       geometry_inputfile = p_geo_poscar
 
+   case('EIn', 'gaussian')
+      geometry_inputfile = p_geo_gaussian
+
    end select
    set = .false.
 end subroutine set_geopref
diff --git a/src/setparam.f90 b/src/setparam.f90
index a0cbe404b..286bde66b 100644
--- a/src/setparam.f90
+++ b/src/setparam.f90
@@ -59,10 +59,11 @@ module xtb_setparam
    logical  :: periodic = .false.
 
 !  Geometry input type
-   integer,parameter :: p_geo_coord  = 1
-   integer,parameter :: p_geo_xmol   = 2
-   integer,parameter :: p_geo_sdf    = 3
-   integer,parameter :: p_geo_poscar = 4
+   integer,parameter :: p_geo_coord     = 1
+   integer,parameter :: p_geo_xmol      = 2
+   integer,parameter :: p_geo_sdf       = 3
+   integer,parameter :: p_geo_poscar    = 4
+   integer,parameter :: p_geo_gaussian  = 5
    integer  :: geometry_inputfile = p_geo_coord
 
 !! ------------------------------------------------------------------------