Merge pull request cryos#92 from psavery/fill-uc

Added fill unit cell and periodic distance calc

avogadro/core/spacegroups.cpp

@@ -14,7 +14,19 @@
 
 ******************************************************************************/
 
+#include <algorithm> // for std::count()
+#include <cassert>
+#include <cctype> // for isdigit()
+#include <iostream>
+
+#include "array.h"
+#include "crystaltools.h"
+#include "molecule.h"
 #include "spacegroupdata.h"
+#include "unitcell.h"
+#include "utilities.h"
+#include "vector.h"
+
 #include "spacegroups.h"
 
 namespace Avogadro {
@@ -125,5 +137,166 @@ const char * SpaceGroups::setting(unsigned short hallNumber)
     return space_group_setting[0];
 }
 
+unsigned short SpaceGroups::transformsCount(unsigned short hallNumber)
+{
+  if (hallNumber <= 530) {
+    std::string s = transformsString(hallNumber);
+    return std::count(s.begin(), s.end(), ' ') + 1;
+  }
+  else {
+    return 0;
+  }
+}
+
+Real readTransformCoordinate(const std::string &coordinate,
+                             const Vector3 &v)
+{
+  // The coordinate should be at least 1 character
+  assert(coordinate.size() != 0);
+
+  Real ret = 0.0;
+  Index i = 0;
+  while (i < coordinate.size()) {
+    bool isNeg = false;
+    if (coordinate[i] == '-') {
+      isNeg = true;
+      ++i;
+      assert(i < coordinate.size());
+    }
+    // We assume we are adding, so no need for a boolean here
+    else if (coordinate[i] == '+') {
+      ++i;
+      assert(i < coordinate.size());
+    }
+
+    // Check to see if we have a digit
+    if (isdigit(coordinate[i])) {
+      // We SHOULD have a fraction. Also, we SHOULD only deal with single
+      // digit numbers. Add assertions to make sure this is the case.
+      assert(i + 2 < coordinate.size());
+      assert(coordinate[i + 1] == '/');
+      assert(isdigit(coordinate[i + 2]));
+      // Assert that this is a single digit number
+      if (coordinate.size() > i + 3)
+        assert(!isdigit(coordinate[i + 3]));
+      // Ancient methods used by our forefathers to cast a char to an int
+      Real numerator = coordinate[i] - '0';
+      Real denominator = coordinate[i + 2] - '0';
+      Real fraction = numerator / denominator;
+      fraction *= (isNeg) ? -1.0 : 1.0;
+
+      ret += fraction;
+      i += 3;
+    }
+    else if (coordinate[i] == 'x') {
+      ret += (isNeg) ? -1.0 * v[0] : v[0];
+      ++i;
+    }
+    else if (coordinate[i] == 'y') {
+      ret += (isNeg) ? -1.0 * v[1] : v[1];
+      ++i;
+    }
+    else if (coordinate[i] == 'z') {
+      ret += (isNeg) ? -1.0 * v[2] : v[2];
+      ++i;
+    }
+    else {
+      std::cerr << "In " << __FUNCTION__ << ", error reading string: '"
+                << coordinate << "'\n";
+      return 0;
+    }
+  }
+  return ret;
+}
+
+Vector3 getSingleTransform(const std::string &transform,
+                           const Vector3 &v)
+{
+  Vector3 ret;
+  std::vector<std::string> coordinates = split(transform, ',');
+
+  // This should be 3 in size. Something very bad happened if it is not.
+  assert(coordinates.size() == 3);
+
+  ret[0] = readTransformCoordinate(coordinates[0], v);
+  ret[1] = readTransformCoordinate(coordinates[1], v);
+  ret[2] = readTransformCoordinate(coordinates[2], v);
+  return ret;
+}
+
+Array<Vector3> SpaceGroups::getTransforms(unsigned short hallNumber,
+                                          const Vector3 &v)
+{
+  if (hallNumber == 0 || hallNumber > 530)
+    return Array<Vector3>();
+
+  Array<Vector3> ret;
+
+  std::string transformsStr = transformsString(hallNumber);
+  // These transforms are separated by spaces
+  std::vector<std::string> transforms = split(transformsStr, ' ');
+
+  for (Index i = 0; i < transforms.size(); ++i)
+    ret.push_back(getSingleTransform(transforms[i], v));
+
+  return ret;
+}
+
+void SpaceGroups::fillUnitCell(Molecule &mol, unsigned short hallNumber,
+                               double cartTol)
+{
+  if (!mol.unitCell())
+    return;
+  UnitCell *uc = mol.unitCell();
+
+  Array<unsigned char> atomicNumbers = mol.atomicNumbers();
+  Array<Vector3> positions = mol.atomPositions3d();
+  Index numAtoms = mol.atomCount();
+
+  // We are going to loop through the original atoms. That is why
+  // we have numAtoms cached instead of using atomCount().
+  for (Index i = 0; i < numAtoms; ++i) {
+    unsigned char atomicNum = atomicNumbers[i];
+    Vector3 pos = uc->toFractional(positions[i]);
+
+    Array<Vector3> newAtoms = getTransforms(hallNumber, pos);
+
+    // We skip 0 because it is the original atom.
+    for (Index j = 1; j < newAtoms.size(); ++j) {
+      // The new atoms are in fractional coordinates. Convert to cartesian.
+      Vector3 newCandidate = uc->toCartesian(newAtoms[j]);
+
+      // If there is already an atom in this location within a
+      // certain tolerance, do not add the atom.
+      bool atomAlreadyPresent = false;
+      for (Index k = 0; k < mol.atomCount(); k++) {
+        // If it does not have the same atomic number, skip over it.
+        if (mol.atomicNumber(k) != atomicNum)
+          continue;
+        Real distance = uc->distance(mol.atomPosition3d(k), newCandidate);
+        if (distance <= cartTol)
+          atomAlreadyPresent = true;
+      }
+
+      // If there is already an atom present here, just continue
+      if (atomAlreadyPresent)
+        continue;
+
+      // If we got this far, add the atom!
+      Atom newAtom = mol.addAtom(atomicNum);
+      newAtom.setPosition3d(newCandidate);
+    }
+  }
+  CrystalTools::wrapAtomsToUnitCell(mol);
+}
+
+const char * SpaceGroups::transformsString(unsigned short hallNumber)
+{
+  if (hallNumber <= 530)
+    return space_group_transforms[hallNumber];
+  else
+    return "";
+}
+
 } // end Core namespace
 } // end Avogadro namespace

