From 209799551a0ed634b1e27456b4ab6e5b83170f2f Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 2 Sep 2023 18:53:04 -0400 Subject: [PATCH] Add orca import (#1326) Add ORCA output file import (porting 1.2 code) * Initial support for geometry import * Add support for reading vibrations and IR * Support open shell calculations and orbital energies Signed-off-by: Geoff Hutchison --- avogadro/core/gaussianset.cpp | 27 +- avogadro/core/gaussianset.h | 24 +- avogadro/qtplugins/surfaces/surfaces.cpp | 143 +++-- avogadro/quantumio/CMakeLists.txt | 2 + avogadro/quantumio/orca.cpp | 762 +++++++++++++++++++++++ avogadro/quantumio/orca.h | 124 ++++ 6 files changed, 983 insertions(+), 99 deletions(-) create mode 100644 avogadro/quantumio/orca.cpp create mode 100644 avogadro/quantumio/orca.h diff --git a/avogadro/core/gaussianset.cpp b/avogadro/core/gaussianset.cpp index 8af5eae3a6..016dee32a2 100644 --- a/avogadro/core/gaussianset.cpp +++ b/avogadro/core/gaussianset.cpp @@ -22,9 +22,7 @@ GaussianSet::GaussianSet() : m_numMOs(0), m_init(false) m_scfType = Rhf; } -GaussianSet::~GaussianSet() -{ -} +GaussianSet::~GaussianSet() {} unsigned int GaussianSet::addBasis(unsigned int atom, orbital type) { @@ -184,12 +182,6 @@ bool GaussianSet::setSpinDensityMatrix(const MatrixX& m) return true; } -bool GaussianSet::generateDensityMatrix() -{ - // FIXME: Finish me! - return true; -} - unsigned int GaussianSet::molecularOrbitalCount(ElectronType type) { size_t index(0); @@ -479,7 +471,7 @@ void GaussianSet::initCalculation() m_init = true; } -bool GaussianSet::generateDensity() +bool GaussianSet::generateDensityMatrix() { if (m_scfType == Unknown) return false; @@ -496,9 +488,11 @@ bool GaussianSet::generateDensity() m_density(jBasis, iBasis) += 2.0 * icoeff * jcoeff; m_density(iBasis, jBasis) = m_density(jBasis, iBasis); } - cout << iBasis << ", " << jBasis << ": " << m_density(iBasis, jBasis) - << endl; + // cout << iBasis << ", " << jBasis << ": " << + // m_density(iBasis, jBasis) + // << endl; break; + case Rohf: // ROHF is handled similarly to UHF case Uhf: for (unsigned int iaMO = 0; iaMO < m_electrons[0]; ++iaMO) { double icoeff = m_moMatrix[0](iBasis, iaMO); @@ -512,8 +506,9 @@ bool GaussianSet::generateDensity() m_density(jBasis, iBasis) += icoeff * jcoeff; m_density(iBasis, jBasis) = m_density(jBasis, iBasis); } - cout << iBasis << ", " << jBasis << ": " << m_density(iBasis, jBasis) - << endl; + // cout << iBasis << ", " << jBasis << ": " << + // m_density(iBasis, jBasis) + // << endl; break; default: cout << "Unhandled scf type:" << m_scfType << endl; @@ -523,7 +518,7 @@ bool GaussianSet::generateDensity() return true; } -bool GaussianSet::generateSpinDensity() +bool GaussianSet::generateSpinDensityMatrix() { if (m_scfType != Uhf) return false; @@ -551,4 +546,4 @@ bool GaussianSet::generateSpinDensity() return true; } -} // End namespace Avogadro +} // namespace Avogadro::Core diff --git a/avogadro/core/gaussianset.h b/avogadro/core/gaussianset.h index 15ba4dd664..76c31994fd 100644 --- a/avogadro/core/gaussianset.h +++ b/avogadro/core/gaussianset.h @@ -173,6 +173,13 @@ class AVOGADROCORE_EXPORT GaussianSet : public BasisSet */ bool generateDensityMatrix(); + /** + * @brief Generate the spin density matrix if we have the required + * information. + * @return True on success, false on failure. + */ + bool generateSpinDensityMatrix(); + /** * @return The number of molecular orbitals in the GaussianSet. */ @@ -357,22 +364,9 @@ class AVOGADROCORE_EXPORT GaussianSet : public BasisSet ScfType m_scfType; std::string m_functionalName; - - /** - * @brief Generate the density matrix if we have the required information. - * @return True on success, false on failure. - */ - bool generateDensity(); - - /** - * @brief Generate the spin density matrix if we have the required - * information. - * @return True on success, false on failure. - */ - bool generateSpinDensity(); }; -} // End Core namespace -} // End Avogadro namespace +} // namespace Core +} // namespace Avogadro #endif diff --git a/avogadro/qtplugins/surfaces/surfaces.cpp b/avogadro/qtplugins/surfaces/surfaces.cpp index 79e878f894..bff82f66f6 100644 --- a/avogadro/qtplugins/surfaces/surfaces.cpp +++ b/avogadro/qtplugins/surfaces/surfaces.cpp @@ -42,15 +42,16 @@ namespace { #include #include #include +#include +#include +#include #include #include #include #include #include #include -#include -#include #include #include #include @@ -81,7 +82,7 @@ Surfaces::Surfaces(QObject* p) : ExtensionPlugin(p), d(new PIMPL()) connect(action, SIGNAL(triggered()), SLOT(surfacesActivated())); connect(&m_displayMeshWatcher, SIGNAL(finished()), SLOT(displayMesh())); connect(&m_performEDTStepWatcher, SIGNAL(finished()), SLOT(performEDTStep())); - + m_actions.push_back(action); // Register quantum file formats @@ -92,6 +93,7 @@ Surfaces::Surfaces(QObject* p) : ExtensionPlugin(p), d(new PIMPL()) Io::FileFormatManager::registerFormat(new QuantumIO::MopacAux); Io::FileFormatManager::registerFormat(new QuantumIO::NWChemJson); Io::FileFormatManager::registerFormat(new QuantumIO::NWChemLog); + Io::FileFormatManager::registerFormat(new QuantumIO::ORCAOutput); } Surfaces::~Surfaces() @@ -157,10 +159,12 @@ void Surfaces::surfacesActivated() } m_dialog->setupSteps(m_molecule->coordinate3dCount()); - const auto identifiers = Calc::ChargeManager::instance().identifiersForMolecule(*m_molecule); + const auto identifiers = + Calc::ChargeManager::instance().identifiersForMolecule(*m_molecule); std::set> chargeModels; - for (const auto &identifier: identifiers) - chargeModels.emplace(Calc::ChargeManager::instance().nameForModel(identifier), identifier); + for (const auto& identifier : identifiers) + chargeModels.emplace( + Calc::ChargeManager::instance().nameForModel(identifier), identifier); m_dialog->setupModels(chargeModels); m_dialog->show(); @@ -179,8 +183,7 @@ float Surfaces::resolution() case SolventExcluded: minimum = 0.1; break; - default: - ; + default:; } r = std::max(minimum, std::min(maximum, r)); @@ -218,7 +221,10 @@ void Surfaces::calculateSurface() } } -float inline square(float x) { return x * x; } +float inline square(float x) +{ + return x * x; +} void Surfaces::calculateEDT() { @@ -237,13 +243,14 @@ void Surfaces::calculateEDT() // first, make a list of all atom positions and radii Array atomPositions = m_molecule->atomPositions3d(); - auto *atoms = new std::vector>(); + auto* atoms = new std::vector>(); double max_radius = probeRadius; QtGui::RWLayerManager layerManager; for (size_t i = 0; i < m_molecule->atomCount(); i++) { if (!layerManager.visible(m_molecule->layer(i))) continue; - auto radius = Core::Elements::radiusVDW(m_molecule->atomicNumber(i)) + probeRadius; + auto radius = + Core::Elements::radiusVDW(m_molecule->atomicNumber(i)) + probeRadius; atoms->emplace_back(atomPositions[i], radius); if (radius > max_radius) max_radius = radius; @@ -257,39 +264,41 @@ void Surfaces::calculateEDT() const Vector3 min = m_cube->min(); // then, for each atom, set cubes around it up to a certain radius - QFuture innerFuture = QtConcurrent::map(*atoms, [=](std::pair &in) { - double startPosX = in.first(0) - in.second; - double endPosX = in.first(0) + in.second; - int startIndexX = (startPosX - min(0)) / res; - int endIndexX = (endPosX - min(0)) / res + 1; - for (int indexX = startIndexX; indexX < endIndexX; indexX++) { - double posX = indexX * res + min(0); - double radiusXsq = square(in.second) - square(posX - in.first(0)); - if (radiusXsq < 0.0) - continue; - double radiusX = sqrt(radiusXsq); - double startPosY = in.first(1) - radiusX; - double endPosY = in.first(1) + radiusX; - int startIndexY = (startPosY - min(1)) / res; - int endIndexY = (endPosY - min(1)) / res + 1; - for (int indexY = startIndexY; indexY < endIndexY; indexY++) { - double posY = indexY * res + min(1); - double lengthXYsq = square(radiusX) - square(posY - in.first(1)); - if (lengthXYsq < 0.0) + QFuture innerFuture = + QtConcurrent::map(*atoms, [=](std::pair& in) { + double startPosX = in.first(0) - in.second; + double endPosX = in.first(0) + in.second; + int startIndexX = (startPosX - min(0)) / res; + int endIndexX = (endPosX - min(0)) / res + 1; + for (int indexX = startIndexX; indexX < endIndexX; indexX++) { + double posX = indexX * res + min(0); + double radiusXsq = square(in.second) - square(posX - in.first(0)); + if (radiusXsq < 0.0) continue; - double lengthXY = sqrt(lengthXYsq); - double startPosZ = in.first(2) - lengthXY; - double endPosZ = in.first(2) + lengthXY; - int startIndexZ = (startPosZ - min(2)) / res; - int endIndexZ = (endPosZ - min(2)) / res + 1; - m_cube->fillStripe(indexX, indexY, startIndexZ, endIndexZ - 1, 1.0f); + double radiusX = sqrt(radiusXsq); + double startPosY = in.first(1) - radiusX; + double endPosY = in.first(1) + radiusX; + int startIndexY = (startPosY - min(1)) / res; + int endIndexY = (endPosY - min(1)) / res + 1; + for (int indexY = startIndexY; indexY < endIndexY; indexY++) { + double posY = indexY * res + min(1); + double lengthXYsq = square(radiusX) - square(posY - in.first(1)); + if (lengthXYsq < 0.0) + continue; + double lengthXY = sqrt(lengthXYsq); + double startPosZ = in.first(2) - lengthXY; + double endPosZ = in.first(2) + lengthXY; + int startIndexZ = (startPosZ - min(2)) / res; + int endIndexZ = (endPosZ - min(2)) / res + 1; + m_cube->fillStripe(indexX, indexY, startIndexZ, endIndexZ - 1, + 1.0f); + } } - } - }); - + }); + innerFuture.waitForFinished(); }); - + // SolventExcluded requires an extra pass if (m_dialog->surfaceType() == SolventExcluded) { m_performEDTStepWatcher.setFuture(future); @@ -303,14 +312,14 @@ void Surfaces::performEDTStep() QFuture future = QtConcurrent::run([=]() { const double probeRadius = 1.4; const double scaledProbeRadius = probeRadius / resolution(); - + // make a list of all "outside" cubes in contact with an "inside" cube // these are the only ones that can be "nearest" to an "inside" cube Array relativePositions; // also make a list of all "inside" cubes - auto *insideIndices = new std::vector; + auto* insideIndices = new std::vector; Vector3i size = m_cube->dimensions(); - relativePositions.reserve(size(0) * size(1) * 4); // O(n^2) + relativePositions.reserve(size(0) * size(1) * 4); // O(n^2) insideIndices->reserve(size(0) * size(1) * size(2)); // O(n^3) for (int z = 0; z < size(2); z++) { int zp = std::max(z - 1, 0); @@ -325,31 +334,27 @@ void Surfaces::performEDTStep() } int xp = std::max(x - 1, 0); int xn = std::min(x + 1, size(0) - 1); - if (m_cube->value(xp, y, z) > 0.0 - || m_cube->value(xn, y, z) > 0.0 - || m_cube->value(x, yp, z) > 0.0 - || m_cube->value(x, yn, z) > 0.0 - || m_cube->value(x, y, zp) > 0.0 - || m_cube->value(x, y, zn) > 0.0 - ) { + if (m_cube->value(xp, y, z) > 0.0 || m_cube->value(xn, y, z) > 0.0 || + m_cube->value(x, yp, z) > 0.0 || m_cube->value(x, yn, z) > 0.0 || + m_cube->value(x, y, zp) > 0.0 || m_cube->value(x, y, zn) > 0.0) { relativePositions.push_back(Vector3(x, y, z)); } } } } - + // pass the list to a NeighborPerceiver so it's faster to look up NeighborPerceiver perceiver(relativePositions, scaledProbeRadius); - + // now, exclude all "inside" cubes too close to any "outside" cube - thread_local Array *neighbors = nullptr; - QFuture innerFuture = QtConcurrent::map(*insideIndices, [=](Vector3i &in) { + thread_local Array* neighbors = nullptr; + QFuture innerFuture = QtConcurrent::map(*insideIndices, [=](Vector3i& in) { Vector3 pos = in.cast(); if (neighbors == nullptr) neighbors = new Array; perceiver.getNeighborsInclusiveInPlace(*neighbors, pos); - for (Index neighbor: *neighbors) { - const Vector3 &npos = relativePositions[neighbor]; + for (Index neighbor : *neighbors) { + const Vector3& npos = relativePositions[neighbor]; float distance = (npos - pos).norm(); if (distance <= scaledProbeRadius) { m_cube->setValue(in(0), in(1), in(2), -1.0f); @@ -357,10 +362,10 @@ void Surfaces::performEDTStep() } } }); - + innerFuture.waitForFinished(); }); - + m_displayMeshWatcher.setFuture(future); } @@ -521,7 +526,8 @@ void Surfaces::displayMesh() m_meshGenerator2 = new QtGui::MeshGenerator; connect(m_meshGenerator2, SIGNAL(finished()), SLOT(meshFinished())); } - m_meshGenerator2->initialize(m_cube, m_mesh2, m_isoValue, m_smoothingPasses, true); + m_meshGenerator2->initialize(m_cube, m_mesh2, m_isoValue, m_smoothingPasses, + true); // Start the mesh generation - this needs an improved mutex with a read lock // to function as expected. Write locks are exclusive, read locks can have @@ -533,7 +539,8 @@ void Surfaces::displayMesh() m_meshesLeft = 2; } -Core::Color3f Surfaces::chargeGradient(double value, double clamp, ColormapType colormap) const +Core::Color3f Surfaces::chargeGradient(double value, double clamp, + ColormapType colormap) const { // okay, typically color scales have blue at the bottom, red at the top. // so we need to invert, so blue is positive charge, red is negative charge. @@ -587,22 +594,22 @@ void Surfaces::colorMeshByPotential() { const auto model = m_dialog->colorModel().toStdString(); const auto colormap = getColormapFromString(m_dialog->colormapName()); - + const auto positionsf = m_mesh1->vertices(); Core::Array positions(positionsf.size()); std::transform(positionsf.begin(), positionsf.end(), positions.begin(), - [](const Vector3f &pos) { return pos.cast(); } - ); - const auto potentials = Calc::ChargeManager::instance().potentials(model, *m_molecule, positions); - + [](const Vector3f& pos) { return pos.cast(); }); + const auto potentials = + Calc::ChargeManager::instance().potentials(model, *m_molecule, positions); + double minPotential = *std::min_element(potentials.begin(), potentials.end()); double maxPotential = *std::max_element(potentials.begin(), potentials.end()); double clamp = std::max(std::abs(minPotential), std::abs(maxPotential)); - + Core::Array colors(positions.size()); for (size_t i = 0; i < potentials.size(); i++) colors[i] = chargeGradient(potentials[i], clamp, colormap); - + m_mesh1->setColors(colors); } @@ -804,4 +811,4 @@ void Surfaces::movieFrame() } } -} // namespace Avogadro +} // namespace Avogadro::QtPlugins diff --git a/avogadro/quantumio/CMakeLists.txt b/avogadro/quantumio/CMakeLists.txt index b8e0f4f8d4..6ae5e26e7d 100644 --- a/avogadro/quantumio/CMakeLists.txt +++ b/avogadro/quantumio/CMakeLists.txt @@ -8,6 +8,7 @@ avogadro_headers(QuantumIO mopacaux.h nwchemjson.h nwchemlog.h + orca.h ) # Source files for our data. @@ -19,6 +20,7 @@ target_sources(QuantumIO PRIVATE mopacaux.cpp nwchemjson.cpp nwchemlog.cpp + orca.cpp ) avogadro_add_library(QuantumIO) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp new file mode 100644 index 0000000000..0a69e2a49f --- /dev/null +++ b/avogadro/quantumio/orca.cpp @@ -0,0 +1,762 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "orca.h" + +#include +#include +#include + +#include +#include +#include + +using std::regex; +using std::string; +using std::vector; + +namespace Avogadro::QuantumIO { + +using Core::Array; +using Core::Atom; +using Core::GaussianSet; + +ORCAOutput::ORCAOutput() {} + +ORCAOutput::~ORCAOutput() {} + +std::vector ORCAOutput::fileExtensions() const +{ + std::vector extensions; + extensions.emplace_back("log"); + extensions.emplace_back("out"); + return extensions; +} + +std::vector ORCAOutput::mimeTypes() const +{ + return std::vector(); +} + +bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) +{ + // Read the log file line by line + auto* basis = new GaussianSet; + + while (!in.eof()) + processLine(in, basis); + + // Set up the molecule + int nAtom = 0; + for (unsigned int i = 0; i < m_atomNums.size(); i++) { + Vector3 pos = m_atomPos[i] * BOHR_TO_ANGSTROM; + molecule.addAtom(static_cast(m_atomNums[nAtom++]), pos); + } + + if (0 == molecule.atomCount()) { + appendError("Could not find any atomic coordinates! Are you sure this is " + "an ORCA output file?"); + return false; + } + + // add the partial charges + if (m_partialCharges.size() > 0) { + for (auto it = m_partialCharges.begin(); it != m_partialCharges.end(); + ++it) { + molecule.setPartialCharges(it->first, it->second); + } + } + + 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(); + + molecule.setBasisSet(basis); + basis->setMolecule(&molecule); + load(basis); + + return true; +} + +void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) +{ + // First truncate the line, remove trailing white space and check + string line; + if (!getline(in, line) || Core::trimmed(line).empty()) + return; + + 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 .... + m_currentMode = Atoms; + m_atomPos.clear(); + m_atomNums.clear(); + m_atomLabel.clear(); + getline(in, key); // skip ----- line + getline(in, key); // column titles + } else if (Core::contains(key, "BASIS SET INFORMATION")) { + if (!Core::contains(key, "AUXILIARY")) { // skip auxiliary basis set infos + m_currentMode = GTO; + getline(in, key); // skip ----- line + + // Number of groups of distinct atoms + getline(in, key); + list = Core::split(key, ' '); + if (list.size() > 3) { + m_nGroups = Core::lexicalCast(list[2]); + } else { + return; + } + getline(in, key); // skip blank line + for (int i = 0; i < m_nGroups; ++i) { + getline(in, key); // skip group information + } + getline(in, key); // skip blank line + for (unsigned int i = 0; i < m_atomNums.size(); ++i) { + getline(in, key); // skip group information + } + + // now skip + // blank line + // ---------------------------- + // # Basis set for element : x + // ---------------------------- + // blank line + for (unsigned int i = 0; i < 6; ++i) { + getline(in, key); + } + } + } else if (Core::contains(key, "TOTAL NUMBER OF BASIS SET")) { + m_currentMode = NotParsing; // no longer reading GTOs + } else if (Core::contains(key, "NUMBER OF CARTESIAN GAUSSIAN BASIS")) { + m_currentMode = NotParsing; // no longer reading GTOs + } else if (Core::contains(key, "Number of Electrons")) { + list = Core::split(key, ' '); + m_electrons = Core::lexicalCast(list[5]); + } else if (Core::contains(key, "ORBITAL ENERGIES")) { + m_currentMode = OrbitalEnergies; + getline(in, key); // skip ------------ + getline(in, key); // check if SPIN UP ORBITALS are present + if (Core::contains(key, "SPIN UP ORBITALS")) { + m_openShell = true; + m_readBeta = false; + } else { + m_openShell = false; + m_readBeta = false; + } + getline(in, key); // skip column titles + } else if (Core::contains(key, "SPIN DOWN ORBITALS")) { + m_currentMode = OrbitalEnergies; + m_openShell = true; + m_readBeta = true; + getline(in, key); // skip column headers + } 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; + unsigned int numColumns, numRows; + numColumns = 0; + numRows = 0; + // parsing a line -- what mode are we in? + + switch (m_currentMode) { + case Atoms: { + if (key.empty()) + break; + list = Core::split(key, ' '); + while (!key.empty()) { + 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); + + m_atomNums.push_back(Core::lexicalCast(list[2])); + m_atomPos.push_back(pos); + m_atomLabel.push_back(Core::trimmed(list[1])); + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + 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 OrbitalEnergies: { + // should start at the first orbital + if (!m_readBeta) + m_orbitalEnergy.clear(); + else + m_betaOrbitalEnergy.clear(); + + if (key.empty()) + break; + list = Core::split(key, ' '); + while (!key.empty()) { + if (list.size() != 4) { + break; + } + + // energy in Hartree + double energy = Core::lexicalCast(list[2]); + if (!m_readBeta) + m_orbitalEnergy.push_back(energy); + else + m_betaOrbitalEnergy.push_back(energy); + + getline(in, key); + key = Core::trimmed(key); + list = Core::split(key, ' '); + } + 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.resize(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_IRintensities[i] = 0.0; + 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() < 7) { + break; + } + // the first entry might be 5 or 6 because of removed rotations / + // translations + int index = Core::lexicalCast(list[0]); + double intensity = Core::lexicalCast(list[3]); + m_IRintensities[index] = 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()) + break; + nGTOs = 0; + list = Core::split(key, ' '); + int nShells; + // init all vectors etc. + m_basisAtomLabel.clear(); + m_orcaNumShells.resize(0); + m_basisFunctions.resize(0); + m_orcaShellTypes.resize(0); + + m_a.resize(0); + m_c.resize(0); + m_shellNums.resize(0); + m_shellTypes.resize(0); + m_shelltoAtom.resize(0); + while (Core::trimmed(list[0]) == "NewGTO") { + m_basisAtomLabel.push_back(Core::trimmed(list[1])); + + getline(in, key); + key = Core::trimmed(key); + + list = Core::split(key, ' '); + + nShells = 0; + m_basisFunctions.push_back( + new std::vector*>); + shellFunctions.resize(0); + shellTypes.resize(0); + while (Core::trimmed(list[0]) != "end;") { + + int nFunc = Core::lexicalCast(Core::trimmed(list[1])); + shellTypes.push_back(orbitalIdx(Core::trimmed(list[0]))); + shellFunctions.push_back(nFunc); + m_basisFunctions.at(nGTOs)->push_back( + new std::vector(nFunc)); + + for (int i = 0; i < nFunc; i++) { + getline(in, key); + key = Core::trimmed(key); + + list = Core::split(key, ' '); + m_basisFunctions.at(nGTOs)->at(nShells)->at(i).x() = + Core::lexicalCast(list[1]); // exponent + m_basisFunctions.at(nGTOs)->at(nShells)->at(i).y() = + Core::lexicalCast(list[2]); // coeff + } + + nShells++; + getline(in, key); + key = Core::trimmed(key); + + list = Core::split(key, ' '); + } + 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; + nGTOs++; + + getline(in, key); + getline(in, key); + getline(in, key); + key = Core::trimmed(key); + + list = Core::split(key, ' '); + if (list.size() == 0) + break; // unexpected structure - suppose no more NewGTOs + } + + // create input for gaussian basisset + int nBasis = nGTOs; + int nAtoms = m_atomLabel.size(); + m_currentAtom = 0; + for (int i = 0; i < nAtoms; i++) { + m_currentAtom++; + for (int j = 0; j < nBasis; j++) { + if (m_atomLabel.at(i) == m_basisAtomLabel.at(j)) { + for (unsigned int k = 0; k < m_orcaNumShells.at(j).size(); k++) { + for (int l = 0; l < m_orcaNumShells.at(j).at(k); l++) { + m_a.push_back(m_basisFunctions.at(j)->at(k)->at(l).x()); + m_c.push_back(m_basisFunctions.at(j)->at(k)->at(l).y()); + } + m_shellNums.push_back(m_orcaNumShells.at(j).at(k)); + m_shellTypes.push_back(m_orcaShellTypes.at(j).at(k)); + m_shelltoAtom.push_back(m_currentAtom); + } + break; + } + } + } + m_currentMode = NotParsing; + break; + } + case MO: { + + m_MOcoeffs.clear(); // if the orbitals were punched multiple times + std::vector orcaOrbitals; + + while (!Core::trimmed(key).empty()) { + // currently reading the sequence number + getline(in, key); // energies + getline(in, key); // symmetries + getline(in, key); // skip ----------- + getline(in, key); // now we've got coefficients + + regex rx("[.][0-9]{6}[0-9-]"); + auto key_begin = std::sregex_iterator(key.begin(), key.end(), rx); + auto key_end = std::sregex_iterator(); + for (std::sregex_iterator i = key_begin; i != key_end; ++i) { + key += i->str() + " "; + } + + list = Core::split(key, ' '); + + numColumns = list.size() - 2; + columns.resize(numColumns); + while (list.size() > 2) { + orcaOrbitals.push_back(list[1]); + for (unsigned int i = 0; i < numColumns; ++i) { + columns[i].push_back(Core::lexicalCast(list[i + 2])); + } + + getline(in, key); + key_begin = std::sregex_iterator(key.begin(), key.end(), rx); + key_end = std::sregex_iterator(); + for (std::sregex_iterator i = key_begin; i != key_end; ++i) { + key += i->str() + " "; + } + + list = Core::split(key, ' '); + if (list.size() != numColumns + 2) + break; + + } // ok, we've finished one batch of MO coeffs + // now reorder the p orbitals from "orcaStyle" (pz, px,py) + // to expected Avogadro (px,py,pz) + int idx = 0; + while (idx < orcaOrbitals.size()) { + if (Core::contains(orcaOrbitals.at(idx), "pz")) { + for (unsigned int i = 0; i < numColumns; i++) { + std::swap(columns[i].at(idx), columns[i].at(idx + 1)); + } + idx++; + for (unsigned int i = 0; i < numColumns; i++) { + std::swap(columns[i].at(idx), columns[i].at(idx + 1)); + } + idx++; + idx++; + } else { + idx++; + } + } + + // Now we need to re-order the MO coeffs, so we insert one MO at a + // time + for (unsigned int i = 0; i < numColumns; ++i) { + numRows = columns[i].size(); + for (unsigned int j = 0; j < numRows; ++j) { + m_MOcoeffs.push_back(columns[i][j]); + } + } + columns.clear(); + orcaOrbitals.clear(); + + } // finished parsing MOs + if (m_MOcoeffs.size() != numRows * numRows) { + m_orcaSuccess = false; + } + m_numBasisFunctions = numRows; + if (m_openShell) { + // TODO: parse both alpha and beta orbitals + + m_BetaMOcoeffs.clear(); // if the orbitals were punched multiple times + getline(in, key); + while (!Core::trimmed(key).empty()) { + // currently reading the sequence number + getline(in, key); // energies + getline(in, key); // symmetries + getline(in, key); // skip ----------- + getline(in, key); // now we've got coefficients + + regex rx("[.][0-9]{6}[0-9-]"); + auto key_begin = std::sregex_iterator(key.begin(), key.end(), rx); + auto key_end = std::sregex_iterator(); + for (std::sregex_iterator i = key_begin; i != key_end; ++i) { + key += i->str() + " "; + } + + list = Core::split(key, ' '); + numColumns = list.size() - 2; + columns.resize(numColumns); + while (list.size() > 2) { + orcaOrbitals.push_back(list[1]); + // columns.resize(numColumns); + for (unsigned int i = 0; i < numColumns; ++i) { + columns[i].push_back(Core::lexicalCast(list[i + 2])); + } + + getline(in, key); + key_begin = std::sregex_iterator(key.begin(), key.end(), rx); + key_end = std::sregex_iterator(); + for (std::sregex_iterator i = key_begin; i != key_end; ++i) { + key += i->str() + " "; + } + list = Core::split(key, ' '); + if (list.size() != numColumns + 2) + break; + + } // ok, we've finished one batch of MO coeffs + // now reorder the p orbitals from "orcaStyle" (pz, px,py) to + // expected Avogadro (px,py,pz) + int idx = 0; + while (idx < orcaOrbitals.size()) { + if (Core::contains(orcaOrbitals.at(idx), "pz")) { + for (unsigned int i = 0; i < numColumns; i++) { + std::swap(columns[i].at(idx), columns[i].at(idx + 1)); + } + idx++; + for (unsigned int i = 0; i < numColumns; i++) { + std::swap(columns[i].at(idx), columns[i].at(idx + 1)); + } + idx++; + idx++; + } else { + idx++; + } + } + + // Now we need to re-order the MO coeffs, so we insert one MO at a + // time + for (unsigned int i = 0; i < numColumns; ++i) { + numRows = columns[i].size(); + for (unsigned int j = 0; j < numRows; ++j) { + + m_BetaMOcoeffs.push_back(columns[i][j]); + } + } + columns.clear(); + orcaOrbitals.clear(); + + if (Core::trimmed(key).empty()) + getline(in, key); // skip the blank line after the MOs + } // finished parsing 2nd. MOs + if (m_MOcoeffs.size() != numRows * numRows) { + m_orcaSuccess = false; + } + m_numBasisFunctions = numRows; + } + + m_currentMode = NotParsing; + break; + } + default:; + } // end switch + } // end if (mode) +} + +void ORCAOutput::load(GaussianSet* basis) +{ + // Now load up our basis set + basis->setElectronCount(m_electrons); + + // Set up the GTO primitive counter, go through the shells and add them + int nGTO = 0; + 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) == GaussianSet::SP) { + // SP orbital type - currently have to unroll into two shells + int tmpGTO = nGTO; + 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, GaussianSet::P); + for (int j = 0; j < m_shellNums.at(i); ++j) { + basis->addGto(p, m_csp.at(nSP), m_a.at(tmpGTO)); + ++tmpGTO; + ++nSP; + } + } else { + int b = basis->addBasis(m_shelltoAtom.at(i) - 1, m_shellTypes.at(i)); + for (int j = 0; j < m_shellNums.at(i); ++j) { + basis->addGto(b, m_c.at(nGTO), m_a.at(nGTO)); + ++nGTO; + } + } + } + + // Now to load in the MO coefficients + if (m_MOcoeffs.size()) + basis->setMolecularOrbitals(m_MOcoeffs); + if (m_BetaMOcoeffs.size()) + basis->setMolecularOrbitals(m_BetaMOcoeffs, Core::BasisSet::Beta); + + if (m_orbitalEnergy.size()) + basis->setMolecularOrbitalEnergy(m_orbitalEnergy); + if (m_betaOrbitalEnergy.size()) + basis->setMolecularOrbitalEnergy(m_betaOrbitalEnergy, Core::BasisSet::Beta); + + // TODO: set orbital symmetries + + m_homo = ceil(m_electrons / 2.0); + basis->generateDensityMatrix(); +} + +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 new file mode 100644 index 0000000000..5334d1de33 --- /dev/null +++ b/avogadro/quantumio/orca.h @@ -0,0 +1,124 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#ifndef AVOGADRO_QUANTUMIO_ORCA_H +#define AVOGADRO_QUANTUMIO_ORCA_H + +#include "avogadroquantumioexport.h" +#include +#include +#include + +#include +#include + +namespace Avogadro { +namespace QuantumIO { + +class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat +{ +public: + ORCAOutput(); + ~ORCAOutput() override; + + Operations supportedOperations() const override + { + return Read | File | Stream | String; + } + + FileFormat* newInstance() const override { return new ORCAOutput; } + std::string identifier() const override { return "Avogadro: Orca"; } + std::string name() const override { return "Orca"; } + std::string description() const override { return "Orca output format."; } + + std::string specificationUrl() const override { return ""; } + + std::vector fileExtensions() const override; + std::vector mimeTypes() const override; + + bool read(std::istream& in, Core::Molecule& molecule) override; + bool write(std::ostream&, const Core::Molecule&) override + { + // Empty, as we do not write out Orca output files. + return false; + } + +private: + void outputAll(); + + void processLine(std::istream& in, Core::GaussianSet* basis); + void load(Core::GaussianSet* basis); + + // OrcaStuff + void orcaWarningMessage(const std::string& m); + Core::GaussianSet::orbital orbitalIdx(std::string txt); + bool m_orcaSuccess; + + std::vector m_atomLabel; + std::vector m_basisAtomLabel; + + std::vector m_atomNums; + std::vector m_atomPos; + + std::vector shellFunctions; + std::vector shellTypes; + std::vector> m_orcaNumShells; + std::vector> m_orcaShellTypes; + int m_nGroups; + + std::vector*>*> m_basisFunctions; + + enum mode + { + Atoms, + GTO, + MO, + OrbitalEnergies, + Charges, + Frequencies, + VibrationalModes, + IR, + Raman, + Electronic, + NMR, + NotParsing, + Unrecognized + }; + + double m_coordFactor; + mode m_currentMode; + int m_electrons; + + bool m_openShell; + bool m_readBeta; + + int m_homo; + + int m_currentAtom; + unsigned int m_numBasisFunctions; + std::vector m_shellTypes; + std::vector m_shellNums; + std::vector m_shelltoAtom; + std::vector m_a; + std::vector m_c; + std::vector m_csp; + std::vector m_orbitalEnergy; + std::vector m_MOcoeffs; + std::vector m_betaOrbitalEnergy; + std::vector m_BetaMOcoeffs; + + std::string m_chargeType; + std::map m_partialCharges; + + Core::Array m_frequencies; + Core::Array m_IRintensities; + Core::Array m_RamanIntensities; + Core::Array> m_vibDisplacements; +}; + +} // namespace QuantumIO +} // namespace Avogadro + +#endif