Added several spglib functions for crystals

Added space group perception, primitive reduction, and cell conventionalization.
In addition, the user can specify the tolerance for all of these features.
All of the features use spglib to perform them. When performing any action,
the program tells the user what the tolerance is currently set to.

avogadro/core/avospglib.cpp

@@ -83,5 +83,155 @@ unsigned short AvoSpglib::getHallNumber(const Molecule &mol, double cartTol)
   return hallNumber;
 }
 
+
+bool AvoSpglib::reduceToPrimitive(Molecule &mol, double cartTol)
+{
+  if (!mol.unitCell())
+    return false;
+
+  const UnitCell *uc = mol.unitCell();
+  Matrix3 cellMat = uc->cellMatrix();
+
+  double lattice[3][3];
+  // Spglib expects column vectors
+  for (Index i = 0; i < 3; ++i) {
+    for (Index j = 0; j < 3; ++j) {
+      lattice[i][j] = cellMat(i,j);
+    }
+  }
+
+  Index numAtoms = mol.atomCount();
+  double (*positions)[3] = new double[numAtoms][3];
+  int *types = new int[numAtoms];
+
+  const Array<unsigned char> &atomicNums = mol.atomicNumbers();
+  const Array<Vector3> &pos = mol.atomPositions3d();
+
+  // Positions need to be in fractional coordinates
+  for (Index i = 0; i < numAtoms; ++i) {
+    Vector3 fracCoords = uc->toFractional(pos[i]);
+    positions[i][0] = fracCoords[0];
+    positions[i][1] = fracCoords[1];
+    positions[i][2] = fracCoords[2];
+    types[i] = atomicNums[i];
+  }
+
+  // Run the spglib algorithm
+  Index newNumAtoms = spg_find_primitive(lattice, positions, types,
+                                         numAtoms, cartTol);
+
+  // If 0 is returned, the algorithm failed.
+  if (newNumAtoms == 0) {
+    delete [] positions;
+    delete [] types;
+    return false;
+  }
+
+  // Let's create a new molecule with the primitive information
+  Molecule newMol;
+
+  // First, we will make the unit cell
+  Matrix3 newCellMat;
+  for (Index i = 0; i < 3; ++i) {
+    for (Index j = 0; j < 3; ++j) {
+      newCellMat(i,j) = lattice[i][j];
+    }
+  }
+  UnitCell *newCell = new UnitCell(newCellMat);
+  newMol.setUnitCell(newCell);
+
+  // Next, add in the atoms
+  for (Index i = 0; i < newNumAtoms; ++i) {
+    Atom newAtom = newMol.addAtom(types[i]);
+    Vector3 newAtomPos(positions[i][0], positions[i][1], positions[i][2]);
+    // We must convert it back to cartesian before adding it
+    newAtom.setPosition3d(newCell->toCartesian(newAtomPos));
+  }
+
+  delete [] positions;
+  delete [] types;
+
+  // Set the new molecule
+  mol = newMol;
+  return true;
+}
+
+bool AvoSpglib::conventionalizeCell(Molecule &mol, double cartTol)
+{
+  if (!mol.unitCell())
+    return false;
+
+  const UnitCell *uc = mol.unitCell();
+  Matrix3 cellMat = uc->cellMatrix();
+
+  double lattice[3][3];
+  // Spglib expects column vectors
+  for (Index i = 0; i < 3; ++i) {
+    for (Index j = 0; j < 3; ++j) {
+      lattice[i][j] = cellMat(i,j);
+    }
+  }
+
+  Index numAtoms = mol.atomCount();
+  // spg_refine_cell() can cause the number of atoms to increase by
+  // as much as 4x. So, we must make these arrays at least 4x the number of
+  // atoms.
+  // See http://atztogo.github.io/spglib/api.html#spg-refine-cell
+  double (*positions)[3] = new double[numAtoms * 4][3];
+  int *types = new int[numAtoms * 4];
+
+  const Array<unsigned char> &atomicNums = mol.atomicNumbers();
+  const Array<Vector3> &pos = mol.atomPositions3d();
+
+  // Positions need to be in fractional coordinates
+  for (Index i = 0; i < numAtoms; ++i) {
+    Vector3 fracCoords = uc->toFractional(pos[i]);
+    positions[i][0] = fracCoords[0];
+    positions[i][1] = fracCoords[1];
+    positions[i][2] = fracCoords[2];
+    types[i] = atomicNums[i];
+  }
+
+  // Run the spglib algorithm
+  Index newNumAtoms = spg_refine_cell(lattice, positions, types,
+                                      numAtoms, cartTol);
+
+  // If 0 is returned, the algorithm failed.
+  if (newNumAtoms == 0) {
+    delete [] positions;
+    delete [] types;
+    return false;
+  }
+
+  // Let's create a new molecule with the conventional information
+  Molecule newMol;
+
+  // First, we will make the unit cell
+  Matrix3 newCellMat;
+  for (Index i = 0; i < 3; ++i) {
+    for (Index j = 0; j < 3; ++j) {
+      newCellMat(i,j) = lattice[i][j];
+    }
+  }
+
+  UnitCell *newCell = new UnitCell(newCellMat);
+  newMol.setUnitCell(newCell);
+
+  // Next, add in the atoms
+  for (Index i = 0; i < newNumAtoms; ++i) {
+    Atom newAtom = newMol.addAtom(types[i]);
+    Vector3 newAtomPos(positions[i][0], positions[i][1], positions[i][2]);
+    // We must convert it back to cartesian before adding it
+    newAtom.setPosition3d(newCell->toCartesian(newAtomPos));
+  }
+
+  delete [] positions;
+  delete [] types;
+
+  // Set the new molecule
+  mol = newMol;
+  return true;
+}
+
 } // end Core namespace
 } // end Avogadro namespace

