From 92943bc6a9177279c783ee23f79ef1db90793e67 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Fri, 1 Sep 2023 21:57:54 -0400 Subject: [PATCH] Add support for reading vibrations and IR Signed-off-by: Geoff Hutchison --- avogadro/quantumio/orca.cpp | 277 ++++++++++++++++++++++++++++++++---- avogadro/quantumio/orca.h | 22 ++- 2 files changed, 268 insertions(+), 31 deletions(-) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 6f581fc8c4..6f1e08ec7c 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -9,8 +9,8 @@ #include #include -#include #include +#include using std::string; using std::vector; @@ -59,12 +59,28 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) return false; } + // add the partial charges + if (m_partialCharges.size() > 0) { + for (auto it = m_partialCharges.begin(); it != m_partialCharges.end(); + ++it) { + + std::cerr << " Adding partial charges of type " << it->first << "\n"; + molecule.setPartialCharges(it->first, it->second); + } + } + + std::cerr << " frequencies " << m_frequencies.size() << "\n"; + std::cerr << " vibDisplacements " << m_vibDisplacements.size() << "\n"; + std::cerr << " IRintensities " << m_IRintensities.size() << "\n"; + if (m_frequencies.size() > 0 && m_frequencies.size() == m_vibDisplacements.size() && - m_frequencies.size() == m_intensities.size()) { + m_frequencies.size() == m_IRintensities.size()) { molecule.setVibrationFrequencies(m_frequencies); - molecule.setVibrationIRIntensities(m_intensities); + molecule.setVibrationIRIntensities(m_IRintensities); molecule.setVibrationLx(m_vibDisplacements); + if (m_RamanIntensities.size()) + molecule.setVibrationRamanIntensities(m_RamanIntensities); } // Do simple bond perception. @@ -88,6 +104,7 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) string key = Core::trimmed(line); vector list; int nGTOs = 0; + float vibScaling = 1.0f; if (Core::contains(key, "CARTESIAN COORDINATES (A.U.)")) { m_coordFactor = 1.; // leave the coords in BOHR .... @@ -137,10 +154,50 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) list = Core::split(key, ' '); m_electrons = Core::lexicalCast(list[5]); } else if (Core::contains(key, "SPIN UP ORBITALS") && !m_openShell) { - m_openShell = true; // not yet implemented + m_openShell = true; // TODO } else if (Core::contains(key, "MOLECULAR ORBITALS")) { m_currentMode = MO; getline(in, key); //------------ + } else if (Core::contains(key, "ATOMIC CHARGES")) { + m_currentMode = Charges; + // figure out what type of charges we have + list = Core::split(key, ' '); + if (list.size() > 2) { + m_chargeType = Core::trimmed(list[0]); // e.g. MULLIKEN or LOEWDIN + } + getline(in, key); // skip ------------ + } else if (Core::contains(key, "VIBRATIONAL FREQUENCIES")) { + m_currentMode = Frequencies; + getline(in, key); // skip ------------ + getline(in, key); // skip blank line + getline(in, key); // scaling factor + // Scaling factor for frequencies = 1.000000000 + list = Core::split(key, ' '); + if (list.size() > 6) + vibScaling = Core::lexicalCast(list[5]); + getline(in, key); // skip blank line + } else if (Core::contains(key, "NORMAL MODES")) { + m_currentMode = VibrationalModes; + + getline(in, key); // skip ------------ + getline(in, key); // skip blank line + getline(in, key); // skip comment + getline(in, key); // skip more comments + getline(in, key); // skip even more comment + getline(in, key); // skip blank line + } else if (Core::contains(key, "IR SPECTRUM")) { + m_currentMode = IR; + getline(in, key); // skip ------------ + getline(in, key); // skip blank line + getline(in, key); // skip column titles + getline(in, key); // skip more column titles + getline(in, key); // skip ------------ + } else if (Core::contains(key, "RAMAN SPECTRUM")) { + m_currentMode = Raman; + getline(in, key); // skip ------------ + getline(in, key); // skip blank line + getline(in, key); // skip column titles + getline(in, key); // skip ------------ } else { vector> columns; @@ -158,9 +215,10 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) if (list.size() < 8) { break; } - Eigen::Vector3d pos(Core::lexicalCast(list[5]) * m_coordFactor, - Core::lexicalCast(list[6]) * m_coordFactor, - Core::lexicalCast(list[7]) * m_coordFactor); + Eigen::Vector3d pos( + Core::lexicalCast(list[5]) * m_coordFactor, + Core::lexicalCast(list[6]) * m_coordFactor, + Core::lexicalCast(list[7]) * m_coordFactor); m_atomNums.push_back(Core::lexicalCast(list[2])); m_atomPos.push_back(pos); @@ -172,6 +230,156 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) m_currentMode = NotParsing; break; } + case Charges: { + // should start at the first atom + if (key.empty()) + break; + + Eigen::MatrixXd charges(m_atomNums.size(), 1); + charges.setZero(); + + list = Core::split(key, ' '); + while (!key.empty()) { + if (list.size() != 4) { + break; + } + // e.g. 0 O : -0.714286 + int atomIndex = Core::lexicalCast(list[0]); + double charge = Core::lexicalCast(list[3]); + charges(atomIndex, 0) = charge; + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + m_partialCharges[m_chargeType] = charges; + m_currentMode = NotParsing; + break; + } + case Frequencies: { + // should start at the first frequency - include zeros + if (key.empty()) + break; + list = Core::split(key, ' '); + while (!key.empty()) { + if (list.size() != 3) { + break; + } + // e.g. 0: 0.00 cm**-1 + double freq = Core::lexicalCast(list[1]); + m_frequencies.push_back(freq); + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + + // okay, now set up the normal mode arrays + m_vibDisplacements.resize(m_frequencies.size()); + m_IRintensities.reserve(m_frequencies.size()); + // we don't bother with Raman, because that's less common + for (unsigned int i = 0; i < m_frequencies.size(); i++) { + m_vibDisplacements[i].resize(m_atomNums.size()); + for (unsigned int j = 0; j < m_atomNums.size(); j++) + m_vibDisplacements[i].push_back(Eigen::Vector3d()); + } + + m_currentMode = NotParsing; + break; + } + case VibrationalModes: { + if (key.empty()) + break; + list = Core::split(key, ' '); + vector modeIndex; + while (!key.empty()) { + // first we get a set of column numbers + // e.g. 1 2 3 4 5 6 7 8 9 10 + modeIndex.clear(); + for (unsigned int i = 0; i < list.size(); i++) { + modeIndex.push_back(Core::lexicalCast(list[i])); + } + // now we read the displacements .. there should be 3N lines + // x,y,z for each atom + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + for (unsigned int i = 0; i < 3*m_atomNums.size(); i++) { + unsigned int atomIndex = i / 3; + unsigned int coordIndex = i % 3; + for (unsigned int j = 0; j < modeIndex.size(); j++) { + m_vibDisplacements[modeIndex[j]][atomIndex][coordIndex] = + Core::lexicalCast(list[j+1]); + } + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + + } + + m_currentMode = NotParsing; + break; + } + case IR: { + if (key.empty()) + break; + list = Core::split(key, ' '); + while (!key.empty()) { + // e.g. 6: 1711.76 0.014736 74.47 0.002686 (-0.021704 0.027180 + // 0.038427) + if (list.size() != 8) { + break; + } + // the first entry might be 5 or 6 because of removed rotations / + // translations + int index = Core::lexicalCast(list[0]); + if (m_IRintensities.size() == 0 && index > 0) { + while (m_IRintensities.size() < index + 1) { + m_IRintensities.push_back(0.0); + } + } + double intensity = Core::lexicalCast(list[3]); + m_IRintensities.push_back(intensity); + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + + m_currentMode = NotParsing; + break; + } + case Raman: { + if (key.empty()) + break; + list = Core::split(key, ' '); + while (!key.empty()) { + // e.g. 6: 76.62 0.000000 0.465517 + if (list.size() != 4) { + break; + } + // the first entry might be 5 or 6 because of removed rotations / + // translations + int index = Core::lexicalCast(list[0]); + if (m_RamanIntensities.size() == 0 && index > 0) { + while (m_RamanIntensities.size() < index) { + m_RamanIntensities.push_back(0.0); + } + } + // index, frequency, activity, depolarization + double activity = Core::lexicalCast(list[2]); + m_RamanIntensities.push_back(activity); + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + + m_currentMode = NotParsing; + break; + } case GTO: { // // should start at the first newGTO if (key.empty()) @@ -228,7 +436,8 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) list = Core::split(key, ' '); } - m_orcaShellTypes.push_back(std::vector(shellTypes.size())); + m_orcaShellTypes.push_back( + std::vector(shellTypes.size())); m_orcaShellTypes.at(nGTOs) = shellTypes; m_orcaNumShells.push_back(std::vector(shellFunctions.size())); m_orcaNumShells.at(nGTOs) = shellFunctions; @@ -439,15 +648,15 @@ void ORCAOutput::load(GaussianSet* basis) int nSP = 0; // number of SP shells for (unsigned int i = 0; i < m_shellTypes.size(); ++i) { // Handle the SP case separately - this should possibly be a distinct type - if (m_shellTypes.at(i) == Core::GaussianSet::SP) { + if (m_shellTypes.at(i) == GaussianSet::SP) { // SP orbital type - currently have to unroll into two shells int tmpGTO = nGTO; - int s = basis->addBasis(m_shelltoAtom.at(i) - 1, Core::GaussianSet::S); + int s = basis->addBasis(m_shelltoAtom.at(i) - 1, GaussianSet::S); for (int j = 0; j < m_shellNums.at(i); ++j) { basis->addGto(s, m_c.at(nGTO), m_a.at(nGTO)); ++nGTO; } - int p = basis->addBasis(m_shelltoAtom.at(i) - 1, Core::GaussianSet::P); + int p = basis->addBasis(m_shelltoAtom.at(i) - 1, GaussianSet::P); for (int j = 0; j < m_shellNums.at(i); ++j) { basis->addGto(p, m_csp.at(nSP), m_a.at(tmpGTO)); ++tmpGTO; @@ -470,21 +679,37 @@ void ORCAOutput::load(GaussianSet* basis) basis->generateDensityMatrix(); } -Core::GaussianSet::orbital ORCAOutput::orbitalIdx(std::string txt) { - if (txt == "S") return Core::GaussianSet::S; - if (txt == "SP") return Core::GaussianSet::SP; - if (txt == "P") return Core::GaussianSet::P; - if (txt == "D") return Core::GaussianSet::D5; //// orca only uses Spherical - 5 d components - if (txt == "D5") return Core::GaussianSet::D5; - if (txt == "F") return Core::GaussianSet::F7; //// orca only uses Spherical - 7 f components - if (txt == "F7") return Core::GaussianSet::F7; - if (txt == "G") return Core::GaussianSet::G9; //// orca only uses Spherical - 9 g components - if (txt == "G9") return Core::GaussianSet::G9; - if (txt == "H") return Core::GaussianSet::H11; //// orca only uses Spherical - 11 h components - if (txt == "H11") return Core::GaussianSet::H11; - if (txt == "I") return Core::GaussianSet::I13; //// orca only uses Spherical - 13 i components - if (txt == "I13") return Core::GaussianSet::I13; - return Core::GaussianSet::UU; +GaussianSet::orbital ORCAOutput::orbitalIdx(std::string txt) +{ + if (txt == "S") + return GaussianSet::S; + if (txt == "SP") + return GaussianSet::SP; + if (txt == "P") + return GaussianSet::P; + if (txt == "D") + return GaussianSet::D5; //// orca only uses Spherical - 5 d components + if (txt == "D5") + return GaussianSet::D5; + if (txt == "F") + return GaussianSet::F7; //// orca only uses Spherical - 7 f components + if (txt == "F7") + return GaussianSet::F7; + if (txt == "G") + return GaussianSet::G9; //// orca only uses Spherical - 9 g components + if (txt == "G9") + return GaussianSet::G9; + if (txt == "H") + return GaussianSet::H11; //// orca only uses Spherical - 11 h + /// components + if (txt == "H11") + return GaussianSet::H11; + if (txt == "I") + return GaussianSet::I13; //// orca only uses Spherical - 13 i + /// components + if (txt == "I13") + return GaussianSet::I13; + return GaussianSet::UU; } } // namespace Avogadro::QuantumIO diff --git a/avogadro/quantumio/orca.h b/avogadro/quantumio/orca.h index fa8e234e21..d22a44d26b 100644 --- a/avogadro/quantumio/orca.h +++ b/avogadro/quantumio/orca.h @@ -11,6 +11,7 @@ #include #include +#include #include namespace Avogadro { @@ -51,7 +52,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat void load(Core::GaussianSet* basis); // OrcaStuff - void orcaWarningMessage(const std::string &m); + void orcaWarningMessage(const std::string& m); Core::GaussianSet::orbital orbitalIdx(std::string txt); bool m_orcaSuccess; @@ -63,17 +64,24 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat std::vector shellFunctions; std::vector shellTypes; - std::vector > m_orcaNumShells; - std::vector > m_orcaShellTypes; + std::vector> m_orcaNumShells; + std::vector> m_orcaShellTypes; int m_nGroups; - std::vector *> *> m_basisFunctions; + std::vector*>*> m_basisFunctions; enum mode { Atoms, GTO, MO, + Charges, + Frequencies, + VibrationalModes, + IR, + Raman, + Electronic, + NMR, NotParsing, Unrecognized }; @@ -98,8 +106,12 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat std::vector m_orbitalEnergy; std::vector m_MOcoeffs; + std::string m_chargeType; + std::map m_partialCharges; + Core::Array m_frequencies; - Core::Array m_intensities; + Core::Array m_IRintensities; + Core::Array m_RamanIntensities; Core::Array> m_vibDisplacements; };