Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

chemicaljson enhancements #1321

Merged
merged 4 commits into from
Jul 26, 2023
Merged

chemicaljson enhancements #1321

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
84c8ed3
Add support for properties, open-shell orbital energies and occ.
ghutchis Jul 26, 2023
49d8293
Add support for orbital energies for fchk and molden
ghutchis Jul 26, 2023
d6d488b
Fix formatting
ghutchis Jul 26, 2023
c4fde81
Fix formatting
ghutchis Jul 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
148 changes: 142 additions & 6 deletions avogadro/io/cjsonformat.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -388,6 +388,7 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
json moCoefficients = orbitals["moCoefficients"];
json moCoefficientsA = orbitals["alphaCoefficients"];
json moCoefficientsB = orbitals["betaCoefficients"];
bool openShell = false;
if (isNumericArray(moCoefficients)) {
std::vector<double> coeffs;
for (auto& moCoefficient : moCoefficients)
Expand All @@ -403,6 +404,7 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
coeffsB.push_back(static_cast<double>(i));
basis->setMolecularOrbitals(coeffsA, BasisSet::Alpha);
basis->setMolecularOrbitals(coeffsB, BasisSet::Beta);
openShell = true;
} else {
std::cout << "No orbital cofficients found!" << std::endl;
}
Expand All @@ -429,11 +431,44 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
coeffsB.push_back(static_cast<double>(i));
basis->setMolecularOrbitals(coeffsA, BasisSet::Alpha, idx);
basis->setMolecularOrbitals(coeffsB, BasisSet::Beta, idx);
openShell = true;
}
}
// Set the first step as active.
basis->setActiveSetStep(0);
}
if (openShell) {
// look for alpha and beta orbital energies
json energiesA = orbitals["alphaEnergies"];
json energiesB = orbitals["betaEnergies"];
// check if they are numeric arrays
if (isNumericArray(energiesA) && isNumericArray(energiesB)) {
std::vector<double> moEnergiesA;
for (auto& i : energiesA)
moEnergiesA.push_back(static_cast<double>(i));
std::vector<double> moEnergiesB;
for (auto& i : energiesB)
moEnergiesB.push_back(static_cast<double>(i));
basis->setMolecularOrbitalEnergy(moEnergiesA, BasisSet::Alpha);
basis->setMolecularOrbitalEnergy(moEnergiesB, BasisSet::Beta);

// look for alpha and beta orbital occupations
json occupationsA = orbitals["alphaOccupations"];
json occupationsB = orbitals["betaOccupations"];
// check if they are numeric arrays
if (isNumericArray(occupationsA) && isNumericArray(occupationsB)) {
std::vector<unsigned char> moOccupationsA;
for (auto& i : occupationsA)
moOccupationsA.push_back(static_cast<unsigned char>(i));
std::vector<unsigned char> moOccupationsB;
for (auto& i : occupationsB)
moOccupationsB.push_back(static_cast<unsigned char>(i));
basis->setMolecularOrbitalOccupancy(moOccupationsA,
BasisSet::Alpha);
basis->setMolecularOrbitalOccupancy(moOccupationsB, BasisSet::Beta);
}
}
}
}
molecule.setBasisSet(basis);
}
Expand Down Expand Up @@ -495,7 +530,59 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
molecule.setData("totalSpinMultiplicity",
static_cast<int>(properties["totalSpinMultiplicity"]));
}
// todo - put everything into the data map
// iterate through everything else
for (auto& element : properties.items()) {
if (element.key() == "totalCharge" ||
element.key() == "totalSpinMultiplicity") {
continue;
}
if (element.value().type() == json::value_t::array) {
// check if it is a numeric array to go into Eigen::MatrixXd
json j = element.value(); // convenience
std::size_t rows = j.size();
MatrixX matrix;
matrix.resize(rows, 1); // default to 1 columns
bool isNumeric = true;

for (std::size_t row = 0; row < j.size(); ++row) {
const auto& jrow = j.at(row);
// check to see if we have a simple vector or a matrix
if (jrow.type() == json::value_t::array) {
matrix.conservativeResize(rows, jrow.size());
for (std::size_t col = 0; col < jrow.size(); ++col) {
const auto& value = jrow.at(col);
if (value.type() == json::value_t::number_float ||
value.type() == json::value_t::number_integer ||
value.type() == json::value_t::number_unsigned)
matrix(row, col) = value.get<double>();
else {
isNumeric = false;
break;
}
}
} else if (jrow.type() == json::value_t::number_float ||
jrow.type() == json::value_t::number_integer ||
jrow.type() == json::value_t::number_unsigned) {
// just a row vector
matrix(row, 0) = jrow.get<double>();
} else {
isNumeric = false;
break;
}
}
if (isNumeric)
molecule.setData(element.key(), matrix);
// TODO: add support for non-numeric arrays
// std::cout << " property: " << element.key() << " = " << matrix
// << " size " << matrix.rows() << 'x' << matrix.cols()
// << std::endl;
} else {
molecule.setData(element.key(), element.value());
// std::cout << " property: " << element.key() << " = "
// << element.value() << " type "
// << element.value().type_name() << std::endl;
}
}
}
}

