Skip to content

Commit

Permalink
Add support for xyz files via ASE
Browse files Browse the repository at this point in the history
  • Loading branch information
gmatteo committed Nov 8, 2023
1 parent 01c39d3 commit d7591d4
Show file tree
Hide file tree
Showing 8 changed files with 89 additions and 48 deletions.
2 changes: 2 additions & 0 deletions abipy/abilab.py
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,8 @@ def _straceback():
# Vasp files.
("POSCAR", Structure),
(".vasp", Structure),
# ASE files
(".xyz", Structure),
#("ZINVCONV.nc", ZinvConvFile),
#("TETRATEST.nc", TetraTestFile),
# QE/CP files
Expand Down
13 changes: 10 additions & 3 deletions abipy/core/structure.py
Original file line number Diff line number Diff line change
Expand Up @@ -290,7 +290,6 @@ def from_file(cls, filepath: str, primitive: bool = False, sort: bool = False) -
raise ValueError("Cannot find structure in Abinit output file `%s`" % filepath)

elif filepath.endswith(".abivars") or filepath.endswith(".ucell"):
#print("in abivars")
with open(filepath, "rt") as fh:
return cls.from_abistring(fh.read())

Expand All @@ -314,6 +313,15 @@ def from_file(cls, filepath: str, primitive: bool = False, sort: bool = False) -

if new.__class__ != cls: new.__class__ = cls

elif filepath.endswith(".xyz"):
# ASE extended xyz format.
try:
from ase.io import read
except ImportError:
raise RuntimeError("ase is required to read xyz files. Use `pip install ase`")
atoms = read(filepath)
return cls.as_structure(atoms)

