<jupyter_start><jupyter_text>Kinetic energy distribution*Note: This notebook can be run locally by cloning the*[Github repository](https://github.com/shirtsgroup/physical_validation).*The notebook is located in* `doc/examples/kinetic_energy_distribution.ipynb`. *Be awarethat probabilistic quantities such as error estimates based on bootstrappingwill differ when repeating the analysis.*<jupyter_code># enable plotting in notebook
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>The results imported here are results from example simulations which arestored in another Python file. In real-world usage, the results wouldeither come from the Python interface of the simulation package, fromflat files containing the results, or from package-specific parsers. See[SimulationData](../simulation_data.rst)for more details.<jupyter_code>from simulation_results import example_simulations

import physical_validation<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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"
)<jupyter_output><empty_output><jupyter_text>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 thecomplete system (i.e. no per-molecule or per-atom information) so we will onlyfill these.<jupyter_code>num_molecules = 900
system_data = physical_validation.data.SystemData(
    # Each water molecule has three atoms
    natoms=num_molecules * 3,
    # Each molecule has three constraints
    nconstraints=num_molecules * 3,
    # In this simulation, translational center of mass motion was removed
    ndof_reduction_tra=3,
    # Rotational center of mass motion was not removed
    ndof_reduction_rot=0,
)<jupyter_output><empty_output><jupyter_text>Next, we will create the ensemble information:<jupyter_code>ensemble_data = physical_validation.data.EnsembleData(
    ensemble="NVT",
    natoms=num_molecules * 3,
    volume=3.01125 ** 3,
    temperature=298.15,
)<jupyter_output><empty_output><jupyter_text>Now we will create the simulation data objects which we will feed to the physical validation tests.<jupyter_code>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"]
    ),
)<jupyter_output><empty_output><jupyter_text>We can now check the velocity-rescale temperature coupling, using firstthe 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).<jupyter_code>physical_validation.kinetic_energy.distribution(
    data=simulation_data_vrescale, strict=True, screen=True
)<jupyter_output>After equilibration, decorrelation and tail pruning, 96.02% (4802 frames) of original Kinetic energy remain.<jupyter_text>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:<jupyter_code># 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
)<jupyter_output>After equilibration, decorrelation and tail pruning, 96.02% (4802 frames) of original Kinetic energy remain.
Kinetic energy distribution check (non-strict)
Analytical distribution (T=298.15 K):
 * mu: 6689.47 kJ/mol
 * sigma: 128.77 kJ/mol
Trajectory:
 * mu: 6688.12 +- 1.93 kJ/mol
   T(mu) = 298.09 +- 0.09 K
 * sigma: 128.70 +- 1.35 kJ/mol
   T(sigma) = 297.98 +- 3.12 K<jupyter_text>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. We can now compare this to the Berendsen temperature coupling, using firstthe strict, then the non-strict kinetic energy distribution check.<jupyter_code>physical_validation.kinetic_energy.distribution(
    data=simulation_data_berendsen, strict=True, screen=True
)<jupyter_output>After equilibration, decorrelation and tail pruning, 95.96% (4799 frames) of original Kinetic energy remain.<jupyter_text>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.<jupyter_code># 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
)<jupyter_output>After equilibration, decorrelation and tail pruning, 95.96% (4799 frames) of original Kinetic energy remain.
Kinetic energy distribution check (non-strict)
Analytical distribution (T=298.15 K):
 * mu: 6689.47 kJ/mol
 * sigma: 128.77 kJ/mol
Trajectory:
 * mu: 6690.21 +- 1.49 kJ/mol
   T(mu) = 298.18 +- 0.07 K
 * sigma: 98.81 +- 1.04 kJ/mol
   T(sigma) = 228.78 +- 2.41 K
{'fraction_non_alphanumeric': 0.06688120889883867, 'fraction_numerical': 0.02672449979012173, 'mean_word_length': 5.492279745685741, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '602505', 'license': 'MIT', 'n_tokens_mistral': 2193, 'n_tokens_neox': 1932, 'n_words': 740, 'path': 'doc/examples/kinetic_energy_distribution.ipynb', 'repo_name': 'shirtsgroup/physical_validation', 'size': '7147'}