Properly outputs Hessians in the format Gaussian wants

src/io/writer/gaussian.f90

@@ -41,14 +41,65 @@ subroutine writeMoleculeGaussianExternal(mol, unit)
 end subroutine writeMoleculeGaussianExternal
 
 
-subroutine writeResultsGaussianExternal(unit, energy, dipole, gradient)
-   integer, intent(in) :: unit
+subroutine writeResultsGaussianExternal(mol, unit, energy, dipole, gradient, hess)
+   implicit none
+   type(TMolecule), intent(in) :: mol
+   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 :: hess(:, :)
+
+   integer :: i,j,count
+   real(wp) :: buf(3)
 
    write(unit, '(4D20.12)') energy, dipole
-   write(unit, '(3D20.12)') gradient
+
+   if (present(gradient)) then
+      write(unit, '(3D20.12)') gradient
+   elseif (present(hess)) then
+      ! In case the gradient is not computed but the Hessian is, we output
+      ! gradients of zero. I don't know what would cause that to be honest...
+      write(unit, '(3D20.12)') (/ 0.0, 0.0, 0.0 /)
+   endif
+
+   if (present(hess)) then
+      ! 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
+      do i=1,6
+         write(unit, '(3D20.12)') (/ 0.0, 0.0, 0.0 /)
+      enddo
+
+      ! DDip(I), I=1,9*NAtoms    3D20.12
+      do i=1,9 * mol%n
+         write(unit, '(3D20.12)') (/ 0.0, 0.0, 0.0 /)
+      enddo
+
+      ! 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
+
+      ! that is, we need to print in group of three numbers, which is why the
+      ! whole buffer / count thing is present
+      count = 1
+      do i=1,3*mol%n
+         do j=1,i
+            buf(count) = hess(i,j)
+
+            if (count .eq. 3) then
+               count = 1
+               write(unit, '(3D20.12)') buf
+            else
+               count = count + 1
+            end if
+
+         enddo
+      enddo
+   endif
 
 end subroutine writeResultsGaussianExternal
 

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
@@ -819,8 +819,13 @@ subroutine xtbMain(env, argParser)
       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(mol, ich, etot, res%dipole, gradient=g, hess=fres%hess)
+      else
+         call writeResultsGaussianExternal(mol, ich, etot, res%dipole, gradient=g)
+      end if
       call close_file(ich)
    end if
 
@@ -883,7 +888,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
 
 
    ! ------------------------------------------------------------------------
@@ -1273,10 +1278,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')
 
@@ -1310,9 +1315,9 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
 
       case('--json')
          call set_write(env,'json','true')
-       
+
       case('--ceasefiles')
-         restart = .false. 
+         restart = .false.
          verbose=.false.
          ceasefiles = .true.
          call set_write(env,'wiberg','false')
@@ -1321,7 +1326,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')
@@ -1456,7 +1461,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)