Merge branch 'master' of github.com:openchemistry/avogadrolibs into a…

…dd-bond-labels

Signed-off-by: Geoff Hutchison <[email protected]>

README.md

@@ -1,7 +1,7 @@
 # ![Avogadro 2][Avogadro2Logo] Avogadro 2
 
 [![Latest Release](https://img.shields.io/github/v/release/openchemistry/avogadrolibs)](https://github.com/OpenChemistry/avogadrolibs/releases) [![BSD License](https://img.shields.io/github/license/openchemistry/avogadrolibs)](https://github.com/OpenChemistry/avogadrolibs/blob/master/LICENSE) [![Build Status](https://img.shields.io/github/actions/workflow/status/openchemistry/avogadrolibs/build_cmake.yml?branch=master)](https://github.com/OpenChemistry/avogadrolibs/actions) [![Codacy Badge](https://app.codacy.com/project/badge/Grade/44bb12662c564ed8a27ee8a7fd89ed50)](https://app.codacy.com/gh/OpenChemistry/avogadrolibs/dashboard?utm_source=gh&utm_medium=referral&utm_content=&utm_campaign=Badge_grade)  
-[![Download Count](https://avogadro.cc/downloads.svg)](https://github.com/OpenChemistry/avogadrolibs/releases) [![Citation Count](https://avogadro.cc/citations.svg)](http://doi.org/10.1186/1758-2946-4-17)  
+[![Download Count](https://avogadro.cc/downloads.svg?readme)](https://github.com/OpenChemistry/avogadrolibs/releases) [![Citation Count](https://avogadro.cc/citations.svg?readme)](http://doi.org/10.1186/1758-2946-4-17)  
 [![PRs Welcome](https://img.shields.io/badge/PRs-welcome-brightgreen.svg?style=flat)](http://makeapullrequest.com) [![GitHub contributors](https://img.shields.io/github/contributors/openchemistry/avogadrolibs.svg?style=flat&color=0bf)](https://github.com/OpenChemistry/avogadrolibs/graphs/contributors)  [![OpenCollective Backers](https://img.shields.io/opencollective/all/open-chemistry)](https://opencollective.com/open-chemistry)
 
 ## Introduction

avogadro/core/atomutilities.h

@@ -8,8 +8,6 @@
 
 #include "avogadrocoreexport.h"
 
-#include "avogadrocoreexport.h"
-
 #include <avogadro/core/molecule.h>
 
 #include <vector>

avogadro/core/coordinateblockgenerator.h

@@ -6,7 +6,7 @@
 #ifndef AVOGADRO_CORE_COORDINATEBLOCKGENERATOR_H
 #define AVOGADRO_CORE_COORDINATEBLOCKGENERATOR_H
 
-#include <avogadrocoreexport.h>
+#include "avogadrocoreexport.h"
 
 #include <sstream>
 #include <string>

avogadro/core/gaussiansettools.cpp

@@ -11,7 +11,6 @@
 
 #include <iostream>
 
-
 using std::vector;
 
 namespace Avogadro::Core {
@@ -22,9 +21,7 @@ GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol)
     m_basis = dynamic_cast<GaussianSet*>(m_molecule->basisSet());
 }
 
-GaussianSetTools::~GaussianSetTools()
-{
-}
+GaussianSetTools::~GaussianSetTools() {}
 
 bool GaussianSetTools::calculateMolecularOrbital(Cube& cube, int moNumber) const
 {
@@ -56,6 +53,12 @@ double GaussianSetTools::calculateMolecularOrbital(const Vector3& position,
 
 bool GaussianSetTools::calculateElectronDensity(Cube& cube) const
 {
+  const MatrixX& matrix = m_basis->densityMatrix();
+  if (matrix.rows() == 0 || matrix.cols() == 0) {
+    // we don't have a density matrix, so generate one
+    m_basis->generateDensityMatrix();
+  }
+
   for (size_t i = 0; i < cube.data()->size(); ++i) {
     Vector3 pos = cube.position(i);
     cube.setValue(i, calculateElectronDensity(pos));
@@ -67,6 +70,7 @@ double GaussianSetTools::calculateElectronDensity(const Vector3& position) const
 {
   const MatrixX& matrix = m_basis->densityMatrix();
   int matrixSize(static_cast<int>(m_basis->moMatrix().rows()));
+
   if (matrix.rows() != matrixSize || matrix.cols() != matrixSize) {
     return 0.0;
   }
@@ -153,7 +157,7 @@ inline vector<double> GaussianSetTools::calculateValues(
   // Calculate the deltas for the position
   for (Index i = 0; i < atomsSize; ++i) {
     deltas.emplace_back(pos -
-                     (m_molecule->atom(i).position3d() * ANGSTROM_TO_BOHR));
+                        (m_molecule->atom(i).position3d() * ANGSTROM_TO_BOHR));
     dr2.push_back(deltas[i].squaredNorm());
   }
 
@@ -246,7 +250,7 @@ inline void GaussianSetTools::pointD(unsigned int moIndex, const Vector3& delta,
        i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
     // Calculate the common factor
     double tmpGTO = exp(-gtoA[i] * dr2);
-    for (double & component : components)
+    for (double& component : components)
       component += gtoCN[cIndex++] * tmpGTO;
   }
 
@@ -279,7 +283,7 @@ inline void GaussianSetTools::pointD5(unsigned int moIndex,
        i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
     // Calculate the common factor
     double tmpGTO = exp(-gtoA[i] * dr2);
-    for (double & component : components)
+    for (double& component : components)
       component += gtoCN[cIndex++] * tmpGTO;
   }
 
@@ -318,7 +322,7 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta,
        i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
     // Calculate the common factor
     double tmpGTO = exp(-gtoA[i] * dr2);
-    for (double & component : components)
+    for (double& component : components)
       component += gtoCN[cIndex++] * tmpGTO;
   }
 
@@ -336,16 +340,7 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta,
   double componentsF[10] = {
     // molden order
     // e.g https://gau2grid.readthedocs.io/en/latest/order.html
-    xxx,
-    yyy,
-    zzz,
-    xyy,
-    xxy,
-    xxz,
-    xzz,
-    yzz,
-    yyz,
-    xyz
+    xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
   };
 
   for (int i = 0; i < 10; ++i)
@@ -371,7 +366,7 @@ inline void GaussianSetTools::pointF7(unsigned int moIndex,
        i < m_basis->gtoIndices()[moIndex + 1]; ++i) {
     // Calculate the common factor
     double tmpGTO = exp(-gtoA[i] * dr2);
-    for (double & component : components)
+    for (double& component : components)
       component += gtoCN[cIndex++] * tmpGTO;
   }
 
@@ -418,4 +413,4 @@ final normalization
     values[baseIndex + i] += components[i] * componentsF[i];
 }
 
-} // End Avogadro namespace
+} // namespace Avogadro::Core

avogadro/core/molecule.cpp

@@ -9,11 +9,13 @@
 #include "color3f.h"
 #include "cube.h"
 #include "elements.h"
+#include "gaussianset.h"
 #include "layermanager.h"
 #include "mdlvalence_p.h"
 #include "mesh.h"
 #include "neighborperceiver.h"
 #include "residue.h"
+#include "slaterset.h"
 #include "unitcell.h"
 
 #include <algorithm>
@@ -75,6 +77,55 @@ Molecule::Molecule(const Molecule& other)
   }
 }
 
+void Molecule::readProperties(const Molecule& other)
+{
+  m_label = other.m_label;
+  m_colors = other.m_colors;
+  // merge data maps by iterating through other's map
+  for (auto it = other.m_data.constBegin(); it != other.m_data.constEnd();
+       ++it) {
+    // even if we have the same key, we want to overwrite
+    m_data.setValue(it->first, it->second);
+  }
+  // merge partial charge maps
+  for (auto it = other.m_partialCharges.cbegin();
+       it != other.m_partialCharges.cend(); ++it) {
+    m_partialCharges[it->first] = it->second;
+  }
+
+  // copy orbital information
+  SlaterSet* slaterSet = dynamic_cast<SlaterSet*>(other.m_basisSet);
+  if (slaterSet != nullptr) {
+    m_basisSet = slaterSet->clone();
+    m_basisSet->setMolecule(this);
+  }
+  GaussianSet* gaussianSet = dynamic_cast<GaussianSet*>(other.m_basisSet);
+  if (gaussianSet != nullptr) {
+    m_basisSet = gaussianSet->clone();
+    m_basisSet->setMolecule(this);
+  }
+
+  // copy over spectra information
+  if (other.m_vibrationFrequencies.size() > 0) {
+    m_vibrationFrequencies = other.m_vibrationFrequencies;
+    m_vibrationIRIntensities = other.m_vibrationIRIntensities;
+    m_vibrationRamanIntensities = other.m_vibrationRamanIntensities;
+    m_vibrationLx = other.m_vibrationLx;
+  }
+
+  // Copy over any meshes
+  for (Index i = 0; i < other.meshCount(); ++i) {
+    Mesh* m = addMesh();
+    *m = *other.mesh(i);
+  }
+
+  // Copy over any cubes
+  for (Index i = 0; i < other.cubeCount(); ++i) {
+    Cube* c = addCube();
+    *c = *other.cube(i);
+  }
+}
+
 Molecule::Molecule(Molecule&& other) noexcept
   : m_data(other.m_data), m_partialCharges(std::move(other.m_partialCharges)),
     m_customElementMap(std::move(other.m_customElementMap)),
@@ -1172,6 +1223,11 @@ bool Molecule::setCoordinate3d(int coord)
   return false;
 }
 
+void Molecule::clearCoordinate3d()
+{
+  m_coordinates3d.clear();
+}
+
 Array<Vector3> Molecule::coordinate3d(int index) const
 {
   return m_coordinates3d[index];

avogadro/core/molecule.h

@@ -73,6 +73,13 @@ class AVOGADROCORE_EXPORT Molecule
   /** Destroys the molecule object. */
   virtual ~Molecule();
 
+  /**
+   * Adds the properties from the supplied
+   * molecule to this molecule. Does not otherwise
+   * modify atoms / bonds / residues, etc.
+   */
+  void readProperties(const Molecule& other);
+
   /** Sets the data value with @p name to @p value. */
   void setData(const std::string& name, const Variant& value);
 
@@ -586,6 +593,11 @@ class AVOGADROCORE_EXPORT Molecule
   Array<Vector3> coordinate3d(int index) const;
   bool setCoordinate3d(const Array<Vector3>& coords, int index);
 
+  /**
+   * Clear coordinate sets (except the default set)
+   */
+  void clearCoordinate3d();
+
   /**
    * Timestep property is used when molecular dynamics trajectories are read
    */

avogadro/core/variant-inline.h

@@ -13,16 +13,28 @@
 namespace Avogadro {
 namespace Core {
 
-inline Variant::Variant() : m_type(Null)
-{
-}
+inline Variant::Variant() : m_type(Null) {}
 
 template <typename T>
 inline Variant::Variant(T v) : m_type(Null)
 {
   setValue(v);
 }
 
+template <>
+inline Variant::Variant(const char* v) : m_type(String)
+{
+  m_value.string = new std::string(v);
+}
+
+template <>
+inline Variant::Variant(const MatrixXf& v) : m_type(Matrix)
+{
+  MatrixX* m = new MatrixX(v.rows(), v.cols());
+  *m = v.cast<double>();
+  m_value.matrix = m;
+}
+
 inline Variant::Variant(const Variant& variant) : m_type(variant.type())
 {
   if (m_type == String)
@@ -444,7 +456,7 @@ inline T Variant::lexical_cast(const std::string& str)
   return value;
 }
 
-} // end Core namespace
-} // end Avogadro namespace
+} // namespace Core
+} // namespace Avogadro
 
 #endif // AVOGADRO_CORE_VARIANT_INLINE_H

avogadro/core/variant.h

@@ -146,8 +146,8 @@ class AVOGADROCORE_EXPORT Variant
   } m_value;
 };
 
-} // end Core namespace
-} // end Avogadro namespace
+} // namespace Core
+} // namespace Avogadro
 
 #include "variant-inline.h"
 

