The equipartition function fails to plot figures for rotational D.O.F in an edge case

[wehs7661]

I was using the package to test the physical validity of a system composed of four vdW sites. The simulation was a 1D metadynamics only biasing the torsion 1-2-3-4. One special thing about the system is that I fixed the lambda state at a fully uncoupled state (state 19 as can be seen in the attached mdp file). The code (check_validity.py) I was using is as follows:

import physical_validation
import numpy as np

if __name__ == "__main__":
    parser = physical_validation.data.GromacsParser()

    res = parser.get_simulation_data(
    mdp='mdout.mdp',
    top='sys2_cccc_lig.top',
    edr='sys2_cccc.edr',
    trr='sys2_cccc_lig.trr'
    )

    res.system.ndof_reduction_tra = 0

    results = physical_validation.kinetic_energy.equipartition(
        data=res,
        strict=True,
        filename='test',
    )

    print(results)

As a result, I got the following error upon the execution of the command python check_validity.py, which seemed that the package failed to plot the histogram for the rotational D.O.F. (The output files corresponding to the total, tra, and rni D.O.F. were still generated before hitting the error.)

Traceback (most recent call last):
  File "check_validity.py", line 19, in <module>
    )
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/physical_validation/kinetic_energy.py", line 250, in equipartition
    temp_unit=data.units.temperature_str
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/physical_validation/util/kinetic_energy.py", line 494, in check_equipartition
    ene_unit=ene_unit, temp_unit=temp_unit))
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/physical_validation/util/kinetic_energy.py", line 975, in test_group
    ene_unit=ene_unit)
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/physical_validation/util/kinetic_energy.py", line 172, in check_distribution
    screen=screen)
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/physical_validation/util/plot.py", line 72, in plot
    _, x, _ = ax.hist(y, r['hist'], **args)
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py", line 1565, in inner
    return func(ax, *map(sanitize_sequence, args), **kwargs)
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/matplotlib/axes/_axes.py", line 6649, in hist
    m, bins = np.histogram(x[i], bins, weights=w[i], **hist_kwargs)
  File "<__array_function__ internals>", line 6, in histogram
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/numpy/lib/histograms.py", line 792, in histogram
    bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
  File "/home/wei-tse/anaconda3/lib/python3.7/site-packages/numpy/lib/histograms.py", line 424, in _get_bin_edges
    raise ValueError('`bins` must be positive, when an integer')
ValueError: `bins` must be positive, when an integer

Here I share the input files for the code: https://drive.google.com/file/d/1SFxH1pL_3TjC3RyCn5EZcZcb8rYpaRKp/view?usp=sharing Note that in the .top and .trr file, the water molecules have been stripped off.

The problem does not happen if the system is fixed at the fully coupled state. Apparently, this is an edge case that fails the equipartition function. I'll look into the problem in the future when I have time, but it would be helpful if someone knows what was probably happening.

[ptmerz]

Hi @wehs7661 , thanks for reporting this. I think this was fixed recently in master, but not yet released. I can't replicate it with the current master version. Could you install the latest development version (you can use pip install git+https://github.com/shirtsgroup/physical-validation.git as described here https://physical-validation.readthedocs.io/en/latest/userguide.html#development-version), and confirm that the problem doesn't arise there?

(Note that the rotational kinetic energy seems to be very off in this case, so the plot that I get from the latest version is not very useful either, but it doesn't crash. I'm not sure off-hand how we could improve that plot, but that's a separate question anyway.)

[wehs7661]

Okay yes, the current version in master solves the problem! I'll close the issue since it has been solved. Thanks!