Adds output of Hessians to xtb when acting as an external program for Gaussian 16

scripts/xtb-gaussian

@@ -1,2 +1,3 @@
 #!/usr/bin/env bash
-2>&1 exec xtb $2 --grad > $4
+flag=$(awk 'NR == 1 {if ($2 == 2) {print "--hess"} else {print "--grad"} }' $2)
+2>&1 exec xtb $2 $flag > $4

src/hessian.F90

@@ -302,6 +302,15 @@
       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
+      ! Store the raw dipole gradient in the intensity output,
+      ! this we require a LHS reallocation of dipt from 3*n to 9*n
+      res%dipt = reshape(dipd, [size(dipd)])
+      return
+   end if
+
    ! prepare all for diag
    ! copy
    k=0

src/io/reader/gaussian.f90

@@ -29,6 +29,7 @@
 
 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 @@
 
    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!"

src/io/writer/gaussian.f90

@@ -41,14 +41,65 @@
 end subroutine writeMoleculeGaussianExternal
 
 
-subroutine writeResultsGaussianExternal(unit, energy, dipole, gradient)
-   integer, intent(in) :: unit
+subroutine writeResultsGaussianExternal(unit, energy, dipole, gradient, hessian, dipgrad)
+   integer, intent(in)  :: unit
    real(wp), intent(in) :: energy
    real(wp), intent(in) :: dipole(:)
-   real(wp), intent(in) :: gradient(:, :)
+   real(wp), intent(in), optional :: gradient(:, :)
+   real(wp), intent(in), optional :: hessian(:, :)
+   real(wp), intent(in), optional :: dipgrad(:)
+
+   integer :: i,j,ij,nat
+
+   if (present(gradient)) then
+      nat = size(gradient, 2)
+   else if (present(hessian)) then
+      nat = size(hessian, 2) / 3
+   else if (present(dipgrad)) then
+      nat = size(dipgrad) / 9
+   else
+      nat = 0
+   end if
 
    write(unit, '(4D20.12)') energy, dipole
-   write(unit, '(3D20.12)') gradient
+
+   if (present(gradient)) then
+      write(unit, '(3D20.12)') gradient
+   else
+      ! In case we have no gradient, but either hessian or dipole gradient
+      ! the gradient entry is padded with zeros
+      if (nat > 0) write(unit, '(3D20.12)') spread(0.0_wp, 1, 3*nat)
+   endif
+
+   ! First, we fake the polarizability and the dipole derivatives, which the
+   ! Gaussian Manual says should be of this form,
+   !
+   ! Polar(I), I=1,6        3D20.12   2 rows of 3 x 0.0
+   write(unit, '(3D20.12)') spread(0.0_wp, 1, 6)
+
+   ! DDip(I), I=1,9*NAtoms  3D20.12   9 x Natoms or 3 natoms rows of 3 x .0.0
+   if (present(dipgrad)) then
+      write(unit, '(3D20.12)') dipgrad
+   else
+      if (nat > 0) write(unit, '(3D20.12)') spread(0.0_wp, 1, 9*nat)
+   end if
+
+   if (present(hessian)) then
+      ! Now we are iterating over Hessian matrix elements. We will
+      ! append those to the output file in the correct format, given in
+      ! the Gaussian manual as
+      !
+      ! FFX(I), I=1,(3*NAtoms*(3*NAtoms+1))/2      3D20.12
+
+      ij = 0
+      do i = 1, size(hessian, 2)
+         do j = 1, i
+            ij = ij + 1
+            write(unit, '(D20.12)', advance='no') hessian(j, i)
+            if (mod(ij, 3) == 0) write(unit, '(a)')
+         end do
+      end do
+   end if
 
 end subroutine writeResultsGaussianExternal
 

src/prog/main.F90

@@ -821,8 +821,14 @@
       else
          cdum = 'xtb-gaussian.EOu'
       end if
+
       call open_file(ich, cdum, 'w')
-      call writeResultsGaussianExternal(ich, etot, res%dipole, g)
+      if (allocated(fres%hess)) then
+         call writeResultsGaussianExternal(ich, etot, res%dipole, gradient=g, &
+            & hessian=fres%hess, dipgrad=fres%dipt)
+      else
+         call writeResultsGaussianExternal(ich, etot, res%dipole, gradient=g)
+      end if
       call close_file(ich)
    end if
 
@@ -854,7 +860,8 @@
       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
@@ -1318,7 +1325,7 @@
          call set_write(env,'json','true')
 
       case('--ceasefiles')
-         restart = .false. 
+         restart = .false.
          verbose=.false.
          ceasefiles = .true.
          call set_write(env,'wiberg','false')

src/set_module.f90

@@ -946,6 +946,9 @@
 
    case('vasp','poscar')
       geometry_inputfile = p_geo_poscar
+
+   case('EIn', 'gaussian')
+      geometry_inputfile = p_geo_gaussian
    end select
    set = .false.
 end subroutine set_geopref

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
 
 !! ------------------------------------------------------------------------