Initial attempt to enforce unit cell molecules to have fractional cen…

…troids in the gamma point for occ cg

include/occ/crystal/crystal.h

@@ -5,9 +5,9 @@
 #include <occ/core/molecule.h>
 #include <occ/crystal/asymmetric_unit.h>
 #include <occ/crystal/hkl.h>
+#include <occ/crystal/site_mapping_table.h>
 #include <occ/crystal/spacegroup.h>
 #include <occ/crystal/unitcell.h>
-#include <occ/crystal/site_mapping_table.h>
 #include <vector>
 
 namespace occ::crystal {
@@ -488,6 +488,13 @@ class Crystal {
 
   /**
    * \brief Specify the behaviour for guessing/finding bonds.
+   *
+   * Controls whether bonds should be automatically determined based on
+   * atomic positions and covalent radii. When enabled, bonds are guessed
+   * using distance criteria. When disabled, no bonds are created.
+   *
+   * \param guess If true, enables automatic bond detection. If false,
+   * disables bond detection. Defaults to true.
    */
   void set_connectivity_criteria(bool guess = true);
 
@@ -501,7 +508,6 @@ class Crystal {
    *
    * \param c The crystal lattice from which to create the primitive
    * supercell.
-   *
    * \param hkl The lattice vector triplet (h, k, l) that defines
    * the primitive supercell.
    *
@@ -510,9 +516,57 @@ class Crystal {
    */
   static Crystal create_primitive_supercell(const Crystal &c, HKL hkl);
 
+  /**
+   * \brief Creates a mapping table for atomic sites in the crystal.
+   *
+   * Generates a table that maps atomic sites in the crystal to their
+   * equivalent positions under symmetry operations. This is useful for
+   * analyzing site equivalences and symmetry relationships.
+   *
+   * \return A SiteMappingTable object containing the atomic site mappings.
+   */
   SiteMappingTable atom_site_mapping_table() const;
+
+  /**
+   * \brief Creates a mapping table for molecular sites in the crystal.
+   *
+   * Generates a table that maps molecular sites in the crystal to their
+   * equivalent positions under symmetry operations. This helps identify
+   * molecular symmetry relationships and equivalent molecular positions.
+   *
+   * \return A SiteMappingTable object containing the molecular site mappings.
+   */
   SiteMappingTable molecule_site_mapping_table() const;
 
+  /**
+   * \brief Checks if molecules are enforced to have centroids in the gamma
+   * point.
+   *
+   * Returns whether unit cell molecules are currently being enforced to have
+   * their geometric centroids in the range [0,1) for each fractional
+   * coordinate.
+   *
+   * \return True if gamma point enforcement is enabled, false otherwise.
+   */
+  inline bool gamma_point_unit_cell_molecules() const {
+    return m_enforce_gamma_point;
+  }
+
+  /**
+   * \brief Controls enforcement of molecule centroids in the gamma point.
+   *
+   * Sets whether unit cell molecules should have their geometric centroids
+   * enforced to lie within the range [0,1) for each fractional coordinate.
+   * When enabled, molecules will be shifted to satisfy this constraint while
+   * maintaining their internal structure and connectivity.
+   *
+   * \param set If true, enforces gamma point centering. If false, allows
+   * molecules to have centroids outside the [0,1) range.
+   */
+  inline void set_gamma_point_unit_cell_molecules(bool set) {
+    m_enforce_gamma_point = set;
+  }
+
 private:
   AsymmetricUnit m_asymmetric_unit;
   SpaceGroup m_space_group;
@@ -523,6 +577,8 @@ class Crystal {
   void update_unit_cell_connectivity() const;
   void update_unit_cell_atoms() const;
 
+  bool m_enforce_gamma_point{false};
+
   mutable std::vector<typename PeriodicBondGraph::VertexDescriptor>
       m_bond_graph_vertices;
   mutable PeriodicBondGraph m_bond_graph;

include/occ/main/occ_cg.h

@@ -18,6 +18,7 @@ struct CGConfig {
   bool write_kmcpp_file{false};
   bool use_xtb{false};
   bool symmetric_solvent_contribution{false};
+  bool gamma_point_molecules{true};
   std::string xtb_solvation_model{"cpcmx"};
   bool list_solvents{false};
   bool crystal_is_atomic{false};

src/crystal/crystal.cpp

@@ -384,7 +384,9 @@ void Crystal::update_unit_cell_molecules() const {
     g.breadth_first_traversal_with_edge(v.first, visitor);
     uc_mol_idx++;
   }
+
   Mat3N cart_pos = to_cartesian(atoms.frac_pos + shifts.cast<double>());
+
   for (size_t i = 0; i < uc_mol_idx; i++) {
     auto idx = atom_indices[i];
 
@@ -401,6 +403,30 @@ void Crystal::update_unit_cell_molecules() const {
     IMat3N shift_hkl = shifts(Eigen::all, idx);
 
     occ::core::Molecule m(atoms.atomic_numbers(idx), cart_pos({0, 1, 2}, idx));
+
+    if (m_enforce_gamma_point) {
+      Vec3 frac_centroid = to_fractional(m.centroid());
+      IVec3 centroid_shift;
+      centroid_shift[0] = -std::floor(frac_centroid[0]);
+      centroid_shift[1] = -std::floor(frac_centroid[1]);
+      centroid_shift[2] = -std::floor(frac_centroid[2]);
+      occ::log::debug("Frac centroid before shift to gamma point: {}",
+                      frac_centroid.transpose());
+      occ::log::debug(
+          "Frac centroid after shift: {}",
+          (frac_centroid + centroid_shift.cast<double>()).transpose());
+
+      if (centroid_shift.squaredNorm() > 0) {
+        // Apply shift to atomic positions
+        Mat3N shifted_cart_pos = cart_pos({0, 1, 2}, idx);
+        shifted_cart_pos += to_cartesian(centroid_shift.cast<double>());
+        m = occ::core::Molecule(atoms.atomic_numbers(idx), shifted_cart_pos);
+
+        // Update unit cell shifts
+        shift_hkl.colwise() += centroid_shift;
+      }
+    }
+
     m.set_unit_cell_idx(Eigen::Map<const IVec>(idx.data(), idx.size()));
     m.set_unit_cell_shift(shift_hkl);
     m.set_asymmetric_unit_idx(atoms.asym_idx(idx));
@@ -717,5 +743,4 @@ SiteMappingTable Crystal::molecule_site_mapping_table() const {
   return SiteMappingTable::build_molecule_table(*this);
 }
 
-
 } // namespace occ::crystal

src/io/load_geometry.cpp

@@ -1,38 +1,33 @@
-#include <occ/io/load_geometry.h>
+#include <occ/io/cifparser.h>
 #include <occ/io/dftb_gen.h>
+#include <occ/io/load_geometry.h>
 #include <occ/io/xyz.h>
-#include <occ/io/cifparser.h>
 
 namespace occ::io {
 
 occ::crystal::Crystal load_crystal(const std::string &filename) {
-    if(CifParser::is_likely_cif_filename(filename)) {
-	occ::io::CifParser parser;
-	return parser.parse_crystal(filename).value();
-    }
-    else if(DftbGenFormat::is_likely_gen_filename(filename)) {
-	DftbGenFormat parser;
-	parser.parse(filename);
-	return parser.crystal().value();
-    }
-    else throw std::runtime_error(
-	fmt::format("Unknown filetype when reading crystal from '{}'", filename)
-    );
+  if (CifParser::is_likely_cif_filename(filename)) {
+    occ::io::CifParser parser;
+    return parser.parse_crystal(filename).value();
+  } else if (DftbGenFormat::is_likely_gen_filename(filename)) {
+    DftbGenFormat parser;
+    parser.parse(filename);
+    return parser.crystal().value();
+  } else
+    throw std::runtime_error(fmt::format(
+        "Unknown filetype when reading crystal from '{}'", filename));
 }
 
 occ::core::Molecule load_molecule(const std::string &filename) {
-    if(XyzFileReader::is_likely_xyz_filename(filename)) {
-	return molecule_from_xyz_file(filename);
-    }
-    else if(DftbGenFormat::is_likely_gen_filename(filename)) {
-	DftbGenFormat parser;
-	parser.parse(filename);
-	return parser.molecule().value();
-    }
-    else throw std::runtime_error(
-	fmt::format("Unknown filetype when reading molecule from '{}'", filename)
-    );
-
+  if (XyzFileReader::is_likely_xyz_filename(filename)) {
+    return molecule_from_xyz_file(filename);
+  } else if (DftbGenFormat::is_likely_gen_filename(filename)) {
+    DftbGenFormat parser;
+    parser.parse(filename);
+    return parser.molecule().value();
+  } else
+    throw std::runtime_error(fmt::format(
+        "Unknown filetype when reading molecule from '{}'", filename));
 }
 
-}
+} // namespace occ::io

src/main/occ_cg.cpp

@@ -110,7 +110,12 @@ CLI::App *add_cg_subcommand(CLI::App &app) {
       "--symmetric-solvent-contribution,--symmetric_solvent_contribution",
       config->symmetric_solvent_contribution,
       "Crystal growth interactions will have permutational symmetry (i.e. A->B "
-      "== B->A)");
+      "== B->A) (default: false)");
+  cg->add_flag(
+      "--gamma-point-molecules,--gamma_point_molecules",
+      config->gamma_point_molecules,
+      "Enforce that the resulting unit cell molecules (e.g. in the net file) "
+      "must have geometric centroids in the range [0,1) (default: true)");
   cg->add_option("--surface-energies", config->max_facets,
                  "Calculate surface energies and write .gmf morphology files");
   cg->add_flag("--list-available-solvents", config->list_solvents,
@@ -290,6 +295,8 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) {
   Crystal c_symm =
       occ::io::load_crystal(config.lattice_settings.crystal_filename);
 
+  c_symm.set_gamma_point_unit_cell_molecules(config.gamma_point_molecules);
+
   if (config.crystal_is_atomic) {
     c_symm.set_connectivity_criteria(false);
   }
@@ -307,7 +314,10 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) {
   if (config.use_xtb)
     opts.use_asymmetric_partition = false;
 
-  occ::log::info("Enforcing asymmetry via partitioning: {}", opts.use_asymmetric_partition);
+  occ::log::info("Enforcing asymmetry via partitioning:   {}",
+                 opts.use_asymmetric_partition);
+  occ::log::info("Enforcing unit cell molecules in gamma: {}",
+                 config.gamma_point_molecules);
 
   opts.wavefunction_choice =
       (config.wavefunction_choice == "gas" ? WavefunctionChoice::GasPhase
@@ -389,14 +399,6 @@ occ::cg::CrystalGrowthResult run_cg_impl(CGConfig const &config) {
 
 occ::cg::CrystalGrowthResult run_cg(CGConfig const &config) {
   occ::cg::CrystalGrowthResult result;
-  std::string basename =
-      fs::path(config.lattice_settings.crystal_filename).stem().string();
-  Crystal c_symm =
-      occ::io::load_crystal(config.lattice_settings.crystal_filename);
-
-  if (config.crystal_is_atomic) {
-    c_symm.set_connectivity_criteria(false);
-  }
 
   if (config.use_xtb) {
     result = run_cg_impl<driver::XTBCrystalGrowthCalculator>(config);