else:
# Invoke pymatgen and change class. Note that AbinitSpacegroup is missing here.
new = super().from_file(filepath, primitive=primitive, sort=sort)
Expand Down Expand Up @@ -2624,7 +2632,7 @@ def diff_structures(structures, fmt="cif", mode="table", headers=(), file=sys.st
print(diff, file=file)

else:
raise ValueError("Unsupported mode: `%s`" % str(mode))
raise ValueError(f"Unsupported {mode=}")


def structure2siesta(structure: Structure, verbose=0) -> str:
Expand Down Expand Up @@ -2686,7 +2694,6 @@ class StructDiff:
"""
Print difference among structures.
"""

def __init__(self, labels: list[str], structures):
"""
Args:
Expand Down
63 changes: 44 additions & 19 deletions abipy/dynamics/analyzer.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
from monty.collections import AttrDict #, dict2namedtuple
from pymatgen.util.string import latexify
from pymatgen.core.lattice import Lattice
from pymatgen.core.composition import Composition
from abipy.core.mixins import TextFile, NotebookWriter
from abipy.core.structure import Structure
from abipy.tools.typing import Figure, PathLike
Expand All @@ -46,6 +47,21 @@
kBoltzEv = 8.617333e-05
#nCar = 56 # FIXME: Harcoded


def common_oxidation_states():
oxi2symbols = {
+1: ["Li", "Na", "K", "Rb", "Cs", "Fr"],
+2: ["Be", "Mg", "Ca", "Sr", "Ba", "Ra"],
}

# map: symbol to oxidation state.
symb2oxi = {}
for oxi, slist in oxi2symbols.items():
for s in slist:
symb2oxi[s] = oxi
return symb2oxi


def read_structure_postac_ucmats(traj_filepath: PathLike, step_skip: int) -> tuple[Structure, np.ndarray, np.ndarray, int]:
"""
Read all configurations from an ASE trajectory file.
Expand Down Expand Up @@ -711,7 +727,7 @@ def get_dw_symbol(self, symbol, t0: float = 0.0, tmax=None, atom_inds=None):
label = r"Linear fit: D={:.2E} cm$^2$/s, $r^2$={:.2f}".format(naive_d, fit.rvalue**2)
return AttrDict(naive_d=naive_d, ts=ts, fit=fit, label=label)

def get_msdtt0_symbol(self, symbol: str, tmax: float, atom_inds=None, nprocs=None) -> Msdtt0:
def get_msdtt0_symbol_tmax(self, symbol: str, tmax: float, atom_inds=None, nprocs=None) -> Msdtt0:
r"""
Calculates the MSD for every possible pair of time points using the formula:
Expand Down Expand Up @@ -1571,10 +1587,10 @@ def iter_mdatc(self):
for itemp, mda in enumerate(self):
yield mda, mda.temperature, self.temp_cmap(float(itemp) / len(self))

def get_msdtt0_symbol(self, symbol: str, tmax: float, atom_inds=None, nprocs=None) -> Msdtt0List:
def get_msdtt0_symbol_tmax(self, symbol: str, tmax: float, atom_inds=None, nprocs=None) -> Msdtt0List:
msdtt0_list = Msdtt0List()
for mda in self:
obj = mda.get_msdtt0_symbol(symbol, tmax, atom_inds=atom_inds, nprocs=nprocs)
obj = mda.get_msdtt0_symbol_tmax(symbol, tmax, atom_inds=atom_inds, nprocs=nprocs)
msdtt0_list.append(obj)

return msdtt0_list
Expand Down Expand Up @@ -1602,7 +1618,7 @@ def _nrows_ncols_nplots(self, size=None):
# data_list = DiffusionDataList()
# for it, ((mda, temp), params) in enumerate(zip(self.iter_mdat(), params_list)):
# p = AttrDict(**params)
# msq_tt0 = mda.get_msdtt0_symbol(symbol, p.tmax, nprocs=nprocs)
# msq_tt0 = mda.get_msdtt0_symbol_tmax(symbol, p.tmax, nprocs=nprocs)
# sigma = msq_tt0.get_sigma_berend(t1=p.t1, t2=p.t2)
# data_t = msq_tt0.get_diffusion_with_sigma(p.fit_time_start, p.fit_time_stop, p.block_size1, p.block_size2, sigma)
# data_list.append(data_t)
Expand Down Expand Up @@ -1781,15 +1797,16 @@ class Entry:

key: str
symbol: str
formula: str
#formula: str
composition: str
temps: np.ndarray
diffusions: np.ndarray
err_diffusions: np.ndarray
volumes: np.ndarray
style: dict

@classmethod
def from_file(cls, filepath: PathLike, key, symbol, formula, style) -> Entry:
def from_file(cls, filepath: PathLike, key, symbol, composition, style) -> Entry:

if str(filepath).endswith(".dat"):
try:
Expand All @@ -1803,7 +1820,7 @@ def from_file(cls, filepath: PathLike, key, symbol, formula, style) -> Entry:
else:
try:
# Read data in CSV format. Assuming header with at least the following entries:
# T,diffusion,err_diffusion,volume,symbol,formula
# T,diffusion,err_diffusion,volume,symbol,composition
df = pd.read_csv(filepath, skipinitialspace=True) #, delim_whitespace=True)

def get_unique(col):
Expand All @@ -1813,7 +1830,7 @@ def get_unique(col):
return v0

symbol = get_unique("symbol")
formula = get_unique("formula")
composition = get_unique("composition")
temps = df["T"].values
volumes = df["volume"].values
diffusions = df["diffusion"].values
Expand All @@ -1824,9 +1841,12 @@ def get_unique(col):
except Exception as exc:
raise RuntimeError(f"Exception while reading {filepath=}") from exc

composition = Composition(composition)
#print(composition)

return cls(key=key,
symbol=symbol,
formula=formula,
composition=composition,
temps=temps,
diffusions=diffusions,
err_diffusions=err_diffusions,
Expand Down Expand Up @@ -1934,7 +1954,7 @@ class ArrheniusPlotter:
from abipy.dynamics.analyzer import ArrheniusPlotter
symbol = "Li"
formula = "c-LLZO"
composition = "c-LLZO"
key_path = {
"matgl-MD": "diffusion_cLLZO-matgl.dat",
Expand All @@ -1948,19 +1968,19 @@ class ArrheniusPlotter:
plotter = ArrheniusPlotter()
for key, path in key_path.items():
plotter.add_entry_from_file(path, key, symbol, formula, style=style_key[key])
plotter.add_entry_from_file(path, key, symbol, composition, style=style_key[key])
# temperature grid refined
thinvt_arange = (0.6, 2.5, 0.01)
xlims, ylims = (0.5, 2.0), (-7, -4)
plotter.plot(thinvt_arange=thinvt_arange,
xlims=xlims, ylims=ylims, text='LLZO cubic',
savefig=None)
xlims=xlims, ylims=ylims, text='LLZO cubic', savefig=None)
"""

def __init__(self):
self.entries = []
self.sym2oxi = common_oxidation_states()

def __iter__(self):
return self.entries.__iter__()
Expand Down Expand Up @@ -2012,11 +2032,11 @@ def get_min_max_temp(self) -> tuple[float, float]:
max_temp = max([e.temperatures.max() for e in self])
return min_temp, max_temp

def add_entry_from_file(self, filepath: PathLike, key: str, symbol: str, formula: str, style=None) -> Entry:
def add_entry_from_file(self, filepath: PathLike, key: str, symbol: str, composition: str, style=None) -> Entry:
"""
"""
style = style or {}
self.append(Entry.from_file(filepath, key, symbol, formula, style))
self.append(Entry.from_file(filepath, key, symbol, composition, style))

@add_fig_kwargs
def plot(self, thinvt_arange=None, what="diffusion", ncar=None, colormap="jet", with_t=True, text=None,
Expand Down Expand Up @@ -2058,17 +2078,22 @@ def plot(self, thinvt_arange=None, what="diffusion", ncar=None, colormap="jet",
if "color" not in my_style and "c" not in my_style:
my_style["color"] = cmap(ie / len(self))

symbol = entry.symbol
composition = entry.composition

diff = entry.get_diffusion_data(fit_thinvt=fit_thinvt)
label = r"D$_{\mathrm{%s}}$ (%s): " % (entry.symbol, entry.key)
label = r"D$_{\mathrm{%s}}$ (%s): " % (symbol, entry.key)
label += r"E$_\mathrm{a}$=" + str('{:.2F}'.format(diff.e_act)) + 'eV'

if what == "diffusion":
data = diff
else:
if ncar is None:
raise ValueError("Computation of conductivity requires the specification of ncar!")
if len(self.symbols()) != 1:
raise ValueError(f"Computation of conductivity requires entries with same symbol while {self.symbols()=}!")
if symbol not in self.symb2oxi:
raise ValueError(f"No entry for {symbol=} found in symb2oxi! Please add the oxistate manually!")
ncar = symb2oxi[symbol] * composition[symbol]
print(self.symb2oxi[symbol], composition[symbol])

sigma_data, tsigma_data = entry.get_conductivity_data(ncar, fit_thinvt=fit_thinvt)
data = dict(sigma=sigma_data, tsigma=tsigma_data)[what]

Expand Down
46 changes: 23 additions & 23 deletions abipy/ml/aseml.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,22 +58,22 @@
# Helper functions
###################

def nprocs_for_ntasks(nprocs, ntasks, title=None) -> int:
"""
Return the number of procs to be used in a multiprocessing Pool.
If negative or None, use half the procs available in the system.
"""
if nprocs is None or nprocs <= 0:
nprocs = max(1, os.cpu_count() // 2)
else:
nprocs = int(nprocs)

nprocs = min(nprocs, ntasks)
if title is not None:
print(title)
print(f"Using multiprocessing pool with {nprocs=} for {ntasks=} ...")

return nprocs
#def nprocs_for_ntasks(nprocs, ntasks, title=None) -> int:
# """
# Return the number of procs to be used in a multiprocessing Pool.
# If negative or None, use half the procs available in the system.
# """
# if nprocs is None or nprocs <= 0:
# nprocs = max(1, os.cpu_count() // 2)
# else:
# nprocs = int(nprocs)
#
# nprocs = min(nprocs, ntasks)
# if title is not None:
# print(title)
# print(f"Using multiprocessing pool with {nprocs=} for {ntasks=} ...")
#
# return nprocs


_CELLPAR_KEYS = ["a", "b", "c", "angle(b,c)", "angle(a,c)", "angle(a,b)"]
Expand Down Expand Up @@ -124,7 +124,7 @@ def fix_atoms(atoms: Atoms,
fix_inds: list[int] | None = None,
fix_symbols: list[str] | None = None) -> None:
"""
Fix atoms by indices and by symbols.
Fix atoms by indices and/or by symbols.
Args:
atoms: ASE atoms
Expand Down Expand Up @@ -347,7 +347,7 @@ def make_square_axes(ax_mat):

class AseResultsComparator(HasPickleIO):
"""
This object allows one to compare energies, forces and stressee computed
This object allows one to compare energies, forces and stresses computed
for the same structure but with different methods e.g. results obtained
with different ML potentials.
"""
Expand Down Expand Up @@ -2028,6 +2028,7 @@ def postprocess_images(self, images):
ef, de = nebtools.get_barrier(fit=False)
neb_data.update(barrier_without_fit=float(ef),
energy_change_without_fit=float(de),
nn_name=str(self.nn_name),
)

self.write_json("neb_data.json", neb_data, info="JSON document with NEB results",
Expand Down Expand Up @@ -2232,15 +2233,15 @@ def run(self) -> None:
self.write_script("ase_gui.sh", text=f"""\
# To visualize the results, use:
ase gui {nebtraj_file}@-{self.nimages}
ase gui "{nebtraj_file}@-{self.nimages}"
# then select `tools->neb` in the gui.
""", info="Shell script to visualize NEB results with ase gui")

self.write_script("ase_nebplot.sh", text=f"""\
# This command create a series of plots showing the progression of the neb relaxation
ase nebplot --share-x --share-y --nimages {self.nimages} {nebtraj_file}
ase nebplot --share-x --share-y --nimages {self.nimages} "{nebtraj_file}"
""", info="Shell script to create a series of plots showing the progression of the neb relaxation")

self._finalize()
Expand Down Expand Up @@ -2621,7 +2622,6 @@ def run(self, nprocs, print_dataframes=True) -> AseResultsComparator:
results_list = [abi_results]

ntasks = len(abi_results)
#nprocs = nprocs_for_ntasks(nprocs, ntasks, title="Begin ML computation")
nprocs = 1

for nn_name in self.nn_names:
Expand Down Expand Up @@ -2775,7 +2775,7 @@ def run(self, steps: int):
self.dyn.run(steps)


class MlCompareNNs(_MlNebBase):
class MlCompareNNs(MlBase):
"""
Compare energies, forces and stresses obtaiend with different ML potentials.
Also profile the time required.
Expand Down Expand Up @@ -2821,7 +2821,7 @@ def to_string(self, verbose=0) -> str:

def run(self, print_dataframes=True) -> AseResultsComparator:
"""
Run calculations
Run calculations.
"""
workdir = self.workdir

Expand Down
5 changes: 5 additions & 0 deletions abipy/ppcodes/oncv_plotter.py
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,11 @@ def _add_rc_vlines_ax(self, ax, with_lloc=False) -> None:
color = "magenta" if self.parser.lloc == 4 else "k"
ax.axvline(self.parser.rc5, lw=2, color=color, ls="--")

def plotly_atan_logders(self, *args, **kwargs):
mpl_fig = self.plot_atan_logders(*args, show=False, **kwargs)
from plotly.tools import mpl_to_plotly
return mpl_to_plotly(mpl_fig)

@add_fig_kwargs
def plot_atan_logders(self, ax=None, with_xlabel=True,
fontsize: int = 8, **kwargs) -> Figure:
Expand Down
4 changes: 2 additions & 2 deletions abipy/scripts/abicomp.py
Original file line number Diff line number Diff line change
Expand Up @@ -229,7 +229,7 @@ def abicomp_xrd(options):
Compare X-ray diffraction plots (requires FILES with structure).
"""
if len(options.paths) < 2:
print("You need more than one structure to compare!")
print("You need more than one structures to compare!")
return 1

structures = [abilab.Structure.from_file(p) for p in options.paths]
Expand Down Expand Up @@ -584,7 +584,7 @@ def abicomp_abiwan(options):


def abicomp_pseudos(options):
""""Compare multiple pseudos Print table to terminal."""
""""Compare multiple pseudos and print table to terminal."""
# Make sure entries in index are unique.
index = [os.path.basename(p) for p in options.paths]
if len(index) != len(set(index)): index = [os.path.relpath(p) for p in options.paths]
Expand Down
2 changes: 2 additions & 0 deletions abipy/scripts/oncv.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,8 @@ def oncv_plot(options) -> int:

out_path = _find_oncv_output(options.filepath)
plotter = OncvPlotter.from_file(out_path)
plotter.plotly_atan_logders().show()
return 0

plotter.expose(slide_mode=options.slide_mode, slide_timeout=options.slide_timeout,
use_web=options.expose_web, verbose=options.verbose)
Expand Down
2 changes: 1 addition & 1 deletion abipy/tools/plotting.py
Original file line number Diff line number Diff line change
Expand Up @@ -1818,7 +1818,7 @@ def wrapper(*args, **kwargs):
================ ====================================================================
title Title of the plot (Default: None).
show True to show the figure (default: True).
hovormode True to show the hover info (default: False)
hovermode True to show the hover info (default: False)
savefig "abc.png" , "abc.jpeg" or "abc.webp" to save the figure to a file.
write_json Write plotly figure to `write_json` JSON file.
Inside jupyter-lab, one can right-click the `write_json` file from
Expand Down

0 comments on commit d7591d4

Please sign in to comment.