avogadro/core/spacegroups.h

@@ -91,6 +91,36 @@ class AVOGADROCORE_EXPORT SpaceGroups
    * space group.
    */
   static const char * setting(unsigned short hallNumber);
+
+  /**
+   * Get the number of transforms for a given hall number.
+   * If an invalid hall number is given, 0 will be returned.
+   */
+  static unsigned short transformsCount(unsigned short hallNumber);
+
+  /**
+   * Get an array of transforms for a given hall number and a vector v.
+   * The vector should be in fractional coordinates.
+   * If an invalid hall number is given, an empty array will be returned.
+   */
+  static Array<Vector3> getTransforms(unsigned short hallNumber,
+                                      const Vector3 &v);
+
+  /**
+   * Fill a crystal with atoms by using transforms from a hall number.
+   * Nothing will be done if the molecule does not have a unit cell.
+   * The cartesian tolerance is used to check if an atom is already
+   * present at that location. If there is another atom within that
+   * distance, the new atom will not be placed there.
+   */
+  static void fillUnitCell(Molecule &mol, unsigned short hallNumber,
+                           double cartTol = 1e-5);
+
+private:
+  /**
+   * Get the transforms string stored in the database.
+   */
+  static const char * transformsString(unsigned short hallNumber);
 };
 
 } // end Core namespace

avogadro/core/unitcell.h

@@ -127,6 +127,23 @@ class AVOGADROCORE_EXPORT UnitCell
   void wrapCartesian(const Vector3 &cart, Vector3 &wrapped) const;
   /** @} */
 
+  /**
+   * Find the minimum fractional image of a fractional vector @a v.
+   * A minimum image has all fractional coordinates between -0.5 and 0.5.
+   */
+  static Vector3 minimumImageFractional(const Vector3 &v);
+
+  /**
+   * Find the minimum image of a Cartesian vector @a v.
+   * A minimum image has all fractional coordinates between -0.5 and 0.5
+   */
+  Vector3 minimumImage(const Vector3 &v) const;
+
+  /**
+   * Find the shortest distance between vectors @a v1 and @a v2.
+   */
+  Real distance(const Vector3 &v1, const Vector3 &v2) const;
+
 private:
   static Real signedAngleRadians(const Vector3 &v1, const Vector3 &v2,
                                  const Vector3 &axis);
@@ -319,6 +336,24 @@ inline void UnitCell::wrapCartesian(const Vector3 &cart, Vector3 &wrapped) const
   toCartesian(wrapped, wrapped);
 }
 
+inline Vector3 UnitCell::minimumImageFractional(const Vector3 &v)
+{
+  Real x = v[0] - rint(v[0]);
+  Real y = v[1] - rint(v[1]);
+  Real z = v[2] - rint(v[2]);
+  return Vector3(x, y, z);
+}
+
+inline Vector3 UnitCell::minimumImage(const Vector3 &v) const
+{
+  return toCartesian(minimumImageFractional(toFractional(v)));
+}
+
+inline Real UnitCell::distance(const Vector3 &v1, const Vector3 &v2) const
+{
+  return std::fabs(minimumImage(v1 - v2).norm());
+}
+
 } // namespace Core
 } // namespace Avogadro
 

avogadro/qtplugins/spacegroup/spacegroup.cpp

@@ -24,12 +24,21 @@
 #include <avogadro/qtgui/molecule.h>
 #include <avogadro/qtgui/rwmolecule.h>
 
+#include <QtWidgets/QAbstractItemView>
 #include <QtWidgets/QAction>
+#include <QtWidgets/QDialogButtonBox>
+#include <QtWidgets/QHeaderView>
 #include <QtWidgets/QInputDialog>
