Merge f394830 into 1abbecd

pymatgen/analysis/diffusion/aimd/clustering.py

@@ -116,7 +116,7 @@ def should_stop(self, old_centroids, centroids, iterations):
             iterations: Number of iterations thus far.
         """
         if iterations > self.max_iterations:
-            warnings.warn("Max iterations %d reached!" % self.max_iterations)
+            warnings.warn(f"Max iterations {self.max_iterations} reached!")
             return True
         if old_centroids is None:
             return False
@@ -192,7 +192,7 @@ def should_stop(self, old_centroids, centroids, iterations):
             iterations: Number of iterations thus far.
         """
         if iterations > self.max_iterations:
-            warnings.warn("Max iterations %d reached!" % self.max_iterations)
+            warnings.warn("Max iterations {self.max_iterations} reached!")
             return True
         if old_centroids is None:
             return False

pymatgen/analysis/diffusion/aimd/pathway.py

@@ -182,9 +182,9 @@ def generate_stable_sites(self, p_ratio=0.25, d_cutoff=1.0):
 
         if nc < nions:
             raise ValueError(
-                "The number of clusters ({}) is smaller than that of "
-                "mobile ions ({})! Please try to decrease either "
-                "'p_ratio' or 'd_cut' values!".format(nc, nions)
+                f"The number of clusters ({nc}) is smaller than that of "
+                f"mobile ions ({nions})! Please try to decrease either "
+                "'p_ratio' or 'd_cut' values!"
             )
 
         # For each low-energy site (cluster centroid), its coordinates are obtained
@@ -227,7 +227,7 @@ def to_chgcar(self, filename="CHGCAR.vasp"):
         """
 
         count = 1
-        VolinAu = self.structure.lattice.volume / 0.5291772083 ** 3
+        #VolinAu = self.structure.lattice.volume / 0.5291772083 ** 3
         symbols = self.structure.symbol_set
         natoms = [str(int(self.structure.composition[symbol])) for symbol in symbols]
         init_fcoords = np.array(self.structure.frac_coords)
@@ -237,21 +237,22 @@ def to_chgcar(self, filename="CHGCAR.vasp"):
             f.write(" 1.00 \n")
 
             for i in range(3):
-                f.write(" {0} {1} {2} \n".format(*self.structure.lattice.matrix[i, :]))
+                v = self.structure.lattice.matrix[i, :]
+                f.write(f" {v[0]} {v[1]} {v[2]} \n")
 
             f.write(" " + " ".join(symbols) + "\n")
             f.write(" " + " ".join(natoms) + "\n")
             f.write("direct\n")
             for fcoord in init_fcoords:
-                f.write(" {0:.8f}  {1:.8f}  {2:.8f} \n".format(*fcoord))
+                f.write(" {fcoord[0]:.8f}  {fcoord[1]:.8f}  {fcoord[2]:.8f} \n")
 
             f.write(" \n")
             f.write(" {0} {1} {2} \n".format(*self.lens))
 
             for i in range(self.lens[2]):
                 for j in range(self.lens[1]):
                     for k in range(self.lens[0]):
-                        f.write(" {0:.10e} ".format(self.Pr[k, j, i] * VolinAu))
+                        f.write(" {self.Pr[k, j, i] * VolinAu:.10e} ")
                         if count % 5 == 0:
                             f.write("\n")
                         count += 1

pymatgen/analysis/diffusion/aimd/rdf.py

@@ -276,7 +276,7 @@ def export_rdf(self, filename: str):
             f.write("\n")
 
             for r, gr in zip(self.interval, self.rdf):
-                f.write(delimiter.join(["%s" % v for v in [r, gr]]))
+                f.write(delimiter.join([f"{v}" for v in [r, gr]]))
                 f.write("\n")
 
 

pymatgen/analysis/diffusion/aimd/van_hove.py

@@ -542,7 +542,7 @@ def plot_rdf_evolution(
         """
         if df is None:
             df = self.get_df(func=EvolutionAnalyzer.rdf, pair=pair)
-        x_label, cb_label = "$r$ ({}-{}) ($\\rm\AA$)".format(*pair), "$g(r)$"
+        x_label, cb_label = f"$r$ ({pair[0]}-{pair[1]}) ($\\rm\AA$)", "$g(r)$"
         p = self.plot_evolution_from_data(df=df, x_label=x_label, cb_label=cb_label, cmap=cmap)
 
         return p
