Add thermal expansion example for LAMMPS

atomistics/calculators/__init__.py

@@ -11,6 +11,7 @@
         calc_energy_and_forces_with_lammps,
         calc_energy_with_lammps,
         calc_forces_with_lammps,
+        calc_molecular_dynamics_thermal_expansion_with_lammps,
         evaluate_with_lammps,
         evaluate_with_lammps_library,
         get_potential_dataframe,

atomistics/calculators/interface.py

@@ -19,13 +19,17 @@ class TaskEnum(StrEnum):
     calc_forces = "calc_forces"
     optimize_positions = "optimize_positions"
     optimize_positions_and_volume = "optimize_positions_and_volume"
+    calc_molecular_dynamics_thermal_expansion = (
+        "calc_molecular_dynamics_thermal_expansion"
+    )
 
 
 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"
+    volume_over_temperature = "calc_molecular_dynamics_thermal_expansion"
 
 
 if TYPE_CHECKING:

atomistics/calculators/lammps/__init__.py

@@ -2,6 +2,7 @@
     calc_energy_and_forces_with_lammps,
     calc_energy_with_lammps,
     calc_forces_with_lammps,
+    calc_molecular_dynamics_thermal_expansion_with_lammps,
     evaluate_with_lammps,
     evaluate_with_lammps_library,
     optimize_positions_and_volume_with_lammps,

atomistics/calculators/lammps/calculator.py

@@ -1,19 +1,24 @@
 from __future__ import annotations
 
