Ase structures

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 2.1.2
+current_version = 2.2.0
 commit = True
 tag = True
 

CITATION.cff

@@ -61,4 +61,4 @@ abstract: >-
   examples illustrate the accuracy and efficiency of
   our implementation.
 license-url: 'https://github.com/ICAMS/calphy/blob/main/LICENSE'
-version: 2.1.2
+version: 2.2.0

calphy/__init__.py

@@ -28,7 +28,7 @@
 from calphy.alchemy import Alchemy
 from calphy.routines import MeltingTemp
 
-__version__ = "2.1.2"
+__version__ = "2.2.0"
 
 def addtest(a,b):
     return a+b

calphy/composition_transformation.py

@@ -27,11 +27,9 @@
 import numpy as np
 import os
 import random
-import pyscal3.core as pc
 from mendeleev import element
 from ase.io import read, write
 from ase.atoms import Atoms
-from pyscal3.core import element_dict
 from calphy.integrators import kb
 
 
@@ -163,9 +161,9 @@ def _log(val):
         entropy_term = kb * np.sum(ents)
         return entropy_term
 
-    def convert_to_pyscal(self):
+    def read_structure(self):
         """
-        Convert a given system to pyscal and give a dict of type mappings
+        Read the input structure and set up the element <-> type mappings
         """
         # Create Z_of_type mapping to properly read LAMMPS data files
         # This ensures atoms are correctly identified by their element
@@ -178,20 +176,19 @@ def convert_to_pyscal(self):
         aseobj = read(
             self.calc.lattice, format="lammps-data", style="atomic", Z_of_type=Z_of_type
         )
-        pstruct = pc.System(aseobj, format="ase")
 
         # here we have to validate the input composition dict; and map it
-        typelist = pstruct.atoms.species
+        typelist = aseobj.get_chemical_symbols()
         types, typecounts = np.unique(typelist, return_counts=True)
         composition = {types[x]: typecounts[x] for x in range(len(types))}
 
         atomsymbols = self.calc.element
         atomtypes = [x + 1 for x in range(len(self.calc.element))]
 
-        self.pyscal_structure = pstruct
+        self.structure = aseobj
         self.typedict = dict(zip(atomsymbols, atomtypes))
         self.reversetypedict = dict(zip(atomtypes, atomsymbols))
-        self.natoms = self.pyscal_structure.natoms
+        self.natoms = len(self.structure)
 
         # Count of actual unique atom types present in the structure
         # This matches what's declared in the LAMMPS data file header
@@ -239,8 +236,9 @@ def mark_atoms(self):
             self.atom_mark.append(False)
 
         # Use species (element symbols) instead of numeric types
-        self.atom_species = self.pyscal_structure.atoms.species
-        self.atom_type = self.pyscal_structure.atoms.types
+        self.atom_species = self.structure.get_chemical_symbols()
+        # placeholder types; update_types assigns the mapping type of every atom
+        self.atom_type = [1] * self.natoms
         self.mappings = [f"{x}-{x}" for x in self.atom_species]
 
     def update_mark_atoms(self):
@@ -378,9 +376,6 @@ def update_types(self):
         for x in range(len(self.atom_type)):
             self.atom_type[x] = self.mappingdict[self.mappings[x]]
 
-        # Update pyscal structure types
-        self.pyscal_structure.atoms.types = self.atom_type
-
     def iselement(self, symbol):
         try:
             _ = element(symbol)
@@ -582,9 +577,8 @@ def write_structure(self, outfilename, for_fe_mode=False):
         from ase.io import write as ase_write
         from ase import Atoms as ASEAtoms
 
-        # Get positions and cell from pyscal structure
-        positions = self.pyscal_structure.atoms.positions
-        cell = self.pyscal_structure.box
+        positions = self.structure.positions
+        cell = self.structure.cell
 
         # Build per-atom type numbers and chemical symbols.
         # for_fe_mode: remap each atom to its TARGET element's typedict entry so
@@ -598,8 +592,8 @@ def write_structure(self, outfilename, for_fe_mode=False):
             invert_typedict = {v: k for k, v in self.typedict.items()}
             symbols = [invert_typedict[t] for t in fe_types]
         else:
-            fe_types = None  # will use pyscal types below
-            symbols = [self.reversetypedict[t] for t in self.pyscal_structure.atoms.types]
+            fe_types = None  # will use the mapping types below
+            symbols = [self.reversetypedict[t] for t in self.atom_type]
 
         ase_atoms = ASEAtoms(symbols=symbols, positions=positions, cell=cell, pbc=True)
 
