From b8a2df95ba9b4f74c6a8ac8bf3eb7d9ac3ab3f66 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 19 Sep 2023 09:25:49 -0400 Subject: [PATCH 01/42] List correct branch to base changes on in CONTRIBUTING.md The default branch is `trunk`, not `main`. --- CONTRIBUTING.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 2af1642e..30a0bdfe 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -6,7 +6,7 @@ the project. ## Base your work off the correct branch -All new work should be based on `main`. +All new work should be based on `trunk`. ## Propose a minimal set of related changes From 596834fdc808822b1dca1d12898ec83b9a38cb1b Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 19 Sep 2023 09:29:25 -0400 Subject: [PATCH 02/42] Update PULL_REQUEST_TEMPLATE.md There is no CONTRIBUTORS.md file --- .github/PULL_REQUEST_TEMPLATE.md | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md index 503b1aac..f0268873 100644 --- a/.github/PULL_REQUEST_TEMPLATE.md +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -17,4 +17,3 @@ Resolves #??? - [ ] I have reviewed the [**Contributor Guidelines**](https://github.com/glotzerlab/hoomd-validation/blob/trunk/CONTRIBUTING.md). - [ ] I agree with the terms of the [**HOOMD-blue Contributor Agreement**](https://github.com/glotzerlab/hoomd-validation/blob/trunk/ContributorAgreement.md). -- [ ] My name is on the list of contributors (`CONTRIBUTORS.md`) in the pull request source branch. From fba606e8aba72e349b051641de5e7a3f0c700e96 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 19 Sep 2023 11:34:47 -0400 Subject: [PATCH 03/42] Add patchy-particle-pressure tests --- hoomd_validation/patchy_particle_pressure.py | 668 +++++++++++++++++++ hoomd_validation/util.py | 42 ++ 2 files changed, 710 insertions(+) create mode 100644 hoomd_validation/patchy_particle_pressure.py diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py new file mode 100644 index 00000000..167a6142 --- /dev/null +++ b/hoomd_validation/patchy_particle_pressure.py @@ -0,0 +1,668 @@ +# Copyright (c) 2022-2023 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +"""Test for consistency between NVT and NPT simulations of patchy particles.""" + +from config import CONFIG +from project_class import Project +from flow import aggregator +import util +import os +import json +import pathlib + +# Run parameters shared between simulations. +# Step counts must be even and a multiple of the log quantity period. +RANDOMIZE_STEPS = 20_000 +EQUILIBRATE_STEPS = 100_000 +RUN_STEPS = 2_500_000 +RESTART_STEPS = RUN_STEPS // 50 +TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS + +WRITE_PERIOD = 1_000 +LOG_PERIOD = {'trajectory': 50_000, 'quantities': 100} +NUM_CPU_RANKS = min(8, CONFIG["max_cores_sim"]) + +WALLTIME_STOP_SECONDS = CONFIG["max_walltime"] * 3600 - 10 * 60 + + +def job_statepoints(): + """list(dict): A list of statepoints for this subproject.""" + num_particles = 64**2 + replicate_indices = range(CONFIG["replicates"]) + # statepoint chosen to be in a dense liquid + # FILL IN STATEPOINT AND JUSTIFICATION + params_list = [(0.7, 0.6, 1.0, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ + for temperature, density, pressure, chi, lambda_ in params_list: + for idx in replicate_indices: + yield ({ + "subproject": "patchy_particle_pressure", + "density": density, + "pressure": pressure, + "chi": chi, + "lambda_": lambda_, + "num_particles": num_particles, + "replicate_idx": idx, + "temperature": temperature, + }) + + +def is_patchy_particle_pressure(job): + """Test if a job is part of the patchy_particle_pressure subproject.""" + return job.statepoint['subproject'] == 'patchy_particle_pressure' + + +partition_jobs_cpu_serial = aggregator.groupsof(num=min( + CONFIG["replicates"], CONFIG["max_cores_submission"]), + sort_by='density', + select=is_patchy_particle_pressure,) + +partition_jobs_cpu_mpi = aggregator.groupsof(num=min( + CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure,) + +partition_jobs_gpu = aggregator.groupsof(num=min(CONFIG["replicates"], + CONFIG["max_gpus_submission"]), + sort_by='density', select=is_patchy_particle_pressure,) + + +@Project.post.isfile('patchy_particle_pressure_initial_state.gsd') +@Project.operation(directives=dict( + executable=CONFIG["executable"], + nranks=util.total_ranks_function(NUM_CPU_RANKS), + walltime=1), + aggregator=partition_jobs_cpu_mpi,) +def patchy_particle_pressure_create_initial_state(*jobs): + """Create initial system configuration.""" + import hoomd + import numpy + import itertools + + communicator = hoomd.communicator.Communicator( + ranks_per_partition=NUM_CPU_RANKS) + job = jobs[communicator.partition] + + if communicator.rank == 0: + print('starting patchy_particle_pressure_create_initial_state:', job) + + num_particles = job.statepoint['num_particles'] + density = job.statepoint['density'] + temperature = job.statepoint['temperature'] + chi = job.statepoint['chi'] + lambda_ = job.statepoint['lambda_'] + + box_volume = num_particles / density + L = box_volume**(1 / 3.) + + N = int(numpy.ceil(num_particles**(1. / 3.))) + x = numpy.linspace(-L / 2, L / 2, N, endpoint=False) + + if x[1] - x[0] < 1.0: + raise RuntimeError('density too high to initialize on square lattice') + + position = list(itertools.product(x, repeat=3))[:num_particles] + + # create snapshot + device = hoomd.device.CPU( + communicator=communicator, + message_filename=job.fn('create_initial_state.log')) + snap = hoomd.Snapshot(communicator) + + if communicator.rank == 0: + snap.particles.N = num_particles + snap.particles.types = ['A'] + snap.configuration.box = [L, L, L, 0, 0, 0] + snap.particles.position[:] = position + snap.particles.typeid[:] = [0] * num_particles + + # Use hard sphere + patches Monte Carlo to randomize the initial + # configuration + mc = hoomd.hpmc.integrate.Sphere(default_d=0.05, default_a=0.1) + diameter = 1.0 + mc.shape['A'] = dict(diameter=diameter, orientable=True) + delta = 2 * numpy.arcsin(numpy.sqrt(chi)) + patch_code = util._single_patch_kern_frenkel_code( + delta, lambda_, diameter, temperature) + r_cut = diameter + diameter * (lambda_ - 1) + patches = hoomd.hpmc.pair.user.CPPPotential( + r_cut=r_cut, code=patch_code, param_array=[]) + mc.pair_potential = patches + + sim = hoomd.Simulation(device=device, seed=util.make_seed(job)) + sim.create_state_from_snapshot(snap) + sim.operations.integrator = mc + + device.notice('Randomizing initial state...') + sim.run(RANDOMIZE_STEPS) + device.notice(f'Move counts: {mc.translate_moves} {mc.rotate_moves}') + device.notice('Done.') + + trajectory_logger = hoomd.logging.Logger(categories=['object']) + trajectory_logger.add(mc, quantities=['type_shapes']) + + hoomd.write.GSD.write( + state=sim.state, + filename=job.fn("patchy_particle_pressure_initial_state.gsd"), + mode='wb', + logger=trajectory_logger, + ) + + if communicator.rank == 0: + print(f'completed patchy_particle_pressure_create_initial_state: ' + f'{job} in {communicator.walltime} s') + + +def make_mc_simulation(job, + device, + initial_state, + sim_mode, + extra_loggables=[]): + """Make a patchy sphere MC Simulation. + + Args: + job (`signac.job.Job`): Signac job object. + + device (`hoomd.device.Device`): Device object. + + initial_state (str): Path to the gsd file to be used as an initial state + for the simulation. + + sim_mode (str): String defining the simulation mode. + + extra_loggables (list[tuple]): List of extra loggables to log to gsd + files. Each tuple is a pair of the instance and the loggable + quantity name. + """ + import hoomd + import numpy + from custom_actions import ComputeDensity + + # integrator and patchy potential + temperature = job.statepoint['temperature'] + chi = job.statepoint['chi'] + lambda_ = job.statepoint['lambda_'] + diameter = 1.0 + mc = hoomd.hpmc.integrate.Sphere(default_d=0.05, default_a=0.1) + mc.shape['A'] = dict(diameter=diameter, orientable=True) + delta = 2 * numpy.arcsin(numpy.sqrt(chi)) + patch_code = util._single_patch_kern_frenkel_code( + delta, lambda_, diameter, temperature) + r_cut = diameter + diameter * (lambda_ - 1) + patches = hoomd.hpmc.pair.user.CPPPotential( + r_cut=r_cut, code=patch_code, param_array=[]) + mc.pair_potential = patches + + # compute the density and pressure + 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']) + for loggable, quantity in extra_loggables: + logger_gsd.add(loggable, quantities=[quantity]) + + trajectory_logger = hoomd.logging.Logger(categories=['object']) + trajectory_logger.add(mc, quantities=['type_shapes']) + + # make simulation + sim = util.make_simulation( + job=job, + device=device, + initial_state=initial_state, + integrator=mc, + sim_mode=sim_mode, + logger=logger_gsd, + table_write_period=WRITE_PERIOD, + trajectory_write_period=LOG_PERIOD['trajectory'], + log_write_period=LOG_PERIOD['quantities'], + log_start_step=RANDOMIZE_STEPS + EQUILIBRATE_STEPS, + trajectory_logger=trajectory_logger, + ) + + sim.operations.computes.append(sdf) + compute_density.attach(sim) + + for loggable, _ in extra_loggables: + # call attach method explicitly so we can access simulation state when + # computing the loggable quantity + if hasattr(loggable, 'attach'): + loggable.attach(sim) + + # move size tuner + move_size_tuner = hoomd.hpmc.tune.MoveSize.scale_solver( + moves=['d', 'a'], + target=0.2, + max_translation_move=0.5, + max_rotation_move=2 * numpy.pi / 4, + trigger=hoomd.trigger.And([ + hoomd.trigger.Periodic(100), + hoomd.trigger.Before(RANDOMIZE_STEPS + EQUILIBRATE_STEPS // 2) + ])) + sim.operations.add(move_size_tuner) + + return sim + + +def run_nvt_sim(job, device, complete_filename): + """Run MC sim in NVT.""" + import hoomd + + sim_mode = 'nvt' + restart_filename = util.get_job_filename(sim_mode, device, 'restart', 'gsd') + if job.isfile(restart_filename): + initial_state = job.fn(restart_filename) + restart = True + else: + initial_state = job.fn('patchy_particle_pressure_initial_state.gsd') + restart = False + + sim = make_mc_simulation(job, + device, + initial_state, + sim_mode, + extra_loggables=[]) + + if not restart: + # equilibrate + device.notice('Equilibrating...') + sim.run(EQUILIBRATE_STEPS // 2) + sim.run(EQUILIBRATE_STEPS // 2) + device.notice('Done.') + + # Print acceptance ratio as measured during the 2nd half of the + # equilibration. + translate_moves = sim.operations.integrator.translate_moves + translate_acceptance = translate_moves[0] / sum(translate_moves) + rotate_moves = sim.operations.integrator.rotate_moves + rotate_acceptance = rotate_moves[0] / sum(rotate_moves) + device.notice(f'Translate move acceptance: {translate_acceptance}') + device.notice( + f'Trial translate move size: {sim.operations.integrator.d["A"]}') + device.notice(f'Rotate move acceptance: {rotate_acceptance}') + device.notice( + f'Trial rotate move size: {sim.operations.integrator.a["A"]}') + + # save move size to a file + if device.communicator.rank == 0: + name = util.get_job_filename(sim_mode, device, 'move_size', 'json') + with open(job.fn(name), 'w') as f: + json.dump( + dict( + d_A=sim.operations.integrator.d["A"], + a_A=sim.operations.integrator.a["A"], + ), f) + else: + device.notice('Restarting...') + # read move size from the file + name = util.get_job_filename(sim_mode, device, 'move_size', 'json') + with open(job.fn(name), 'r') as f: + data = json.load(f) + + sim.operations.integrator.d["A"] = data['d_A'] + sim.operations.integrator.a["A"] = data['a_A'] + mcd = sim.operations.integrator.d["A"] + device.notice(f'Restored translate trial move size: {mcd}') + mca = sim.operations.integrator.a["A"] + device.notice(f'Restored rotate trial move size: {mca}') + + # run + device.notice('Running...') + + util.run_up_to_walltime(sim=sim, + end_step=TOTAL_STEPS, + steps=RESTART_STEPS, + walltime_stop=WALLTIME_STOP_SECONDS) + + hoomd.write.GSD.write(state=sim.state, + filename=job.fn(restart_filename), + mode='wb') + + if sim.timestep == TOTAL_STEPS: + pathlib.Path(job.fn(complete_filename)).touch() + device.notice('Done.') + else: + device.notice('Ending run early due to walltime limits at:' + f'{device.communicator.walltime}') + + +def run_npt_sim(job, device, complete_filename): + """Run MC sim in NPT.""" + import hoomd + + # device + sim_mode = 'npt' + restart_filename = util.get_job_filename(sim_mode, device, 'restart', 'gsd') + if job.isfile(restart_filename): + initial_state = job.fn(restart_filename) + restart = True + else: + initial_state = job.fn('patchy_particle_pressure_initial_state.gsd') + restart = False + + # box updates + boxmc = hoomd.hpmc.update.BoxMC(betaP=job.statepoint.pressure, + trigger=hoomd.trigger.Periodic(1)) + boxmc.volume = dict(weight=1.0, mode='ln', delta=1e-6) + + # simulation + sim = make_mc_simulation(job, + device, + initial_state, + sim_mode, + extra_loggables=[ + (boxmc, 'volume_moves'), + ]) + + sim.operations.add(boxmc) + + boxmc_tuner = hoomd.hpmc.tune.BoxMCMoveSize.scale_solver( + trigger=hoomd.trigger.And([ + hoomd.trigger.Periodic(400), + hoomd.trigger.Before(RANDOMIZE_STEPS + EQUILIBRATE_STEPS // 2) + ]), + boxmc=boxmc, + moves=['volume'], + target=0.5) + sim.operations.add(boxmc_tuner) + + if not restart: + # equilibrate + device.notice('Equilibrating...') + sim.run(EQUILIBRATE_STEPS // 2) + sim.run(EQUILIBRATE_STEPS // 2) + device.notice('Done.') + + # Print acceptance ratio as measured during the 2nd half of the + # equilibration. + translate_moves = sim.operations.integrator.translate_moves + translate_acceptance = translate_moves[0] / sum(translate_moves) + device.notice(f'Translate move acceptance: {translate_acceptance}') + device.notice( + f'Translate trial move size: {sim.operations.integrator.d["A"]}') + rotate_moves = sim.operations.integrator.rotate_moves + rotate_acceptance = rotate_moves[0] / sum(rotate_moves) + device.notice(f'Rotate move acceptance: {rotate_acceptance}') + device.notice( + f'Rotate trial move size: {sim.operations.integrator.a["A"]}') + + volume_moves = boxmc.volume_moves + volume_acceptance = volume_moves[0] / sum(volume_moves) + device.notice(f'Volume move acceptance: {volume_acceptance}') + device.notice(f'Volume move size: {boxmc.volume["delta"]}') + + # save move sizes to a file + if device.communicator.rank == 0: + name = util.get_job_filename(sim_mode, device, 'move_size', 'json') + with open(job.fn(name), 'w') as f: + json.dump( + dict(d_A=sim.operations.integrator.d["A"], + a_A=sim.operations.integrator.a["A"], + volume_delta=boxmc.volume['delta']), f) + else: + device.notice('Restarting...') + # read move size from the file + name = util.get_job_filename(sim_mode, device, 'move_size', 'json') + with open(job.fn(name), 'r') as f: + data = json.load(f) + + sim.operations.integrator.d["A"] = data['d_A'] + sim.operations.integrator.a["A"] = data['a_A'] + mcd = sim.operations.integrator.d["A"] + device.notice(f'Restored translate trial move size: {mcd}') + mca = sim.operations.integrator.a["A"] + device.notice(f'Restored rotate trial move size: {mca}') + boxmc.volume = dict(weight=1.0, mode='ln', delta=data['volume_delta']) + device.notice(f'Restored volume move size: {boxmc.volume["delta"]}') + + # run + device.notice('Running...') + util.run_up_to_walltime(sim=sim, + end_step=TOTAL_STEPS, + steps=100_000, + walltime_stop=WALLTIME_STOP_SECONDS) + + hoomd.write.GSD.write(state=sim.state, + filename=job.fn(restart_filename), + mode='wb') + + if sim.timestep == TOTAL_STEPS: + pathlib.Path(job.fn(complete_filename)).touch() + device.notice('Done.') + else: + device.notice('Ending run early due to walltime limits at:' + f'{device.communicator.walltime}') + + +sampling_jobs = [] +job_definitions = [ + { + 'mode': 'nvt', + 'device_name': 'cpu', + 'ranks_per_partition': NUM_CPU_RANKS, + 'aggregator': partition_jobs_cpu_mpi + }, + { + 'mode': 'npt', + 'device_name': 'cpu', + 'ranks_per_partition': NUM_CPU_RANKS, + 'aggregator': partition_jobs_cpu_mpi + }, +] + +job_definitions = [] +if CONFIG["enable_gpu"]: + job_definitions.extend([ + { + 'mode': 'nvt', + 'device_name': 'gpu', + 'ranks_per_partition': 1, + 'aggregator': partition_jobs_gpu + }, + { + 'mode': 'npt', + 'device_name': 'gpu', + 'ranks_per_partition': 1, + 'aggregator': partition_jobs_gpu + }, + ]) + + +def add_sampling_job(mode, device_name, ranks_per_partition, aggregator): + """Add a sampling job to the workflow.""" + directives = dict(walltime=CONFIG["max_walltime"], + executable=CONFIG["executable"], + nranks=util.total_ranks_function(ranks_per_partition)) + + if device_name == 'gpu': + directives['ngpu'] = directives['nranks'] + + @Project.pre.after(patchy_particle_pressure_create_initial_state) + @Project.post.isfile(f'{mode}_{device_name}_complete') + @Project.operation(name=f'patchy_particle_pressure_{mode}_{device_name}', + directives=directives, + aggregator=aggregator) + def sampling_operation(*jobs): + """Perform sampling simulation given the definition.""" + import hoomd + + communicator = hoomd.communicator.Communicator( + ranks_per_partition=ranks_per_partition) + job = jobs[communicator.partition] + + if communicator.rank == 0: + print(f'starting patchy_particle_pressure_{mode}_{device_name}:', job) + + if device_name == 'gpu': + device_cls = hoomd.device.GPU + elif device_name == 'cpu': + device_cls = hoomd.device.CPU + + device = device_cls( + communicator=communicator, + message_filename=job.fn(f'{mode}_{device_name}.log')) + + globals().get(f'run_{mode}_sim')( + job, device, complete_filename=f'{mode}_{device_name}_complete') + + if communicator.rank == 0: + print(f'completed patchy_particle_pressure_{mode}_{device_name} ' + f'{job} in {communicator.walltime} s') + + sampling_jobs.append(sampling_operation) + + +for definition in job_definitions: + add_sampling_job(**definition) + + +@Project.pre(is_patchy_particle_pressure) +@Project.pre.after(*sampling_jobs) +@Project.post.true('patchy_particle_pressure_analysis_complete') +@Project.operation(directives=dict(walltime=CONFIG['short_walltime'], + 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 + import matplotlib.figure + matplotlib.style.use('fivethirtyeight') + + print('starting patchy_particle_pressure_analyze:', job) + + sim_modes = [] + for _ensemble in ['nvt', 'npt']: + for _device in ['cpu', 'gpu']: + if job.isfile(f'{_ensemble}_{_device}_quantities.gsd'): + sim_modes.append(f'{_ensemble}_{_device}') + + util._sort_sim_modes(sim_modes) + + timesteps = {} + pressures = {} + 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'] + + pressures[sim_mode] = log_traj['log/hpmc/compute/SDF/betaP'] + + densities[sim_mode] = log_traj[ + 'log/custom_actions/ComputeDensity/density'] + + # save averages + for mode in sim_modes: + job.document[mode] = dict(pressure=float(numpy.mean(pressures[mode])), + density=float(numpy.mean(densities[mode]))) + + # Plot results + fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') + ax = fig.add_subplot(2, 1, 1) + util.plot_timeseries(ax=ax, + timesteps=timesteps, + data=densities, + ylabel=r"$\rho$", + expected=job.sp.density, + max_points=500) + ax.legend() + + ax = fig.add_subplot(2, 1, 2) + util.plot_timeseries(ax=ax, + timesteps=timesteps, + data=pressures, + ylabel=r"$\beta P$", + expected=job.sp.pressure, + max_points=500) + + fig.suptitle(f"$\\rho={job.statepoint.density}$, " + f"$N={job.statepoint.num_particles}$, " + f"replicate={job.statepoint.replicate_idx}") + fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight') + + job.document['patchy_particle_pressure_analysis_complete'] = True + + +@Project.pre( + lambda *jobs: util.true_all(*jobs, key='patchy_particle_pressure_analysis_complete')) +@Project.post(lambda *jobs: util.true_all( + *jobs, key='patchy_particle_pressure_compare_modes_complete')) +@Project.operation(directives=dict(executable=CONFIG["executable"]), + aggregator=aggregator.groupby( + key=['density', 'num_particles'], + sort_by='replicate_idx', + select=is_patchy_particle_pressure)) +def patchy_particle_pressure_compare_modes(*jobs): + """Compares the tested simulation modes.""" + import numpy + import matplotlib + import matplotlib.style + import matplotlib.figure + matplotlib.style.use('fivethirtyeight') + + print('starting patchy_particle_pressure_compare_modes:', jobs[0]) + + sim_modes = [] + for _ensemble in ['nvt', 'npt']: + for _device in ['cpu', 'gpu']: + if jobs[0].isfile(f'{_ensemble}_{_device}_quantities.gsd'): + sim_modes.append(f'{_ensemble}_{_device}') + + util._sort_sim_modes(sim_modes) + + quantity_names = ['density', 'pressure'] + + labels = { + 'density': r'$\frac{\rho_\mathrm{sample} - \rho}{\rho} \cdot 1000$', + 'pressure': r'$\frac{P_\mathrm{sample} - P}{P} \cdot 1000$', + } + + # grab the common statepoint parameters + set_density = jobs[0].sp.density + set_pressure = jobs[0].sp.pressure + num_particles = jobs[0].sp.num_particles + + quantity_reference = dict(density=set_density, pressure=set_pressure) + + fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') + fig.suptitle(f"$\\rho={set_density}$, $N={num_particles}$") + + for i, quantity_name in enumerate(quantity_names): + ax = fig.add_subplot(2, 1, i + 1) + + # organize data from jobs + quantities = {mode: [] for mode in sim_modes} + for jb in jobs: + for mode in sim_modes: + quantities[mode].append( + getattr(getattr(jb.doc, mode), quantity_name)) + + if quantity_reference[quantity_name] is not None: + reference = quantity_reference[quantity_name] + else: + avg_value = { + mode: numpy.mean(quantities[mode]) for mode in sim_modes + } + reference = numpy.mean([avg_value[mode] for mode in sim_modes]) + + avg_quantity, stderr_quantity = util.plot_vs_expected( + ax=ax, + values=quantities, + ylabel=labels[quantity_name], + expected=reference, + relative_scale=1000, + separate_nvt_npt=True) + + filename = f'patchy_particle_pressure_compare_density{round(set_density, 2)}.svg' + fig.savefig(os.path.join(jobs[0]._project.path, filename), + bbox_inches='tight') + + for job in jobs: + job.document['patchy_particle_pressure_compare_modes_complete'] = True diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index 770ec375..9a39e058 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -299,3 +299,45 @@ 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 _single_patch_kern_frenkel_code(delta_rad, lambda_, sigma, kT): + """Generate code for JIT compilation of Kern-Frenkel potential. + + Args: + delta_rad (float): Half-opening angle of patchy interaction in radians + + lambda_ (float): range of patchy interaction relative to hard core + diameter + + sigma (float): Diameter of hard core of particles + + kT (float): Temperature; sets the energy scale + + The terminology (e.g., `ehat`) comes from the "Modelling Patchy Particles" + HOOMD-blue tutorial. + + """ + patch_code = f""" + const float delta = {delta_rad}; + const float lambda = {lambda_:f}; + const float sigma = {sigma:f}; // hard core diameter + const float kT = {kT:f}; + const vec3 ehat_particle_reference_frame(1, 0, 0); + vec3 ehat_i = rotate(q_i, ehat_particle_reference_frame); + vec3 ehat_j = rotate(q_j, ehat_particle_reference_frame); + vec3 r_hat_ij = r_ij / sqrtf(dot(r_ij, r_ij)); + bool patch_on_i_is_aligned_with_r_ij = dot(ehat_i, r_hat_ij) >= cos(delta); + bool patch_on_j_is_aligned_with_r_ji = dot(ehat_j, -r_hat_ij) >= cos(delta); + float rsq = dot(r_ij, r_ij); + if (patch_on_i_is_aligned_with_r_ij + && patch_on_j_is_aligned_with_r_ji + && dot(r_ij, r_ij) < lambda*sigma*lambda*sigma) + {{ + return -1 / kT; + }} + else + {{ + return 0.0; + }} + """ From 91d4f865257bf900cf1356195ab345d67224209b Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 19 Sep 2023 11:38:55 -0400 Subject: [PATCH 04/42] Add patchy_particle_pressure to init and project files --- hoomd_validation/init.py | 2 ++ hoomd_validation/project.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/hoomd_validation/init.py b/hoomd_validation/init.py index 9166d174..a4f4b203 100644 --- a/hoomd_validation/init.py +++ b/hoomd_validation/init.py @@ -14,6 +14,7 @@ import hard_disk import hard_sphere import simple_polygon +import patchy_particle_pressure subprojects = [ alj_2d, @@ -22,6 +23,7 @@ hard_disk, hard_sphere, simple_polygon, + patchy_particle_pressure ] project = signac.init_project(path=config.project_root) diff --git a/hoomd_validation/project.py b/hoomd_validation/project.py index 54dd00db..c5adff66 100644 --- a/hoomd_validation/project.py +++ b/hoomd_validation/project.py @@ -14,6 +14,7 @@ import hard_disk import hard_sphere import simple_polygon +import patchy_particle_pressure # use srun on delta (mpiexec fails on multiple nodes) flow.environments.xsede.DeltaEnvironment.mpi_cmd = "srun" @@ -25,6 +26,7 @@ "hard_disk", "hard_sphere", "simple_polygon", + "patchy_particle_pressure", ] if __name__ == "__main__": From c9b49dcc7361cb294983a5432e0c1b804710e3a5 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 17:04:24 -0400 Subject: [PATCH 05/42] Fix typos --- hoomd_validation/init.py | 12 ++++++------ hoomd_validation/patchy_particle_pressure.py | 4 ++-- hoomd_validation/util.py | 1 + 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/hoomd_validation/init.py b/hoomd_validation/init.py index a4f4b203..c7bd2774 100644 --- a/hoomd_validation/init.py +++ b/hoomd_validation/init.py @@ -17,12 +17,12 @@ import patchy_particle_pressure subprojects = [ - alj_2d, - lj_fluid, - lj_union, - hard_disk, - hard_sphere, - simple_polygon, +# alj_2d, +# lj_fluid, +# lj_union, +# hard_disk, +# hard_sphere, +# simple_polygon, patchy_particle_pressure ] diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 167a6142..c728382f 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -122,7 +122,7 @@ def patchy_particle_pressure_create_initial_state(*jobs): diameter = 1.0 mc.shape['A'] = dict(diameter=diameter, orientable=True) delta = 2 * numpy.arcsin(numpy.sqrt(chi)) - patch_code = util._single_patch_kern_frenkel_code( + patch_code = util._single_patch_kern_frenkel_code( delta, lambda_, diameter, temperature) r_cut = diameter + diameter * (lambda_ - 1) patches = hoomd.hpmc.pair.user.CPPPotential( @@ -454,7 +454,7 @@ def run_npt_sim(job, device, complete_filename): }, ] -job_definitions = [] +#job_definitions = [] if CONFIG["enable_gpu"]: job_definitions.extend([ { diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index 9a39e058..54999e06 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -341,3 +341,4 @@ def _single_patch_kern_frenkel_code(delta_rad, lambda_, sigma, kT): return 0.0; }} """ + return patch_code From 82be5bb245dfc7a6d23783cd967bac4831aeb818 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 18:03:26 -0400 Subject: [PATCH 06/42] Add pressure --- hoomd_validation/patchy_particle_pressure.py | 4 ++-- hoomd_validation/util.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index c728382f..7be4e37a 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -15,7 +15,7 @@ # Step counts must be even and a multiple of the log quantity period. RANDOMIZE_STEPS = 20_000 EQUILIBRATE_STEPS = 100_000 -RUN_STEPS = 2_500_000 +RUN_STEPS = 1_000_000 RESTART_STEPS = RUN_STEPS // 50 TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS @@ -32,7 +32,7 @@ def job_statepoints(): replicate_indices = range(CONFIG["replicates"]) # statepoint chosen to be in a dense liquid # FILL IN STATEPOINT AND JUSTIFICATION - params_list = [(0.7, 0.6, 1.0, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ + params_list = [(0.7, 0.6, -1.2085475196748359, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: yield ({ diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index 54999e06..9385fc53 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -28,7 +28,7 @@ def get_job_filename(sim_mode, device, name, type): return f"{sim_mode}_{suffix}_{name}.{type}" -def run_up_to_walltime(sim, end_step, steps, walltime_stop): +def run_up_to_walltime(sim, end_step, steps, walltime_stop, minutes_per_run=None): """Run a simulation, stopping early if a walltime limit is reached. Args: From 387ecde6d949b3b4ef416e945c6d46bc2206273b Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 18:15:33 -0400 Subject: [PATCH 07/42] test --- hoomd_validation/patchy_particle_pressure.py | 1 + 1 file changed, 1 insertion(+) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 7be4e37a..1c91580e 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -318,6 +318,7 @@ def run_nvt_sim(job, device, complete_filename): steps=RESTART_STEPS, walltime_stop=WALLTIME_STOP_SECONDS) + hoomd.write.GSD.write(state=sim.state, filename=job.fn(restart_filename), mode='wb') From 2a3868545293f5017c0bc8120dd9ec0a8a658b70 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 18:16:25 -0400 Subject: [PATCH 08/42] undo test --- hoomd_validation/patchy_particle_pressure.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 1c91580e..7be4e37a 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -318,7 +318,6 @@ def run_nvt_sim(job, device, complete_filename): steps=RESTART_STEPS, walltime_stop=WALLTIME_STOP_SECONDS) - hoomd.write.GSD.write(state=sim.state, filename=job.fn(restart_filename), mode='wb') From 495aff900d18260a8fc587f3e9381eb3296e12fc Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 18:17:04 -0400 Subject: [PATCH 09/42] test --- hoomd_validation/patchy_particle_pressure.py | 1 + 1 file changed, 1 insertion(+) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 7be4e37a..1c91580e 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -318,6 +318,7 @@ def run_nvt_sim(job, device, complete_filename): steps=RESTART_STEPS, walltime_stop=WALLTIME_STOP_SECONDS) + hoomd.write.GSD.write(state=sim.state, filename=job.fn(restart_filename), mode='wb') From 3f4d5e69d683872bd652f52ad9bc9ffa237db55a Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 21 Sep 2023 18:17:32 -0400 Subject: [PATCH 10/42] Undo test --- hoomd_validation/patchy_particle_pressure.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 1c91580e..7be4e37a 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -318,7 +318,6 @@ def run_nvt_sim(job, device, complete_filename): steps=RESTART_STEPS, walltime_stop=WALLTIME_STOP_SECONDS) - hoomd.write.GSD.write(state=sim.state, filename=job.fn(restart_filename), mode='wb') From e4ba892187a745a97fcd92416faf0e483113e4e4 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Fri, 22 Sep 2023 12:59:24 -0400 Subject: [PATCH 11/42] Change min number of cpu ranks --- hoomd_validation/patchy_particle_pressure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 7be4e37a..56a3fa91 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -21,7 +21,7 @@ WRITE_PERIOD = 1_000 LOG_PERIOD = {'trajectory': 50_000, 'quantities': 100} -NUM_CPU_RANKS = min(8, CONFIG["max_cores_sim"]) +NUM_CPU_RANKS = min(9, CONFIG["max_cores_sim"]) WALLTIME_STOP_SECONDS = CONFIG["max_walltime"] * 3600 - 10 * 60 From 4fd71403de805bef0964fa333ee8d4e0e35b2876 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Fri, 22 Sep 2023 21:48:09 -0400 Subject: [PATCH 12/42] Add contributors line back to PR template --- .github/PULL_REQUEST_TEMPLATE.md | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md index f0268873..503b1aac 100644 --- a/.github/PULL_REQUEST_TEMPLATE.md +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -17,3 +17,4 @@ Resolves #??? - [ ] I have reviewed the [**Contributor Guidelines**](https://github.com/glotzerlab/hoomd-validation/blob/trunk/CONTRIBUTING.md). - [ ] I agree with the terms of the [**HOOMD-blue Contributor Agreement**](https://github.com/glotzerlab/hoomd-validation/blob/trunk/ContributorAgreement.md). +- [ ] My name is on the list of contributors (`CONTRIBUTORS.md`) in the pull request source branch. From 6531fadb00ad65574423ec001a22bd425f1010bc Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Fri, 22 Sep 2023 21:48:46 -0400 Subject: [PATCH 13/42] Add mailmap --- .mailmap | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 .mailmap diff --git a/.mailmap b/.mailmap new file mode 100644 index 00000000..36fb78bf --- /dev/null +++ b/.mailmap @@ -0,0 +1,3 @@ +Tim Moore Tim Moore +Tommy Waltmann Tommy Waltmann <53307607+tommy-waltmann@users.noreply.github.com> +Tommy Waltmann tommy-waltmann <53307607+tommy-waltmann@users.noreply.github.com> From 9c4933c1d8bed8a9945730b5e155fa5b23a97648 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Fri, 22 Sep 2023 22:05:45 -0400 Subject: [PATCH 14/42] Add CONTRIBUTORS.md --- CONTRIBUTORS.md | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 CONTRIBUTORS.md diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md new file mode 100644 index 00000000..0c0a8e12 --- /dev/null +++ b/CONTRIBUTORS.md @@ -0,0 +1,8 @@ +# Contributors + +The following people have contributed to HOOMD-validation: + +* Joshua A. Anderson, University of Michigan +* Tommy Waltmann, University of Michigan +* Tim Moore, University of Michigan +* Brandon Butler, University of Michigan From e847901275c7282f1313a62562af5d5820641f53 Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Mon, 25 Sep 2023 09:20:41 -0400 Subject: [PATCH 15/42] Update CONTRIBUTORS.md --- CONTRIBUTORS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 0c0a8e12..0500dcaf 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -1,6 +1,6 @@ # Contributors -The following people have contributed to HOOMD-validation: +The following people have contributed to hoomd-validation: * Joshua A. Anderson, University of Michigan * Tommy Waltmann, University of Michigan From cab505548279c1032d8f68203e13444a27bec266 Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Mon, 25 Sep 2023 09:28:54 -0400 Subject: [PATCH 16/42] Fix dependabot configuration. --- .github/dependabot.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 8d64d26f..aa3347c1 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -2,7 +2,7 @@ version: 2 updates: - package-ecosystem: "github-actions" directory: "/" - target-branch: main + target-branch: trunk schedule: interval: "monthly" pull-request-branch-name: @@ -11,7 +11,7 @@ updates: - package-ecosystem: "pip" directory: ".github/workflows" - target-branch: main + target-branch: trunk schedule: interval: "monthly" pull-request-branch-name: From 29ba09a1e34fc6204c2de10a89470a58acc84eed Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Mon, 25 Sep 2023 09:29:14 -0400 Subject: [PATCH 17/42] Ignore code-workspace files. --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 243a6810..be5f0b2c 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,4 @@ signac_project_document.json .signac_sp_cache.json.gz __pycache__ .signac +*.code-workspace From 3da5766b7e30fc88fe698f35c875974d0578754c Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Mon, 25 Sep 2023 09:29:25 -0400 Subject: [PATCH 18/42] Update actions used in CI tests. --- .github/workflows/CI.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index fd8855b9..b2e33f65 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,12 +22,12 @@ jobs: name: flow-status runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3.5.2 + - uses: actions/checkout@v4.0.0 - name: Set up Python - uses: actions/setup-python@v4.6.1 + uses: actions/setup-python@v4.7.0 with: - python-version: '3.10' - - uses: actions/cache@v3.3.1 + python-version: '3.11' + - uses: actions/cache@v3.3.2 with: path: ~/.cache/pip key: ${{ runner.os }}-pip-${{ hashFiles('.github/workflows/requirements-test.txt') }} From cfcd54afa30313bb2e9faa4c86ae84d4e3d98c58 Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Mon, 25 Sep 2023 09:30:35 -0400 Subject: [PATCH 19/42] Bump python packages used in testing. --- .github/workflows/requirements-test.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/requirements-test.txt b/.github/workflows/requirements-test.txt index 9b59446f..f4473234 100644 --- a/.github/workflows/requirements-test.txt +++ b/.github/workflows/requirements-test.txt @@ -1,5 +1,5 @@ -signac==2.0.0 -signac-flow==0.25.1 -numpy==1.24.3 -PyYAML==6.0 -gsd==2.9.0 +signac==2.1.0 +signac-flow==0.26.1 +numpy==1.26.0 +PyYAML==6.0.1 +gsd==3.1.1 From 82745310ac8a7360218fac642b04afc134354539 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 25 Sep 2023 14:22:36 +0000 Subject: [PATCH 20/42] Bump actions/checkout from 4.0.0 to 4.1.0 Bumps [actions/checkout](https://github.com/actions/checkout) from 4.0.0 to 4.1.0. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v4.0.0...v4.1.0) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index b2e33f65..1ab79d8d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,7 +22,7 @@ jobs: name: flow-status runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4.0.0 + - uses: actions/checkout@v4.1.0 - name: Set up Python uses: actions/setup-python@v4.7.0 with: From f307199dd9978188f524dc6e99365af34a3b2b03 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Mon, 25 Sep 2023 17:24:08 -0400 Subject: [PATCH 21/42] Update statepoint based on NVT sims --- hoomd_validation/patchy_particle_pressure.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 56a3fa91..5548e076 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -31,8 +31,9 @@ def job_statepoints(): num_particles = 64**2 replicate_indices = range(CONFIG["replicates"]) # statepoint chosen to be in a dense liquid - # FILL IN STATEPOINT AND JUSTIFICATION - params_list = [(0.7, 0.6, -1.2085475196748359, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ + # nvt simulations at density = 0.6 yielded a measured pressure of + # -1.218 +/- 0.002 (mean +/- std. error of means) over 32 replicas + params_list = [(0.7, 0.6, -1.218, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: yield ({ From a74838a20e890e0a0bf5397babf790073ad7c036 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 1 Oct 2023 00:07:32 +0000 Subject: [PATCH 22/42] Bump gsd from 3.1.1 to 3.2.0 in /.github/workflows Bumps [gsd](https://github.com/glotzerlab/gsd) from 3.1.1 to 3.2.0. - [Release notes](https://github.com/glotzerlab/gsd/releases) - [Changelog](https://github.com/glotzerlab/gsd/blob/trunk-patch/CHANGELOG.rst) - [Commits](https://github.com/glotzerlab/gsd/compare/v3.1.1...v3.2.0) --- updated-dependencies: - dependency-name: gsd dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/requirements-test.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/requirements-test.txt b/.github/workflows/requirements-test.txt index f4473234..e7a2e2b0 100644 --- a/.github/workflows/requirements-test.txt +++ b/.github/workflows/requirements-test.txt @@ -2,4 +2,4 @@ signac==2.1.0 signac-flow==0.26.1 numpy==1.26.0 PyYAML==6.0.1 -gsd==3.1.1 +gsd==3.2.0 From 6bc2cbeb48a428ee0e9252c3674e1195d373c90f Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 4 Oct 2023 12:47:51 -0400 Subject: [PATCH 23/42] Update dashboard --- dashboard.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/dashboard.py b/dashboard.py index ab309f7b..6237daf0 100644 --- a/dashboard.py +++ b/dashboard.py @@ -32,11 +32,21 @@ def job_title(self, job): f"rho={job.statepoint.density}" elif (job.statepoint.subproject == 'hard_disk' or job.statepoint.subproject == 'hard_sphere' - or job.statepoint.subproject == 'simple_polygon'): + or job.statepoint.subproject == 'simple_polygon' + or job.statepoint.subproject == 'patchy_particle_pressure'): return f"{job.statepoint.subproject}: rho={job.statepoint.density}" else: raise RuntimeError("Unexpected job") + def job_sorter(self, job): + return ( + job.sp.density, + job.sp.pressure, + job.sp.temperature, + job.sp.chi, + job.sp.replicate_idx, + ) + if __name__ == "__main__": ValidationDashboard(modules=modules).main() From 2293b2e46e1dd5483f3feb3a2a2c23d15ccf7d71 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 4 Oct 2023 12:48:40 -0400 Subject: [PATCH 24/42] Update job initialization and aggregation --- hoomd_validation/patchy_particle_pressure.py | 106 +++++++++++++------ 1 file changed, 73 insertions(+), 33 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 56a3fa91..8ae3de72 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -20,7 +20,7 @@ TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS WRITE_PERIOD = 1_000 -LOG_PERIOD = {'trajectory': 50_000, 'quantities': 100} +LOG_PERIOD = {'trajectory': 50_000, 'quantities': 500} NUM_CPU_RANKS = min(9, CONFIG["max_cores_sim"]) WALLTIME_STOP_SECONDS = CONFIG["max_walltime"] * 3600 - 10 * 60 @@ -28,11 +28,27 @@ def job_statepoints(): """list(dict): A list of statepoints for this subproject.""" - num_particles = 64**2 + num_particles = 8**3 replicate_indices = range(CONFIG["replicates"]) # statepoint chosen to be in a dense liquid - # FILL IN STATEPOINT AND JUSTIFICATION - params_list = [(0.7, 0.6, -1.2085475196748359, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ + # nvt simulations at density = 0.95 yielded a measured pressure of + # 6.415 +/- 0.003 (mean +/- std. error of means) over 32 replicas + params_list = [ + (1000.0, 0.95, 12.12, 0.7, 1.5), + (10.0, 0.95, 11.83, 0.7, 1.5), + (1.0, 0.95, 8.27, 0.7, 1.5), + # next 3 are from 10.1063/1.3054361 + # pressure from NVT simulations + (0.5714, 0.8, -5.871736803897682, 1.0, 1.5), + (1.0, 0.8, -1.0235450460184001, 1.0, 1.5), + (3.0, 0.8, 3.7415391338219885, 1.0, 1.5), + # next 2 are from 10.1063/1.3054361 + # pressure from webplotdigitizer fig 7 + (1.0, 0.8, 1.70178459919393820, 1.0, 1.5 + ), # pressure = pressure = pressure/kT + (3.0, 0.8, 13.549010655204555, 1.0, 1.5), # pressure = pressure + (3.0, 0.8, 4.516336885068185, 1.0, 1.5), # pressure = pressure/kT + ] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: yield ({ @@ -52,27 +68,47 @@ def is_patchy_particle_pressure(job): return job.statepoint['subproject'] == 'patchy_particle_pressure' -partition_jobs_cpu_serial = aggregator.groupsof(num=min( - CONFIG["replicates"], CONFIG["max_cores_submission"]), - sort_by='density', - select=is_patchy_particle_pressure,) - -partition_jobs_cpu_mpi = aggregator.groupsof(num=min( - CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), - sort_by='density', - select=is_patchy_particle_pressure,) - -partition_jobs_gpu = aggregator.groupsof(num=min(CONFIG["replicates"], - CONFIG["max_gpus_submission"]), - sort_by='density', select=is_patchy_particle_pressure,) +def is_patchy_particle_pressure_positive_pressure(job): + """Test if a job is part of the patchy_particle_pressure subproject.""" + return job.statepoint[ + 'subproject'] == 'patchy_particle_pressure' and job.statepoint[ + 'pressure'] > 0.0 + + +partition_jobs_cpu_serial = aggregator.groupsof( + num=min(CONFIG["replicates"], CONFIG["max_cores_submission"]), + sort_by='density', + select=is_patchy_particle_pressure, +) + +partition_jobs_cpu_mpi_nvt = aggregator.groupsof( + num=min(CONFIG["replicates"], + CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure, +) + +partition_jobs_cpu_mpi_npt = aggregator.groupsof( + num=min(CONFIG["replicates"], + CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure_positive_pressure, +) + +partition_jobs_gpu = aggregator.groupsof( + num=min(CONFIG["replicates"], CONFIG["max_gpus_submission"]), + sort_by='density', + select=is_patchy_particle_pressure, +) @Project.post.isfile('patchy_particle_pressure_initial_state.gsd') -@Project.operation(directives=dict( - executable=CONFIG["executable"], - nranks=util.total_ranks_function(NUM_CPU_RANKS), - walltime=1), - aggregator=partition_jobs_cpu_mpi,) +@Project.operation( + directives=dict(executable=CONFIG["executable"], + nranks=util.total_ranks_function(NUM_CPU_RANKS), + walltime=1), + aggregator=partition_jobs_cpu_mpi_nvt, +) def patchy_particle_pressure_create_initial_state(*jobs): """Create initial system configuration.""" import hoomd @@ -122,11 +158,12 @@ def patchy_particle_pressure_create_initial_state(*jobs): diameter = 1.0 mc.shape['A'] = dict(diameter=diameter, orientable=True) delta = 2 * numpy.arcsin(numpy.sqrt(chi)) - patch_code = util._single_patch_kern_frenkel_code( - delta, lambda_, diameter, temperature) + patch_code = util._single_patch_kern_frenkel_code(delta, lambda_, diameter, + temperature) r_cut = diameter + diameter * (lambda_ - 1) - patches = hoomd.hpmc.pair.user.CPPPotential( - r_cut=r_cut, code=patch_code, param_array=[]) + patches = hoomd.hpmc.pair.user.CPPPotential(r_cut=r_cut, + code=patch_code, + param_array=[]) mc.pair_potential = patches sim = hoomd.Simulation(device=device, seed=util.make_seed(job)) @@ -186,11 +223,12 @@ def make_mc_simulation(job, mc = hoomd.hpmc.integrate.Sphere(default_d=0.05, default_a=0.1) mc.shape['A'] = dict(diameter=diameter, orientable=True) delta = 2 * numpy.arcsin(numpy.sqrt(chi)) - patch_code = util._single_patch_kern_frenkel_code( - delta, lambda_, diameter, temperature) + patch_code = util._single_patch_kern_frenkel_code(delta, lambda_, diameter, + temperature) r_cut = diameter + diameter * (lambda_ - 1) - patches = hoomd.hpmc.pair.user.CPPPotential( - r_cut=r_cut, code=patch_code, param_array=[]) + patches = hoomd.hpmc.pair.user.CPPPotential(r_cut=r_cut, + code=patch_code, + param_array=[]) mc.pair_potential = patches # compute the density and pressure @@ -444,13 +482,13 @@ def run_npt_sim(job, device, complete_filename): 'mode': 'nvt', 'device_name': 'cpu', 'ranks_per_partition': NUM_CPU_RANKS, - 'aggregator': partition_jobs_cpu_mpi + 'aggregator': partition_jobs_cpu_mpi_nvt }, { 'mode': 'npt', 'device_name': 'cpu', 'ranks_per_partition': NUM_CPU_RANKS, - 'aggregator': partition_jobs_cpu_mpi + 'aggregator': partition_jobs_cpu_mpi_npt }, ] @@ -483,6 +521,7 @@ def add_sampling_job(mode, device_name, ranks_per_partition, aggregator): @Project.pre.after(patchy_particle_pressure_create_initial_state) @Project.post.isfile(f'{mode}_{device_name}_complete') + #@Project.post(lambda j: j.isfile(f'{mode}_{device_name}_complete') or j.sp.pressure <= 0) @Project.operation(name=f'patchy_particle_pressure_{mode}_{device_name}', directives=directives, aggregator=aggregator) @@ -495,7 +534,8 @@ def sampling_operation(*jobs): job = jobs[communicator.partition] if communicator.rank == 0: - print(f'starting patchy_particle_pressure_{mode}_{device_name}:', job) + print(f'starting patchy_particle_pressure_{mode}_{device_name}:', + job) if device_name == 'gpu': device_cls = hoomd.device.GPU From 1de4ca54ca3ae2d58f9f521c7db49894ab7d9b7e Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 4 Oct 2023 12:49:43 -0400 Subject: [PATCH 25/42] Update plotting --- hoomd_validation/patchy_particle_pressure.py | 58 +++++++++++++++----- hoomd_validation/util.py | 41 ++++++++++---- 2 files changed, 74 insertions(+), 25 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 8ae3de72..9b5d99bb 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -605,7 +605,7 @@ def patchy_particle_pressure_analyze(job): # Plot results fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') - ax = fig.add_subplot(2, 1, 1) + ax = fig.add_subplot(2, 2, 1) util.plot_timeseries(ax=ax, timesteps=timesteps, data=densities, @@ -614,31 +614,54 @@ def patchy_particle_pressure_analyze(job): max_points=500) ax.legend() - ax = fig.add_subplot(2, 1, 2) + ax_distribution = fig.add_subplot(2, 2, 2, sharey=ax) + util.plot_distribution( + ax_distribution, + {k: v for k, v in densities.items() if not k.startswith('nvt')}, + #densities, + r'', + expected=job.sp.density, + bins=50, + plot_rotated=True, + ) + + ax = fig.add_subplot(2, 2, 3) util.plot_timeseries(ax=ax, timesteps=timesteps, data=pressures, ylabel=r"$\beta P$", expected=job.sp.pressure, - max_points=500) + max_points=500,) + ax_distribution = fig.add_subplot(2, 2, 4, sharey=ax) + util.plot_distribution( + ax_distribution, + pressures, + r'', + expected=job.sp.pressure, + bins=50, + plot_rotated=True, + ) fig.suptitle(f"$\\rho={job.statepoint.density}$, " f"$N={job.statepoint.num_particles}$, " + f"T={job.statepoint.temperature}, " + f"$\\chi={job.statepoint.chi}$, " f"replicate={job.statepoint.replicate_idx}") - fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight') + fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight', transparent=False) job.document['patchy_particle_pressure_analysis_complete'] = True -@Project.pre( - lambda *jobs: util.true_all(*jobs, key='patchy_particle_pressure_analysis_complete')) +@Project.pre(lambda *jobs: util.true_all( + *jobs, key='patchy_particle_pressure_analysis_complete')) @Project.post(lambda *jobs: util.true_all( *jobs, key='patchy_particle_pressure_compare_modes_complete')) -@Project.operation(directives=dict(executable=CONFIG["executable"]), - aggregator=aggregator.groupby( - key=['density', 'num_particles'], - sort_by='replicate_idx', - select=is_patchy_particle_pressure)) +@Project.operation( + directives=dict(executable=CONFIG["executable"]), + aggregator=aggregator.groupby( + key=['pressure', 'density', 'temperature', 'chi', 'num_particles'], + sort_by='replicate_idx', + select=is_patchy_particle_pressure)) def patchy_particle_pressure_compare_modes(*jobs): """Compares the tested simulation modes.""" import numpy @@ -667,12 +690,15 @@ def patchy_particle_pressure_compare_modes(*jobs): # grab the common statepoint parameters set_density = jobs[0].sp.density set_pressure = jobs[0].sp.pressure + set_temperature = jobs[0].sp.temperature + set_chi = jobs[0].sp.chi num_particles = jobs[0].sp.num_particles quantity_reference = dict(density=set_density, pressure=set_pressure) fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') - fig.suptitle(f"$\\rho={set_density}$, $N={num_particles}$") + fig.suptitle( + f"$\\rho={set_density}$, $N={num_particles}$, $T={set_temperature}$, $\\chi={set_chi}$") for i, quantity_name in enumerate(quantity_names): ax = fig.add_subplot(2, 1, i + 1) @@ -700,9 +726,13 @@ def patchy_particle_pressure_compare_modes(*jobs): relative_scale=1000, separate_nvt_npt=True) - filename = f'patchy_particle_pressure_compare_density{round(set_density, 2)}.svg' + filename = f'patchy_particle_pressure_compare_' + filename += f'density{round(set_density, 2)}_' + filename += f'temperature{round(set_temperature, 4)}_' + filename += f'pressure{round(set_pressure, 3)}_' + filename += f'chi{round(set_chi, 2)}.svg' fig.savefig(os.path.join(jobs[0]._project.path, filename), - bbox_inches='tight') + bbox_inches='tight', transparent=False) for job in jobs: job.document['patchy_particle_pressure_compare_modes_complete'] = True diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index 9385fc53..8fc92a52 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -152,7 +152,7 @@ def make_seed(job, sim_mode=None): return int(signac.job.calc_id(statepoint), 16) & 0xffff -def plot_distribution(ax, data, xlabel, expected=None, bins=100): +def plot_distribution(ax, data, xlabel, expected=None, bins=100, plot_rotated=False): """Plot distributions.""" import numpy @@ -168,30 +168,49 @@ def plot_distribution(ax, data, xlabel, expected=None, bins=100): bins=bins, range=(range_min, range_max), density=True) + bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) if numpy.all(data_arr == data_arr[0]): histogram[:] = 0 max_density_histogram = max(max_density_histogram, numpy.max(histogram)) - ax.plot(bin_edges[:-1], histogram, label=mode) + if plot_rotated: + XX, YY = histogram, bin_centers #bin_edges[:-1] + ax.set_ylabel(xlabel) + ax.set_xlabel('probability density') + else: + XX, YY = bin_edges[:-1], histogram + ax.set_xlabel(xlabel) + ax.set_ylabel('probability density') + ax.plot(XX, YY, label=mode) - ax.set_xlabel(xlabel) - ax.set_ylabel('probability density') if callable(expected): - ax.plot(bin_edges[:-1], - expected(bin_edges[:-1]), + if plot_rotated: + #Y, X = bin_edges[:-1], expected(bin_edges[:-1]) + Y, X = bin_centers, expected(bin_centers) + else: + X, Y = bin_edges[:-1], expected(bin_edges[:-1]) + ax.plot(X, + Y, linestyle='dashed', color='k', label='expected') elif expected is not None: - ax.vlines(x=expected, - ymin=0, - ymax=max_density_histogram, - linestyles='dashed', - colors='k') + if plot_rotated: + ax.hlines(y=expected, + xmin=0, + xmax=max_density_histogram, + linestyles='dashed', + colors='k') + else: + ax.vlines(x=expected, + ymin=0, + ymax=max_density_histogram, + linestyles='dashed', + colors='k') def plot_vs_expected(ax, From e00350dee1fa584ad8da3e03def869288ab88588 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 4 Oct 2023 12:50:43 -0400 Subject: [PATCH 26/42] Changes during testing --- dashboard.py | 15 +++++- hoomd_validation/patchy_particle_pressure.py | 50 +++++++++++++++----- 2 files changed, 52 insertions(+), 13 deletions(-) diff --git a/dashboard.py b/dashboard.py index ab309f7b..a403fd8b 100644 --- a/dashboard.py +++ b/dashboard.py @@ -32,11 +32,24 @@ def job_title(self, job): f"rho={job.statepoint.density}" elif (job.statepoint.subproject == 'hard_disk' or job.statepoint.subproject == 'hard_sphere' - or job.statepoint.subproject == 'simple_polygon'): + or job.statepoint.subproject == 'simple_polygon' + or job.statepoint.subproject == 'patchy_particle_pressure'): return f"{job.statepoint.subproject}: rho={job.statepoint.density}" else: raise RuntimeError("Unexpected job") + def job_sorter(self, job): + if job.statepoint.subproject == 'patchy_particle_pressure': + return ( + job.sp.density, + job.sp.pressure, + job.sp.temperature, + job.sp.chi, + job.sp.replicate_idx, + ) + else: + return job.statepoint.num_particles + if __name__ == "__main__": ValidationDashboard(modules=modules).main() diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 5548e076..2c1e5c6f 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -20,7 +20,7 @@ TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS WRITE_PERIOD = 1_000 -LOG_PERIOD = {'trajectory': 50_000, 'quantities': 100} +LOG_PERIOD = {'trajectory': 50_000, 'quantities': 500} NUM_CPU_RANKS = min(9, CONFIG["max_cores_sim"]) WALLTIME_STOP_SECONDS = CONFIG["max_walltime"] * 3600 - 10 * 60 @@ -28,12 +28,26 @@ def job_statepoints(): """list(dict): A list of statepoints for this subproject.""" - num_particles = 64**2 + num_particles = 8**3 replicate_indices = range(CONFIG["replicates"]) # statepoint chosen to be in a dense liquid - # nvt simulations at density = 0.6 yielded a measured pressure of - # -1.218 +/- 0.002 (mean +/- std. error of means) over 32 replicas - params_list = [(0.7, 0.6, -1.218, 0.7, 1.5)] # kT, rho, pressure, chi, lambda_ + # nvt simulations at density = 0.95 yielded a measured pressure of + # 6.415 +/- 0.003 (mean +/- std. error of means) over 32 replicas + params_list = [ + (10.0, 0.95, 11.83, 0.7, 1.5), + (1.0, 0.95, 8.27, 0.7, 1.5), + (0.6, 0.95, 8.27, 0.7, 1.5), + # next 3 are from 10.1063/1.3054361 + # pressure from NVT simulations + (0.5714, 0.8, -5.871736803897682, 1.0, 1.5), + (1.0, 0.8, -1.0235450460184001, 1.0, 1.5), + (3.0, 0.8, 3.7415391338219885, 1.0, 1.5), + # next 2 are from 10.1063/1.3054361 + # pressure from webplotdigitizer fig 7 + (1.0, 0.8, 1.70178459919393820, 1.0, 1.5), # pressure = pressure = pressure/kT + (3.0, 0.8, 13.549010655204555, 1.0, 1.5), # pressure = pressure + (3.0, 0.8, 4.516336885068185, 1.0, 1.5), # pressure = pressure/kT + ] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: yield ({ @@ -52,17 +66,27 @@ def is_patchy_particle_pressure(job): """Test if a job is part of the patchy_particle_pressure subproject.""" return job.statepoint['subproject'] == 'patchy_particle_pressure' +def is_patchy_particle_pressure_positive_pressure(job): + """Test if a job is part of the patchy_particle_pressure subproject.""" + return job.statepoint['subproject'] == 'patchy_particle_pressure' and job.statepoint['pressure'] > 0.0 + + partition_jobs_cpu_serial = aggregator.groupsof(num=min( CONFIG["replicates"], CONFIG["max_cores_submission"]), sort_by='density', select=is_patchy_particle_pressure,) -partition_jobs_cpu_mpi = aggregator.groupsof(num=min( +partition_jobs_cpu_mpi_nvt = aggregator.groupsof(num=min( CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), sort_by='density', select=is_patchy_particle_pressure,) +partition_jobs_cpu_mpi_npt = aggregator.groupsof(num=min( + CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure_positive_pressure,) + partition_jobs_gpu = aggregator.groupsof(num=min(CONFIG["replicates"], CONFIG["max_gpus_submission"]), sort_by='density', select=is_patchy_particle_pressure,) @@ -73,7 +97,7 @@ def is_patchy_particle_pressure(job): executable=CONFIG["executable"], nranks=util.total_ranks_function(NUM_CPU_RANKS), walltime=1), - aggregator=partition_jobs_cpu_mpi,) + aggregator=partition_jobs_cpu_mpi_nvt,) def patchy_particle_pressure_create_initial_state(*jobs): """Create initial system configuration.""" import hoomd @@ -445,13 +469,13 @@ def run_npt_sim(job, device, complete_filename): 'mode': 'nvt', 'device_name': 'cpu', 'ranks_per_partition': NUM_CPU_RANKS, - 'aggregator': partition_jobs_cpu_mpi + 'aggregator': partition_jobs_cpu_mpi_nvt }, { 'mode': 'npt', 'device_name': 'cpu', 'ranks_per_partition': NUM_CPU_RANKS, - 'aggregator': partition_jobs_cpu_mpi + 'aggregator': partition_jobs_cpu_mpi_npt }, ] @@ -585,6 +609,7 @@ def patchy_particle_pressure_analyze(job): fig.suptitle(f"$\\rho={job.statepoint.density}$, " f"$N={job.statepoint.num_particles}$, " + f"temperature={job.statepoint.temperature}, " f"replicate={job.statepoint.replicate_idx}") fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight') @@ -597,7 +622,7 @@ def patchy_particle_pressure_analyze(job): *jobs, key='patchy_particle_pressure_compare_modes_complete')) @Project.operation(directives=dict(executable=CONFIG["executable"]), aggregator=aggregator.groupby( - key=['density', 'num_particles'], + key=['pressure', 'density', 'temperature', 'chi', 'num_particles'], sort_by='replicate_idx', select=is_patchy_particle_pressure)) def patchy_particle_pressure_compare_modes(*jobs): @@ -628,12 +653,13 @@ def patchy_particle_pressure_compare_modes(*jobs): # grab the common statepoint parameters set_density = jobs[0].sp.density set_pressure = jobs[0].sp.pressure + set_temperature = jobs[0].sp.temperature num_particles = jobs[0].sp.num_particles quantity_reference = dict(density=set_density, pressure=set_pressure) fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') - fig.suptitle(f"$\\rho={set_density}$, $N={num_particles}$") + fig.suptitle(f"$\\rho={set_density}$, $N={num_particles}$, $T={set_temperature}$") for i, quantity_name in enumerate(quantity_names): ax = fig.add_subplot(2, 1, i + 1) @@ -661,7 +687,7 @@ def patchy_particle_pressure_compare_modes(*jobs): relative_scale=1000, separate_nvt_npt=True) - filename = f'patchy_particle_pressure_compare_density{round(set_density, 2)}.svg' + filename = f'patchy_particle_pressure_compare_density{round(set_density, 2)}_temperature{round(set_temperature, 4)}.svg' fig.savefig(os.path.join(jobs[0]._project.path, filename), bbox_inches='tight') From 28d5a8e98af62110d55ab2e09cc591ebc5261265 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 4 Oct 2023 17:54:56 -0400 Subject: [PATCH 27/42] Update statepoints based on nvt simulations --- hoomd_validation/patchy_particle_pressure.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index f595ea3e..d208f593 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -34,20 +34,14 @@ def job_statepoints(): # nvt simulations at density = 0.95 yielded a measured pressure of # 6.415 +/- 0.003 (mean +/- std. error of means) over 32 replicas params_list = [ - (10.0, 0.95, 11.83, 0.7, 1.5), - (1.0, 0.95, 8.27, 0.7, 1.5), - (0.6, 0.95, 8.27, 0.7, 1.5), + (10.0, 0.95, 11.942862877236088, 0.7, 1.5), + (1.0, 0.95, 10.147994162256277, 0.7, 1.5), + (0.6, 0.95, 8.926992055668622, 0.7, 1.5), # next 3 are from 10.1063/1.3054361 # pressure from NVT simulations - (0.5714, 0.8, -5.871736803897682, 1.0, 1.5), - (1.0, 0.8, -1.0235450460184001, 1.0, 1.5), - (3.0, 0.8, 3.7415391338219885, 1.0, 1.5), - # next 2 are from 10.1063/1.3054361 - # pressure from webplotdigitizer fig 7 - (1.0, 0.8, 1.70178459919393820, 1.0, 1.5 - ), # pressure = pressure = pressure/kT - (3.0, 0.8, 13.549010655204555, 1.0, 1.5), # pressure = pressure - (3.0, 0.8, 4.516336885068185, 1.0, 1.5), # pressure = pressure/kT + (0.5714, 0.8, -0.22167766330477995, 1.0, 1.5), + (1.0, 0.8, 2.2911186181, 1.0, 1.5), + (3.0, 0.8, 4.8393938019, 1.0, 1.5), ] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: From 760c3395a25ad0b9264ff7ca5f74f482abf8d483 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Sun, 8 Oct 2023 18:38:44 -0400 Subject: [PATCH 28/42] Allow floats for short job walltime in config --- hoomd_validation/config_parser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd_validation/config_parser.py b/hoomd_validation/config_parser.py index 47ba41fe..d5f27777 100644 --- a/hoomd_validation/config_parser.py +++ b/hoomd_validation/config_parser.py @@ -34,7 +34,7 @@ def __init__(self, config_file_path=DEFAULT_CONFIG_PATH): config.get("max_cores_submission", 16)) self["max_gpus_submission"] = int(config.get("max_gpus_submission", 1)) self["max_walltime"] = int(config.get("max_walltime", 24)) - self["short_walltime"] = int(config.get("short_walltime", 2)) + self["short_walltime"] = float(config.get("short_walltime", 2)) self["replicates"] = int(config.get("replicates", 32)) self["enable_llvm"] = bool(config.get("enable_llvm", True)) self["enable_gpu"] = bool(config.get("enable_gpu", True)) From f0e6510a0d796776e136d2f20c8b4a943a174f1b Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Sun, 8 Oct 2023 18:39:30 -0400 Subject: [PATCH 29/42] Update patchy particle pressure. --- hoomd_validation/patchy_particle_pressure.py | 22 +++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index d208f593..4b85f592 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -14,34 +14,34 @@ # Run parameters shared between simulations. # Step counts must be even and a multiple of the log quantity period. RANDOMIZE_STEPS = 20_000 -EQUILIBRATE_STEPS = 100_000 -RUN_STEPS = 1_000_000 +EQUILIBRATE_STEPS = 200_000 +RUN_STEPS = 2_000_000 RESTART_STEPS = RUN_STEPS // 50 TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS WRITE_PERIOD = 1_000 LOG_PERIOD = {'trajectory': 50_000, 'quantities': 500} -NUM_CPU_RANKS = min(9, CONFIG["max_cores_sim"]) +NUM_CPU_RANKS = min(16, CONFIG["max_cores_sim"]) WALLTIME_STOP_SECONDS = CONFIG["max_walltime"] * 3600 - 10 * 60 def job_statepoints(): """list(dict): A list of statepoints for this subproject.""" - num_particles = 8**3 + num_particles = 16**3 replicate_indices = range(CONFIG["replicates"]) # statepoint chosen to be in a dense liquid # nvt simulations at density = 0.95 yielded a measured pressure of # 6.415 +/- 0.003 (mean +/- std. error of means) over 32 replicas params_list = [ - (10.0, 0.95, 11.942862877236088, 0.7, 1.5), - (1.0, 0.95, 10.147994162256277, 0.7, 1.5), - (0.6, 0.95, 8.926992055668622, 0.7, 1.5), + (10.0, 0.95, 11.951630531873338, 0.7, 1.5), + (1.0, 0.95, 10.208625410213362, 0.7, 1.5), + (0.6, 0.95, 8.927827449359, 0.7, 1.5), # next 3 are from 10.1063/1.3054361 # pressure from NVT simulations (0.5714, 0.8, -0.22167766330477995, 1.0, 1.5), - (1.0, 0.8, 2.2911186181, 1.0, 1.5), - (3.0, 0.8, 4.8393938019, 1.0, 1.5), + (1.0, 0.8, 2.2766339608381325, 1.0, 1.5), + (3.0, 0.8, 4.837436833423719, 1.0, 1.5), ] # kT, rho, pressure, chi, lambda_ for temperature, density, pressure, chi, lambda_ in params_list: for idx in replicate_indices: @@ -710,12 +710,14 @@ def patchy_particle_pressure_compare_modes(*jobs): expected=reference, relative_scale=1000, separate_nvt_npt=True) + ax.axhline(0.0, c='k', ls='--') filename = f'patchy_particle_pressure_compare_' filename += f'density{round(set_density, 2)}_' filename += f'temperature{round(set_temperature, 4)}_' filename += f'pressure{round(set_pressure, 3)}_' - filename += f'chi{round(set_chi, 2)}.svg' + filename += f'chi{round(set_chi, 2)}_' + filename += f'{num_particles}particles.svg' fig.savefig(os.path.join(jobs[0]._project.path, filename), bbox_inches='tight', transparent=False) From 308755e0be24fd128306412b06c0ea43416a3d80 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Mon, 9 Oct 2023 13:57:07 -0400 Subject: [PATCH 30/42] Parse max_walltime as a float --- hoomd_validation/config_parser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd_validation/config_parser.py b/hoomd_validation/config_parser.py index d5f27777..9639a579 100644 --- a/hoomd_validation/config_parser.py +++ b/hoomd_validation/config_parser.py @@ -33,7 +33,7 @@ def __init__(self, config_file_path=DEFAULT_CONFIG_PATH): self["max_cores_submission"] = int( config.get("max_cores_submission", 16)) self["max_gpus_submission"] = int(config.get("max_gpus_submission", 1)) - self["max_walltime"] = int(config.get("max_walltime", 24)) + self["max_walltime"] = float(config.get("max_walltime", 24)) self["short_walltime"] = float(config.get("short_walltime", 2)) self["replicates"] = int(config.get("replicates", 32)) self["enable_llvm"] = bool(config.get("enable_llvm", True)) From acb65fad779cfc05eedb85312f7799e9e1b7f45c Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Mon, 9 Oct 2023 14:48:54 -0400 Subject: [PATCH 31/42] Address PR comments --- hoomd_validation/init.py | 14 +- hoomd_validation/patchy_particle_pressure.py | 209 ++++++++++++++----- hoomd_validation/util.py | 80 ++----- 3 files changed, 179 insertions(+), 124 deletions(-) diff --git a/hoomd_validation/init.py b/hoomd_validation/init.py index c7bd2774..0fe01589 100644 --- a/hoomd_validation/init.py +++ b/hoomd_validation/init.py @@ -17,13 +17,13 @@ import patchy_particle_pressure subprojects = [ -# alj_2d, -# lj_fluid, -# lj_union, -# hard_disk, -# hard_sphere, -# simple_polygon, - patchy_particle_pressure + # alj_2d, + # lj_fluid, + # lj_union, + # hard_disk, + # hard_sphere, + # simple_polygon, + patchy_particle_pressure, ] project = signac.init_project(path=config.project_root) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 4b85f592..6304249f 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -30,20 +30,24 @@ def job_statepoints(): """list(dict): A list of statepoints for this subproject.""" num_particles = 16**3 replicate_indices = range(CONFIG["replicates"]) - # statepoint chosen to be in a dense liquid - # nvt simulations at density = 0.95 yielded a measured pressure of - # 6.415 +/- 0.003 (mean +/- std. error of means) over 32 replicas params_list = [ - (10.0, 0.95, 11.951630531873338, 0.7, 1.5), - (1.0, 0.95, 10.208625410213362, 0.7, 1.5), - (0.6, 0.95, 8.927827449359, 0.7, 1.5), - # next 3 are from 10.1063/1.3054361 - # pressure from NVT simulations - (0.5714, 0.8, -0.22167766330477995, 1.0, 1.5), - (1.0, 0.8, 2.2766339608381325, 1.0, 1.5), - (3.0, 0.8, 4.837436833423719, 1.0, 1.5), - ] # kT, rho, pressure, chi, lambda_ - for temperature, density, pressure, chi, lambda_ in params_list: + # kern-frenkel statepoints. densities/temperatures chosen to prorduce + # a dense liquid. pressures measure from NVT simulations. + (10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), + (1.0, 0.95, 10.208625410213362, 0.7, 1.5, 1.0), + (0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), + # hard sphere + square well statepoints, from 10.1063/1.3054361 + # pressure from NVT simulations, NOT from values in paper + (0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), + (1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), + (3.0, 0.8, 4.837436833423719, 1.0, 1.5, 1.0), + # hard sphere + square well + repulsive shoulder statepoints. + # temperatures/densities from initial tests, pressures from NVT + # simulations. + (3.0, 0.7, -5.0, 1.0, 1.0, -1.0), + (1.0, 0.7, -5.0, 1.0, 1.0, -1.0), + ] # kT, rho, pressure, chi, lambda_, long_range_interaction_scale_factor + for temperature, density, pressure, chi, lambda_, lrisf in params_list: for idx in replicate_indices: yield ({ "subproject": "patchy_particle_pressure", @@ -54,6 +58,7 @@ def job_statepoints(): "num_particles": num_particles, "replicate_idx": idx, "temperature": temperature, + "long_range_interaction_scale_factor": lrisf, }) @@ -61,11 +66,6 @@ def is_patchy_particle_pressure(job): """Test if a job is part of the patchy_particle_pressure subproject.""" return job.statepoint['subproject'] == 'patchy_particle_pressure' -def is_patchy_particle_pressure_positive_pressure(job): - """Test if a job is part of the patchy_particle_pressure subproject.""" - return job.statepoint['subproject'] == 'patchy_particle_pressure' and job.statepoint['pressure'] > 0.0 - - def is_patchy_particle_pressure_positive_pressure(job): """Test if a job is part of the patchy_particle_pressure subproject.""" @@ -73,27 +73,94 @@ def is_patchy_particle_pressure_positive_pressure(job): 'subproject'] == 'patchy_particle_pressure' and job.statepoint[ 'pressure'] > 0.0 -partition_jobs_cpu_mpi_nvt = aggregator.groupsof(num=min( - CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), - sort_by='density', - select=is_patchy_particle_pressure,) -partition_jobs_cpu_mpi_npt = aggregator.groupsof(num=min( - CONFIG["replicates"], CONFIG["max_cores_submission"] // NUM_CPU_RANKS), - sort_by='density', - select=is_patchy_particle_pressure_positive_pressure,) +partition_jobs_cpu_mpi_nvt = aggregator.groupsof( + num=min(CONFIG["replicates"], + CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure, +) + +partition_jobs_cpu_mpi_npt = aggregator.groupsof( + num=min(CONFIG["replicates"], + CONFIG["max_cores_submission"] // NUM_CPU_RANKS), + sort_by='density', + select=is_patchy_particle_pressure_positive_pressure, +) + +partition_jobs_gpu = aggregator.groupsof( + num=min(CONFIG["replicates"], CONFIG["max_gpus_submission"]), + sort_by='density', + select=is_patchy_particle_pressure, +) + + +def _single_patch_kern_frenkel_code(delta_rad, sq_well_lambda, sigma, kT, + long_range_interaction_scale_factor): + """Generate code for JIT compilation of Kern-Frenkel potential. + + Args: + delta_rad (float): Half-opening angle of patchy interaction in radians + + sq_well_lambda (float): range of patchy interaction relative to hard + core diameter + + sigma (float): Diameter of hard core of particles + + kT (float): Temperature; sets the energy scale -partition_jobs_gpu = aggregator.groupsof(num=min(CONFIG["replicates"], - CONFIG["max_gpus_submission"]), - sort_by='density', select=is_patchy_particle_pressure,) + long_range_interaction_scale_factor (float): factor to multiply the + pair potential by when the interparticle separation is between + sq_well_lambda/2*sigma and sq_well_lambda*sigma. + + The terminology (e.g., `ehat`) comes from the "Modelling Patchy Particles" + HOOMD-blue tutorial. + + """ + patch_code = f""" + const float delta = {delta_rad}; + const float lambda = {sq_well_lambda:f}; + const float sigma = {sigma:f}; // hard core diameter + const float kT = {kT:f}; + const float long_range_interaction_scale_factor = { + long_range_interaction_scale_factor:f}; + const vec3 ehat_particle_reference_frame(1, 0, 0); + vec3 ehat_i = rotate(q_i, ehat_particle_reference_frame); + vec3 ehat_j = rotate(q_j, ehat_particle_reference_frame); + vec3 r_hat_ij = r_ij / sqrtf(dot(r_ij, r_ij)); + bool patch_on_i_is_aligned_with_r_ij = dot(ehat_i, r_hat_ij) >= cos(delta); + bool patch_on_j_is_aligned_with_r_ji = dot(ehat_j, -r_hat_ij) >= cos(delta); + float rsq = dot(r_ij, r_ij); + if (patch_on_i_is_aligned_with_r_ij && patch_on_j_is_aligned_with_r_ji) + {{ + if (rsq < lambda / 2 * sigma * lambda / 2 * sigma) + {{ + return -1 / kT; + }} + else if (rsq < lambda * sigma * lambda * sigma) + {{ + return -1 / kT * long_range_interaction_scale_factor; + }} + else + {{ + return 0.0; + }} + }} + else + {{ + return 0.0; + }} + """ + return patch_code @Project.post.isfile('patchy_particle_pressure_initial_state.gsd') -@Project.operation(directives=dict( - executable=CONFIG["executable"], - nranks=util.total_ranks_function(NUM_CPU_RANKS), - walltime=1), - aggregator=partition_jobs_cpu_mpi_nvt,) +@Project.operation( + directives=dict(executable=CONFIG["executable"], + nranks=util.total_ranks_function(NUM_CPU_RANKS), + walltime=1), + aggregator=partition_jobs_cpu_mpi_nvt, +) def patchy_particle_pressure_create_initial_state(*jobs): """Create initial system configuration.""" import hoomd @@ -112,6 +179,8 @@ def patchy_particle_pressure_create_initial_state(*jobs): temperature = job.statepoint['temperature'] chi = job.statepoint['chi'] lambda_ = job.statepoint['lambda_'] + long_range_interaction_scale_factor = job.statepoint[ + 'long_range_interaction_scale_factor'] box_volume = num_particles / density L = box_volume**(1 / 3.) @@ -143,8 +212,13 @@ def patchy_particle_pressure_create_initial_state(*jobs): diameter = 1.0 mc.shape['A'] = dict(diameter=diameter, orientable=True) delta = 2 * numpy.arcsin(numpy.sqrt(chi)) - patch_code = util._single_patch_kern_frenkel_code(delta, lambda_, diameter, - temperature) + patch_code = _single_patch_kern_frenkel_code( + delta, + lambda_, + diameter, + temperature, + long_range_interaction_scale_factor, + ) r_cut = diameter + diameter * (lambda_ - 1) patches = hoomd.hpmc.pair.user.CPPPotential(r_cut=r_cut, code=patch_code, @@ -204,12 +278,19 @@ def make_mc_simulation(job, temperature = job.statepoint['temperature'] chi = job.statepoint['chi'] lambda_ = job.statepoint['lambda_'] + long_range_interaction_scale_factor = job.statepoint[ + 'long_range_interaction_scale_factor'] diameter = 1.0 mc = hoomd.hpmc.integrate.Sphere(default_d=0.05, default_a=0.1) mc.shape['A'] = dict(diameter=diameter, orientable=True) delta = 2 * numpy.arcsin(numpy.sqrt(chi)) - patch_code = util._single_patch_kern_frenkel_code(delta, lambda_, diameter, - temperature) + patch_code = _single_patch_kern_frenkel_code( + delta, + lambda_, + diameter, + temperature, + long_range_interaction_scale_factor, + ) r_cut = diameter + diameter * (lambda_ - 1) patches = hoomd.hpmc.pair.user.CPPPotential(r_cut=r_cut, code=patch_code, @@ -477,7 +558,6 @@ def run_npt_sim(job, device, complete_filename): }, ] -#job_definitions = [] if CONFIG["enable_gpu"]: job_definitions.extend([ { @@ -506,7 +586,6 @@ def add_sampling_job(mode, device_name, ranks_per_partition, aggregator): @Project.pre.after(patchy_particle_pressure_create_initial_state) @Project.post.isfile(f'{mode}_{device_name}_complete') - #@Project.post(lambda j: j.isfile(f'{mode}_{device_name}_complete') or j.sp.pressure <= 0) @Project.operation(name=f'patchy_particle_pressure_{mode}_{device_name}', directives=directives, aggregator=aggregator) @@ -603,7 +682,6 @@ def patchy_particle_pressure_analyze(job): util.plot_distribution( ax_distribution, {k: v for k, v in densities.items() if not k.startswith('nvt')}, - #densities, r'', expected=job.sp.density, bins=50, @@ -611,12 +689,14 @@ def patchy_particle_pressure_analyze(job): ) ax = fig.add_subplot(2, 2, 3) - util.plot_timeseries(ax=ax, - timesteps=timesteps, - data=pressures, - ylabel=r"$\beta P$", - expected=job.sp.pressure, - max_points=500,) + util.plot_timeseries( + ax=ax, + timesteps=timesteps, + data=pressures, + ylabel=r"$\beta P$", + expected=job.sp.pressure, + max_points=500, + ) ax_distribution = fig.add_subplot(2, 2, 4, sharey=ax) util.plot_distribution( ax_distribution, @@ -632,7 +712,9 @@ def patchy_particle_pressure_analyze(job): f"T={job.statepoint.temperature}, " f"$\\chi={job.statepoint.chi}$, " f"replicate={job.statepoint.replicate_idx}") - fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight', transparent=False) + fig.savefig(job.fn('nvt_npt_plots.svg'), + bbox_inches='tight', + transparent=False) job.document['patchy_particle_pressure_analysis_complete'] = True @@ -641,12 +723,18 @@ def patchy_particle_pressure_analyze(job): *jobs, key='patchy_particle_pressure_analysis_complete')) @Project.post(lambda *jobs: util.true_all( *jobs, key='patchy_particle_pressure_compare_modes_complete')) -@Project.operation( - directives=dict(executable=CONFIG["executable"]), - aggregator=aggregator.groupby( - key=['pressure', 'density', 'temperature', 'chi', 'num_particles'], - sort_by='replicate_idx', - select=is_patchy_particle_pressure)) +@Project.operation(directives=dict(executable=CONFIG["executable"]), + aggregator=aggregator.groupby( + key=[ + 'pressure', + 'density', + 'temperature', + 'chi', + 'num_particles', + 'long_range_interaction_scale_factor', + ], + sort_by='replicate_idx', + select=is_patchy_particle_pressure)) def patchy_particle_pressure_compare_modes(*jobs): """Compares the tested simulation modes.""" import numpy @@ -678,12 +766,16 @@ def patchy_particle_pressure_compare_modes(*jobs): set_temperature = jobs[0].sp.temperature set_chi = jobs[0].sp.chi num_particles = jobs[0].sp.num_particles + lrisf = jobs[0].sp.long_range_interaction_scale_factor quantity_reference = dict(density=set_density, pressure=set_pressure) fig = matplotlib.figure.Figure(figsize=(10, 10 / 1.618 * 2), layout='tight') - fig.suptitle( - f"$\\rho={set_density}$, $N={num_particles}$, $T={set_temperature}$, $\\chi={set_chi}$") + fig.suptitle(f"$\\rho={set_density}$, " + f"$N={num_particles}$, " + f"T={set_temperature}, " + f"$\\chi={set_chi}$, " + f"$\\varepsilon_2/\\varepsilon_1 = {lrisf}") for i, quantity_name in enumerate(quantity_names): ax = fig.add_subplot(2, 1, i + 1) @@ -712,14 +804,15 @@ def patchy_particle_pressure_compare_modes(*jobs): separate_nvt_npt=True) ax.axhline(0.0, c='k', ls='--') - filename = f'patchy_particle_pressure_compare_' + filename = 'patchy_particle_pressure_compare_' filename += f'density{round(set_density, 2)}_' filename += f'temperature{round(set_temperature, 4)}_' filename += f'pressure{round(set_pressure, 3)}_' filename += f'chi{round(set_chi, 2)}_' filename += f'{num_particles}particles.svg' fig.savefig(os.path.join(jobs[0]._project.path, filename), - bbox_inches='tight', transparent=False) + bbox_inches='tight', + transparent=False) for job in jobs: job.document['patchy_particle_pressure_compare_modes_complete'] = True diff --git a/hoomd_validation/util.py b/hoomd_validation/util.py index 8fc92a52..8c6676b8 100644 --- a/hoomd_validation/util.py +++ b/hoomd_validation/util.py @@ -28,7 +28,7 @@ def get_job_filename(sim_mode, device, name, type): return f"{sim_mode}_{suffix}_{name}.{type}" -def run_up_to_walltime(sim, end_step, steps, walltime_stop, minutes_per_run=None): +def run_up_to_walltime(sim, end_step, steps, walltime_stop): """Run a simulation, stopping early if a walltime limit is reached. Args: @@ -152,7 +152,12 @@ def make_seed(job, sim_mode=None): return int(signac.job.calc_id(statepoint), 16) & 0xffff -def plot_distribution(ax, data, xlabel, expected=None, bins=100, plot_rotated=False): +def plot_distribution(ax, + data, + independent_variable_label, + expected=None, + bins=100, + plot_rotated=False): """Plot distributions.""" import numpy @@ -176,27 +181,27 @@ def plot_distribution(ax, data, xlabel, expected=None, bins=100, plot_rotated=Fa max_density_histogram = max(max_density_histogram, numpy.max(histogram)) if plot_rotated: - XX, YY = histogram, bin_centers #bin_edges[:-1] - ax.set_ylabel(xlabel) + ax.set_ylabel(independent_variable_label) ax.set_xlabel('probability density') + ax.plot(histogram, bin_centers, label=mode) else: - XX, YY = bin_edges[:-1], histogram - ax.set_xlabel(xlabel) + ax.set_xlabel(independent_variable_label) ax.set_ylabel('probability density') - ax.plot(XX, YY, label=mode) - + ax.plot(bin_centers, histogram, label=mode) if callable(expected): if plot_rotated: - #Y, X = bin_edges[:-1], expected(bin_edges[:-1]) - Y, X = bin_centers, expected(bin_centers) + ax.plot(expected(bin_centers), + bin_centers, + linestyle='dashed', + color='k', + label='expected') else: - X, Y = bin_edges[:-1], expected(bin_edges[:-1]) - ax.plot(X, - Y, - linestyle='dashed', - color='k', - label='expected') + ax.plot(bin_centers, + expected(bin_centers), + linestyle='dashed', + color='k', + label='expected') elif expected is not None: if plot_rotated: @@ -318,46 +323,3 @@ 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 _single_patch_kern_frenkel_code(delta_rad, lambda_, sigma, kT): - """Generate code for JIT compilation of Kern-Frenkel potential. - - Args: - delta_rad (float): Half-opening angle of patchy interaction in radians - - lambda_ (float): range of patchy interaction relative to hard core - diameter - - sigma (float): Diameter of hard core of particles - - kT (float): Temperature; sets the energy scale - - The terminology (e.g., `ehat`) comes from the "Modelling Patchy Particles" - HOOMD-blue tutorial. - - """ - patch_code = f""" - const float delta = {delta_rad}; - const float lambda = {lambda_:f}; - const float sigma = {sigma:f}; // hard core diameter - const float kT = {kT:f}; - const vec3 ehat_particle_reference_frame(1, 0, 0); - vec3 ehat_i = rotate(q_i, ehat_particle_reference_frame); - vec3 ehat_j = rotate(q_j, ehat_particle_reference_frame); - vec3 r_hat_ij = r_ij / sqrtf(dot(r_ij, r_ij)); - bool patch_on_i_is_aligned_with_r_ij = dot(ehat_i, r_hat_ij) >= cos(delta); - bool patch_on_j_is_aligned_with_r_ji = dot(ehat_j, -r_hat_ij) >= cos(delta); - float rsq = dot(r_ij, r_ij); - if (patch_on_i_is_aligned_with_r_ij - && patch_on_j_is_aligned_with_r_ji - && dot(r_ij, r_ij) < lambda*sigma*lambda*sigma) - {{ - return -1 / kT; - }} - else - {{ - return 0.0; - }} - """ - return patch_code From 065a5b2654abdd275810b15f087e0dcf9729d69c Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 10 Oct 2023 09:02:57 -0400 Subject: [PATCH 32/42] Update square well + repulsive shoulder pressure --- hoomd_validation/patchy_particle_pressure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 6304249f..faabd559 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -44,7 +44,7 @@ def job_statepoints(): # hard sphere + square well + repulsive shoulder statepoints. # temperatures/densities from initial tests, pressures from NVT # simulations. - (3.0, 0.7, -5.0, 1.0, 1.0, -1.0), + (3.0, 0.7, 4.00804, 1.0, 1.0, -1.0), (1.0, 0.7, -5.0, 1.0, 1.0, -1.0), ] # kT, rho, pressure, chi, lambda_, long_range_interaction_scale_factor for temperature, density, pressure, chi, lambda_, lrisf in params_list: From d12c1dc88b825d18b9deb40df1c0cd2c39fda27b Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 10 Oct 2023 09:13:34 -0400 Subject: [PATCH 33/42] Update filenames --- hoomd_validation/patchy_particle_pressure.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index faabd559..d573fee6 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -809,7 +809,8 @@ def patchy_particle_pressure_compare_modes(*jobs): filename += f'temperature{round(set_temperature, 4)}_' filename += f'pressure{round(set_pressure, 3)}_' filename += f'chi{round(set_chi, 2)}_' - filename += f'{num_particles}particles.svg' + filename += f'{num_particles}particles_' + filename += f"$epsilonrepulsive{lrisf}.svg" fig.savefig(os.path.join(jobs[0]._project.path, filename), bbox_inches='tight', transparent=False) From 4d0022ca53a9fe017d081fcd457a3560c1dbb60e Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 10 Oct 2023 09:17:01 -0400 Subject: [PATCH 34/42] Update filenames --- hoomd_validation/patchy_particle_pressure.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index d573fee6..374af7b7 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -810,7 +810,7 @@ def patchy_particle_pressure_compare_modes(*jobs): filename += f'pressure{round(set_pressure, 3)}_' filename += f'chi{round(set_chi, 2)}_' filename += f'{num_particles}particles_' - filename += f"$epsilonrepulsive{lrisf}.svg" + filename += f"epsilonrepulsive{lrisf}.svg" fig.savefig(os.path.join(jobs[0]._project.path, filename), bbox_inches='tight', transparent=False) From 7d2fe5f9ff542eba32d91069644e4b3c4ea7c204 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 10 Oct 2023 09:26:21 -0400 Subject: [PATCH 35/42] Fix figure title --- hoomd_validation/patchy_particle_pressure.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 374af7b7..f64879fb 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -775,7 +775,8 @@ def patchy_particle_pressure_compare_modes(*jobs): f"$N={num_particles}$, " f"T={set_temperature}, " f"$\\chi={set_chi}$, " - f"$\\varepsilon_2/\\varepsilon_1 = {lrisf}") + "$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" + f"$={lrisf}$") for i, quantity_name in enumerate(quantity_names): ax = fig.add_subplot(2, 1, i + 1) From 133e718a62905cd7208491103a2979a875821509 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Tue, 10 Oct 2023 11:06:15 -0400 Subject: [PATCH 36/42] Fix figure titles --- hoomd_validation/patchy_particle_pressure.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index f64879fb..52eab704 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -707,11 +707,13 @@ def patchy_particle_pressure_analyze(job): plot_rotated=True, ) - fig.suptitle(f"$\\rho={job.statepoint.density}$, " - f"$N={job.statepoint.num_particles}$, " - f"T={job.statepoint.temperature}, " - f"$\\chi={job.statepoint.chi}$, " - f"replicate={job.statepoint.replicate_idx}") + fig.suptitle(f"$\\rho={job.sp.density}$, " + f"$N={job.sp.num_particles}$, " + f"T={job.sp.temperature}, " + f"$\\chi={job.sp.chi}$, " + f"replicate={job.statepoint.replicate_idx}, " + "$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" + f"$={job.sp.long_range_interaction_scale_factor}$") fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight', transparent=False) From 042f44356eb1cb73a487e2851db506338bfce1f7 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 19 Oct 2023 13:22:52 -0400 Subject: [PATCH 37/42] Update statepoints and run lengths --- hoomd_validation/patchy_particle_pressure.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 52eab704..fb0c9072 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -15,7 +15,7 @@ # Step counts must be even and a multiple of the log quantity period. RANDOMIZE_STEPS = 20_000 EQUILIBRATE_STEPS = 200_000 -RUN_STEPS = 2_000_000 +RUN_STEPS = 500_000 RESTART_STEPS = RUN_STEPS // 50 TOTAL_STEPS = RANDOMIZE_STEPS + EQUILIBRATE_STEPS + RUN_STEPS @@ -33,19 +33,19 @@ def job_statepoints(): params_list = [ # kern-frenkel statepoints. densities/temperatures chosen to prorduce # a dense liquid. pressures measure from NVT simulations. - (10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), + #(10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), (1.0, 0.95, 10.208625410213362, 0.7, 1.5, 1.0), - (0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), + #(0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), # hard sphere + square well statepoints, from 10.1063/1.3054361 # pressure from NVT simulations, NOT from values in paper - (0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), - (1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), + #(0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), + #(1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), (3.0, 0.8, 4.837436833423719, 1.0, 1.5, 1.0), # hard sphere + square well + repulsive shoulder statepoints. # temperatures/densities from initial tests, pressures from NVT # simulations. (3.0, 0.7, 4.00804, 1.0, 1.0, -1.0), - (1.0, 0.7, -5.0, 1.0, 1.0, -1.0), + #(1.0, 0.7, -5.0, 1.0, 1.0, -1.0), ] # kT, rho, pressure, chi, lambda_, long_range_interaction_scale_factor for temperature, density, pressure, chi, lambda_, lrisf in params_list: for idx in replicate_indices: @@ -559,6 +559,7 @@ def run_npt_sim(job, device, complete_filename): ] if CONFIG["enable_gpu"]: + job_definitions = [] job_definitions.extend([ { 'mode': 'nvt', From 67d01fe7ee19b15d5160814f63041fa7af0bad78 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 19 Oct 2023 16:22:03 -0400 Subject: [PATCH 38/42] Add docstring to job_sorter --- dashboard.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dashboard.py b/dashboard.py index a403fd8b..f6485981 100644 --- a/dashboard.py +++ b/dashboard.py @@ -39,6 +39,7 @@ def job_title(self, job): raise RuntimeError("Unexpected job") def job_sorter(self, job): + """Sort jobs.""" if job.statepoint.subproject == 'patchy_particle_pressure': return ( job.sp.density, From 4388fbb8fabd6e6a0d7a0498ae6e715a093af696 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 19 Oct 2023 16:23:16 -0400 Subject: [PATCH 39/42] Uncomment subprojects in init.py --- hoomd_validation/init.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/hoomd_validation/init.py b/hoomd_validation/init.py index 0fe01589..79ee0357 100644 --- a/hoomd_validation/init.py +++ b/hoomd_validation/init.py @@ -17,12 +17,12 @@ import patchy_particle_pressure subprojects = [ - # alj_2d, - # lj_fluid, - # lj_union, - # hard_disk, - # hard_sphere, - # simple_polygon, + alj_2d, + lj_fluid, + lj_union, + hard_disk, + hard_sphere, + simple_polygon, patchy_particle_pressure, ] From 5327c4771bbd3838616c1bbe6470d60d7c66f042 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 19 Oct 2023 16:26:26 -0400 Subject: [PATCH 40/42] flake8 --- hoomd_validation/patchy_particle_pressure.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index fb0c9072..30c142f7 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -33,19 +33,19 @@ def job_statepoints(): params_list = [ # kern-frenkel statepoints. densities/temperatures chosen to prorduce # a dense liquid. pressures measure from NVT simulations. - #(10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), + # (10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), (1.0, 0.95, 10.208625410213362, 0.7, 1.5, 1.0), - #(0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), + # (0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), # hard sphere + square well statepoints, from 10.1063/1.3054361 # pressure from NVT simulations, NOT from values in paper - #(0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), - #(1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), + # (0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), + # (1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), (3.0, 0.8, 4.837436833423719, 1.0, 1.5, 1.0), # hard sphere + square well + repulsive shoulder statepoints. # temperatures/densities from initial tests, pressures from NVT # simulations. (3.0, 0.7, 4.00804, 1.0, 1.0, -1.0), - #(1.0, 0.7, -5.0, 1.0, 1.0, -1.0), + # (1.0, 0.7, -5.0, 1.0, 1.0, -1.0), ] # kT, rho, pressure, chi, lambda_, long_range_interaction_scale_factor for temperature, density, pressure, chi, lambda_, lrisf in params_list: for idx in replicate_indices: @@ -713,7 +713,7 @@ def patchy_particle_pressure_analyze(job): f"T={job.sp.temperature}, " f"$\\chi={job.sp.chi}$, " f"replicate={job.statepoint.replicate_idx}, " - "$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" + r"$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" f"$={job.sp.long_range_interaction_scale_factor}$") fig.savefig(job.fn('nvt_npt_plots.svg'), bbox_inches='tight', @@ -778,7 +778,7 @@ def patchy_particle_pressure_compare_modes(*jobs): f"$N={num_particles}$, " f"T={set_temperature}, " f"$\\chi={set_chi}$, " - "$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" + r"$\\varepsilon_{\mathrm{rep}}/\\varepsilon_{\mathrm{att}}$" f"$={lrisf}$") for i, quantity_name in enumerate(quantity_names): From feb285841cdb0a4596cddd2739ec0de93d9341b7 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Thu, 19 Oct 2023 16:39:22 -0400 Subject: [PATCH 41/42] Remove calls to communicator.walltime on rank 0 --- hoomd_validation/patchy_particle_pressure.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 30c142f7..7a2549ee 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -245,8 +245,7 @@ def patchy_particle_pressure_create_initial_state(*jobs): ) if communicator.rank == 0: - print(f'completed patchy_particle_pressure_create_initial_state: ' - f'{job} in {communicator.walltime} s') + print(f'completed patchy_particle_pressure_create_initial_state: {job}') def make_mc_simulation(job, @@ -430,8 +429,7 @@ def run_nvt_sim(job, device, complete_filename): pathlib.Path(job.fn(complete_filename)).touch() device.notice('Done.') else: - device.notice('Ending run early due to walltime limits at:' - f'{device.communicator.walltime}') + device.notice(f'Ending {job} run early due to walltime limits.') def run_npt_sim(job, device, complete_filename): @@ -538,8 +536,7 @@ def run_npt_sim(job, device, complete_filename): pathlib.Path(job.fn(complete_filename)).touch() device.notice('Done.') else: - device.notice('Ending run early due to walltime limits at:' - f'{device.communicator.walltime}') + device.notice(f'Ending {job} run early due to walltime limits.') sampling_jobs = [] @@ -616,7 +613,7 @@ def sampling_operation(*jobs): if communicator.rank == 0: print(f'completed patchy_particle_pressure_{mode}_{device_name} ' - f'{job} in {communicator.walltime} s') + f'{job}') sampling_jobs.append(sampling_operation) From c4f8f6ab7a6d628fa632150a209950be26959709 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Sun, 22 Oct 2023 23:44:26 -0400 Subject: [PATCH 42/42] Remove unused state points. --- hoomd_validation/patchy_particle_pressure.py | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/hoomd_validation/patchy_particle_pressure.py b/hoomd_validation/patchy_particle_pressure.py index 7a2549ee..3a960c6c 100644 --- a/hoomd_validation/patchy_particle_pressure.py +++ b/hoomd_validation/patchy_particle_pressure.py @@ -31,21 +31,16 @@ def job_statepoints(): num_particles = 16**3 replicate_indices = range(CONFIG["replicates"]) params_list = [ - # kern-frenkel statepoints. densities/temperatures chosen to prorduce - # a dense liquid. pressures measure from NVT simulations. - # (10.0, 0.95, 11.951630531873338, 0.7, 1.5, 1.0), + # kern-frenkel. density and temperature chosen to produce + # a dense liquid. pressure measured from NVT simulations on CPU. (1.0, 0.95, 10.208625410213362, 0.7, 1.5, 1.0), - # (0.6, 0.95, 8.927827449359, 0.7, 1.5, 1.0), - # hard sphere + square well statepoints, from 10.1063/1.3054361 - # pressure from NVT simulations, NOT from values in paper - # (0.5714, 0.8, -0.2692376274894095, 1.0, 1.5, 1.0), - # (1.0, 0.8, 2.2766339608381325, 1.0, 1.5, 1.0), + # hard sphere + square well, from 10.1063/1.3054361 + # pressure from NVT simulations (on CPU), NOT from values in paper (3.0, 0.8, 4.837436833423719, 1.0, 1.5, 1.0), - # hard sphere + square well + repulsive shoulder statepoints. - # temperatures/densities from initial tests, pressures from NVT - # simulations. + # hard sphere + square well + repulsive shoulder. + # temperatures/densities for dense liquid based on initial tests. + # pressure measured from NVT simulations on CPU. (3.0, 0.7, 4.00804, 1.0, 1.0, -1.0), - # (1.0, 0.7, -5.0, 1.0, 1.0, -1.0), ] # kT, rho, pressure, chi, lambda_, long_range_interaction_scale_factor for temperature, density, pressure, chi, lambda_, lrisf in params_list: for idx in replicate_indices: