Skip to content

Commit

Permalink
Merge pull request #1380 from ghutchis/reading-fchk-vibrations
Browse files Browse the repository at this point in the history 
Reading Gaussian fchk vibrations when present
  • Loading branch information
@ghutchis
ghutchis authored Oct 19, 2023
2 parents 545f9e4 + 4f0291a commit 62eecfc
Show file tree
Hide file tree
Showing 2 changed files with 73 additions and 26 deletions.
85 changes: 61 additions & 24 deletions avogadro/quantumio/gaussianfchk.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -84,14 +95,15 @@ void GaussianFchk::processLine(std::istream& in)

string tmp = line.substr(43);
vector<string> list = Core::split(tmp, ' ');
std::vector<double> tmpVec;

// Big switch statement checking for various things we are interested in
if (Core::contains(key, "RHF")) {
m_scftype = Rhf;
} else if (Core::contains(key, "UHF")) {
m_scftype = Uhf;
} else if (key == "Number of atoms" && list.size() > 1) {
cout << "Number of atoms = " << Core::lexicalCast<int>(list[1]) << endl;
m_numAtoms = Core::lexicalCast<int>(list[1]);
} else if (key == "Charge" && list.size() > 1) {
m_charge = Core::lexicalCast<signed char>(list[1]);
} else if (key == "Multiplicity" && list.size() > 1) {
Expand All @@ -104,13 +116,11 @@ void GaussianFchk::processLine(std::istream& in)
m_electronsBeta = Core::lexicalCast<int>(list[1]);
} else if (key == "Number of basis functions" && list.size() > 1) {
m_numBasisFunctions = Core::lexicalCast<int>(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<int>(list[2]));
if (static_cast<int>(m_aNums.size()) != Core::lexicalCast<int>(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) {
Expand All @@ -134,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<int>(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<int>(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<double>();
Expand All @@ -152,37 +163,63 @@ void GaussianFchk::processLine(std::istream& in)
}

m_betaOrbitalEnergy = readArrayD(in, Core::lexicalCast<int>(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<int>(list[2]), 16);
if (static_cast<int>(m_MOcoeffs.size()) ==
Core::lexicalCast<int>(list[2]))
cout << "MO coefficients, n = " << m_MOcoeffs.size() << endl;
} else if (m_scftype == Uhf) {
m_alphaMOcoeffs = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
if (static_cast<int>(m_alphaMOcoeffs.size()) ==
Core::lexicalCast<int>(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<int>(list[2]), 16);
if (static_cast<int>(m_betaMOcoeffs.size()) ==
Core::lexicalCast<int>(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<int>(list[2]), 16))
cout << "SCF density matrix read in " << m_density.rows() << endl;
else
if (!readDensityMatrix(in, Core::lexicalCast<int>(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<int>(list[2]), 16))
cout << "SCF spin density matrix read in " << m_spinDensity.rows()
<< endl;
else
if (!readSpinDensityMatrix(in, Core::lexicalCast<int>(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<int>(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
tmpVec = readArrayD(in, Core::lexicalCast<int>(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(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(tmpVec[i]);
}
// now check if we have Raman intensities
if (tmpVec[threeN + m_normalModes] != 0.0) {
for (unsigned int i = threeN + m_normalModes;
i < threeN + 2 * m_normalModes; ++i) {
m_RamanIntensities.push_back(tmpVec[i]);
}
}
} else if (key == "Vib-Modes" && list.size() > 2) {
tmpVec = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
m_vibDisplacements.clear();
if (tmpVec.size() == m_numAtoms * 3 * m_normalModes) {
for (unsigned int i = 0; i < m_normalModes; ++i) {
Core::Array<Vector3> mode;
for (unsigned int j = 0; j < m_numAtoms; ++j) {
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);
}
}
}
}

Expand Down
14 changes: 12 additions & 2 deletions avogadro/quantumio/gaussianfchk.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
#define AVOGADRO_QUANTUMIO_GAUSSIANFCHK_H

#include "avogadroquantumioexport.h"

#include <avogadro/core/array.h>
#include <avogadro/core/gaussianset.h>
#include <avogadro/io/fileformat.h>

Expand Down Expand Up @@ -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;
Expand All @@ -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<double> m_frequencies;
Core::Array<double> m_IRintensities;
Core::Array<double> m_RamanIntensities;
Core::Array<Core::Array<Vector3>> m_vibDisplacements;
};
}
}

} // namespace QuantumIO
} // namespace Avogadro

#endif

0 comments on commit 62eecfc

Please sign in to comment.