avogadro/core/avospglib.h

@@ -42,7 +42,29 @@ class AVOGADROCORE_EXPORT AvoSpglib
    * @return The Hall number for the crystal.
    */
   static unsigned short getHallNumber(const Molecule &mol,
-                                      double cartTol = 0.05);
+                                      double cartTol = 1e-5);
+
+  /**
+   * Use spglib to reduce the crystal to a primitive cell. Unless the molecule
+   * is missing its unit cell, it will be edited by spglib.
+   *
+   * @param mol The molecule to be reduced to its primitive cell.
+   * @param cartTol The cartesian tolerance for spglib.
+   * @return False if the molecule has no unit cell or if the
+                   spglib algorithm failed. True otherwise.
+   */
+  static bool reduceToPrimitive(Molecule &mol, double cartTol = 1e-5);
+
+  /**
+   * Use spglib to refine the crystal to its conventional cell. Unless the
+   * molecule is missing its unit cell, it will be edited by spglib.
+   *
+   * @param mol The molecule to be conventionalized.
+   * @param cartTol The cartesian tolerance for spglib.
+   * @return False if the molecule has no unit cell or if the
+   *         spglib algorithm failed. True otherwise.
+   */
+  static bool conventionalizeCell(Molecule &mol, double cartTol = 1e-5);
 
 };
 

avogadro/qtplugins/CMakeLists.txt

@@ -112,6 +112,7 @@ add_subdirectory(openbabel)
 add_subdirectory(quantumoutput)
 add_subdirectory(playertool)
 add_subdirectory(povray)
+add_subdirectory(spacegroup)
 add_subdirectory(spectra)
 
 if(USE_MOLEQUEUE)

avogadro/qtplugins/crystal/crystal.cpp

@@ -55,45 +55,45 @@ Crystal::Crystal(QObject *parent_) :
   connect(m_importCrystalClipboardAction, SIGNAL(triggered()),
           SLOT(importCrystalClipboard()));
   m_actions.push_back(m_importCrystalClipboardAction);
-  m_importCrystalClipboardAction->setProperty("menu priority", 50);
+  m_importCrystalClipboardAction->setProperty("menu priority", 220);
 
   // this will be changed when the molecule is set:
   m_toggleUnitCellAction->setText(tr("Toggle Unit Cell"));
   connect(m_toggleUnitCellAction, SIGNAL(triggered()), SLOT(toggleUnitCell()));
   m_actions.push_back(m_toggleUnitCellAction);
-  m_toggleUnitCellAction->setProperty("menu priority", -1);
+  m_toggleUnitCellAction->setProperty("menu priority", 210);
 
   m_editUnitCellAction->setText(tr("Edit Unit Cell..."));
   connect(m_editUnitCellAction, SIGNAL(triggered()), SLOT(editUnitCell()));
   m_actions.push_back(m_editUnitCellAction);
-  m_editUnitCellAction->setProperty("menu priority", -50);
+  m_editUnitCellAction->setProperty("menu priority", 190);
 
   m_wrapAtomsToCellAction->setText(tr("&Wrap Atoms to Unit Cell"));
   connect(m_wrapAtomsToCellAction, SIGNAL(triggered()),
           SLOT(wrapAtomsToCell()));
   m_actions.push_back(m_wrapAtomsToCellAction);
-  m_wrapAtomsToCellAction->setProperty("menu priority", -200);
+  m_wrapAtomsToCellAction->setProperty("menu priority", 180);
 
   m_standardOrientationAction->setText(tr("Rotate to Standard &Orientation"));
   connect(m_standardOrientationAction, SIGNAL(triggered()),
           SLOT(standardOrientation()));
   m_actions.push_back(m_standardOrientationAction);
-  m_standardOrientationAction->setProperty("menu priority", -250);
+  m_standardOrientationAction->setProperty("menu priority", 170);
 
   m_scaleVolumeAction->setText(tr("Scale Cell &Volume"));
   connect(m_scaleVolumeAction, SIGNAL(triggered()), SLOT(scaleVolume()));
   m_actions.push_back(m_scaleVolumeAction);
-  m_scaleVolumeAction->setProperty("menu priority", -275);
+  m_scaleVolumeAction->setProperty("menu priority", 160);
 
   m_buildSupercellAction->setText(tr("Build &Supercell"));
   connect(m_buildSupercellAction, SIGNAL(triggered()), SLOT(buildSupercell()));
   m_actions.push_back(m_buildSupercellAction);
-  m_buildSupercellAction->setProperty("menu priority", -300);
+  m_buildSupercellAction->setProperty("menu priority", 150);
 
   m_niggliReduceAction->setText(tr("Reduce Cell (&Niggli)"));
   connect(m_niggliReduceAction, SIGNAL(triggered()), SLOT(niggliReduce()));
   m_actions.push_back(m_niggliReduceAction);
-  m_niggliReduceAction->setProperty("menu priority", -350);
+  m_niggliReduceAction->setProperty("menu priority", 140);
 
   updateActions();
 }

avogadro/qtplugins/spacegroup/CMakeLists.txt

@@ -0,0 +1,17 @@
+set(spacegroup_srcs
+  spacegroup.cpp
+)
+
+set(spacegroup_uis
+
+)
+
+avogadro_plugin(SpaceGroup
+  "Space group features for crystals."
+  ExtensionPlugin
+  spacegroup.h
+  SpaceGroup
+  "${spacegroup_srcs}"
+  "${spacegroup_uis}"
+)
+