Calculate density matrix when needed for electron density surfaces

avogadro/core/gaussiansettools.cpp

@@ -55,7 +55,13 @@ 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 +73,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;
   }

avogadro/core/internalcoordinates.cpp

@@ -0,0 +1,108 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#include "internalcoordinates.h"
+#include "matrix.h"
+
+#include <cmath>
+
+namespace Avogadro::Core {
+
+Array<Vector3> internalToCartesian(
+  const Molecule& molecule, const Array<InternalCoordinate>& internalCoords)
+{
+  Array<Vector3> coords(molecule.atomCount());
+  Vector3 ab;
+  Vector3 bc;
+  Vector3 n;
+  Matrix3 m;
+
+  for (Index i = 0; i < molecule.atomCount()++ i) {
+    Real sinTheta, cosTheta, sinPhi, cosPhi;
+    Real length = internalCoords[i].length;
+    Real angle = internalCoords[i].angle;
+    Real dihedral = internalCoords[i].dihedral;
+
+    switch (i) {
+      case 0:
+        coords[i] = Vector3(0.0, 0.0, 0.0);
+        break;
+      case 1:
+        coords[i] = Vector3(length, 0.0, 0.0);
+        ab = Vector3(1.0, 0.0, 0.0); // normalized
+        break;
+      case 2:
+        sinTheta = std::sin(angle*DEG_TO_RAD);
+        cosTheta = std::cos(angle*DEG_TO_RAD);
+        coords[i] = Vector3(coords[i - 1] + length * cosTheta,
+                            length * sinTheta, 0.0);
+        bc = (coords[i] - coords[i - 1]) / length;
+        break;
+      default: 
+        // NeRF formula
+        // see J. Comp. Chem. Vol. 26, No. 10, p. 1063-1068 (2005)
+        // https://doi.org/10.1002/jcc.20237
+        sinTheta = std::sin(internalCoords[i].angle);
+        cosTheta = std::cos(internalCoords[i].angle);
+        sinPhi = std::sin(internalCoords[i].dihedral);
+        cosPhi = std::cos(internalCoords[i].dihedral);
+
+        n = (ab.cross(bc)).normalized();
+
+        // D2 in the paper nomenclature (page 1066)
+        // D2 = (RcosTheta, R*cosPhi*sinTheta, R*sinPhi*sinTheta)
+        coords[i] = Vector3(length*cosTheta, R*sinTheta*cosPhi, R*sinTheta*sinPhi);
+        m.col(0) = bc;
+        m.col(1) = n.cross(bc);
+        m.col(2) = n;
+        coords[i] = m * coords[i] + coords[i - 1];
+
+        // set up the vectors for the next iteration
+        ab = bc;
+        // we know the length, so we don't need .normalized()
+        //  .. save ourself a square root
+        bc = (coords[i] - coords[i - 1]) / length;
+        break;
+    }
+  }
+
+  return coords;
+}
+
+Array<InternalCoordinate> cartesianToInternal(const Molecule& molecule)
+{
+  Array<InternalCoordinate> internalCoords(molecule.atomCount());
+  /*
+  for (Index i = 0; i < molecule.numAtoms(); ++i) {
+
+    Vector3 a = molecule.atom(i).pos();
+    Vector3 b = molecule.atom(j).pos();
+    Vector3 c = molecule.atom(k).pos();
+    Vector3 ab = b - a;
+    Vector3 bc = c - b;
+    Vector3 ac = c - a;
+    Real lengthAB = ab.length();
+    Real lengthBC = bc.length();
+    Real lengthAC = ac.length();
+    Real angleAB = std::acos(ab.dot(bc) / (lengthAB * lengthBC));
+    Real angleAC = std::acos(ac.dot(bc) / (lengthAC * lengthBC));
+    Real angleBC = std::acos(bc.dot(ab) / (lengthBC * lengthAB));
+    Real dihedral =
+      std::acos((ab.dot(ac) * std::sin(angleAB) * std::sin(angleAC)) /
+                (lengthAB * lengthAC));
+
+    InternalCoordinate coord;
+    coord.a = i;
+    coord.b = j;
+    coord.c = k;
+    coord.length = lengthAB;
+    coord.angle = angleAB;
+    coord.dihedral = dihedral;
+    internalCoords.append(coord);
+  */
+  return internalCoords;
+}
+
+} // end namespace Avogadro::Core

avogadro/core/internalcoordinates.h

@@ -0,0 +1,35 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#ifndef AVOGADRO_CORE_INTERNALCOORDINATES_H
+#define AVOGADRO_CORE_INTERNALCOORDINATES_H
+
+#include "array.h"
+#include "avogadrocore.h"
+#include "vector.h"
+
+namespace Avogadro {
+namespace Core {
+
+/** A simple struct to define internal / z-matrix coordinates. */
+struct InternalCoordinate
+{
+  Index a;
+  Index b;
+  Index c;
+  Real length;
+  Real angle;
+  Real dihedral;
+};
+
+Array<Vector3> internalToCartesian(
+  const Molecule& molecule, const Array<InternalCoordinate>& internalCoords);
+
+Array<InternalCoordinate> cartesianToInternal(const Molecule& molecule);
+
+} // namespace Core
+} // namespace Avogadro
+
+#endif // AVOGADRO_CORE_INTERNALCOORDINATES_H

avogadro/qtplugins/CMakeLists.txt

