Add plot_lruj.py example

abipy/data/lruj_data/lruj.out

@@ -0,0 +1,83 @@
+
+.Version 9.11.6 of LRUJ 
+.(MPI version, prepared for a x86_64_linux_intel18.0 computer) 
+
+.Copyright (C) 1998-2022 ABINIT group . 
+ LRUJ comes with ABSOLUTELY NO WARRANTY.
+ It is free software, and you are welcome to redistribute it
+ under certain conditions (GNU General Public License,
+ see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt).
+
+ ABINIT is a project of the Universite Catholique de Louvain,
+ Corning Inc. and other collaborators, see ~abinit/doc/developers/contributors.txt .
+ Please read https://docs.abinit.org/theory/acknowledgments for suggested
+ acknowledgments of the ABINIT effort.
+ For more information, see https://www.abinit.org .
+
+.Starting date : Wed 29 Nov 2023.
+- ( at 15h33 )
+
+ Number of perturbations detected:  6
+ Including unperturbed state, we have  7 data points.
+ Hunds J determination, implemented by L. MacEnulty August 2021
+ Maximum degree of polynomials analyzed:  5
+ NOTE: Unlike the ujdet utility, lruj treats the 
+ response functions as scalars, not matrices!
+ See lruj tutorial for more information.
+
+
+*************************************************************************************************
+**************************************  Linear Response J  **************************************
+
+ Total number of atoms:    4
+ Index of perturbed atom:    1
+ Value of macro_uj:     4
+ Value of dmatpuopt:     3
+ Mixing constant factored out of Chi0:   0.450
+ Percentage of AE orbital within the PAW sphere of perturbed subspace:     90.46723%
+
+  Perturbations         Magnetizations
+ --------------- -----------------------------
+    beta [eV]     Unscreened      Screened
+ --------------- -----------------------------
+  -0.1500000676   1.3322568566   1.3664401054
+  -0.1000000451   1.3250447778   1.3485578255
+  -0.0500000225   1.3176770671   1.3298301923
+   0.0000000000   1.3101482092   1.3101482092
+   0.0500000225   1.3024514902   1.2893790961
+   0.1000000451   1.2945812190   1.2673571836
+   0.1500000676   1.2865306151   1.2438737113
+                                                                       RMS Errors
+                                                         ---------------------------------------
+ Regression   Chi0 [eV^-1]   Chi [eV^-1]      J [eV]    | Chi0 [eV^-1]  Chi [eV^-1]     J [eV]
+--------------------------------------------------------|---------------------------------------
+ Linear:        -0.3386211    -0.4075366     0.4993858  |    0.0006971    0.0020836    0.1364740
+ Quadratic:     -0.3386211    -0.4075366     0.4993858  |    0.0000137    0.0001529    0.0322580
+ Cubic:         -0.3383005    -0.4039797     0.4805799  |    0.0000002    0.0000120    0.0086817
+ Degree 4 :     -0.3383005    -0.4039797     0.4805799  |    0.0000001    0.0000009    0.0025307
+ Degree 5 :     -0.3383000    -0.4040198     0.4808302  |    0.0000001    0.0000000    0.0011046
+*************************************************************************************************
+*************************************************************************************************
+
+Linear Response UJ (LRUJ) program complete. Live long and prosper. ~LMac
+
+
+--- !LRUJ_Abipy_Plots
+natom: 4
+ndata: 7
+pawujat: 1
+macro_uj: 4
+diem_token: diemixmag
+diem:   4.50000000E-01
+chi0_coefficients_degree1: [  1.30981289E+00,  -1.52379516E-01, ]
+chi_coefficients_degree1: [  1.30794090E+00,  -4.07536646E-01, ]
+chi0_coefficients_degree2: [  1.31014818E+00,  -1.52379516E-01,  -3.35292571E-02, ]
+chi_coefficients_degree2: [  1.31016241E+00,  -4.07536646E-01,  -2.22150364E-01, ]
+chi0_coefficients_degree3: [  1.31014818E+00,  -1.52235216E-01,  -3.35292571E-02,  -8.24572223E-03, ]
+chi_coefficients_degree3: [  1.31016241E+00,  -4.03979653E-01,  -2.22150364E-01,  -2.03256604E-01, ]
+chi0_coefficients_degree4: [  1.31014812E+00,  -1.52235216E-01,  -3.35038810E-02,  -8.24572223E-03,  -1.06049292E-03, ]
+chi_coefficients_degree4: [  1.31014819E+00,  -4.03979653E-01,  -2.16856067E-01,  -2.03256604E-01,  -2.21254002E-01, ]
+chi0_coefficients_degree5: [  1.31014812E+00,  -1.52234989E-01,  -3.35038810E-02,  -8.28815298E-03,  -1.06049292E-03,   1.45476754E-03, ]
+chi_coefficients_degree5: [  1.31014819E+00,  -4.04019800E-01,  -2.16856067E-01,  -1.95748183E-01,  -2.21254002E-01,  -2.57431330E-01, ]
+...
+

