From 76faf341a9f43b8cf06cef8fe209a33c0e54cfc5 Mon Sep 17 00:00:00 2001 From: BelizSertcan <48097107+BelizSertcan@users.noreply.github.com> Date: Tue, 9 Sep 2025 23:55:42 +0200 Subject: [PATCH] bugfix in force constant calculation --- src/common/physcon.F | 2 +- src/motion/vibrational_analysis.F | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/common/physcon.F b/src/common/physcon.F index 033ffc146b..6b7fba46df 100644 --- a/src/common/physcon.F +++ b/src/common/physcon.F @@ -77,7 +77,7 @@ MODULE physcon PUBLIC :: a_bohr, a_fine, a_mass, angstrom, atm, bar, bohr, boltzmann, c_light_au, & debye, e_charge, e_gfactor, e_mass, evolt, femtoseconds, h_bar, & hertz, joule, kcalmol, kelvin, kjmol, massunit, mu_perm, n_avogadro, newton, & - p_mass, pascal, picoseconds, seconds, vibfac, wavenumbers + p_mass, pascal, picoseconds, seconds, vibfac, wavenumbers, c_light PUBLIC :: write_physcon diff --git a/src/motion/vibrational_analysis.F b/src/motion/vibrational_analysis.F index 0f73d74143..f699ec03d5 100644 --- a/src/motion/vibrational_analysis.F +++ b/src/motion/vibrational_analysis.F @@ -72,7 +72,7 @@ MODULE vibrational_analysis USE particle_types, ONLY: particle_type USE physcon, ONLY: & a_bohr, angstrom, bohr, boltzmann, debye, e_mass, h_bar, hertz, kelvin, kjmol, massunit, & - n_avogadro, pascal, vibfac, wavenumbers + n_avogadro, pascal, vibfac, wavenumbers, c_light, joule USE replica_methods, ONLY: rep_env_calc_e_f,& rep_env_create USE replica_types, ONLY: rep_env_release,& @@ -632,8 +632,6 @@ CONTAINS rmass(i) = norm/massunit ! Renormalize displacements and convert in Angstrom D(:, i) = SQRT(norm)*D(:, i) - ! Force constants - konst(i) = SIGN(1.0_dp, H_eigval2(i))*2.0_dp*pi**2*(ABS(H_eigval2(i))/massunit)**2*rmass(i) IF (calc_intens) THEN D_deriv = 0._dp DO j = 1, nvib @@ -671,6 +669,8 @@ CONTAINS END IF ! Convert frequencies to cm^-1 H_eigval2(i) = SIGN(1.0_dp, H_eigval2(i))*SQRT(ABS(H_eigval2(i))*massunit)*vibfac/1000.0_dp + ! Force constant in au, conversion to mdyne/A is 15.57 + konst(i) = SIGN(1.0_dp, H_eigval2(i))*rmass(i)*massunit*(2.0_dp*pi*c_light*100*ABS(H_eigval2(i))*h_bar/joule)**2 END DO IF (calc_intens) THEN IF (iounit > 0) THEN @@ -881,7 +881,7 @@ CONTAINS WRITE (UNIT=iw, FMT="(T2,'VIB|Red.Masses (a.u.)',3(1X,F12.6,8X))") & (m(icol), icol=from, to) WRITE (UNIT=iw, FMT="(T2,'VIB|Frc consts (a.u.)',3(1X,F12.6,8X))") & - (k(icol), icol=from, to) + (k(icol), icol=from, to) WRITE (UNIT=iw, FMT="(T2,' ATOM',2X,'EL',7X,3(4X,' X ',1X,' Y ',1X,' Z '))") DO iatom = 1, natom, 3 katom = iatom/3