Fix S24 weas visualisation: triplets -> only adsorbed molecule on surface

Author: joehart2001

mlip_testing/analysis/surfaces/S24/analyse_S24.py

@@ -49,14 +49,37 @@ def system_names() -> list:
     """
     system_names = []
     for model_name in MODELS:
-        for system_path in (CALC_PATH / model_name).glob("*.xyz"):
-            if system_path.name == "s24_mol_surface_atoms.xyz":
-                continue
-            system_names.append(system_path.stem)
-        break
+        model_dir = CALC_PATH / model_name
+        if model_dir.exists():
+            for system_path in sorted(model_dir.glob("*.xyz")):
+                mol_surface = read(system_path)
+                if "system_name" in mol_surface.info:
+                    system_names.append(mol_surface.info["system_name"])
+            break
     return system_names
 
 
+def sys_ids() -> list:
+    """
+    Get list of system IDs.
+
+    Returns
+    -------
+    list
+        List of all system IDs.
+    """
+    sys_ids = []
+    for model_name in MODELS:
+        model_dir = CALC_PATH / model_name
+        if model_dir.exists():
+            for system_path in sorted(model_dir.glob("*.xyz")):
+                mol_surface = read(system_path)
+                if "sys_id" in mol_surface.info:
+                    sys_ids.append(mol_surface.info["sys_id"])
+            break
+    return sys_ids
+
+
 @pytest.fixture
 @plot_parity(
     filename=OUT_PATH / "figure_adsorption_energies.json",
@@ -65,6 +88,7 @@ def system_names() -> list:
     y_label="Reference adsorption energy / eV",
     hoverdata={
         "System": system_names(),
+        "Sys ID": sys_ids(),
     },
 )
 def adsorption_energies() -> dict[str, list]:
@@ -80,39 +104,26 @@ def adsorption_energies() -> dict[str, list]:
     ref_stored = False
 
     for model_name in MODELS:
-        for system_path in sorted((CALC_PATH / model_name).glob("*.xyz")):
-            if system_path.name == "s24_mol_surface_atoms.xyz":
-                continue
-
-            structs = read(system_path, index=":")
-            if len(structs) != 3:
-                continue
-
-            surface, mol_surface, molecule = structs
-
-            # Get predicted energies
-            surface_e = surface.get_potential_energy()
-            mol_surf_e = mol_surface.get_potential_energy()
-            molecule_e = molecule.get_potential_energy()
-            pred_ads_energy = compute_adsorption_energy(
-                surface_e, mol_surf_e, molecule_e
-            )
+        model_dir = CALC_PATH / model_name
+        if not model_dir.exists():
+            results[model_name] = []
+            continue
+            
+        for system_path in sorted(model_dir.glob("*.xyz")):
+            mol_surface = read(system_path)
+            
+            # Get pre-calculated adsorption energies
+            pred_ads_energy = mol_surface.info["adsorption_energy"]
             results[model_name].append(pred_ads_energy)
-
-            # Get reference energies (only store once)
+            
             if not ref_stored:
-                ref_surface_e = surface.info["ref_energy"]
-                ref_mol_surf_e = mol_surface.info["ref_energy"]
-                ref_molecule_e = molecule.info["ref_energy"]
-                ref_ads_energy = compute_adsorption_energy(
-                    ref_surface_e, ref_mol_surf_e, ref_molecule_e
-                )
+                ref_ads_energy = mol_surface.info["ref_adsorption_energy"]
                 results["ref"].append(ref_ads_energy)
 
-            # Write structures in order
+            # Write molecule-surface structure to app data
             structs_dir = OUT_PATH / model_name
             structs_dir.mkdir(parents=True, exist_ok=True)
-            write(structs_dir / f"{system_path.stem}.xyz", structs)
+            write(structs_dir / f"{system_path.stem}.xyz", mol_surface)
 
         ref_stored = True
     return results

mlip_testing/app/surfaces/S24/app_S24.py

@@ -35,8 +35,7 @@ def register_callbacks(self) -> None:
         # Assets dir will be parent directory
         structs = [
             f"assets/surfaces/S24/{list(MODELS.keys())[0]}/{struct_file.stem}.xyz"
-            for struct_file in structs_dir.glob("*.xyz")
-            if struct_file.name != "s24_mol_surface_atoms.xyz"
+            for struct_file in sorted(structs_dir.glob("*.xyz"))
         ]
 
         plot_from_table_column(
@@ -79,7 +78,7 @@ def get_app() -> S24App:
 
 if __name__ == "__main__":
     # Create Dash app
-    full_app = Dash(__name__, assets_folder=DATA_PATH.parent)
+    full_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent)
 
     # Construct layout and register callbacks
     s24_app = get_app()

mlip_testing/calcs/surfaces/S24/calc_S24.py

@@ -78,52 +78,53 @@ def run(self):
         data = get_benchmark_data("s24.zip") / "s24/s24_data.extxyz"
         atoms_list = read(data, ":")
 
-        mol_surface_list = []
-
         for i in tqdm(
             range(0, len(atoms_list), 3),
             desc=f"Processing Triplets for model: {self.model_name}",
         ):
             surface = atoms_list[i]
             mol_surface = atoms_list[i + 1]
             molecule = atoms_list[i + 2]
-            surface_formula = surface.get_chemical_formula()
-            molecule_formula = molecule.get_chemical_formula()
+            
+            # Create sys_id
+            sys_id = f"{(i // 3) + 1:03d}"
 
             # Store reference energies
             surface.info["ref_energy"] = surface.get_potential_energy()
             mol_surface.info["ref_energy"] = mol_surface.get_potential_energy()
             molecule.info["ref_energy"] = molecule.get_potential_energy()
 
             # Store system information
+            surface_formula = surface.get_chemical_formula()
+            molecule_formula = molecule.get_chemical_formula()
             system_name = f"{surface_formula}-{molecule_formula}"
-            surface.info["system_name"] = system_name
+            
+            mol_surface.info["sys_id"] = sys_id
             mol_surface.info["system_name"] = system_name
-            molecule.info["system_name"] = system_name
-
-            # Store structure types
-            surface.info["structure_type"] = "surface"
-            mol_surface.info["structure_type"] = "mol_surface"
-            molecule.info["structure_type"] = "molecule"
 
             # Evaluate with the model
             triplet = [surface, mol_surface, molecule]
             self.evaluate_energies(triplet, calc)
 
-            mol_surface_list.append(mol_surface)
-
-        # Write structures to output
-        write_dir = OUT_PATH / self.model_name
-        write_dir.mkdir(parents=True, exist_ok=True)
-        write(write_dir / "s24_mol_surface_atoms.xyz", mol_surface_list)
-
-        # Write all structures organized by system
-        for i in range(0, len(atoms_list), 3):
-            surface = atoms_list[i]
-            mol_surface = atoms_list[i + 1]
-            molecule = atoms_list[i + 2]
-            system_name = surface.info["system_name"]
-            write(write_dir / f"{system_name}.xyz", [surface, mol_surface, molecule])
+            # Calculate and store adsorption energies
+            surface_e = surface.get_potential_energy()
+            mol_surf_e = mol_surface.get_potential_energy()
+            molecule_e = molecule.get_potential_energy()
+            pred_ads_energy = self.compute_adsorption_energy(surface_e, mol_surf_e, molecule_e)
+            
+            ref_surface_e = surface.info["ref_energy"]
+            ref_mol_surf_e = mol_surface.info["ref_energy"]
+            ref_molecule_e = molecule.info["ref_energy"]
+            ref_ads_energy = self.compute_adsorption_energy(ref_surface_e, ref_mol_surf_e, ref_molecule_e)
+            
+            # Store adsorption energies in mol_surface
+            mol_surface.info["adsorption_energy"] = pred_ads_energy
+            mol_surface.info["ref_adsorption_energy"] = ref_ads_energy
+
+            # Save individual molecule-surface structure using sys_id
+            write_dir = OUT_PATH / self.model_name
+            write_dir.mkdir(parents=True, exist_ok=True)
+            write(write_dir / f"{sys_id}.xyz", mol_surface)
 
 
 def build_project(repro: bool = False) -> None: