From a07a8b12c87673860a378b75a03337cebec0e4f1 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 18 Oct 2023 17:45:20 -0400 Subject: [PATCH 1/3] Reading Gaussian fchk vibrations when present Signed-off-by: Geoff Hutchison --- avogadro/quantumio/gaussianfchk.cpp | 55 ++++++++++++++++++++++++++++- avogadro/quantumio/gaussianfchk.h | 14 ++++++-- 2 files changed, 66 insertions(+), 3 deletions(-) diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp index c58f0ce3f3..9f748a1269 100644 --- a/avogadro/quantumio/gaussianfchk.cpp +++ b/avogadro/quantumio/gaussianfchk.cpp @@ -57,6 +57,17 @@ bool GaussianFchk::read(std::istream& in, Core::Molecule& molecule) m_aPos[i + 1] * BOHR_TO_ANGSTROM, m_aPos[i + 2] * BOHR_TO_ANGSTROM)); } + + if (m_frequencies.size() > 0 && + m_frequencies.size() == m_vibDisplacements.size() && + m_frequencies.size() == m_IRintensities.size()) { + molecule.setVibrationFrequencies(m_frequencies); + molecule.setVibrationIRIntensities(m_IRintensities); + molecule.setVibrationLx(m_vibDisplacements); + if (m_RamanIntensities.size()) + molecule.setVibrationRamanIntensities(m_RamanIntensities); + } + // Do simple bond perception. molecule.perceiveBondsSimple(); molecule.perceiveBondOrders(); @@ -91,7 +102,7 @@ void GaussianFchk::processLine(std::istream& in) } else if (Core::contains(key, "UHF")) { m_scftype = Uhf; } else if (key == "Number of atoms" && list.size() > 1) { - cout << "Number of atoms = " << Core::lexicalCast(list[1]) << endl; + m_numAtoms = Core::lexicalCast(list[1]); } else if (key == "Charge" && list.size() > 1) { m_charge = Core::lexicalCast(list[1]); } else if (key == "Multiplicity" && list.size() > 1) { @@ -183,6 +194,48 @@ void GaussianFchk::processLine(std::istream& in) << endl; else cout << "Error reading in the SCF spin density matrix.\n"; + } else if (key == "Number of Normal Modes" && list.size() > 1) { + m_normalModes = Core::lexicalCast(list[1]); + } else if (key == "Vib-E2" && list.size() > 2) { + m_frequencies.clear(); + m_IRintensities.clear(); + m_RamanIntensities.clear(); + + unsigned threeN = m_numAtoms * 3; // degrees of freedom + std::vector tmp = + readArrayD(in, Core::lexicalCast(list[2]), 16); + + // read in the first 3N-6 elements as frequencies + for (unsigned int i = 0; i < m_normalModes; ++i) { + m_frequencies.push_back(tmp[i]); + } + // skip to after threeN elements then read IR intensities + for (unsigned int i = threeN; i < threeN + m_normalModes; ++i) { + m_IRintensities.push_back(tmp[i]); + } + // now check if we have Raman intensities + if (tmp[threeN + m_normalModes] != 0.0) { + for (unsigned int i = threeN + m_normalModes; + i < threeN + 2 * m_normalModes; ++i) { + m_RamanIntensities.push_back(tmp[i]); + } + } + } else if (key == "Vib-Modes" && list.size() > 2) { + std::vector tmp = readArrayD(in, Core::lexicalCast(list[2]), 16); + m_vibDisplacements.clear(); + if (tmp.size() == m_numAtoms * 3 * m_normalModes) { + for (unsigned int i = 0; i < m_normalModes; ++i) { + Core::Array mode; + for (unsigned int j = 0; j < m_numAtoms; ++j) { + Vector3 v(tmp[i * m_numAtoms * 3 + j * 3], + tmp[i * m_numAtoms * 3 + j * 3 + 1], + tmp[i * m_numAtoms * 3 + j * 3 + 2]); + mode.push_back(v); + } + m_vibDisplacements.push_back(mode); + } + } + cout << "Read " << m_vibDisplacements.size() << " vibrational modes\n"; } } diff --git a/avogadro/quantumio/gaussianfchk.h b/avogadro/quantumio/gaussianfchk.h index e2c3c4abea..f232c2ffbd 100644 --- a/avogadro/quantumio/gaussianfchk.h +++ b/avogadro/quantumio/gaussianfchk.h @@ -7,6 +7,8 @@ #define AVOGADRO_QUANTUMIO_GAUSSIANFCHK_H #include "avogadroquantumioexport.h" + +#include #include #include @@ -67,6 +69,8 @@ class AVOGADROQUANTUMIO_EXPORT GaussianFchk : public Io::FileFormat int m_electrons; int m_electronsAlpha; int m_electronsBeta; + int m_normalModes; + int m_numAtoms; unsigned char m_spin; signed char m_charge; unsigned int m_numBasisFunctions; @@ -87,8 +91,14 @@ class AVOGADROQUANTUMIO_EXPORT GaussianFchk : public Io::FileFormat MatrixX m_density; /// Total density matrix MatrixX m_spinDensity; /// Spin density matrix Core::ScfType m_scftype; + + Core::Array m_frequencies; + Core::Array m_IRintensities; + Core::Array m_RamanIntensities; + Core::Array> m_vibDisplacements; }; -} -} + +} // namespace QuantumIO +} // namespace Avogadro #endif From 66e0bfc3fee9856570dc69ede33e8c1bd38cbdb0 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 18 Oct 2023 17:47:27 -0400 Subject: [PATCH 2/3] Fix formatting Signed-off-by: Geoff Hutchison --- avogadro/quantumio/gaussianfchk.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp index 9f748a1269..2a2d1c792c 100644 --- a/avogadro/quantumio/gaussianfchk.cpp +++ b/avogadro/quantumio/gaussianfchk.cpp @@ -221,7 +221,8 @@ void GaussianFchk::processLine(std::istream& in) } } } else if (key == "Vib-Modes" && list.size() > 2) { - std::vector tmp = readArrayD(in, Core::lexicalCast(list[2]), 16); + std::vector tmp = + readArrayD(in, Core::lexicalCast(list[2]), 16); m_vibDisplacements.clear(); if (tmp.size() == m_numAtoms * 3 * m_normalModes) { for (unsigned int i = 0; i < m_normalModes; ++i) { From 4f0291a17df3881e14df1b2a95cc9f03dad5f9c7 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 18 Oct 2023 20:06:17 -0400 Subject: [PATCH 3/3] Fix code quality (too many tmp variables) Signed-off-by: Geoff Hutchison --- avogadro/quantumio/gaussianfchk.cpp | 55 ++++++++++------------------- 1 file changed, 19 insertions(+), 36 deletions(-) diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp index 2a2d1c792c..2daf34a2f5 100644 --- a/avogadro/quantumio/gaussianfchk.cpp +++ b/avogadro/quantumio/gaussianfchk.cpp @@ -95,6 +95,7 @@ void GaussianFchk::processLine(std::istream& in) string tmp = line.substr(43); vector list = Core::split(tmp, ' '); + std::vector tmpVec; // Big switch statement checking for various things we are interested in if (Core::contains(key, "RHF")) { @@ -115,13 +116,11 @@ void GaussianFchk::processLine(std::istream& in) m_electronsBeta = Core::lexicalCast(list[1]); } else if (key == "Number of basis functions" && list.size() > 1) { m_numBasisFunctions = Core::lexicalCast(list[1]); - cout << "Number of basis functions = " << m_numBasisFunctions << endl; + // cout << "Number of basis functions = " << m_numBasisFunctions << endl; } else if (key == "Atomic numbers" && list.size() > 2) { m_aNums = readArrayI(in, Core::lexicalCast(list[2])); if (static_cast(m_aNums.size()) != Core::lexicalCast(list[2])) cout << "Reading atomic numbers failed.\n"; - else - cout << "Reading atomic numbers succeeded.\n"; } // Now we get to the meat of it - coordinates of the atoms else if (key == "Current cartesian coordinates" && list.size() > 2) { @@ -145,15 +144,16 @@ void GaussianFchk::processLine(std::istream& in) } else if (key == "Alpha Orbital Energies") { if (m_scftype == Rhf) { m_orbitalEnergy = readArrayD(in, Core::lexicalCast(list[2]), 16); - cout << "MO energies, n = " << m_orbitalEnergy.size() << endl; + // cout << "MO energies, n = " << m_orbitalEnergy.size() << endl; } else if (m_scftype == Uhf) { m_alphaOrbitalEnergy = readArrayD(in, Core::lexicalCast(list[2]), 16); - cout << "Alpha MO energies, n = " << m_alphaOrbitalEnergy.size() << endl; + // cout << "Alpha MO energies, n = " << m_alphaOrbitalEnergy.size() << + // endl; } } else if (key == "Beta Orbital Energies") { if (m_scftype != Uhf) { - cout << "UHF detected. Reassigning Alpha properties." << endl; + // cout << "UHF detected. Reassigning Alpha properties." << endl; m_scftype = Uhf; m_alphaOrbitalEnergy = m_orbitalEnergy; m_orbitalEnergy = vector(); @@ -163,36 +163,22 @@ void GaussianFchk::processLine(std::istream& in) } m_betaOrbitalEnergy = readArrayD(in, Core::lexicalCast(list[2]), 16); - cout << "Beta MO energies, n = " << m_betaOrbitalEnergy.size() << endl; + // cout << "Beta MO energies, n = " << m_betaOrbitalEnergy.size() << endl; } else if (key == "Alpha MO coefficients" && list.size() > 2) { if (m_scftype == Rhf) { m_MOcoeffs = readArrayD(in, Core::lexicalCast(list[2]), 16); - if (static_cast(m_MOcoeffs.size()) == - Core::lexicalCast(list[2])) - cout << "MO coefficients, n = " << m_MOcoeffs.size() << endl; } else if (m_scftype == Uhf) { m_alphaMOcoeffs = readArrayD(in, Core::lexicalCast(list[2]), 16); - if (static_cast(m_alphaMOcoeffs.size()) == - Core::lexicalCast(list[2])) - cout << "Alpha MO coefficients, n = " << m_alphaMOcoeffs.size() << endl; } else { cout << "Error, alpha MO coefficients, n = " << m_MOcoeffs.size() << endl; } } else if (key == "Beta MO coefficients" && list.size() > 2) { m_betaMOcoeffs = readArrayD(in, Core::lexicalCast(list[2]), 16); - if (static_cast(m_betaMOcoeffs.size()) == - Core::lexicalCast(list[2])) - cout << "Beta MO coefficients, n = " << m_betaMOcoeffs.size() << endl; } else if (key == "Total SCF Density" && list.size() > 2) { - if (readDensityMatrix(in, Core::lexicalCast(list[2]), 16)) - cout << "SCF density matrix read in " << m_density.rows() << endl; - else + if (!readDensityMatrix(in, Core::lexicalCast(list[2]), 16)) cout << "Error reading in the SCF density matrix.\n"; } else if (key == "Spin SCF Density" && list.size() > 2) { - if (readSpinDensityMatrix(in, Core::lexicalCast(list[2]), 16)) - cout << "SCF spin density matrix read in " << m_spinDensity.rows() - << endl; - else + if (!readSpinDensityMatrix(in, Core::lexicalCast(list[2]), 16)) cout << "Error reading in the SCF spin density matrix.\n"; } else if (key == "Number of Normal Modes" && list.size() > 1) { m_normalModes = Core::lexicalCast(list[1]); @@ -202,41 +188,38 @@ void GaussianFchk::processLine(std::istream& in) m_RamanIntensities.clear(); unsigned threeN = m_numAtoms * 3; // degrees of freedom - std::vector tmp = - readArrayD(in, Core::lexicalCast(list[2]), 16); + tmpVec = readArrayD(in, Core::lexicalCast(list[2]), 16); // read in the first 3N-6 elements as frequencies for (unsigned int i = 0; i < m_normalModes; ++i) { - m_frequencies.push_back(tmp[i]); + m_frequencies.push_back(tmpVec[i]); } // skip to after threeN elements then read IR intensities for (unsigned int i = threeN; i < threeN + m_normalModes; ++i) { - m_IRintensities.push_back(tmp[i]); + m_IRintensities.push_back(tmpVec[i]); } // now check if we have Raman intensities - if (tmp[threeN + m_normalModes] != 0.0) { + if (tmpVec[threeN + m_normalModes] != 0.0) { for (unsigned int i = threeN + m_normalModes; i < threeN + 2 * m_normalModes; ++i) { - m_RamanIntensities.push_back(tmp[i]); + m_RamanIntensities.push_back(tmpVec[i]); } } } else if (key == "Vib-Modes" && list.size() > 2) { - std::vector tmp = - readArrayD(in, Core::lexicalCast(list[2]), 16); + tmpVec = readArrayD(in, Core::lexicalCast(list[2]), 16); m_vibDisplacements.clear(); - if (tmp.size() == m_numAtoms * 3 * m_normalModes) { + if (tmpVec.size() == m_numAtoms * 3 * m_normalModes) { for (unsigned int i = 0; i < m_normalModes; ++i) { Core::Array mode; for (unsigned int j = 0; j < m_numAtoms; ++j) { - Vector3 v(tmp[i * m_numAtoms * 3 + j * 3], - tmp[i * m_numAtoms * 3 + j * 3 + 1], - tmp[i * m_numAtoms * 3 + j * 3 + 2]); + Vector3 v(tmpVec[i * m_numAtoms * 3 + j * 3], + tmpVec[i * m_numAtoms * 3 + j * 3 + 1], + tmpVec[i * m_numAtoms * 3 + j * 3 + 2]); mode.push_back(v); } m_vibDisplacements.push_back(mode); } } - cout << "Read " << m_vibDisplacements.size() << " vibrational modes\n"; } }