@@ -112,6 +112,7 @@ add_subdirectory(forcefield)
 add_subdirectory(hydrogens)
 add_subdirectory(importpqr)
 add_subdirectory(insertdna)
+add_subdirectory(insertpeptide)
 add_subdirectory(label)
 add_subdirectory(lineformatinput)
 add_subdirectory(manipulator)

avogadro/qtplugins/insertpeptide/CMakeLists.txt

@@ -0,0 +1,32 @@
+include(ExternalProject)
+
+avogadro_plugin(InsertPeptide
+  "Insert oligopeptide chains"
+  ExtensionPlugin
+  insertpeptide.h
+  InsertPeptide
+  insertpeptide.cpp
+  insertpeptidedialog.ui
+)
+
+# Install the fragments
+set(_fragments "${AvogadroLibs_SOURCE_DIR}/../fragments")
+
+# Look in parallel directory for the molecule fragment repository
+if(NOT EXISTS "${_fragments}")
+# download molecules...
+ExternalProject_Add(fragments
+    GIT_REPOSITORY https://github.com/openchemistry/fragments 
+    GIT_TAG main
+    # or https://github.com/OpenChemistry/molecules/archive/refs/heads/master.zip
+    SOURCE_DIR "${AvogadroLibs_SOURCE_DIR}/../fragments"
+    CONFIGURE_COMMAND ""
+    BUILD_COMMAND ""
+    INSTALL_COMMAND ""
+)
+endif()
+
+install(DIRECTORY "${AvogadroLibs_SOURCE_DIR}/../fragments"
+  DESTINATION "${INSTALL_DATA_DIR}/avogadro2"
+  PATTERN ".git" EXCLUDE
+)

avogadro/qtplugins/insertpeptide/insertpeptide.cpp

@@ -0,0 +1,80 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#include "insertpeptide.h"
+#include "ui_insertpeptidedialog.h"
+
+#include <avogadro/qtgui/molecule.h>
+#include <avogadro/qtgui/rwmolecule.h>
+
+#include <QtCore/QDebug>
+
+#include <QtWidgets/QAction>
+
+using Avogadro::QtGui::Molecule;
+
+namespace Avogadro::QtPlugins {
+
+  class InsertPeptideDialog : public QDialog, public Ui::InsertPeptideDialog
+    {
+    public:
+    InsertPeptideDialog(QWidget *parent=0) : QDialog(parent) {
+        setWindowFlags(Qt::Dialog | Qt::Tool);
+        setupUi(this);
+      }
+    };
+
+InsertPeptide::InsertPeptide(QObject* parent_)
+  : Avogadro::QtGui::ExtensionPlugin(parent_), m_dialog(nullptr)
+{
+  auto* action = new QAction(tr("Peptide…"), this);
+  connect(action, SIGNAL(triggered()), SLOT(showDialog()));
+  m_actions.append(action);
+}
+
+InsertPeptide::~InsertPeptide()
+{
+}
+
+QList<QAction*> InsertPeptide::actions() const
+{
+  return m_actions;
+}
+
+QStringList InsertPeptide::menuPath(QAction* action) const
+{
+  return QStringList() << tr("&Build") << tr("&Insert");
+}
+
+void InsertPeptide::setMolecule(QtGui::Molecule* mol)
+{
+  m_molecule = mol;
+}
+
+void InsertPeptide::showDialog()
+{
+  if (m_molecule == nullptr)
+    return;
+
+  if (m_dialog == nullptr) {
+    m_dialog = new InsertPeptideDialog(qobject_cast<QWidget*>(parent()));
+  }
+
+  m_dialog->show();
+}
+
+void InsertPeptide::performInsert(const QString& sequence)
+{
+  if (m_molecule == nullptr)
+    return;
+
+  // read the file into the new fragment
+  Avogadro::QtGui::Molecule newMol(m_molecule->parent());
+
+  //m_molecule->undoMolecule()->appendMolecule(newMol, tr("Insert Peptide"));
+  emit requestActiveTool("Manipulator");
+}
+
+} // namespace Avogadro

avogadro/qtplugins/insertpeptide/insertpeptide.h

@@ -0,0 +1,50 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#ifndef AVOGADRO_QTPLUGINS_INSERTPEPTIDE_H
+#define AVOGADRO_QTPLUGINS_INSERTPEPTIDE_H
+
+#include <avogadro/qtgui/extensionplugin.h>
+
+namespace Avogadro {
+namespace QtPlugins {
+
+class InsertPeptideDialog;
+
+/**
+ * @brief Insert oligopeptide sequences
+ */
+class InsertPeptide : public QtGui::ExtensionPlugin
+{
+  Q_OBJECT
+public:
+  explicit InsertPeptide(QObject* parent_ = nullptr);
+  ~InsertPeptide() override;
+
+  QString name() const override { return tr("InsertPeptide"); }
+  QString description() const override
+  {
+    return tr("Insert oligopeptide sequences.");
+  }
+  QList<QAction*> actions() const override;
+  QStringList menuPath(QAction*) const override;
+
+public slots:
+  void setMolecule(QtGui::Molecule*) override;
+
+private slots:
+  void showDialog();
+  void performInsert(const QString& sequence);
+
+private:
+  QList<QAction*> m_actions;
+  QtGui::Molecule* m_molecule;
+  InsertPeptideDialog* m_dialog;
+};
+
+} // namespace QtPlugins
+} // namespace Avogadro
+
+#endif // AVOGADRO_QTPLUGINS_INSERTPEPTIDE_H