diff --git a/README.md b/README.md index 6d1feb742e..723d098439 100644 --- a/README.md +++ b/README.md @@ -30,6 +30,8 @@ types of particle simulations. Additional information and publications. - [HOOMD-blue benchmark scripts](https://github.com/glotzerlab/hoomd-benchmarks): Scripts to evaluate the performance of HOOMD-blue simulations. +- [HOOMD-blue validation tests](https://github.com/glotzerlab/hoomd-validation): + Scripts to validate that HOOMD-blue performs accurate simulations. ## Related tools diff --git a/hoomd/conftest.py b/hoomd/conftest.py index 69bcd8fe33..91a35960a3 100644 --- a/hoomd/conftest.py +++ b/hoomd/conftest.py @@ -14,8 +14,6 @@ import atexit import os import numpy -import math -import warnings from hoomd.logging import LoggerCategories from hoomd.snapshot import Snapshot from hoomd import Simulation @@ -548,108 +546,6 @@ def autotuned_kernel_parameter_check(instance, activate, all_optional=False): assert instance.kernel_parameters == initial_kernel_parameters -class BlockAverage: - """Block average method for estimating standard deviation of the mean. - - Args: - data: List of values - """ - - def __init__(self, data): - # round down to the nearest power of 2 - N = 2**int(math.log(len(data)) / math.log(2)) - if N != len(data): - warnings.warn( - "Ignoring some data. Data array should be a power of 2.") - - block_sizes = [] - block_mean = [] - block_variance = [] - - # take means of blocks and the mean/variance of all blocks, growing - # blocks by factors of 2 - block_size = 1 - while block_size <= N // 8: - num_blocks = N // block_size - block_data = numpy.zeros(num_blocks) - - for i in range(0, num_blocks): - start = i * block_size - end = start + block_size - block_data[i] = numpy.mean(data[start:end]) - - block_mean.append(numpy.mean(block_data)) - block_variance.append(numpy.var(block_data) / (num_blocks - 1)) - - block_sizes.append(block_size) - block_size *= 2 - - self._block_mean = numpy.array(block_mean) - self._block_variance = numpy.array(block_variance) - self._block_sizes = numpy.array(block_sizes) - self.data = numpy.array(data) - - # check for a plateau in the relative error before the last data point - block_relative_error = numpy.sqrt(self._block_variance) / numpy.fabs( - self._block_mean) - relative_error_derivative = (numpy.diff(block_relative_error) - / numpy.diff(self._block_sizes)) - if numpy.all(relative_error_derivative > 0): - warnings.warn("Block averaging failed to plateau, run longer") - - def get_hierarchical_errors(self): - """Get details on the hierarchical errors.""" - return (self._block_sizes, self._block_mean, self._block_variance) - - @property - def standard_deviation(self): - """float: The error estimate on the mean.""" - if numpy.all(self.data == self.data[0]): - return 0 - - return numpy.sqrt(numpy.max(self._block_variance)) - - @property - def mean(self): - """float: The mean.""" - return self._block_mean[-1] - - @property - def relative_error(self): - """float: The relative error.""" - return self.standard_deviation / numpy.fabs(self.mean) - - def assert_close(self, - reference_mean, - reference_deviation, - z=6, - max_relative_error=0.02): - """Assert that the distribution is constent with a given reference. - - Also assert that the relative error of the distribution is small. - Otherwise, test runs with massive fluctuations would likely lead to - passing tests. - - Args: - reference_mean: Known good mean value - reference_deviation: Standard deviation of the known good value - z: Number of standard deviations - max_relative_error: Maximum relative error to allow - """ - sample_mean = self.mean - sample_deviation = self.standard_deviation - - assert sample_deviation / sample_mean <= max_relative_error - - # compare if 0 is within the confidence interval around the difference - # of the means - deviation_diff = ((sample_deviation**2 - + reference_deviation**2)**(1 / 2.)) - mean_diff = math.fabs(sample_mean - reference_mean) - deviation_allowed = z * deviation_diff - assert mean_diff <= deviation_allowed - - class ListWriter(hoomd.custom.Action): """Log a single quantity to a list. diff --git a/hoomd/hpmc/pytest/CMakeLists.txt b/hoomd/hpmc/pytest/CMakeLists.txt index c61eb5f3c5..2b11dcf68e 100644 --- a/hoomd/hpmc/pytest/CMakeLists.txt +++ b/hoomd/hpmc/pytest/CMakeLists.txt @@ -20,7 +20,6 @@ set(files __init__.py test_small_box_2d.py test_small_box_3d.py conftest.py - test_nec.py depletants_sphere.py ) diff --git a/hoomd/hpmc/pytest/test_nec.py b/hoomd/hpmc/pytest/test_nec.py deleted file mode 100644 index d040bda598..0000000000 --- a/hoomd/hpmc/pytest/test_nec.py +++ /dev/null @@ -1,112 +0,0 @@ -# Copyright (c) 2009-2022 The Regents of the University of Michigan. -# Part of HOOMD-blue, released under the BSD 3-Clause License. - -"""Test hoomd.hpmc.nec.integrate.Sphere.""" - -import hoomd -import hoomd.hpmc.tune -from hoomd.conftest import operation_pickling_check - -import pytest -import math - -# from test_sphere_eos.py for hoomd-2.9 -# see for example Guang-Wen Wu and Richard J. Sadus, doi:10.1002.aic10233 -_phi_p_ref = [(0.29054, 0.1), (0.91912, 0.2), (2.2768, 0.3), (5.29102, 0.4), - (8.06553, 0.45), (9.98979, 0.475)] -rel_err_cs = 0.0015 - - -@pytest.mark.parametrize("betap,phi", _phi_p_ref) -@pytest.mark.serial -@pytest.mark.validate -def test_sphere_eos_nec(betap, phi, simulation_factory, - lattice_snapshot_factory): - """Test that NEC runs and computes the pressure correctly.""" - n = 7 - - v_particle = 4 / 3 * math.pi * (0.5)**3 - lattice_length = (v_particle / phi)**(1. / 3) - snap = lattice_snapshot_factory(n=n, a=lattice_length) - - sim = simulation_factory(snap) - sim.seed = 123456 - sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT=1) - - mc = hoomd.hpmc.nec.integrate.Sphere(default_d=0.05, update_fraction=0.05) - mc.shape['A'] = dict(diameter=1) - mc.chain_time = 0.05 - sim.operations.integrator = mc - - triggerTune = hoomd.trigger.Periodic(50, 0) - tune_nec_d = hoomd.hpmc.tune.MoveSize.scale_solver( - triggerTune, - moves=['d'], - target=0.10, - tol=0.001, - max_translation_move=0.15) - sim.operations.tuners.append(tune_nec_d) - - tune_nec_ct = hoomd.hpmc.nec.tune.ChainTime.scale_solver(triggerTune, - target=20.0, - tol=1.0, - gamma=20.0) - sim.operations.tuners.append(tune_nec_ct) - - # equilibrate - sim.run(1000) - sim.operations.tuners.clear() - - pressures = [] - mc.nselect = 50 - for i in range(100): - sim.run(1) - pressures.append(mc.virial_pressure) - - mean = sum(pressures) / len(pressures) - variance = sum(p**2 for p in pressures) / len(pressures) - mean**2 - std_dev = variance**0.5 - error = std_dev / (len(pressures) - 1)**0.5 - - # confidence interval, 0.95 quantile of the normal distribution - ci = 1.96 - - assert math.isclose(betap, mean, abs_tol=ci * (error + rel_err_cs * betap)) - assert error < 1 - - mcCounts = mc.counters - necCounts = mc.nec_counters - - assert mcCounts.overlap_errors == 0 - assert necCounts.overlap_errors == 0 - - -nec_test_parameters = [ - ( - hoomd.hpmc.nec.integrate.Sphere, - dict(default_d=0.2, chain_time=0.75, update_fraction=0.125, nselect=2), - dict(diameter=1), - ), - ( - hoomd.hpmc.nec.integrate.ConvexPolyhedron, - dict(default_d=0.2, - default_a=0.1, - chain_probability=0.5, - chain_time=0.75, - update_fraction=0.125, - nselect=2), - dict(vertices=[[1, 1, 1], [1, 1, -1], [1, -1, 1], [1, -1, -1], - [-1, 1, 1], [-1, 1, -1], [-1, -1, 1], [-1, -1, -1]]), - ), -] - - -@pytest.mark.serial -@pytest.mark.parametrize('integrator_cls, integrator_args, shape', - nec_test_parameters) -def test_pickling(integrator_cls, integrator_args, shape, simulation_factory, - two_particle_snapshot_factory): - mc = integrator_cls(**integrator_args) - mc.shape["A"] = shape - sim = simulation_factory(two_particle_snapshot_factory(L=1000)) - operation_pickling_check(mc, sim) diff --git a/hoomd/md/pytest/CMakeLists.txt b/hoomd/md/pytest/CMakeLists.txt index 0bec714e23..508fb8add3 100644 --- a/hoomd/md/pytest/CMakeLists.txt +++ b/hoomd/md/pytest/CMakeLists.txt @@ -19,7 +19,6 @@ set(files __init__.py test_dihedral.py test_flags.py test_kernel_parameters.py - test_lj_equation_of_state.py test_potential.py test_pppm_coulomb.py test_manifolds.py diff --git a/hoomd/md/pytest/test_alj.py b/hoomd/md/pytest/test_alj.py index 77c5ee2a9a..4db282bb26 100644 --- a/hoomd/md/pytest/test_alj.py +++ b/hoomd/md/pytest/test_alj.py @@ -2,7 +2,6 @@ # Part of HOOMD-blue, released under the BSD 3-Clause License. import numpy as np -import numpy.testing as npt import pytest import hoomd @@ -10,125 +9,6 @@ from hoomd import md -@pytest.mark.validate -def test_conservation(simulation_factory, lattice_snapshot_factory): - # For test, use a unit area hexagon. - particle_vertices = np.array([[6.20403239e-01, 0.00000000e+00, 0], - [3.10201620e-01, 5.37284966e-01, 0], - [-3.10201620e-01, 5.37284966e-01, 0], - [-6.20403239e-01, 7.59774841e-17, 0], - [-3.10201620e-01, -5.37284966e-01, 0], - [3.10201620e-01, -5.37284966e-01, 0]]) - area = 1.0 - circumcircle_radius = 0.6204032392788702 - incircle_radius = 0.5372849659264116 - num_vertices = len(particle_vertices) - - circumcircle_diameter = 2 * circumcircle_radius - - # Just initialize in a simple cubic lattice. - sim = simulation_factory( - lattice_snapshot_factory(a=4 * circumcircle_diameter, - n=10, - dimensions=2)) - sim.seed = 175 - - # Initialize moments of inertia since original simulation was HPMC. - mass = area - # https://math.stackexchange.com/questions/2004798/moment-of-inertia-for-a-n-sided-regular-polygon # noqa - moment_inertia = (mass * circumcircle_diameter**2 / 6) - moment_inertia *= (1 + 2 * np.cos(np.pi / num_vertices)**2) - - with sim.state.cpu_local_snapshot as snapshot: - snapshot.particles.mass[:] = mass - snapshot.particles.moment_inertia[:] = np.array([0, 0, moment_inertia]) - # Not sure if this should be incircle or circumcircle; - # probably doesn't matter based on current usage, but may - # matter in the future for the potential if it's modified - # to actually use diameter. - snapshot.particles.diameter[:] = circumcircle_diameter - kT = 0.3 - sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT) - - # Create box resize updater - packing_fraction = 0.4 - final_area = area * sim.state.N_particles / packing_fraction - L_final = np.sqrt(final_area) - final_box = hoomd.Box.square(L_final) - - n_compression_start = int(1e4) - n_compression_end = int(1e5) - n_compression_total = n_compression_end - n_compression_start - - box_resize = hoomd.update.BoxResize( - box1=sim.state.box, - box2=final_box, - trigger=int(n_compression_total / 10000), - variant=hoomd.variant.Ramp(0, 1, n_compression_start, - n_compression_total), - filter=hoomd.filter.All()) - sim.operations += box_resize - - # Define forces and methods - r_cut_scale = 1.3 - kernel_scale = (1 / np.cos(np.pi / num_vertices)) - incircle_diameter = 2 * incircle_radius - r_cut_set = incircle_diameter * kernel_scale * r_cut_scale - - alj = md.pair.aniso.ALJ(default_r_cut=r_cut_set, nlist=md.nlist.Cell(0.4)) - - alj.shape["A"] = { - "vertices": particle_vertices, - "faces": [], - "rounding_radii": 0 - } - - alpha = 0 # Make it WCA-only (no attraction) - eps_att = 1.0 - alj.params[("A", "A")] = { - "epsilon": eps_att, - "sigma_i": incircle_diameter, - "sigma_j": incircle_diameter, - "alpha": alpha - } - - nve = md.methods.NVE(filter=hoomd.filter.All()) - integrator = md.Integrator(dt=1e-4, - forces=[alj], - methods=[nve], - integrate_rotational_dof=True) - sim.operations.integrator = integrator - - # Compress box - sim.run(n_compression_end) - - thermo = md.compute.ThermodynamicQuantities(hoomd.filter.All()) - sim.operations += thermo - - # Reset velocities after the compression, and equilibriate - sim.state.thermalize_particle_momenta(hoomd.filter.All(), kT) - sim.run(1000) - - # run sim and get values back - w = hoomd.conftest.ManyListWriter([(thermo, 'potential_energy'), - (thermo, 'kinetic_energy'), - (integrator, 'linear_momentum')]) - writer = hoomd.write.CustomWriter(action=w, trigger=1) - sim.operations.writers.append(writer) - sim.run(1000) - pe, ke, momentum = w.data - total_energies = np.array(pe) + np.array(ke) - - # Ensure energy conservation up to the 3 digit per-particle. - npt.assert_allclose(total_energies, - total_energies[0], - atol=0.03 * sim.state.N_particles) - - # Test momentum conservation. - p_magnitude = np.linalg.norm(momentum, axis=-1) - npt.assert_allclose(p_magnitude, p_magnitude[0], atol=1e-13) - - def test_type_shapes(simulation_factory, two_particle_snapshot_factory): alj = md.pair.aniso.ALJ(md.nlist.Cell(buffer=0.1)) sim = simulation_factory(two_particle_snapshot_factory(d=2.0)) diff --git a/hoomd/md/pytest/test_lj_equation_of_state.py b/hoomd/md/pytest/test_lj_equation_of_state.py deleted file mode 100644 index 7092f28a68..0000000000 --- a/hoomd/md/pytest/test_lj_equation_of_state.py +++ /dev/null @@ -1,127 +0,0 @@ -# Copyright (c) 2009-2022 The Regents of the University of Michigan. -# Part of HOOMD-blue, released under the BSD 3-Clause License. - -import hoomd -import pytest -import hoomd.conftest -import numpy - -# Selected state point(s) in the high density fluid. These state point(s) are in -# single phase regions of the phase diagram and suitable for NVT/NPT and MC -# cross validation. -# T_star, rho_star, mean_U_ref, sigma_U_ref, mean_P_ref, sigma_P_ref, -# log_period, equilibration_steps, run_steps -statepoints = [ - (1.4, 0.9, -4.6622, 0.0006089, 6.6462, 0.00328, 64, 2**10, 2**16), -] - - -@pytest.mark.validate -@pytest.mark.parametrize( - 'T_star, rho_star, mean_U_ref, sigma_U_ref, mean_P_ref, sigma_P_ref,' - 'log_period, equilibration_steps, run_steps', statepoints) -@pytest.mark.parametrize('method_name', ['Langevin', 'NVT', 'NPT']) -def test_lj_equation_of_state( - T_star, - rho_star, - mean_U_ref, - sigma_U_ref, - mean_P_ref, - sigma_P_ref, - log_period, - equilibration_steps, - run_steps, - method_name, - fcc_snapshot_factory, - simulation_factory, - device, -): - # construct the system at the given density - n = 6 - if device.communicator.num_ranks > 1: - # MPI tests need a box large enough to decompose - n = 8 - N = n**3 * 4 - V = N / rho_star - L = V**(1 / 3) - r_cut = 2.5 - a = L / n - - snap = fcc_snapshot_factory(n=n, a=a) - sim = simulation_factory(snap) - sim.seed = 10 - - # set the simulation parameters - integrator = hoomd.md.Integrator(dt=0.005) - lj = hoomd.md.pair.LJ(nlist=hoomd.md.nlist.Cell(buffer=0.4), - default_r_cut=r_cut, - mode='shift') - lj.params.default = {'sigma': 1, 'epsilon': 1} - integrator.forces.append(lj) - if method_name == 'NVT': - method = hoomd.md.methods.NVT(filter=hoomd.filter.All(), - kT=T_star, - tau=0.1) - elif method_name == 'Langevin': - method = hoomd.md.methods.Langevin(filter=hoomd.filter.All(), kT=T_star) - elif method_name == 'NPT': - method = hoomd.md.methods.NPT(filter=hoomd.filter.All(), - kT=T_star, - tau=0.1, - S=mean_P_ref, - tauS=0.5, - couple='xyz') - integrator.methods.append(method) - sim.operations.integrator = integrator - - # equilibrate the simulation - sim.run(0) - sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=T_star) - if method_name == 'NVT': - method.thermalize_thermostat_dof() - elif method_name == 'NPT': - method.thermalize_thermostat_and_barostat_dof() - sim.run(equilibration_steps) - - # log energy and pressure - thermo = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All()) - sim.operations.computes.append(thermo) - energy_log = hoomd.conftest.ListWriter(thermo, 'potential_energy') - pressure_log = hoomd.conftest.ListWriter(thermo, 'pressure') - volume_log = hoomd.conftest.ListWriter(thermo, 'volume') - log_trigger = hoomd.trigger.Periodic(log_period) - sim.operations.writers.append( - hoomd.write.CustomWriter(action=energy_log, trigger=log_trigger)) - sim.operations.writers.append( - hoomd.write.CustomWriter(action=pressure_log, trigger=log_trigger)) - sim.operations.writers.append( - hoomd.write.CustomWriter(action=volume_log, trigger=log_trigger)) - - sim.always_compute_pressure = True - sim.run(run_steps) - - # compute the average and error - energy = hoomd.conftest.BlockAverage(numpy.array(energy_log.data) / N) - pressure = hoomd.conftest.BlockAverage(pressure_log.data) - rho = hoomd.conftest.BlockAverage(N / numpy.array(volume_log.data)) - - # Useful information to know when the test fails - print('U_ref = ', mean_U_ref, '+/-', sigma_U_ref) - print('U = ', energy.mean, '+/-', energy.standard_deviation, '(', - energy.relative_error * 100, '%)') - print('P_ref = ', mean_P_ref, '+/-', sigma_P_ref) - print('P = ', pressure.mean, '+/-', pressure.standard_deviation, - pressure.standard_deviation, '(', pressure.relative_error * 100, '%)') - print('rho = ', rho.mean, '+/-', rho.standard_deviation) - - print(f'Statepoint entry: {T_star:0.4}, {rho_star:0.4}, ' - f'{energy.mean:0.5}, {energy.standard_deviation:0.4}, ' - f'{pressure.mean:0.5}, {pressure.standard_deviation:0.4}') - - energy.assert_close(mean_U_ref, sigma_U_ref) - - if method_name == 'NVT' or method_name == 'Langevin': - pressure.assert_close(mean_P_ref, sigma_P_ref) - - if method_name == 'NPT': - rho.assert_close(rho_star, 0) diff --git a/sphinx-doc/index.rst b/sphinx-doc/index.rst index d698613315..a035512df2 100644 --- a/sphinx-doc/index.rst +++ b/sphinx-doc/index.rst @@ -51,6 +51,8 @@ Resources Additional information and publications. - `HOOMD-blue benchmark scripts `_: Scripts to evaluate the performance of HOOMD-blue simulations. +- `HOOMD-blue validation tests `_: + Scripts to validate that HOOMD-blue performs accurate simulations. Related tools ============= diff --git a/sphinx-doc/testing.rst b/sphinx-doc/testing.rst index 5c9ebfbd2b..67fe6ae17d 100644 --- a/sphinx-doc/testing.rst +++ b/sphinx-doc/testing.rst @@ -19,12 +19,12 @@ APIs to provide inputs and check for correct outputs. For example, test that the hard sphere HPMC simulation executes for several steps. System integration tests may take several seconds. -**Validation tests** rigorously check that HOOMD simulations sample the correct -statistical ensembles. For example, validate the a Lennard-Jones simulation at a -given density matches the pressure in the NIST reference. Validation tests -should run long enough to ensure reasonable sampling, but not too long. These -test run in a CI environment on every pull request. Individual validation tests -should execute in less than 10 minutes. +**Long running tests** check for correct behavior, but require up to a minute to execute. Mark long +running tests with the ``validate`` label. + +**Validation tests** rigorously check that simulations sample the correct statistical ensembles. +These tests take hours to execute on many CPU cores or GPUs. Find HOOMD's validation tests in the +hoomd-validation_ repository. Requirements ------------ @@ -93,22 +93,22 @@ this file when a test fails on rank 1. This will result in an ``MPI_ABORT`` on r .. _OpenMPI: https://www.open-mpi.org/ -Running validation tests ------------------------- +Executing long runing tests +--------------------------- -Longer running validation tests do not execute by default. Run these with the ``--validate`` command -line option to pytest:: +Longer running tests do not execute by default. Run these with the ``--validate`` command line +option to pytest:: $ python3 -m pytest build/hoomd --validate -m validate $ mpirun -n 2 hoomd/pytest/pytest-openmpi.sh build/hoomd -v -x -ra --validate -m validate .. note:: - The ``-m validate`` option selects *only* the validation tests. + The ``-m validate`` option selects *only* the long running tests. .. note:: - To run validation tests on an installed ``hoomd`` package, you need to specify additional + To run long running tests on an installed ``hoomd`` package, you need to specify additional options:: python3 -m pytest --pyargs hoomd -p hoomd.pytest_plugin_validate -m validate --validate @@ -125,4 +125,7 @@ time. Add any new ``test_*.py`` files to the list in the corresponding ``CMakeLists.txt`` file. Only add C++ tests for classes that have no Python interface or otherwise require low level testing. -If you are unsure, please check with the lead developers prior to adding new C++ tests. +If you are unsure, please check with the lead developers prior to adding new C++ tests. Add +new validation tests to the hoomd-validation_ repository. + +.. _hoomd-validation: https://github.com/glotzerlab/hoomd-validation/