From 6eaad7e0974eb9a55b9b4f66175fecacc3f26dc7 Mon Sep 17 00:00:00 2001
From: Geoff Hutchison <geoff.hutchison@gmail.com>
Date: Tue, 17 Oct 2023 12:37:16 -0400
Subject: [PATCH] Make sure "fill unit cell" generates copies on the unit cell
 too

Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
---
 avogadro/core/spacegroups.cpp | 88 +++++++++++++++++++++++++++++++----
 1 file changed, 79 insertions(+), 9 deletions(-)

diff --git a/avogadro/core/spacegroups.cpp b/avogadro/core/spacegroups.cpp
index bf67444de8..f8a32a23e6 100644
--- a/avogadro/core/spacegroups.cpp
+++ b/avogadro/core/spacegroups.cpp
@@ -20,13 +20,9 @@
 
 namespace Avogadro::Core {
 
-SpaceGroups::SpaceGroups()
-{
-}
+SpaceGroups::SpaceGroups() {}
 
-SpaceGroups::~SpaceGroups()
-{
-}
+SpaceGroups::~SpaceGroups() {}
 
 unsigned short SpaceGroups::hallNumber(const std::string& spaceGroup)
 {
@@ -64,7 +60,6 @@ unsigned short SpaceGroups::hallNumber(const std::string& spaceGroup)
   return hall; // can't find anything
 }
 
-
 CrystalSystem SpaceGroups::crystalSystem(unsigned short hallNumber)
 {
   if (hallNumber == 1 || hallNumber == 2)
@@ -254,7 +249,7 @@ 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;
@@ -305,7 +300,82 @@ void SpaceGroups::fillUnitCell(Molecule& mol, unsigned short hallNumber,
       newAtom.setPosition3d(newCandidate);
     }
   }
+  // @todo make this optional
   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
+  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 +429,4 @@ const char* SpaceGroups::transformsString(unsigned short hallNumber)
     return "";
 }
 
-} // end Avogadro namespace
+} // namespace Avogadro::Core