diff --git a/.github/workflows/requirements-test.txt b/.github/workflows/requirements-test.txt index ac142c32..781c7c7f 100644 --- a/.github/workflows/requirements-test.txt +++ b/.github/workflows/requirements-test.txt @@ -1,5 +1,6 @@ -signac==2.1.0 -signac-flow==0.26.1 +h5py==3.10.0 +gsd==3.2.0 numpy==1.26.1 PyYAML==6.0.1 -gsd==3.2.0 +signac==2.1.0 +signac-flow==0.26.1 diff --git a/hoomd_validation/alj_2d.py b/hoomd_validation/alj_2d.py index 97a2be53..484e541e 100644 --- a/hoomd_validation/alj_2d.py +++ b/hoomd_validation/alj_2d.py @@ -349,7 +349,6 @@ def alj_2d_nve_md_job(*jobs): aggregator=analysis_aggregator) def alj_2d_conservation_analyze(*jobs): """Analyze the output of NVE simulations and inspect conservation.""" - import gsd.hoomd import math import matplotlib import matplotlib.style @@ -359,7 +358,7 @@ def alj_2d_conservation_analyze(*jobs): print('starting alj_2d_conservation_analyze:', jobs[0]) sim_modes = ['nve_md_cpu'] - if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.h5')): sim_modes.extend(['nve_md_gpu']) timesteps = [] @@ -372,20 +371,22 @@ def alj_2d_conservation_analyze(*jobs): job_linear_momentum = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - job_timesteps[sim_mode] = log_traj['configuration/step'] + job_timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] job_energies[sim_mode] = ( log_traj[ - 'log/md/compute/ThermodynamicQuantities/potential_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/potential_energy'] + log_traj[ - 'log/md/compute/ThermodynamicQuantities/kinetic_energy']) + 'hoomd-data/md/compute/ThermodynamicQuantities/kinetic_energy'] + ) job_energies[sim_mode] = ( job_energies[sim_mode] - job_energies[sim_mode][0]) / job.statepoint["num_particles"] - momentum_vector = log_traj['log/md/Integrator/linear_momentum'] + momentum_vector = log_traj[ + 'hoomd-data/md/Integrator/linear_momentum'] job_linear_momentum[sim_mode] = [ math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) / job.statepoint["num_particles"] for v in momentum_vector diff --git a/hoomd_validation/hard_disk.py b/hoomd_validation/hard_disk.py index 1a895b26..a14555b6 100644 --- a/hoomd_validation/hard_disk.py +++ b/hoomd_validation/hard_disk.py @@ -164,13 +164,12 @@ def make_mc_simulation(job, compute_density = ComputeDensity() sdf = hoomd.hpmc.compute.SDF(xmax=0.02, dx=1e-4) - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, quantities=['translate_moves']) - logger_gsd.add(sdf, quantities=['betaP']) - logger_gsd.add(compute_density, quantities=['density']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, quantities=['translate_moves']) + logger.add(sdf, quantities=['betaP']) + logger.add(compute_density, quantities=['density']) for loggable, quantity in extra_loggables: - logger_gsd.add(loggable, quantities=[quantity]) + logger.add(loggable, quantities=[quantity]) # make simulation sim = util.make_simulation(job=job, @@ -178,7 +177,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -398,14 +397,12 @@ def run_nec_sim(job, device, complete_filename): # compute the density compute_density = ComputeDensity() - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, - quantities=[ - 'translate_moves', 'particles_per_chain', - 'virial_pressure' - ]) - logger_gsd.add(compute_density, quantities=['density']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, + quantities=[ + 'translate_moves', 'particles_per_chain', 'virial_pressure' + ]) + logger.add(compute_density, quantities=['density']) # make simulation sim = util.make_simulation(job=job, @@ -413,7 +410,7 @@ def run_nec_sim(job, device, complete_filename): initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -588,7 +585,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def hard_disk_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -603,7 +599,7 @@ def hard_disk_analyze(job): 'npt_cpu', ] - if os.path.exists(job.fn('nvt_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_gpu_quantities.h5')): sim_modes.extend(['nvt_gpu']) util._sort_sim_modes(sim_modes) @@ -613,18 +609,18 @@ def hard_disk_analyze(job): densities = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - timesteps[sim_mode] = log_traj['configuration/step'] + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] if 'nec' in sim_mode: pressures[sim_mode] = log_traj[ - 'log/hpmc/nec/integrate/Sphere/virial_pressure'] + 'hoomd-data/hpmc/nec/integrate/Sphere/virial_pressure'] else: - pressures[sim_mode] = log_traj['log/hpmc/compute/SDF/betaP'] + pressures[sim_mode] = log_traj['hoomd-data/hpmc/compute/SDF/betaP'] densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] # save averages for mode in sim_modes: @@ -683,7 +679,7 @@ def hard_disk_compare_modes(*jobs): 'npt_cpu', ] - if os.path.exists(jobs[0].fn('nvt_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_gpu_quantities.h5')): sim_modes.extend(['nvt_gpu']) util._sort_sim_modes(sim_modes) diff --git a/hoomd_validation/hard_sphere.py b/hoomd_validation/hard_sphere.py index 01327f49..72cc0865 100644 --- a/hoomd_validation/hard_sphere.py +++ b/hoomd_validation/hard_sphere.py @@ -159,13 +159,12 @@ def make_mc_simulation(job, compute_density = ComputeDensity() sdf = hoomd.hpmc.compute.SDF(xmax=0.02, dx=1e-4) - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, quantities=['translate_moves']) - logger_gsd.add(sdf, quantities=['betaP']) - logger_gsd.add(compute_density, quantities=['density']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, quantities=['translate_moves']) + logger.add(sdf, quantities=['betaP']) + logger.add(compute_density, quantities=['density']) for loggable, quantity in extra_loggables: - logger_gsd.add(loggable, quantities=[quantity]) + logger.add(loggable, quantities=[quantity]) # make simulation sim = util.make_simulation(job=job, @@ -173,7 +172,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -310,14 +309,12 @@ def run_nec_sim(job, device, complete_filename): # compute the density compute_density = ComputeDensity() - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, - quantities=[ - 'translate_moves', 'particles_per_chain', - 'virial_pressure' - ]) - logger_gsd.add(compute_density, quantities=['density']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, + quantities=[ + 'translate_moves', 'particles_per_chain', 'virial_pressure' + ]) + logger.add(compute_density, quantities=['density']) # make simulation sim = util.make_simulation(job=job, @@ -325,7 +322,7 @@ def run_nec_sim(job, device, complete_filename): initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -467,7 +464,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def hard_sphere_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -482,7 +478,7 @@ def hard_sphere_analyze(job): 'npt_cpu', ] - if os.path.exists(job.fn('nvt_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_gpu_quantities.h5')): sim_modes.extend(['nvt_gpu']) util._sort_sim_modes(sim_modes) @@ -492,17 +488,17 @@ def hard_sphere_analyze(job): densities = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) - timesteps[sim_mode] = log_traj['configuration/step'] + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] if 'nec' in sim_mode: pressures[sim_mode] = log_traj[ - 'log/hpmc/nec/integrate/Sphere/virial_pressure'] + 'hoomd-data/hpmc/nec/integrate/Sphere/virial_pressure'] else: - pressures[sim_mode] = log_traj['log/hpmc/compute/SDF/betaP'] + pressures[sim_mode] = log_traj['hoomd-data/hpmc/compute/SDF/betaP'] densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] # save averages for mode in sim_modes: @@ -561,7 +557,7 @@ def hard_sphere_compare_modes(*jobs): 'npt_cpu', ] - if os.path.exists(jobs[0].fn('nvt_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_gpu_quantities.h5')): sim_modes.extend(['nvt_gpu']) util._sort_sim_modes(sim_modes) diff --git a/hoomd_validation/lj_fluid.py b/hoomd_validation/lj_fluid.py index e51ba2e7..755130ae 100644 --- a/hoomd_validation/lj_fluid.py +++ b/hoomd_validation/lj_fluid.py @@ -524,13 +524,12 @@ def make_mc_simulation(job, # compute the density compute_density = ComputeDensity() - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(lj_jit_potential, quantities=['energy']) - logger_gsd.add(mc, quantities=['translate_moves']) - logger_gsd.add(compute_density) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(lj_jit_potential, quantities=['energy']) + logger.add(mc, quantities=['translate_moves']) + logger.add(compute_density) for loggable in extra_loggables: - logger_gsd.add(loggable) + logger.add(loggable) # make simulation sim = util.make_simulation(job=job, @@ -538,7 +537,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -561,8 +560,7 @@ def _compute_virial_pressure(): return job.statepoint.num_particles * job.statepoint.kT / V + w / (3 * V) - logger_gsd[('custom', 'virial_pressure')] = (_compute_virial_pressure, - 'scalar') + logger[('custom', 'virial_pressure')] = (_compute_virial_pressure, 'scalar') # move size tuner mstuner = hpmc.tune.MoveSize.scale_solver( @@ -832,7 +830,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def lj_fluid_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import math import matplotlib @@ -849,7 +846,7 @@ def lj_fluid_analyze(job): 'npt_bussi_md_cpu', ] - if os.path.exists(job.fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -857,10 +854,10 @@ def lj_fluid_analyze(job): 'npt_bussi_md_gpu', ]) - if os.path.exists(job.fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(job.fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -872,30 +869,31 @@ def lj_fluid_analyze(job): linear_momentum = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - timesteps[sim_mode] = log_traj['configuration/step'] + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] if 'md' in sim_mode: energies[sim_mode] = log_traj[ - 'log/md/compute/ThermodynamicQuantities/potential_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/potential_energy'] else: energies[sim_mode] = log_traj[ - 'log/hpmc/pair/user/CPPPotential/energy'] * job.statepoint.kT + 'hoomd-data/hpmc/pair/user/CPPPotential/energy'] * job.statepoint.kT energies[sim_mode] /= job.statepoint.num_particles if 'md' in sim_mode: pressures[sim_mode] = log_traj[ - 'log/md/compute/ThermodynamicQuantities/pressure'] + 'hoomd-data/md/compute/ThermodynamicQuantities/pressure'] else: - pressures[sim_mode] = log_traj['log/custom/virial_pressure'] + pressures[sim_mode] = log_traj['hoomd-data/custom/virial_pressure'] densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] if 'md' in sim_mode and 'langevin' not in sim_mode: - momentum_vector = log_traj['log/md/Integrator/linear_momentum'] + momentum_vector = log_traj[ + 'hoomd-data/md/Integrator/linear_momentum'] linear_momentum[sim_mode] = [ math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) for v in momentum_vector ] @@ -984,7 +982,7 @@ def lj_fluid_compare_modes(*jobs): 'npt_bussi_md_cpu', ] - if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -992,10 +990,10 @@ def lj_fluid_compare_modes(*jobs): 'npt_bussi_md_gpu', ]) - if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -1085,7 +1083,6 @@ def lj_fluid_compare_modes(*jobs): aggregator=analysis_aggregator) def lj_fluid_distribution_analyze(*jobs): """Checks that MD follows the correct KE distribution.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -1102,7 +1099,7 @@ def lj_fluid_distribution_analyze(*jobs): 'npt_bussi_md_cpu', ] - if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -1110,10 +1107,10 @@ def lj_fluid_distribution_analyze(*jobs): 'npt_bussi_md_gpu', ]) - if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -1143,12 +1140,12 @@ def lj_fluid_distribution_analyze(*jobs): else: n_dof = num_particles * 3 - 3 - print('Reading' + job.fn(sim_mode + '_quantities.gsd')) - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + print('Reading' + job.fn(sim_mode + '_quantities.h5')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) if 'md' in sim_mode: ke = log_traj[ - 'log/md/compute/ThermodynamicQuantities/kinetic_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/kinetic_energy'] ke_means_expected[sim_mode].append( numpy.mean(ke) - 1 / 2 * n_dof * kT) ke_sigmas_expected[sim_mode].append( @@ -1161,23 +1158,28 @@ def lj_fluid_distribution_analyze(*jobs): if 'md' in sim_mode: potential_energy_samples[sim_mode].extend( - list(log_traj['log/md/compute/ThermodynamicQuantities' - '/potential_energy'])) + list( + log_traj['hoomd-data/md/compute/ThermodynamicQuantities' + '/potential_energy'])) else: potential_energy_samples[sim_mode].extend( - list(log_traj['log/hpmc/pair/user/CPPPotential/energy'] - * job.statepoint.kT)) + list( + log_traj['hoomd-data/hpmc/pair/user/CPPPotential/energy'] + * job.statepoint.kT)) if 'md' in sim_mode: pressure_samples[sim_mode].extend( list(log_traj[ - 'log/md/compute/ThermodynamicQuantities/pressure'])) + 'hoomd-data/md/compute/ThermodynamicQuantities/pressure'] + )) else: pressure_samples[sim_mode].extend( - list(log_traj['log/custom/virial_pressure'])) + list(log_traj['hoomd-data/custom/virial_pressure'])) density_samples[sim_mode].extend( - list(log_traj['log/custom_actions/ComputeDensity/density'])) + list( + log_traj['hoomd-data/custom_actions/ComputeDensity/density'] + )) ax = fig.add_subplot(2, 2, 1) util.plot_vs_expected(ax, ke_means_expected, '$ - 1/2 N_{dof} k T$') @@ -1370,7 +1372,6 @@ def lj_fluid_nve_md_job(*jobs): aggregator=nve_analysis_aggregator) def lj_fluid_conservation_analyze(*jobs): """Analyze the output of NVE simulations and inspect conservation.""" - import gsd.hoomd import numpy import math import matplotlib @@ -1381,7 +1382,7 @@ def lj_fluid_conservation_analyze(*jobs): print('starting lj_fluid_conservation_analyze:', jobs[0]) sim_modes = ['nve_md_cpu'] - if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.h5')): sim_modes.extend(['nve_md_gpu']) timesteps = [] @@ -1394,20 +1395,22 @@ def lj_fluid_conservation_analyze(*jobs): job_linear_momentum = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - job_timesteps[sim_mode] = log_traj['configuration/step'] + job_timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] job_energies[sim_mode] = ( log_traj[ - 'log/md/compute/ThermodynamicQuantities/potential_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/potential_energy'] + log_traj[ - 'log/md/compute/ThermodynamicQuantities/kinetic_energy']) + 'hoomd-data/md/compute/ThermodynamicQuantities/kinetic_energy'] + ) job_energies[sim_mode] = ( job_energies[sim_mode] - job_energies[sim_mode][0]) / job.statepoint["num_particles"] - momentum_vector = log_traj['log/md/Integrator/linear_momentum'] + momentum_vector = log_traj[ + 'hoomd-data/md/Integrator/linear_momentum'] job_linear_momentum[sim_mode] = [ math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) / job.statepoint["num_particles"] for v in momentum_vector diff --git a/hoomd_validation/lj_union.py b/hoomd_validation/lj_union.py index 64b0d218..aa2af4e1 100644 --- a/hoomd_validation/lj_union.py +++ b/hoomd_validation/lj_union.py @@ -567,13 +567,12 @@ def make_mc_simulation(job, # compute the density compute_density = ComputeDensity() - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(lj_jit_potential, quantities=['energy']) - logger_gsd.add(mc, quantities=['translate_moves', 'rotate_moves']) - logger_gsd.add(compute_density) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(lj_jit_potential, quantities=['energy']) + logger.add(mc, quantities=['translate_moves', 'rotate_moves']) + logger.add(compute_density) for loggable in extra_loggables: - logger_gsd.add(loggable) + logger.add(loggable) # make simulation sim = util.make_simulation(job=job, @@ -581,7 +580,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -866,7 +865,8 @@ def sampling_operation(*jobs): device = device_cls(communicator=communicator, message_filename=util.get_message_filename( - job, f'{mode}_mc_{device_name}.log')) + job, f'{mode}_mc_{device_name}.log'), + num_cpu_threads=1) globals().get(f'run_{mode}_mc_sim')( job, device, complete_filename=f'{mode}_mc_{device_name}_complete') @@ -890,7 +890,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def lj_union_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import math import matplotlib @@ -907,7 +906,7 @@ def lj_union_analyze(job): 'npt_bussi_md_cpu', ] - if os.path.exists(job.fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -915,10 +914,10 @@ def lj_union_analyze(job): 'npt_bussi_md_gpu', ]) - if os.path.exists(job.fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(job.fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(job.fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -930,31 +929,32 @@ def lj_union_analyze(job): linear_momentum = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - timesteps[sim_mode] = log_traj['configuration/step'] + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] if 'md' in sim_mode: energies[sim_mode] = log_traj[ - 'log/md/compute/ThermodynamicQuantities/potential_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/potential_energy'] else: energies[sim_mode] = ( - log_traj['log/hpmc/pair/user/CPPPotentialUnion/energy'] + log_traj['hoomd-data/hpmc/pair/user/CPPPotentialUnion/energy'] * job.statepoint.kT) energies[sim_mode] /= job.statepoint.num_particles if 'md' in sim_mode: pressures[sim_mode] = log_traj[ - 'log/md/compute/ThermodynamicQuantities/pressure'] + 'hoomd-data/md/compute/ThermodynamicQuantities/pressure'] else: pressures[sim_mode] = numpy.full(len(energies[sim_mode]), numpy.nan) densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] if 'md' in sim_mode and 'langevin' not in sim_mode: - momentum_vector = log_traj['log/md/Integrator/linear_momentum'] + momentum_vector = log_traj[ + 'hoomd-data/md/Integrator/linear_momentum'] linear_momentum[sim_mode] = [ math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) for v in momentum_vector ] @@ -1040,7 +1040,7 @@ def lj_union_compare_modes(*jobs): 'npt_bussi_md_cpu', ] - if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -1048,10 +1048,10 @@ def lj_union_compare_modes(*jobs): 'npt_bussi_md_gpu', ]) - if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -1128,7 +1128,6 @@ def lj_union_compare_modes(*jobs): aggregator=analysis_aggregator) def lj_union_distribution_analyze(*jobs): """Checks that MD follows the correct KE distribution.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -1145,7 +1144,7 @@ def lj_union_distribution_analyze(*jobs): 'npt_bussi_md_cpu', ] - if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_langevin_md_gpu_quantities.h5')): sim_modes.extend([ 'nvt_langevin_md_gpu', 'nvt_mttk_md_gpu', @@ -1153,10 +1152,10 @@ def lj_union_distribution_analyze(*jobs): 'npt_bussi_md_gpu', ]) - if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_cpu_quantities.h5')): sim_modes.extend(['nvt_mc_cpu', 'npt_mc_cpu']) - if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nvt_mc_gpu_quantities.h5')): sim_modes.extend(['nvt_mc_gpu']) util._sort_sim_modes(sim_modes) @@ -1194,13 +1193,13 @@ def lj_union_distribution_analyze(*jobs): n_rotate_dof = num_particles * 3 - print('Reading' + job.fn(sim_mode + '_quantities.gsd')) - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + print('Reading' + job.fn(sim_mode + '_quantities.h5')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) if 'md' in sim_mode: # https://doi.org/10.1371/journal.pone.0202764 ke_translate = log_traj[ - 'log/md/compute/ThermodynamicQuantities/' + 'hoomd-data/md/compute/ThermodynamicQuantities/' 'translational_kinetic_energy'] ke_translate_means_expected[sim_mode].append( numpy.mean(ke_translate) - 1 / 2 * n_translate_dof * kT) @@ -1209,8 +1208,9 @@ def lj_union_distribution_analyze(*jobs): - 1 / math.sqrt(2) * math.sqrt(n_translate_dof) * kT) ke_translate_samples[sim_mode].extend(ke_translate) - ke_rotate = log_traj['log/md/compute/ThermodynamicQuantities/' - 'rotational_kinetic_energy'] + ke_rotate = log_traj[ + 'hoomd-data/md/compute/ThermodynamicQuantities/' + 'rotational_kinetic_energy'] ke_rotate_means_expected[sim_mode].append( numpy.mean(ke_rotate) - 1 / 2 * n_rotate_dof * kT) ke_rotate_sigmas_expected[sim_mode].append( @@ -1225,21 +1225,21 @@ def lj_union_distribution_analyze(*jobs): if 'md' in sim_mode: potential_energy_samples[sim_mode].extend( - log_traj['log/md/compute/ThermodynamicQuantities' + log_traj['hoomd-data/md/compute/ThermodynamicQuantities' '/potential_energy']) else: - potential_energy_samples[sim_mode].extend( - log_traj['log/hpmc/pair/user/CPPPotentialUnion/energy'] - * job.statepoint.kT) + potential_energy_samples[sim_mode].extend(log_traj[ + 'hoomd-data/hpmc/pair/user/CPPPotentialUnion/energy'] + * job.statepoint.kT) if 'md' in sim_mode: - pressure_samples[sim_mode].extend( - log_traj['log/md/compute/ThermodynamicQuantities/pressure']) + pressure_samples[sim_mode].extend(log_traj[ + 'hoomd-data/md/compute/ThermodynamicQuantities/pressure']) else: pressure_samples[sim_mode].extend([job.statepoint.pressure]) density_samples[sim_mode].extend( - log_traj['log/custom_actions/ComputeDensity/density']) + log_traj['hoomd-data/custom_actions/ComputeDensity/density']) ax = fig.add_subplot(4, 2, 1) util.plot_vs_expected(ax, ke_translate_means_expected, @@ -1450,7 +1450,6 @@ def lj_union_nve_md_job(*jobs): aggregator=nve_analysis_aggregator) def lj_union_conservation_analyze(*jobs): """Analyze the output of NVE simulations and inspect conservation.""" - import gsd.hoomd import math import matplotlib import matplotlib.style @@ -1460,7 +1459,7 @@ def lj_union_conservation_analyze(*jobs): print('starting lj_union_conservation_analyze:', jobs[0]) sim_modes = ['nve_md_cpu'] - if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.gsd')): + if os.path.exists(jobs[0].fn('nve_md_gpu_quantities.h5')): sim_modes.extend(['nve_md_gpu']) timesteps = [] @@ -1473,20 +1472,22 @@ def lj_union_conservation_analyze(*jobs): job_linear_momentum = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - job_timesteps[sim_mode] = log_traj['configuration/step'] + job_timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] job_energies[sim_mode] = ( log_traj[ - 'log/md/compute/ThermodynamicQuantities/potential_energy'] + 'hoomd-data/md/compute/ThermodynamicQuantities/potential_energy'] + log_traj[ - 'log/md/compute/ThermodynamicQuantities/kinetic_energy']) + 'hoomd-data/md/compute/ThermodynamicQuantities/kinetic_energy'] + ) job_energies[sim_mode] = ( job_energies[sim_mode] - job_energies[sim_mode][0]) / job.statepoint["num_particles"] - momentum_vector = log_traj['log/md/Integrator/linear_momentum'] + momentum_vector = log_traj[ + 'hoomd-data/md/Integrator/linear_momentum'] job_linear_momentum[sim_mode] = [ math.sqrt(v[0]**2 + v[1]**2 + v[2]**2) / job.statepoint["num_particles"] for v in momentum_vector diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index a75e2193..05dcbca4 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -290,14 +290,13 @@ def make_mc_simulation(job, compute_density = ComputeDensity() sdf = hoomd.hpmc.compute.SDF(xmax=0.02, dx=1e-4) - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, quantities=['translate_moves']) - logger_gsd.add(sdf, quantities=['betaP']) - logger_gsd.add(compute_density, quantities=['density']) - logger_gsd.add(patches, quantities=['energy']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, quantities=['translate_moves']) + logger.add(sdf, quantities=['betaP']) + logger.add(compute_density, quantities=['density']) + logger.add(patches, quantities=['energy']) for loggable, quantity in extra_loggables: - logger_gsd.add(loggable, quantities=[quantity]) + logger.add(loggable, quantities=[quantity]) trajectory_logger = hoomd.logging.Logger(categories=['object']) trajectory_logger.add(mc, quantities=['type_shapes']) @@ -309,7 +308,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -594,7 +593,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def patchy_particle_pressure_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -605,7 +603,7 @@ def patchy_particle_pressure_analyze(job): sim_modes = [] for _ensemble in ['nvt', 'npt']: - if job.isfile(f'{_ensemble}_cpu_quantities.gsd'): + if job.isfile(f'{_ensemble}_cpu_quantities.h5'): sim_modes.append(f'{_ensemble}_cpu') util._sort_sim_modes(sim_modes) @@ -615,14 +613,14 @@ def patchy_particle_pressure_analyze(job): densities = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - timesteps[sim_mode] = log_traj['configuration/step'] + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] - pressures[sim_mode] = log_traj['log/hpmc/compute/SDF/betaP'] + pressures[sim_mode] = log_traj['hoomd-data/hpmc/compute/SDF/betaP'] densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] # save averages for mode in sim_modes: @@ -711,7 +709,7 @@ def patchy_particle_pressure_compare_modes(*jobs): sim_modes = [] for _ensemble in ['nvt', 'npt']: - if jobs[0].isfile(f'{_ensemble}_cpu_quantities.gsd'): + if jobs[0].isfile(f'{_ensemble}_cpu_quantities.h5'): sim_modes.append(f'{_ensemble}_cpu') util._sort_sim_modes(sim_modes) diff --git a/hoomd_validation/simple_polygon.py b/hoomd_validation/simple_polygon.py index 1284edc2..f698f9ea 100644 --- a/hoomd_validation/simple_polygon.py +++ b/hoomd_validation/simple_polygon.py @@ -177,13 +177,12 @@ def make_mc_simulation(job, compute_density = ComputeDensity() sdf = hoomd.hpmc.compute.SDF(xmax=0.02, dx=1e-4) - # log to gsd - logger_gsd = hoomd.logging.Logger(categories=['scalar', 'sequence']) - logger_gsd.add(mc, quantities=['translate_moves']) - logger_gsd.add(sdf, quantities=['betaP']) - logger_gsd.add(compute_density, quantities=['density']) + logger = hoomd.logging.Logger(categories=['scalar', 'sequence']) + logger.add(mc, quantities=['translate_moves']) + logger.add(sdf, quantities=['betaP']) + logger.add(compute_density, quantities=['density']) for loggable, quantity in extra_loggables: - logger_gsd.add(loggable, quantities=[quantity]) + logger.add(loggable, quantities=[quantity]) trajectory_logger = hoomd.logging.Logger(categories=['object']) trajectory_logger.add(mc, quantities=['type_shapes']) @@ -195,7 +194,7 @@ def make_mc_simulation(job, initial_state=initial_state, integrator=mc, sim_mode=sim_mode, - logger=logger_gsd, + logger=logger, table_write_period=WRITE_PERIOD, trajectory_write_period=LOG_PERIOD['trajectory'], log_write_period=LOG_PERIOD['quantities'], @@ -480,7 +479,6 @@ def sampling_operation(*jobs): executable=CONFIG["executable"])) def simple_polygon_analyze(job): """Analyze the output of all simulation modes.""" - import gsd.hoomd import numpy import matplotlib import matplotlib.style @@ -491,7 +489,7 @@ def simple_polygon_analyze(job): sim_modes = [] for _ensemble in ['nvt', 'npt']: - if job.isfile(f'{_ensemble}_cpu_quantities.gsd'): + if job.isfile(f'{_ensemble}_cpu_quantities.h5'): sim_modes.append(f'{_ensemble}_cpu') util._sort_sim_modes(sim_modes) @@ -501,14 +499,14 @@ def simple_polygon_analyze(job): densities = {} for sim_mode in sim_modes: - log_traj = gsd.hoomd.read_log(job.fn(sim_mode + '_quantities.gsd')) + log_traj = util.read_log(job.fn(sim_mode + '_quantities.h5')) - timesteps[sim_mode] = log_traj['configuration/step'] + timesteps[sim_mode] = log_traj['hoomd-data/Simulation/timestep'] - pressures[sim_mode] = log_traj['log/hpmc/compute/SDF/betaP'] + pressures[sim_mode] = log_traj['hoomd-data/hpmc/compute/SDF/betaP'] densities[sim_mode] = log_traj[ - 'log/custom_actions/ComputeDensity/density'] + 'hoomd-data/custom_actions/ComputeDensity/density'] # save averages for mode in sim_modes: @@ -563,7 +561,7 @@ def simple_polygon_compare_modes(*jobs): sim_modes = [] for _ensemble in ['nvt', 'npt']: - if jobs[0].isfile(f'{_ensemble}_cpu_quantities.gsd'): + if jobs[0].isfile(f'{_ensemble}_cpu_quantities.h5'): sim_modes.append(f'{_ensemble}_cpu') util._sort_sim_modes(sim_modes) diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index f442eb7f..fbb4b3f2 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -6,6 +6,7 @@ import numpy import signac import os +import h5py def true_all(*jobs, key): @@ -121,7 +122,7 @@ def make_simulation( output=file) sim.operations.add(table_writer) - # write particle trajectory to gsd file + # write particle trajectory to a gsd file trajectory_writer = hoomd.write.GSD( filename=job.fn(get_job_filename(sim_mode, device, 'trajectory', 'gsd')), @@ -133,16 +134,16 @@ def make_simulation( mode='ab') sim.operations.add(trajectory_writer) - # write logged quantities to gsd file - quantity_writer = hoomd.write.GSD( - filter=hoomd.filter.Null(), - filename=job.fn(get_job_filename(sim_mode, device, 'quantities', - 'gsd')), + # write logged quantities to h5 file + logger.add(sim, quantities=['timestep']) + + quantity_writer = hoomd.write.HDF5Log( + filename=job.fn(get_job_filename(sim_mode, device, 'quantities', 'h5')), trigger=hoomd.trigger.And([ hoomd.trigger.Periodic(log_write_period), hoomd.trigger.After(log_start_step) ]), - mode='ab', + mode='w', logger=logger) sim.operations.add(quantity_writer) @@ -333,3 +334,13 @@ def plot_timeseries(ax, def _sort_sim_modes(sim_modes): """Sort simulation modes for comparison.""" sim_modes.sort(key=lambda x: ('nvt' in x or 'nec' in x, 'md' in x, x)) + + +def read_log(filename): + """Read a HDF5 log as a dictionary of logged quantities.""" + with h5py.File(mode='r', name=filename) as f: + keys = [] + f.visit(lambda name: keys.append(name)) + result = {key: numpy.array(f[key]) for key in keys} + + return result diff --git a/setup.cfg b/setup.cfg index 1618dcd6..6bbeada2 100644 --- a/setup.cfg +++ b/setup.cfg @@ -12,8 +12,8 @@ ignore = D107 exclude = .git, __pycache__ -max_line_length = 80 -max_doc_length = 80 +max_line_length = 88 +max_doc_length = 88 hang_closing = True docstring-convention=google rst-roles =