metastable node detection and removal for BinaryStrategy

examples/TernaryExamples.ipynb

@@ -34,8 +34,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: total: 32.8 s\n",
-      "Wall time: 32.9 s\n"
+      "CPU times: total: 39.7 s\n",
+      "Wall time: 40.4 s\n"
      ]
     },
     {

pycalphad/mapping/plotting.py

@@ -420,7 +420,7 @@ def plot_ternary(strategy: TernaryStrategy, x: v.StateVariable = None, y: v.Stat
     if ax is None:
         fig, ax = plt.subplots(1, 1, subplot_kw={'projection': "triangular"})
 
-    plot_binary(strategy, x, y, ax, tielines=tielines, label_node=label_nodes, legend_generator=legend_generator, tieline_color=tieline_color, tie_triangle_color=tie_triangle_color, *args, **kwargs)
+    plot_binary(strategy, x, y, ax, tielines=tielines, label_nodes=label_nodes, legend_generator=legend_generator, tieline_color=tieline_color, tie_triangle_color=tie_triangle_color, *args, **kwargs)
     ax.set_xlim([0,1])
     ax.set_ylim([0,1])
 

pycalphad/mapping/strategy/binary_strategy.py

@@ -207,3 +207,17 @@ def _determine_start_direction(self, node: Node, exit_point: Point, proposed_dir
                 return exit_point, axis_var, d, av_delta
 
         return None
+
+    def _process_new_node(self, zpf_line: ZPFLine, new_node: Node):
+        """
+        Post global eq check after finding node to ensure it's not a metastable node
+        """
+        _log.info("Checking if new node is metastable")
+        cs_result = zeq._find_global_min_cs(new_node, system_info=self.system_info, pdens=self.GLOBAL_MIN_PDENS, tol=self.GLOBAL_MIN_TOL, num_candidates=self.GLOBAL_MIN_NUM_CANDIDATES)
+        if cs_result is None:
+            _log.info("Global eq check on new node passed")
+            super()._process_new_node(zpf_line, new_node)
+        else:
+            _log.info("Global eq check failed. New node is metastable. Removing current zpf line.")
+            self.zpf_lines.pop(-1)
+

pycalphad/mapping/strategy/ternary_strategy.py

@@ -295,7 +295,7 @@ def _process_new_node(self, zpf_line: ZPFLine, new_node: Node):
         # Do a global equilibrium check before attempting to add the new node
         global_eq_check, test_point = self._check_full_global_equilibrium(new_node, add_global_point_if_false = False)
         if test_point is None:
-            _log.info("Global eq check could not be determined")
+            _log.info("Global eq check on new node could not be determined")
 
             zpf_line.current_delta *= self.DELTA_SCALE
             if zpf_line.current_delta / self.axis_delta[zpf_line.axis_var] < self.MIN_DELTA_RATIO:
@@ -310,12 +310,12 @@ def _process_new_node(self, zpf_line: ZPFLine, new_node: Node):
 
         else:
             if global_eq_check:
-                _log.info("Global eq check passed")
+                _log.info("Global eq check on new node passed")
                 return super()._process_new_node(zpf_line, new_node)
             else:
                 zpf_line_phases = zpf_line.stable_phases_with_multiplicity
                 test_node_phases = test_point.stable_phases_with_multiplicity
-                _log.info(f"Global eq check failed. Comparing zpf line {zpf_line_phases} with node phases {test_node_phases}")
+                _log.info(f"Global eq check failed. New node is metastable. Comparing zpf line {zpf_line_phases} with node phases {test_node_phases}")
                 if len(set(test_node_phases) - set(zpf_line_phases)) == 1:
                     _log.info("Node can be added as a proper node")
                     new_node = self._create_node_from_point(test_point, zpf_line.points[-1], None, None)
@@ -324,7 +324,7 @@ def _process_new_node(self, zpf_line: ZPFLine, new_node: Node):
                     zpf_line.current_delta *= self.DELTA_SCALE
                     if zpf_line.current_delta / self.axis_delta[zpf_line.axis_var] < self.MIN_DELTA_RATIO:
                         zpf_line.status = ZPFState.FAILED
-                        _log.info("Node is added as a starting point")
+                        _log.info("Node is added as a starting point. Removing current zpf line.")
                         new_node = self._create_node_from_point(test_point, None, None, None)
                         if not self.node_queue.add_node(new_node):
                             _log.info(f"Node {new_node.fixed_phases}, {new_node.free_phases} has already been added")

pycalphad/mapping/zpf_equilibrium.py

@@ -70,15 +70,15 @@ def update_equilibrium_with_new_conditions(point: Point, new_conditions: dict[v.
     new_point = Point(new_conditions, np.array(results.chemical_potentials), [cs for cs in comp_sets if cs.fixed], [cs for cs in comp_sets if not cs.fixed])
     return new_point, orig_cs
 
-def find_global_min_point(point: Point, system_info: dict, pdens = 500, tol = 1e-5, num_candidates = 1):
+def _find_global_min_cs(point: Point, system_info: dict, pdens = 500, tol = 1e-5, num_candidates = 1):
     """
+    Finds potential composition set for global min check
+
     For each possible phase:
         1. Sample DOF and find CS that maximizes driving force
         2. Create a DormantPhase with CS and compute driving force with potentials at equilibrium
             Or check the top N CS that maximizes driving force and compute driving force and take max
-        3. If driving force is positive, then new phase is stable
-        4. Check that the new CS doesn"t match with a currently stable CS
-        4. Hope that this works on miscibility gaps
+        3. If driving force is positive, then new phase may be stable in when attempting to find a new node
 
     Parameters
     ----------
@@ -98,11 +98,11 @@ def find_global_min_point(point: Point, system_info: dict, pdens = 500, tol = 1e
 
     Returns
     -------
-    new_point : Point with added composition set
+    (cs, dG) : (potential new composition set, driving force of composition set)
 
     or
 
-    None : if point is global min or global min was unable to be found
+    None : if no composition set was found below the current chemical potential hyperplane
     """
     dbf = system_info["dbf"]
     comps = system_info["comps"]
@@ -148,6 +148,49 @@ def find_global_min_point(point: Point, system_info: dict, pdens = 500, tol = 1e
     if dG < tol:
         return None
     else:
+        return cs, dG
+
+def find_global_min_point(point: Point, system_info: dict, pdens = 500, tol = 1e-5, num_candidates = 1):
+    """
+    Checks global min on current point and attempts to find a new point 
+    with a new composition set if current point is not global min
+
+    1. Find potential new composition set for global min
+    2. If no composition set is found, then return
+    3. If composition set is found, check that the new CS doesn't match
+       with a currently stable CS
+    4. Hope that this works on miscibility gaps
+
+    Parameters
+    ----------
+    point : Point
+        Point to check global min
+    system_info : dict
+        Dictionary containing information for pycalphad.calculate
+    pdens : int
+        Sampling density
+    tol : float
+        Tolerance for whether a new CS is considered stable
+    num_candidates : int
+        Number of candidate CS to check driving force
+        To avoid redundant calculations, this will only check driving force
+        for unique phases. So setting this to a high number will not significantly
+        affect performance
+
+    Returns
+    -------
+    new_point : Point with added composition set
+
+    or
+
+    None : if point is global min or global min was unable to be found
+    """
+    min_cs_result = _find_global_min_cs(point, system_info, pdens, tol, num_candidates)
+    if min_cs_result is None:
+        return None
+    else:
+        cs, dG = min_cs_result
+
         _log.info(f'Global min potentially detected. {point.stable_phases} + {cs.phase_record.phase_name} with dG = {dG}')
         if _detect_degenerate_phase(point, cs):
             new_point = Point(point.global_conditions, point.chemical_potentials, point.fixed_composition_sets, point.free_composition_sets)