@@ -611,7 +605,7 @@ def write_structure(self, outfilename, for_fe_mode=False):
             lines = f.readlines()
 
         # Choose which per-atom type list to use
-        atom_types_to_write = fe_types if for_fe_mode else list(self.pyscal_structure.atoms.types)
+        atom_types_to_write = fe_types if for_fe_mode else list(self.atom_type)
 
         # Find the Atoms section and replace type numbers
         # Support different ASE formats: "Atoms # atomic", "Atoms # full", or just "Atoms"
@@ -635,7 +629,7 @@ def write_structure(self, outfilename, for_fe_mode=False):
                         break
 
         # Verify all atoms were updated
-        expected_atoms = len(self.pyscal_structure.atoms.types)
+        expected_atoms = len(self.atom_type)
         if atom_idx != expected_atoms:
             raise RuntimeError(
                 f"Failed to update all atoms in {outfilename}. "
@@ -666,7 +660,7 @@ def prepare_mappings(self):
         self.unique_mappings = []
 
         self.get_composition_transformation()
-        self.convert_to_pyscal()
+        self.read_structure()
 
         self.mark_atoms()
         self.update_mark_atoms()

calphy/helpers.py

@@ -32,8 +32,8 @@
 
 from ase.io import read, write
 
-import pyscal3.core as pc
-from pyscal3.trajectory import Trajectory
+import pyscal3
+from pyscal3 import Trajectory
 
 from calphy.runner import (
     ExecutableRunner,
@@ -343,10 +343,8 @@ def read_data(lmp, file):
 def get_structures(file, species, index=None):
     traj = Trajectory(file)
     if index is None:
-        aseobjs = traj[:].to_ase(species=species)
-    else:
-        aseobjs = traj[index].to_ase(species=species)
-    return aseobjs
+        return traj[:].to_atoms(species=species)
+    return traj[index].to_atoms(species=species)
 
 
 def remap_box(lmp, x, y, z):
@@ -394,14 +392,13 @@ def compute_msd(lmp, options):
 
 
 def find_solid_fraction(file):
-    sys = pc.System(file)
+    atoms = Trajectory(file)[0].to_atoms()[0]
     try:
-        sys.find.neighbors(method="cutoff", cutoff=0)
+        pyscal3.find_neighbors(atoms, method="cutoff", cutoff=0)
     except RuntimeError:
-        sys.find.neighbors(method="cutoff", cutoff=5.0)
-    sys.find.solids(cluster=False)
-    solids = np.sum(sys.atoms.solid)
-    return solids
+        pyscal3.find_neighbors(atoms, method="cutoff", cutoff=5.0)
+    pyscal3.find_solids(atoms, cluster=False)
+    return np.sum(atoms.arrays["pyscal_solid"])
 
 
 def write_data(lmp, file):

calphy/input.py

@@ -48,12 +48,17 @@
 import itertools
 import os
 import warnings
-from pyscal3 import System
-from pyscal3.core import structure_dict, element_dict, _make_crystal
 from ase.io import read, write
 import shutil
 
-__version__ = "2.1.2"
+from calphy.structures import (
+    BUILTIN_LATTICES,
+    canonical_lattice,
+    default_lattice,
+    make_lattice,
+)
+
+__version__ = "2.2.0"
 
 
 def _check_equal(val):
@@ -560,6 +565,8 @@ class Calculation(_StrictInput, title="Main input class"):
 
     reference_phase: Annotated[str, Field(default="")]
     lattice_constant: Annotated[float, Field(default=0)]
+    # c/a of the built-in hcp lattice; None gives the ideal sqrt(8/3)
+    c_over_a: Annotated[Union[float, None], Field(default=None, gt=0)]
     repeat: Annotated[
         conlist(int, min_length=3, max_length=3), Field(default=[1, 1, 1])
     ]
@@ -868,71 +875,38 @@ def _validate_all(self) -> "Input":
         write_structure_file = False
         rename_structure_file = False
 
-        if self.lattice == "":
-            # fetch from dict
+        if self.lattice == "" or canonical_lattice(self.lattice) is not None:
             if len(self.element) > 1:
                 raise ValueError(
                     "Cannot create lattice for more than one element, provide a lammps-data file explicitly"
                 )
-            if self.element[0] in element_dict.keys():
-                self.lattice = element_dict[self.element[0]]["structure"]
-                self.lattice_constant = element_dict[self.element[0]][
-                    "lattice_constant"
-                ]
-            else:
-                raise ValueError(
-                    "Could not find structure, please provide lattice and lattice_constant explicitely"
-                )
-
-            if self.repeat == [1, 1, 1]:
-                self.repeat = [5, 5, 5]
-
-            structure = _make_crystal(
-                self.lattice.lower(),
-                lattice_constant=self.lattice_constant,
-                repetitions=self.repeat,
-                element=self.element,
-            )
-            structure = structure.write.ase()
-
-            # extract composition
-            types, typecounts = np.unique(
-                structure.get_chemical_symbols(), return_counts=True
-            )
-
-            for c, t in enumerate(types):
-                self._element_dict[t]["count"] = typecounts[c]
-                self._element_dict[t]["composition"] = typecounts[c] / np.sum(
-                    typecounts
-                )
+            default = default_lattice(self.element[0])
 
-            self._natoms = len(structure)
-            self._original_lattice = self.lattice.lower()
-            write_structure_file = True
+            if self.lattice == "":
+                # ground-state lattice of the element
+                if default is None:
+                    raise ValueError(
+                        "Could not find structure, please provide lattice and lattice_constant explicitely"
+                    )
+                self.lattice, self.lattice_constant = default
+                if self.repeat == [1, 1, 1]:
+                    self.repeat = [5, 5, 5]
+            elif self.lattice_constant == 0:
+                if default is None:
+                    raise ValueError("Please provide lattice_constant!")
+                self.lattice_constant = default[1]
 
-        elif self.lattice.lower() in structure_dict.keys():
-            if len(self.element) > 1:
+            if self.c_over_a is not None and canonical_lattice(self.lattice) != "hcp":
                 raise ValueError(
-                    "Cannot create lattice for more than one element, provide a lammps-data file explicitly"
+                    f"c_over_a only applies to the hcp lattice, not {self.lattice}"
                 )
-
-            # this is a valid structure
-            if self.lattice_constant == 0:
-                # we try try to get lattice_constant
-                if self.element[0] in element_dict.keys():
-                    self.lattice_constant = element_dict[self.element[0]][
-                        "lattice_constant"
-                    ]
-                else:
-                    raise ValueError("Please provide lattice_constant!")
-            # now create lattice
-            structure = _make_crystal(
-                self.lattice.lower(),
-                lattice_constant=self.lattice_constant,
-                repetitions=self.repeat,
-                element=self.element,
+            structure = make_lattice(
+                self.lattice,
+                self.element[0],
+                self.lattice_constant,
+                self.repeat,
+                c_over_a=self.c_over_a,
             )
-            structure = structure.write.ase()
 
             # extract composition
             types, typecounts = np.unique(
@@ -945,14 +919,13 @@ def _validate_all(self) -> "Input":
                     typecounts
                 )
 
-            # concdict_counts = {str(t): typecounts[c] for c, t in enumerate(types)}
-            # concdict_frac = {str(t): typecounts[c]/np.sum(typecounts) for c, t in enumerate(types)}
-            # self._composition = concdict_frac
-            # self._composition_counts = concdict_counts
             self._natoms = len(structure)
-            self._original_lattice = self.lattice.lower()
+            self._original_lattice = canonical_lattice(self.lattice)
             write_structure_file = True
 
+        elif self.c_over_a is not None:
+            raise ValueError("c_over_a only applies to the built-in hcp lattice")
+
         elif self.lattice.split("-")[0] == "mp":
             # confirm here that API key exists
             if not self.materials_project.api_key:
@@ -1012,7 +985,11 @@ def _validate_all(self) -> "Input":
         else:
             # this is a file
             if not os.path.exists(self.lattice):
-                raise ValueError(f"File {self.lattice} could not be found")
+                raise ValueError(
+                    f"lattice '{self.lattice}' is neither a built-in lattice "
+                    f"({', '.join(BUILTIN_LATTICES)}), a materials project id "
+                    "nor an existing lammps-data file"
+                )
             if self.file_format == "lammps-data":
                 # create atomic numbers for proper reading
                 Z_of_type = dict(

calphy/integrators.py

@@ -42,7 +42,6 @@
 from numpy import trapezoid as trapz
 
 from tqdm import tqdm
-import pyscal3.core as pc
 from ase.io import read
 
 # Constants

calphy/postprocessing.py

@@ -374,7 +374,19 @@ def gather_results(
     if not extract_phase_prefix:
         del datadict["phase_name"]
 
-    df = pd.DataFrame(data=datadict)
+    # pandas >= 3 stores text columns with a string dtype that turns None into
+    # NaN; keep them as object so a missing value (error_code of a successful
+    # run, say) stays None
+    df = pd.DataFrame(
+        data={
+            key: (
+                pd.Series(val, dtype=object)
+                if any(isinstance(v, str) for v in val)
+                else val
+            )
+            for key, val in datadict.items()
+        }
+    )
     return df