add fix for compscale

Author: srmnitc

calphy/composition_transformation.py

@@ -172,7 +172,17 @@ def convert_to_pyscal(self):
         """
         Convert a given system to pyscal and give a dict of type mappings
         """
-        aseobj = read(self.calc.lattice, format="lammps-data", style="atomic")
+        # Create Z_of_type mapping to properly read LAMMPS data files
+        # This ensures atoms are correctly identified by their element
+        Z_of_type = dict(
+            [
+                (count + 1, element(el).atomic_number)
+                for count, el in enumerate(self.calc.element)
+            ]
+        )
+        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
@@ -191,7 +201,6 @@ def convert_to_pyscal(self):
         self.actual_species = len(self.typedict)
         self.new_species = len(self.output_chemical_composition) - len(self.typedict)
         self.maxtype = self.actual_species + 1  # + self.new_species
-        # print(self.typedict)
 
     def get_composition_transformation(self):
         """
@@ -215,26 +224,28 @@ def get_composition_transformation(self):
         self.to_remove = to_remove
         self.to_add = to_add
 
-    def get_random_index_of_species(self, species):
+    def get_random_index_of_species(self, species_name):
         """
-        Get a random index of a given species
+        Get a random index of a given species by element name
         """
-        ids = [count for count, x in enumerate(self.atom_type) if x == species]
+        ids = [count for count, x in enumerate(self.atom_species) if x == species_name]
         return ids[np.random.randint(0, len(ids))]
 
     def mark_atoms(self):
         for i in range(self.natoms):
             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.mappings = [f"{x}-{x}" for x in self.atom_type]
+        self.mappings = [f"{x}-{x}" for x in self.atom_species]
 
     def update_mark_atoms(self):
         self.marked_atoms = []
         for key, val in self.to_remove.items():
-            # print(f"Element {key}, count {val}")
+            # key is the element name (e.g., "Mg")
             for i in range(100000):
-                rint = self.get_random_index_of_species(self.typedict[key])
+                rint = self.get_random_index_of_species(key)
                 if rint not in self.marked_atoms:
                     self.atom_mark[rint] = True
                     self.marked_atoms.append(rint)
@@ -257,22 +268,20 @@ def update_typedicts(self):
 
     def compute_possible_mappings(self):
         self.possible_mappings = []
-        # Now make a list of possible mappings
+        # Now make a list of possible mappings using element names
         for key1, val1 in self.to_remove.items():
             for key2, val2 in self.to_add.items():
                 mapping = f"{key1}-{key2}"
                 if mapping not in self.restrictions:
-                    self.possible_mappings.append(
-                        f"{self.typedict[key1]}-{self.typedict[key2]}"
-                    )
+                    self.possible_mappings.append(mapping)
 
     def update_mappings(self):
         marked_atoms = self.marked_atoms.copy()
         for key, val in self.to_add.items():
             # now get all
 
             # we to see if we can get val number of atoms from marked ones
-            if val < len(marked_atoms):
+            if val > len(marked_atoms):
                 raise ValueError(
                     f"Not enough atoms to choose {val} from {len(marked_atoms)} not possible"
                 )
@@ -284,8 +293,8 @@ def update_mappings(self):
                 to_del = []
                 for x in range(len(self.marked_atoms)):
                     random_choice = np.random.choice(marked_atoms)
-                    # find corresponding mappiong
-                    mapping = f"{self.atom_type[random_choice]}-{self.typedict[key]}"
+                    # find corresponding mapping using species name
+                    mapping = f"{self.atom_species[random_choice]}-{key}"
                     if mapping in self.possible_mappings:
                         # this is a valid choice
                         self.mappings[random_choice] = mapping
@@ -314,12 +323,8 @@ def update_mappings(self):
             mapsplit = mapping.split("-")
             if not mapsplit[0] == mapsplit[1]:
                 transformation_dict = {}
-                transformation_dict["primary_element"] = self.reversetypedict[
-                    int(mapsplit[0])
-                ]
-                transformation_dict["secondary_element"] = self.reversetypedict[
-                    int(mapsplit[1])
-                ]
+                transformation_dict["primary_element"] = mapsplit[0]
+                transformation_dict["secondary_element"] = mapsplit[1]
                 transformation_dict["count"] = self.unique_mapping_counts[count]
                 self.transformation_list.append(transformation_dict)
 
@@ -333,25 +338,23 @@ def prepare_pair_lists(self):
         self.pair_list_new = []
         for mapping in self.unique_mappings:
             map_split = mapping.split("-")
-            # conserved atom
+            # conserved atom - mappings now use element names directly
             if map_split[0] == map_split[1]:
-                self.pair_list_old.append(self.reversetypedict[int(map_split[0])])
-                self.pair_list_new.append(self.reversetypedict[int(map_split[0])])
+                self.pair_list_old.append(map_split[0])
+                self.pair_list_new.append(map_split[0])
             else:
-                self.pair_list_old.append(self.reversetypedict[int(map_split[0])])
-                self.pair_list_new.append(self.reversetypedict[int(map_split[1])])
+                self.pair_list_old.append(map_split[0])
+                self.pair_list_new.append(map_split[1])
         self.new_atomtype = np.array(range(len(self.unique_mappings))) + 1
         self.mappingdict = dict(zip(self.unique_mappings, self.new_atomtype))
 
     def update_types(self):
+        # Update atom_type based on mapping to new types
         for x in range(len(self.atom_type)):
             self.atom_type[x] = self.mappingdict[self.mappings[x]]
 
-        # smartify these loops
-        # npyscal = len(self.pyscal_structure.atoms.types)
+        # Update pyscal structure types
         self.pyscal_structure.atoms.types = self.atom_type
-        # for count in range(npyscal)):
-        #    self.pyscal_structure.atoms.types[count] = self.atom_type[count]
 
     def iselement(self, symbol):
         try:
@@ -392,7 +395,7 @@ def update_pair_coeff(self, pair_coeff):
                 # If element_group matches our pair_list_old, replace with pair_list_new
                 # Otherwise replace with pair_list_old (for the old/reference command)
                 if element_group == self.pair_list_old or set(element_group) == set(
-                    self.element
+                    self.calc.element
                 ):
                     # This needs special handling - we'll mark position for later
                     result_parts.append("__ELEMENTS__")

tests/test_composition_transformation_bug.py

@@ -0,0 +1,182 @@
+"""
+Test for composition transformation bug fix.
+
+This test verifies that composition_transformation.py correctly handles
+LAMMPS data files where all atoms are of a single type that needs to be
+partially transformed to another element.
+
+Bug: When reading LAMMPS data files, the code wasn't providing Z_of_type
+mapping, causing ASE to incorrectly identify atom types. Additionally, the
+code was using pyscal's numeric types instead of element species names.
+"""
+
+import pytest
+import numpy as np
+import os
+import tempfile
+from calphy.composition_transformation import CompositionTransformation
+from calphy.phase_diagram import SimpleCalculation
+
+
+@pytest.fixture
+def mg_structure_file():
+    """Create a temporary LAMMPS data file with only Mg atoms (all type 2)."""
+    with tempfile.TemporaryDirectory() as tmpdir:
+        filepath = os.path.join(tmpdir, "Mg.lammps.data")
+        
+        # Simple HCP Mg structure - all atoms are type 2
+        natoms = 2048
+        content = f"""(written by ASE)
+
+{natoms} atoms
+2 atom types
+
+0.0                   25.68  xlo xhi
+0.0      44.479064738368763  ylo yhi
+0.0                41.93544  zlo zhi
+
+Atoms # atomic
+
+"""
+        
+        # Generate atom positions for a simple HCP lattice
+        nx, ny, nz = 8, 8, 16  # Should give 8*8*16*2 = 2048 atoms
+        lx, ly, lz = 25.68, 44.479064738368763, 41.93544
+        
+        atom_id = 1
+        for ix in range(nx):
+            for iy in range(ny):
+                for iz in range(nz):
+                    for sublat in range(2):
+                        x = (ix + sublat * 0.5) * lx / nx
+                        y = (iy + sublat * 0.5) * ly / ny
+                        z = (iz + sublat * 0.5) * lz / nz
+                        content += f"  {atom_id:4d}   2  {x:20.15f}  {y:20.15f}  {z:20.15f}\n"
+                        atom_id += 1
+        
+        with open(filepath, 'w') as f:
+            f.write(content)
+        
+        yield filepath
+
+
+def test_composition_transformation_single_type_file(mg_structure_file):
+    """
+    Test composition transformation when starting structure has all atoms
+    of a single type (type 2 = Mg) and we want to add Al atoms.
+    
+    Input: 2048 Mg atoms (all type 2 in LAMMPS file)
+    Output: 102 Al + 1946 Mg atoms
+    
+    This test reproduces the bug where the transformation fails with
+    "ValueError: high <= 0" due to incorrect type/species handling.
+    """
+    # Create a simple calculation object
+    calc = SimpleCalculation(
+        lattice=mg_structure_file,
+        element=["Al", "Mg"],
+        input_chemical_composition={"Al": 0, "Mg": 2048},
+        output_chemical_composition={"Al": 102, "Mg": 1946}
+    )
+    
+    # Create composition transformation
+    comp = CompositionTransformation(calc)
+    
+    # Verify basic properties
+    assert comp.natoms == 2048
+    assert comp.typedict == {"Al": 1, "Mg": 2}
+    
+    # Verify input/output compositions
+    assert comp.input_chemical_composition == {"Al": 0, "Mg": 2048}
+    assert comp.output_chemical_composition == {"Al": 102, "Mg": 1946}
+    
+    # Verify transformation requirements
+    assert comp.to_remove == {"Mg": 102}
+    assert comp.to_add == {"Al": 102}
+    
+    # Check that we have exactly one transformation: Mg -> Al
+    assert len(comp.transformation_list) == 1
+    assert comp.transformation_list[0]['primary_element'] == 'Mg'
+    assert comp.transformation_list[0]['secondary_element'] == 'Al'
+    assert comp.transformation_list[0]['count'] == 102
+    
+    # Verify pair lists
+    assert comp.pair_list_old == ['Mg', 'Mg']
+    assert comp.pair_list_new == ['Al', 'Mg']
+    
+    # Test pair_coeff update
+    pair_coeff_input = "pair_coeff * * /path/to/potential.yace Al Mg"
+    pc_old, pc_new = comp.update_pair_coeff(pair_coeff_input)
+    
+    # Verify pair coeffs are different and properly formatted
+    # pc_old should have Mg because all atoms start as Mg
+    # pc_new should have both Al and Mg after transformation
+    assert "Mg" in pc_old
+    assert "Al" in pc_new and "Mg" in pc_new
+    
+    # Write structure and verify
+    output_structure = os.path.join(os.path.dirname(mg_structure_file), "transformed.lammps.data")
+    comp.write_structure(output_structure)
+    
+    assert os.path.exists(output_structure)
+    
+    # Read and verify the written structure
+    with open(output_structure, 'r') as f:
+        content = f.read()
+        assert "2048 atoms" in content
+    
+    # Verify entropy contribution calculation
+    entropy = comp.entropy_contribution
+    assert isinstance(entropy, float)
+    assert entropy < 0  # Should be negative for this transformation
+
+
+def test_composition_transformation_with_pace_potential(mg_structure_file):
+    """
+    Test with PACE potential format (as in the original bug report).
+    Verifies that pair_coeff commands are correctly generated for PACE potentials.
+    """
+    # Create a simple calculation object
+    calc = SimpleCalculation(
+        lattice=mg_structure_file,
+        element=["Al", "Mg"],
+        input_chemical_composition={"Al": 0, "Mg": 2048},
+        output_chemical_composition={"Al": 102, "Mg": 1946}
+    )
+    
+    # Create composition transformation
+    comp = CompositionTransformation(calc)
+    
+    # Test PACE potential format
+    pair_coeff_input = "pair_coeff * * /home/user/AlMgZn.yace Al Mg"
+    pc_old, pc_new = comp.update_pair_coeff(pair_coeff_input)
+    
+    # Both should contain the potential file path
+    assert "/home/user/AlMgZn.yace" in pc_old
+    assert "/home/user/AlMgZn.yace" in pc_new
+    
+    # Verify element specifications
+    assert "Mg" in pc_old
+    assert "Al" in pc_new and "Mg" in pc_new
+
+
+def test_composition_transformation_swap_types(mg_structure_file):
+    """
+    Test get_swap_types method returns correct type mappings.
+    """
+    # Create a simple calculation object
+    calc = SimpleCalculation(
+        lattice=mg_structure_file,
+        element=["Al", "Mg"],
+        input_chemical_composition={"Al": 0, "Mg": 2048},
+        output_chemical_composition={"Al": 102, "Mg": 1946}
+    )
+    
+    # Create composition transformation
+    comp = CompositionTransformation(calc)
+    
+    # Get swap types
+    swap_types = comp.get_swap_types()
+    
+    assert isinstance(swap_types, list)
+    assert len(swap_types) == 2