<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Kinetic energy distribution
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:  In this example, we will check the distributions of two simulations of 900 water molecules simulated under NVT conditions at 298.15K.
One simulation used velocity-rescale temperature coupling, the other used Berendsen temperature coupling.
Step3:  First, we will create the object containing system information.
The example system used in both simulations consists of 900 water molecules which are fully constrained.
For the kinetic energy distribution check, we only need information about the
complete system (i.e. no per-molecule or per-atom information) so we will only
fill these.
Step4:  Next, we will create the ensemble information
Step5:  Now we will create the simulation data objects which we will feed to the physical validation tests.
Step6:  We can now check the velocity-rescale temperature coupling, using first
the strict, then the non-strict kinetic energy distribution check.
We are using screen=True to display a result plot on screen (see argument filename to print that same plot to file).
Step7:  The strict test prints and returns a p-value which indicates that the null hypothesis stands with high confidence. The printed figure confirms that the sampled and the analytical distribution are very similar.
The non-strict test confirms this result
Step8:  The non-strict test calculates the temperature of the mean and the width of the distribution. Analytically, we are expecting both to be close to 298.15K. The result confirms that the simulation behaves as expected, with both the calculated mean and variance being within one standard deviation of the analytical value. The returned tuple is the number of standard deviations the calculated mean and variance differ from the expected value - less than a standard deviations for the mean, and only about 0.05 standard deviations for the variance. Note that any deviation within about 2-3 standard deviations is usually considered as statistically insignificant.
Also note that the error estimates were computed using bootstrapping and are hence a probabilistic quantity which will slightly differ when repeating the analysis.
We can now compare this to the Berendsen temperature coupling, using first
the strict, then the non-strict kinetic energy distribution check.
Step9:  The calculated p-value means that the trajectory is very unklikely to come from a Maxwell-Boltzmann distribution at 298.15K. Commonly, one would consider a p-value of less than 0.05 (5% confidence interval) a reasonable cut-off, which means that the test results make it easy to reject the null hypothesis that the kinetic energy samples the correct distribution. This does not come as a suprise, the Berendsen thermostat is known to keep the correct mean temperature, but not to sample the right distribution.
<ASSISTANT_TASK:>
Python Code:
# enable plotting in notebook
%matplotlib notebook
Explanation: Kinetic energy distribution
Note: This notebook can be run locally by cloning the
Github repository.
The notebook is located in doc/examples/kinetic_energy_distribution.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_vrescale = example_simulations.get(
    "900 water molecules, NVT at 298K with v-rescale thermostat"
)
simulation_berendsen = example_simulations.get(
    "900 water molecules, NVT at 298K with Berendsen thermostat"
)
Explanation: In this example, we will check the distributions of two simulations of 900 water molecules simulated under NVT conditions at 298.15K.
One simulation used velocity-rescale temperature coupling, the other used Berendsen temperature coupling.
End of explanation
num_molecules = 900
# Each water molecule has three atoms
num_atoms_per_molecule = 3
# Each fully constrained water molecule has three constraints
num_constraints_per_molecule = 3
# In this simulation, translational center of mass motion was removed
num_constrained_translational_dof = 3
# Rotational center of mass motion was not removed
num_constrained_rotational_dof = 0
system_data = physical_validation.data.SystemData(
    natoms=num_molecules * num_atoms_per_molecule,
    nconstraints=num_molecules * num_constraints_per_molecule,
    ndof_reduction_tra=num_constrained_translational_dof,
    ndof_reduction_rot=num_constrained_rotational_dof,
)
Explanation: First, we will create the object containing system information.
The example system used in both simulations consists of 900 water molecules which are fully constrained.
For the kinetic energy distribution check, we only need information about the
complete system (i.e. no per-molecule or per-atom information) so we will only
fill these.
End of explanation
ensemble_data = physical_validation.data.EnsembleData(
    ensemble="NVT",
    natoms=num_molecules * 3,
    volume=3.01125 ** 3,
    temperature=298.15,
)
Explanation: Next, we will create the ensemble information:
End of explanation
simulation_data_vrescale = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    system=system_data,
    ensemble=ensemble_data,
    observables=physical_validation.data.ObservableData(
        # This test requires only the kinetic energy
        kinetic_energy=simulation_vrescale["kinetic energy"]
    ),
)
simulation_data_berendsen = physical_validation.data.SimulationData(
    # Example simulations were performed using GROMACS
    units=physical_validation.data.UnitData.units("GROMACS"),
    system=system_data,
    ensemble=ensemble_data,
    observables=physical_validation.data.ObservableData(
        # This test requires only the kinetic energy
        kinetic_energy=simulation_berendsen["kinetic energy"]
    ),
)
Explanation: Now we will create the simulation data objects which we will feed to the physical validation tests.
End of explanation
physical_validation.kinetic_energy.distribution(
    data=simulation_data_vrescale, strict=True, screen=True
)
Explanation: We can now check the velocity-rescale temperature coupling, using first
the strict, then the non-strict kinetic energy distribution check.
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
# We turn plotting off here (`screen=False`), because the plot is 
# identical to the one in the strict test
physical_validation.kinetic_energy.distribution(
    data=simulation_data_vrescale, strict=False, screen=False
)
Explanation: The strict test prints and returns a p-value which indicates that the null hypothesis stands with high confidence. The printed figure confirms that the sampled and the analytical distribution are very similar.
The non-strict test confirms this result:
End of explanation
physical_validation.kinetic_energy.distribution(
    data=simulation_data_berendsen, strict=True, screen=True
)
Explanation: The non-strict test calculates the temperature of the mean and the width of the distribution. Analytically, we are expecting both to be close to 298.15K. The result confirms that the simulation behaves as expected, with both the calculated mean and variance being within one standard deviation of the analytical value. The returned tuple is the number of standard deviations the calculated mean and variance differ from the expected value - less than a standard deviations for the mean, and only about 0.05 standard deviations for the variance. Note that any deviation within about 2-3 standard deviations is usually considered as statistically insignificant.
Also note that the error estimates were computed using bootstrapping and are hence a probabilistic quantity which will slightly differ when repeating the analysis.
We can now compare this to the Berendsen temperature coupling, using first
the strict, then the non-strict kinetic energy distribution check.
End of explanation
# We turn plotting off here (`screen=False`), because the plot is 
# identical to the one in the strict test
physical_validation.kinetic_energy.distribution(
    data=simulation_data_berendsen, strict=False, screen=False
)
Explanation: The calculated p-value means that the trajectory is very unklikely to come from a Maxwell-Boltzmann distribution at 298.15K. Commonly, one would consider a p-value of less than 0.05 (5% confidence interval) a reasonable cut-off, which means that the test results make it easy to reject the null hypothesis that the kinetic energy samples the correct distribution. This does not come as a suprise, the Berendsen thermostat is known to keep the correct mean temperature, but not to sample the right distribution.
End of explanation
<END_TASK>