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.

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

Initial support for total charge and spin through CJSON. #1058

Merged
merged 1 commit into from
Jul 27, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
17 changes: 17 additions & 0 deletions avogadro/core/molecule.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -321,6 +321,23 @@ const Array<signed char>& Molecule::formalCharges() const
return m_formalCharges;
}

signed char Molecule::totalCharge() const
{
signed char charge = 0;

// check the data map first
if (m_data.hasValue("totalCharge")) {
charge = m_data.value("totalCharge").toInt();
}

if (m_formalCharges.size() > 0) {
for (Index i = 0; i < m_formalCharges.size(); ++i)
charge += m_formalCharges[i];
return charge;
}
return charge; // should be zero
}

Array<Vector3ub>& Molecule::colors()
{
return m_colors;
Expand Down
11 changes: 10 additions & 1 deletion avogadro/core/molecule.h
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,15 @@ class AVOGADROCORE_EXPORT Molecule
/** \overload */
const Array<signed char>& formalCharges() const;

/**
* Get the total charge on the molecule.
* The method will first check to see if a total charge has been set. If not,
* it will calculate the total charge from the formal charges (if set).
* If neither has been set, it will assume the total charge is zero.
* @return The total charge of the molecule.
*/
signed char totalCharge() const;

/**
* Get the formal charge for the requested atom.
* @param atomId The index of the atom.
Expand All @@ -159,7 +168,7 @@ class AVOGADROCORE_EXPORT Molecule
*/
bool setFormalCharge(Index atomId, signed char charge);

/** Returns a vector of colors for the atoms in the moleucle. */
/** \returns a vector of colors for the atoms in the moleucle. */
Array<Vector3ub>& colors();

/** \overload */
Expand Down
69 changes: 47 additions & 22 deletions avogadro/io/cjsonformat.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
// This represents our minimal spec for a molecule - atoms that have an
// atomic number.
if (isNumericArray(atomicNumbers) && atomicNumbers.size() > 0) {
for (auto & atomicNumber : atomicNumbers)
for (auto& atomicNumber : atomicNumbers)
molecule.addAtom(atomicNumber);
} else {
appendError("Malformed array for in atoms.elements.number");
Expand Down Expand Up @@ -367,21 +367,21 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
json occupations = orbitals["occupations"];
if (isNumericArray(occupations)) {
std::vector<unsigned char> occs;
for (auto & occupation : occupations)
for (auto& occupation : occupations)
occs.push_back(static_cast<unsigned char>(occupation));
basis->setMolecularOrbitalOccupancy(occupations);
}
json energies = orbitals["energies"];
if (isNumericArray(energies)) {
std::vector<double> energyArray;
for (auto & energie : energies)
for (auto& energie : energies)
energyArray.push_back(static_cast<double>(energie));
basis->setMolecularOrbitalEnergy(energyArray);
}
json numbers = orbitals["numbers"];
if (isNumericArray(numbers)) {
std::vector<unsigned int> numArray;
for (auto & number : numbers)
for (auto& number : numbers)
numArray.push_back(static_cast<unsigned int>(number));
basis->setMolecularOrbitalNumber(numArray);
}
Expand All @@ -390,16 +390,16 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
json moCoefficientsB = orbitals["betaCoefficients"];
if (isNumericArray(moCoefficients)) {
std::vector<double> coeffs;
for (auto & moCoefficient : moCoefficients)
for (auto& moCoefficient : moCoefficients)
coeffs.push_back(static_cast<double>(moCoefficient));
basis->setMolecularOrbitals(coeffs);
} else if (isNumericArray(moCoefficientsA) &&
isNumericArray(moCoefficientsB)) {
std::vector<double> coeffsA;
for (auto & i : moCoefficientsA)
for (auto& i : moCoefficientsA)
coeffsA.push_back(static_cast<double>(i));
std::vector<double> coeffsB;
for (auto & i : moCoefficientsB)
for (auto& i : moCoefficientsB)
coeffsB.push_back(static_cast<double>(i));
basis->setMolecularOrbitals(coeffsA, BasisSet::Alpha);
basis->setMolecularOrbitals(coeffsB, BasisSet::Beta);
Expand All @@ -416,16 +416,16 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
moCoefficientsB = orbSets[idx]["betaCoefficients"];
if (isNumericArray(moCoefficients)) {
std::vector<double> coeffs;
for (auto & moCoefficient : moCoefficients)
for (auto& moCoefficient : moCoefficients)
coeffs.push_back(static_cast<double>(moCoefficient));
basis->setMolecularOrbitals(coeffs, BasisSet::Paired, idx);
} else if (isNumericArray(moCoefficientsA) &&
isNumericArray(moCoefficientsB)) {
std::vector<double> coeffsA;
for (auto & i : moCoefficientsA)
for (auto& i : moCoefficientsA)
coeffsA.push_back(static_cast<double>(i));
std::vector<double> coeffsB;
for (auto & i : moCoefficientsB)
for (auto& i : moCoefficientsB)
coeffsB.push_back(static_cast<double>(i));
basis->setMolecularOrbitals(coeffsA, BasisSet::Alpha, idx);
basis->setMolecularOrbitals(coeffsB, BasisSet::Beta, idx);
Expand All @@ -444,23 +444,23 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
json frequencies = vibrations["frequencies"];
if (isNumericArray(frequencies)) {
Array<double> freqs;
for (auto & frequencie : frequencies) {
for (auto& frequencie : frequencies) {
freqs.push_back(static_cast<double>(frequencie));
}
molecule.setVibrationFrequencies(freqs);
}
json intensities = vibrations["intensities"];
if (isNumericArray(intensities)) {
Array<double> intens;
for (auto & intensitie : intensities) {
for (auto& intensitie : intensities) {
intens.push_back(static_cast<double>(intensitie));
}
molecule.setVibrationIRIntensities(intens);
}
json raman = vibrations["ramanIntensities"];
if (isNumericArray(raman)) {
Array<double> intens;
for (auto & i : raman) {
for (auto& i : raman) {
intens.push_back(static_cast<double>(i));
}
molecule.setVibrationRamanIntensities(intens);
Expand All @@ -473,7 +473,7 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
Array<Vector3> mode;
mode.resize(arr.size() / 3);
double* ptr = &mode[0][0];
for (auto & j : arr) {
for (auto& j : arr) {
*(ptr++) = static_cast<double>(j);
}
disps.push_back(mode);
Expand All @@ -483,6 +483,22 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule)
}
}

// properties
if (jsonRoot.find("properties") != jsonRoot.end()) {
json properties = jsonRoot["properties"];
if (properties.is_object()) {
if (properties.find("totalCharge") != properties.end()) {
molecule.setData("totalCharge",
static_cast<int>(properties["totalCharge"]));
}
if (properties.find("totalSpinMultiplicity") != properties.end()) {
molecule.setData("totalSpinMultiplicity",
static_cast<int>(properties["totalSpinMultiplicity"]));
}
// todo - put everything into the data map
}
}

if (jsonRoot.find("layer") != jsonRoot.end()) {
auto names = LayerManager::getMoleculeInfo(&molecule);
json visible = jsonRoot["layer"]["visible"];
Expand Down Expand Up @@ -537,6 +553,15 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
root["inchi"] = molecule.data("inchi").toString().c_str();
}

json properties;
properties["totalCharge"] = molecule.totalCharge();
if (molecule.data("totalSpinMultiplicity").type() == Variant::Int)
properties["totalSpinMultiplicity"] =
molecule.data("totalSpinMultiplicity").toInt();
else
properties["totalSpinMultiplicity"] = 1; // assume singlet if not specified
root["properties"] = properties;

if (molecule.unitCell()) {
json unitCell;
unitCell["a"] = molecule.unitCell()->a();
Expand Down Expand Up @@ -605,7 +630,7 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)

auto atomIndices = gaussian->atomIndices();
json shellToAtomMap;
for (unsigned int & atomIndice : atomIndices)
for (unsigned int& atomIndice : atomIndices)
shellToAtomMap.push_back(atomIndice);
basis["shellToAtomMap"] = shellToAtomMap;

Expand Down Expand Up @@ -648,22 +673,22 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
auto energies = gaussian->moEnergy();
if (energies.size() > 0) {
json energyData;
for (double & energie : energies) {
for (double& energie : energies) {
energyData.push_back(energie);
}
root["orbitals"]["energies"] = energyData;
}
auto occ = gaussian->moOccupancy();
if (occ.size() > 0) {
json occData;
for (unsigned char & it : occ)
for (unsigned char& it : occ)
occData.push_back(static_cast<int>(it));
root["orbitals"]["occupations"] = occData;
}
auto num = gaussian->moNumber();
if (num.size() > 0) {
json numData;
for (unsigned int & it : num)
for (unsigned int& it : num)
numData.push_back(it);
root["orbitals"]["numbers"] = numData;
}
Expand Down Expand Up @@ -726,7 +751,7 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
if (molecule.atomPositions3d().size() == molecule.atomCount()) {
// everything gets real-space Cartesians
json coords3d;
for (const auto & it : molecule.atomPositions3d()) {
for (const auto& it : molecule.atomPositions3d()) {
coords3d.push_back(it.x());
coords3d.push_back(it.y());
coords3d.push_back(it.z());
Expand All @@ -739,7 +764,7 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
Array<Vector3> fcoords;
CrystalTools::fractionalCoordinates(
*molecule.unitCell(), molecule.atomPositions3d(), fcoords);
for (auto & fcoord : fcoords) {
for (auto& fcoord : fcoords) {
coordsFractional.push_back(fcoord.x());
coordsFractional.push_back(fcoord.y());
coordsFractional.push_back(fcoord.z());
Expand All @@ -751,7 +776,7 @@ bool CjsonFormat::write(std::ostream& file, const Molecule& molecule)
// 2d positions:
if (molecule.atomPositions2d().size() == molecule.atomCount()) {
json coords2d;
for (const auto & it : molecule.atomPositions2d()) {
for (const auto& it : molecule.atomPositions2d()) {
coords2d.push_back(it.x());
coords2d.push_back(it.y());
}
Expand Down Expand Up @@ -903,4 +928,4 @@ vector<std::string> CjsonFormat::mimeTypes() const
return mime;
}

} // namespace Avogadro
} // namespace Avogadro::Io