Fix imaginary frequencies

src/hessian.F90

@@ -86,11 +86,12 @@ subroutine numhess( &
    real(wp) :: sum1,sum2,trdip(3),dipole(3)
    real(wp) :: trpol(3),sl(3,3)
    integer  :: n3,i,j,k,ic,jc,ia,ja,ii,jj,info,lwork,a,b,ri,rj
-   integer  :: nread,kend,lowmode
+   integer  :: nread,lowmode,nskip
    integer  :: nonfrozh
    integer  :: fixmode
+   integer, allocatable :: skiplist(:)
    integer, allocatable :: nb(:,:)
-   integer, allocatable :: indx(:),molvec(:),izero(:)
+   integer, allocatable :: indx(:),molvec(:)
    real(wp),allocatable :: bond(:,:)
 
 !$ integer  :: nproc
@@ -129,7 +130,7 @@ subroutine numhess( &
    allocate(hss(n3*(n3+1)/2),hsb(n3*(n3+1)/2),h(n3,n3),htb(n3,n3),hbias(n3,n3), &
       & gl(3,mol%n),isqm(n3),xyzsave(3,mol%n),dipd(3,n3), &
       & pold(n3),nb(20,mol%n),indx(mol%n),molvec(mol%n),bond(mol%n,mol%n), &
-      & freq_scal(n3),fc_tb(n3),fc_bias(n3),amass(n3), h_dummy(n3,n3), izero(n3))
+      & freq_scal(n3),fc_tb(n3),fc_bias(n3),amass(n3), h_dummy(n3,n3))
 
    if (set%elprop == p_elprop_alpha) then
       allocate(dalphadr(6,n3), source = 0.0_wp)
@@ -412,14 +413,9 @@ subroutine numhess( &
    else
       write(env%unit,'(1x,a)') 'projected vibrational frequencies (cm⁻¹)'
    endif
-   k=0
    do i=1,n3
       ! Eigenvalues in atomic units, convert to wavenumbers
       res%freq(i)=autorcm*sign(sqrt(abs(res%freq(i))),res%freq(i))/sqrt(amutoau)
-      if(abs(res%freq(i)).lt.0.01_wp) then
-         k=k+1
-         izero(k)=i
-      endif
    enddo
 
    ! scale frequencies
@@ -442,42 +438,26 @@ subroutine numhess( &
    if (mol%n > 1) then
       h = 0.0_wp
       isqm = 0.0_wp
-      kend=0
-      if (freezeset%n == 0) then
-         kend=6
-         if(res%linear)then
-            kend=5
-            do i=1,kend
-               izero(i)=i
-            enddo
-            res%freq(1:kend)=0
-         endif
-         do k=1,kend
-            h(1:n3,k)=res%hess(1:n3,izero(k))
-            isqm(  k)=res%freq(izero(k))
-         enddo
-      else if (freezeset%n <= 2) then
-         ! for systems with one fixed atom, there should be 2 and 3 degrees of freedom for linear and non-linear systems, respectively
-         ! for systems with two fixed atoms, there should be 0 and 1 degrees of freedom for linear and non-linear systems, respectively
-         ! for linear systems with more than two fixed atoms, there should be 0 degrees of freedom
-         ! 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"
-         ! for three atom systems we assume that the plane was constructed (or linear system is used)
-      endif
-      j=kend
-      do k=1,n3
-         if(abs(res%freq(k)).gt.0.01_wp)then
-            j=j+1
-            if(j.gt.n3) then
-               call env%error('internal error while sorting hessian', source)
-               return
-            end if
-            h(1:n3,j)=res%hess(1:n3,k)
-            isqm(  j)=res%freq(   k)
+      j = 0
+      nskip = 0
+      allocate(skiplist(n3), source = 0)
+      do k=1, n3
+         if (abs(res%freq(k)) > 0.05_wp) then
+            j = j + 1
+            h(1:n3,j) = res%hess(1:n3,k)
+            isqm(  j) = res%freq(   k)
+         else
+            nskip = nskip + 1
+            skiplist(nskip) = k
          endif
       enddo
+      do while (j /= n3)
+         j = j + 1
+         h(1:n3,j) = res%hess(1:n3,skiplist(nskip))
+         isqm(  j) = 0.0_wp
+         nskip = nskip - 1
+      end do
+      deallocate(skiplist)
    end if
    res%hess = h
    res%freq = isqm

src/main/property.F90

@@ -1199,7 +1199,7 @@ subroutine print_thermo(iunit, nat, nvib_in, at, xyz, freq, etot, htot, gtot, ni
       real(wp) xx(10), sthr, temp, scale_factor
       real(wp) aa, bb, cc, vibthr, ithr
       real(wp) escf, symnum, wt, avmom, diff
-      real(wp) :: omega, maxfreq, fswitch, lnq_r, lnq_v
+      real(wp) :: omega, fswitch, lnq_r, lnq_v
       real(wp), allocatable :: et(:), ht(:), gt(:), ts(:)
       integer nn, nvib, i, j, k, n, nvib_theo, isthr
       integer, intent(out) :: nimag
@@ -1373,10 +1373,9 @@ subroutine print_thermo_sthr_lnq(iunit, nvib, vibs, avmom, sthr, temp)
       real(wp), intent(in) :: temp
 
       integer :: i
-      real(wp) :: maxfreq, omega, lnq_r, lnq_v, fswitch
+      real(wp) :: omega, lnq_r, lnq_v, fswitch
 
       write (iunit, '(a)')
-      maxfreq = max(300.0_wp, chg_inverted(0.99_wp, sthr))
       write (iunit, '(a8,a14,a12,10x,a12,10x,a12)') &
          "mode", "ω/cm⁻¹", "ln{qvib}", "ln{qrot}", "ln{qtot}"
       write (iunit, '(3x,72("-"))')
@@ -1385,7 +1384,6 @@ subroutine print_thermo_sthr_lnq(iunit, nvib, vibs, avmom, sthr, temp)
          lnq_r = lnqvib(temp, omega)
          lnq_v = lnqrot(temp, omega, avmom)
          fswitch = 1.0_wp - chg_switching(omega, sthr)
-         if (omega > maxfreq) exit
          write (iunit, '(i8,f10.2,2(f12.5,1x,"(",f6.2,"%)"),f12.5)') &
             i, omega, lnq_v, (1.0_wp - fswitch) * 100, &
             lnq_r, fswitch * 100, (1.0_wp - fswitch) * lnq_v + fswitch * lnq_r
@@ -1408,15 +1406,14 @@ subroutine print_thermo_sthr_ts(iunit, nvib, vibs, avmom_si, sthr_rcm, temp)
       real(wp), intent(in) :: temp       !< temperature
 
       integer :: i
-      real(wp) :: maxfreq, omega, s_r, s_v, fswitch
+      real(wp) :: omega, s_r, s_v, fswitch
       real(wp) :: beta, xxmom, e, ewj, mu, RT, sthr, avmom
       beta = 1.0_wp / kB / temp ! beta in 1/Eh
       sthr = sthr_rcm * rcmtoau ! sthr in Eh
       RT = kb * temp * autokcal ! RT in kcal/mol for printout
       avmom = avmom_si * kgtome * aatoau**2 * 1.0e+20_wp ! in me·α²
 
       write (iunit, '(a)')
-      maxfreq = max(300.0_wp, chg_inverted(0.99_wp, sthr_rcm))
       write (iunit, '(a8,a14,1x,a27,a27,a12)') &
          "mode", "ω/cm⁻¹", "T·S(HO)/kcal·mol⁻¹", "T·S(FR)/kcal·mol⁻¹", "T·S(vib)"
       write (iunit, '(3x,72("-"))')
@@ -1446,7 +1443,6 @@ subroutine print_thermo_sthr_ts(iunit, nvib, vibs, avmom_si, sthr_rcm, temp)
          end if
          ! Head-Gordon weighting
          fswitch = 1.0_wp - chg_switching(omega, sthr)
-         if (omega > maxfreq * rcmtoau) exit
          write (iunit, '(i8,f10.2,2(f12.5,1x,"(",f6.2,"%)"),f12.5)') &
             i, omega * autorcm, -RT * s_v, (1.0_wp - fswitch) * 100, &
             -RT * s_r, fswitch * 100, -RT * ((1.0_wp - fswitch) * s_v + fswitch * s_r)