From 482072bc89098d9bd8a34f1ec41078d41cb1ed5a Mon Sep 17 00:00:00 2001 From: Cyrille Lavigne Date: Sun, 4 Oct 2020 08:42:30 -0400 Subject: [PATCH] Properly outputs Hessians in the format Gaussian wants --- src/io/writer/gaussian.f90 | 59 +++++++++++++++++++++++++++++++++++--- src/prog/main.F90 | 29 +++++++++++-------- 2 files changed, 72 insertions(+), 16 deletions(-) diff --git a/src/io/writer/gaussian.f90 b/src/io/writer/gaussian.f90 index 544fa08c4..18363687d 100644 --- a/src/io/writer/gaussian.f90 +++ b/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 diff --git a/src/prog/main.F90 b/src/prog/main.F90 index e2b01ba63..0238828d9 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 @@ -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)