+#include <QtWidgets/QLayout>
 #include <QtWidgets/QMessageBox>
+#include <QtWidgets/QScrollBar>
+#include <QtWidgets/QTableView>
+#include <QtWidgets/QVBoxLayout>
 
 #include <QtCore/QStringList>
 
+#include <QtGui/QStandardItemModel>
+
 #include <sstream>
 
 using Avogadro::Core::AvoSpglib;
@@ -49,6 +58,7 @@ SpaceGroup::SpaceGroup(QObject *parent_) :
   m_reduceToPrimitiveAction(new QAction(this)),
   m_conventionalizeCellAction(new QAction(this)),
   m_symmetrizeAction(new QAction(this)),
+  m_fillUnitCellAction(new QAction(this)),
   m_setToleranceAction(new QAction(this))
 {
   m_perceiveSpaceGroupAction->setText(tr("Perceive Space Group"));
@@ -74,6 +84,11 @@ SpaceGroup::SpaceGroup(QObject *parent_) :
   m_actions.push_back(m_symmetrizeAction);
   m_symmetrizeAction->setProperty("menu priority", 60);
 
+  m_fillUnitCellAction->setText(tr("Fill Unit Cell"));
+  connect(m_fillUnitCellAction, SIGNAL(triggered()), SLOT(fillUnitCell()));
+  m_actions.push_back(m_fillUnitCellAction);
+  m_fillUnitCellAction->setProperty("menu priority", 50);
+
   m_setToleranceAction->setText(tr("Set Tolerance"));
   connect(m_setToleranceAction, SIGNAL(triggered()), SLOT(setTolerance()));
   m_actions.push_back(m_setToleranceAction);
@@ -258,6 +273,56 @@ void SpaceGroup::symmetrize()
   }
 }
 
+void SpaceGroup::fillUnitCell()
+{
+  QStandardItemModel spacegroups;
+  QStringList modelHeader;
+  modelHeader << tr("International")
+              << tr("Hall")
+              << tr("Hermann-Mauguin");
+  spacegroups.setHorizontalHeaderLabels(modelHeader);
+  for (unsigned int i = 1; i <= 530; ++i) {
+    QList<QStandardItem*> row;
+    row << new QStandardItem(QString::number(
+                                 Core::SpaceGroups::internationalNumber(i)))
+        << new QStandardItem(QString(Core::SpaceGroups::hallSymbol(i)))
+        << new QStandardItem(QString(
+                                 Core::SpaceGroups::internationalShort(i)));
+    spacegroups.appendRow(row);
+  }
+
+  QDialog dialog;
+  dialog.setLayout(new QVBoxLayout);
+  dialog.setWindowTitle(tr("Select Space Group"));
+  QTableView *view = new QTableView;
+  view->setSelectionBehavior(QAbstractItemView::SelectRows);
+  view->setSelectionMode(QAbstractItemView::SingleSelection);
+  view->setCornerButtonEnabled(false);
+  view->setVerticalScrollMode(QAbstractItemView::ScrollPerPixel);
+  view->setHorizontalScrollBarPolicy(Qt::ScrollBarAlwaysOff);
+  view->verticalHeader()->hide();
+  view->setModel(&spacegroups);
+  dialog.layout()->addWidget(view);
+  view->selectRow(0);
+  view->resizeColumnsToContents();
+  view->resizeRowsToContents();
+  view->setMinimumWidth(view->horizontalHeader()->length()
+                        + view->verticalScrollBar()->sizeHint().width());
+  connect(view, SIGNAL(activated(QModelIndex)), &dialog, SLOT(accept()));
+  QDialogButtonBox *buttons =
+    new QDialogButtonBox(QDialogButtonBox::Ok | QDialogButtonBox::Cancel);
+  connect(buttons, SIGNAL(accepted()), &dialog, SLOT(accept()));
+  connect(buttons, SIGNAL(rejected()), &dialog, SLOT(reject()));
+  dialog.layout()->addWidget(buttons);
+  if (dialog.exec() != QDialog::Accepted)
+    return;
+
+  unsigned short hallNumber = view->currentIndex().row() + 1;
+
+  Core::SpaceGroups::fillUnitCell(*m_molecule, hallNumber, m_spgTol);
+  m_molecule->emitChanged(Molecule::Added | Molecule::Atoms);
+}
+
 void SpaceGroup::setTolerance()
 {
   bool ok;

avogadro/qtplugins/spacegroup/spacegroup.h

@@ -49,6 +49,7 @@ private slots:
   void reduceToPrimitive();
   void conventionalizeCell();
   void symmetrize();
+  void fillUnitCell();
   void setTolerance();
 
 private:
@@ -60,6 +61,7 @@ private slots:
   QAction *m_reduceToPrimitiveAction;
   QAction *m_conventionalizeCellAction;
   QAction *m_symmetrizeAction;
+  QAction *m_fillUnitCellAction;
   QAction *m_setToleranceAction;
 };
 

tests/core/CMakeLists.txt

@@ -16,7 +16,7 @@ set(tests
   Molecule
   Mutex
   RingPerceiver
-  Spglib
+  Spacegroup
   Utilities
   UnitCell
   Variant