Added in ideal gas potential and ideal gas test system for the MC bar…

…ostat.

chiron/potential.py

@@ -63,6 +63,72 @@ def compute_pairlist(self, positions, cutoff) -> jnp.array:
         return distance[interacting_mask], displacement_vectors[interacting_mask], pairs
 
 
+class IdealGasPotential(NeuralNetworkPotential):
+    def __init__(
+        self,
+        topology: Topology,
+    ):
+        """
+        Initialize the Ideal Gas potential.
+
+        Parameters
+        ----------
+        topology : Topology
+            The topology of the system
+
+        """
+
+        if not isinstance(topology, Topology):
+            if not isinstance(topology, property):
+                if topology is not None:
+                    raise TypeError(
+                        f"Topology must be a Topology object or None, type(topology) = {type(topology)}"
+                    )
+
+        self.topology = topology
+
+    def compute_energy(self, positions: jnp.array, nbr_list=None, debug_mode=False):
+        """
+        Compute the energy for an ideal gas, which is always 0.
+
+        Parameters
+        ----------
+        positions : jnp.array
+            The positions of the particles in the system
+        nbr_list : NeighborList, default=None
+            Instance of a neighbor list or pair list class to use.
+            If None, an unoptimized N^2 pairlist will be used without PBC conditions.
+        Returns
+        -------
+        potential_energy : float
+            The total potential energy of the system.
+
+        """
+        # Compute the pair distances and displacement vectors
+
+        return 0.0
+
+    def compute_force(self, positions: jnp.array, nbr_list=None) -> jnp.array:
+        """
+        Compute the  force for ideal gas particles, which is always 0.
+
+        Parameters
+        ----------
+        positions : jnp.array
+            The positions of the particles in the system
+        nbr_list : NeighborList, optional
+            Instance of the neighborlist class to use. By default, set to None, which will use an N^2 pairlist
+
+        Returns
+        -------
+        force : jnp.array
+            The forces on the particles in the system
+
+        """
+
+        return 0.0
+
+
 class LJPotential(NeuralNetworkPotential):
     def __init__(
         self,

chiron/tests/test_testsystems.py

@@ -1,3 +1,13 @@
+import pytest
+
+
+@pytest.fixture(scope="session")
+def prep_temp_dir(tmpdir_factory):
+    """Create a temporary directory for the test."""
+    tmpdir = tmpdir_factory.mktemp("test_testsystems")
+    return tmpdir
+
+
 def compute_openmm_reference_energy(testsystem, positions):
     from openmm import unit
     from openmm.app import Simulation
@@ -203,3 +213,115 @@ def test_LJ_fluid():
         assert jnp.isclose(
             e_chiron_energy, e_openmm_energy.value_in_unit_system(unit.md_unit_system)
         ), "Chiron LJ fluid energy does not match openmm"
+
+
+def test_ideal_gas(prep_temp_dir):
+    from openmmtools.testsystems import IdealGas
+    from openmm import unit
+
+    # Use the LennardJonesFluid example from openmmtools to initialize particle positions and topology
+    # For this example, the topology provides the masses for the particles
+    # The default LennardJonesFluid example considers the system to be Argon with 39.9 amu
+    n_particles = 216
+    temperature = 298 * unit.kelvin
+    pressure = 1 * unit.atmosphere
+    mass = unit.Quantity(39.9, unit.gram / unit.mole)
+
+    ideal_gas = IdealGas(
+        nparticles=n_particles, temperature=temperature, pressure=pressure
+    )
+
+    from chiron.potential import IdealGasPotential
+    import jax.numpy as jnp
+
+    #
+    cutoff = 0.0 * unit.nanometer
+    ideal_gas_potential = IdealGasPotential(ideal_gas.topology)
+
+    from chiron.states import SamplerState, ThermodynamicState
+
+    # define the thermodynamic state
+    thermodynamic_state = ThermodynamicState(
+        potential=ideal_gas_potential,
+        temperature=temperature,
+        pressure=pressure,
+    )
+
+    # define the sampler state
+    sampler_state = SamplerState(
+        x0=ideal_gas.positions,
+        box_vectors=ideal_gas.system.getDefaultPeriodicBoxVectors(),
+    )
+
+    from chiron.neighbors import PairList, OrthogonalPeriodicSpace
+
+    # define the pair list for an orthogonal periodic space
+    # since particles are non-interacting, this will not really do much
+    # but will appropriately wrap particles in space
+    nbr_list = PairList(OrthogonalPeriodicSpace(), cutoff=cutoff)
+    nbr_list.build_from_state(sampler_state)
+
+    from chiron.reporters import SimulationReporter
+
+    # initialize a reporter to save the simulation data
+    filename = f"{prep_temp_dir}/test_ideal_gas_vol.h5"
+
+    reporter1 = SimulationReporter(filename, ideal_gas.topology, 100)
+
+    from chiron.mcmc import (
+        MetropolisDisplacementMove,
+        MoveSet,
+        MCMCSampler,
+        MCBarostatMove,
+    )
+
+    mc_disp_move = MetropolisDisplacementMove(
+        seed=1234,
+        displacement_sigma=0.1 * unit.nanometer,
+        nr_of_moves=10,
+    )
+
+    mc_barostat_move = MCBarostatMove(
+        seed=1234,
+        volume_max_scale=0.2,
+        nr_of_moves=100,
+        adjust_box_scaling=True,
+        adjust_frequency=50,
+        simulation_reporter=reporter1,
+    )
+    move_set = MoveSet(
+        [
+            ("MetropolisDisplacementMove", mc_disp_move),
+            ("MCBarostatMove", mc_barostat_move),
+        ]
+    )
+
+    sampler = MCMCSampler(move_set, sampler_state, thermodynamic_state)
+    sampler.run(n_iterations=30, nbr_list=nbr_list)  # how many times to repeat
+
+    import h5py
+
+    with h5py.File(filename, "r") as f:
+        volume = f["volume"][:]
+        steps = f["step"][:]
+
+    # get expectations
+    ideal_volume = ideal_gas.get_volume_expectation(thermodynamic_state)
+    ideal_volume_std = ideal_gas.get_volume_standard_deviation(thermodynamic_state)
+
+    volume_mean = jnp.mean(jnp.array(volume)) * unit.nanometer**3
+    volume_std = jnp.std(jnp.array(volume)) * unit.nanometer**3
+
+    ideal_density = mass * n_particles / unit.AVOGADRO_CONSTANT_NA / ideal_volume
+    measured_density = mass * n_particles / unit.AVOGADRO_CONSTANT_NA / volume_mean
+
+    assert jnp.isclose(
+        ideal_density.value_in_unit(unit.kilogram / unit.meter**3),
+        measured_density.value_in_unit(unit.kilogram / unit.meter**3),
+        atol=1e-1,
+    )
+    # see if within 5% of ideal volume
+    assert abs(ideal_volume - volume_mean) / ideal_volume < 0.05
+
+    # see if within 10% of the ideal standard deviation of the volume
+    assert abs(ideal_volume_std - volume_std) / ideal_volume_std < 0.1