Auto rdf

.gitignore

@@ -1,64 +1,18 @@
 # Output files from MD runs
 *.gsd
 *.txt
+*.swp
+notes
+phantomwalk.egg-info/
 
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[codz]
 *$py.class
 
-# Distribution / packaging
-dist/
-*.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/
+phantomwalk/examples/rdf.csv
+phantomwalk/examples/trajectory.gsd
+phantomwalk/examples/log.txt
 
 # Jupyter Notebook
 .ipynb_checkpoints
@@ -67,126 +21,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__/

phantomwalk/examples/1-bead-spring-dpd.ipynb

@@ -73,36 +73,101 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "last_frame, s = create_polymer_system_dpd(\n",
-    "    num_pol=100,\n",
-    "    num_mon=100,\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",
-    "    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=10,\n",
-    "    log_file_name='log.txt',\n",
-    "    log_write_freq=10)\n",
+    "nruns = 5\n",
+    "times = []\n",
+    "rs = []\n",
     "\n",
-    "print(f\"Finished in time = {s:.2f}s\")"
+    "density=1.4 #number density. 1.4 is close-packed, so things get weird above it!\n",
+    "npoly=100\n",
+    "nmono=10\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 #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, e_cut = create_polymer_system_dpd(\n",
+    "        bond_l=bond_l,\n",
+    "        num_pol=npoly,\n",
+    "        num_mon=nmono,\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=loop_timeout,\n",
+    "        gsd_file_name='trajectory.gsd',\n",
+    "        gsd_write_freq=50,\n",
+    "        log_file_name='log.txt',\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",
+    "        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(\"\\nN={}: {:.2f} ({:.2f})s , r = {:.2f} ({:.2f})\".format(N,np.average(times),np.std(times), np.average(rs), np.std(rs)))\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d4262eec-2f8d-44d9-b7f6-bdfd0c55329b",
+   "id": "31c44f45-e350-431e-bb38-a8ec845aceae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rdf_data = np.genfromtxt(\"rdf.csv\", delimiter=\",\")\n",
+    "plt.plot(rdf_data[:, 0], rdf_data[:, 1])\n",
+    "plt.title(\"Radial Distribution Function (of the last run)\")\n",
+    "plt.xlabel(\"$r$\")\n",
+    "plt.ylabel(\"$g(r)$\")\n",
+    "plt.show()\n",
+    "#print(rdf_data)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5f279783-91f5-4560-8d38-c8be52c09443",
+   "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",
+    "x_values = range(slice_idx, len(pe[:short_idx]))\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.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",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "59980585-3703-42c8-8c71-81e76eefcbb3",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -124,7 +189,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.13"
+   "version": "3.12.13"
   }
  },
  "nbformat": 4,