@@ -569,7 +569,7 @@ def plot_atomic_evolution(
         if df is None:
             df = self.get_df(func=EvolutionAnalyzer.atom_dist, specie=specie, direction=direction)
         x_label, cb_label = (
-            "Atomic distribution along {} ".format(direction),
+            f"Atomic distribution along {direction} ",
             "Probability",
         )
         p = self.plot_evolution_from_data(df=df, x_label=x_label, cb_label=cb_label, cmap=cmap)

pymatgen/analysis/diffusion/analyzer.py

@@ -463,7 +463,7 @@ def get_framework_rms_plot(self, plt=None, granularity=200, matching_s=None):
         plt.plot(plot_dt, rms[:, 0], label="RMS")
         plt.plot(plot_dt, rms[:, 1], label="max")
         plt.legend(loc="best")
-        plt.xlabel("Timestep ({})".format(unit))
+        plt.xlabel(f"Timestep ({unit})")
         plt.ylabel("normalized distance")
         plt.tight_layout()
         return plt
@@ -499,7 +499,7 @@ def get_msd_plot(self, plt=None, mode="specie"):
         elif mode == "sites":
             for i, site in enumerate(self.structure):
                 sd = self.sq_disp_ions[i, :]
-                plt.plot(plot_dt, sd, label="%s - %d" % (site.specie.__str__(), i))
+                plt.plot(plot_dt, sd, label=f"{site.specie.__str__()} - {i}")
             plt.legend(loc=2, prop={"size": 20})
         elif mode == "mscd":
             plt.plot(plot_dt, self.mscd, "r")
@@ -512,7 +512,7 @@ def get_msd_plot(self, plt=None, mode="specie"):
             plt.plot(plot_dt, self.msd_components[:, 2], "b")
             plt.legend(["Overall", "a", "b", "c"], loc=2, prop={"size": 20})
 
-        plt.xlabel("Timestep ({})".format(unit))
+        plt.xlabel(f"Timestep ({unit})")
         if mode == "mscd":
             plt.ylabel("MSCD ($\\AA^2$)")
         else:
@@ -551,7 +551,7 @@ def export_msdt(self, filename):
             f.write(delimiter.join(["t", "MSD", "MSD_a", "MSD_b", "MSD_c", "MSCD"]))
             f.write("\n")
             for dt, msd, msdc, mscd in zip(self.dt, self.msd, self.msd_components, self.mscd):
-                f.write(delimiter.join(["%s" % v for v in [dt, msd] + list(msdc) + [mscd]]))
+                f.write(delimiter.join([f"{v}" for v in [dt, msd] + list(msdc) + [mscd]]))
                 f.write("\n")
 
     @classmethod
@@ -918,7 +918,7 @@ def get_arrhenius_plot(temps, diffusivities, diffusivity_errors=None, **kwargs):
     plt.text(
         0.6,
         0.85,
-        "E$_a$ = {:.0f} meV".format(Ea * 1000),
+        f"E$_a$ = {Ea * 1000:.0f} meV",
         fontsize=30,
         transform=plt.axes().transAxes,
     )

pymatgen/analysis/diffusion/neb/full_path_mapper.py

@@ -358,7 +358,7 @@ def assign_cost_to_graph(self, cost_keys=["hop_distance"]):
             cost_val = np.prod([v[ik] for ik in cost_keys])
             self.add_data_to_similar_edges(k, {"cost": cost_val})
 
-    def get_path(self, max_val=100000):
+    def get_path(self, max_val=100000, flip_hops=True):
         """
         obtain a pathway through the material using hops that are in the current graph
         Basic idea:
@@ -367,6 +367,10 @@ def get_path(self, max_val=100000):
             or any other neighboring UC not containing p1.
         Args:
             max_val: Filter the graph by a cost
+            flip_hops: If true, hops in paths returned will be flipped so
+                isites and esites match to form a coherent path.
+                If false, hops will retain their original orientation
+                from the migration graph.
         Returns:
             Generator for List of Dicts:
             Each dict contains the information of a hop
@@ -421,7 +425,10 @@ def get_path(self, max_val=100000):
                         found_ += 1
                 if found_ != 1:
                     raise RuntimeError("More than one edge matched in original graph.")
-            yield u, path_hops
+            if flip_hops is True:  # flip hops in path to form coherent pathway
+                yield u, order_path(path_hops, u)
+            else:
+                yield u, path_hops
 
     def get_summary_dict(self, added_keys: List[str] = None) -> dict:
         """
@@ -635,12 +642,7 @@ def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None):
             isym = self.symm_structure.wyckoff_symbols[migration_hop.iindex]
             esym = self.symm_structure.wyckoff_symbols[migration_hop.eindex]
             mask_out.write_file(
-                "{}_{}_{}_tot({:0.2f}).vasp".format(
-                    mask_file_seedname,
-                    isym,
-                    esym,
-                    mask_out.data[self.potential_data_key].sum(),
-                )
+                f"{mask_file_seedname}_{isym}_{esym}_tot({mask_out.data[self.potential_data_key].sum():0.2f}).vasp"
             )
 
         return (
@@ -768,6 +770,62 @@ def get_hop_site_sequence(hop_list: List[Dict], start_u: Union[int, str], key: s
     return site_seq
 
 
+def order_path(hop_list: List[Dict], start_u: Union[int, str]) -> List[Dict]:
+    """
+    Takes a list of hop dictionaries and flips hops (switches isite and esite)
+    as needed to form a coherent path / sequence of sites according to
+    get_hop_site_sequence().
+    For example if hop_list = [{iindex:0, eindex:1, etc.}, {iindex:0, eindex:1, etc.}]
+    then the output is [{iindex:0, eindex:1, etc.}, {iindex:1, eindex:0, etc.}] so that
+    the following hop iindex matches the previous hop's eindex.
+    Args:
+        hop_list: a list of the data on a sequence of hops
+        start_u: the site index of the starting sites
+    Returns:
+        a list of the data on a sequence of hops with hops in coherent orientation
+    """
+    seq = get_hop_site_sequence(hop_list, start_u)
+
+    ordered_path = []
+    for n, hop in zip(seq[:-1], hop_list):
+        if n == hop["iindex"]:  # don't flip hop
+            ordered_path.append(hop)
+        else:
+            # create flipped hop
+            fh = MigrationHop(
+                isite=hop["hop"].esite,
+                esite=hop["hop"].isite,
+                symm_structure=hop["hop"].symm_structure,
+                host_symm_struct=None,
+                symprec=hop["hop"].symprec,
+            )
+            # must manually set iindex and eindex
+            fh.iindex = hop["hop"].eindex
+            fh.eindex = hop["hop"].iindex
+            fhd = {
+                "to_jimage": tuple(-1 * i for i in hop["to_jimage"]),
+                "ipos": fh.isite.frac_coords,
+                "epos": fh.esite.frac_coords,
+                "ipos_cart": fh.isite.coords,
+                "epos_cart": fh.esite.coords,
+                "hop": fh,
+                "hop_label": hop["hop_label"],
+                "iindex": hop["eindex"],
+                "eindex": hop["iindex"],
+                "hop_distance": fh.length,
+            }
+            # flip any data that is in a list to match flipped hop orientation
+            for k in hop.keys():
+                if k not in fhd.keys():
+                    if isinstance(hop[k], list):
+                        fhd[k] = hop[k][::-1]
+                    else:
+                        fhd[k] = hop[k]
+            ordered_path.append(fhd)
+
+    return ordered_path
+
+
 """
 Note the current pathway algorithm no longer needs supercells but the following
 functions might still be useful for other applications

pymatgen/analysis/diffusion/neb/pathfinder.py

@@ -354,23 +354,13 @@ def __init__(
 
     def __repr__(self):
         return (
-            "Path of %.4f A from %s [%.3f, %.3f, %.3f] "
-            "(ind: %d, Wyckoff: %s) to %s [%.3f, %.3f, %.3f] (ind: %d, Wyckoff: %s)"
-            % (
-                self.length,
-                self.isite.specie,
-                self.isite.frac_coords[0],
-                self.isite.frac_coords[1],
-                self.isite.frac_coords[2],
-                self.iindex,
-                self.symm_structure.wyckoff_symbols[self.iindex],
-                self.esite.specie,
-                self.esite.frac_coords[0],
-                self.esite.frac_coords[1],
-                self.esite.frac_coords[2],
-                self.eindex,
-                self.symm_structure.wyckoff_symbols[self.eindex],
-            )
+            f"Path of {self.length:.4f} A from {self.isite.specie} "
+            f"[{self.isite.frac_coords[0]:.3f}, {self.isite.frac_coords[1]:.3f}, "
+            f"{self.isite.frac_coords[2]:.3f}] (ind: {self.iindex}, "
+            f"Wyckoff: {self.symm_structure.wyckoff_symbols[self.iindex]}) "
+            f"to {self.esite.specie} [{self.esite.frac_coords[0]:.3f}, "
+            f"{self.esite.frac_coords[1]:.3f}, {self.esite.frac_coords[2]:.3f}] "
+            f"(ind: {self.eindex}, Wyckoff: {self.symm_structure.wyckoff_symbols[self.eindex]})"
         )
 
     @property

pymatgen/analysis/diffusion/neb/tests/test_full_path_mapper.py

@@ -14,6 +14,7 @@
     MigrationGraph,
     MigrationHop,
     get_hop_site_sequence,
+    order_path,
 )
 from pymatgen.analysis.diffusion.neb.periodic_dijkstra import _get_adjacency_with_images
 
@@ -163,24 +164,24 @@ def test_periodic_dijkstra(self):
 
     def test_get_path(self):
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
-        paths = [*self.fpm_li.get_path()]
+        paths = [*self.fpm_li.get_path(flip_hops=False)]
         p_strings = {"->".join(map(str, get_hop_site_sequence(ipath, start_u=u))) for u, ipath in paths}
         self.assertIn("5->7->5", p_strings)
         # convert each pathway to a string representation
-        paths = [*self.fpm_li.get_path(max_val=2.0)]
+        paths = [*self.fpm_li.get_path(max_val=2.0, flip_hops=False)]
         p_strings = {"->".join(map(str, get_hop_site_sequence(ipath, start_u=u))) for u, ipath in paths}
 
         # After checking trimming the graph more hops are needed for the same path
         self.assertIn("5->3->7->2->5", p_strings)
 
         self.fpm_mg.assign_cost_to_graph()  # use 'hop_distance'
-        paths = [*self.fpm_mg.get_path()]
+        paths = [*self.fpm_mg.get_path(flip_hops=False)]
         p_strings = {"->".join(map(str, get_hop_site_sequence(ipath, start_u=u))) for u, ipath in paths}
         self.assertIn("1->0->1", p_strings)
 
     def test_get_key_in_path(self):
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
-        paths = [*self.fpm_li.get_path()]
+        paths = [*self.fpm_li.get_path(flip_hops=False)]
         hop_seq_info = [get_hop_site_sequence(ipath, start_u=u, key="hop_distance") for u, ipath in paths]
 
         hop_distances = {}
@@ -200,6 +201,27 @@ def test_not_matching_first(self):
         for u, v, d in fpm_lmo.m_graph.graph.edges(data=True):
             self.assertIn(d["hop"].eindex, {0, 1})
 
+    def test_order_path(self):
+        # add list data to migration graph - to test if list data is flipped
+        for n, hop_d in self.fpm_li.unique_hops.items():
+            data = {"data": [hop_d["iindex"], hop_d["eindex"]]}
+            self.fpm_li.add_data_to_similar_edges(n, data, hop_d["hop"])
+
+        self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
+
+        for n, hop_list in [*self.fpm_li.get_path(flip_hops=False)]:
+            ordered = order_path(hop_list, n)
+
+            # check if isites and esites are oriented to form a coherent path
+            for h1, h2 in zip(ordered[:-1], ordered[1:]):
+                self.assertEqual(h1["eindex"], h2["iindex"])
+
+            # if hop was flipped, check list data was flipped too
+            for h, ord_h in zip(hop_list, ordered):
+                if h["iindex"] != ord_h["iindex"]:  # check only if hop was flipped
+                    self.assertEqual(h["data"][0], ord_h["data"][-1])
+                    self.assertEqual(h["data"][-1], ord_h["data"][0])
+
 
 class ChargeBarrierGraphTest(unittest.TestCase):
     def setUp(self):