avogadro/io/cjsonformat.cpp

@@ -214,21 +214,6 @@ bool CjsonFormat::deserialize(std::istream& file, Molecule& molecule,
     }
   }
 
-  // Partial charges are optional, but if present should be loaded.
-  json partialCharges = atoms["partialCharges"];
-  if (partialCharges.is_object()) {
-    // keys are types, values are arrays of charges
-    for (auto& kv : partialCharges.items()) {
-      MatrixX charges(atomCount, 1);
-      if (isNumericArray(kv.value()) && kv.value().size() == atomCount) {
-        for (size_t i = 0; i < kv.value().size(); ++i) {
-          charges(i, 0) = kv.value()[i];
-        }
-        molecule.setPartialCharges(kv.key(), charges);
-      }
-    }
-  }
-
   // Bonds are optional, but if present should be loaded.
   json bonds = jsonRoot["bonds"];
   if (bonds.is_object() && isNumericArray(bonds["connections"]["index"])) {
@@ -617,6 +602,21 @@ bool CjsonFormat::deserialize(std::istream& file, Molecule& molecule,
     }
   }
 
+  // Partial charges are optional, but if present should be loaded.
+  json partialCharges = atoms["partialCharges"];
+  if (partialCharges.is_object()) {
+    // keys are types, values are arrays of charges
+    for (auto& kv : partialCharges.items()) {
+      MatrixX charges(atomCount, 1);
+      if (isNumericArray(kv.value()) && kv.value().size() == atomCount) {
+        for (size_t i = 0; i < kv.value().size(); ++i) {
+          charges(i, 0) = kv.value()[i];
+        }
+        molecule.setPartialCharges(kv.key(), charges);
+      }
+    }
+  }
+
   if (jsonRoot.find("layer") != jsonRoot.end()) {
     auto names = LayerManager::getMoleculeInfo(&molecule);
     json visible = jsonRoot["layer"]["visible"];
@@ -923,6 +923,20 @@ bool CjsonFormat::serialize(std::ostream& file, const Molecule& molecule,
     if (hasCustomColors)
       root["atoms"]["colors"] = colors;
 
+    // check for partial charges
+    auto partialCharges = molecule.partialChargeTypes();
+    if (!partialCharges.empty()) {
+      // add them to the atoms object
+      for (const auto& type : partialCharges) {
+        MatrixX chargesMatrix = molecule.partialCharges(type);
+        json charges;
+        for (Index i = 0; i < molecule.atomCount(); ++i) {
+          charges.push_back(chargesMatrix(i, 0));
+        }
+        root["atoms"]["partialCharges"][type] = charges;
+      }
+    }
+
     // 3d positions:
     if (molecule.atomPositions3d().size() == molecule.atomCount()) {
       // everything gets real-space Cartesians