Implement structure optimization protocol

atomistics/calculators/ase.py

@@ -7,17 +7,42 @@
 if TYPE_CHECKING:
     from ase import Atoms
     from ase.calculators.calculator import Calculator as ASECalculator
+    from ase.optimize.optimize import Optimizer
     from atomistics.calculators.interface import TaskName
 
 
+def ase_optimize_structure(
+    structure, ase_calculator, ase_optimizer, ase_optimizer_kwargs
+):
+    structure_optimized = structure.copy()
+    structure_optimized.calc = ase_calculator
+    ase_optimizer_obj = ase_optimizer(structure_optimized)
+    ase_optimizer_obj.run(**ase_optimizer_kwargs)
+    return structure_optimized
+
+
 @as_task_dict_evaluator
 def evaluate_with_ase(
-    structure: Atoms, tasks: list[TaskName], ase_calculator: ASECalculator
+    structure: Atoms,
+    tasks: list[TaskName],
+    ase_calculator: ASECalculator,
+    ase_optimizer: Optimizer = None,
+    ase_optimizer_kwargs: dict = None,
 ):
-    structure.calc = ase_calculator
     results = {}
-    if "calc_energy" in tasks:
-        results["energy"] = structure.get_potential_energy()
-    if "calc_forces" in tasks:
-        results["forces"] = structure.get_forces()
+    if "optimize_positions" in tasks:
+        results["structure_with_optimized_positions"] = ase_optimize_structure(
+            structure=structure,
+            ase_calculator=ase_calculator,
+            ase_optimizer=ase_optimizer,
+            ase_optimizer_kwargs=ase_optimizer_kwargs,
+        )
+    elif "calc_energy" in tasks or "calc_forces" in tasks:
+        structure.calc = ase_calculator
+        if "calc_energy" in tasks:
+            results["energy"] = structure.get_potential_energy()
+        if "calc_forces" in tasks:
+            results["forces"] = structure.get_forces()
+    else:
+        raise ValueError("The ASE calculator does not implement:", tasks)
     return results

atomistics/calculators/interface.py

@@ -17,11 +17,15 @@ def __str__(self):
 class TaskEnum(StrEnum):
     calc_energy = "calc_energy"
     calc_forces = "calc_forces"
+    optimize_positions = "optimize_positions"
+    optimize_positions_and_volume = "optimize_positions_and_volume"
 
 
 class TaskOutputEnum(Enum):
     energy = "calc_energy"
     forces = "calc_forces"
+    structure_with_optimized_positions = "optimize_positions"
+    structure_with_optimized_positions_and_volume = "optimize_positions_and_volume"
 
 
 if TYPE_CHECKING:

atomistics/calculators/lammps.py

@@ -17,26 +17,72 @@
     from atomistics.calculators.interface import TaskName
 
 