abipy/electrons/lruj.py

@@ -109,7 +109,7 @@ def from_file(cls, filepath: PathLike):
           pert_name = r"$\alpha$"
           metric = r"N $(n^{\uparrow} + n^{\downarrow})$"
           parname = 'U'
-        
+
         chi0_coefficients = {}
         chi_coefficients = {}
         for k, v in data.items():
@@ -145,7 +145,7 @@ def find(header, dtype=None):
         -0.1500000676   8.6964981922   8.6520722003
 
        -OR-
- 
+
         Perturbations         Magnetizations
         --------------- -----------------------------
            beta [eV]     Unscreened      Screened
@@ -189,7 +189,9 @@ def find(header, dtype=None):
 #===============================================================================================================
 #===============================================================================================================
     @add_fig_kwargs
-    def plot(self, ax=None, degrees="all", inset=True, insetdegree=1, insetlocale="lower left", ptcolor0='k', ptcolor='k', gradcolor1='#3575D5',gradcolor2='#FDAE7B', ptitle="default", fontsize=12, **kwargs) -> Figure:
+    def plot(self, ax=None, degrees="all", inset=True, insetdegree=1, insetlocale="lower left",
+             ptcolor0='k', ptcolor='k', gradcolor1='#3575D5',gradcolor2='#FDAE7B',
+             ptitle="default", fontsize=12, **kwargs) -> Figure:
         """
         Plot
 
@@ -220,7 +222,7 @@ def plot(self, ax=None, degrees="all", inset=True, insetdegree=1, insetlocale="l
         xstart, xstop = 1.1 * self.alphas.min(), 1.1 * self.alphas.max()
         xs = np.arange(xstart, xstop, step=0.01)
 
-        #Prepare colors and coefficients for polynomials the use wants to plot
+        # Prepare colors and coefficients for polynomials the use wants to plot
         if degrees == "all":
           degrees = self.chi0_coefficients.keys()
 
@@ -288,6 +290,7 @@ def dfvalue(keywo):
           parambox.patch.set_edgecolor(insetcolor)
           parambox.patch.set_facecolor('white')
           ax.add_artist(parambox)
+
         return fig
 
 

abipy/examples/plot/plot_lruj.py

@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+r"""
+Plot the LRUJ results
+=====================
+
+This example shows how to parse the output file produced by lruj and plot the results
+See also https://docs.abinit.org/tutorial/lruj/#43-execution-of-the-lruj-post-processinng-utility
+"""
+from abipy.electrons.lruj import LrujAnalyzer, LrujResults
+import abipy.data as abidata
+import os
+
+# Initialize LrujResults from the main output file of lruj
+outfile = abidata.ref_file("lruj_data/lruj.out")
+lr = LrujResults.from_file(outfile)
+
+#%%
+# Plot the fits.
+
+lr.plot(degrees="all", insetdegree=4, ptcolor0='blue',
+         ptitle="Hello World", fontsize=9)
+
+#filepaths = [
+#    "tlruj_2.o_DS1_LRUJ.nc",
+#    "tlruj_2.o_DS2_LRUJ.nc",
+#    "tlruj_2.o_DS3_LRUJ.nc",
+#    "tlruj_2.o_DS4_LRUJ.nc",
+#]
+#
+#root = "tutorial_tlruj_1-tlruj_2-tlruj_3/"
+#filepaths = [os.path.join(root, p) for p in filepaths]
+#
+#lruj = LrujAnalyzer(verbose=1)
+#lruj.add_ncpaths("foo", filepaths)
+#lruj.add_ncpaths("bar", filepaths)
+
+#print(lruj)
+#lruj.plot()

abipy/ml/aseml.py

@@ -198,8 +198,101 @@ def diff_two_structures(label1, structure1, label2, structure2, fmt, file=sys.st
         print(l1.ljust(pad), " | ", l2, file=file)
 
 
+class AseTrajectoryPlotter:
+    """
+    Plot an ASE trajectory with matplotlib.
+    """
+
+    def __init__(self, traj: Trajectory):
+        self.traj = traj
+
+    @classmethod
+    def from_file(cls, filepath: PathLike) -> AseTrajectoryPlotter:
+        """Initialize object from file."""
+        return cls(read(filepath, index=":"))
+
+    def __str__(self) -> str:
+        return self.to_string()
+
+    def to_string(self, verbose=0) -> str:
+        """String representation with verbosity level verbose."""
+        lines = [f"ASE trajectory with {len(self.traj)} configuration(s)."]
+        app = lines.append
+        if len(self.traj) == 1:
+            first = AseResults.from_atoms(self.traj[0])
+            app(first.to_string(verbose=verbose))
+        else:
+            first, last = AseResults.from_atoms(self.traj[0]), AseResults.from_atoms(self.traj[-1])
+            raise NotImplementedError()
+
+        return "\n".join(lines)
+
+    #@add_fig_kwargs
+    #def plot_lattice(self, what_list=("abc", "angles", "volume"), ax_list=None,
+    #                 fontsize=8, xlims=None, **kwargs) -> Figure:
+    #    """
+    #    Plot lattice lengths/angles/volume as a function the trajectory index.
+    #    """
+    #    energies = [ene=float(atoms.get_potential_energy()) for atoms in self.traj]
+    #
+    #    stress_voigt = atoms.get_stress()
+    #    forces=atoms.get_forces(),
+    #    try:
+    #        magmoms = atoms.get_magnetic_moments()
+    #    except PropertyNotImplementedError:
+    #        magmoms = None
+
+    #@add_fig_kwargs
+    #def plot_lattice(self, what_list=("abc", "angles", "volume"), ax_list=None,
+    #                 fontsize=8, xlims=None, **kwargs) -> Figure:
+    #    """
+    #    Plot lattice lengths/angles/volume as a function the trajectory index.
+
+    #    Args:
+    #        what_list: List of strings specifying the quantities to plot. Default all
+    #        ax_list: List of axis or None if a new figure should be created.
+    #        fontsize: fontsize for legends and titles
+    #        xlims: Set the data limits for the x-axis. Accept tuple e.g. ``(left, right)``
+    #               or scalar e.g. ``left``. If left (right) is None, default values are used.
+    #    """
+    #    what_list = list_strings(what_list)
+    #    ax_list, fig, plt = get_axarray_fig_plt(ax_list, nrows=1, ncols=len(what_list),
+    #                                            sharex=True, sharey=False, squeeze=False)
+    #    markers = ["o", "^", "v"]
+
+    #    if "abc" in what_list:
+    #        # plot lattice parameters.
+    #        for i, label in enumerate(["a", "b", "c"]):
+    #            ax.plot(self.times, [lattice.abc[i] for lattice in self.lattices],
+    #                    label=label, marker=markers[i])
+    #        ax.set_ylabel("abc (A)")
+
+    #    if "angles" in what_list:
+    #        # plot lattice angles.
+    #        for i, label in enumerate(["alpha", "beta", "gamma"]):
+    #            ax.plot(self.times, [lattice.angles[i] for lattice in self.lattices],
+    #                    label=label, marker=markers[i])
+    #        ax.set_ylabel(r"$\alpha\beta\gamma$ (degree)")
+
+    #    if "volume" in what_list:
+    #        # plot lattice volume.
+    #        marker = "o"
+    #        ax.plot(self.times, [lattice.volume for lattice in self.lattices],
+    #                label="Volume", marker=marker)
+    #        ax.set_ylabel(r'$V\, (A^3)$')
+
+    #    for ix, ax in enumerate(ax_list):
+    #        set_axlims(ax, xlims, "x")
+    #        if ix == len(ax_list) - 1:
+    #            ax.set_xlabel('t (ps)', fontsize=fontsize)
+    #        ax.legend(loc="best", shadow=True, fontsize=fontsize)
+
+    #    return fig
+
+
+
 @dataclasses.dataclass
