Changeset 2390e6 for src/Actions/AtomAction/SaturateAction.cpp

Timestamp:

Apr 28, 2021, 10:02:49 PM (5 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

Candidate_v1.7.0, stable

Children:

e0e77e

Parents:

87c1cc

git-author:

Frederik Heber <frederik.heber@…> (04/17/21 09:41:33)

git-committer:

Frederik Heber <frederik.heber@…> (04/28/21 22:02:49)

Message:

SaturateAction may use outer shell information.

• when considering the outer shell, then we look at both occupied and unoccupied orbitals. The ideal polygon will be created for the number of all orbitals, but only the unoccupied will be allowed to be taken by either present bonds or saturation hydrogens.
• SphericalPointDistribution allows to occupy less than present points.
• I have tested this manually by adding oxygen, saturation, adding carbon, saturating again, then turning oxygen into nitrogen and saturating. Results were always exactly as hoped for.
• DOCU: Added note to action in userguide.

File:

• src/Actions/AtomAction/SaturateAction.cpp (modified) (2 diffs)

src/Actions/AtomAction/SaturateAction.cpp

@@ -86 +86 @@
     const BondList& ListOfBonds = _atom->getListOfBonds();
     SphericalPointDistribution PointSphere(typical_distance);
+    int num_orbitals_outer_shell;
+    int num_unoccupied_orbitals;
+    if (params.useOuterShell.get()) {
+      /*
+       * Number of valence electrons (getValence()) + number of unoccupied
+       * orbitals (getNoValenceOrbitals()) is the number of electrons able
+       * to fit into the outer shell. Hence, half the number gives the number
+       * of orbitals in the outer shell
+       */
+      num_orbitals_outer_shell =
+          (((int)_atom->getType()->getValence())+_atom->getType()->getNoValenceOrbitals())/2;
+      num_unoccupied_orbitals = _atom->getType()->getNoValenceOrbitals();
+    } else {
+      /**
+       * Here, same number of points on the polyhedra and available ones.
+       */
+      num_orbitals_outer_shell = num_unoccupied_orbitals = _atom->getType()->getNoValenceOrbitals();
+    }
     if (ListOfBonds.size() == 0) {
-      vacant_positions = PointSphere.getSimplePolygon(_atom->getType()->getNoValenceOrbitals());
+      vacant_positions = PointSphere.getSimplePolygon(num_orbitals_outer_shell);
+
+      // if less available then present in the ideal polygon, remove some
+      if (num_unoccupied_orbitals < vacant_positions.size())
+        vacant_positions.resize(num_unoccupied_orbitals);
+
       LOG(3, "DEBUG: Using ideal positions as " << vacant_positions);
     } else {
-      // get ideal polygon and currently occupied positions
-      const SphericalPointDistribution::Polygon_t ideal_positions =
-          PointSphere.getSimplePolygon(_atom->getType()->getNoValenceOrbitals());
-      LOG(3, "DEBUG: ideal positions are " << ideal_positions);
       SphericalPointDistribution::WeightedPolygon_t current_positions;
       for (BondList::const_iterator bonditer = ListOfBonds.begin();
@@ -104 +123 @@
       LOG(3, "DEBUG: current occupied positions are " << current_positions);
 
+      // get ideal polygon and currently occupied positions
+      const SphericalPointDistribution::Polygon_t ideal_positions =
+          PointSphere.getSimplePolygon(num_orbitals_outer_shell);
+      LOG(3, "DEBUG: ideal positions are " << ideal_positions);
+
       // find the best matching rotated polygon
       vacant_positions = PointSphere.getRemainingPoints(
           current_positions,
-          _atom->getType()->getNoValenceOrbitals());
+          num_orbitals_outer_shell,
+          num_unoccupied_orbitals);
       LOG(3, "DEBUG: Resulting vacant positions are " << vacant_positions);