diff --git a/.gitignore b/.gitignore index 0971d20..de85771 100644 --- a/.gitignore +++ b/.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__/ diff --git a/phantomwalk/examples/1-bead-spring-dpd.ipynb b/phantomwalk/examples/1-bead-spring-dpd.ipynb index c2fe0e8..34a9447 100644 --- a/phantomwalk/examples/1-bead-spring-dpd.ipynb +++ b/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, diff --git a/phantomwalk/examples/2-dpd-energy-gsd.ipynb b/phantomwalk/examples/2-dpd-energy-gsd.ipynb index 449a253..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." ] }, { @@ -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=0.9\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, e_cut = 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", @@ -97,13 +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= \",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], 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(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", @@ -113,39 +103,39 @@ ] }, { - "cell_type": "raw", - "id": "674da1d0-3767-4475-8091-3a66564779b4", + "cell_type": "code", + "execution_count": null, + "id": "0b947636", "metadata": {}, + "outputs": [], "source": [ + "plt.legend()\n", "plt.plot(bonds, label=\"bond energy\")\n", - "plt.title(\"Energy vs Time\")\n", + "plt.title(\"Bond 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\")" + "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, - "id": "0b947636", + "id": "133f7da1-415e-49d6-9731-b7f642a62469", "metadata": {}, "outputs": [], "source": [ - "bonds = log[\"mdbondHarmonicenergy\"]\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()" + "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()" ] }, { "cell_type": "code", "execution_count": null, - "id": "133f7da1-415e-49d6-9731-b7f642a62469", + "id": "315bac34-3982-4a53-8466-3ff095590187", "metadata": {}, "outputs": [], "source": [] @@ -167,7 +157,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.13" + "version": "3.12.13" } }, "nbformat": 4, diff --git a/phantomwalk/examples/3-dpd-to-lj-wca.ipynb b/phantomwalk/examples/3-dpd-to-lj-wca.ipynb index f54dfaf..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,12 +23,11 @@ "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", "import hoomd \n", - "from cmeutils.sampling import is_equilibrated\n", "import time\n", "import freud\n", "import matplotlib_inline\n", @@ -46,209 +47,48 @@ "num_pol=100\n", "num_mon=100\n", "N = num_pol*num_mon\n", - "density=0.9\n", - "A=1000\n", - "r_cut=1.13\n", - "min_pair_dist=1.05\n", - "dpd_final_frame, s = create_polymer_system_dpd(\n", + "density=1.4\n", + "A=50000\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=20000,\n", - " r_cut=r_cut,\n", " A=A,\n", - " gamma=800,\n", - " sim_seed=1234,\n", - " np_seed=1234,\n", - " energy=True,\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", + " write = True,\n", " min_pair_dist=min_pair_dist,\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", - " random_seed=24,\n", - " dt=0.001,\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=1.2,\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,\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": { @@ -267,7 +107,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 c07092d..cd3d1fd 100644 --- a/phantomwalk/lib/create_system_dpd.py +++ b/phantomwalk/lib/create_system_dpd.py @@ -1,36 +1,44 @@ import numpy as np -import freud -import gsd, gsd.hoomd import hoomd import time -from phantomwalk.lib.dpd_utils import initialize_snapshot_rand_walk,check_bond_length_equilibration,check_inter_particle_distance,add_hoomd_writers,simulation_energy_end +#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): + ''' + 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] + 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, + min_pair_dist=0.80, + energy_scaling= 1, + 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. @@ -43,21 +51,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 @@ -65,10 +73,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' @@ -90,11 +99,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, @@ -104,7 +109,6 @@ def create_polymer_system_dpd( ) 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) @@ -119,52 +123,54 @@ def create_polymer_system_dpd( DPD = hoomd.md.pair.DPD(nlist, default_r_cut=r_cut, kT=kT) DPD.params[('A', 'A')] = dict(A=A, gamma=gamma) 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 per particle= {:.2f}, max per bond= {:.2f}".format(maxPerParticle, maxPerBond)) if write: - 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 DPD.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(), get_close(rdf), loop_timeout + simulation.run(sim_steps_incr) + for writer in simulation.operations.writers: + if hasattr(writer, "flush"): + writer.flush() + + 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") + 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"): + writer.flush() + + 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(), get_close(rdf), 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) - 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, maxPerParticle diff --git a/phantomwalk/lib/dpd_utils.py b/phantomwalk/lib/dpd_utils.py index 86983d3..ce27992 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): ''' @@ -11,36 +9,27 @@ def initialize_snapshot_rand_walk(num_pol, num_mon, density, bond_length=1.0, se ''' 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 @@ -48,51 +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] - + 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", @@ -128,6 +75,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=100, r_max=2.0) + rdf_calc = FreudRDFCalc(sim, rdf) + + gsd_logger = hoomd.logging.Logger( categories=["scalar", "string", "sequence"] ) @@ -174,54 +137,160 @@ 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, 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. +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 ---------- - energy_idx : int, default -1 - Number of initial simulation steps to discard before - performing equilibration analysis. Default is to return the last frame. + 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 ------- - float, energy of last frame(s) of dpd simulation - + hoomd.Simulation + HOOMD simulation object after short equilibration run. """ - 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 + + 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) -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 + ''' 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) ''' - 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 + + # Integrator + integrator_lj = hoomd.md.Integrator(dt=dt) + integrator_lj.forces = forces - return pair_energy + integrator_lj.methods.append( + hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All()) + ) -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 + # 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 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