Patchy particle pressure #66

tcmoore3 · 2023-10-08T22:56:26Z

Description

Validate pressure computation for discontinuous pair potentials simulated with HPMC.

I also changed the comparison plots to show the distributions of sampled pressures and densities.

Motivation and context

The pressure pressure computation for systems with discontinuous pair potentials in HPMC is currently not tested for validity but should be.

I will add another validation test for the hard sphere + square well + square shoulder system we discussed, but I wanted to get this PR opened because the code will not change very much.

Resolves #43.

How has this been tested?

We discovered a conceptual error in the implementation of the SDF calculation for systems with pair potentials. We implemented a fix at glotzerlab/hoomd-blue@8d7cb2c. That commit is the version of HOOMD-blue I used to generate the results below.

I ran NVT simulations of 4,096 particles interacting either through a hard sphere + square well potential, or through a Kern–Frenkel potential, i.e., a hard sphere + square well model with an angular component to the pair potential. I ran NVT simulations, measured the pressure for the given density and temperature, and then ran NPT simulations at the corresponding pressure. The agreement between the pressure measured in NVT simulations and the density measured at the corresponding NPT simulations is shown in the figures below. These results are from 16 replicates at each statepoint. \chi = 1.0 corresponds to the hard sphere + square well state points, and \chi = 0.7 is the Kern–Frenkel model.

Checklist:

I have reviewed the Contributor Guidelines.
I agree with the terms of the HOOMD-blue Contributor Agreement.
My name is on the list of contributors (CONTRIBUTORS.md) in the pull request source branch.

…d-validation into patchy-particle-pressure

joaander · 2023-10-09T15:10:18Z

With short simulations, I tested:

sigma = 1.0
lambda_ = 1.5
epsilon = 1/3.0

with patch code:

    evaluate_square_well = f'''float rsq = dot(r_ij, r_ij);
                               if (rsq < {(lambda_ * sigma/2)**2}f)
                                   return param_array[0];
                               else if (rsq < {(lambda_ * sigma)**2}f)
                                   return param_array[1];
                               else
                                   return 0.0f;
                            '''
    square_well = hoomd.hpmc.pair.user.CPPPotential(r_cut=lambda_ * sigma,
                                                    code=evaluate_square_well,
                                                    param_array=[-epsilon, epsilon])

at a variety of pressures (2-10) and densities (around ~0.7).

joaander

I see you have GPU jobs definitions in the workflow. Do you plan to validate that these work properly before merging?

joaander · 2023-10-09T15:13:17Z

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))


Suggested change

self["max_walltime"] = int(config.get("max_walltime", 24))

self["max_walltime"] = float(config.get("max_walltime", 24))

For consistency with the below change.

joaander · 2023-10-09T15:13:57Z

hoomd_validation/init.py

+#    alj_2d,
+#    lj_fluid,
+#    lj_union,
+#    hard_disk,
+#    hard_sphere,
+#    simple_polygon,


Suggested change

# alj_2d,

# lj_fluid,

# lj_union,

# hard_disk,

# hard_sphere,

# simple_polygon,

alj_2d,

lj_fluid,

lj_union,

hard_disk,

hard_sphere,

simple_polygon,

joaander · 2023-10-09T15:14:09Z

hoomd_validation/init.py

+#    hard_disk,
+#    hard_sphere,
+#    simple_polygon,
+    patchy_particle_pressure


Suggested change

patchy_particle_pressure

patchy_particle_pressure,

joaander · 2023-10-09T15:14:56Z

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):


Suggested change

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):

unused.

joaander · 2023-10-09T15:16:09Z

hoomd_validation/util.py


        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]


Remove the need for the variables by placing the call to plot inside the conditional.

Alternately, replace meaningless XX and YY variable names with something meaningful.

joaander · 2023-10-09T15:19:51Z

hoomd_validation/util.py

+            ax.set_ylabel(xlabel)
+            ax.set_xlabel('probability density')
+        else:
+            XX, YY = bin_edges[:-1], histogram


The unrotated plot should also use bin_centers.

joaander · 2023-10-09T15:20:08Z

hoomd_validation/util.py

+            #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])


The unrotated plot should also use bin_centers.

joaander · 2023-10-09T15:20:57Z

hoomd_validation/util.py

-                expected(bin_edges[:-1]),
+        if plot_rotated:
+            #Y, X = bin_edges[:-1], expected(bin_edges[:-1])
+            Y, X = bin_centers, expected(bin_centers)


Remove the need for the variables by placing the call to plot inside the conditional.

Alternately, name these variables with a PEP8 compliant name.

joaander · 2023-10-09T15:21:51Z

hoomd_validation/util.py

@@ -299,3 +318,46 @@ 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):


Move this function to patchy_particle_pressure.py.

joaander · 2023-10-09T15:32:26Z

hoomd_validation/patchy_particle_pressure.py

+
+    @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)


Suggested change

#@Project.post(lambda j: j.isfile(f'{mode}_{device_name}_complete') or j.sp.pressure <= 0)

tcmoore3 · 2023-10-10T13:47:37Z

Results for the hard sphere + square well + repulsive shoulder potential. It shows similar quantitative agreement NPT and NVT simulations as the Kern–Frenkel and hard sphere + square well systems.