-LAMMPS_INPUT_TEMPLATE = """\
+LAMMPS_STATIC_RUN_INPUT_TEMPLATE = """\
 thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol
 thermo_modify format float %20.15g
 thermo 100
 run 0"""
 
 
+LAMMPS_MINIMIZE_POSITIONS_INPUT_TEMPLATE = """\
+thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol
+thermo_modify format float %20.15g
+thermo 10
+min_style cg
+minimize 0.0 0.0001 100000 10000000"""
+
+
+LAMMPS_MINIMIZE_POSITIONS_AND_VOLUME_INPUT_TEMPLATE = """\
+fix ensemble all box/relax iso 0.0
+thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol
+thermo_modify format float %20.15g
+thermo 10
+min_style cg
+minimize 0.0 0.0001 100000 10000000"""
+
+
 @as_task_dict_evaluator
 def evaluate_with_lammps_library(
     structure: Atoms,
     tasks: list[TaskName],
     potential_dataframe: DataFrame,
     lmp: LammpsASELibrary,
 ):
-    lmp = _run_simulation(structure, potential_dataframe, LAMMPS_INPUT_TEMPLATE, lmp)
     results = {}
-    if "calc_energy" in tasks:
-        results["energy"] = lmp.interactive_energy_pot_getter()
-    if "calc_forces" in tasks:
-        results["forces"] = lmp.interactive_forces_getter()
+    if "optimize_positions_and_volume" in tasks:
+        lmp = _run_simulation(
+            structure=structure,
+            potential_dataframe=potential_dataframe,
+            input_template=LAMMPS_MINIMIZE_POSITIONS_AND_VOLUME_INPUT_TEMPLATE,
+            lmp=lmp,
+        )
+        structure_copy = structure.copy()
+        structure_copy.set_cell(lmp.interactive_cells_getter(), scale_atoms=True)
+        structure_copy.positions = lmp.interactive_positions_getter()
+        results["structure_with_optimized_positions_and_volume"] = structure_copy
+    elif "optimize_positions" in tasks:
+        lmp = _run_simulation(
+            structure=structure,
+            potential_dataframe=potential_dataframe,
+            input_template=LAMMPS_MINIMIZE_POSITIONS_INPUT_TEMPLATE,
+            lmp=lmp,
+        )
+        structure_copy = structure.copy()
+        structure_copy.positions = lmp.interactive_positions_getter()
+        results["structure_with_optimized_positions"] = structure_copy
+    elif "calc_energy" in tasks or "calc_forces" in tasks:
+        lmp = _run_simulation(
+            structure=structure,
+            potential_dataframe=potential_dataframe,
+            input_template=LAMMPS_STATIC_RUN_INPUT_TEMPLATE,
+            lmp=lmp,
+        )
+        if "calc_energy" in tasks:
+            results["energy"] = lmp.interactive_energy_pot_getter()
+        if "calc_forces" in tasks:
+            results["forces"] = lmp.interactive_forces_getter()
+    else:
+        raise ValueError("The LAMMPS calculator does not implement:", tasks)
     lmp.interactive_lib_command("clear")
     return results
 

atomistics/calculators/wrapper.py

@@ -33,6 +33,8 @@ def _convert_task_dict(old_task_dict: dict[TaskName, dict[str, Atoms]]) -> TaskD
     """
     task_dict = {}
     for method_name, subdict in old_task_dict.items():
+        if not isinstance(subdict, dict):
+            subdict = {"label_hidden": subdict}
         for label, structure in subdict.items():
             try:
                 task_dict[label][1].append(method_name)
