From f05c511fdccb746ea2cc8605633744c6c49be7c6 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:13:34 +0100 Subject: [PATCH 1/7] Do not drop high freqs Signed-off-by: Igor S. Gerasimov --- src/main/property.F90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/main/property.F90 b/src/main/property.F90 index ab22c3717..46b3c970c 100644 --- a/src/main/property.F90 +++ b/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,7 +1406,7 @@ 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 @@ -1416,7 +1414,6 @@ subroutine print_thermo_sthr_ts(iunit, nvib, vibs, avmom_si, sthr_rcm, temp) 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) From e888a059f353c725d2d2dd4495a8673f37d8ede1 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:15:32 +0100 Subject: [PATCH 2/7] Small reformat Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index 98b4ee547..f479aed80 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -467,15 +467,15 @@ subroutine numhess( & ! 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 + do k=1, n3 + if (abs(res%freq(k)) > 0.01_wp) then + j = j + 1 + if(j > 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) + h(1:n3,j) = res%hess(1:n3,k) + isqm( j) = res%freq( k) endif enddo end if From ea74337656530d7ef4564fa86606b1948b87d4aa Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:15:51 +0100 Subject: [PATCH 3/7] Update low freq threshold Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index f479aed80..5565cafb3 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -468,7 +468,7 @@ subroutine numhess( & endif j=kend do k=1, n3 - if (abs(res%freq(k)) > 0.01_wp) then + if (abs(res%freq(k)) > 0.05_wp) then j = j + 1 if(j > n3) then call env%error('internal error while sorting hessian', source) From 0f6f80c0aec504812a48d24cf42255a6e71efd33 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:18:13 +0100 Subject: [PATCH 4/7] Sort all freqs Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 28 ++-------------------------- 1 file changed, 2 insertions(+), 26 deletions(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index 5565cafb3..2ce54be0c 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -86,7 +86,7 @@ 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 integer :: nonfrozh integer :: fixmode integer, allocatable :: nb(:,:) @@ -442,31 +442,7 @@ 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 + j = 0 do k=1, n3 if (abs(res%freq(k)) > 0.05_wp) then j = j + 1 From 7168387ee62a83c2a0e83bf66792917313c1504e Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:21:52 +0100 Subject: [PATCH 5/7] Sort up to 6 zero freqs Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index 2ce54be0c..1e4a09a89 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -89,6 +89,7 @@ subroutine numhess( & integer :: nread,lowmode integer :: nonfrozh integer :: fixmode + integer :: skiplist(6), nskip integer, allocatable :: nb(:,:) integer, allocatable :: indx(:),molvec(:),izero(:) real(wp),allocatable :: bond(:,:) @@ -443,17 +444,27 @@ subroutine numhess( & h = 0.0_wp isqm = 0.0_wp j = 0 + nskip = 0 do k=1, n3 if (abs(res%freq(k)) > 0.05_wp) then j = j + 1 - if(j > n3) then + h(1:n3,j) = res%hess(1:n3,k) + isqm( j) = res%freq( k) + else + nskip = nskip + 1 + if(nskip > 6) 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) + 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 end if res%hess = h res%freq = isqm From 821d7fb40cf3bcd7e02983c0022fd97d4e76e1f3 Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:24:20 +0100 Subject: [PATCH 6/7] Allow n3 zero freqs (for frozen atoms) Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index 1e4a09a89..fe01bacbb 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -86,10 +86,10 @@ 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,lowmode + integer :: nread,lowmode,nskip integer :: nonfrozh integer :: fixmode - integer :: skiplist(6), nskip + integer, allocatable :: skiplist(:) integer, allocatable :: nb(:,:) integer, allocatable :: indx(:),molvec(:),izero(:) real(wp),allocatable :: bond(:,:) @@ -445,6 +445,7 @@ subroutine numhess( & isqm = 0.0_wp 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 @@ -452,10 +453,6 @@ subroutine numhess( & isqm( j) = res%freq( k) else nskip = nskip + 1 - if(nskip > 6) then - call env%error('internal error while sorting hessian', source) - return - end if skiplist(nskip) = k endif enddo @@ -465,6 +462,7 @@ subroutine numhess( & isqm( j) = 0.0_wp nskip = nskip - 1 end do + deallocate(skiplist) end if res%hess = h res%freq = isqm From 0a1c9e191d5f035150a2ab2df7b2632e3d872a7a Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" Date: Wed, 5 Mar 2025 05:51:49 +0100 Subject: [PATCH 7/7] Remove izero Signed-off-by: Igor S. Gerasimov --- src/hessian.F90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/hessian.F90 b/src/hessian.F90 index fe01bacbb..d36f0942f 100644 --- a/src/hessian.F90 +++ b/src/hessian.F90 @@ -91,7 +91,7 @@ subroutine numhess( & integer :: fixmode integer, allocatable :: skiplist(:) integer, allocatable :: nb(:,:) - integer, allocatable :: indx(:),molvec(:),izero(:) + integer, allocatable :: indx(:),molvec(:) real(wp),allocatable :: bond(:,:) !$ integer :: nproc @@ -130,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) @@ -413,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