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