-class AseResults:
+class AseResults(HasPickleIO):
     """
     Container with the results produced by an ASE calculator.
     """
@@ -262,12 +355,14 @@ def to_string(self, verbose: int = 0) -> str:
         lines = []; app = lines.append
 
         app(f"Energy: {self.ene} (eV)")
-        app(f"Pressure: {self.pressure} ")
+        app(f"Pressure: {self.pressure} (Gpa)")
+
         fstats = self.get_fstats()
         for k, v in fstats.items():
-            app(f"{k} = {v}")
-        #app('Stress tensor:', r.stress)
-        if verbose:
+            app(f"{k} = {v} (eV/Ang)")
+
+        #if verbose:
+        if True:
             app('Forces (eV/Ang):')
             positions = self.atoms.get_positions()
             df = pd.DataFrame(dict(
@@ -278,7 +373,15 @@ def to_string(self, verbose: int = 0) -> str:
                 fy=self.forces[:,1],
                 fz=self.forces[:,2],
             ))
-            app(str(df))
+            app(df.to_string())
+
+        #if self.magmoms is not None:
+        #    for ia, (atom, magmoms) in enumerate(zip(self.atoms, self.magmoms)):
+        #        print(atom, magmoms)
+
+        app('Stress tensor:')
+        for row in self.strees:
+            app(str(row))
 
         return "\n".join(lines)
 
