From cfa136a22cb73fb764a6b99b3b5d9653b45c33de Mon Sep 17 00:00:00 2001 From: Rebecca Stephan Date: Mon, 12 Aug 2024 13:13:10 +0200 Subject: [PATCH 1/2] Add Zndraw in tutorials --- doc/tutorials/constant_pH/constant_pH.ipynb | 44 ++++- .../electrodes/electrodes_part2.ipynb | 66 +++++-- .../ferrofluid/ferrofluid_part1.ipynb | 174 ++---------------- .../ferrofluid/ferrofluid_part2.ipynb | 85 +-------- .../ferrofluid/ferrofluid_part3.ipynb | 45 ++++- .../langevin_dynamics/langevin_dynamics.ipynb | 74 +++++++- .../raspberry_electrophoresis.ipynb | 52 ++++++ 7 files changed, 277 insertions(+), 263 deletions(-) diff --git a/doc/tutorials/constant_pH/constant_pH.ipynb b/doc/tutorials/constant_pH/constant_pH.ipynb index e147a6c3bb..8dabb3dc2b 100644 --- a/doc/tutorials/constant_pH/constant_pH.ipynb +++ b/doc/tutorials/constant_pH/constant_pH.ipynb @@ -153,6 +153,7 @@ "import espressomd.electrostatics\n", "import espressomd.reaction_methods\n", "import espressomd.polymer\n", + "import espressomd.zn\n", "from espressomd.interactions import HarmonicBond" ] }, @@ -569,7 +570,7 @@ "\n", "# add thermostat and short integration to let the system relax further\n", "system.thermostat.set_langevin(kT=KT_REDUCED, gamma=1.0, seed=7)\n", - "system.integrator.run(steps=1000)\n", + "system.integrator.run(1000)\n", "\n", "if USE_ELECTROSTATICS:\n", " COULOMB_PREFACTOR=BJERRUM_LENGTH_REDUCED * KT_REDUCED\n", @@ -739,6 +740,40 @@ "outputs": [], "source": [] }, + { + "cell_type": "markdown", + "id": "cd417295-bab3-46a5-a574-fc60d06371a5", + "metadata": {}, + "source": [ + "We will initialize the zndraw visualizer to visualize the simulation for increasing pH value:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04dd303c", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "color = {TYPES[\"HA\"]: \"#7fc454\", #green\n", + " TYPES[\"A\"]: \"#225204\", #dark green\n", + " TYPES[\"B\"]: \"#fca000\", #orange\n", + " TYPES[\"Na\"]: \"#ff0000\", #red\n", + " TYPES[\"Cl\"]: \"#030ffc\" #blue\n", + " }\n", + "radii = {TYPES[\"HA\"]: 2,\n", + " TYPES[\"A\"]: 2,\n", + " TYPES[\"B\"]: 2,\n", + " TYPES[\"Na\"]: 2,\n", + " TYPES[\"Cl\"]: 2\n", + " }\n", + "\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", + "vis.update()" + ] + }, { "cell_type": "markdown", "id": "07daa583", @@ -778,11 +813,14 @@ "outputs": [], "source": [ "# SOLUTION CELL\n", - "def perform_sampling(type_A, num_samples, num_As:np.ndarray, reaction_steps, \n", + "def perform_sampling(type_A, num_samples, num_As:np.ndarray, reaction_steps,\n", " prob_integration=0.5, integration_steps=1000):\n", " for i in range(num_samples):\n", " if USE_WCA and np.random.random() < prob_integration:\n", - " system.integrator.run(integration_steps)\n", + " for _ in range(integration_steps):\n", + " system.integrator.run(1)\n", + " global vis\n", + " vis.update()\n", " # we should do at least one reaction attempt per reactive particle\n", " RE.reaction(steps=reaction_steps)\n", " num_As[i] = system.number_of_particles(type=type_A)" diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb index 4cae3ebcd5..1aaf34d86e 100644 --- a/doc/tutorials/electrodes/electrodes_part2.ipynb +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -214,6 +214,7 @@ "import espressomd.observables\n", "import espressomd.accumulators\n", "import espressomd.shapes\n", + "import espressomd.zn\n", "\n", "espressomd.assert_features(['WCA', 'ELECTROSTATICS'])\n", "rng = np.random.default_rng(42)\n", @@ -476,20 +477,20 @@ "source": [ "# SOLUTION CELL\n", "offset = LJ_SIGMA # avoid unfavorable overlap at close distance to the walls\n", - "init_part_btw_z1 = offset \n", + "init_part_btw_z1 = offset\n", "init_part_btw_z2 = box_l_z - offset\n", - "ion_pos = np.empty((3), dtype=float)\n", + "ion_pos = np.empty((3,), dtype=float)\n", "\n", - "for i in range (N_IONPAIRS):\n", - " ion_pos[0] = rng.random(1) * system.box_l[0]\n", - " ion_pos[1] = rng.random(1) * system.box_l[1]\n", - " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", + "for i in range(N_IONPAIRS):\n", + " ion_pos[0] = rng.random(1)[0] * system.box_l[0]\n", + " ion_pos[1] = rng.random(1)[0] * system.box_l[1]\n", + " ion_pos[2] = rng.random(1)[0] * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", " system.part.add(pos=ion_pos, type=types[\"Cation\"], q=charges[\"Cation\"])\n", - " \n", - "for i in range (N_IONPAIRS):\n", - " ion_pos[0] = rng.random(1) * system.box_l[0]\n", - " ion_pos[1] = rng.random(1) * system.box_l[1]\n", - " ion_pos[2] = rng.random(1) * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", + "\n", + "for i in range(N_IONPAIRS):\n", + " ion_pos[0] = rng.random(1)[0] * system.box_l[0]\n", + " ion_pos[1] = rng.random(1)[0] * system.box_l[1]\n", + " ion_pos[2] = rng.random(1)[0] * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1\n", " system.part.add(pos=ion_pos, type=types[\"Anion\"], q=charges[\"Anion\"])" ] }, @@ -1008,7 +1009,40 @@ "\n", "With the above knowledge, we can now assess the \n", "differential capacitance of the system, by changing the applied voltage\n", - "difference and determining the corresponding surface charge density." + "difference and determining the corresponding surface charge density.\n", + "We will use the ZnDraw visualizer to visualize our system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef5d908f-a7f0-4845-9555-34822128c141", + "metadata": {}, + "outputs": [], + "source": [ + "color = {\n", + " types[\"Cation\"]: \"#ff0000\", #red\n", + " types[\"Anion\"]: \"#030ffc\" #blue\n", + " }\n", + " \n", + "radii = {\n", + " types[\"Cation\"]: LJ_SIGMA,\n", + " types[\"Anion\"]: LJ_SIGMA\n", + " }\n", + "\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", + "#vis.draw_constraints([floor, ceiling])\n", + "\n", + "#note: The system displayed is the enlarged system as you can see the ELG-GAP along the nonperiodic direction. The particles are only in the smaller subsystem.\n", + "#note: The particles are shown bigger for visualization purpose." + ] + }, + { + "cell_type": "markdown", + "id": "21da58e8-e666-4d06-8538-a8d6af521148", + "metadata": {}, + "source": [ + "Do the sampling from high to low potential:" ] }, { @@ -1029,7 +1063,9 @@ "for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):\n", " system.auto_update_accumulators.clear()\n", " system.electrostatics.solver = setup_electrostatic_solver(potential_diff)\n", - " system.integrator.run(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE)\n", + " for i in range(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE//50):\n", + " system.integrator.run(50)\n", + " vis.update()\n", " sigmas = []\n", "\n", " for tm in range(N_SAMPLES_CAP):\n", @@ -1039,7 +1075,9 @@ " system.auto_update_accumulators.clear()\n", " system.auto_update_accumulators.add(density_accumulator_cation)\n", " system.auto_update_accumulators.add(density_accumulator_anion)\n", - " system.integrator.run(STEPS_PER_SAMPLE)\n", + " for j in range(int(STEPS_PER_SAMPLE/50)):\n", + " system.integrator.run(50)\n", + " vis.update()\n", "\n", " cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]\n", " anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]\n", diff --git a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb index d65111cce2..2bd4d2ad56 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb @@ -257,6 +257,7 @@ "import espressomd.magnetostatics\n", "import espressomd.cluster_analysis\n", "import espressomd.pair_criteria\n", + "import espressomd.zn\n", "\n", "espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])\n", "\n", @@ -607,7 +608,7 @@ "id": "a6c4c46b", "metadata": {}, "source": [ - "## Sampling" + "## Sampling and cluster analysis" ] }, { @@ -615,6 +616,8 @@ "id": "af0906f5", "metadata": {}, "source": [ + "To quantify the number of clusters and their respective sizes, we now want to perform a cluster analysis.\n", + "For that we can use ESPREsSo's [cluster analysis class](https://espressomd.github.io/doc/analysis.html?highlight=cluster%20analysis#cluster-analysis).\n", "The system will be sampled over 100 loops." ] }, @@ -628,159 +631,6 @@ "LOOPS = 100" ] }, - { - "cell_type": "markdown", - "id": "516c2f33", - "metadata": {}, - "source": [ - "As the system is two dimensional, we can simply do a scatter plot to get a visual representation of a system state. To get a better insight of how a ferrofluid system develops during time we will create a video of the development of our system during the sampling. If you only want to sample the system simply go to [Sampling without animation](#Sampling-without-animation)" - ] - }, - { - "cell_type": "markdown", - "id": "60dc1653", - "metadata": {}, - "source": [ - "### Sampling with animation" - ] - }, - { - "cell_type": "markdown", - "id": "5f7fe3a7", - "metadata": {}, - "source": [ - "To get an animation of the system development we have to create a function which will save the video and embed it in an html string." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ed0ee691", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.animation as animation\n", - "import tempfile\n", - "import base64\n", - "\n", - "VIDEO_TAG = \"\"\"\"\"\"\n", - "\n", - "\n", - "def anim_to_html(anim):\n", - " if not hasattr(anim, '_encoded_video'):\n", - " with tempfile.NamedTemporaryFile(suffix='.mp4') as f:\n", - " anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])\n", - " with open(f.name, \"rb\") as g:\n", - " video = g.read()\n", - " anim._encoded_video = base64.b64encode(video).decode('ascii')\n", - " plt.close(anim._fig)\n", - " return VIDEO_TAG.format(anim._encoded_video)\n", - "\n", - "\n", - "animation.Animation._repr_html_ = anim_to_html\n", - "\n", - "\n", - "def init():\n", - " # Set x and y range\n", - " ax.set_ylim(0, BOX_SIZE)\n", - " ax.set_xlim(0, BOX_SIZE)\n", - " x_data, y_data = [], []\n", - " part.set_data(x_data, y_data)\n", - " return part," - ] - }, - { - "cell_type": "markdown", - "id": "89d9844a", - "metadata": {}, - "source": [ - "## Exercise:\n", - "\n", - "In the following an animation loop is defined, however it is incomplete.\n", - "Extend the code such that in every loop the system is integrated for 100 steps.\n", - "Afterwards `x_data` and `y_data` have to be populated by the folded $x$- and $y$- positions of the particles.\n", - "\n", - "(You may copy and paste the incomplete code template to the empty cell below.)\n", - "\n", - "```python\n", - "def run(i):\n", - " # < excercise >\n", - " \n", - " # Save current system state as a plot\n", - " x_data, y_data = # < excercise >\n", - " ax.figure.canvas.draw()\n", - " part.set_data(x_data, y_data)\n", - " print(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', end='\\r')\n", - " return part,\n", - "```" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "91004d62", - "metadata": {}, - "outputs": [], - "source": [ - "# SOLUTION CELL\n", - "def run(i):\n", - " system.integrator.run(100)\n", - "\n", - " # Save current system state as a plot\n", - " x_data, y_data = system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1]\n", - " ax.figure.canvas.draw()\n", - " part.set_data(x_data, y_data)\n", - " print(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', end='\\r')\n", - " return part," - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "95caca0c", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "id": "0f7e9093", - "metadata": {}, - "source": [ - "Now we use the animation class of matplotlib to save snapshots of the system as frames of a video which is then displayed after the sampling is finished. Between two frames are 100 integration steps.\n", - "\n", - "In the video chain-like and ring-like clusters should be visible, as well as some isolated monomers." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3d11e830", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(10, 10))\n", - "part, = ax.plot([], [], 'o')\n", - "\n", - "animation.FuncAnimation(fig, run, frames=LOOPS, blit=True, interval=0, repeat=False, init_func=init)" - ] - }, - { - "cell_type": "markdown", - "id": "a7ab793c", - "metadata": {}, - "source": [ - "## Cluster analysis\n", - "\n", - "To quantify the number of clusters and their respective sizes, we now want to perform a cluster analysis.\n", - "For that we can use ESPREsSo's [cluster analysis class](https://espressomd.github.io/doc/analysis.html?highlight=cluster%20analysis#cluster-analysis)." - ] - }, { "cell_type": "markdown", "id": "2f6d201a", @@ -847,15 +697,19 @@ "id": "6b8d8e7f", "metadata": {}, "source": [ - "### Sampling without animation" + "We initialize an instance of the ZnDraw visualizer to animate the sampling. As our system is constrained to 2D, the particles are on a single face of the cubic system box." ] }, { - "cell_type": "markdown", + "cell_type": "code", "id": "2098e033", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "The following code just samples the system and does a cluster analysis every loops (100 by default) simulation steps." + "color = {0: \"#7fc454\"} \n", + "radii = {0: LJ_SIGMA/2}\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)" ] }, { @@ -878,6 +732,7 @@ " for c in cluster_structure.clusters:\n", " cluster_sizes.append(# < excercise >)\n", " system.integrator.run(100)\n", + " vis.update()\n", "```" ] }, @@ -899,7 +754,8 @@ " n_clusters.append(len(cluster_structure.clusters))\n", " for c in cluster_structure.clusters:\n", " cluster_sizes.append(c[1].size())\n", - " system.integrator.run(100)" + " system.integrator.run(100)\n", + " vis.update()" ] }, { @@ -915,7 +771,7 @@ "id": "0042827e", "metadata": {}, "source": [ - "You may want to get a visualization of the current state of the system. For that we plot the particle positions folded to the simulation box using matplotlib." + "You may want to get a 2D visualization of the current state of the system. For that we plot the particle positions folded to the simulation box using matplotlib." ] }, { diff --git a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb index 453d5f9e97..e13eb759f9 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb @@ -59,6 +59,7 @@ "source": [ "import espressomd\n", "import espressomd.magnetostatics\n", + "import espressomd.zn\n", "\n", "espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])\n", "\n", @@ -305,95 +306,27 @@ "plt.show()" ] }, - { - "cell_type": "markdown", - "id": "2f52d782", - "metadata": {}, - "source": [ - "## Video of the development of the system" - ] - }, { "cell_type": "markdown", "id": "eb06ab55", "metadata": {}, "source": [ - "You may want to get an insight of how the system develops in time. Thus we now create a function which will save a video and embed it in an html string to create a video of the systems development " + "You may want to get an insight of how the system develops in time. We use the ZnDraw visualiser:" ] }, { "cell_type": "code", "execution_count": null, - "id": "e25ba03b", + "id": "a6807ebd-2f97-476c-94e1-ea5fe1d383a4", "metadata": {}, "outputs": [], "source": [ - "import matplotlib.pyplot as plt\n", - "import matplotlib.animation as animation\n", - "import tempfile\n", - "import base64\n", - "\n", - "VIDEO_TAG = \"\"\"\"\"\"\n", - "\n", - "\n", - "def anim_to_html(anim):\n", - " if not hasattr(anim, '_encoded_video'):\n", - " with tempfile.NamedTemporaryFile(suffix='.mp4') as f:\n", - " anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])\n", - " with open(f.name, \"rb\") as g:\n", - " video = g.read()\n", - " anim._encoded_video = base64.b64encode(video).decode('ascii')\n", - " plt.close(anim._fig)\n", - " return VIDEO_TAG.format(anim._encoded_video)\n", - "\n", - "\n", - "animation.Animation._repr_html_ = anim_to_html\n", - "\n", - "\n", - "def init():\n", - " # Set x and y range\n", - " ax.set_ylim(0, box_size)\n", - " ax.set_xlim(0, box_size)\n", - " xdata, ydata = [], []\n", - " part.set_data(xdata, ydata)\n", - " return part,\n", - "\n", - "\n", - "def run(i):\n", + "color = {0: \"#7fc454\"} \n", + "radii = {0: LJ_SIGMA/2}\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", + "for _ in range(100):\n", " system.integrator.run(50)\n", - "\n", - " # Save current system state as a plot\n", - " xdata, ydata = particles.pos_folded[:, 0], particles.pos_folded[:, 1]\n", - " ax.figure.canvas.draw()\n", - " part.set_data(xdata, ydata)\n", - " print(f'progress: {(i + 1):3.0f}%', end='\\r')\n", - " return part," - ] - }, - { - "cell_type": "markdown", - "id": "25b6af88", - "metadata": {}, - "source": [ - "We now can start the sampling over the animation class of matplotlib" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8c9fe146", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(10, 10))\n", - "part, = ax.plot([], [], 'o')\n", - "\n", - "animation.FuncAnimation(fig, run, frames=100, blit=True, interval=0, repeat=False, init_func=init)" + " vis.update()" ] }, { @@ -401,7 +334,7 @@ "id": "278cfc7f", "metadata": {}, "source": [ - "In the visualization video we can see that the single chains break and connect to each other during time. Also some monomers are present which break from and connect to chains. If you want to have some more frames, i.e. a longer video, just adjust the frames parameter in FuncAnimation." + "In the visualization video we can see that the single chains break and connect to each other during time. Also some monomers are present which break from and connect to chains." ] }, { diff --git a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb index b7f896f46a..8bce1733db 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb @@ -134,6 +134,7 @@ "source": [ "import espressomd\n", "import espressomd.magnetostatics\n", + "import espressomd.zn\n", "\n", "espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])\n", "\n", @@ -536,6 +537,14 @@ "As in **part II** we use the same system for every value of the Langevin parameter $\\alpha$. Thus we use that the system is already pre-equilibrated from the previous run so we save some equilibration time. For scientific purposes one would use a new system for every value for the Langevin parameter to ensure that the systems are independent and no correlation effects are measured. Also one would perform more than just one simulation for each value of $\\alpha$ to increase the precision of the results." ] }, + { + "cell_type": "markdown", + "id": "395a4fc0-baaf-48d6-8b88-b8ba42408bfe", + "metadata": {}, + "source": [ + "Additionally, we will visualize the sampling of the magnetization curve using the ZnDraw visualizer. Note, that the magnetic field is increasing during the video." + ] + }, { "cell_type": "code", "execution_count": null, @@ -579,7 +588,8 @@ " magnetizations[ndx] = np.mean(magn_temp)\n", "\n", " # remove constraint\n", - " system.constraints.clear()" + " if not alpha == alphas[-1]:\n", + " system.constraints.clear()" ] }, { @@ -668,6 +678,39 @@ "At sufficiently high volume fractions $\\phi$ and dipolar interaction parameters $\\lambda$ these clusters can be so rigid, that simulations with normal methods are impossible as the relaxation times exceeds normal simulation times by far, resulting in strongly correlated configurations and thus measurements." ] }, + { + "cell_type": "markdown", + "id": "d2333fe6-c427-48b6-92f0-d946f4084b30", + "metadata": {}, + "source": [ + "## Visualization" + ] + }, + { + "cell_type": "markdown", + "id": "7fd720f4-858e-4d2e-8d4a-23644ab996e3", + "metadata": {}, + "source": [ + "Finally, we use ZnDraw to visualize the evolving system for high $\\alpha$. We just need to continue running our system as it is already warmed up for high magnetic field." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e1d8ade-efa9-4d73-a328-27b0e63993b4", + "metadata": {}, + "outputs": [], + "source": [ + "# initializing zndraw visualizer\n", + "color = {0: \"#7fc454\"} \n", + "radii = {0: lj_sigma/2}\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", + "\n", + "for _ in range(400):\n", + " system.integrator.run(5)\n", + " vis.update()" + ] + }, { "cell_type": "markdown", "id": "81e06b9b", diff --git a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb index 804725d179..5f6471c40b 100644 --- a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb +++ b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb @@ -159,14 +159,6 @@ "## 2. Simulating Brownian motion" ] }, - { - "cell_type": "markdown", - "id": "c38456fa", - "metadata": {}, - "source": [ - "We will simulate the diffusion of a single particle that is coupled to an implicit solvent." - ] - }, { "cell_type": "code", "execution_count": null, @@ -181,6 +173,7 @@ "import espressomd\n", "import espressomd.accumulators\n", "import espressomd.observables\n", + "import espressomd.zn\n", "\n", "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n", "\n", @@ -188,15 +181,76 @@ "KT = 1.1\n", "STEPS = 1000000\n", "\n", + "gammas = [1.0, 2.0, 4.0, 10.0]" + ] + }, + { + "cell_type": "markdown", + "id": "c38456fa", + "metadata": {}, + "source": [ + "We will simulate the diffusion of a single particle that is coupled to an implicit solvent. First, we have to setup our system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45a6ad79-fae4-4bd7-b041-012b193b83bf", + "metadata": {}, + "outputs": [], + "source": [ "# System setup\n", "system = espressomd.System(box_l=[16] * 3)\n", "system.time_step = 0.01\n", "system.cell_system.skin = 0.4\n", "\n", - "particle = system.part.add(pos=[0, 0, 0])\n", + "particle = system.part.add(pos=[0, 0, 0])" + ] + }, + { + "cell_type": "markdown", + "id": "3d594176-bd64-4ce8-9e25-f0efefbf996c", + "metadata": {}, + "source": [ + "In order to get an idea of the Brownian motion, we will visualize our system evolving for a few steps with the ZnDraw visualizer:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80788a7f-4bf2-4fa6-9292-5de8b2a926ff", + "metadata": {}, + "outputs": [], + "source": [ + "system.thermostat.set_langevin(kT=KT, gamma=gammas[0], seed=42)\n", "\n", + "# Setup ZnDraw visualizer\n", + "color = {0: \"#7fc454\"} \n", + "radii = {0: 1}\n", + "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", + "\n", + "system.integrator.run(1000)\n", + "for _ in range(200):\n", + " system.integrator.run(50)\n", + " vis.update()" + ] + }, + { + "cell_type": "markdown", + "id": "933cea77-4bca-47a6-8a82-cc13d07ec673", + "metadata": {}, + "source": [ + "Now we can do the actual simulation where we simulate for different friction coefficients and sample our observables:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf997a0a", + "metadata": {}, + "outputs": [], + "source": [ "# Run for different friction coefficients\n", - "gammas = [1.0, 2.0, 4.0, 10.0]\n", "tau_results = []\n", "msd_results = []\n", "vacf_results = []\n", diff --git a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb index 1d91441855..a72c35294e 100644 --- a/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb +++ b/doc/tutorials/raspberry_electrophoresis/raspberry_electrophoresis.ipynb @@ -71,6 +71,7 @@ "import espressomd.interactions\n", "import espressomd.electrostatics\n", "import espressomd.lb\n", + "import espressomd.zn\n", "\n", "import sys\n", "import tqdm\n", @@ -683,6 +684,56 @@ "system.thermostat.set_lb(LB_fluid=lb, seed=123, gamma=20.0)" ] }, + { + "cell_type": "markdown", + "id": "a7388457-8c1b-4839-bfc6-e34778ae9511", + "metadata": {}, + "source": [ + "## Visualization" + ] + }, + { + "cell_type": "markdown", + "id": "f3deb299-de31-4901-8f65-0ff06112232c", + "metadata": {}, + "source": [ + "We use the ZnDraw visualizer and run our simulation for a few steps in order to get an idea what is happening in our system:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c24f4ee4-2a6a-46ff-ba89-4f6c7f7fb227", + "metadata": {}, + "outputs": [], + "source": [ + "# Initializing ZnDraw visualizer\n", + "colors = {\n", + " TYPE_CENTRAL: \"#fca000\", #orange\n", + " TYPE_SURFACE: \"#7fc454\", #green\n", + " TYPE_CATIONS: \"#ff0000\", #red\n", + " TYPE_ANIONS: \"#030ffc\" #blue\n", + " }\n", + "radii = {\n", + " TYPE_CENTRAL: sig_ss/2,\n", + " TYPE_SURFACE: sig_ss/2,\n", + " TYPE_CATIONS: sig_ss/2,\n", + " TYPE_ANIONS: sig_ss/2\n", + " }\n", + "arrows_config = {'colormap': [[-0.5, 0.9, 0.5], [-0.0, 0.9, 0.5]],\n", + " 'normalize': True,\n", + " 'colorrange': [0, 1],\n", + " 'scale_vector_thickness': True,\n", + " 'opacity': 1.0}\n", + "\n", + "lbfield = espressomd.zn.LBField(system, step_x=2, step_y=2, step_z=5, scale=0.3)\n", + "vis = espressomd.zn.Visualizer(system, colors=colors, radii=radii, folded=True, vector_field=lbfield)\n", + "\n", + "for _ in range(50):\n", + " system.integrator.run(10)\n", + " vis.update()" + ] + }, { "cell_type": "markdown", "id": "f848c501", @@ -709,6 +760,7 @@ "with open('posVsTime.dat', 'w') as f: # file where the raspberry trajectory will be written to\n", " for _ in tqdm.tqdm(range(num_iterations)):\n", " system.integrator.run(num_steps_per_iteration)\n", + " vis.update()\n", " pos = central_part.pos - initial_pos\n", " f.write(f\"{system.time:.2f} {pos[0]:.4f} {pos[1]:.4f} {pos[2]:.4f}\\n\")\n", "\n", From efacee2da2bc083b56911d3e310fd45b69cc5398 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 23 Aug 2024 19:48:16 +0200 Subject: [PATCH 2/2] Edit and format tutorials --- doc/tutorials/constant_pH/constant_pH.ipynb | 4 +- .../electrodes/electrodes_part2.ipynb | 7 +- .../ferrofluid/ferrofluid_part1.ipynb | 74 +++++++++---------- .../ferrofluid/ferrofluid_part2.ipynb | 22 +++--- .../ferrofluid/ferrofluid_part3.ipynb | 8 +- .../langevin_dynamics/langevin_dynamics.ipynb | 2 +- .../widom_insertion/widom_insertion.ipynb | 12 ++- 7 files changed, 69 insertions(+), 60 deletions(-) diff --git a/doc/tutorials/constant_pH/constant_pH.ipynb b/doc/tutorials/constant_pH/constant_pH.ipynb index 8dabb3dc2b..0c52dde352 100644 --- a/doc/tutorials/constant_pH/constant_pH.ipynb +++ b/doc/tutorials/constant_pH/constant_pH.ipynb @@ -745,7 +745,7 @@ "id": "cd417295-bab3-46a5-a574-fc60d06371a5", "metadata": {}, "source": [ - "We will initialize the zndraw visualizer to visualize the simulation for increasing pH value:" + "We will initialize ZnDraw to visualize the simulation for increasing pH values:" ] }, { @@ -753,7 +753,7 @@ "execution_count": null, "id": "04dd303c", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ diff --git a/doc/tutorials/electrodes/electrodes_part2.ipynb b/doc/tutorials/electrodes/electrodes_part2.ipynb index 1aaf34d86e..e7007c62e3 100644 --- a/doc/tutorials/electrodes/electrodes_part2.ipynb +++ b/doc/tutorials/electrodes/electrodes_part2.ipynb @@ -1010,7 +1010,7 @@ "With the above knowledge, we can now assess the \n", "differential capacitance of the system, by changing the applied voltage\n", "difference and determining the corresponding surface charge density.\n", - "We will use the ZnDraw visualizer to visualize our system:" + "We will use ZnDraw to visualize our system:" ] }, { @@ -1033,8 +1033,9 @@ "vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)\n", "#vis.draw_constraints([floor, ceiling])\n", "\n", - "#note: The system displayed is the enlarged system as you can see the ELG-GAP along the nonperiodic direction. The particles are only in the smaller subsystem.\n", - "#note: The particles are shown bigger for visualization purpose." + "# note: you may need to zoom out since the ELC gap region takes a significant portion of\n", + "# the box along the non-periodic direction. The particles are only in the smaller subsystem.\n", + "# note: The particles are shown bigger for visualization purpose." ] }, { diff --git a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb index 2bd4d2ad56..72d116a517 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part1.ipynb @@ -15,15 +15,13 @@ "source": [ "## Table of Contents\n", "1. [Introduction](#Introduction)\n", - "2. [The Model](#The-Model)\n", + "2. [The model](#The-model)\n", "3. [Structure of this tutorial](#Structure-of-this-tutorial)\n", - "4. [Compiling ESPResSo for this Tutorial](#Compiling-ESPResSo-for-this-Tutorial)\n", - "5. [A Monolayer-Ferrofluid System in ESPResSo](#A-Monolayer-Ferrofluid-System-in-ESPResSo)\n", + "4. [Compiling ESPResSo for this tutorial](#Compiling-ESPResSo-for-this-tutorial)\n", + "5. [A monolayer-ferrofluid system in ESPResSo](#A-monolayer-ferrofluid-system-in-ESPResSo)\n", " 1. [Setup](#Setup)\n", - " 2. [Sampling](#Sampling)\n", - " 3. [Sampling with animation](#Sampling-with-animation)\n", - " 4. [Sampling without animation](#Sampling-without-animation)\n", - " 5. [Cluster distribution](#Cluster-distribution)" + " 2. [Sampling and cluster analysis](#Sampling-and-cluster-analysis)\n", + " 3. [Cluster distribution](#Cluster-distribution)" ] }, { @@ -54,10 +52,8 @@ "metadata": {}, "source": [ "
\n", - "\n", - "
\n", - "
Figure 1: Schematic representation of electrostatically stabilization (picture top) and steric stabilization (picture bottom) [3]
\n", - "
\n", + "\n", + "
Figure 1: Schematic representation of electrostatically stabilization (picture top) and steric stabilization (picture bottom) [3]
\n", "
" ] }, @@ -67,10 +63,8 @@ "metadata": {}, "source": [ "
\n", - "ferrofluid on glass plate under which a strong magnet is placed\n", - "
\n", - "
Figure 2: Real Ferrofluid exposed to an external magnetic field (neodymium magnet) [4]
\n", - "
\n", + "ferrofluid on glass plate under which a strong magnet is placed\n", + "
Figure 2: Real ferrofluid exposed to an external magnetic field (neodymium magnet) [4]
\n", "
" ] }, @@ -79,7 +73,7 @@ "id": "095339c1", "metadata": {}, "source": [ - "## The Model" + "## The model" ] }, { @@ -124,7 +118,7 @@ " \\lambda = \\frac{\\tilde{u}_{\\text{DD}}}{u_T} = \\gamma \\frac{\\mu^2}{k_{\\text{B}}T\\sigma^3}\n", "\\end{equation}\n", "\n", - "where $u_\\mathrm{T} = k_{\\text{B}}T$ is the thermal energy and $\\tilde{u}_{DD}$ is the absolute value of the dipole-dipole interaction energy at close contact (cc) and head-to-tail configuration (htt) (see figure 4) per particle, i.e. in formulas it reads\n", + "where $u_\\mathrm{T} = k_{\\text{B}}T$ is the thermal energy and $\\tilde{u}_{DD}$ is the absolute value of the dipole-dipole interaction energy at close contact (cc) and head-to-tail configuration (htt) (see figure 3) per particle, i.e. in formulas it reads\n", "\n", "\\begin{equation}\n", " \\tilde{u}_{\\text{DD}} = \\frac{ \\left| u_{\\text{DD}}^{\\text{htt, cc}} \\right| }{2}\n", @@ -143,11 +137,9 @@ "id": "62b95d9c", "metadata": {}, "source": [ - "
\n", - "schematic representation of head to tail configuration\n", - "
\n", - "
Figure 4: Schematic representation of the head-to-tail configuration of two magnetic particles at close contact.
\n", - "
\n", + "
\n", + "schematic representation of head to tail configuration\n", + "
Figure 3: Schematic representation of the head-to-tail configuration of two magnetic particles at close contact.
\n", "
" ] }, @@ -182,7 +174,7 @@ "id": "b6d5a65e", "metadata": {}, "source": [ - "## Compiling ESPResSo for this Tutorial" + "## Compiling ESPResSo for this tutorial" ] }, { @@ -218,7 +210,7 @@ "id": "a0bffbeb", "metadata": {}, "source": [ - "## A Monolayer-Ferrofluid System in ESPResSo" + "## A monolayer-ferrofluid system in ESPResSo" ] }, { @@ -324,7 +316,8 @@ "id": "136aa645", "metadata": {}, "source": [ - "## Exercise:\n", + "**Exercise:**\n", + "\n", "How large does `BOX_SIZE` have to be for a system of `N_PART` particles with a volume (area) fraction `PHI`?\n", "Define `BOX_SIZE`." ] @@ -384,7 +377,8 @@ "outputs": [], "source": [ "# Lennard-Jones interaction\n", - "system.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA, cutoff=LJ_CUT, shift=\"auto\")" + "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", + " epsilon=LJ_EPSILON, sigma=LJ_SIGMA, cutoff=LJ_CUT, shift=\"auto\")" ] }, { @@ -394,7 +388,7 @@ "source": [ "Now we generate random positions and orientations of the particles and their dipole moments. \n", "\n", - "**Hint:**\n", + "Hint:\n", "It should be noted that we seed the random number generator of numpy. Thus the initial configuration of our system is the same every time this script will be executed. You can change it to another one to simulate with a different initial configuration." ] }, @@ -403,9 +397,11 @@ "id": "c49729f7", "metadata": {}, "source": [ - "## Exercise:\n", + "**Exercise:**\n", + "\n", "How does one set up randomly oriented dipole moments?\n", - "*Hint*: Think of the way that different methods could introduce a bias in the distribution of the orientations.\n", + "\n", + "Hint: Think of the way that different methods could introduce a bias in the distribution of the orientations.\n", "\n", "Create a variable `dip` as a `N_PART x 3` numpy array, which contains the randomly distributed dipole moments." ] @@ -523,7 +519,7 @@ "For the simulation of our system we choose the velocity Verlet integrator.\n", "After that we set up the thermostat which is, in our case, a Langevin thermostat to simulate in an NVT ensemble.\n", "\n", - "**Hint:**\n", + "Hint:\n", "It should be noted that we seed the Langevin thermostat, thus the time evolution of the system is partly predefined.\n", "Partly because of the numeric accuracy and the automatic tuning algorithms of Dipolar P3M and DLC where the resulting parameters are slightly different every time.\n", "You can change the seed to get a guaranteed different time evolution." @@ -562,7 +558,7 @@ "execution_count": null, "id": "fd12b1f6", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ @@ -592,7 +588,7 @@ "execution_count": null, "id": "95c4a8b6", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ @@ -608,7 +604,7 @@ "id": "a6c4c46b", "metadata": {}, "source": [ - "## Sampling and cluster analysis" + "### Sampling and cluster analysis" ] }, { @@ -617,7 +613,7 @@ "metadata": {}, "source": [ "To quantify the number of clusters and their respective sizes, we now want to perform a cluster analysis.\n", - "For that we can use ESPREsSo's [cluster analysis class](https://espressomd.github.io/doc/analysis.html?highlight=cluster%20analysis#cluster-analysis).\n", + "For that we can use ESPResSo's [cluster analysis class](https://espressomd.github.io/doc/analysis.html#cluster-analysis).\n", "The system will be sampled over 100 loops." ] }, @@ -636,7 +632,7 @@ "id": "2f6d201a", "metadata": {}, "source": [ - "## Exercise:\n", + "**Exercise:**\n", "\n", "Setup a cluster analysis object (`ClusterStructure` class) and assign its instance to the variable `cluster_structure`.\n", "As criterion for the cluster analysis use a distance criterion where particles are assumed to be\n", @@ -717,7 +713,7 @@ "id": "98b99dd6", "metadata": {}, "source": [ - "## Exercise:\n", + "**Exercise:**\n", "\n", "Write an integration loop which runs a cluster analysis on the system, saving the number of clusters `n_clusters` and the size distribution `cluster_sizes`.\n", "Take the following as a starting point:\n", @@ -803,7 +799,7 @@ "id": "4c0828ea", "metadata": {}, "source": [ - "## Cluster distribution" + "### Cluster distribution" ] }, { @@ -819,12 +815,12 @@ "id": "1c6be57c", "metadata": {}, "source": [ - "## Exercise:\n", + "**Exercise:**\n", "\n", "Use `numpy` to calculate a histogram of the cluster sizes and assign it to the variable `size_dist`.\n", "Take only clusters up to a size of 19 particles into account.\n", "\n", - "*Hint*: In order not to count clusters with size 20 or more, one may include an additional bin containing these.\n", + "Hint: In order not to count clusters with size 20 or more, one may include an additional bin containing these.\n", "The reason for that is that `numpy` defines the histogram bins as half-open intervals with the open border at the upper bin edge.\n", "Consequently clusters of larger sizes are attributed to the last bin.\n", "By not using the last bin in the plot below, these clusters can effectively be neglected." diff --git a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb index e13eb759f9..80dcc881bb 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part2.ipynb @@ -167,7 +167,8 @@ "pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))\n", "\n", "# Add particles\n", - "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)], dip=dip, fix=N_PART * [(False, False, True)])\n", + "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)],\n", + " dip=dip, fix=N_PART * [(False, False, True)])\n", "\n", "# Remove overlap between particles by means of the steepest descent method\n", "MASS = 1.0\n", @@ -270,7 +271,7 @@ "execution_count": null, "id": "7750295c", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ @@ -471,7 +472,8 @@ "pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))\n", "\n", "# Add particles\n", - "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)], dip=dip, fix=N_PART * [(False, False, True)])\n", + "particles = system.part.add(pos=pos, rotation=N_PART * [(True, True, True)],\n", + " dip=dip, fix=N_PART * [(False, False, True)])\n", "\n", "# Remove overlap between particles by means of the steepest descent method\n", "system.integrator.set_steepest_descent(f_max=0, gamma=0.1, max_displacement=0.05)\n", @@ -803,7 +805,7 @@ "outputs": [], "source": [ "def dL_dB(alpha):\n", - " return (1. / (alpha**2.) - 1. / ((np.sinh(alpha))**2.)) * dipm / (KT)" + " return (1. / (alpha**2) - 1. / np.sinh(alpha)**2) * dipm / (KT)" ] }, { @@ -884,8 +886,8 @@ "plt.xlabel(r'$\\alpha$', fontsize=20)\n", "plt.ylabel(r'$M^*$', fontsize=20)\n", "plt.plot(x, L(x), label='Langevin function')\n", - "plt.plot(x, magnetization_approx_perp(x), label='q2D approximation $\\perp$')\n", - "plt.plot(x, magnetization_approx_para(x), label='q2D approximation $\\parallel$')\n", + "plt.plot(x, magnetization_approx_perp(x), label=r'q2D approximation $\\perp$')\n", + "plt.plot(x, magnetization_approx_para(x), label=r'q2D approximation $\\parallel$')\n", "# < exercise >\n", "plt.legend(fontsize=20)\n", "plt.show()\n", @@ -907,10 +909,10 @@ "plt.xlabel(r'$\\alpha$', fontsize=20)\n", "plt.ylabel(r'$M^*$', fontsize=20)\n", "plt.plot(x, L(x), label='Langevin function')\n", - "plt.plot(x, magnetization_approx_perp(x), label='q2D approximation $\\perp$')\n", - "plt.plot(x, magnetization_approx_para(x), label='q2D approximation $\\parallel$')\n", - "plt.plot(alphas, magnetization_perp / N_PART, 'o', label='simulation results $\\perp$')\n", - "plt.plot(alphas, magnetization_para / N_PART, 'o', label='simulation results $\\parallel$')\n", + "plt.plot(x, magnetization_approx_perp(x), label=r'q2D approximation $\\perp$')\n", + "plt.plot(x, magnetization_approx_para(x), label=r'q2D approximation $\\parallel$')\n", + "plt.plot(alphas, magnetization_perp / N_PART, 'o', label=r'simulation results $\\perp$')\n", + "plt.plot(alphas, magnetization_para / N_PART, 'o', label=r'simulation results $\\parallel$')\n", "plt.legend(fontsize=20)\n", "plt.show()" ] diff --git a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb index 8bce1733db..b3be7c3c32 100644 --- a/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb +++ b/doc/tutorials/ferrofluid/ferrofluid_part3.ipynb @@ -273,7 +273,7 @@ "execution_count": null, "id": "89966f9b", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ @@ -334,11 +334,11 @@ "execution_count": null, "id": "bd24f8eb", "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [], "source": [ - "# initialize array for hold the sampled dipole moments\n", + "# initialize array that will hold the sampled dipole moments\n", "dipms = np.full((loops, 3), np.nan)\n", "\n", "# sample dipole moment\n", @@ -542,7 +542,7 @@ "id": "395a4fc0-baaf-48d6-8b88-b8ba42408bfe", "metadata": {}, "source": [ - "Additionally, we will visualize the sampling of the magnetization curve using the ZnDraw visualizer. Note, that the magnetic field is increasing during the video." + "Additionally, we will visualize the sampling of the magnetization curve using ZnDraw. Note, that the magnetic field is increasing during the video." ] }, { diff --git a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb index 5f6471c40b..6c602341e1 100644 --- a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb +++ b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb @@ -240,7 +240,7 @@ "id": "933cea77-4bca-47a6-8a82-cc13d07ec673", "metadata": {}, "source": [ - "Now we can do the actual simulation where we simulate for different friction coefficients and sample our observables:" + "Now we can do the actual production run where we simulate for different friction coefficients and sample our observables:" ] }, { diff --git a/doc/tutorials/widom_insertion/widom_insertion.ipynb b/doc/tutorials/widom_insertion/widom_insertion.ipynb index 04e9140afd..cf06a7c507 100644 --- a/doc/tutorials/widom_insertion/widom_insertion.ipynb +++ b/doc/tutorials/widom_insertion/widom_insertion.ipynb @@ -195,6 +195,7 @@ "Now we are ready to set up the system. Because we will need to rescale the simulation box, we will initially only set up the short-range interactions.\n", "\n", "**Exercise:**\n", + "\n", "* Set up a system with box length $10\\,\\sigma$, integrator time step $\\Delta t=0.01 \\,\\tau$ ($\\tau$ is the Lennard-Jones timescale, i.e. equal to 1 in the employed unit system) and Verlet skin width of $\\delta = 0.4\\sigma$\n", "* Add a total of ``N_ION_PAIRS`` ion pairs at random positions to the simulation box\n", "* Add a WCA interaction (purely repulsive Lennard-Jones interaction) between the particles\n", @@ -244,10 +245,12 @@ "metadata": {}, "source": [ "**Exercise:**\n", + "\n", "* Implement a function ``system_setup(c_salt_SI)`` that takes the desired salt concentration ``c_salt_SI`` in SI-units (mol/l) as an argument and rescales the simulation box accordingly\n", "* Afterwards, the function should also add the P3M actor to account for electrostatic interactions\n", "\n", "**Hints:**\n", + "\n", "* Use the prefactor ``PREF`` defined above to convert the concentration from SI-units to the unit system employed in the simulation\n", "* An accuracy of $10^{-3}$ for the P$^3$M should be enough for this exercise" ] @@ -290,9 +293,11 @@ "metadata": {}, "source": [ "**Exercise:**\n", + "\n", "* Implement a function ``warmup()`` that removes potential overlaps between particles with 10000 steps of the steepest descent integrator and performs a warmup integration with 100 loops of 100 simulation steps\n", "\n", "**Hints:**\n", + "\n", "* keep in mind that after the overlaps have been removed, the integrator should be switched\n", " back to Velocity-Verlet with the Langevin thermostat (using ``GAMMA`` and ``KT``)\n", " in order to perform a warmup integration\n", @@ -340,10 +345,11 @@ "metadata": {}, "source": [ "**Exercise:**\n", - "* Add the functionality to perform Widom insertion moves. Make sure that you always insert an ion pair in order to conserve the system electroneutrality.\n", "\n", + "* Add the functionality to perform Widom insertion moves. Make sure that you always insert an ion pair in order to conserve the system electroneutrality.\n", "\n", "**Hints:**\n", + "\n", "* Refer to the documentation to find out how to set up Widom particle insertion ([Widom Insertion](https://espressomd.github.io/doc/reaction_methods.html#widom-insertion-for-homogeneous-systems))" ] }, @@ -377,6 +383,7 @@ "metadata": {}, "source": [ "**Exercise:**\n", + "\n", "* Implement a function ``calculate_excess_chemical_potential(n_md_steps, n_mc_steps, sample_size)``\n", " that performs ``n_mc_steps`` Widom particle insertions every ``n_md_steps`` MD steps in a loop\n", " that repeats ``sample_size`` times and returns the chemical potential and its error as a tuple via\n", @@ -471,6 +478,7 @@ "a logarithmically scaled x-axis.\n", "\n", "**Exercise:**\n", + "\n", "* Explain the observed behaviour qualitatively\n", "* How do you expect the excess chemical potential to behave in the limit $c_{\\mathrm{salt}}\\rightarrow 0$ mol/l?" ] @@ -524,9 +532,11 @@ "$B \\simeq 1.0$ and $C = 0.2$ in reduced units for NaCl in an aqueous solution at 25°C.\n", "\n", "**Exercise:**\n", + "\n", "* Compare your simulation results with the Debye-Hückel theory and the Davies equation by plotting all three curves in a single plot\n", "\n", "**Hint:**\n", + "\n", "* re-use the code from the previous figure\n", "* create a logarithmic concentration range from $10^{-4}$ to $10^{0}$ mol/l with ``np.logspace(-4, 0.0, num=500, base=10)`` to plot the analytical solutions" ]