Expand Down Expand Up @@ -558,6 +645,35 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
// or approximate from formal charges and # of electrons
properties["totalCharge"] = molecule.totalCharge();
properties["totalSpinMultiplicity"] = molecule.totalSpinMultiplicity();
// loop through all other properties
const auto map = molecule.dataMap();
for (const auto& element : map) {
if (element.first == "name" || element.first == "inchi")
continue;

if (element.second.type() == Variant::String)
properties[element.first] = element.second.toString().c_str();
else if (element.second.type() == Variant::Double)
properties[element.first] = element.second.toDouble();
else if (element.second.type() == Variant::Float)
properties[element.first] = element.second.toFloat();
else if (element.second.type() == Variant::Int)
properties[element.first] = element.second.toInt();
else if (element.second.type() == Variant::Bool)
properties[element.first] = element.second.toBool();
else if (element.second.type() == Variant::Matrix) {
MatrixX m = element.second.toMatrix();
json matrix;
for (int i = 0; i < m.rows(); ++i) {
json row;
for (int j = 0; j < m.cols(); ++j) {
row.push_back(m(i, j));
}
matrix.push_back(row);
}
properties[element.first] = matrix;
}
}
root["properties"] = properties;

if (molecule.unitCell()) {
Expand Down Expand Up @@ -616,7 +732,8 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
}
basis["shellTypes"] = shellTypes;

// This bit is slightly tricky, map from our index to primitives per shell.
// This bit is slightly tricky, map from our index to primitives per
// shell.
if (gaussian->gtoIndices().size() && gaussian->atomIndices().size()) {
auto gtoIndices = gaussian->gtoIndices();
auto gtoA = gaussian->gtoA();
Expand Down Expand Up @@ -646,8 +763,8 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
root["basisSet"] = basis;
}

// Now get the MO matrix, potentially other things. Need to get a handle on
// when we have just one (paired), or two (alpha and beta) to write.
// Now get the MO matrix, potentially other things. Need to get a handle
// on when we have just one (paired), or two (alpha and beta) to write.
auto moMatrix = gaussian->moMatrix();
auto betaMatrix = gaussian->moMatrix(BasisSet::Beta);
json moCoefficients;
Expand All @@ -674,14 +791,33 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
for (double& energie : energies) {
energyData.push_back(energie);
}
root["orbitals"]["energies"] = energyData;