@@ -325,19 +428,6 @@ def diff_stats(xs, ys):
     )
 
 
-def make_square_axes(ax_mat):
-    """
-    Make an axes square in screen units.
-    Should be called after plotting.
-    """
-    return
-    for ax in ax_mat.flat:
-        #ax.set_aspect(1 / ax.get_data_ratio())
-        #ax.set(adjustable='box', aspect='equal')
-        ax.set(adjustable='datalim', aspect='equal')
-    #ax.set_aspect(1 / ax.get_data_ratio())
-
-
 class AseResultsComparator(HasPickleIO):
     """
     This object allows one to compare energies, forces and stresses computed
@@ -595,7 +685,6 @@ def plot_energies(self, fontsize=8, **kwargs):
                ax.set_title(f"{key1}/{key2} MAE: {stats.MAE:.6f}", fontsize=fontsize)
 
         if "title" not in kwargs: fig.suptitle(f"Energies in eV for {self.structure.latex_formula}")
-        #make_square_axes(ax_mat)
         return fig
 
     @add_fig_kwargs
@@ -628,7 +717,6 @@ def plot_forces(self, fontsize=8, **kwargs):
                 ax.set_title(f"{key1}/{key2} MAE: {stats.MAE:.6f}", fontsize=fontsize)
 
         if "title" not in kwargs: fig.suptitle(f"Cartesian forces in ev/Ang for {self.structure.latex_formula}")
-        #make_square_axes(ax_mat)
         return fig
 
     @add_fig_kwargs
@@ -661,7 +749,6 @@ def plot_stresses(self, fontsize=6, **kwargs):
                 ax.set_title(f"{key1}/{key2} MAE: {stats.MAE:.6f}", fontsize=fontsize)
 
         if "title" not in kwargs: fig.suptitle(f"Stresses in (eV/Ang^2) for {self.structure.latex_formula}")
-        #make_square_axes(ax_mat)
         return fig
 
     @add_fig_kwargs
@@ -2787,6 +2874,29 @@ def run(self, steps: int):
         self.dyn.run(steps)
 
 
+class GsMl(MlBase):
+    """
+    Single point calculation of energy, forces and stress with ML potential.
+    """
+    def __init__(self, atoms, nn_name, verbose, workdir, prefix=None):
+        super().__init__(workdir, prefix)
+        self.atoms = atoms
+        self.nn_name = nn_name
+        self.verbose = verbose
+
+    def run(self):
+        calc = CalcBuilder(self.nn_name).get_calculator()
+        self.atoms.calc = calc
+        res = AseResults.from_atoms(self.atoms)
+        print(res.to_string(verbose=self.verbose))
+
+        # Write ASE traj file with results.
+        with open(self.workdir / "gs.traj", "wb") as fd:
+            write_traj(fd, [self.atoms])
+
+        return 0
+
+
 class MlCompareNNs(MlBase):
     """
     Compare energies, forces and stresses obtained with different ML potentials.