+import numpy as np
 from typing import TYPE_CHECKING
 
 from pylammpsmpi import LammpsASELibrary
 
 from atomistics.calculators.wrapper import as_task_dict_evaluator
 from atomistics.calculators.lammps.helpers import (
     lammps_run,
+    lammps_thermal_expansion_loop,
     lammps_shutdown,
     template_render_minimize,
     template_render_run,
 )
 from atomistics.calculators.lammps.commands import (
     LAMMPS_THERMO_STYLE,
     LAMMPS_THERMO,
+    LAMMPS_ENSEMBLE_NPT,
+    LAMMPS_VELOCITY,
+    LAMMPS_TIMESTEP,
     LAMMPS_MINIMIZE,
     LAMMPS_RUN,
     LAMMPS_MINIMIZE_VOLUME,
@@ -26,7 +31,9 @@
     from atomistics.calculators.interface import TaskName
 
 
-def calc_energy_with_lammps(structure: Atoms, potential_dataframe: DataFrame, lmp=None):
+def calc_energy_with_lammps(
+    structure: Atoms, potential_dataframe: DataFrame, lmp=None, **kwargs
+):
     template_str = LAMMPS_THERMO_STYLE + "\n" + LAMMPS_THERMO + "\n" + LAMMPS_RUN
     lmp_instance = lammps_run(
         structure=structure,
@@ -37,13 +44,16 @@ def calc_energy_with_lammps(structure: Atoms, potential_dataframe: DataFrame, lm
             run=0,
         ),
         lmp=lmp,
+        **kwargs,
     )
     energy_pot = lmp_instance.interactive_energy_pot_getter()
     lammps_shutdown(lmp_instance=lmp_instance, close_instance=lmp is None)
     return energy_pot
 
 
-def calc_forces_with_lammps(structure: Atoms, potential_dataframe: DataFrame, lmp=None):
+def calc_forces_with_lammps(
+    structure: Atoms, potential_dataframe: DataFrame, lmp=None, **kwargs
+):
     template_str = LAMMPS_THERMO_STYLE + "\n" + LAMMPS_THERMO + "\n" + LAMMPS_RUN
     lmp_instance = lammps_run(
         structure=structure,
@@ -54,14 +64,15 @@ def calc_forces_with_lammps(structure: Atoms, potential_dataframe: DataFrame, lm
             run=0,
         ),
         lmp=lmp,
+        **kwargs,
     )
     forces = lmp_instance.interactive_forces_getter()
     lammps_shutdown(lmp_instance=lmp_instance, close_instance=lmp is None)
     return forces
 
 
 def calc_energy_and_forces_with_lammps(
-    structure, potential_dataframe: DataFrame, lmp=None
+    structure, potential_dataframe: DataFrame, lmp=None, **kwargs
 ):
     template_str = LAMMPS_THERMO_STYLE + "\n" + LAMMPS_THERMO + "\n" + LAMMPS_RUN
     lmp_instance = lammps_run(
@@ -73,6 +84,7 @@ def calc_energy_and_forces_with_lammps(
             run=0,
         ),
         lmp=lmp,
+        **kwargs,
     )
     energy_pot = lmp_instance.interactive_energy_pot_getter()
     forces = lmp_instance.interactive_forces_getter()
@@ -90,6 +102,7 @@ def optimize_positions_and_volume_with_lammps(
     maxeval=10000000,
     thermo=10,
     lmp=None,
+    **kwargs,
 ):
     template_str = (
         LAMMPS_MINIMIZE_VOLUME
@@ -113,6 +126,7 @@ def optimize_positions_and_volume_with_lammps(
             thermo=thermo,
         ),
         lmp=lmp,
+        **kwargs,
     )
     structure_copy = structure.copy()
     structure_copy.set_cell(lmp_instance.interactive_cells_getter(), scale_atoms=True)
@@ -131,6 +145,7 @@ def optimize_positions_with_lammps(
     maxeval=10000000,
     thermo=10,
     lmp=None,
+    **kwargs,
 ):
     template_str = LAMMPS_THERMO_STYLE + "\n" + LAMMPS_THERMO + "\n" + LAMMPS_MINIMIZE
     lmp_instance = lammps_run(
@@ -146,13 +161,65 @@ def optimize_positions_with_lammps(
             thermo=thermo,
         ),
         lmp=lmp,
+        **kwargs,
     )
     structure_copy = structure.copy()
     structure_copy.positions = lmp_instance.interactive_positions_getter()
     lammps_shutdown(lmp_instance=lmp_instance, close_instance=lmp is None)
     return structure_copy
 
 
+def calc_molecular_dynamics_thermal_expansion_with_lammps(
+    structure,
+    potential_dataframe,
+    Tstart=15,
+    Tstop=1500,
+    Tstep=5,
+    Tdamp=0.1,
+    run=100,
+    thermo=100,
+    timestep=0.001,
+    Pstart=0.0,
+    Pstop=0.0,
+    Pdamp=1.0,
+    seed=4928459,
+    dist="gaussian",
+    lmp=None,
+    **kwargs,
+):
+    init_str = (
+        LAMMPS_THERMO_STYLE
+        + "\n"
+        + LAMMPS_TIMESTEP
+        + "\n"
+        + LAMMPS_THERMO
+        + "\n"
+        + LAMMPS_VELOCITY
+        + "\n"
+    )
+    run_str = LAMMPS_ENSEMBLE_NPT + "\n" + LAMMPS_RUN
+    temperature_lst = np.arange(Tstart, Tstop + Tstep, Tstep).tolist()
+    volume_md_lst = lammps_thermal_expansion_loop(
+        structure=structure,
+        potential_dataframe=potential_dataframe,
+        init_str=init_str,
+        run_str=run_str,
+        temperature_lst=temperature_lst,
+        run=run,
+        thermo=thermo,
+        timestep=timestep,
+        Tdamp=Tdamp,
+        Pstart=Pstart,
+        Pstop=Pstop,
+        Pdamp=Pdamp,
+        seed=seed,
+        dist=dist,
+        lmp=lmp,
+        **kwargs,
+    )
+    return temperature_lst, volume_md_lst
+
+
 @as_task_dict_evaluator
 def evaluate_with_lammps_library(
     structure: Atoms,
@@ -178,6 +245,17 @@ def evaluate_with_lammps_library(
             lmp=lmp,
             **lmp_optimizer_kwargs,
         )
+    elif "calc_molecular_dynamics_thermal_expansion" in tasks:
+        (
+            temperature_lst,
+            volume_md_lst,
+        ) = calc_molecular_dynamics_thermal_expansion_with_lammps(
+            structure=structure,
+            potential_dataframe=potential_dataframe,
+            lmp=lmp,
+            **lmp_optimizer_kwargs,
+        )
+        results["volume_over_temperature"] = (temperature_lst, volume_md_lst)
     elif "calc_energy" in tasks and "calc_forces" in tasks:
         results["energy"], results["forces"] = calc_energy_and_forces_with_lammps(
             structure=structure,

atomistics/calculators/lammps/commands.py

@@ -15,3 +15,12 @@
 
 
 LAMMPS_MINIMIZE_VOLUME = "fix ensemble all box/relax iso 0.0"
+
+
+LAMMPS_TIMESTEP = "timestep {{timestep}}"
+
+
+LAMMPS_VELOCITY = "velocity all create $(2 * {{ temp }}) {{seed}} dist {{dist}}"
+
+
+LAMMPS_ENSEMBLE_NPT = "fix ensemble all npt temp {{Tstart}} {{Tstop}} {{Tdamp}} iso {{Pstart}} {{Pstop}} {{Pdamp}}"

atomistics/calculators/lammps/helpers.py

@@ -1,19 +1,8 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
-
 from jinja2 import Template
-import pandas
 from pylammpsmpi import LammpsASELibrary
 
-from atomistics.calculators.wrapper import as_task_dict_evaluator
-from atomistics.calculators.lammps.commands import (
-    LAMMPS_THERMO_STYLE,
-    LAMMPS_THERMO,
-    LAMMPS_MINIMIZE,
-    LAMMPS_RUN,
-    LAMMPS_MINIMIZE_VOLUME,
-)
 from atomistics.calculators.lammps.potential import validate_potential_dataframe
 
 
@@ -47,12 +36,12 @@ def template_render_run(
     )
 
 
-def lammps_run(structure, potential_dataframe, input_template, lmp=None):
+def lammps_run(structure, potential_dataframe, input_template, lmp=None, **kwargs):
     potential_dataframe = validate_potential_dataframe(
         potential_dataframe=potential_dataframe
     )
     if lmp is None:
-        lmp = LammpsASELibrary()
+        lmp = LammpsASELibrary(**kwargs)
 
     # write structure to LAMMPS
     lmp.interactive_structure_setter(
@@ -75,6 +64,56 @@ def lammps_run(structure, potential_dataframe, input_template, lmp=None):
     return lmp
 
 
+def lammps_thermal_expansion_loop(
+    structure,
+    potential_dataframe,
+    init_str,
+    run_str,
+    temperature_lst,
+    run=100,
+    thermo=100,
+    timestep=0.001,
+    Tdamp=0.1,
+    Pstart=0.0,
+    Pstop=0.0,
+    Pdamp=1.0,
+    seed=4928459,
+    dist="gaussian",
+    lmp=None,
+    **kwargs,
+):
+    lmp_instance = lammps_run(
+        structure=structure,
+        potential_dataframe=potential_dataframe,
+        input_template=Template(init_str).render(
+            thermo=thermo,
+            temp=temperature_lst[0],
+            timestep=timestep,
+            seed=seed,
+            dist=dist,
+        ),
+        lmp=lmp,
+        **kwargs,
+    )
+
+    volume_md_lst = []
+    for temp in temperature_lst:
+        run_str_rendered = Template(run_str).render(
+            run=run,
+            Tstart=temp - 5,
+            Tstop=temp,
+            Tdamp=Tdamp,
+            Pstart=Pstart,
+            Pstop=Pstop,
+            Pdamp=Pdamp,
+        )
+        for l in run_str_rendered.split("\n"):
+            lmp_instance.interactive_lib_command(l)
+        volume_md_lst.append(lmp_instance.interactive_volume_getter())
+    lammps_shutdown(lmp_instance=lmp_instance, close_instance=lmp is None)
+    return volume_md_lst
+
+
 def lammps_shutdown(lmp_instance, close_instance=True):
     lmp_instance.interactive_lib_command("clear")
     if close_instance:

atomistics/workflows/__init__.py

@@ -1,6 +1,9 @@
 from atomistics.workflows.elastic.workflow import ElasticMatrixWorkflow
 from atomistics.workflows.evcurve.workflow import EnergyVolumeCurveWorkflow
 from atomistics.workflows.langevin.workflow import LangevinWorkflow
+from atomistics.workflows.molecular_dynamics.workflow import (
+    calc_molecular_dynamics_thermal_expansion,
+)
 from atomistics.workflows.structure_optimization.workflow import (
     optimize_positions,
     optimize_positions_and_volume,

atomistics/workflows/molecular_dynamics/workflow.py

@@ -0,0 +1,2 @@
+def calc_molecular_dynamics_thermal_expansion(structure):
+    return {"calc_molecular_dynamics_thermal_expansion": structure}