auto betaEnergies = gaussian->moEnergy(BasisSet::Beta);
if (betaEnergies.size() > 0) {
json betaEnergyData;
for (double& energie : betaEnergies) {
betaEnergyData.push_back(energie);
}
root["orbitals"]["alphaEnergies"] = energyData;
root["orbitals"]["betaEnergies"] = betaEnergyData;
} else
root["orbitals"]["energies"] = energyData;
}
auto occ = gaussian->moOccupancy();
if (occ.size() > 0) {
json occData;
for (unsigned char& it : occ)
occData.push_back(static_cast<int>(it));
root["orbitals"]["occupations"] = occData;

auto betaOcc = gaussian->moOccupancy(BasisSet::Beta);
if (betaOcc.size() > 0) {
json betaOccData;
for (unsigned char& it : betaOcc)
betaOccData.push_back(static_cast<int>(it));
root["orbitals"]["alphaOccupations"] = occData;
root["orbitals"]["betaOccupations"] = betaOccData;
} else
root["orbitals"]["occupations"] = occData;
}
auto num = gaussian->moNumber();
if (num.size() > 0) {
Expand Down
28 changes: 14 additions & 14 deletions avogadro/quantumio/gaussianfchk.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,27 +11,23 @@

#include <iostream>

using std::vector;
using std::string;
using std::cout;
using std::endl;
using std::string;
using std::vector;

namespace Avogadro::QuantumIO {

using Core::Atom;
using Core::BasisSet;
using Core::GaussianSet;
using Core::Rhf;
using Core::Uhf;
using Core::Rohf;
using Core::Uhf;

GaussianFchk::GaussianFchk() : m_scftype(Rhf)
{
}
GaussianFchk::GaussianFchk() : m_scftype(Rhf) {}

GaussianFchk::~GaussianFchk()
{
}
GaussianFchk::~GaussianFchk() {}

std::vector<std::string> GaussianFchk::fileExtensions() const
{
Expand Down Expand Up @@ -282,6 +278,10 @@ void GaussianFchk::load(GaussianSet* basis)
basis->setDensityMatrix(m_density);
if (m_spinDensity.rows())
basis->setSpinDensityMatrix(m_spinDensity);
if (m_alphaOrbitalEnergy.size())
basis->setMolecularOrbitalEnergy(m_alphaOrbitalEnergy, BasisSet::Alpha);
if (m_betaOrbitalEnergy.size())
basis->setMolecularOrbitalEnergy(m_betaOrbitalEnergy, BasisSet::Beta);
} else {
cout << "Basis set is not valid!\n";
}
Expand All @@ -303,7 +303,7 @@ vector<int> GaussianFchk::readArrayI(std::istream& in, unsigned int n)
return tmp;

vector<string> list = Core::split(line, ' ');
for (auto & i : list) {
for (auto& i : list) {
if (tmp.size() >= n) {
cout << "Too many variables read in. File may be inconsistent. "
<< tmp.size() << " of " << n << endl;
Expand Down Expand Up @@ -338,7 +338,7 @@ vector<double> GaussianFchk::readArrayD(std::istream& in, unsigned int n,

if (width == 0) { // we can split by spaces
vector<string> list = Core::split(line, ' ');
for (auto & i : list) {
for (auto& i : list) {
if (tmp.size() >= n) {
cout << "Too many variables read in. File may be inconsistent. "
<< tmp.size() << " of " << n << endl;
Expand Down Expand Up @@ -395,7 +395,7 @@ bool GaussianFchk::readDensityMatrix(std::istream& in, unsigned int n,

if (width == 0) { // we can split by spaces
vector<string> list = Core::split(line, ' ');
for (auto & k : list) {
for (auto& k : list) {
if (cnt >= n) {
cout << "Too many variables read in. File may be inconsistent. "
<< cnt << " of " << n << endl;
Expand Down Expand Up @@ -471,7 +471,7 @@ bool GaussianFchk::readSpinDensityMatrix(std::istream& in, unsigned int n,

if (width == 0) { // we can split by spaces
vector<string> list = Core::split(line, ' ');
for (auto & k : list) {
for (auto& k : list) {
if (cnt >= n) {
cout << "Too many variables read in. File may be inconsistent. "
<< cnt << " of " << n << endl;
Expand Down Expand Up @@ -566,4 +566,4 @@ void GaussianFchk::outputAll()
cout << endl << endl;
}
}
}
} // namespace Avogadro::QuantumIO
15 changes: 9 additions & 6 deletions avogadro/quantumio/molden.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,10 +11,10 @@

#include <iostream>

using std::vector;
using std::string;
using std::cout;
using std::endl;
using std::string;
using std::vector;

namespace Avogadro::QuantumIO {

Expand All @@ -26,9 +26,7 @@ MoldenFile::MoldenFile()
{
}

MoldenFile::~MoldenFile()
{
}
MoldenFile::~MoldenFile() {}

std::vector<std::string> MoldenFile::fileExtensions() const
{
Expand Down Expand Up @@ -161,6 +159,9 @@ void MoldenFile::processLine(std::istream& in)
list = Core::split(line, ' ');
if (Core::contains(line, "Occup"))
m_electrons += Core::lexicalCast<int>(list[1]);
else if (Core::contains(line, "Ene"))
m_orbitalEnergy.push_back(Core::lexicalCast<double>(list[1]));
// TODO: track alpha beta spin
}

// Parse the molecular orbital coefficients.
Expand Down Expand Up @@ -224,6 +225,8 @@ void MoldenFile::load(GaussianSet* basis)
// Now to load in the MO coefficients
if (m_MOcoeffs.size())
basis->setMolecularOrbitals(m_MOcoeffs);
if (m_orbitalEnergy.size())
basis->setMolecularOrbitalEnergy(m_orbitalEnergy);
}

void MoldenFile::outputAll()
Expand All @@ -238,4 +241,4 @@ void MoldenFile::outputAll()
cout << m_MOcoeff << "\t";
cout << endl;
}
}
} // namespace Avogadro::QuantumIO
Loading