feat: add dice input generation and player routines for NVT/NPT simulations

diceplayer/config/dice_config.py

@@ -29,6 +29,13 @@ class DiceConfig(BaseModel):
     upbuf: int = Field(
         360, description="Buffer size for the potential energy calculations"
     )
+    irdf: int = Field(
+        0,
+        description="Flag for calculating radial distribution functions (0: no, 1: yes)",
+    )
+    vstep: int = Field(
+        5000, description="Frequency of volume change moves in NPT simulations"
+    )
     combrule: Literal["+", "*"] = Field(
         "*", description="Combination rule for the Lennard-Jones potential"
     )

diceplayer/config/player_config.py

@@ -20,7 +20,20 @@ class RoutineType(str, Enum):
 
 class PlayerConfig(BaseModel):
     """
-    Data Transfer Object for the player configuration.
+    Configuration for DICEPlayer simulations.
+
+    Attributes:
+        type: Type of simulation to perform (charge, geometry, or both).
+        max_cyc: Maximum number of cycles for the geometry optimization.
+        switch_cyc: Cycle at which to switch from charge to geometry optimization (if type is "both").
+        mem: Memory configuration for QM calculations.
+        nprocs: Number of processors to use for QM calculations.
+        ncores: Number of cores to use for QM calculations.
+        dice: Configuration parameters specific to DICE simulations.
+        gaussian: Configuration parameters specific to Gaussian calculations.
+        altsteps: Number of steps for the alternate simulation (default: 20000).
+        geoms_file: File name for the geometries output (default: "geoms.xyz").
+        simulation_dir: Directory name for the simulation files (default: "simfiles").
     """
 
     model_config = ConfigDict(

diceplayer/dice/__init__.py

@@ -0,0 +1,174 @@
+"""
+DICE Monte Carlo Simulation Interface
+=====================================
+
+This package provides utilities for configuring and running simulations with
+the DICE Monte Carlo molecular simulation program.
+
+DICE performs statistical sampling of molecular systems using the Metropolis
+Monte Carlo algorithm. Simulations are defined by text input files containing
+keywords that control the thermodynamic ensemble, system composition, and
+simulation parameters.
+
+---------------------------------------------------------------------
+
+Simulation Ensembles
+--------------------
+
+DICE supports multiple statistical ensembles.
+
+NVT
+    Canonical ensemble where the following properties remain constant:
+
+    - N: number of molecules
+    - V: system volume
+    - T: temperature
+
+    The system density is fixed and the simulation box volume does not change
+    during the simulation.
+
+NPT
+    Isothermal–isobaric ensemble where the following properties remain constant:
+
+    - N: number of molecules
+    - P: pressure
+    - T: temperature
+
+    The simulation box volume is allowed to fluctuate in order to maintain the
+    target pressure.
+
+---------------------------------------------------------------------
+
+Simulation Stages
+-----------------
+
+Simulations are typically executed in multiple stages.
+
+Thermalization (TER)
+    Initial phase where the system relaxes to the desired thermodynamic
+    conditions. Molecular configurations stabilize and the system reaches
+    equilibrium.
+
+    During this stage statistical properties are **not accumulated**.
+
+Production / Equilibration (EQ)
+    Main sampling phase after the system has equilibrated.
+
+    Statistical properties such as energies, densities, and radial
+    distribution functions are collected and configurations may be saved
+    for later analysis.
+
+---------------------------------------------------------------------
+
+Typical Simulation Pipeline
+---------------------------
+
+Two common execution workflows are used.
+
+NVT Simulation
+    Used when the system density is known.
+
+    1. NVT.ter  → thermalization at constant density
+    2. NVT.eq   → production sampling
+
+NPT Simulation
+    Used when the equilibrium density is unknown.
+
+    1. NVT.ter  → initial thermalization at approximate density
+    2. NPT.ter  → pressure relaxation (volume adjustment)
+    3. NPT.eq   → production sampling at target pressure
+
+---------------------------------------------------------------------
+
+DICE Input Keywords
+-------------------
+
+The following keywords are used in the generated input files.
+
+title
+    Descriptive title printed in the simulation output.
+
+ncores
+    Number of CPU cores used by the DICE executable.
+
+ljname
+    File containing Lennard-Jones parameters and molecular topology.
+
+outname
+    Prefix used for simulation output files.
+
+nmol
+    Number of molecules of each species in the system.
+
+dens
+    System density (g/cm³). Used only in NVT simulations or for
+    initialization of NPT runs.
+
+press
+    Target pressure used in NPT simulations.
+
+temp
+    Simulation temperature.
+
+nstep
+    Number of Monte Carlo cycles executed in the simulation stage.
+
+init
+    Defines how the simulation initializes molecular coordinates.
+
+    yes
+        Random initial configuration.
+
+    no
+        Continue from a previous configuration.
+
+    yesreadxyz
+        Read coordinates from a previously saved XYZ configuration.
+
+vstep
+    Frequency of volume-change moves in NPT simulations.
+
+mstop
+    Molecule displacement control flag used internally by DICE.
+
+accum
+    Enables or disables accumulation of statistical averages.
+
+iprint
+    Frequency of simulation information printed to the output.
+
+isave
+    Frequency at which configurations are written to trajectory files.
+
+irdf
+    Controls calculation of radial distribution functions.
+
+seed
+    Random number generator seed used by the Monte Carlo algorithm.
+
+upbuf
+    Buffer size parameter used internally by DICE during thermalization.
+
+---------------------------------------------------------------------
+
+Output Files
+------------
+
+Important output files produced during the simulation include:
+
+phb.xyz
+    XYZ trajectory containing sampled molecular configurations.
+
+last.xyz
+    Final configuration of the simulation, often used as the starting
+    configuration for the next simulation cycle.
+
+---------------------------------------------------------------------
+
+References
+----------
+
+DICE is a Monte Carlo molecular simulation program developed primarily
+by researchers at the University of São Paulo (USP) for studying liquids,
+solutions, and solvation phenomena.
+"""

diceplayer/dice/dice_input.py

@@ -0,0 +1,186 @@
+from diceplayer.config import PlayerConfig
+
+from typing_extensions import Self
+
+import random
+from abc import ABC
+from dataclasses import dataclass, fields
+from typing import Any, Sequence, TextIO
+
+
+@dataclass(slots=True)
+class BaseConfig(ABC):
+    ncores: int
+    ljname: str
+    outname: str
+    nmol: Sequence[int]
+    temp: float
+    seed: int
+    isave: int
+    press: float = 1.0
+
+    def write_dice_config(self, io_writer: TextIO) -> None:
+        for field in fields(self):
+            key = field.name
+            value = getattr(self, key)
+
+            if value is None:
+                continue
+
+            io_writer.write(f"{key} = {self._serialize_value(value)}\n")
+
+        io_writer.write("$end\n")
+
+    @classmethod
+    def from_config(cls, config: PlayerConfig, **kwargs) -> Self:
+        base_fields = cls._extract_base_fields(config)
+        return cls(**base_fields, **kwargs)
+
+    @staticmethod
+    def _extract_base_fields(config: PlayerConfig) -> dict[str, Any]:
+        seed: int
+        if config.dice.seed is not None:
+            seed = config.dice.seed
+        else:
+            seed = random.randint(0, 2**32 - 1)
+
+        return dict(
+            ncores=config.ncores,
+            ljname=config.dice.ljname,
+            outname=config.dice.outname,
+            nmol=config.dice.nmol,
+            temp=config.dice.temp,
+            seed=seed,
+            isave=config.dice.isave,
+            press=config.dice.press,
+        )
+
+    @staticmethod
+    def _get_nstep(config: PlayerConfig, idx: int) -> int:
+        if len(config.dice.nstep) > idx:
+            return config.dice.nstep[idx]
+        return config.dice.nstep[-1]
+
+    @staticmethod
+    def _serialize_value(value: Any) -> str:
+        if value is None:
+            raise ValueError("DICE configuration cannot serialize None values")
+
+        if isinstance(value, bool):
+            return "yes" if value else "no"
+
+        if isinstance(value, (list, tuple)):
+            return " ".join(str(v) for v in value)
+
+        return str(value)
+
+
+# -----------------------------------------------------
+# NVT BASE
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NVTConfig(BaseConfig):
+    title: str = "Diceplayer Run - NVT"
+    dens: float = 0.0
+    nstep: int = 0
+    vstep: int = 0
+
+
+# -----------------------------------------------------
+# NVT THERMALIZATION
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NVTTerConfig(NVTConfig):
+    title: str = "Diceplayer Run - NVT Thermalization"
+    upbuf: int = 360
+
+    @classmethod
+    def from_config(cls, config: PlayerConfig, **kwargs) -> Self:
+        return super(NVTTerConfig, cls).from_config(
+            config,
+            dens=config.dice.dens,
+            nstep=cls._get_nstep(config, 0),
+            upbuf=config.dice.upbuf,
+            vstep=0,
+            **kwargs,
+        )
+
+
+# -----------------------------------------------------
+# NVT PRODUCTION
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NVTEqConfig(NVTConfig):
+    title: str = "Diceplayer Run - NVT Production"
+    irdf: int = 0
+
+    @classmethod
+    def from_config(cls, config: PlayerConfig, **kwargs) -> Self:
+        return super(NVTEqConfig, cls).from_config(
+            config,
+            dens=config.dice.dens,
+            nstep=cls._get_nstep(config, 1),
+            irdf=config.dice.irdf,
+            vstep=0,
+            **kwargs,
+        )
+
+
+# -----------------------------------------------------
+# NPT BASE
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NPTConfig(BaseConfig):
+    title: str = "Diceplayer Run - NPT"
+    nstep: int = 0
+    vstep: int = 5000
+
+
+# -----------------------------------------------------
+# NPT THERMALIZATION
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NPTTerConfig(NPTConfig):
+    title: str = "Diceplayer Run - NPT Thermalization"
+    dens: float | None = None
+
+    @classmethod
+    def from_config(cls, config: PlayerConfig, **kwargs) -> Self:
+        return super(NPTTerConfig, cls).from_config(
+            config,
+            dens=config.dice.dens,
+            nstep=cls._get_nstep(config, 1),
+            vstep=config.dice.vstep,
+            **kwargs,
+        )
+
+
+# -----------------------------------------------------
+# NPT PRODUCTION
+# -----------------------------------------------------
+
+
+@dataclass(slots=True)
+class NPTEqConfig(NPTConfig):
+    title: str = "Diceplayer Run - NPT Production"
+    dens: float | None = None
+
+    @classmethod
+    def from_config(cls, config: PlayerConfig, **kwargs) -> Self:
+        return super(NPTEqConfig, cls).from_config(
+            config,
+            dens=config.dice.dens,
+            nstep=cls._get_nstep(config, 2),
+            vstep=config.dice.vstep,
+            **kwargs,
+        )

diceplayer/player.py

@@ -1,4 +1,4 @@
-from diceplayer.config.player_config import PlayerConfig
+from diceplayer.config.player_config import PlayerConfig, RoutineType
 from diceplayer.logger import logger
 from diceplayer.state.state_handler import StateHandler
 from diceplayer.state.state_model import StateModel
@@ -14,17 +14,16 @@ class PlayerFlags(TypedDict):
 class Player:
     def __init__(self, config: PlayerConfig):
         self.config = config
+        self._state_handler = StateHandler(config.simulation_dir)
 
     def play(self, **flags: Unpack[PlayerFlags]):
-        state_handler = StateHandler(self.config.simulation_dir)
-
         if not flags["continuation"]:
             logger.info(
                 "Continuation flag is not set. Starting a new simulation and deleting any existing state."
             )
-            state_handler.delete()
+            self._state_handler.delete()
 
-        state = state_handler.get(self.config, force=flags["force"])
+        state = self._state_handler.get(self.config, force=flags["force"])
 
         if state is None:
             state = StateModel.from_config(self.config)
@@ -35,7 +34,36 @@ class Player:
             logger.info(
                 f"Starting cycle {state.current_cycle + 1} of {self.config.max_cyc}."
             )
+
+            step_directory = self.config.simulation_dir / f"{state.current_cycle::02d}"
+            if not step_directory.exists():
+                step_directory.mkdir(parents=True)
+
+            current_routine = self._fetch_current_routine(state.current_cycle)
+            if current_routine == RoutineType.CHARGE:
+                self._charge_opt_routine(state)
+            elif current_routine == RoutineType.GEOMETRY:
+                self._geometry_opt_routine(state)
+            else:
+                logger.error(f"Invalid routine type: {current_routine}")
+                return
+
             state.current_cycle += 1
-            state_handler.save(state)
+            self._state_handler.save(state)
 
         logger.info("Reached maximum number of cycles. Simulation complete.")
+
+    def _fetch_current_routine(self, current_cycle: int) -> RoutineType:
+        if self.config.type != RoutineType.BOTH:
+            return self.config.type
+
+        if current_cycle < self.config.switch_cyc:
+            return RoutineType.CHARGE
+
+        return RoutineType.GEOMETRY
+
+    def _charge_opt_routine(self, state: StateModel) -> None:
+        pass
+
+    def _geometry_opt_routine(self, state: StateModel) -> None:
+        pass

tests/config/test_player_config.py

@@ -1,56 +1,62 @@
 from diceplayer.config.dice_config import DiceConfig
 from diceplayer.config.gaussian_config import GaussianConfig
-from diceplayer.config.player_config import PlayerConfig
+from diceplayer.config.player_config import PlayerConfig, RoutineType
 
 import pytest
 
+from typing import Any
 
-def get_config_dict():
+
+class TestPlayerConfig:
+    @pytest.fixture
+    def dice_payload(self) -> dict[str, Any]:
         return {
-        "opt": True,
-        "mem": 12,
-        "maxcyc": 100,
-        "nprocs": 4,
-        "ncores": 4,
-        "dice": {
             "ljname": "test",
             "outname": "test",
             "dens": 1.0,
             "nmol": [1],
             "nstep": [1, 1],
-        },
+        }
-        "gaussian": {
+
+    @pytest.fixture
+    def gaussian_payload(self) -> dict[str, Any]:
+        return {
             "level": "test",
             "qmprog": "g16",
             "keywords": "test",
-        },
         }
 
-
-class TestPlayerConfig:
     @pytest.fixture
-    def dice_config(self):
+    def player_payload(
-        return DiceConfig(
+        self, dice_payload: dict[str, Any], gaussian_payload: dict[str, Any]
-            ljname="test",
+    ) -> dict[str, Any]:
-            outname="test",
+        return {
-            dens=1.0,
+            "type": "both",
-            nmol=[1],
+            "mem": 12,
-            nstep=[1, 1],
+            "max_cyc": 100,
-        )
+            "switch_cyc": 50,
+            "nprocs": 4,
+            "ncores": 4,
+            "dice": dice_payload,
+            "gaussian": gaussian_payload,
+        }
 
     @pytest.fixture
-    def gaussian_config(self):
+    def dice_config(self, dice_payload: dict[str, Any]) -> DiceConfig:
-        return GaussianConfig(
+        return DiceConfig.model_validate(dice_payload)
-            level="test",
-            qmprog="g16",
-            keywords="test",
-        )
 
-    def test_class_instantiation(self, dice_config, gaussian_config):
+    @pytest.fixture
+    def gaussian_config(self, gaussian_payload: dict[str, Any]):
+        return GaussianConfig.model_validate(gaussian_payload)
+
+    def test_class_instantiation(
+        self, dice_config: DiceConfig, gaussian_config: GaussianConfig
+    ):
         player_dto = PlayerConfig(
-            opt=True,
+            type=RoutineType.BOTH,
             mem=12,
-            maxcyc=100,
+            max_cyc=100,
+            switch_cyc=50,
             nprocs=4,
             ncores=4,
             dice=dice_config,
@@ -61,22 +67,25 @@ class TestPlayerConfig:
         assert isinstance(player_dto.dice, DiceConfig)
         assert isinstance(player_dto.gaussian, GaussianConfig)
 
-    def test_min_altsteps(self, dice_config, gaussian_config):
+    def test_min_altsteps(
+        self, dice_config: DiceConfig, gaussian_config: GaussianConfig
+    ):
         player_dto = PlayerConfig(
-            opt=True,
+            type=RoutineType.BOTH,
             mem=12,
-            maxcyc=100,
+            max_cyc=100,
+            switch_cyc=50,
             nprocs=4,
             ncores=4,
-            altsteps=100,
+            altsteps=0,
             dice=dice_config,
             gaussian=gaussian_config,
         )
 
         assert player_dto.altsteps == 20000
 
-    def test_from_dict(self):
+    def test_from_dict(self, player_payload: dict[str, Any]):
-        player_dto = PlayerConfig.model_validate(get_config_dict())
+        player_dto = PlayerConfig.model_validate(player_payload)
 
         assert isinstance(player_dto, PlayerConfig)
         assert isinstance(player_dto.dice, DiceConfig)

tests/dice/test_dice_input.py

@@ -0,0 +1,65 @@
+from pathlib import Path
+
+from diceplayer.config import PlayerConfig
+from diceplayer.dice.dice_input import NVTEqConfig, NVTTerConfig, NPTTerConfig, NPTEqConfig
+
+import pytest
+
+
+class TestDiceInput:
+    @pytest.fixture
+    def player_config(self) -> PlayerConfig:
+        return PlayerConfig.model_validate(
+            {
+                "type": "both",
+                "mem": 12,
+                "max_cyc": 100,
+                "switch_cyc": 50,
+                "nprocs": 4,
+                "ncores": 4,
+                "dice": {
+                    "ljname": "test",
+                    "outname": "test",
+                    "dens": 1.0,
+                    "nmol": [1],
+                    "nstep": [1, 1],
+                },
+                "gaussian": {
+                    "level": "test",
+                    "qmprog": "g16",
+                    "keywords": "test",
+                },
+            }
+        )
+
+    def test_generate_nvt_ter_input(self, player_config: PlayerConfig):
+        dice_input = NVTTerConfig.from_config(player_config)
+
+        assert isinstance(dice_input, NVTTerConfig)
+
+    def test_generate_nvt_eq_input(self, player_config: PlayerConfig):
+        dice_input = NVTEqConfig.from_config(player_config)
+
+        assert isinstance(dice_input, NVTEqConfig)
+
+    def test_generate_npt_ter_input(self, player_config: PlayerConfig):
+        dice_input = NPTTerConfig.from_config(player_config)
+
+        assert isinstance(dice_input, NPTTerConfig)
+
+    def test_generate_npt_eq_input(self, player_config: PlayerConfig):
+        dice_input = NPTEqConfig.from_config(player_config)
+
+        assert isinstance(dice_input, NPTEqConfig)
+
+    def test_write_dice_config(self, player_config: PlayerConfig, tmp_path: Path):
+        dice_input = NVTTerConfig.from_config(player_config)
+
+        output_file = tmp_path / "nvt_ter.inp"
+
+        with open(output_file, "w") as file:
+            dice_input.write_dice_config(file)
+
+            assert output_file.exists()
+
+        print(output_file.read_text())

tests/state/test_state_handler.py

@@ -12,9 +12,10 @@ class TestStateHandler:
     @pytest.fixture
     def player_config(self) -> PlayerConfig:
         return PlayerConfig(
-            opt=True,
+            type="both",
             mem=12,
-            maxcyc=100,
+            max_cyc=100,
+            switch_cyc=50,
             nprocs=4,
             ncores=4,
             dice=DiceConfig(
@@ -87,7 +88,7 @@ class TestStateHandler:
 
         state_handler.save(state)
 
-        different_config = player_config.model_copy(update={"opt": False})
+        different_config = player_config.model_copy(update={"max_cyc": 200})
 
         retrieved_state = state_handler.get(different_config)
 
@@ -106,7 +107,7 @@ class TestStateHandler:
 
         state_handler.save(state)
 
-        different_config = player_config.model_copy(update={"opt": False})
+        different_config = player_config.model_copy(update={"max_cyc": 200})
 
         retrieved_state = state_handler.get(different_config, force=True)