Skip to content

Commit

Permalink
Merge pull request #1797 from ghutchis/orca-electronic-spectra
Browse files Browse the repository at this point in the history
Parse Orca electronic spectra
  • Loading branch information
ghutchis authored Nov 18, 2024
2 parents 5c130f5 + 6c3a053 commit 519658f
Show file tree
Hide file tree
Showing 4 changed files with 131 additions and 0 deletions.
12 changes: 12 additions & 0 deletions avogadro/quantumio/gaussianfchk.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,13 @@ bool GaussianFchk::read(std::istream& in, Core::Molecule& molecule)
molecule.setVibrationRamanIntensities(m_RamanIntensities);
}

// set the total charge
molecule.setData("totalCharge", m_charge);
// set the spin multiplicity
molecule.setData("totalSpinMultiplicity", m_spin);
// dipole moment
molecule.setData("dipoleMoment", m_dipoleMoment);

// Do simple bond perception.
molecule.perceiveBondsSimple();
molecule.perceiveBondOrders();
Expand Down Expand Up @@ -108,6 +115,11 @@ void GaussianFchk::processLine(std::istream& in)
m_charge = Core::lexicalCast<signed char>(list[1]);
} else if (key == "Multiplicity" && list.size() > 1) {
m_spin = Core::lexicalCast<char>(list[1]);
} else if (key == "Dipole Moment" && list.size() > 2) {
vector<double> dipole = readArrayD(in, Core::lexicalCast<int>(list[2]));
m_dipoleMoment = Vector3(dipole[0], dipole[1], dipole[2]);
// convert from au
m_dipoleMoment *= 2.541746;
} else if (key == "Number of electrons" && list.size() > 1) {
m_electrons = Core::lexicalCast<int>(list[1]);
} else if (key == "Number of alpha electrons" && list.size() > 1) {
Expand Down
1 change: 1 addition & 0 deletions avogadro/quantumio/gaussianfchk.h
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ class AVOGADROQUANTUMIO_EXPORT GaussianFchk : public Io::FileFormat
std::vector<double> m_betaMOcoeffs;
MatrixX m_density; /// Total density matrix
MatrixX m_spinDensity; /// Spin density matrix
Vector3 m_dipoleMoment;
Core::ScfType m_scftype;

Core::Array<double> m_frequencies;
Expand Down
109 changes: 109 additions & 0 deletions avogadro/quantumio/orca.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,27 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
molecule.setVibrationRamanIntensities(m_RamanIntensities);
}

if (m_electronicTransitions.size() > 0 &&
m_electronicTransitions.size() == m_electronicIntensities.size()) {
MatrixX electronicData(m_electronicTransitions.size(), 2);
for (size_t i = 0; i < m_electronicTransitions.size(); ++i) {
electronicData(i, 0) = m_electronicTransitions[i];
electronicData(i, 1) = m_electronicIntensities[i];
}
molecule.setSpectra("Electronic", electronicData);
std::cout << "UV/Vis data found." << electronicData.rows() << std::endl;
if (m_electronicRotations.size() == m_electronicTransitions.size()) {
MatrixX electronicRotations(m_electronicTransitions.size(), 2);
for (size_t i = 0; i < m_electronicTransitions.size(); ++i) {
electronicRotations(i, 0) = m_electronicTransitions[i];
electronicRotations(i, 1) = m_electronicRotations[i];
}
molecule.setSpectra("CircularDichroism", electronicRotations);
std::cout << "CircularDichroism data found." << electronicRotations.rows()
<< std::endl;
}
}

// guess bonds and bond orders
molecule.perceiveBondsSimple();
molecule.perceiveBondOrders();
Expand All @@ -102,6 +123,10 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
basis->setMolecule(&molecule);
load(basis);

molecule.setData("totalCharge", m_charge);
molecule.setData("totalSpinMultiplicity", m_spin);
molecule.setData("dipoleMoment", m_dipoleMoment);

return true;
}