@@ -72,10 +74,13 @@ def evaluate_with_calculator(
             output = calculate(structure, tasks, *calculate_args, **calculate_kwargs)
             for task_name in tasks:
                 result_name = TaskOutputEnum(task_name).name
-                try:
-                    results_dict[result_name][label] = output[result_name]
-                except KeyError:
-                    results_dict[result_name] = {label: output[result_name]}
+                if label != "label_hidden":
+                    try:
+                        results_dict[result_name][label] = output[result_name]
+                    except KeyError:
+                        results_dict[result_name] = {label: output[result_name]}
+                else:
+                    results_dict[result_name] = output[result_name]
         return results_dict
 
     return evaluate_with_calculator

atomistics/workflows/structure_optimization/workflow.py

@@ -0,0 +1,6 @@
+def optimize_positions_and_volume(structure):
+    return {"optimize_positions_and_volume": structure}
+
+
+def optimize_positions(structure):
+    return {"optimize_positions": structure}

tests/test_elastic_emt.py

@@ -9,15 +9,15 @@
 
 class TestElastic(unittest.TestCase):
     def test_calc_elastic(self):
-        calculator = calculator = ElasticMatrixWorkflow(
+        calculator = ElasticMatrixWorkflow(
             structure=bulk("Al", a=4.0, cubic=True),
             num_of_point=5,
             eps_range=0.005,
             sqrt_eta=True,
             fit_order=2
         )
-        structure_dict = calculator.generate_structures()
-        result_dict = evaluate_with_ase(structure_dict, ase_calculator=EMT())
+        task_dict = calculator.generate_structures()
+        result_dict = evaluate_with_ase(task_dict=task_dict, ase_calculator=EMT())
         elastic_dict = calculator.analyse_structures(output_dict=result_dict)
         self.assertTrue(np.isclose(elastic_dict["C"][0, 0], 52.62435421))
         self.assertTrue(np.isclose(elastic_dict["C"][0, 1], 32.6743838))

tests/test_elastic_gpaw.py

@@ -25,9 +25,9 @@ def test_calc_elastic(self):
             sqrt_eta=True,
             fit_order=2
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_ase(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             ase_calculator=GPAW(
                 xc="PBE",
                 mode=PW(300),

tests/test_elastic_lammps.py

@@ -5,6 +5,7 @@
 import unittest
 
 from atomistics.workflows.elastic.workflow import ElasticMatrixWorkflow
+from atomistics.workflows.structure_optimization.workflow import optimize_positions_and_volume
 
 try:
     from atomistics.calculators.lammps import (
@@ -23,25 +24,30 @@ class TestElastic(unittest.TestCase):
     def test_calc_elastic(self):
         potential = '1999--Mishin-Y--Al--LAMMPS--ipr1'
         resource_path = os.path.join(os.path.dirname(__file__), "static", "lammps")
-        structure = bulk("Al", a=4.05, cubic=True)
+        structure = bulk("Al", cubic=True)
         df_pot = get_potential_dataframe(
             structure=structure,
             resource_path=resource_path
         )
         df_pot_selected = df_pot[df_pot.Name == potential].iloc[0]
+        task_dict = optimize_positions_and_volume(structure=structure)
+        result_dict = evaluate_with_lammps(
+            task_dict=task_dict,
+            potential_dataframe=df_pot_selected,
+        )
         calculator = ElasticMatrixWorkflow(
-            structure=structure,
+            structure=result_dict["structure_with_optimized_positions_and_volume"],
             num_of_point=5,
             eps_range=0.005,
             sqrt_eta=True,
             fit_order=2
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_lammps(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             potential_dataframe=df_pot_selected,
         )
         elastic_dict = calculator.analyse_structures(output_dict=result_dict)
         self.assertTrue(np.isclose(elastic_dict["C"][0, 0], 114.10393023))
         self.assertTrue(np.isclose(elastic_dict["C"][0, 1], 60.51098897))
-        self.assertTrue(np.isclose(elastic_dict["C"][3, 3], 51.23931149))
+        self.assertTrue(np.isclose(elastic_dict["C"][3, 3], 51.23853765))

tests/test_evcurve_abinit.py

@@ -43,9 +43,9 @@ def test_calc_evcurve(self):
             axes=['x', 'y', 'z'],
             strains=None,
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_ase(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             ase_calculator=Abinit(
                 label='abinit_evcurve',
                 nbands=32,

tests/test_evcurve_emt.py

@@ -18,8 +18,8 @@ def test_calc_evcurve(self):
             axes=['x', 'y', 'z'],
             strains=None,
         )
-        structure_dict = calculator.generate_structures()
-        result_dict = evaluate_with_ase(task_dict=structure_dict, ase_calculator=EMT())
+        task_dict = calculator.generate_structures()
+        result_dict = evaluate_with_ase(task_dict=task_dict, ase_calculator=EMT())
         fit_dict = calculator.analyse_structures(output_dict=result_dict)
         self.assertTrue(np.isclose(fit_dict['volume_eq'], 63.72615218844302))
         self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 39.544084907317895))

tests/test_evcurve_gpaw.py

@@ -28,9 +28,9 @@ def test_calc_evcurve(self):
             axes=['x', 'y', 'z'],
             strains=None,
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_ase(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             ase_calculator=GPAW(
                 xc="PBE",
                 mode=PW(300),

tests/test_evcurve_lammps.py

@@ -5,6 +5,7 @@
 import unittest
 
 from atomistics.workflows.evcurve.workflow import EnergyVolumeCurveWorkflow
+from atomistics.workflows.structure_optimization.workflow import optimize_positions_and_volume
 
 
 try:
@@ -24,14 +25,19 @@ class TestEvCurve(unittest.TestCase):
     def test_calc_evcurve(self):
         potential = '1999--Mishin-Y--Al--LAMMPS--ipr1'
         resource_path = os.path.join(os.path.dirname(__file__), "static", "lammps")
-        structure = bulk("Al", a=4.05, cubic=True)
+        structure = bulk("Al", cubic=True)
         df_pot = get_potential_dataframe(
             structure=structure,
             resource_path=resource_path
         )
         df_pot_selected = df_pot[df_pot.Name == potential].iloc[0]
+        task_dict = optimize_positions_and_volume(structure=structure)
+        result_dict = evaluate_with_lammps(
+            task_dict=task_dict,
+            potential_dataframe=df_pot_selected,
+        )
         calculator = EnergyVolumeCurveWorkflow(
-            structure=structure,
+            structure=result_dict["structure_with_optimized_positions_and_volume"],
             num_points=11,
             fit_type='polynomial',
             fit_order=3,
@@ -47,4 +53,4 @@ def test_calc_evcurve(self):
         fit_dict = calculator.analyse_structures(output_dict=result_dict)
         self.assertTrue(np.isclose(fit_dict['volume_eq'], 66.43019853103964))
         self.assertTrue(np.isclose(fit_dict['bulkmodul_eq'], 77.7250135953191))
-        self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.279502459079921))
+        self.assertTrue(np.isclose(fit_dict['b_prime_eq'], 1.2795467367276832))

tests/test_evcurve_qe.py

@@ -46,9 +46,9 @@ def test_calc_evcurve(self):
             axes=['x', 'y', 'z'],
             strains=None,
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_ase(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             ase_calculator=Espresso(
                 pseudopotentials=pseudopotentials,
                 tstress=True,

tests/test_evcurve_siesta.py

@@ -47,9 +47,9 @@ def test_calc_evcurve(self):
             axes=['x', 'y', 'z'],
             strains=None,
         )
-        structure_dict = calculator.generate_structures()
+        task_dict = calculator.generate_structures()
         result_dict = evaluate_with_ase(
-            task_dict=structure_dict,
+            task_dict=task_dict,
             ase_calculator=Siesta(
                 label="siesta",
                 xc="PBE",

tests/test_optimize_positions_emt.py

@@ -0,0 +1,31 @@
+import numpy as np
+from ase.build import bulk
+from ase.calculators.emt import EMT
+from ase.optimize import BFGS
+import unittest
+
+from atomistics.calculators.ase import evaluate_with_ase
+from atomistics.workflows.structure_optimization.workflow import optimize_positions
+
+
+class TestOptimizePositions(unittest.TestCase):
+    def test_optimize_positions(self):
+        structure = bulk("Al", a=4.0, cubic=True)
+        positions_before_displacement = structure.positions.copy()
+        structure.positions[0] += [0.01, 0.01, 0.01]
+        task_dict = optimize_positions(structure=structure)
+        result_dict = evaluate_with_ase(
+            task_dict=task_dict,
+            ase_calculator=EMT(),
+            ase_optimizer=BFGS,
+            ase_optimizer_kwargs={"fmax": 0.00005}
+        )
+        structure_optimized = result_dict["structure_with_optimized_positions"]
+        self.assertTrue(
+            all(np.isclose(
+                positions_before_displacement,
+                structure_optimized.positions-structure_optimized.positions[0],
+                atol=1e-5
+            ).flatten())
+        )
+