<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Ensemble check
Note
Step1:  The results imported here are the time series of kinetic and potential energy from example simulations, which are
stored in another Python file. In real-world usage, the results would
either come from the Python interface of the simulation package, from
flat files containing the results, or from package-specific parsers. See
SimulationData
for more details.
Step2:  Check NVT simulations
To check the configurational quantities in NVT, two (otherwise identical)
simulations run at different temperatures are required.
We start by loading the first NVT simulation of 900 water molecules, which
was performed at 298.15K using velocity-rescale temperature coupling.
Step3:  It is not trivial to decide at which temperature to perform a second simulation.
The best results are achieved when the two simulations are close enough to have
good overlap between the distributions, while keeping them far enough apart to
be able to distinguish the physical difference between the distributions from the
statistical error present in simulations.
physical_validation offers functionality to compute a rule-of-thumb estimate of
the optimal interval in state point between two functions.
We will now use our first simulation result to get an estimate of where a second
simulation would optimally be located
Step4:  The second simulation available in our example set was performed at 308.15K, which
is reasonably close to the estimate calculated above. Let's load these results
Step5:  Using both simulation data objects, we can now check the ensemble sampled
by our simulations.
We are using screen=True to display a result plot on screen.
See argument filename to print that same plot to file.
Step6:  By default, the ensemble check is estimating the distance in temperature between the two sampled ensembles using a maximum likelihood approach. This distance estimate is expected to be close to the true value. As a rule of thumb, if the true interval is not within about 2-3 standard deviations of the estimated interval, the trajectory is unlikely to have been sampled from the expected ensemble. The quantiles (number of standard deviations) of difference between the true value and the estimate is returned from the test as a machine-readable test result.
Note that in order to print the plot, the line is also linearly fitted to the simulations. This leads to a slightly different estimate, which explains the difference between the quantiles printed in the plot and in the terminal. As the maximum likelihood estimate is considered to be more exact, it's value is reported on the terminal and used as a return value.
We will now repeat this analysis using the same system, but a simulation which was performed using Berendsen pressure coupling. This temperature coupling method was found not to sample the expected ensemble. We will use this example to illustrate that the physical_validation checks are able to pick up this discrepancy in sampling.
Since the simulated system was identical to the one first analyzed, we will simply replace the observable trajectory in our simulation data objects
Step7:  The check is confirming that the sampled ensemble using the Berendsen thermostat is not behaving as expected when changing the temperature. The reported estimated temperature interval is around 15 standard deviations from the true value, which makes it easy to reject the hypothesis that the potential energy was sampled from the correct ensemble.
Check NPT simulations
To check the sampled ensemble of the configurational quantities in NPT, we again need two otherwise identical simulations which were performed at slightly different state points (target temperature and / or pressure). The checks can be performed using identical pressure and different temperatures, which will test whether the sampled ensembles exhibit the expected temperature dependence. The checks can also be performed using identical temperature and different pressures, which will in turn test the pressure dependence of the sampled ensemble. Finally, we can also use two simulations which differs in both the target temperature and pressure, combining the two tests into one. Here, we will showcase the last option for a system of 900 water molecules, sampled using velocity-rescale temperature coupling and Parrinello-Rahman pressure coupling. These coupling algorithms were analytically shown to sample the correct distribution, so we will check whether the simulated results fulfill this expectation.
We will start by loading a first simulation, performed at 298.15 K and 1 bar.
Step8:  As in the NVT case, we can use this simulation to have physical_validation suggesting a state point to perform the second simulation in.
Step9:  The rule of thumb suggests that a second state point with a temperature difference of about 7.8 K and a pressure difference of about 315 bar would be optimal. The second simulation which is available in our example set was performed at 308.15 K and 101 bar, so at a distance of 10 K and 100 bar. According to the physical_validation estimate, the pressure distance should be a bit further to have optimal error recognition. The check will, however, not be invalid with this choice of state points.
Step10:  Using both simulation data objects, we can now check the ensemble sampled by our simulations. Note that plotting is not available for NPT simulations which differ in both temperature and pressure, since the 2-dimensional plot would be very hard to interpret.
Step11:  The ensemble check now prints both the estimated temperature and pressure intervals. We note that in both cases, the true value is within less than a standard deviation, which means that the null hypothesis of sampling the expected ensemble stands.
It's worth noting that the true pressure difference is given as 98.3 bar rather than 100 bar. When checking simulations which differ in both their pressure and temperature, the pressure interval can only be approximated since the temperature and pressure are not perfectly separable in the NPT partition function. Please refer to Merz & Shirts 2018, eq (18), for details.
Check $\mu$VT simulations
To check the sampled ensemble of the configurational quantities in $\mu$VT, we again need two otherwise identical simulations which were performed at slightly different state points (target temperature and / or chemical potential). The checks can be performed using identical chemical potential and different temperatures, which will test whether the sampled ensembles exhibit the expected temperature dependence. The checks can also be performed using identical temperature and different chemical potentials, which will in turn test the dependence of the sampled ensemble on the chemical potential. Finally, we can also use two simulations which differ in both the target temperature and chemical potentials, combining the two tests into one.
Here, we will showcase checking the dependence on the chemical potential for a system of difluoromethane vapor, sampled using GCMC.
We will start by loading a first simulation, performed at a temperature of 300K and chemical potential of -37 kJ/mol.
Step12:  As for the other ensembles, we can use this simulation to have physical_validation suggesting a state point to perform the second simulation in.
Step13:  The rule of thumb suggests that a second state point with chemical potential difference of about 0.6 kJ/mol would be optimal. For the temperature, the suggestion of a second state point varies depending on whether a 1D fit (simulations varying only in their target temperature) or a 2D fit (simulations varying in both their target temperature and the chemical potential) is to be performed.
The second simulation available in the example set was performed at 300K and -37.0 kJ/mol. The chemical potential difference (0.5 kJ/mol) is close to the value suggested by physical_validation, so we should get close to optimal error recognition when testing for the dependence on the chemical potential.
We proceed by loading the second simulation
Step14:  Using both simulation data objects, we can now check the ensemble sampled by our simulations. We are using screen=True to display a result plot on screen. See argument filename to print that same plot to file.
<ASSISTANT_TASK:>
Python Code:
# enable plotting in notebook
%matplotlib notebook
Explanation: Ensemble check
Note: This notebook can be run locally by cloning the
Github repository.
The notebook is located in doc/examples/ensemble_check.ipynb. Be aware
that probabilistic quantities such as error estimates based on bootstrapping
will differ when repeating the analysis.
End of explanation
from simulation_results import example_simulations
import physical_validation
Explanation: The results imported here are the time series of kinetic and potential energy from example simulations, which are
stored in another Python file. In real-world usage, the results would
either come from the Python interface of the simulation package, from
flat files containing the results, or from package-specific parsers. See
SimulationData
for more details.
End of explanation
simulation_nvt_vrescale_low = example_simulations.get(
    "900 water molecules, NVT at 298K with v-rescale thermostat"
)
num_molecules = 900
simulation_data_nvt_low = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="NVT",
        natoms=num_molecules * 3,
        volume=3.01125 ** 3,
        temperature=298.15,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires only the potential energy
        potential_energy=simulation_nvt_vrescale_low["potential energy"]
    ),
)
Explanation: Check NVT simulations
To check the configurational quantities in NVT, two (otherwise identical)
simulations run at different temperatures are required.
We start by loading the first NVT simulation of 900 water molecules, which
was performed at 298.15K using velocity-rescale temperature coupling.
End of explanation
physical_validation.ensemble.estimate_interval(
    data=simulation_data_nvt_low,
)
Explanation: It is not trivial to decide at which temperature to perform a second simulation.
The best results are achieved when the two simulations are close enough to have
good overlap between the distributions, while keeping them far enough apart to
be able to distinguish the physical difference between the distributions from the
statistical error present in simulations.
physical_validation offers functionality to compute a rule-of-thumb estimate of
the optimal interval in state point between two functions.
We will now use our first simulation result to get an estimate of where a second
simulation would optimally be located:
End of explanation
simulation_nvt_vrescale_high = example_simulations.get(
    "900 water molecules, NVT at 308K with v-rescale thermostat"
)
simulation_data_nvt_high = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="NVT",
        natoms=num_molecules * 3,
        volume=3.01125 ** 3,
        temperature=308.15,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires only the potential energy
        potential_energy=simulation_nvt_vrescale_high["potential energy"]
    ),
)
Explanation: The second simulation available in our example set was performed at 308.15K, which
is reasonably close to the estimate calculated above. Let's load these results:
End of explanation
physical_validation.ensemble.check(
    data_sim_one=simulation_data_nvt_low,
    data_sim_two=simulation_data_nvt_high,
    screen=True,
)
Explanation: Using both simulation data objects, we can now check the ensemble sampled
by our simulations.
We are using screen=True to display a result plot on screen.
See argument filename to print that same plot to file.
End of explanation
simulation_nvt_berendsen_low = example_simulations.get(
    "900 water molecules, NVT at 298K with Berendsen thermostat"
)
simulation_data_nvt_low.observables = physical_validation.data.ObservableData(
    potential_energy=simulation_nvt_berendsen_low["potential energy"]
)
simulation_nvt_berendsen_high = example_simulations.get(
    "900 water molecules, NVT at 308K with Berendsen thermostat"
)
simulation_data_nvt_high.observables = physical_validation.data.ObservableData(
    potential_energy=simulation_nvt_berendsen_high["potential energy"]
)
physical_validation.ensemble.check(
    data_sim_one=simulation_data_nvt_low,
    data_sim_two=simulation_data_nvt_high,
    screen=True,
)
Explanation: By default, the ensemble check is estimating the distance in temperature between the two sampled ensembles using a maximum likelihood approach. This distance estimate is expected to be close to the true value. As a rule of thumb, if the true interval is not within about 2-3 standard deviations of the estimated interval, the trajectory is unlikely to have been sampled from the expected ensemble. The quantiles (number of standard deviations) of difference between the true value and the estimate is returned from the test as a machine-readable test result.
Note that in order to print the plot, the line is also linearly fitted to the simulations. This leads to a slightly different estimate, which explains the difference between the quantiles printed in the plot and in the terminal. As the maximum likelihood estimate is considered to be more exact, it's value is reported on the terminal and used as a return value.
We will now repeat this analysis using the same system, but a simulation which was performed using Berendsen pressure coupling. This temperature coupling method was found not to sample the expected ensemble. We will use this example to illustrate that the physical_validation checks are able to pick up this discrepancy in sampling.
Since the simulated system was identical to the one first analyzed, we will simply replace the observable trajectory in our simulation data objects:
End of explanation
simulation_npt_low = example_simulations.get(
    "900 water molecules, NPT at 298K and 1bar, using v-rescale and Parrinello-Rahman"
)
num_molecules = 900
simulation_data_npt_low = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="NPT",
        natoms=num_molecules * 3,
        pressure=1.0,
        temperature=298.15,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires the potential energy and the volume
        potential_energy=simulation_npt_low["potential energy"],
        volume=simulation_npt_low["volume"],
    ),
)
Explanation: The check is confirming that the sampled ensemble using the Berendsen thermostat is not behaving as expected when changing the temperature. The reported estimated temperature interval is around 15 standard deviations from the true value, which makes it easy to reject the hypothesis that the potential energy was sampled from the correct ensemble.
Check NPT simulations
To check the sampled ensemble of the configurational quantities in NPT, we again need two otherwise identical simulations which were performed at slightly different state points (target temperature and / or pressure). The checks can be performed using identical pressure and different temperatures, which will test whether the sampled ensembles exhibit the expected temperature dependence. The checks can also be performed using identical temperature and different pressures, which will in turn test the pressure dependence of the sampled ensemble. Finally, we can also use two simulations which differs in both the target temperature and pressure, combining the two tests into one. Here, we will showcase the last option for a system of 900 water molecules, sampled using velocity-rescale temperature coupling and Parrinello-Rahman pressure coupling. These coupling algorithms were analytically shown to sample the correct distribution, so we will check whether the simulated results fulfill this expectation.
We will start by loading a first simulation, performed at 298.15 K and 1 bar.
End of explanation
physical_validation.ensemble.estimate_interval(data=simulation_data_npt_low)
Explanation: As in the NVT case, we can use this simulation to have physical_validation suggesting a state point to perform the second simulation in.
End of explanation
simulation_npt_high = example_simulations.get(
    "900 water molecules, NPT at 308K and 101bar, using v-rescale and Parrinello-Rahman"
)
num_molecules = 900
simulation_data_npt_high = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="NPT",
        natoms=num_molecules * 3,
        pressure=101.0,
        temperature=308.15,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires the potential energy and the volume
        potential_energy=simulation_npt_high["potential energy"],
        volume=simulation_npt_high["volume"],
    ),
)
Explanation: The rule of thumb suggests that a second state point with a temperature difference of about 7.8 K and a pressure difference of about 315 bar would be optimal. The second simulation which is available in our example set was performed at 308.15 K and 101 bar, so at a distance of 10 K and 100 bar. According to the physical_validation estimate, the pressure distance should be a bit further to have optimal error recognition. The check will, however, not be invalid with this choice of state points.
End of explanation
physical_validation.ensemble.check(
    data_sim_one=simulation_data_npt_low,
    data_sim_two=simulation_data_npt_high,
)
Explanation: Using both simulation data objects, we can now check the ensemble sampled by our simulations. Note that plotting is not available for NPT simulations which differ in both temperature and pressure, since the 2-dimensional plot would be very hard to interpret.
End of explanation
simulation_muvt_low = example_simulations.get(
    "GCMC Difluoromethane vapor, muVT at 300K and -37.5kJ/mol"
)
simulation_data_muvt_low = physical_validation.data.SimulationData(
    # Example simulations were performed in units compatible with GROMACS units
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="muVT",
        mu=-37.5,
        volume=512,
        temperature=300,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires the potential energy and the volume
        potential_energy=simulation_muvt_low["potential energy"],
        number_of_species=simulation_muvt_low["number of species"],
    ),
)
Explanation: The ensemble check now prints both the estimated temperature and pressure intervals. We note that in both cases, the true value is within less than a standard deviation, which means that the null hypothesis of sampling the expected ensemble stands.
It's worth noting that the true pressure difference is given as 98.3 bar rather than 100 bar. When checking simulations which differ in both their pressure and temperature, the pressure interval can only be approximated since the temperature and pressure are not perfectly separable in the NPT partition function. Please refer to Merz & Shirts 2018, eq (18), for details.
Check $\mu$VT simulations
To check the sampled ensemble of the configurational quantities in $\mu$VT, we again need two otherwise identical simulations which were performed at slightly different state points (target temperature and / or chemical potential). The checks can be performed using identical chemical potential and different temperatures, which will test whether the sampled ensembles exhibit the expected temperature dependence. The checks can also be performed using identical temperature and different chemical potentials, which will in turn test the dependence of the sampled ensemble on the chemical potential. Finally, we can also use two simulations which differ in both the target temperature and chemical potentials, combining the two tests into one.
Here, we will showcase checking the dependence on the chemical potential for a system of difluoromethane vapor, sampled using GCMC.
We will start by loading a first simulation, performed at a temperature of 300K and chemical potential of -37 kJ/mol.
End of explanation
physical_validation.ensemble.estimate_interval(data=simulation_data_muvt_low)
Explanation: As for the other ensembles, we can use this simulation to have physical_validation suggesting a state point to perform the second simulation in.
End of explanation
simulation_muvt_high = example_simulations.get(
    "GCMC Difluoromethane vapor, muVT at 300K and -37.0kJ/mol"
)
simulation_data_muvt_high = physical_validation.data.SimulationData(
    # Example simulations were performed in units compatible with GROMACS units
    units=physical_validation.data.UnitData.units("GROMACS"),
    ensemble=physical_validation.data.EnsembleData(
        ensemble="muVT",
        mu=-37.0,
        volume=512,
        temperature=300,
    ),
    observables=physical_validation.data.ObservableData(
        # This test requires the potential energy and the volume
        potential_energy=simulation_muvt_high["potential energy"],
        number_of_species=simulation_muvt_high["number of species"],
    ),
)
Explanation: The rule of thumb suggests that a second state point with chemical potential difference of about 0.6 kJ/mol would be optimal. For the temperature, the suggestion of a second state point varies depending on whether a 1D fit (simulations varying only in their target temperature) or a 2D fit (simulations varying in both their target temperature and the chemical potential) is to be performed.
The second simulation available in the example set was performed at 300K and -37.0 kJ/mol. The chemical potential difference (0.5 kJ/mol) is close to the value suggested by physical_validation, so we should get close to optimal error recognition when testing for the dependence on the chemical potential.
We proceed by loading the second simulation:
End of explanation
physical_validation.ensemble.check(
    data_sim_one=simulation_data_muvt_low,
    data_sim_two=simulation_data_muvt_high,
    screen=True
)
Explanation: Using both simulation data objects, we can now check the ensemble sampled by our simulations. We are using screen=True to display a result plot on screen. See argument filename to print that same plot to file.
End of explanation
<END_TASK>