From 0fe09c5d6f9fc0be47109afbeea9474ca7a42e19 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Wed, 10 Jun 2026 18:21:54 -0600 Subject: [PATCH 01/18] Auto rdf --- .gitignore | 5 +- phantomwalk/examples/1-bead-spring-dpd.ipynb | 1091 +++++++++++++++++- phantomwalk/lib/create_system_dpd.py | 5 +- phantomwalk/lib/dpd_utils.py | 21 + 4 files changed, 1108 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index b7faf40..e0bf55e 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,10 @@ __pycache__/ *.py[codz] *$py.class +phantomwalk/examples/rdf.csv +phantomwalk/examples/trajectory.gsd +phantomwalk/examples/log.txt + # C extensions *.so @@ -14,7 +18,6 @@ dist/ downloads/ eggs/ .eggs/ -lib/ lib64/ parts/ sdist/ diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 71d3526..427f783 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -88,7 +88,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 8, "id": "62050dc2-fe73-4d06-8490-32aa8e86daed", "metadata": {}, "outputs": [ @@ -96,26 +96,1085 @@ "name": "stdout", "output_type": "stream", "text": [ - "10000\n", + "1000\n", "\n", - "Running with A=800, gamma=800, k=20000, num_pol=100, num_mon=100\n", - "Total build time: 0.0018270720001964946\n", - "Total build and simulation time: 1.680786765999983\n", - "Finished in time = 1.68s\n" + "Running with A=5000, gamma=1000, k=20000, num_pol=100, num_mon=10\n", + "Total build time: 0.0008774581365287304\n", + "Total build and simulation time: 0.2862441670149565\n", + "Finished in time = 0.29s\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2026-06-10T18:19:24.795131\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.9, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "First nozero RDF value: [0.83999997 0.02113025]\n" ] } ], "source": [ "last_frame, s = create_polymer_system_dpd(\n", " num_pol=100,\n", - " num_mon=100,\n", + " num_mon=10,\n", " density=0.8,\n", " k=20000,\n", " bond_l=1.0,\n", " r_cut=1.15,\n", " kT=1.0,\n", - " A=800,\n", - " gamma=800,\n", + " A=5000,\n", + " gamma=1000,\n", " dt=0.001,\n", " sim_seed=1234,\n", " np_seed=1234,\n", @@ -125,17 +1184,25 @@ " min_pair_dist=1.05,\n", " write=True,\n", " gsd_file_name='trajectory.gsd',\n", - " gsd_write_freq=10,\n", + " gsd_write_freq=100,\n", " log_file_name='log.txt',\n", " log_write_freq=10)\n", "\n", - "print(f\"Finished in time = {s:.2f}s\")" + "print(f\"Finished in time = {s:.2f}s\")\n", + "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", + "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", + "plt.title(\"Radial Distribution Function\")\n", + "plt.xlabel(\"$r$\")\n", + "plt.ylabel(\"$g(r)$\")\n", + "plt.show()\n", + "b = (rdf_data[:,1] !=0).argmax()\n", + "print(\"First nozero RDF value:\", rdf_data[b])" ] }, { "cell_type": "code", "execution_count": null, - "id": "d4262eec-2f8d-44d9-b7f6-bdfd0c55329b", + "id": "31c44f45-e350-431e-bb38-a8ec845aceae", "metadata": {}, "outputs": [], "source": [] diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index 06174e3..d052253 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -121,7 +121,7 @@ def create_polymer_system_dpd( integrator.forces.append(DPD) if write: - add_hoomd_writers( + rdf = add_hoomd_writers( simulation, gsd_file_name, gsd_write_freq, @@ -167,4 +167,7 @@ def create_polymer_system_dpd( end_time = time.perf_counter() total_time = end_time - start_time print("Total build and simulation time:", end_time - start_time) + np.savetxt( + "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)" + ) return snap, total_time diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index f96d09f..2711441 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -128,6 +128,22 @@ def add_hoomd_writers( and does not return a value. """ + + class FreudRDFCalc(hoomd.custom.Action): + """Compute RDF periodically as the simulation progresses.""" + + def __init__(self, sim, rdf): + self._sim = sim + self._rdf = rdf + + def act(self, timestep): + snap = self._sim.state.get_snapshot() + self._rdf.compute(system=snap, reset=True) + + rdf = freud.density.RDF(bins=50, r_max=4) + rdf_calc = FreudRDFCalc(sim, rdf) + + gsd_logger = hoomd.logging.Logger( categories=["scalar", "string", "sequence"] ) @@ -174,8 +190,13 @@ def add_hoomd_writers( logger=logger, max_header_len=None, ) + + rdf_action = hoomd.write.CustomWriter(action=rdf_calc, trigger=log_trigger) + sim.operations.writers.append(rdf_action) + sim.operations.writers.append(gsd_writer) sim.operations.writers.append(table_file) + return rdf def check_pair_energy(energy_idx=-1, log_file_name="log.txt"): """Check whether the pair interaction energy has equilibrated. From 394b620d4239c6ca752400d7608af0b48a61ae8a Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Wed, 10 Jun 2026 18:22:42 -0600 Subject: [PATCH 02/18] no notebook output --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 1100 +----------------- 1 file changed, 4 insertions(+), 1096 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 427f783..d77b905 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -21,30 +21,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "88429372", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", - "\n", - "****** PyMBAR will use 64-bit JAX! *******\n", - "* JAX is currently set to 32-bit bitsize *\n", - "* which is its default. *\n", - "* *\n", - "* PyMBAR requires 64-bit mode and WILL *\n", - "* enable JAX's 64-bit mode when called. *\n", - "* *\n", - "* This MAY cause problems with other *\n", - "* Uses of JAX in the same code. *\n", - "******************************************\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "import sys\n", "import os\n", @@ -88,1082 +68,10 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "62050dc2-fe73-4d06-8490-32aa8e86daed", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000\n", - "\n", - "Running with A=5000, gamma=1000, k=20000, num_pol=100, num_mon=10\n", - "Total build time: 0.0008774581365287304\n", - "Total build and simulation time: 0.2862441670149565\n", - "Finished in time = 0.29s\n" - ] - }, - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " 2026-06-10T18:19:24.795131\n", - " image/svg+xml\n", - " \n", - " \n", - " Matplotlib v3.10.9, https://matplotlib.org/\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "\n" - ], - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "First nozero RDF value: [0.83999997 0.02113025]\n" - ] - } - ], + "outputs": [], "source": [ "last_frame, s = create_polymer_system_dpd(\n", " num_pol=100,\n", From 2c63717049793fe0709bab6e317dc746938675cb Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Fri, 12 Jun 2026 15:17:07 -0600 Subject: [PATCH 03/18] output shortest distance --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 73 +++++++++++--------- phantomwalk/examples/3-dpd-to-lj-wca.ipynb | 71 +++---------------- phantomwalk/lib/create_system_dpd.py | 6 +- 3 files changed, 54 insertions(+), 96 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index d77b905..8aa787d 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -73,38 +73,47 @@ "metadata": {}, "outputs": [], "source": [ - "last_frame, s = create_polymer_system_dpd(\n", - " num_pol=100,\n", - " num_mon=10,\n", - " density=0.8,\n", - " k=20000,\n", - " bond_l=1.0,\n", - " r_cut=1.15,\n", - " kT=1.0,\n", - " A=5000,\n", - " gamma=1000,\n", - " dt=0.001,\n", - " sim_seed=1234,\n", - " np_seed=1234,\n", - " sim_steps_incr=100,\n", - " loop_timeout=60,\n", - " energy=True,\n", - " min_pair_dist=1.05,\n", - " write=True,\n", - " gsd_file_name='trajectory.gsd',\n", - " gsd_write_freq=100,\n", - " log_file_name='log.txt',\n", - " log_write_freq=10)\n", - "\n", - "print(f\"Finished in time = {s:.2f}s\")\n", - "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", - "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", - "plt.title(\"Radial Distribution Function\")\n", - "plt.xlabel(\"$r$\")\n", - "plt.ylabel(\"$g(r)$\")\n", - "plt.show()\n", - "b = (rdf_data[:,1] !=0).argmax()\n", - "print(\"First nozero RDF value:\", rdf_data[b])" + "nruns = 5\n", + "times = []\n", + "rs = []\n", + "for i in range(nruns):\n", + " seed = np.random.randint(50000)\n", + " last_frame, s = create_polymer_system_dpd(\n", + " num_pol=100,\n", + " num_mon=10,\n", + " density=1.3,\n", + " k=5000,\n", + " bond_l=1.0,\n", + " r_cut=1.1,\n", + " kT=1.0,\n", + " A=5000,\n", + " gamma=1000,\n", + " dt=0.002,\n", + " sim_seed=seed,\n", + " np_seed=seed,\n", + " sim_steps_incr=100,\n", + " loop_timeout=60,\n", + " energy=True,\n", + " min_pair_dist=.8,\n", + " write=True,\n", + " gsd_file_name='trajectory.gsd',\n", + " gsd_write_freq=100,\n", + " log_file_name='log.txt',\n", + " log_write_freq=10)\n", + " \n", + " #print(f\"Finished in time = {s:.2f}s\")\n", + " rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", + " #plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", + " #plt.title(\"Radial Distribution Function\")\n", + " #plt.xlabel(\"$r$\")\n", + " #plt.ylabel(\"$g(r)$\")\n", + " #lt.show()\n", + " b = (rdf_data[:,1] !=0).argmax()\n", + " #print(\"First nozero RDF value:\", rdf_data[b])\n", + " times.append(s)\n", + " rs.append(rdf_data[b][0])\n", + " print(\"{:.2f}s, r = {:.2f} g(r) = {:.4f}\".format(s,rdf_data[b][0],rdf_data[b][1]))\n", + "print(\"\\n\\n{:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(np.average(times),np.std(times), np.average(rs), np.std(rs)))" ] }, { diff --git a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb index cc410a7..5ece95f 100644 --- a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb +++ b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb @@ -11,30 +11,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "88429372", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", - "\n", - "****** PyMBAR will use 64-bit JAX! *******\n", - "* JAX is currently set to 32-bit bitsize *\n", - "* which is its default. *\n", - "* *\n", - "* PyMBAR requires 64-bit mode and WILL *\n", - "* enable JAX's 64-bit mode when called. *\n", - "* *\n", - "* This MAY cause problems with other *\n", - "* Uses of JAX in the same code. *\n", - "******************************************\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "import sys\n", "import os\n", @@ -58,30 +38,17 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "1be674a6-7412-4f5b-aad2-24cb3cc0f869", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "10000\n", - "\n", - "Running with A=1000, gamma=800, k=20000, num_pol=100, num_mon=100\n", - "Total build time: 0.002074115000141319\n", - "Total build and simulation time: 2.389522753000165\n", - "Finished in time = 2.39s\n" - ] - } - ], + "outputs": [], "source": [ "num_pol=100\n", "num_mon=100\n", "N = num_pol*num_mon\n", - "density=0.9\n", + "density=0.8\n", "A=1000\n", - "r_cut=1.13\n", + "r_cut=1.1\n", "min_pair_dist=1.05\n", "dpd_final_frame, s = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", @@ -90,7 +57,7 @@ " k=20000,\n", " r_cut=r_cut,\n", " A=A,\n", - " gamma=800,\n", + " gamma=1000,\n", " sim_seed=1234,\n", " np_seed=1234,\n", " energy=True,\n", @@ -109,7 +76,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "6139f233-10e2-4288-8cc0-cc78717d9f35", "metadata": {}, "outputs": [], @@ -267,28 +234,10 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "48f19fa1-6ec6-4244-a23c-fd9d774da3f1", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "LJ simulation finished.\n" - ] - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "run_lj_simulation(dpd_final_frame=dpd_final_frame)" ] diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index d052253..51d0fe7 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -90,7 +90,7 @@ def create_polymer_system_dpd( execution time of the DPD workflow, build + simulation wall time ''' - print(num_pol*num_mon) + #print(num_pol*num_mon) print(f"\nRunning with A={A}, gamma={gamma}, k={k}, " f"num_pol={num_pol}, num_mon={num_mon}") start_time = time.perf_counter() @@ -104,7 +104,7 @@ def create_polymer_system_dpd( ) build_stop = time.perf_counter() - print("Total build time: ", build_stop-start_time) + #print("Total build time: ", build_stop-start_time) harmonic = hoomd.md.bond.Harmonic() harmonic.params["b"] = dict(r0=bond_l, k=k) integrator = hoomd.md.Integrator(dt=dt) @@ -166,7 +166,7 @@ def create_polymer_system_dpd( end_time = time.perf_counter() total_time = end_time - start_time - print("Total build and simulation time:", end_time - start_time) + #print("Total build and simulation time:", end_time - start_time) np.savetxt( "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)" ) From c7e40f72ceb59da6cb01ee8c56762c6daa7b29cd Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Fri, 12 Jun 2026 15:45:28 -0600 Subject: [PATCH 04/18] pseudocode for AA implementation --- phantomwalk/lib/dpd_utils.py | 41 ++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 2711441..28ff1e9 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -5,6 +5,47 @@ import time from cmeutils.sampling import is_equilibrated +def initialize_from_topology(bond_list, density, box, bond_length=1.0, seed=1234): + ''' + Currently expects bond_list to be list of ordered tuples with indices from 0 to N-1, each tuple describing a bond + + ''' + rng = np.random.default_rng(seed) + N = np.max(bond_list[:,1])+1 #assumes ordered tuples + + #Overall algorithm: + # give all particles initial random coords + positions = rng.uniform(0, np.max(box), size=(N, 3)) + # Generate a uniform-sphere-delta for each bond + # Cumulative sum of deltas gives positions for a chain off of a branch (or origin) + thetas = rng.uniform(0,2*np.pi,size=N) + phis = np.arccos(rng.uniform(-1,1,size=N) + x = np.sin(phis)*np.cos(thetas) + y = np.sin(phis)*np.sin(thetas) + z = np.cos(phis) + deltas = np.stack([x,y,z],axis=2) * bond_length + displacements = np.cumsum(deltas, axis=1) + #positions_view[:, 1:, :] = starts[:, None, :] + displacements #deprecated + #TODO: loop over bonds and add displacements here + + #pbc + #TODO Wrap coordinates along each axis separately + positions %= L + positions -= L/2 + + #TODO: Fold in PR 83 mupt, mupt.builder.all_atom_dpd functionality. + + frame = gsd.hoomd.Frame() + frame.particles.types = ['A'] #update with all types + frame.particles.N = N + frame.particles.position = positions + frame.bonds.N = len(bonds_list) + frame.bonds.group = bond_list + frame.bonds.types = ['b'] #update with all types + frame.configuration.box = box + + return frame + def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, seed=1234): ''' Create a HOOMD snapshot of a cubic box with the number density given by input parameters. Configure particles using a random walk. From 5bf067aae2d9a68d711504a1abad466a6e857fe2 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Fri, 12 Jun 2026 18:01:59 -0600 Subject: [PATCH 05/18] sketching pieces --- phantomwalk/lib/dpd_utils.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 28ff1e9..5d4e16d 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -5,6 +5,24 @@ import time from cmeutils.sampling import is_equilibrated + + +#https://github.com/joelaforet/mupt/blob/issue-77-aa-dpd-builder/mupt/builders/all_atom_dpd.py +class _ParameterTables: + """HOOMD-ready bonded and vdW parameter tables.""" + + bond_params: dict[str, dict[str, float]] = field(default_factory=dict) + angle_params: dict[str, dict[str, float]] = field(default_factory=dict) + dihedral_params: dict[str, dict[str, float]] = field(default_factory=dict) + improper_params: dict[str, dict[str, float]] = field(default_factory=dict) + bond_type_by_group: dict[tuple[int, int], str] = field(default_factory=dict) + angle_type_by_group: dict[tuple[int, int, int], str] = field(default_factory=dict) + dihedral_type_by_group: dict[tuple[int, int, int, int], list[str]] = field(default_factory=dict) + improper_type_by_group: dict[tuple[int, int, int, int], list[str]] = field(default_factory=dict) + atom_epsilons: dict[int, float] = field(default_factory=dict) + atom_types_by_global: dict[int, str] = field(default_factory=dict) + epsilon_by_type: dict[str, float] = field(default_factory=dict) + def initialize_from_topology(bond_list, density, box, bond_length=1.0, seed=1234): ''' Currently expects bond_list to be list of ordered tuples with indices from 0 to N-1, each tuple describing a bond From 90c7d6214efa6a1a11b10b19e5b761f51e537b31 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Fri, 12 Jun 2026 21:55:00 -0600 Subject: [PATCH 06/18] Dialed in some fast big simulations at 1.3 density --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 2781 +++++++++++++++++- phantomwalk/examples/2-dpd-energy-gsd.ipynb | 2713 ++++++++++++++++- phantomwalk/lib/create_system_dpd.py | 4 +- phantomwalk/lib/dpd_utils.py | 2 +- 4 files changed, 5458 insertions(+), 42 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 5f6fc60..fabd1be 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -21,7 +21,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "88429372", "metadata": {}, "outputs": [], @@ -68,59 +68,2790 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "id": "62050dc2-fe73-4d06-8490-32aa8e86daed", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.68s, r = 0.91 g(r) = 0.0591\n", + "0.66s, r = 0.91 g(r) = 0.0296\n", + "0.66s, r = 0.91 g(r) = 0.0444\n", + "0.67s, r = 0.91 g(r) = 0.0296\n", + "0.65s, r = 0.89 g(r) = 0.0077\n", + "\n", + "\n", + "0.66 (0.01)s , r = 0.91 (0.01)\n" + ] + } + ], "source": [ "nruns = 5\n", "times = []\n", "rs = []\n", + "\n", + "npoly=10\n", + "nmono=100\n", + "N=npoly*nmono\n", + "A=50000\n", + "k=50000\n", + "dt=0.001\n", + "r_cut=1.01\n", + "min_pair_dist=0.988\n", + "gamma = 1200\n", + "density=1.3\n", + "\n", "for i in range(nruns):\n", " seed = np.random.randint(50000)\n", " last_frame, s = create_polymer_system_dpd(\n", - " num_pol=100,\n", - " num_mon=10,\n", - " density=1.3,\n", - " k=5000,\n", " bond_l=1.0,\n", - " r_cut=1.1,\n", + " num_pol=npoly,\n", + " num_mon=nmono,\n", " kT=1.0,\n", - " A=5000,\n", - " gamma=1000,\n", - " dt=0.002,\n", " sim_seed=seed,\n", " np_seed=seed,\n", " sim_steps_incr=100,\n", - " loop_timeout=60,\n", + " loop_timeout=600,\n", " energy=True,\n", - " min_pair_dist=.8,\n", " write=True,\n", " gsd_file_name='trajectory.gsd',\n", " gsd_write_freq=100,\n", " log_file_name='log.txt',\n", - " log_write_freq=10)\n", - " \n", - " #print(f\"Finished in time = {s:.2f}s\")\n", + " log_write_freq=10,\n", + " density=density,\n", + " \n", + " dt=dt,\n", + " r_cut=r_cut,\n", + " min_pair_dist=min_pair_dist,\n", + " A=A,\n", + " k=k,\n", + " gamma=gamma)\n", " rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", - " #plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", - " #plt.title(\"Radial Distribution Function\")\n", - " #plt.xlabel(\"$r$\")\n", - " #plt.ylabel(\"$g(r)$\")\n", - " #lt.show()\n", " b = (rdf_data[:,1] !=0).argmax()\n", - " #print(\"First nozero RDF value:\", rdf_data[b])\n", " times.append(s)\n", " rs.append(rdf_data[b][0])\n", " print(\"{:.2f}s, r = {:.2f} g(r) = {:.4f}\".format(s,rdf_data[b][0],rdf_data[b][1]))\n", - "print(\"\\n\\n{:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(np.average(times),np.std(times), np.average(rs), np.std(rs)))" + "print(\"\\n\\n{:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(np.average(times),np.std(times), np.average(rs), np.std(rs)))\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "31c44f45-e350-431e-bb38-a8ec845aceae", "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2026-06-12T21:53:09.622157\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.9, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[9.99999978e-03 0.00000000e+00]\n", + " [2.99999993e-02 0.00000000e+00]\n", + " [4.99999970e-02 0.00000000e+00]\n", + " [7.00000003e-02 0.00000000e+00]\n", + " [8.99999961e-02 0.00000000e+00]\n", + " [1.09999999e-01 0.00000000e+00]\n", + " [1.29999995e-01 0.00000000e+00]\n", + " [1.50000006e-01 0.00000000e+00]\n", + " [1.69999987e-01 0.00000000e+00]\n", + " [1.89999998e-01 0.00000000e+00]\n", + " [2.09999993e-01 0.00000000e+00]\n", + " [2.29999989e-01 0.00000000e+00]\n", + " [2.50000000e-01 0.00000000e+00]\n", + " [2.69999981e-01 0.00000000e+00]\n", + " [2.89999992e-01 0.00000000e+00]\n", + " [3.10000002e-01 0.00000000e+00]\n", + " [3.29999983e-01 0.00000000e+00]\n", + " [3.49999994e-01 0.00000000e+00]\n", + " [3.70000005e-01 0.00000000e+00]\n", + " [3.89999986e-01 0.00000000e+00]\n", + " [4.09999967e-01 0.00000000e+00]\n", + " [4.30000007e-01 0.00000000e+00]\n", + " [4.49999988e-01 0.00000000e+00]\n", + " [4.69999969e-01 0.00000000e+00]\n", + " [4.90000010e-01 0.00000000e+00]\n", + " [5.09999990e-01 0.00000000e+00]\n", + " [5.29999971e-01 0.00000000e+00]\n", + " [5.49999952e-01 0.00000000e+00]\n", + " [5.69999993e-01 0.00000000e+00]\n", + " [5.89999974e-01 0.00000000e+00]\n", + " [6.10000014e-01 0.00000000e+00]\n", + " [6.29999995e-01 0.00000000e+00]\n", + " [6.49999976e-01 0.00000000e+00]\n", + " [6.69999957e-01 0.00000000e+00]\n", + " [6.89999998e-01 0.00000000e+00]\n", + " [7.09999979e-01 0.00000000e+00]\n", + " [7.30000019e-01 0.00000000e+00]\n", + " [7.50000000e-01 0.00000000e+00]\n", + " [7.69999981e-01 0.00000000e+00]\n", + " [7.89999962e-01 0.00000000e+00]\n", + " [8.09999943e-01 0.00000000e+00]\n", + " [8.29999983e-01 0.00000000e+00]\n", + " [8.49999964e-01 0.00000000e+00]\n", + " [8.70000005e-01 0.00000000e+00]\n", + " [8.89999986e-01 7.72768073e-03]\n", + " [9.09999967e-01 2.21752264e-02]\n", + " [9.29999948e-01 3.18475366e-01]\n", + " [9.49999988e-01 2.19071674e+00]\n", + " [9.69999969e-01 6.31044579e+00]\n", + " [9.90000010e-01 9.31815243e+00]\n", + " [1.00999999e+00 5.99454355e+00]\n", + " [1.02999997e+00 2.33675718e+00]\n", + " [1.04999995e+00 1.23811281e+00]\n", + " [1.06999993e+00 1.01582968e+00]\n", + " [1.08999991e+00 8.55239391e-01]\n", + " [1.11000001e+00 7.35277355e-01]\n", + " [1.13000000e+00 6.13604367e-01]\n", + " [1.14999998e+00 5.92449248e-01]\n", + " [1.16999996e+00 5.50008595e-01]\n", + " [1.18999994e+00 5.10064900e-01]\n", + " [1.21000004e+00 4.51530993e-01]\n", + " [1.23000002e+00 4.89567488e-01]\n", + " [1.25000000e+00 4.27016288e-01]\n", + " [1.26999998e+00 4.59215075e-01]\n", + " [1.28999996e+00 5.03942072e-01]\n", + " [1.30999994e+00 4.35167283e-01]\n", + " [1.32999992e+00 3.94493490e-01]\n", + " [1.34999990e+00 5.10520041e-01]\n", + " [1.37000000e+00 5.54431260e-01]\n", + " [1.38999999e+00 5.19582152e-01]\n", + " [1.40999997e+00 5.11105120e-01]\n", + " [1.42999995e+00 5.29835641e-01]\n", + " [1.44999993e+00 5.91018200e-01]\n", + " [1.47000003e+00 6.06200635e-01]\n", + " [1.49000001e+00 6.86542988e-01]\n", + " [1.50999999e+00 6.92638040e-01]\n", + " [1.52999997e+00 6.27580464e-01]\n", + " [1.54999995e+00 6.98118865e-01]\n", + " [1.56999993e+00 7.89712429e-01]\n", + " [1.58999991e+00 7.62705624e-01]\n", + " [1.60999990e+00 9.11538541e-01]\n", + " [1.63000000e+00 1.03215218e+00]\n", + " [1.64999998e+00 1.16467202e+00]\n", + " [1.66999996e+00 1.46617603e+00]\n", + " [1.68999994e+00 1.50240195e+00]\n", + " [1.70999992e+00 1.47793388e+00]\n", + " [1.73000002e+00 1.47463214e+00]\n", + " [1.75000000e+00 1.33518815e+00]\n", + " [1.76999998e+00 1.34231174e+00]\n", + " [1.78999996e+00 1.24561477e+00]\n", + " [1.80999994e+00 1.21823835e+00]\n", + " [1.82999992e+00 1.19907069e+00]\n", + " [1.84999990e+00 1.19831800e+00]\n", + " [1.87000000e+00 1.33037162e+00]\n", + " [1.88999999e+00 1.32635713e+00]\n", + " [1.90999997e+00 1.49168754e+00]\n", + " [1.92999995e+00 1.49872899e+00]\n", + " [1.94999993e+00 1.51322103e+00]\n", + " [1.96999991e+00 1.38643694e+00]\n", + " [1.99000001e+00 1.11911285e+00]]\n" + ] + } + ], + "source": [ + "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", + "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", + "plt.title(\"Radial Distribution Function\")\n", + "plt.xlabel(\"$r$\")\n", + "plt.ylabel(\"$g(r)$\")\n", + "plt.show()\n", + "print(rdf_data)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "5f279783-91f5-4560-8d38-c8be52c09443", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total steps 13200\n", + "DPD Energy Cutoff= 65.82111560065738\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2026-06-12T21:53:10.528353\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.9, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from dpd_utils import calculate_pair_energy\n", + "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=npoly,num_mon=nmono,density=density)\n", + "log = np.genfromtxt(\"log.txt\", names=True)\n", + "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", + "pairs = log[\"mdpairDPDenergy\"]\n", + "print(\"Total steps\",len(pe)*100+100)\n", + "print(\"DPD Energy Cutoff= \",U/N)\n", + "\n", + "slice_idx = 100\n", + "short_idx= None\n", + "x_values = range(slice_idx, len(pe[:short_idx]))\n", + "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", + "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", + "plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", + "plt.title(\"Energy vs Time\")\n", + "plt.xlabel(\"Frame (time/log_freq)\")\n", + "plt.ylabel(\"Energy/N\")\n", + "\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59980585-3703-42c8-8c71-81e76eefcbb3", + "metadata": {}, "outputs": [], "source": [] } @@ -141,7 +2872,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.13" + "version": "3.12.13" } }, "nbformat": 4, diff --git a/phantomwalk/examples/2-dpd-energy-gsd.ipynb b/phantomwalk/examples/2-dpd-energy-gsd.ipynb index 449a253..bb5ebc0 100644 --- a/phantomwalk/examples/2-dpd-energy-gsd.ipynb +++ b/phantomwalk/examples/2-dpd-energy-gsd.ipynb @@ -11,10 +11,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "88429372", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + "\n", + "****** PyMBAR will use 64-bit JAX! *******\n", + "* JAX is currently set to 32-bit bitsize *\n", + "* which is its default. *\n", + "* *\n", + "* PyMBAR requires 64-bit mode and WILL *\n", + "* enable JAX's 64-bit mode when called. *\n", + "* *\n", + "* This MAY cause problems with other *\n", + "* Uses of JAX in the same code. *\n", + "******************************************\n", + "\n" + ] + } + ], "source": [ "import sys\n", "import os\n", @@ -40,15 +60,26 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "1be674a6-7412-4f5b-aad2-24cb3cc0f869", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Running with A=800, gamma=800, k=20000, num_pol=100, num_mon=100\n", + "Simulation timed out\n", + "Finished in time = 0.00s\n" + ] + } + ], "source": [ "num_pol=100\n", "num_mon=100\n", "N = num_pol*num_mon\n", - "density=0.9\n", + "density=1.3\n", "A=800\n", "r_cut=1.13\n", "min_pair_dist=1.05\n", @@ -70,10 +101,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "id": "47e9b287-e7bd-4fbb-bfb4-b27c59206eb1", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "124469.44818043037\n" + ] + } + ], "source": [ "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=num_pol,num_mon=num_mon,density=density)\n", "print(U)" @@ -89,21 +128,1414 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "id": "1b04fd76-9337-421e-bcfc-2b351c73e700", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total steps 104200\n", + "DPD Energy Cutoff= 12.446944818043036\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2026-06-12T18:11:55.976843\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.9, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "log = np.genfromtxt(\"log.txt\", names=True)\n", "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", "pairs = log[\"mdpairDPDenergy\"]\n", "print(\"Total steps\",len(pe)*100+100)\n", + "print(\"DPD Energy Cutoff= \",U/N)\n", + "\n", "slice_idx = 0\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", - "plt.plot(x_values,pe[slice_idx:short_idx], label=\"potential energy\")\n", - "plt.plot(x_values,pairs[slice_idx:short_idx], label=\"DPD pair energy\")\n", - "plt.hlines(U,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", + "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", + "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", + "plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", "plt.title(\"Energy vs Time\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy\")\n", @@ -144,9 +1576,1262 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "id": "133f7da1-415e-49d6-9731-b7f642a62469", "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2026-06-12T21:33:21.305317\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.9, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[9.99999978e-03 0.00000000e+00]\n", + " [2.99999993e-02 0.00000000e+00]\n", + " [4.99999970e-02 0.00000000e+00]\n", + " [7.00000003e-02 0.00000000e+00]\n", + " [8.99999961e-02 0.00000000e+00]\n", + " [1.09999999e-01 0.00000000e+00]\n", + " [1.29999995e-01 0.00000000e+00]\n", + " [1.50000006e-01 0.00000000e+00]\n", + " [1.69999987e-01 0.00000000e+00]\n", + " [1.89999998e-01 0.00000000e+00]\n", + " [2.09999993e-01 0.00000000e+00]\n", + " [2.29999989e-01 0.00000000e+00]\n", + " [2.50000000e-01 0.00000000e+00]\n", + " [2.69999981e-01 0.00000000e+00]\n", + " [2.89999992e-01 0.00000000e+00]\n", + " [3.10000002e-01 0.00000000e+00]\n", + " [3.29999983e-01 0.00000000e+00]\n", + " [3.49999994e-01 0.00000000e+00]\n", + " [3.70000005e-01 0.00000000e+00]\n", + " [3.89999986e-01 0.00000000e+00]\n", + " [4.09999967e-01 0.00000000e+00]\n", + " [4.30000007e-01 0.00000000e+00]\n", + " [4.49999988e-01 0.00000000e+00]\n", + " [4.69999969e-01 0.00000000e+00]\n", + " [4.90000010e-01 0.00000000e+00]\n", + " [5.09999990e-01 0.00000000e+00]\n", + " [5.29999971e-01 0.00000000e+00]\n", + " [5.49999952e-01 0.00000000e+00]\n", + " [5.69999993e-01 0.00000000e+00]\n", + " [5.89999974e-01 0.00000000e+00]\n", + " [6.10000014e-01 0.00000000e+00]\n", + " [6.29999995e-01 0.00000000e+00]\n", + " [6.49999976e-01 0.00000000e+00]\n", + " [6.69999957e-01 0.00000000e+00]\n", + " [6.89999998e-01 0.00000000e+00]\n", + " [7.09999979e-01 0.00000000e+00]\n", + " [7.30000019e-01 0.00000000e+00]\n", + " [7.50000000e-01 0.00000000e+00]\n", + " [7.69999981e-01 0.00000000e+00]\n", + " [7.89999962e-01 0.00000000e+00]\n", + " [8.09999943e-01 0.00000000e+00]\n", + " [8.29999983e-01 0.00000000e+00]\n", + " [8.49999964e-01 0.00000000e+00]\n", + " [8.70000005e-01 1.61740478e-04]\n", + " [8.89999986e-01 1.85464346e-03]\n", + " [9.09999967e-01 3.42237689e-02]\n", + " [9.29999948e-01 3.75659347e-01]\n", + " [9.49999988e-01 2.12913251e+00]\n", + " [9.69999969e-01 6.38623619e+00]\n", + " [9.90000010e-01 9.32189941e+00]\n", + " [1.00999999e+00 5.86049223e+00]\n", + " [1.02999997e+00 2.40495610e+00]\n", + " [1.04999995e+00 1.30290556e+00]\n", + " [1.06999993e+00 9.48677957e-01]\n", + " [1.08999991e+00 7.72342980e-01]\n", + " [1.11000001e+00 6.71785176e-01]\n", + " [1.13000000e+00 6.23863041e-01]\n", + " [1.14999998e+00 5.86108208e-01]\n", + " [1.16999996e+00 5.55777013e-01]\n", + " [1.18999994e+00 5.20482302e-01]\n", + " [1.21000004e+00 5.19804120e-01]\n", + " [1.23000002e+00 4.97619092e-01]\n", + " [1.25000000e+00 4.86955285e-01]\n", + " [1.26999998e+00 4.79405373e-01]\n", + " [1.28999996e+00 4.77862120e-01]\n", + " [1.30999994e+00 4.69837993e-01]\n", + " [1.32999992e+00 4.74741757e-01]\n", + " [1.34999990e+00 4.84658182e-01]\n", + " [1.37000000e+00 4.91095632e-01]\n", + " [1.38999999e+00 4.97499943e-01]\n", + " [1.40999997e+00 5.17201483e-01]\n", + " [1.42999995e+00 5.36301434e-01]\n", + " [1.44999993e+00 5.55993795e-01]\n", + " [1.47000003e+00 5.77137053e-01]\n", + " [1.49000001e+00 6.08183384e-01]\n", + " [1.50999999e+00 6.34086013e-01]\n", + " [1.52999997e+00 6.68373168e-01]\n", + " [1.54999995e+00 7.17355371e-01]\n", + " [1.56999993e+00 7.84075201e-01]\n", + " [1.58999991e+00 8.83406878e-01]\n", + " [1.60999990e+00 9.75062847e-01]\n", + " [1.63000000e+00 1.07528138e+00]\n", + " [1.64999998e+00 1.19554257e+00]\n", + " [1.66999996e+00 1.33810508e+00]\n", + " [1.68999994e+00 1.46330953e+00]\n", + " [1.70999992e+00 1.51521719e+00]\n", + " [1.73000002e+00 1.46107221e+00]\n", + " [1.75000000e+00 1.35471618e+00]\n", + " [1.76999998e+00 1.29108131e+00]\n", + " [1.78999996e+00 1.24045670e+00]\n", + " [1.80999994e+00 1.21756566e+00]\n", + " [1.82999992e+00 1.23899102e+00]\n", + " [1.84999990e+00 1.25472844e+00]\n", + " [1.87000000e+00 1.29916048e+00]\n", + " [1.88999999e+00 1.36729598e+00]\n", + " [1.90999997e+00 1.43487251e+00]\n", + " [1.92999995e+00 1.48630548e+00]\n", + " [1.94999993e+00 1.46616626e+00]\n", + " [1.96999991e+00 1.35239911e+00]\n", + " [1.99000001e+00 1.14362824e+00]]\n" + ] + } + ], + "source": [ + "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", + "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", + "plt.title(\"Radial Distribution Function\")\n", + "plt.xlabel(\"$r$\")\n", + "plt.ylabel(\"$g(r)$\")\n", + "plt.show()\n", + "print(rdf_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "315bac34-3982-4a53-8466-3ff095590187", + "metadata": {}, "outputs": [], "source": [] } @@ -167,7 +2852,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.13" + "version": "3.12.13" } }, "nbformat": 4, diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index 6273bc2..7f30f32 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -91,8 +91,8 @@ def create_polymer_system_dpd( ''' #print(num_pol*num_mon) - print(f"\nRunning with A={A}, gamma={gamma}, k={k}, " - f"num_pol={num_pol}, num_mon={num_mon}") + #print(f"\nRunning with A={A}, gamma={gamma}, k={k}, " + # f"num_pol={num_pol}, num_mon={num_mon}") start_time = time.perf_counter() frame = initialize_snapshot_rand_walk( diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 47458de..87d9b66 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -140,7 +140,7 @@ def act(self, timestep): snap = self._sim.state.get_snapshot() self._rdf.compute(system=snap, reset=True) - rdf = freud.density.RDF(bins=50, r_max=4) + rdf = freud.density.RDF(bins=100, r_max=2.0) rdf_calc = FreudRDFCalc(sim, rdf) From 4a74c328f37d9d95e1701359b1201a22f300761d Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Fri, 12 Jun 2026 22:04:02 -0600 Subject: [PATCH 07/18] clean up gitignore, no notebook output --- .gitignore | 1 - phantomwalk/examples/1-bead-spring-dpd.ipynb | 2700 +----------------- phantomwalk/examples/2-dpd-energy-gsd.ipynb | 2687 +---------------- 3 files changed, 19 insertions(+), 5369 deletions(-) diff --git a/.gitignore b/.gitignore index 6366fcc..6cf64c2 100644 --- a/.gitignore +++ b/.gitignore @@ -16,7 +16,6 @@ phantomwalk/examples/log.txt # Distribution / packaging dist/ -<<<<<<< HEAD downloads/ eggs/ .eggs/ diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index fabd1be..3c7134e 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -21,7 +21,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "88429372", "metadata": {}, "outputs": [], @@ -68,25 +68,10 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "62050dc2-fe73-4d06-8490-32aa8e86daed", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "0.68s, r = 0.91 g(r) = 0.0591\n", - "0.66s, r = 0.91 g(r) = 0.0296\n", - "0.66s, r = 0.91 g(r) = 0.0444\n", - "0.67s, r = 0.91 g(r) = 0.0296\n", - "0.65s, r = 0.89 g(r) = 0.0077\n", - "\n", - "\n", - "0.66 (0.01)s , r = 0.91 (0.01)\n" - ] - } - ], + "outputs": [], "source": [ "nruns = 5\n", "times = []\n", @@ -138,1247 +123,10 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "31c44f45-e350-431e-bb38-a8ec845aceae", "metadata": {}, - "outputs": [ - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " 2026-06-12T21:53:09.622157\n", - " image/svg+xml\n", - " \n", - " \n", - " Matplotlib v3.10.9, https://matplotlib.org/\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "\n" - ], - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[9.99999978e-03 0.00000000e+00]\n", - " [2.99999993e-02 0.00000000e+00]\n", - " [4.99999970e-02 0.00000000e+00]\n", - " [7.00000003e-02 0.00000000e+00]\n", - " [8.99999961e-02 0.00000000e+00]\n", - " [1.09999999e-01 0.00000000e+00]\n", - " [1.29999995e-01 0.00000000e+00]\n", - " [1.50000006e-01 0.00000000e+00]\n", - " [1.69999987e-01 0.00000000e+00]\n", - " [1.89999998e-01 0.00000000e+00]\n", - " [2.09999993e-01 0.00000000e+00]\n", - " [2.29999989e-01 0.00000000e+00]\n", - " [2.50000000e-01 0.00000000e+00]\n", - " [2.69999981e-01 0.00000000e+00]\n", - " [2.89999992e-01 0.00000000e+00]\n", - " [3.10000002e-01 0.00000000e+00]\n", - " [3.29999983e-01 0.00000000e+00]\n", - " [3.49999994e-01 0.00000000e+00]\n", - " [3.70000005e-01 0.00000000e+00]\n", - " [3.89999986e-01 0.00000000e+00]\n", - " [4.09999967e-01 0.00000000e+00]\n", - " [4.30000007e-01 0.00000000e+00]\n", - " [4.49999988e-01 0.00000000e+00]\n", - " [4.69999969e-01 0.00000000e+00]\n", - " [4.90000010e-01 0.00000000e+00]\n", - " [5.09999990e-01 0.00000000e+00]\n", - " [5.29999971e-01 0.00000000e+00]\n", - " [5.49999952e-01 0.00000000e+00]\n", - " [5.69999993e-01 0.00000000e+00]\n", - " [5.89999974e-01 0.00000000e+00]\n", - " [6.10000014e-01 0.00000000e+00]\n", - " [6.29999995e-01 0.00000000e+00]\n", - " [6.49999976e-01 0.00000000e+00]\n", - " [6.69999957e-01 0.00000000e+00]\n", - " [6.89999998e-01 0.00000000e+00]\n", - " [7.09999979e-01 0.00000000e+00]\n", - " [7.30000019e-01 0.00000000e+00]\n", - " [7.50000000e-01 0.00000000e+00]\n", - " [7.69999981e-01 0.00000000e+00]\n", - " [7.89999962e-01 0.00000000e+00]\n", - " [8.09999943e-01 0.00000000e+00]\n", - " [8.29999983e-01 0.00000000e+00]\n", - " [8.49999964e-01 0.00000000e+00]\n", - " [8.70000005e-01 0.00000000e+00]\n", - " [8.89999986e-01 7.72768073e-03]\n", - " [9.09999967e-01 2.21752264e-02]\n", - " [9.29999948e-01 3.18475366e-01]\n", - " [9.49999988e-01 2.19071674e+00]\n", - " [9.69999969e-01 6.31044579e+00]\n", - " [9.90000010e-01 9.31815243e+00]\n", - " [1.00999999e+00 5.99454355e+00]\n", - " [1.02999997e+00 2.33675718e+00]\n", - " [1.04999995e+00 1.23811281e+00]\n", - " [1.06999993e+00 1.01582968e+00]\n", - " [1.08999991e+00 8.55239391e-01]\n", - " [1.11000001e+00 7.35277355e-01]\n", - " [1.13000000e+00 6.13604367e-01]\n", - " [1.14999998e+00 5.92449248e-01]\n", - " [1.16999996e+00 5.50008595e-01]\n", - " [1.18999994e+00 5.10064900e-01]\n", - " [1.21000004e+00 4.51530993e-01]\n", - " [1.23000002e+00 4.89567488e-01]\n", - " [1.25000000e+00 4.27016288e-01]\n", - " [1.26999998e+00 4.59215075e-01]\n", - " [1.28999996e+00 5.03942072e-01]\n", - " [1.30999994e+00 4.35167283e-01]\n", - " [1.32999992e+00 3.94493490e-01]\n", - " [1.34999990e+00 5.10520041e-01]\n", - " [1.37000000e+00 5.54431260e-01]\n", - " [1.38999999e+00 5.19582152e-01]\n", - " [1.40999997e+00 5.11105120e-01]\n", - " [1.42999995e+00 5.29835641e-01]\n", - " [1.44999993e+00 5.91018200e-01]\n", - " [1.47000003e+00 6.06200635e-01]\n", - " [1.49000001e+00 6.86542988e-01]\n", - " [1.50999999e+00 6.92638040e-01]\n", - " [1.52999997e+00 6.27580464e-01]\n", - " [1.54999995e+00 6.98118865e-01]\n", - " [1.56999993e+00 7.89712429e-01]\n", - " [1.58999991e+00 7.62705624e-01]\n", - " [1.60999990e+00 9.11538541e-01]\n", - " [1.63000000e+00 1.03215218e+00]\n", - " [1.64999998e+00 1.16467202e+00]\n", - " [1.66999996e+00 1.46617603e+00]\n", - " [1.68999994e+00 1.50240195e+00]\n", - " [1.70999992e+00 1.47793388e+00]\n", - " [1.73000002e+00 1.47463214e+00]\n", - " [1.75000000e+00 1.33518815e+00]\n", - " [1.76999998e+00 1.34231174e+00]\n", - " [1.78999996e+00 1.24561477e+00]\n", - " [1.80999994e+00 1.21823835e+00]\n", - " [1.82999992e+00 1.19907069e+00]\n", - " [1.84999990e+00 1.19831800e+00]\n", - " [1.87000000e+00 1.33037162e+00]\n", - " [1.88999999e+00 1.32635713e+00]\n", - " [1.90999997e+00 1.49168754e+00]\n", - " [1.92999995e+00 1.49872899e+00]\n", - " [1.94999993e+00 1.51322103e+00]\n", - " [1.96999991e+00 1.38643694e+00]\n", - " [1.99000001e+00 1.11911285e+00]]\n" - ] - } - ], + "outputs": [], "source": [ "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", @@ -1386,1445 +134,15 @@ "plt.xlabel(\"$r$\")\n", "plt.ylabel(\"$g(r)$\")\n", "plt.show()\n", - "print(rdf_data)" + "#print(rdf_data)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "5f279783-91f5-4560-8d38-c8be52c09443", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Total steps 13200\n", - "DPD Energy Cutoff= 65.82111560065738\n" - ] - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " 2026-06-12T21:53:10.528353\n", - " image/svg+xml\n", - " \n", - " \n", - " Matplotlib v3.10.9, https://matplotlib.org/\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "\n" - ], - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "from dpd_utils import calculate_pair_energy\n", "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=npoly,num_mon=nmono,density=density)\n", @@ -2834,7 +152,7 @@ "print(\"Total steps\",len(pe)*100+100)\n", "print(\"DPD Energy Cutoff= \",U/N)\n", "\n", - "slice_idx = 100\n", + "slice_idx = 20\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", diff --git a/phantomwalk/examples/2-dpd-energy-gsd.ipynb b/phantomwalk/examples/2-dpd-energy-gsd.ipynb index bb5ebc0..515a004 100644 --- a/phantomwalk/examples/2-dpd-energy-gsd.ipynb +++ b/phantomwalk/examples/2-dpd-energy-gsd.ipynb @@ -11,30 +11,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "88429372", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", - "\n", - "****** PyMBAR will use 64-bit JAX! *******\n", - "* JAX is currently set to 32-bit bitsize *\n", - "* which is its default. *\n", - "* *\n", - "* PyMBAR requires 64-bit mode and WILL *\n", - "* enable JAX's 64-bit mode when called. *\n", - "* *\n", - "* This MAY cause problems with other *\n", - "* Uses of JAX in the same code. *\n", - "******************************************\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "import sys\n", "import os\n", @@ -60,21 +40,10 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "1be674a6-7412-4f5b-aad2-24cb3cc0f869", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Running with A=800, gamma=800, k=20000, num_pol=100, num_mon=100\n", - "Simulation timed out\n", - "Finished in time = 0.00s\n" - ] - } - ], + "outputs": [], "source": [ "num_pol=100\n", "num_mon=100\n", @@ -101,18 +70,10 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "47e9b287-e7bd-4fbb-bfb4-b27c59206eb1", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "124469.44818043037\n" - ] - } - ], + "outputs": [], "source": [ "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=num_pol,num_mon=num_mon,density=density)\n", "print(U)" @@ -128,1401 +89,10 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "id": "1b04fd76-9337-421e-bcfc-2b351c73e700", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Total steps 104200\n", - "DPD Energy Cutoff= 12.446944818043036\n" - ] - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " 2026-06-12T18:11:55.976843\n", - " image/svg+xml\n", - " \n", - " \n", - " Matplotlib v3.10.9, https://matplotlib.org/\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "\n" - ], - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "log = np.genfromtxt(\"log.txt\", names=True)\n", "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", @@ -1576,1247 +146,10 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "id": "133f7da1-415e-49d6-9731-b7f642a62469", "metadata": {}, - "outputs": [ - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " 2026-06-12T21:33:21.305317\n", - " image/svg+xml\n", - " \n", - " \n", - " Matplotlib v3.10.9, https://matplotlib.org/\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "\n" - ], - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[9.99999978e-03 0.00000000e+00]\n", - " [2.99999993e-02 0.00000000e+00]\n", - " [4.99999970e-02 0.00000000e+00]\n", - " [7.00000003e-02 0.00000000e+00]\n", - " [8.99999961e-02 0.00000000e+00]\n", - " [1.09999999e-01 0.00000000e+00]\n", - " [1.29999995e-01 0.00000000e+00]\n", - " [1.50000006e-01 0.00000000e+00]\n", - " [1.69999987e-01 0.00000000e+00]\n", - " [1.89999998e-01 0.00000000e+00]\n", - " [2.09999993e-01 0.00000000e+00]\n", - " [2.29999989e-01 0.00000000e+00]\n", - " [2.50000000e-01 0.00000000e+00]\n", - " [2.69999981e-01 0.00000000e+00]\n", - " [2.89999992e-01 0.00000000e+00]\n", - " [3.10000002e-01 0.00000000e+00]\n", - " [3.29999983e-01 0.00000000e+00]\n", - " [3.49999994e-01 0.00000000e+00]\n", - " [3.70000005e-01 0.00000000e+00]\n", - " [3.89999986e-01 0.00000000e+00]\n", - " [4.09999967e-01 0.00000000e+00]\n", - " [4.30000007e-01 0.00000000e+00]\n", - " [4.49999988e-01 0.00000000e+00]\n", - " [4.69999969e-01 0.00000000e+00]\n", - " [4.90000010e-01 0.00000000e+00]\n", - " [5.09999990e-01 0.00000000e+00]\n", - " [5.29999971e-01 0.00000000e+00]\n", - " [5.49999952e-01 0.00000000e+00]\n", - " [5.69999993e-01 0.00000000e+00]\n", - " [5.89999974e-01 0.00000000e+00]\n", - " [6.10000014e-01 0.00000000e+00]\n", - " [6.29999995e-01 0.00000000e+00]\n", - " [6.49999976e-01 0.00000000e+00]\n", - " [6.69999957e-01 0.00000000e+00]\n", - " [6.89999998e-01 0.00000000e+00]\n", - " [7.09999979e-01 0.00000000e+00]\n", - " [7.30000019e-01 0.00000000e+00]\n", - " [7.50000000e-01 0.00000000e+00]\n", - " [7.69999981e-01 0.00000000e+00]\n", - " [7.89999962e-01 0.00000000e+00]\n", - " [8.09999943e-01 0.00000000e+00]\n", - " [8.29999983e-01 0.00000000e+00]\n", - " [8.49999964e-01 0.00000000e+00]\n", - " [8.70000005e-01 1.61740478e-04]\n", - " [8.89999986e-01 1.85464346e-03]\n", - " [9.09999967e-01 3.42237689e-02]\n", - " [9.29999948e-01 3.75659347e-01]\n", - " [9.49999988e-01 2.12913251e+00]\n", - " [9.69999969e-01 6.38623619e+00]\n", - " [9.90000010e-01 9.32189941e+00]\n", - " [1.00999999e+00 5.86049223e+00]\n", - " [1.02999997e+00 2.40495610e+00]\n", - " [1.04999995e+00 1.30290556e+00]\n", - " [1.06999993e+00 9.48677957e-01]\n", - " [1.08999991e+00 7.72342980e-01]\n", - " [1.11000001e+00 6.71785176e-01]\n", - " [1.13000000e+00 6.23863041e-01]\n", - " [1.14999998e+00 5.86108208e-01]\n", - " [1.16999996e+00 5.55777013e-01]\n", - " [1.18999994e+00 5.20482302e-01]\n", - " [1.21000004e+00 5.19804120e-01]\n", - " [1.23000002e+00 4.97619092e-01]\n", - " [1.25000000e+00 4.86955285e-01]\n", - " [1.26999998e+00 4.79405373e-01]\n", - " [1.28999996e+00 4.77862120e-01]\n", - " [1.30999994e+00 4.69837993e-01]\n", - " [1.32999992e+00 4.74741757e-01]\n", - " [1.34999990e+00 4.84658182e-01]\n", - " [1.37000000e+00 4.91095632e-01]\n", - " [1.38999999e+00 4.97499943e-01]\n", - " [1.40999997e+00 5.17201483e-01]\n", - " [1.42999995e+00 5.36301434e-01]\n", - " [1.44999993e+00 5.55993795e-01]\n", - " [1.47000003e+00 5.77137053e-01]\n", - " [1.49000001e+00 6.08183384e-01]\n", - " [1.50999999e+00 6.34086013e-01]\n", - " [1.52999997e+00 6.68373168e-01]\n", - " [1.54999995e+00 7.17355371e-01]\n", - " [1.56999993e+00 7.84075201e-01]\n", - " [1.58999991e+00 8.83406878e-01]\n", - " [1.60999990e+00 9.75062847e-01]\n", - " [1.63000000e+00 1.07528138e+00]\n", - " [1.64999998e+00 1.19554257e+00]\n", - " [1.66999996e+00 1.33810508e+00]\n", - " [1.68999994e+00 1.46330953e+00]\n", - " [1.70999992e+00 1.51521719e+00]\n", - " [1.73000002e+00 1.46107221e+00]\n", - " [1.75000000e+00 1.35471618e+00]\n", - " [1.76999998e+00 1.29108131e+00]\n", - " [1.78999996e+00 1.24045670e+00]\n", - " [1.80999994e+00 1.21756566e+00]\n", - " [1.82999992e+00 1.23899102e+00]\n", - " [1.84999990e+00 1.25472844e+00]\n", - " [1.87000000e+00 1.29916048e+00]\n", - " [1.88999999e+00 1.36729598e+00]\n", - " [1.90999997e+00 1.43487251e+00]\n", - " [1.92999995e+00 1.48630548e+00]\n", - " [1.94999993e+00 1.46616626e+00]\n", - " [1.96999991e+00 1.35239911e+00]\n", - " [1.99000001e+00 1.14362824e+00]]\n" - ] - } - ], + "outputs": [], "source": [ "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", From 0da18203a8ce24e1b7c3c608903d640b53632397 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Sat, 13 Jun 2026 11:40:51 -0600 Subject: [PATCH 08/18] Fiddling with 1.4 densities - can run acceptably fast but have not been able to continue LJ system yet. --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 12 +-- phantomwalk/examples/3-dpd-to-lj-wca.ipynb | 89 +++++++++++++++----- 2 files changed, 77 insertions(+), 24 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 3c7134e..6065423 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -83,10 +83,12 @@ "A=50000\n", "k=50000\n", "dt=0.001\n", + "#r_cut=1.01\n", + "#min_pair_dist=0.988\n", "r_cut=1.01\n", - "min_pair_dist=0.988\n", + "min_pair_dist=.97\n", "gamma = 1200\n", - "density=1.3\n", + "density=1.4\n", "\n", "for i in range(nruns):\n", " seed = np.random.randint(50000)\n", @@ -102,9 +104,9 @@ " energy=True,\n", " write=True,\n", " gsd_file_name='trajectory.gsd',\n", - " gsd_write_freq=100,\n", + " gsd_write_freq=20,\n", " log_file_name='log.txt',\n", - " log_write_freq=10,\n", + " log_write_freq=1,\n", " density=density,\n", " \n", " dt=dt,\n", @@ -152,7 +154,7 @@ "print(\"Total steps\",len(pe)*100+100)\n", "print(\"DPD Energy Cutoff= \",U/N)\n", "\n", - "slice_idx = 20\n", + "slice_idx = 40\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", diff --git a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb index 7033544..f1ba1c6 100644 --- a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb +++ b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb @@ -11,10 +11,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "88429372", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", + "\n", + "****** PyMBAR will use 64-bit JAX! *******\n", + "* JAX is currently set to 32-bit bitsize *\n", + "* which is its default. *\n", + "* *\n", + "* PyMBAR requires 64-bit mode and WILL *\n", + "* enable JAX's 64-bit mode when called. *\n", + "* *\n", + "* This MAY cause problems with other *\n", + "* Uses of JAX in the same code. *\n", + "******************************************\n", + "\n" + ] + } + ], "source": [ "import sys\n", "import os\n", @@ -38,28 +58,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "id": "1be674a6-7412-4f5b-aad2-24cb3cc0f869", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Finished in time = 0.30s\n" + ] + } + ], "source": [ - "num_pol=100\n", + "num_pol=10\n", "num_mon=100\n", "N = num_pol*num_mon\n", "density=0.8\n", - "A=1000\n", - "r_cut=1.1\n", - "min_pair_dist=1.05\n", + "A=50000\n", + "r_cut=1.2\n", + "min_pair_dist=1.1\n", + "seed = 122\n", "dpd_final_frame, s = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", " num_mon=num_mon,\n", " density=density,\n", - " k=20000,\n", + " k=50000,\n", " r_cut=r_cut,\n", " A=A,\n", - " gamma=1000,\n", - " sim_seed=1234,\n", - " np_seed=1234,\n", + " gamma=1200,\n", + " sim_seed=seed,\n", + " np_seed=seed,\n", " energy=True,\n", " min_pair_dist=min_pair_dist,\n", ")\n", @@ -76,20 +105,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "id": "6139f233-10e2-4288-8cc0-cc78717d9f35", "metadata": {}, "outputs": [], "source": [ "def run_lj_simulation(\n", " dpd_final_frame,\n", - " random_seed=24,\n", - " dt=0.001,\n", + " random_seed=25,\n", + " dt=0.0005,\n", " lj_epsilon=1.0,\n", " lj_sigma=1.0,\n", " lj_r_cut=1.2,\n", " fene_k=30,\n", - " fene_r0=1.05,\n", + " fene_r0=1.,\n", " fene_epsilon=1.0,\n", " fene_sigma=1.0,\n", " fene_delta=0,\n", @@ -234,10 +263,32 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 25, "id": "48f19fa1-6ec6-4244-a23c-fd9d774da3f1", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "**ERROR**: bond.fene: bond out of bounds\n", + "\n" + ] + }, + { + "ename": "RuntimeError", + "evalue": "Error in bond calculation", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mRuntimeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[25]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m run_lj_simulation(dpd_final_frame=dpd_final_frame)\n", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[24]\u001b[39m\u001b[32m, line 140\u001b[39m, in \u001b[36mrun_lj_simulation\u001b[39m\u001b[34m(dpd_final_frame, random_seed, dt, lj_epsilon, lj_sigma, lj_r_cut, fene_k, fene_r0, fene_epsilon, fene_sigma, fene_delta, angle_k, angle_t0, dihedral_k, dihedral_d, dihedral_n, dihedral_phi0)\u001b[39m\n\u001b[32m 136\u001b[39m add_hoomd_writers(sim=LJ_sim)\n\u001b[32m 137\u001b[39m \n\u001b[32m 138\u001b[39m \u001b[38;5;66;03m# Run short equilibration\u001b[39;00m\n\u001b[32m 139\u001b[39m LJ_sim.run(\u001b[32m0\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m140\u001b[39m LJ_sim.run(\u001b[32m100\u001b[39m)\n\u001b[32m 141\u001b[39m \n\u001b[32m 142\u001b[39m \u001b[38;5;66;03m# Flush outputs\u001b[39;00m\n\u001b[32m 143\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m writer \u001b[38;5;28;01min\u001b[39;00m LJ_sim.operations.writers:\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/miniforge3/envs/phantomwalk/lib/python3.12/site-packages/hoomd/simulation.py:561\u001b[39m, in \u001b[36mSimulation.run\u001b[39m\u001b[34m(self, steps, write_at_start)\u001b[39m\n\u001b[32m 558\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m steps_int < \u001b[32m0\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m steps_int > TIMESTEP_MAX - \u001b[32m1\u001b[39m:\n\u001b[32m 559\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33msteps must be in the range [0, \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mTIMESTEP_MAX\u001b[38;5;250m \u001b[39m-\u001b[38;5;250m \u001b[39m\u001b[32m1\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m]\u001b[39m\u001b[33m\"\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m561\u001b[39m \u001b[30;43mself\u001b[39;49m\u001b[30;43m.\u001b[39;49m\u001b[30;43m_cpp_sys\u001b[39;49m\u001b[30;43m.\u001b[39;49m\u001b[30;43mrun\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43msteps_int\u001b[39;49m\u001b[30;43m,\u001b[39;49m\u001b[30;43m \u001b[39;49m\u001b[30;43mwrite_at_start\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n", + "\u001b[31mRuntimeError\u001b[39m: Error in bond calculation" + ] + } + ], "source": [ "run_lj_simulation(dpd_final_frame=dpd_final_frame)" ] @@ -267,7 +318,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.13" + "version": "3.12.13" } }, "nbformat": 4, From 313cb790c7913eaf4cecf3b53cd2381dfb008986 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Mon, 15 Jun 2026 08:11:46 -0600 Subject: [PATCH 09/18] More minimal .gitignore --- .gitignore | 189 +---------------------------------------------------- 1 file changed, 1 insertion(+), 188 deletions(-) diff --git a/.gitignore b/.gitignore index 6cf64c2..042f5a6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ # Output files from MD runs *.gsd *.txt +notes # Byte-compiled / optimized / DLL files __pycache__/ @@ -11,71 +12,6 @@ phantomwalk/examples/rdf.csv phantomwalk/examples/trajectory.gsd phantomwalk/examples/log.txt -# C extensions -*.so - -# Distribution / packaging -dist/ -downloads/ -eggs/ -.eggs/ -lib64/ -parts/ -sdist/ -var/ -wheels/ -share/python-wheels/ -*.egg-info/ - -# PyInstaller -# Usually these files are written by a python script from a template -# before PyInstaller builds the exe, so as to inject date/other infos into it. -*.manifest -*.spec - -# Installer logs -pip-log.txt -pip-delete-this-directory.txt - -# Unit test / coverage reports -htmlcov/ -.tox/ -.nox/ -.coverage -.coverage.* -.cache -nosetests.xml -coverage.xml -*.cover -*.py.cover -.hypothesis/ -.pytest_cache/ -cover/ - -# Translations -*.mo -*.pot - -# Django stuff: -*.log -local_settings.py -db.sqlite3 -db.sqlite3-journal - -# Flask stuff: -instance/ -.webassets-cache - -# Scrapy stuff: -.scrapy - -# Sphinx documentation -docs/_build/ - -# PyBuilder -.pybuilder/ -target/ - # Jupyter Notebook .ipynb_checkpoints @@ -83,126 +19,3 @@ target/ profile_default/ ipython_config.py -# pyenv -# For a library or package, you might want to ignore these files since the code is -# intended to run in multiple environments; otherwise, check them in: -# .python-version - -# pipenv -# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. -# However, in case of collaboration, if having platform-specific dependencies or dependencies -# having no cross-platform support, pipenv may install dependencies that don't work, or not -# install all needed dependencies. -#Pipfile.lock - -# UV -# Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control. -# This is especially recommended for binary packages to ensure reproducibility, and is more -# commonly ignored for libraries. -#uv.lock - -# poetry -# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. -# This is especially recommended for binary packages to ensure reproducibility, and is more -# commonly ignored for libraries. -# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control -#poetry.lock -#poetry.toml - -# pdm -# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. -# pdm recommends including project-wide configuration in pdm.toml, but excluding .pdm-python. -# https://pdm-project.org/en/latest/usage/project/#working-with-version-control -#pdm.lock -#pdm.toml -.pdm-python -.pdm-build/ - -# pixi -# Similar to Pipfile.lock, it is generally recommended to include pixi.lock in version control. -#pixi.lock -# Pixi creates a virtual environment in the .pixi directory, just like venv module creates one -# in the .venv directory. It is recommended not to include this directory in version control. -.pixi - -# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm -__pypackages__/ - -# Celery stuff -celerybeat-schedule -celerybeat.pid - -# SageMath parsed files -*.sage.py - -# Environments -.env -.envrc -.venv -env/ -venv/ -ENV/ -env.bak/ -venv.bak/ - -# Spyder project settings -.spyderproject -.spyproject - -# Rope project settings -.ropeproject - -# mkdocs documentation -/site - -# mypy -.mypy_cache/ -.dmypy.json -dmypy.json - -# Pyre type checker -.pyre/ - -# pytype static type analyzer -.pytype/ - -# Cython debug symbols -cython_debug/ - -# PyCharm -# JetBrains specific template is maintained in a separate JetBrains.gitignore that can -# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore -# and can be added to the global gitignore or merged into this file. For a more nuclear -# option (not recommended) you can uncomment the following to ignore the entire idea folder. -#.idea/ - -# Abstra -# Abstra is an AI-powered process automation framework. -# Ignore directories containing user credentials, local state, and settings. -# Learn more at https://abstra.io/docs -.abstra/ - -# Visual Studio Code -# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore -# that can be found at https://github.com/github/gitignore/blob/main/Global/VisualStudioCode.gitignore -# and can be added to the global gitignore or merged into this file. However, if you prefer, -# you could uncomment the following to ignore the entire vscode folder -# .vscode/ - -# Ruff stuff: -.ruff_cache/ - -# PyPI configuration file -.pypirc - -# Cursor -# Cursor is an AI-powered code editor. `.cursorignore` specifies files/directories to -# exclude from AI features like autocomplete and code analysis. Recommended for sensitive data -# refer to https://docs.cursor.com/context/ignore-files -.cursorignore -.cursorindexingignore - -# Marimo -marimo/_static/ -marimo/_lsp/ -__marimo__/ From 728c36757ffadd9272704171eeb6a59d540aff1c Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Mon, 15 Jun 2026 09:32:51 -0600 Subject: [PATCH 10/18] ignore vim swap files --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 042f5a6..1a4f566 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ # Output files from MD runs *.gsd *.txt +*.swp notes # Byte-compiled / optimized / DLL files From e86ba8d6d55622221b6313b2e571c422abdb2618 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Mon, 15 Jun 2026 14:04:49 -0600 Subject: [PATCH 11/18] Broke things. --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 18 ++-- phantomwalk/lib/create_system_dpd.py | 91 ++++++++++---------- phantomwalk/lib/dpd_utils.py | 74 +--------------- 3 files changed, 55 insertions(+), 128 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 6065423..8a235c1 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -83,12 +83,10 @@ "A=50000\n", "k=50000\n", "dt=0.001\n", - "#r_cut=1.01\n", - "#min_pair_dist=0.988\n", - "r_cut=1.01\n", - "min_pair_dist=.97\n", + "r_cut=1.15\n", + "min_pair_dist=.90\n", "gamma = 1200\n", - "density=1.4\n", + "density=0.2\n", "\n", "for i in range(nruns):\n", " seed = np.random.randint(50000)\n", @@ -100,13 +98,11 @@ " sim_seed=seed,\n", " np_seed=seed,\n", " sim_steps_incr=100,\n", - " loop_timeout=600,\n", - " energy=True,\n", - " write=True,\n", + " loop_timeout=60,\n", " gsd_file_name='trajectory.gsd',\n", " gsd_write_freq=20,\n", " log_file_name='log.txt',\n", - " log_write_freq=1,\n", + " log_write_freq=10,\n", " density=density,\n", " \n", " dt=dt,\n", @@ -114,7 +110,9 @@ " min_pair_dist=min_pair_dist,\n", " A=A,\n", " k=k,\n", - " gamma=gamma)\n", + " gamma=gamma,\n", + " energy_scaling = 1.,\n", + " )\n", " rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", " b = (rdf_data[:,1] !=0).argmax()\n", " times.append(s)\n", diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index 7f30f32..cf79d1b 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -7,6 +7,10 @@ from dpd_utils import initialize_snapshot_rand_walk,check_bond_length_equilibration,check_inter_particle_distance,add_hoomd_writers,check_pair_energy,simulation_energy_end +def get_close(rdf): + b =(rdf.rdf[:,1] !=0).argmax() + return rdf.bin_centers[b][0] + def create_polymer_system_dpd( num_pol, num_mon, @@ -24,13 +28,14 @@ def create_polymer_system_dpd( loop_timeout=60, energy=True, min_pair_dist=1.05, + energy_scaling= 5, + bond_tolerance = 0.05, write=True, gsd_file_name='trajectory.gsd', gsd_write_freq=10, log_file_name='log.txt', log_write_freq=10 ): - ''' Initialize a polymer system in a cubic box using a random walk and a HOOMD simulation with DPD forces. @@ -90,11 +95,7 @@ def create_polymer_system_dpd( execution time of the DPD workflow, build + simulation wall time ''' - #print(num_pol*num_mon) - #print(f"\nRunning with A={A}, gamma={gamma}, k={k}, " - # f"num_pol={num_pol}, num_mon={num_mon}") start_time = time.perf_counter() - frame = initialize_snapshot_rand_walk( num_mon=num_mon, num_pol=num_pol, @@ -103,8 +104,6 @@ def create_polymer_system_dpd( seed=np_seed ) - build_stop = time.perf_counter() - #print("Total build time: ", build_stop-start_time) harmonic = hoomd.md.bond.Harmonic() harmonic.params["b"] = dict(r0=bond_l, k=k) integrator = hoomd.md.Integrator(dt=dt) @@ -117,58 +116,58 @@ def create_polymer_system_dpd( nlist = hoomd.md.nlist.Cell(buffer=0.4,exclusions=['bond']) simulation.operations.nlist = nlist DPD = hoomd.md.pair.DPD(nlist, default_r_cut=r_cut, kT=kT) + DPDc = hoomd.md.pair.DPDConservative(nlist, default_r_cut=r_cut) DPD.params[('A', 'A')] = dict(A=A, gamma=gamma) + DPDc.params[('A', 'A')] = dict(A=A) integrator.forces.append(DPD) + + N = num_mon*num_pol + maxPerParticle = A*( (min_pair_dist*min_pair_dist)/(2*r_cut) - min_pair_dist + r_cut/2) + maxPerParticle *= density*density*energy_scaling + maxPerBond = k*bond_tolerance*bond_tolerance/2 + print("max bond = {:.2f}, max per Particle = {:.2f}".format(maxPerBond, maxPerParticle)) if write: - rdf = add_hoomd_writers( - simulation, - gsd_file_name, - gsd_write_freq, - log_file_name, - log_write_freq - ) - + rdf,thermo = add_hoomd_writers( simulation, gsd_file_name, gsd_write_freq, log_file_name,log_write_freq ) simulation.run(1) for writer in simulation.operations.writers: if hasattr(writer, "flush"): writer.flush() - if energy: - while not simulation_energy_end( - A=A, - r=min_pair_dist, - r_cut=r_cut, - num_pol=num_pol, - num_mon=num_mon, - density=density, - log_file_name=log_file_name, - ): - check_time = time.perf_counter() - if (check_time-start_time) > loop_timeout: - print("Simulation timed out") - return simulation.state.get_snapshot(), 0 - simulation.run(sim_steps_incr) - for writer in simulation.operations.writers: - if hasattr(writer, "flush"): - writer.flush() - else: - while not check_inter_particle_distance(snap,minimum_distance=min_pair_dist): - check_time = time.perf_counter() - if (check_time-start_time) > loop_timeout: - print("Simulation timed out") - return snap,0 - simulation.run(sim_steps_incr) - for writer in simulation.operations.writers: - if hasattr(writer, "flush"): - writer.flush() - snap=simulation.state.get_snapshot() + while DPDc.energy/N > maxPerParticle: + check_time = time.perf_counter() + if (check_time-start_time) > loop_timeout: + print("Simulation timed out in energy") + return simulation.state.get_snapshot(), loop_timeout + simulation.run(sim_steps_incr) + for writer in simulation.operations.writers: + if hasattr(writer, "flush"): + writer.flush() + + while harmonic.energy/N > maxPerBond:#TODO normalize + check_time = time.perf_counter() + if (check_time-start_time) > loop_timeout: + print("Simulation timed out in bond energy") + return simulation.state.get_snapshot(), loop_timeout + simulation.run(sim_steps_incr) + for writer in simulation.operations.writers: + if hasattr(writer, "flush"): + writer.flush() + + while get_close(rdf) < min_pair_dist:#TODO RDF + check_time = time.perf_counter() + if (check_time-start_time) > loop_timeout: + print("Simulation timed out in rdf polish") + return simulation.state.get_snapshot(), loop_timeout + simulation.run(sim_steps_incr) + for writer in simulation.operations.writers: + if hasattr(writer, "flush"): + writer.flush() end_time = time.perf_counter() total_time = end_time - start_time - #print("Total build and simulation time:", end_time - start_time) + print("Total build and simulation time:", end_time - start_time) np.savetxt( "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)" ) - #return snap, total_time return simulation.state.get_snapshot(), total_time diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 4a70606..a5e1540 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -5,101 +5,33 @@ import time from cmeutils.sampling import is_equilibrated - - -#https://github.com/joelaforet/mupt/blob/issue-77-aa-dpd-builder/mupt/builders/all_atom_dpd.py -class _ParameterTables: - """HOOMD-ready bonded and vdW parameter tables.""" - - bond_params: dict[str, dict[str, float]] = field(default_factory=dict) - angle_params: dict[str, dict[str, float]] = field(default_factory=dict) - dihedral_params: dict[str, dict[str, float]] = field(default_factory=dict) - improper_params: dict[str, dict[str, float]] = field(default_factory=dict) - bond_type_by_group: dict[tuple[int, int], str] = field(default_factory=dict) - angle_type_by_group: dict[tuple[int, int, int], str] = field(default_factory=dict) - dihedral_type_by_group: dict[tuple[int, int, int, int], list[str]] = field(default_factory=dict) - improper_type_by_group: dict[tuple[int, int, int, int], list[str]] = field(default_factory=dict) - atom_epsilons: dict[int, float] = field(default_factory=dict) - atom_types_by_global: dict[int, str] = field(default_factory=dict) - epsilon_by_type: dict[str, float] = field(default_factory=dict) - -def initialize_from_topology(bond_list, density, box, bond_length=1.0, seed=1234): - ''' - Currently expects bond_list to be list of ordered tuples with indices from 0 to N-1, each tuple describing a bond - - ''' - rng = np.random.default_rng(seed) - N = np.max(bond_list[:,1])+1 #assumes ordered tuples - - #Overall algorithm: - # give all particles initial random coords - positions = rng.uniform(0, np.max(box), size=(N, 3)) - # Generate a uniform-sphere-delta for each bond - # Cumulative sum of deltas gives positions for a chain off of a branch (or origin) - thetas = rng.uniform(0,2*np.pi,size=N) - phis = np.arccos(rng.uniform(-1,1,size=N) - x = np.sin(phis)*np.cos(thetas) - y = np.sin(phis)*np.sin(thetas) - z = np.cos(phis) - deltas = np.stack([x,y,z],axis=2) * bond_length - displacements = np.cumsum(deltas, axis=1) - #positions_view[:, 1:, :] = starts[:, None, :] + displacements #deprecated - #TODO: loop over bonds and add displacements here - - #pbc - #TODO Wrap coordinates along each axis separately - positions %= L - positions -= L/2 - - #TODO: Fold in PR 83 mupt, mupt.builder.all_atom_dpd functionality. - - frame = gsd.hoomd.Frame() - frame.particles.types = ['A'] #update with all types - frame.particles.N = N - frame.particles.position = positions - frame.bonds.N = len(bonds_list) - frame.bonds.group = bond_list - frame.bonds.types = ['b'] #update with all types - frame.configuration.box = box - - return frame - def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, seed=1234): ''' Create a HOOMD snapshot of a cubic box with the number density given by input parameters. Configure particles using a random walk. ''' rng = np.random.default_rng(seed) - N = num_pol * num_mon L = np.cbrt(N / density) - positions = np.empty((N, 3)) starts = rng.uniform(0, L, size=(num_pol, 3)) - thetas = rng.uniform(0,2*np.pi,size=(num_pol,num_mon-1)) phis = np.arccos(rng.uniform(-1,1,size=(num_pol,num_mon-1))) x = np.sin(phis)*np.cos(thetas) y = np.sin(phis)*np.sin(thetas) z = np.cos(phis) - deltas = np.stack([x,y,z],axis=2) * bond_length displacements = np.cumsum(deltas, axis=1) - positions_view = positions.reshape(num_pol, num_mon, 3) positions_view[:, 0, :] = starts positions_view[:, 1:, :] = starts[:, None, :] + displacements - - #pbc positions %= L - positions -= L/2 - + positions -= L/2 #TODO: use box in flowerMD indices = np.arange(N).reshape(num_pol, num_mon) bonds = np.column_stack([ indices[:, :-1].ravel(), indices[:, 1:].ravel() ]) - frame = gsd.hoomd.Frame() frame.particles.types = ['A'] frame.particles.N = N @@ -107,8 +39,6 @@ def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, se frame.bonds.N = len(bonds) frame.bonds.group = bonds frame.bonds.types = ['b'] - frame.configuration.box = [L, L, L, 0, 0, 0] - return frame def check_bond_length_equilibration(snap, num_mon, num_pol, max_bond_length=1.1, min_bond_length=0.95): @@ -255,7 +185,7 @@ def act(self, timestep): sim.operations.writers.append(gsd_writer) sim.operations.writers.append(table_file) - return rdf + return rdf, thermo_props def check_pair_energy(energy_idx=-1, log_file_name="log.txt"): """Check whether the pair interaction energy has equilibrated. From 30fd9a62911f8b55490bd51ae058cd92f6db67d5 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Mon, 15 Jun 2026 15:11:03 -0600 Subject: [PATCH 12/18] Updated Stop Criteria! --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 25 +++--- phantomwalk/lib/create_system_dpd.py | 30 +++---- phantomwalk/lib/dpd_utils.py | 92 +------------------- 3 files changed, 25 insertions(+), 122 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 8a235c1..7691510 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -77,20 +77,21 @@ "times = []\n", "rs = []\n", "\n", - "npoly=10\n", - "nmono=100\n", + "density=1.1\n", + "npoly=100\n", + "nmono=10\n", "N=npoly*nmono\n", "A=50000\n", "k=50000\n", "dt=0.001\n", - "r_cut=1.15\n", - "min_pair_dist=.90\n", "gamma = 1200\n", - "density=0.2\n", + "r_cut=1.01\n", + "min_pair_dist=.95\n", + "es = 0.05\n", "\n", "for i in range(nruns):\n", " seed = np.random.randint(50000)\n", - " last_frame, s = create_polymer_system_dpd(\n", + " last_frame, closest, s = create_polymer_system_dpd(\n", " bond_l=1.0,\n", " num_pol=npoly,\n", " num_mon=nmono,\n", @@ -111,13 +112,11 @@ " A=A,\n", " k=k,\n", " gamma=gamma,\n", - " energy_scaling = 1.,\n", + " energy_scaling = es,\n", " )\n", - " rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", - " b = (rdf_data[:,1] !=0).argmax()\n", " times.append(s)\n", - " rs.append(rdf_data[b][0])\n", - " print(\"{:.2f}s, r = {:.2f} g(r) = {:.4f}\".format(s,rdf_data[b][0],rdf_data[b][1]))\n", + " rs.append(closest)\n", + " print(\"{:.2f}s, r = {:.2f} \".format(s,closest))\n", "print(\"\\n\\n{:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(np.average(times),np.std(times), np.average(rs), np.std(rs)))\n" ] }, @@ -144,20 +143,16 @@ "metadata": {}, "outputs": [], "source": [ - "from dpd_utils import calculate_pair_energy\n", - "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=npoly,num_mon=nmono,density=density)\n", "log = np.genfromtxt(\"log.txt\", names=True)\n", "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", "pairs = log[\"mdpairDPDenergy\"]\n", "print(\"Total steps\",len(pe)*100+100)\n", - "print(\"DPD Energy Cutoff= \",U/N)\n", "\n", "slice_idx = 40\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", - "plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", "plt.title(\"Energy vs Time\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy/N\")\n", diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index cf79d1b..b5d7499 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -1,15 +1,13 @@ import numpy as np -import freud -import gsd, gsd.hoomd import hoomd import time -from dpd_utils import initialize_snapshot_rand_walk,check_bond_length_equilibration,check_inter_particle_distance,add_hoomd_writers,check_pair_energy,simulation_energy_end +from dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers def get_close(rdf): - b =(rdf.rdf[:,1] !=0).argmax() - return rdf.bin_centers[b][0] + b =(rdf.rdf !=0).argmax() + return rdf.bin_centers[b] def create_polymer_system_dpd( num_pol, @@ -104,6 +102,7 @@ def create_polymer_system_dpd( seed=np_seed ) + build_stop = time.perf_counter() harmonic = hoomd.md.bond.Harmonic() harmonic.params["b"] = dict(r0=bond_l, k=k) integrator = hoomd.md.Integrator(dt=dt) @@ -116,25 +115,25 @@ def create_polymer_system_dpd( nlist = hoomd.md.nlist.Cell(buffer=0.4,exclusions=['bond']) simulation.operations.nlist = nlist DPD = hoomd.md.pair.DPD(nlist, default_r_cut=r_cut, kT=kT) - DPDc = hoomd.md.pair.DPDConservative(nlist, default_r_cut=r_cut) DPD.params[('A', 'A')] = dict(A=A, gamma=gamma) - DPDc.params[('A', 'A')] = dict(A=A) integrator.forces.append(DPD) N = num_mon*num_pol maxPerParticle = A*( (min_pair_dist*min_pair_dist)/(2*r_cut) - min_pair_dist + r_cut/2) maxPerParticle *= density*density*energy_scaling maxPerBond = k*bond_tolerance*bond_tolerance/2 - print("max bond = {:.2f}, max per Particle = {:.2f}".format(maxPerBond, maxPerParticle)) + print("max per particle= {:.2f}, max per bond= {:.2f}".format(maxPerParticle, maxPerBond)) if write: rdf,thermo = add_hoomd_writers( simulation, gsd_file_name, gsd_write_freq, log_file_name,log_write_freq ) + simulation.run(1) + for writer in simulation.operations.writers: if hasattr(writer, "flush"): writer.flush() - while DPDc.energy/N > maxPerParticle: + while DPD.energy/N > maxPerParticle: check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in energy") @@ -144,7 +143,7 @@ def create_polymer_system_dpd( if hasattr(writer, "flush"): writer.flush() - while harmonic.energy/N > maxPerBond:#TODO normalize + while harmonic.energy/frame.bonds.N > maxPerBond: check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in bond energy") @@ -154,20 +153,19 @@ def create_polymer_system_dpd( if hasattr(writer, "flush"): writer.flush() - while get_close(rdf) < min_pair_dist:#TODO RDF + closest = get_close(rdf) + while closest < min_pair_dist: check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in rdf polish") return simulation.state.get_snapshot(), loop_timeout simulation.run(sim_steps_incr) + closest = get_close(rdf) for writer in simulation.operations.writers: if hasattr(writer, "flush"): writer.flush() end_time = time.perf_counter() total_time = end_time - start_time - print("Total build and simulation time:", end_time - start_time) - np.savetxt( - "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)" - ) - return simulation.state.get_snapshot(), total_time + np.savetxt( "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)") + return simulation.state.get_snapshot(), closest, total_time diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index a5e1540..3be7077 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -2,8 +2,6 @@ import freud import gsd, gsd.hoomd import hoomd -import time -from cmeutils.sampling import is_equilibrated def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, seed=1234): ''' @@ -39,49 +37,9 @@ def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, se frame.bonds.N = len(bonds) frame.bonds.group = bonds frame.bonds.types = ['b'] + frame.configuration.box = [L,L,L,0,0,0] return frame -def check_bond_length_equilibration(snap, num_mon, num_pol, max_bond_length=1.1, min_bond_length=0.95): - ''' - Check the bond distances. - - ''' - frame_ds = [] - for j in range(num_pol): - idx = j*num_mon - d1 = snap.particles.position[idx:idx+num_mon-1] - snap.particles.position[idx+1:idx+num_mon] - L = snap.configuration.box[0] - d1 -= L*np.round(d1/L) - bond_l = np.linalg.norm(d1,axis=1) - frame_ds.append(bond_l) - max_frame_bond_l = np.max(np.array(frame_ds)) - min_frame_bond_l = np.min(np.array(frame_ds)) - print("max: ",max_frame_bond_l," min: ",min_frame_bond_l) - if max_frame_bond_l <= max_bond_length and min_frame_bond_l >= min_bond_length: - print("Bonds relaxed.") - return True - if max_frame_bond_l > max_bond_length or min_frame_bond_l < min_bond_length: - return False - -def check_inter_particle_distance(snap, minimum_distance=0.95): - ''' - Check particle separations. - - ''' - positions = snap.particles.position - box = snap.configuration.box - aq = freud.locality.AABBQuery(box,positions) - aq_query = aq.query( - query_points=positions, - query_args=dict(r_min=0.0, r_max=minimum_distance, exclude_ii=True), - ) - nlist = aq_query.toNeighborList() - if len(nlist)==0: - print("Inter-particle separation reached.") - return True - else: - return False - def add_hoomd_writers( sim, gsd_file_name="trajectory.gsd", @@ -187,51 +145,3 @@ def act(self, timestep): sim.operations.writers.append(table_file) return rdf, thermo_props -def check_pair_energy(energy_idx=-1, log_file_name="log.txt"): - """Check whether the pair interaction energy has equilibrated. - - Pair energies are read from the HOOMD log file and analyzed - using pymbar timeseries equilibration detection. - - Parameters - ---------- - energy_idx : int, default -1 - Number of initial simulation steps to discard before - performing equilibration analysis. Default is to return the last frame. - - Returns - ------- - float, energy of last frame(s) of dpd simulation - - """ - log = np.genfromtxt(log_file_name, names=True) - pairs = log["mdpairDPDenergy"] - if pairs.size > 1: - return np.mean(pairs[energy_idx:]) - elif pairs.size == 1: - return pairs - -def calculate_pair_energy(A,r,r_cut,num_pol,num_mon,density): - ''' - Calculate the minimum energy for the conservative force to reach at the given radius. - energy for each pair in the system - ''' - density_scaling = (density/1.414)**2 - constant = (1/2)*A*r_cut - U = ((A*(r**2))/(2*r_cut)) - (A*r) + constant - pair_energy = (13*U*num_pol*num_mon*density_scaling)/2 - - return pair_energy - -def simulation_energy_end(A,r,r_cut,num_pol,num_mon,density,energy_idx=-1,log_file_name='log.txt'): - ''' - Calculate the minimum energy for the conservative force to reach at the given radius. - energy for each pair in the system - ''' - U_goal = calculate_pair_energy(A=A,r=r,r_cut=r_cut,num_pol=num_pol,num_mon=num_mon,density=density) - last_U = check_pair_energy(energy_idx=energy_idx,log_file_name=log_file_name) - if last_U <= U_goal: - return True - else: - return False - From 1533ff633c226ebdc38cbb9f127796c00ab0791b Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Mon, 15 Jun 2026 16:23:31 -0600 Subject: [PATCH 13/18] Fix examples 2, 3 to work with updated stop criteria. TODO: need to add some function documentation. --- phantomwalk/examples/2-dpd-energy-gsd.ipynb | 33 +++------ phantomwalk/examples/3-dpd-to-lj-wca.ipynb | 78 ++++----------------- 2 files changed, 24 insertions(+), 87 deletions(-) diff --git a/phantomwalk/examples/2-dpd-energy-gsd.ipynb b/phantomwalk/examples/2-dpd-energy-gsd.ipynb index 515a004..4267829 100644 --- a/phantomwalk/examples/2-dpd-energy-gsd.ipynb +++ b/phantomwalk/examples/2-dpd-energy-gsd.ipynb @@ -21,11 +21,9 @@ "sys.path.append('../lib/')\n", "import create_system_dpd\n", "from create_system_dpd import create_polymer_system_dpd\n", - "from dpd_utils import calculate_pair_energy\n", "import matplotlib\n", "import numpy as np \n", "import gsd, gsd.hoomd \n", - "from cmeutils.sampling import is_equilibrated\n", "import hoomd \n", "import time\n", "import freud\n", @@ -48,37 +46,26 @@ "num_pol=100\n", "num_mon=100\n", "N = num_pol*num_mon\n", - "density=1.3\n", - "A=800\n", - "r_cut=1.13\n", - "min_pair_dist=1.05\n", - "last_dpd_frame, s = create_polymer_system_dpd(\n", + "density=1.1\n", + "A=50000\n", + "r_cut=1.05\n", + "min_pair_dist=0.9\n", + "last_dpd_frame, closest, s = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", " num_mon=num_mon,\n", " density=density,\n", - " k=20000,\n", + " k=50000,\n", " r_cut=r_cut,\n", " A=A,\n", - " gamma=800,\n", + " gamma=1200,\n", " sim_seed=1234,\n", " np_seed=1234,\n", - " energy=True,\n", " min_pair_dist=min_pair_dist,\n", + " energy_scaling = 4,\n", ")\n", "print(f\"Finished in time = {s:.2f}s\")" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "47e9b287-e7bd-4fbb-bfb4-b27c59206eb1", - "metadata": {}, - "outputs": [], - "source": [ - "U = calculate_pair_energy(A=A,r_cut=r_cut,r=min_pair_dist,num_pol=num_pol,num_mon=num_mon,density=density)\n", - "print(U)" - ] - }, { "cell_type": "markdown", "id": "8b99be6a-bab7-453e-94bd-627766707010", @@ -98,14 +85,14 @@ "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", "pairs = log[\"mdpairDPDenergy\"]\n", "print(\"Total steps\",len(pe)*100+100)\n", - "print(\"DPD Energy Cutoff= \",U/N)\n", + "#print(\"DPD Energy Cutoff= \",U/N)\n", "\n", "slice_idx = 0\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", - "plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", + "#plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", "plt.title(\"Energy vs Time\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy\")\n", diff --git a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb index f1ba1c6..4061d95 100644 --- a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb +++ b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb @@ -11,30 +11,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "88429372", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Warning on use of the timeseries module: If the inherent timescales of the system are long compared to those being analyzed, this statistical inefficiency may be an underestimate. The estimate presumes the use of many statistically independent samples. Tests should be performed to assess whether this condition is satisfied. Be cautious in the interpretation of the data.\n", - "\n", - "****** PyMBAR will use 64-bit JAX! *******\n", - "* JAX is currently set to 32-bit bitsize *\n", - "* which is its default. *\n", - "* *\n", - "* PyMBAR requires 64-bit mode and WILL *\n", - "* enable JAX's 64-bit mode when called. *\n", - "* *\n", - "* This MAY cause problems with other *\n", - "* Uses of JAX in the same code. *\n", - "******************************************\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "import sys\n", "import os\n", @@ -46,7 +26,6 @@ "import numpy as np \n", "import gsd, gsd.hoomd \n", "import hoomd \n", - "from cmeutils.sampling import is_equilibrated\n", "import time\n", "import freud\n", "import matplotlib_inline\n", @@ -58,28 +37,20 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": null, "id": "1be674a6-7412-4f5b-aad2-24cb3cc0f869", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Finished in time = 0.30s\n" - ] - } - ], + "outputs": [], "source": [ "num_pol=10\n", "num_mon=100\n", "N = num_pol*num_mon\n", "density=0.8\n", "A=50000\n", - "r_cut=1.2\n", - "min_pair_dist=1.1\n", + "r_cut=1.1\n", + "min_pair_dist=0.95\n", "seed = 122\n", - "dpd_final_frame, s = create_polymer_system_dpd(\n", + "dpd_final_frame, closest, s = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", " num_mon=num_mon,\n", " density=density,\n", @@ -91,6 +62,7 @@ " np_seed=seed,\n", " energy=True,\n", " min_pair_dist=min_pair_dist,\n", + " energy_scaling = .1\n", ")\n", "print(f\"Finished in time = {s:.2f}s\")" ] @@ -105,7 +77,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": null, "id": "6139f233-10e2-4288-8cc0-cc78717d9f35", "metadata": {}, "outputs": [], @@ -116,12 +88,12 @@ " dt=0.0005,\n", " lj_epsilon=1.0,\n", " lj_sigma=1.0,\n", - " lj_r_cut=1.2,\n", + " lj_r_cut=2.5,\n", " fene_k=30,\n", - " fene_r0=1.,\n", + " fene_r0=1.05,\n", " fene_epsilon=1.0,\n", " fene_sigma=1.0,\n", - " fene_delta=0,\n", + " fene_delta=0.05,\n", " angle_k=3.0,\n", " angle_t0=1.0,\n", " dihedral_k=3.0,\n", @@ -263,32 +235,10 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": null, "id": "48f19fa1-6ec6-4244-a23c-fd9d774da3f1", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "**ERROR**: bond.fene: bond out of bounds\n", - "\n" - ] - }, - { - "ename": "RuntimeError", - "evalue": "Error in bond calculation", - "output_type": "error", - "traceback": [ - "\u001b[31m---------------------------------------------------------------------------\u001b[39m", - "\u001b[31mRuntimeError\u001b[39m Traceback (most recent call last)", - "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[25]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m run_lj_simulation(dpd_final_frame=dpd_final_frame)\n", - "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[24]\u001b[39m\u001b[32m, line 140\u001b[39m, in \u001b[36mrun_lj_simulation\u001b[39m\u001b[34m(dpd_final_frame, random_seed, dt, lj_epsilon, lj_sigma, lj_r_cut, fene_k, fene_r0, fene_epsilon, fene_sigma, fene_delta, angle_k, angle_t0, dihedral_k, dihedral_d, dihedral_n, dihedral_phi0)\u001b[39m\n\u001b[32m 136\u001b[39m add_hoomd_writers(sim=LJ_sim)\n\u001b[32m 137\u001b[39m \n\u001b[32m 138\u001b[39m \u001b[38;5;66;03m# Run short equilibration\u001b[39;00m\n\u001b[32m 139\u001b[39m LJ_sim.run(\u001b[32m0\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m140\u001b[39m LJ_sim.run(\u001b[32m100\u001b[39m)\n\u001b[32m 141\u001b[39m \n\u001b[32m 142\u001b[39m \u001b[38;5;66;03m# Flush outputs\u001b[39;00m\n\u001b[32m 143\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m writer \u001b[38;5;28;01min\u001b[39;00m LJ_sim.operations.writers:\n", - "\u001b[36mFile \u001b[39m\u001b[32m~/miniforge3/envs/phantomwalk/lib/python3.12/site-packages/hoomd/simulation.py:561\u001b[39m, in \u001b[36mSimulation.run\u001b[39m\u001b[34m(self, steps, write_at_start)\u001b[39m\n\u001b[32m 558\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m steps_int < \u001b[32m0\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m steps_int > TIMESTEP_MAX - \u001b[32m1\u001b[39m:\n\u001b[32m 559\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[33msteps must be in the range [0, \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mTIMESTEP_MAX\u001b[38;5;250m \u001b[39m-\u001b[38;5;250m \u001b[39m\u001b[32m1\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m]\u001b[39m\u001b[33m\"\u001b[39m)\n\u001b[32m--> \u001b[39m\u001b[32m561\u001b[39m \u001b[30;43mself\u001b[39;49m\u001b[30;43m.\u001b[39;49m\u001b[30;43m_cpp_sys\u001b[39;49m\u001b[30;43m.\u001b[39;49m\u001b[30;43mrun\u001b[39;49m\u001b[30;43m(\u001b[39;49m\u001b[30;43msteps_int\u001b[39;49m\u001b[30;43m,\u001b[39;49m\u001b[30;43m \u001b[39;49m\u001b[30;43mwrite_at_start\u001b[39;49m\u001b[30;43m)\u001b[39;49m\n", - "\u001b[31mRuntimeError\u001b[39m: Error in bond calculation" - ] - } - ], + "outputs": [], "source": [ "run_lj_simulation(dpd_final_frame=dpd_final_frame)" ] From 7b3158620099c81c66564babdcf6c270cda04289 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Tue, 16 Jun 2026 10:10:15 -0600 Subject: [PATCH 14/18] Cleaned up notebook examples, moved lj sim into utils, updated docstrings --- phantomwalk/examples/1-bead-spring-dpd.ipynb | 48 +++-- phantomwalk/examples/2-dpd-energy-gsd.ipynb | 29 +-- phantomwalk/examples/3-dpd-to-lj-wca.ipynb | 207 +++---------------- phantomwalk/lib/create_system_dpd.py | 49 +++-- phantomwalk/lib/dpd_utils.py | 149 +++++++++++++ 5 files changed, 232 insertions(+), 250 deletions(-) diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index 7691510..34a9447 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/phantomwalk/examples/1-bead-spring-dpd.ipynb @@ -77,31 +77,34 @@ "times = []\n", "rs = []\n", "\n", - "density=1.1\n", + "density=1.4 #number density. 1.4 is close-packed, so things get weird above it!\n", "npoly=100\n", "nmono=10\n", - "N=npoly*nmono\n", - "A=50000\n", + "A=50000 #inter-particle repulsion A and bond-stiffness k seem to work well if they're about the same magnitude\n", "k=50000\n", - "dt=0.001\n", - "gamma = 1200\n", - "r_cut=1.01\n", - "min_pair_dist=.95\n", - "es = 0.05\n", - "\n", + "dt=0.001 #step sizes larger than 0.002 seem hard to keep numerically stable\n", + "gamma = 1200 #drag terms larger than 1200 seem to blow up\n", + "r_cut=1.01 #smaller r_cut means smaller neighbor lists, means faster, but larger helps push particles further apart\n", + "min_pair_dist=.80 #this code attempts to run until no two particles are within this distance. \n", + "#a good heuristic is to start with min_pair_dist around 0.7 or 0.8, and creep it up to see how high you can get it.\n", + "loop_timeout = 600 #seconds after which we just stop and give you the last config\n", + "es = 1 #this scaling factor is used to scale the per-particle-energy: make it a fraction to run longer, a large integer to accept higher-energy configurations\n", + "bond_l = 1.0 #bond length\n", + "bond_tolerance = 0.05 #largest (r-r_0) for bonds that we're willing to tolerate, on average\n", + "N=npoly*nmono\n", "for i in range(nruns):\n", " seed = np.random.randint(50000)\n", - " last_frame, closest, s = create_polymer_system_dpd(\n", - " bond_l=1.0,\n", + " last_frame, closest, s, e_cut = create_polymer_system_dpd(\n", + " bond_l=bond_l,\n", " num_pol=npoly,\n", " num_mon=nmono,\n", - " kT=1.0,\n", + " kT=1.0, #For this paper, keep kT const\n", " sim_seed=seed,\n", " np_seed=seed,\n", " sim_steps_incr=100,\n", - " loop_timeout=60,\n", + " loop_timeout=loop_timeout,\n", " gsd_file_name='trajectory.gsd',\n", - " gsd_write_freq=20,\n", + " gsd_write_freq=50,\n", " log_file_name='log.txt',\n", " log_write_freq=10,\n", " density=density,\n", @@ -113,11 +116,12 @@ " k=k,\n", " gamma=gamma,\n", " energy_scaling = es,\n", + " bond_tolerance = bond_tolerance\n", " )\n", " times.append(s)\n", " rs.append(closest)\n", " print(\"{:.2f}s, r = {:.2f} \".format(s,closest))\n", - "print(\"\\n\\n{:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(np.average(times),np.std(times), np.average(rs), np.std(rs)))\n" + "print(\"\\nN={}: {:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(N,np.average(times),np.std(times), np.average(rs), np.std(rs)))\n" ] }, { @@ -129,7 +133,7 @@ "source": [ "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n", "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n", - "plt.title(\"Radial Distribution Function\")\n", + "plt.title(\"Radial Distribution Function (of the last run)\")\n", "plt.xlabel(\"$r$\")\n", "plt.ylabel(\"$g(r)$\")\n", "plt.show()\n", @@ -143,20 +147,20 @@ "metadata": {}, "outputs": [], "source": [ + "slice_idx = 10\n", + "short_idx= None\n", + "\n", "log = np.genfromtxt(\"log.txt\", names=True)\n", "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", "pairs = log[\"mdpairDPDenergy\"]\n", "print(\"Total steps\",len(pe)*100+100)\n", - "\n", - "slice_idx = 40\n", - "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", - "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", + "plt.plot(x_values,(pe[slice_idx:short_idx] - pairs[slice_idx:short_idx])/N, label=\"bond energy\")\n", "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", - "plt.title(\"Energy vs Time\")\n", + "plt.hlines(e_cut,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", + "plt.title(\"Energy vs Time (of the last run)\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy/N\")\n", - "\n", "plt.legend()" ] }, diff --git a/phantomwalk/examples/2-dpd-energy-gsd.ipynb b/phantomwalk/examples/2-dpd-energy-gsd.ipynb index 4267829..d8af00c 100644 --- a/phantomwalk/examples/2-dpd-energy-gsd.ipynb +++ b/phantomwalk/examples/2-dpd-energy-gsd.ipynb @@ -6,7 +6,7 @@ "metadata": {}, "source": [ "## Initializing polymer systems using a DPD potential\n", - "This notebook walks through the PhantomWalk functions for packing linear polymers in a box. The polymers are first placed in a cubic box using a random walk. Then a short HOOMD simulation is run with the soft force potential of Dissipative Particle Dynamics. The simulation ends when the pair energy from the DPD potential reaches a stable state, as checked with an autocorrelation function." + "This notebook walks through the PhantomWalk functions for packing linear polymers in a box. The polymers are first placed in a cubic box using a random walk. Then a short HOOMD simulation is run with the soft force potential of Dissipative Particle Dynamics. The simulation ends when the pair energy from the DPD potential reaches a stable state." ] }, { @@ -50,7 +50,7 @@ "A=50000\n", "r_cut=1.05\n", "min_pair_dist=0.9\n", - "last_dpd_frame, closest, s = create_polymer_system_dpd(\n", + "last_dpd_frame, closest, s, e_cut = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", " num_mon=num_mon,\n", " density=density,\n", @@ -84,15 +84,16 @@ "log = np.genfromtxt(\"log.txt\", names=True)\n", "pe = log[\"mdcomputeThermodynamicQuantitiespotential_energy\"]\n", "pairs = log[\"mdpairDPDenergy\"]\n", + "bonds = log[\"mdbondHarmonicenergy\"]\n", "print(\"Total steps\",len(pe)*100+100)\n", - "#print(\"DPD Energy Cutoff= \",U/N)\n", + "print(\"DPD Energy Cutoff= \",e_cut)\n", "\n", "slice_idx = 0\n", "short_idx= None\n", "x_values = range(slice_idx, len(pe[:short_idx]))\n", "plt.plot(x_values,pe[slice_idx:short_idx]/N, label=\"potential energy\")\n", "plt.plot(x_values,pairs[slice_idx:short_idx]/N, label=\"DPD pair energy\")\n", - "#plt.hlines(U/N,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", + "plt.hlines(e_cut,xmin=slice_idx,xmax=len(pairs[:short_idx]),color=\"blue\",linestyle=\"--\",label=\"Calculated Pair Energy Cutoff\")\n", "plt.title(\"Energy vs Time\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy\")\n", @@ -101,20 +102,6 @@ "#plt.savefig(\"10-10mers-energy-cut\")" ] }, - { - "cell_type": "raw", - "id": "674da1d0-3767-4475-8091-3a66564779b4", - "metadata": {}, - "source": [ - "plt.plot(bonds, label=\"bond energy\")\n", - "plt.title(\"Energy vs Time\")\n", - "plt.xlabel(\"Frame (time/log_freq)\")\n", - "plt.ylabel(\"Energy\")\n", - "\n", - "plt.legend()\n", - "#plt.savefig(\"100-100mers-new-energy-cut-0.9p.png\")" - ] - }, { "cell_type": "code", "execution_count": null, @@ -122,12 +109,11 @@ "metadata": {}, "outputs": [], "source": [ - "bonds = log[\"mdbondHarmonicenergy\"]\n", + "plt.legend()\n", "plt.plot(bonds, label=\"bond energy\")\n", "plt.title(\"Bond Energy vs Time\")\n", "plt.xlabel(\"Frame (time/log_freq)\")\n", "plt.ylabel(\"Energy\")\n", - "\n", "plt.legend()" ] }, @@ -143,8 +129,7 @@ "plt.title(\"Radial Distribution Function\")\n", "plt.xlabel(\"$r$\")\n", "plt.ylabel(\"$g(r)$\")\n", - "plt.show()\n", - "print(rdf_data)" + "plt.show()" ] }, { diff --git a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb index 4061d95..5400cfa 100644 --- a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb +++ b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb @@ -5,8 +5,10 @@ "id": "9ba5e733", "metadata": {}, "source": [ - "## Initializing polymer systems using a DPD potential\n", - "This notebook walks through the PhantomWalk functions for packing linear polymers in a box. The polymers are first placed in a cubic box using a random walk. Then a short HOOMD simulation is run with the soft force potential of Dissipative Particle Dynamics. The simulation ends when the pair energy from the DPD potential reaches a stable state, as checked with an autocorrelation function." + "## Testing \"hard\" systems after DPD initialization\n", + "This notebook takes the output of a DPD simulation and uses it as input in a lennard-jones (LJ) simulation that has much stronger repulsions between particles, and bonds that cannot extend beyond the fene_r0 length. \n", + "\n", + "If the LJ simulation runs, it is evidence that we have a very stable configuration from DPD. It's kindof anticlimactic." ] }, { @@ -21,7 +23,7 @@ "sys.path.append('../lib/')\n", "import create_system_dpd\n", "from create_system_dpd import create_polymer_system_dpd\n", - "from dpd_utils import add_hoomd_writers\n", + "from dpd_utils import add_hoomd_writers, run_lj_simulation\n", "import matplotlib\n", "import numpy as np \n", "import gsd, gsd.hoomd \n", @@ -42,214 +44,51 @@ "metadata": {}, "outputs": [], "source": [ - "num_pol=10\n", + "num_pol=100\n", "num_mon=100\n", "N = num_pol*num_mon\n", - "density=0.8\n", + "density=1.4\n", "A=50000\n", - "r_cut=1.1\n", - "min_pair_dist=0.95\n", - "seed = 122\n", - "dpd_final_frame, closest, s = create_polymer_system_dpd(\n", + "k=50000\n", + "r_cut=1.01\n", + "min_pair_dist=0.80\n", + "seed = 123\n", + "dpd_final_frame, closest, s, e_cut = create_polymer_system_dpd(\n", " num_pol=num_pol,\n", " num_mon=num_mon,\n", " density=density,\n", - " k=50000,\n", - " r_cut=r_cut,\n", " A=A,\n", + " k=k,\n", + " r_cut=r_cut,\n", + " bond_l = 0.95,\n", " gamma=1200,\n", " sim_seed=seed,\n", " np_seed=seed,\n", - " energy=True,\n", + " write = True,\n", " min_pair_dist=min_pair_dist,\n", - " energy_scaling = .1\n", + " energy_scaling = .5,\n", + " bond_tolerance = 0.04,\n", ")\n", - "print(f\"Finished in time = {s:.2f}s\")" - ] - }, - { - "cell_type": "markdown", - "id": "041801b6-5c19-42b6-a597-66dd526adeae", - "metadata": {}, - "source": [ - "## Run a Lennard Jones Simulation from the Final DPD Frame" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6139f233-10e2-4288-8cc0-cc78717d9f35", - "metadata": {}, - "outputs": [], - "source": [ - "def run_lj_simulation(\n", - " dpd_final_frame,\n", + "print(f\"Finished in time = {s:.2f}s, r_min = {closest:.2f}\")\n", + "lj = run_lj_simulation(dpd_final_frame=dpd_final_frame,\n", + " fene_r0=1.15,\n", + " fene_delta=0.0,\n", " random_seed=25,\n", " dt=0.0005,\n", " lj_epsilon=1.0,\n", " lj_sigma=1.0,\n", " lj_r_cut=2.5,\n", " fene_k=30,\n", - " fene_r0=1.05,\n", " fene_epsilon=1.0,\n", " fene_sigma=1.0,\n", - " fene_delta=0.05,\n", " angle_k=3.0,\n", " angle_t0=1.0,\n", " dihedral_k=3.0,\n", " dihedral_d=-1,\n", " dihedral_n=3,\n", " dihedral_phi0=0\n", - "):\n", - " \"\"\"Run an LJ + FENE + angle + dihedral equilibration simulation in HOOMD-blue.\n", - "\n", - " Parameters\n", - " ----------\n", - " dpd_final_frame : gsd.hoomd.Frame\n", - " Initial configuration used to start the LJ simulation.\n", - "\n", - " random_seed : int, optional, default 24\n", - " Random seed for reproducibility.\n", - "\n", - " dt : float, optional, default 0.001\n", - " Integration timestep.\n", - "\n", - " lj_epsilon : float, optional, default 1.0\n", - " Lennard-Jones interaction strength.\n", - "\n", - " lj_sigma : float, optional, default 1.0\n", - " Lennard-Jones particle size parameter.\n", - "\n", - " lj_r_cut : float, optional, default 1.2\n", - " Lennard-Jones cutoff radius.\n", - "\n", - " fene_k : float, optional, default 30\n", - " FENE bond spring constant.\n", - "\n", - " fene_r0 : float, optional, default 1.05\n", - " FENE maximum bond extension parameter.\n", - "\n", - " fene_epsilon : float, optional, default 1.0\n", - " FENE-WCA epsilon parameter.\n", - "\n", - " fene_sigma : float, optional, default 1.0\n", - " FENE-WCA sigma parameter.\n", - "\n", - " fene_delta : float, optional, default 0\n", - " FENE potential shift parameter.\n", - "\n", - " angle_k : float, optional, default 3.0\n", - " Harmonic angle force constant.\n", - "\n", - " angle_t0 : float, optional, default 1.0\n", - " Equilibrium bond angle (radians).\n", - "\n", - " dihedral_k : float, optional, default 3.0\n", - " Dihedral force constant.\n", - "\n", - " dihedral_d : int, optional, default -1\n", - " Dihedral sign parameter.\n", - "\n", - " dihedral_n : int, optional, default 3\n", - " Dihedral periodicity.\n", - "\n", - " dihedral_phi0 : float, optional, default 0\n", - " Dihedral phase offset (radians).\n", - "\n", - " Returns\n", - " -------\n", - " hoomd.Simulation\n", - " HOOMD simulation object after short equilibration run.\n", - " \"\"\"\n", - "\n", - " forces = []\n", - "\n", - " # Pair force (LJ)\n", - " nlist = hoomd.md.nlist.Cell(buffer=0.40, exclusions=[\"bond\"])\n", - " lj = hoomd.md.pair.LJ(nlist=nlist)\n", - " lj.params[('A', 'A')] = dict(epsilon=lj_epsilon, sigma=lj_sigma)\n", - " lj.r_cut[('A', 'A')] = lj_r_cut\n", - " forces.append(lj)\n", - "\n", - " # FENE bonds\n", - " fene_bond = hoomd.md.bond.FENEWCA()\n", - " fene_bond.params['b'] = dict(\n", - " k=fene_k,\n", - " r0=fene_r0,\n", - " epsilon=fene_epsilon,\n", - " sigma=fene_sigma,\n", - " delta=fene_delta,\n", - " )\n", - " forces.append(fene_bond)\n", - " \n", - " ''' TODO add angles and dihedrals back into frame generation\n", - " # Angle potential\n", - " harmonic_angle = hoomd.md.angle.Harmonic()\n", - " harmonic_angle.params[\"A-A-A\"] = dict(k=angle_k, t0=angle_t0)\n", - " forces.append(harmonic_angle)\n", - "\n", - " # Dihedral potential\n", - " dihedral = hoomd.md.dihedral.Periodic()\n", - " dihedral.params[\"A-A-A-A\"] = dict(\n", - " k=dihedral_k,\n", - " d=dihedral_d,\n", - " n=dihedral_n,\n", - " phi0=dihedral_phi0\n", - " )\n", - " forces.append(dihedral)\n", - " '''\n", - " \n", - " # Integrator\n", - " integrator_lj = hoomd.md.Integrator(dt=dt)\n", - " integrator_lj.forces = forces\n", - "\n", - " integrator_lj.methods.append(\n", - " hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())\n", - " )\n", - "\n", - " # Simulation\n", - " LJ_sim = hoomd.Simulation(\n", - " device=hoomd.device.auto_select(),\n", - " seed=random_seed\n", - " )\n", - "\n", - " LJ_sim.create_state_from_snapshot(snapshot=dpd_final_frame)\n", - " LJ_sim.operations.integrator = integrator_lj\n", - "\n", - " # Use your shared writer setup\n", - " add_hoomd_writers(sim=LJ_sim)\n", - "\n", - " # Run short equilibration\n", - " LJ_sim.run(0)\n", - " LJ_sim.run(100)\n", - "\n", - " # Flush outputs\n", - " for writer in LJ_sim.operations.writers:\n", - " if hasattr(writer, \"flush\"):\n", - " writer.flush()\n", - "\n", - " print(\"LJ simulation finished.\")\n", - "\n", - " return LJ_sim" + " )" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "48f19fa1-6ec6-4244-a23c-fd9d774da3f1", - "metadata": {}, - "outputs": [], - "source": [ - "run_lj_simulation(dpd_final_frame=dpd_final_frame)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0cbf4dee-ccad-42b5-bde0-538965f8ae8f", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index b5d7499..6d77d7d 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -6,6 +6,11 @@ def get_close(rdf): + ''' + Find closest separation between two particles from first nonzero bin of the rdf + + returns value of the bin center with the nonzero rdf + ''' b =(rdf.rdf !=0).argmax() return rdf.bin_centers[b] @@ -13,20 +18,19 @@ def create_polymer_system_dpd( num_pol, num_mon, density, - k=20000, + A=50000, + k=50000, bond_l=1.0, - r_cut=1.15, + r_cut=1.01, kT=1.0, - A=800, - gamma=800, + gamma=1200, dt=0.001, sim_seed=1234, np_seed=1234, sim_steps_incr=100, loop_timeout=60, - energy=True, - min_pair_dist=1.05, - energy_scaling= 5, + min_pair_dist=0.80, + energy_scaling= 1, bond_tolerance = 0.05, write=True, gsd_file_name='trajectory.gsd', @@ -46,21 +50,21 @@ def create_polymer_system_dpd( length of polymers in system density : float, required number density to initalize the system - k : int, default 20000 + A : float, default 50000 + DPD force parameter + k : int, default 50000 spring constant for harmonic bonds bond_l : float, default 1.0 harmonic bond rest length - r_cut : float, default 1.15 + r_cut : float, default 1.01 cutoff pair distance for neighbor list kT : float, default 1.0 temperature of thermostat - A : float, default 1000 - DPD force parameter - gamma : float, default 800 + gamma : float, default 1200 DPD drag parameter (mass/time) dt : float, default 0.001 timestep for HOOMD simulation - sim_seed : int, default 123 + sim_seed : int, default 1234 seed for the HOOMD simulation state np_seed : int, default 1234 seed for random number generator in random walk @@ -68,10 +72,11 @@ def create_polymer_system_dpd( the number of steps to run in a loop before checking simulation end criteria loop_timeout : int, default 60 seconds time out to manually end the simulation before it reaches the cutoff, meant to prevent large file creation - energy : bool, default True - trigger to use energy cutoff instead of manually building neighbor list - min_pair_dist : float, default 1.05 - condition for ending the soft push simulation + min_pair_dist : float, default 0.8 + run until no two particles are within this distance + energy_scaling : float, default 1 + scaling factor to manually adjust per-particle energy cutoff stop criteria + Fractions (0.5, 0.2) will lower the threshold (longer sims), and large numbers (10,15) will shorten simulations. write : bool, True trigger for writing out gsd and log files gsd_file_name : str, default 'trajectory.gsd' @@ -122,7 +127,7 @@ def create_polymer_system_dpd( maxPerParticle = A*( (min_pair_dist*min_pair_dist)/(2*r_cut) - min_pair_dist + r_cut/2) maxPerParticle *= density*density*energy_scaling maxPerBond = k*bond_tolerance*bond_tolerance/2 - print("max per particle= {:.2f}, max per bond= {:.2f}".format(maxPerParticle, maxPerBond)) + #print("max per particle= {:.2f}, max per bond= {:.2f}".format(maxPerParticle, maxPerBond)) if write: rdf,thermo = add_hoomd_writers( simulation, gsd_file_name, gsd_write_freq, log_file_name,log_write_freq ) @@ -137,7 +142,7 @@ def create_polymer_system_dpd( check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in energy") - return simulation.state.get_snapshot(), loop_timeout + return simulation.state.get_snapshot(), get_close(rdf), loop_timeout simulation.run(sim_steps_incr) for writer in simulation.operations.writers: if hasattr(writer, "flush"): @@ -147,7 +152,7 @@ def create_polymer_system_dpd( check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in bond energy") - return simulation.state.get_snapshot(), loop_timeout + return simulation.state.get_snapshot(), get_close(rdf), loop_timeout simulation.run(sim_steps_incr) for writer in simulation.operations.writers: if hasattr(writer, "flush"): @@ -158,7 +163,7 @@ def create_polymer_system_dpd( check_time = time.perf_counter() if (check_time-start_time) > loop_timeout: print("Simulation timed out in rdf polish") - return simulation.state.get_snapshot(), loop_timeout + return simulation.state.get_snapshot(), get_close(rdf), loop_timeout simulation.run(sim_steps_incr) closest = get_close(rdf) for writer in simulation.operations.writers: @@ -168,4 +173,4 @@ def create_polymer_system_dpd( end_time = time.perf_counter() total_time = end_time - start_time np.savetxt( "rdf.csv", np.vstack((rdf.bin_centers, rdf.rdf)).T, delimiter=",", header="r, g(r)") - return simulation.state.get_snapshot(), closest, total_time + return simulation.state.get_snapshot(), closest, total_time, maxPerParticle diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 3be7077..ce27992 100644 --- a/phantomwalk/lib/dpd_utils.py +++ b/phantomwalk/lib/dpd_utils.py @@ -145,3 +145,152 @@ def act(self, timestep): sim.operations.writers.append(table_file) return rdf, thermo_props +def run_lj_simulation( + dpd_final_frame, + random_seed=25, + dt=0.0005, + lj_epsilon=1.0, + lj_sigma=1.0, + lj_r_cut=2.5, + fene_k=30, + fene_r0=1.01, + fene_epsilon=1.0, + fene_sigma=1.0, + fene_delta=0.05, + angle_k=3.0, + angle_t0=1.0, + dihedral_k=3.0, + dihedral_d=-1, + dihedral_n=3, + dihedral_phi0=0 +): + """Run an LJ + FENE + angle + dihedral equilibration simulation in HOOMD-blue. + + Parameters + ---------- + dpd_final_frame : gsd.hoomd.Frame + Initial configuration used to start the LJ simulation. + + random_seed : int, optional, default 24 + Random seed for reproducibility. + + dt : float, optional, default 0.001 + Integration timestep. + + lj_epsilon : float, optional, default 1.0 + Lennard-Jones interaction strength. + + lj_sigma : float, optional, default 1.0 + Lennard-Jones particle size parameter. + + lj_r_cut : float, optional, default 1.2 + Lennard-Jones cutoff radius. + + fene_k : float, optional, default 30 + FENE bond spring constant. + + fene_r0 : float, optional, default 1.05 + FENE maximum bond extension parameter. + + fene_epsilon : float, optional, default 1.0 + FENE-WCA epsilon parameter. + + fene_sigma : float, optional, default 1.0 + FENE-WCA sigma parameter. + + fene_delta : float, optional, default 0 + FENE potential shift parameter. + + angle_k : float, optional, default 3.0 + Harmonic angle force constant. + + angle_t0 : float, optional, default 1.0 + Equilibrium bond angle (radians). + + dihedral_k : float, optional, default 3.0 + Dihedral force constant. + + dihedral_d : int, optional, default -1 + Dihedral sign parameter. + + dihedral_n : int, optional, default 3 + Dihedral periodicity. + + dihedral_phi0 : float, optional, default 0 + Dihedral phase offset (radians). + + Returns + ------- + hoomd.Simulation + HOOMD simulation object after short equilibration run. + """ + + forces = [] + + # Pair force (LJ) + nlist = hoomd.md.nlist.Cell(buffer=0.40, exclusions=["bond"]) + lj = hoomd.md.pair.LJ(nlist=nlist) + lj.params[('A', 'A')] = dict(epsilon=lj_epsilon, sigma=lj_sigma) + lj.r_cut[('A', 'A')] = lj_r_cut + forces.append(lj) + + # FENE bonds + fene_bond = hoomd.md.bond.FENEWCA() + fene_bond.params['b'] = dict( + k=fene_k, + r0=fene_r0, + epsilon=fene_epsilon, + sigma=fene_sigma, + delta=fene_delta, + ) + forces.append(fene_bond) + + ''' TODO add angles and dihedrals back into frame generation + # Angle potential + harmonic_angle = hoomd.md.angle.Harmonic() + harmonic_angle.params["A-A-A"] = dict(k=angle_k, t0=angle_t0) + forces.append(harmonic_angle) + + # Dihedral potential + dihedral = hoomd.md.dihedral.Periodic() + dihedral.params["A-A-A-A"] = dict( + k=dihedral_k, + d=dihedral_d, + n=dihedral_n, + phi0=dihedral_phi0 + ) + forces.append(dihedral) + ''' + + # Integrator + integrator_lj = hoomd.md.Integrator(dt=dt) + integrator_lj.forces = forces + + integrator_lj.methods.append( + hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All()) + ) + + # Simulation + LJ_sim = hoomd.Simulation( + device=hoomd.device.auto_select(), + seed=random_seed + ) + + LJ_sim.create_state_from_snapshot(snapshot=dpd_final_frame) + LJ_sim.operations.integrator = integrator_lj + + # Use your shared writer setup + add_hoomd_writers(sim=LJ_sim) + + # Run short equilibration + LJ_sim.run(0) + LJ_sim.run(100) + + # Flush outputs + for writer in LJ_sim.operations.writers: + if hasattr(writer, "flush"): + writer.flush() + + print("LJ simulation finished.") + + return LJ_sim From b320ad777560c3d2f23b41d8122fc070f7d5f42a Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Tue, 16 Jun 2026 10:15:31 -0600 Subject: [PATCH 15/18] examples weren't finding libs, but does this break cedar's tests? --- phantomwalk/lib/create_system_dpd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index cd3d1fd..575afe8 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -2,8 +2,8 @@ import hoomd import time -#from dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers -from phantomwalk.lib.dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers +from dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers +#from phantomwalk.lib.dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers def get_close(rdf): From aa2dde40315a5392e7e8cc67768b7745c6a590d9 Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Tue, 16 Jun 2026 10:22:26 -0600 Subject: [PATCH 16/18] installed via pip --- phantomwalk/lib/create_system_dpd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/phantomwalk/lib/create_system_dpd.py b/phantomwalk/lib/create_system_dpd.py index 575afe8..cd3d1fd 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -2,8 +2,8 @@ import hoomd import time -from dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers -#from phantomwalk.lib.dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers +#from dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers +from phantomwalk.lib.dpd_utils import initialize_snapshot_rand_walk,add_hoomd_writers def get_close(rdf): From e13f4b33a33b04b0a81b9a1b73d0e5c078eb0b6f Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Tue, 16 Jun 2026 10:22:56 -0600 Subject: [PATCH 17/18] Fix path problem --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 1a4f566..de85771 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ *.txt *.swp notes +phantomwalk.egg-info/ # Byte-compiled / optimized / DLL files __pycache__/ From 29aea2edf1e8bf46f48d228153abcd3c7eaec49a Mon Sep 17 00:00:00 2001 From: Eric Jankowski Date: Tue, 16 Jun 2026 10:33:41 -0600 Subject: [PATCH 18/18] pass tests --- phantomwalk/tests/test_create_system_dpd.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/phantomwalk/tests/test_create_system_dpd.py b/phantomwalk/tests/test_create_system_dpd.py index 3ab658a..e4635cf 100644 --- a/phantomwalk/tests/test_create_system_dpd.py +++ b/phantomwalk/tests/test_create_system_dpd.py @@ -14,17 +14,20 @@ def rm_files(*files): pass def test_creation(): - snap, time = dpd.create_polymer_system_dpd(num_pol=5, num_mon=10, density=0.5) - assert time > 0 + loop_timeout = 60 + snap, closest, time, e_cut = dpd.create_polymer_system_dpd(num_pol=5, num_mon=10, density=0.5, loop_timeout = loop_timeout) + assert time < loop_timeout def test_custom_log_files(): # remove files that might've been output by other tests so that file loading # can fail if it's wrong - rm_files('log.txt', 'trajectory.gsd') + rm_files('log.txt', 'trajectory.gsd', 'rdf.csv') time_string = datetime.datetime.now().strftime("%Y-%m-%d-%H-%M-%S") gsd_file_name = f"{time_string}.gsd" log_file_name = f"{time_string}.txt" + rdf_file_name = f"rdf.csv" + s = dpd.create_polymer_system_dpd( num_pol=5, num_mon=10, @@ -35,8 +38,9 @@ def test_custom_log_files(): gsd_exists = os.path.isfile(gsd_file_name) log_exists = os.path.isfile(log_file_name) + rdf_exists = os.path.isfile(rdf_file_name) # clean up after ourselves - rm_files(gsd_file_name, log_file_name) + rm_files(gsd_file_name, log_file_name, rdf_file_name) assert gsd_exists and log_exists