joaander

More fixes to prevent deadlock.

joaander · 2023-10-10T19:15:23Z

hoomd_validation/patchy_particle_pressure.py

+            print(f'completed patchy_particle_pressure_{mode}_{device_name} '
+                  f'{job} in {communicator.walltime} s')


Suggested change

print(f'completed patchy_particle_pressure_{mode}_{device_name} '

f'{job} in {communicator.walltime} s')

print(f'completed patchy_particle_pressure_{mode}_{device_name}: {job}')

See #67.

joaander · 2023-10-10T19:16:02Z

hoomd_validation/patchy_particle_pressure.py

+        print(f'completed patchy_particle_pressure_create_initial_state: '
+              f'{job} in {communicator.walltime} s')


Suggested change

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}')

See #67.

tcmoore3 · 2023-10-19T20:32:03Z

GPU results from 8 replicates at each of 3 statepoints, 1 KF, 1 HSSW, 1 HSSW with a repulsive shoulder. Note that I took the average pressure from the CPU NVT runs, so the deviation from the expected pressure in the 'nvt_gpu' runs are nonzero.

joaander

Thanks! There is one minor change to make, then this is ready to merge.

joaander · 2023-10-20T16:33:17Z

hoomd_validation/patchy_particle_pressure.py

+        # (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, 4.00804, 1.0, 1.0, -1.0),
+        # (1.0, 0.7, -5.0, 1.0, 1.0, -1.0),


Suggested change

# (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, 4.00804, 1.0, 1.0, -1.0),

# (1.0, 0.7, -5.0, 1.0, 1.0, -1.0),

(1.0, 0.95, 10.208625410213362, 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

(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),

Remove unused state points.

joaander

Thanks!

tcmoore3 added 18 commits September 19, 2023 11:34

Add patchy-particle-pressure tests

fba606e

Add patchy_particle_pressure to init and project files

91d4f86

Fix typos

c9b49dc

Add pressure

82be5bb

test

387ecde

undo test

2a38685

test

495aff9

Undo test

3f4d5e6

Change min number of cpu ranks

e4ba892

Update statepoint based on NVT sims

f307199

Update dashboard

6bc2cbe

Update job initialization and aggregation

2293b2e

Update plotting

1de4ca5

Changes during testing

e00350d

Merge branch 'patchy-particle-pressure' of github.com:glotzerlab/hoom…

47c8463

…d-validation into patchy-particle-pressure

Update statepoints based on nvt simulations

28d5a8e

Allow floats for short job walltime in config

760c339

Update patchy particle pressure.

f0e6510

tcmoore3 requested review from a team and joaander and removed request for a team October 8, 2023 22:56

joaander requested changes Oct 9, 2023

View reviewed changes

tcmoore3 added 6 commits October 9, 2023 13:57

Parse max_walltime as a float

308755e

Address PR comments

acb65fa

Update square well + repulsive shoulder pressure

065a5b2

Update filenames

d12c1dc

Update filenames

4d0022c

Fix figure title

7d2fe5f

Fix figure titles

133e718

joaander requested changes Oct 10, 2023

View reviewed changes

Update statepoints and run lengths

042f443

tcmoore3 mentioned this pull request Oct 19, 2023

Fix SDF computation with patchy potentials. glotzerlab/hoomd-blue#1636

Merged

3 tasks

tcmoore3 added 3 commits October 19, 2023 16:22

Add docstring to job_sorter

67d01fe

Uncomment subprojects in init.py

4388fbb

flake8

5327c47

tcmoore3 requested a review from joaander October 19, 2023 20:34

Remove calls to communicator.walltime on rank 0

feb2858

joaander requested changes Oct 20, 2023

View reviewed changes

Remove unused state points.

c4f8f6a

joaander approved these changes Oct 23, 2023

View reviewed changes

joaander merged commit b6ab253 into trunk Oct 23, 2023
2 checks passed

joaander deleted the patchy-particle-pressure branch October 23, 2023 13:41

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Patchy particle pressure #66

Patchy particle pressure #66

tcmoore3 commented Oct 8, 2023 •

edited

Loading

joaander commented Oct 9, 2023

joaander left a comment

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

joaander Oct 9, 2023

tcmoore3 commented Oct 10, 2023

joaander left a comment

joaander Oct 10, 2023

joaander Oct 10, 2023

tcmoore3 commented Oct 19, 2023

joaander left a comment

joaander Oct 20, 2023

joaander left a comment

	self["max_walltime"] = int(config.get("max_walltime", 24))
	self["max_walltime"] = float(config.get("max_walltime", 24))

	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):

		print(f'completed patchy_particle_pressure_{mode}_{device_name} '
		f'{job} in {communicator.walltime} s')

Patchy particle pressure #66

Patchy particle pressure #66

Conversation

tcmoore3 commented Oct 8, 2023 • edited Loading

Description

Motivation and context

How has this been tested?

Checklist:

joaander commented Oct 9, 2023

joaander left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tcmoore3 commented Oct 10, 2023

joaander left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tcmoore3 commented Oct 19, 2023

joaander left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

joaander left a comment

Choose a reason for hiding this comment

tcmoore3 commented Oct 8, 2023 •

edited

Loading