Merge pull request #1375 from ghutchis/fill-unit-cell

Make sure "fill unit cell" generates copies on the unit cell too

avogadro/core/spacegroups.cpp

@@ -20,51 +20,49 @@
 
 namespace Avogadro::Core {
 
-SpaceGroups::SpaceGroups()
-{
-}
+SpaceGroups::SpaceGroups() {}
 
-SpaceGroups::~SpaceGroups()
-{
-}
+SpaceGroups::~SpaceGroups() {}
 
 unsigned short SpaceGroups::hallNumber(const std::string& spaceGroup)
 {
   unsigned short hall = 0; // can't find anything
-  const unsigned short hall_count = 530;
+  const unsigned short hall_count = 531; // 530 but first one is empty
+  // some files use " instead of = for the space group symbol
+  std::string sg = spaceGroup;
+  std::replace(sg.begin(), sg.end(), '"', '=');
 
   // space_group_hall_symbol
   for (unsigned short i = 0; i < hall_count; ++i) {
-    if (spaceGroup == space_group_hall_symbol[i]) {
+    if (sg == space_group_hall_symbol[i]) {
       return i; // found a match
     }
   }
 
   // space_group_international
   for (unsigned short i = 0; i < hall_count; ++i) {
-    if (spaceGroup == space_group_international[i]) {
+    if (sg == space_group_international[i]) {
       return i; // found a match
     }
   }
 
   // space_group_international_short
   for (unsigned short i = 0; i < hall_count; ++i) {
-    if (spaceGroup == space_group_international_short[i]) {
+    if (sg == space_group_international_short[i]) {
       return i; // found a match
     }
   }
 
   // space_group_international_full
   for (unsigned short i = 0; i < hall_count; ++i) {
-    if (spaceGroup == space_group_international_full[i]) {
+    if (sg == space_group_international_full[i]) {
       return i; // found a match
     }
   }
 
   return hall; // can't find anything
 }
 
-
 CrystalSystem SpaceGroups::crystalSystem(unsigned short hallNumber)
 {
   if (hallNumber == 1 || hallNumber == 2)
@@ -254,14 +252,14 @@ Array<Vector3> SpaceGroups::getTransforms(unsigned short hallNumber,
   // These transforms are separated by spaces
   std::vector<std::string> transforms = split(transformsStr, ' ');
 
-  for (auto & transform : transforms)
+  for (auto& transform : transforms)
     ret.push_back(getSingleTransform(transform, v));
 
   return ret;
 }
 
 void SpaceGroups::fillUnitCell(Molecule& mol, unsigned short hallNumber,
-                               double cartTol)
+                               double cartTol, bool wrapToCell, bool allCopies)
 {
   if (!mol.unitCell())
     return;
@@ -305,7 +303,86 @@ void SpaceGroups::fillUnitCell(Molecule& mol, unsigned short hallNumber,
       newAtom.setPosition3d(newCandidate);
     }
   }
-  CrystalTools::wrapAtomsToUnitCell(mol);
+
+  if (wrapToCell)
+    CrystalTools::wrapAtomsToUnitCell(mol);
+
+  // Now we need to generate any copies on the unit boundary
+  // We need to loop through all the atoms again
+  // if a fractional coordinate contains 0.0, we need to generate a copy
+  // of the atom at 1.0
+  if (!allCopies)
+    return;
+
+  atomicNumbers = mol.atomicNumbers();
+  positions = mol.atomPositions3d();
+  numAtoms = mol.atomCount();
+  for (Index i = 0; i < numAtoms; ++i) {
+    unsigned char atomicNum = atomicNumbers[i];
+    Vector3 pos = uc->toFractional(positions[i]);
+
+    Array<Vector3> newAtoms;
+    // We need to check each coordinate to see if it is 0.0 or 1.0
+    // .. if so, we need to generate the symmetric copy
+    for (Index j = 0; j < 3; ++j) {
+      Vector3 newPos = pos;
+      if (fabs(pos[j]) < 0.001) {
+        newPos[j] = 1.0;
+        newAtoms.push_back(newPos);
+      }
+      if (fabs(pos[j] - 1.0) < 0.001) {
+        newPos[j] = 0.0;
+        newAtoms.push_back(newPos);
+      }
+
+      for (Index k = j + 1; k < 3; ++k) {
+        if (fabs(pos[k]) < 0.001) {
+          newPos[k] = 1.0;
+          newAtoms.push_back(newPos);
+          newPos[k] = 0.0; // revert for the other coords
+        }
+        if (fabs(pos[k] - 1.0) < 0.001) {
+          newPos[k] = 0.0;
+          newAtoms.push_back(newPos);
+          newPos[k] = 1.0; // revert
+        }
+      }
+    }
+    // finally, check if all coordinates are 0.0
+    if (fabs(pos[0]) < 0.001 && fabs(pos[1]) < 0.001 && fabs(pos[2]) < 0.001)
+      newAtoms.push_back(Vector3(1.0, 1.0, 1.0));
+    // or 1.0 .. unlikely, but possible
+    if (fabs(pos[0] - 1.0) < 0.001 && fabs(pos[1] - 1.0) < 0.001 &&
+        fabs(pos[2] - 1.0) < 0.001)
+      newAtoms.push_back(Vector3(0.0, 0.0, 0.0));
+
+    for (Index j = 0; 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 = (mol.atomPosition3d(k) - newCandidate).norm();
+
+        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);
+    }
+  }
 }
 
 void SpaceGroups::reduceToAsymmetricUnit(Molecule& mol,
@@ -359,4 +436,4 @@ const char* SpaceGroups::transformsString(unsigned short hallNumber)
     return "";
 }
 
-} // end Avogadro namespace
+} // namespace Avogadro::Core

avogadro/core/spacegroups.h

@@ -42,7 +42,8 @@ class AVOGADROCORE_EXPORT SpaceGroups
   ~SpaceGroups();
 
   /**
-   * @return The hall number of the matching space group string or 0 if not found
+   * @return The hall number of the matching space group string or 0 if not
+   * found
    */
   static unsigned short hallNumber(const std::string& spaceGroup);
 
@@ -65,8 +66,8 @@ class AVOGADROCORE_EXPORT SpaceGroups
   static const char* schoenflies(unsigned short hallNumber);
 
   /**
-   * @return the Hall symbol for a given hall number. '=' is used instead of '"'.
-   * If an invalid hall number is given, an empty string will be returned.
+   * @return the Hall symbol for a given hall number. '=' is used instead of
+   * '"'. If an invalid hall number is given, an empty string will be returned.
    */
   static const char* hallSymbol(unsigned short hallNumber);
 
@@ -118,7 +119,8 @@ class AVOGADROCORE_EXPORT SpaceGroups
    * distance, the new atom will not be placed there.
    */
   static void fillUnitCell(Molecule& mol, unsigned short hallNumber,
-                           double cartTol = 1e-5);
+                           double cartTol = 1e-5, bool wrapToCell = true,
+                           bool allCopies = true);
 
   /**
    * Reduce a cell to its asymmetric unit.
@@ -137,7 +139,7 @@ class AVOGADROCORE_EXPORT SpaceGroups
   static const char* transformsString(unsigned short hallNumber);
 };
 
-} // end Core namespace
-} // end Avogadro namespace
+} // namespace Core
+} // namespace Avogadro
 
 #endif // AVOGADRO_CORE_SPACE_GROUPS_H