Expand Down Expand Up @@ -157,16 +182,45 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
getline(in, key);
}
}
} else if (Core::contains(key, "Total Charge")) {
list = Core::split(key, ' ');
if (list.size() > 4)
m_charge = Core::lexicalCast<int>(list[4]);
} else if (Core::contains(key, "Multiplicity")) {
list = Core::split(key, ' ');
if (list.size() > 3)
m_spin = Core::lexicalCast<int>(list[3]);
} 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<int>(list[5]);
} else if (Core::contains(key, "Total Dipole Moment")) {
list = Core::split(key, ' ');
m_dipoleMoment = Eigen::Vector3d(Core::lexicalCast<double>(list[4]),
Core::lexicalCast<double>(list[5]),
Core::lexicalCast<double>(list[6]));
// convert from atomic units to Debye
// e.g. https://en.wikipedia.org/wiki/Debye
m_dipoleMoment *= 2.54174628;
} else if (Core::contains(key, "Mayer bond orders")) {
m_currentMode = BondOrders;
// starts at the next line
} else if (Core::contains(
key,
"ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS")) {
m_currentMode = Electronic;
for (int i = 0; i < 4; ++i) {
getline(in, key); // skip header
}
// starts at the next line
} else if (Core::contains(key, "CD SPECTRUM")) {
m_currentMode = ECD;
for (int i = 0; i < 4; ++i) {
getline(in, key); // skip header
}
} else if (Core::contains(key, "ORBITAL ENERGIES")) {
m_currentMode = OrbitalEnergies;
getline(in, key); // skip ------------
Expand Down Expand Up @@ -509,6 +563,61 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
m_currentMode = NotParsing;
break;
}
case Electronic: {
if (key.empty())
break;
list = Core::split(key, ' ');
double wavenumbers;
while (!key.empty()) {
// should have 8 columns
if (list.size() != 8) {
getline(in, key);
key = Core::trimmed(key);
list = Core::split(key, ' ');
continue; // skip any spin-forbidden transitions
}

wavenumbers = Core::lexicalCast<double>(list[1]);
// convert to eV
m_electronicTransitions.push_back(wavenumbers / 8065.544);
m_electronicIntensities.push_back(Core::lexicalCast<double>(list[3]));

getline(in, key);
key = Core::trimmed(key);
list = Core::split(key, ' ');
if (list.size() < 2)
break; // hit the blank line
}
m_currentMode = NotParsing;
}
case ECD: {
if (key.empty())
break;
list = Core::split(key, ' ');

double wavenumbers;
while (!key.empty()) {
// should have 7 columns
if (list.size() != 7) {
getline(in, key);
key = Core::trimmed(key);
list = Core::split(key, ' ');
continue; // skip any spin-forbidden transitions
}

wavenumbers = Core::lexicalCast<double>(list[1]);
// convert to eV
m_electronicTransitions.push_back(wavenumbers / 8065.544);
m_electronicRotations.push_back(Core::lexicalCast<double>(list[3]));

getline(in, key);
key = Core::trimmed(key);
list = Core::split(key, ' ');
if (list.size() < 2)
break; // hit the blank line
}
m_currentMode = NotParsing;
}
case GTO: {
// // should start at the first newGTO
if (key.empty())
Expand Down
9 changes: 9 additions & 0 deletions avogadro/quantumio/orca.h
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
std::vector<int> m_atomNums;
std::vector<Eigen::Vector3d> m_atomPos;

Eigen::Vector3d m_dipoleMoment;

std::vector<std::vector<int>> m_bondOrders;

std::vector<int> shellFunctions;
Expand All @@ -84,6 +86,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
IR,
Raman,
Electronic,
ECD, // electronic circular dichroism
NMR,
BondOrders,
NotParsing,
Expand All @@ -98,6 +101,8 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
bool m_readBeta;

int m_homo;
int m_charge;
int m_spin;

int m_currentAtom;
unsigned int m_numBasisFunctions;
Expand All @@ -119,6 +124,10 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
Core::Array<double> m_IRintensities;
Core::Array<double> m_RamanIntensities;
Core::Array<Core::Array<Vector3>> m_vibDisplacements;

Core::Array<double> m_electronicTransitions; // in eV
Core::Array<double> m_electronicIntensities;
Core::Array<double> m_electronicRotations; // for CD
};

} // namespace QuantumIO
Expand Down

0 comments on commit 519658f

Please sign in to comment.