Implementation of O1NumHess algorithm

README.md

@@ -166,6 +166,7 @@ features reality:
 - P. Pracht ([@pprcht](https://github.com/pprcht))
 - F. Pultar ([@pultar](https://github.com/pultar))
 - J. Seibert ([@liljay42](https://github.com/liljay42))
+- L. M. Seidler ([@lmseidler](https://github.com/lmseidler))
 - P. Shushkov
 - S. Spicher ([@sespic](https://github.com/sespic))
 - M. Stahn ([@MtoLStoN](https://github.com/mtolston))
@@ -245,6 +246,10 @@ for ONIOM refer to:
 - C. Plett, A. Katbashev, S. Ehlert, S. Grimme, M. Bursch, *Phys. Chem. Chem. Phys.*, **2023**, 25, 17860-17868.
   DOI: [10.1039/D3CP02178E](https://doi.org/10.1039/D3CP02178E)
 
+for O1NumHess refer to:
+- B. Wang, S. Luo, Z. Wang, W. Liu, *J. Chem. Theory Comput.* **2025**, 21, 10893-10909.
+  DOI: 10.1021/acs.jctc.5c01354
+
 All references are available in [bibtex format](./assets/references.bib).
 
 

assets/references.bib

@@ -203,3 +203,14 @@ @article{grimme_non-self-consistent_2023
 	pages = {124111},
 	file = {Full Text PDF:/Users/marcelmueller/Zotero/storage/LRJJ6T9I/Grimme et al. - 2023 - A non-self-consistent tight-binding electronic str.pdf:application/pdf;Snapshot:/Users/marcelmueller/Zotero/storage/FAT76IUH/A-non-self-consistent-tight-binding-electronic.html:text/html},
 }
+
+@article{wang2025o1numhess,
+  title={O1NumHess: a fast and accurate seminumerical Hessian algorithm using only O (1) gradients},
+  author={Wang, Bo and Luo, Shaohang and Wang, Zikuan and Liu, Wenjian},
+  journal={J{.} Chem{.} Theory Comput{.}},
+  volume={21},
+  number={21},
+  pages={10893--10909},
+  year={2025},
+  publisher={ACS Publications}
+}

man/xtb.1.adoc

@@ -234,6 +234,9 @@ NOTE: You can only select *one* runtyp, only the first runtyp will be used
     perform a biased numerical hessian calculation on an `ancopt(3)` optimized
     geometry
 
+*--o1nh*::
+    perform the numerical hessian calculation using the ODLR approximation (O1NumHess)
+
 *--md*::
     molecular dynamics simulation on start geometry
 

src/CMakeLists.txt

@@ -110,6 +110,7 @@ list(APPEND srcs
   "${dir}/modef.f90"
   "${dir}/model_hessian.f90"
   "${dir}/neighbor.f90"
+  "${dir}/o1numhess.f90"
   "${dir}/onetri.f90"
   "${dir}/oniom.f90"
   "${dir}/optimizer.f90"

src/hessian.F90

@@ -113,6 +113,9 @@ subroutine numhess( &
    real(wp),allocatable :: amass_au(:), amass_amu(:)
    real(wp) :: asq, gamsq
 
+   ! final residual error in approximating gradients using odlr hessian
+   real(wp) :: final_err
+
    type(TMolecule) :: tmol
 
    logical :: ex,rd
@@ -139,6 +142,7 @@ subroutine numhess( &
    rd=.false.
    xyzsave = mol%xyz
 
+   ! TODO: need to figure out why this is done
    step=0.0001_wp
    call rotmol(mol%n,mol%xyz,step,2.*step,3.*step)
 
@@ -185,7 +189,13 @@ subroutine numhess( &
 !! ========================================================================
 !  Hessian part -----------------------------------------------------------
 
-   if(freezeset%n.gt.0) then
+   if (set%o1numhess) then
+      h = 0.0_wp
+      call calc%odlrhessian(env, mol, chk0, step, h, final_err)
+      call env%check(exitRun)
+      if (exitRun) return
+      write(env%unit, '(A,1X,ES12.4)') "Error norm for predicted gradient (ODLR Hessian):", final_err
+   else if(freezeset%n.gt.0) then
       ! for frozfc of about 10 the frozen modes
       ! approach 5000 cm-1, i.e., come too close to
       ! the real ones
@@ -212,7 +222,6 @@ subroutine numhess( &
       dipd = 0.0_wp
       pold = 0.0_wp
       call calc%hessian(env, mol, chk0, indx(:nonfrozh), step, h, dipd)
-
    else
 !! ------------------------------------------------------------------------
 !  normal case
@@ -461,7 +470,8 @@ subroutine numhess( &
          ! for non-linear systems unless one fixes three atoms defines plane, 1 degree of freedom will exist, otherwise there should be 0 degrees of freedom
          ! anyway, the check here will become more complex and therefore it is not impemented
          ! NOTE: it is not necessary lowest N frequencies
-         error stop "not implemented"
+         call env%error("check for <=2 frozen atoms not implemented", source)
+         return
          ! for three atom systems we assume that the plane was constructed (or linear system is used)
       endif
       j=kend
@@ -498,57 +508,60 @@ subroutine numhess( &
       res%rmass(i)= 1.0_wp / xsum
    enddo
 
-   !--- IR intensity ---! (holds in a similar fashion also for Raman)
-   !  1. res%hess corresponds to the orthonormal eigenvectors of the mass-weighted Hessian
-   !     matrix (-> normal modes of vibration). By mass-weighting the Hessian matrix,
-   !     the normal modes are transformed into the mass-weighted space ("Q basis"), and
-   !     have the units [sqrt(mass) * length]
-   !     To obtain purely cartesian coordinates (-> transforming back into the Cartesian space),
-   !     the mass-weighted normal modes have to be divided by the square root of the mass of the respective atom.
-   !
-   !  2. res%hess(j,i) is the matrix which transforms a derivative with
-   !     respect to the j-th cartesian coordinate ("dipd") into a derivative with
-   !     respect to the i-th normal coordinate.
-   !
-   !  3. amass_au(j) = 1/sqrt(m(j)); m(j) is given in atomic units (a.u.).
-   !
-   !  4. matmul(D x H) = U
-   !
-   !  5. D = dipd(3,n3); H = res%hess(n3:n3); U = Matrix with dipol derivatives
-   !                                              in x, y and z direction per mode
-   !
-   ! Generally nice reads for understanding the necessity of mass-weighting:
-   ! 1) https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/
-   !         Advanced_Theoretical_Chemistry_(Simons)/
-   !         03%3A_Characteristics_of_Energy_Surfaces/3.02%3A_Normal_Modes_of_Vibration
-   ! 2) https://www.cup.uni-muenchen.de/ch/compchem/G98vib.pdf
-
-   do i = 1, n3
-      do k = 1, 3
-         sum2 = 0.0_wp
-         do j = 1, n3
-            sum2 = sum2 + dipd(k,j)*(res%hess(j,i)*amass_au(j))
-         end do
-         trdip(k) = sum2
-      end do
-      res%dipt(i) = autokmmol*(trdip(1)**2+trdip(2)**2+trdip(3)**2)
-   end do
-   ! Raman activities (for intensities, see "write_tm_vibspectrum")
-   if (set%elprop == p_elprop_alpha) then
+   ! odlr hessian check, cannot calculate IR intensities
+   if (.not. set%o1numhess) then
+      !--- IR intensity ---! (holds in a similar fashion also for Raman)
+      !  1. res%hess corresponds to the orthonormal eigenvectors of the mass-weighted Hessian
+      !     matrix (-> normal modes of vibration). By mass-weighting the Hessian matrix,
+      !     the normal modes are transformed into the mass-weighted space ("Q basis"), and
+      !     have the units [sqrt(mass) * length]
+      !     To obtain purely cartesian coordinates (-> transforming back into the Cartesian space),
+      !     the mass-weighted normal modes have to be divided by the square root of the mass of the respective atom.
+      !
+      !  2. res%hess(j,i) is the matrix which transforms a derivative with
+      !     respect to the j-th cartesian coordinate ("dipd") into a derivative with
+      !     respect to the i-th normal coordinate.
+      !
+      !  3. amass_au(j) = 1/sqrt(m(j)); m(j) is given in atomic units (a.u.).
+      !
+      !  4. matmul(D x H) = U
+      !
+      !  5. D = dipd(3,n3); H = res%hess(n3:n3); U = Matrix with dipol derivatives
+      !                                              in x, y and z direction per mode
+      !
+      ! Generally nice reads for understanding the necessity of mass-weighting:
+      ! 1) https://chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/
+      !         Advanced_Theoretical_Chemistry_(Simons)/
+      !         03%3A_Characteristics_of_Energy_Surfaces/3.02%3A_Normal_Modes_of_Vibration
+      ! 2) https://www.cup.uni-muenchen.de/ch/compchem/G98vib.pdf
+
       do i = 1, n3
-         do k = 1,6
+         do k = 1, 3
             sum2 = 0.0_wp
             do j = 1, n3
-               sum2 = sum2 + (res%hess(j,i)*amass_au(j))*dalphadr(k,j)
+               sum2 = sum2 + dipd(k,j)*(res%hess(j,i)*amass_au(j))
+            end do
+            trdip(k) = sum2
+         end do
+         res%dipt(i) = autokmmol*(trdip(1)**2+trdip(2)**2+trdip(3)**2)
+      end do
+      ! Raman activities (for intensities, see "write_tm_vibspectrum")
+      if (set%elprop == p_elprop_alpha) then
+         do i = 1, n3
+            do k = 1,6
+               sum2 = 0.0_wp
+               do j = 1, n3
+                  sum2 = sum2 + (res%hess(j,i)*amass_au(j))*dalphadr(k,j)
+               enddo
+               dalphadq(k,i) = sum2
             enddo
-            dalphadq(k,i) = sum2
+            asq = (dalphadq(1,i)+dalphadq(3,i)+dalphadq(6,i))**2 / 9.0_wp
+            gamsq = ( (dalphadq(1,i)-dalphadq(3,i))**2 + (dalphadq(3,i)-dalphadq(6,i))**2 + (dalphadq(6,i)-dalphadq(1,i))**2 &
+               & + 6.0_wp*(dalphadq(2,i)**2 + dalphadq(5,i)**2 + dalphadq(4,i)**2) )*0.5_wp
+            res%polt(i) = (45.0_wp*asq + 7.0_wp*gamsq)
+            res%polt(i) = res%polt(i) * autoaa4byamu()
          enddo
-         asq = (dalphadq(1,i)+dalphadq(3,i)+dalphadq(6,i))**2 / 9.0_wp
-         gamsq = ( (dalphadq(1,i)-dalphadq(3,i))**2 + (dalphadq(3,i)-dalphadq(6,i))**2 + (dalphadq(6,i)-dalphadq(1,i))**2 &
-            & + 6.0_wp*(dalphadq(2,i)**2 + dalphadq(5,i)**2 + dalphadq(4,i)**2) )*0.5_wp
-         res%polt(i) = (45.0_wp*asq + 7.0_wp*gamsq)
-         res%polt(i) = res%polt(i) * autoaa4byamu()
-      enddo
+      end if
    end if
 
 end subroutine numhess

src/main/property.F90

@@ -652,25 +652,26 @@ subroutine main_freq &
 
       write (iunit, '(1x,a)') 'reduced masses (amu)'
       write (iunit, '(6(i5,'':'',f10.2))') (i, res%rmass(i), i=1, res%n3)
-      write (iunit, '(1x,a)') 'IR intensities (km·mol⁻¹)'
-      write (iunit, '(6(i5,'':'',f10.2))') (i, res%dipt(i), i=1, res%n3)
-      write (iunit, '(1x,a)') 'Raman intensities (Ä⁴*amu⁻¹)'
-      write (iunit, '(6(i5,'':'',f10.2))') (i, res%polt(i), i=1, res%n3)
-
-      call open_file(ifile, 'vibspectrum', 'w')
-      if (set%elprop == p_elprop_alpha) then
-         call write_tm_vibspectrum(ifile, res%n3, res%freq, res%dipt, res%polt,&
-                                           set%ptbsetup%raman_temp, set%ptbsetup%raman_lambda)
-      else
-         call write_tm_vibspectrum(ifile, res%n3, res%freq, res%dipt, res%polt)
+      if (.not. set%o1numhess) then
+         write (iunit, '(1x,a)') 'IR intensities (km·mol⁻¹)'
+         write (iunit, '(6(i5,'':'',f10.2))') (i, res%dipt(i), i=1, res%n3)
+         write (iunit, '(1x,a)') 'Raman intensities (Ä⁴*amu⁻¹)'
+         write (iunit, '(6(i5,'':'',f10.2))') (i, res%polt(i), i=1, res%n3)
+         call open_file(ifile, 'vibspectrum', 'w')
+         if (set%elprop == p_elprop_alpha) then
+            call write_tm_vibspectrum(ifile, res%n3, res%freq, res%dipt, res%polt,&
+            set%ptbsetup%raman_temp, set%ptbsetup%raman_lambda)
+         else
+            call write_tm_vibspectrum(ifile, res%n3, res%freq, res%dipt, res%polt)
+         end if
+         call close_file(ifile)
+         write (iunit, '(1x,a)') 'output can be read by thermo (or use thermo option).'
+         write (iunit, '(1x,a)') 'writing <g98.out> molden fake output.'
+         write (iunit, '(1x,a)') &
+         & 'recommended (thermochemical) frequency scaling factor: 1.0'
+         call g98fake2('g98.out', mol%n, mol%at, mol%xyz, res%freq, res%rmass, res%dipt, res%hess)
       end if
-      call close_file(ifile)
 
-      write (iunit, '(1x,a)') 'output can be read by thermo (or use thermo option).'
-      write (iunit, '(1x,a)') 'writing <g98.out> molden fake output.'
-      write (iunit, '(1x,a)') &
-         & 'recommended (thermochemical) frequency scaling factor: 1.0'
-      call g98fake2('g98.out', mol%n, mol%at, mol%xyz, res%freq, res%rmass, res%dipt, res%hess)
 
       if (set%pr_nmtm) then
          call open_file(ifile, "vib_normal_modes", 'w')

src/meson.build

@@ -108,6 +108,7 @@ srcs += files(
   'modef.f90',
   'model_hessian.f90',
   'neighbor.f90',
+  'o1numhess.f90',
   'onetri.f90',
   'oniom.f90',
   'optimizer.f90',