Skip to content

Commit

Permalink
tutorial 04: part 3: Green-Kubo diffusion coefficient
Browse files Browse the repository at this point in the history
  • Loading branch information
@jngrad
jngrad committed Oct 6, 2020
1 parent 15956f9 commit 8841251
Show file tree
Hide file tree
Showing 3 changed files with 215 additions and 40 deletions.
2 changes: 1 addition & 1 deletion doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -307,7 +307,7 @@
"source": [
"plt.figure(figsize=(10, 8))\n",
"plt.xlabel(r'$\\gamma$')\n",
"plt.ylabel('Diffusion coefficient')\n",
"plt.ylabel('Diffusion coefficient [$\\sigma^2/t$]')\n",
"gammas = np.linspace(0.9 * min(lb_gammas), 1.1 * max(lb_gammas), 50)\n",
"plt.plot(gammas, lbf.kT / gammas, '-', label=r'$k_BT\\gamma^{-1}$')\n",
"plt.plot(lb_gammas, diffusion_results, 'o', label='D')\n",
Expand Down
247 changes: 208 additions & 39 deletions doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,14 +48,15 @@
"\n",
"which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description.\n",
"\n",
"Furthermore we want to compute the diffusion constant of the polymer for different numbers of monomers.\n",
"For this purpose we can again use the multiple tau correlator. The following script computes the mean\n",
"squared displacement for the center of mass of the polymer as well as the average hydrodynamic radius\n",
"$R_h$, end-to-end distance $R_F$ and radius of gyration $R_g$.\n",
"The first task is to compute the average hydrodynamic radius $R_h$, end-to-end distance $R_F$\n",
"and radius of gyration $R_g$ for different polymer lengths.\n",
"\n",
"How do $R_h$, $R_g$, $R_F$ and the diffusion coefficient $D$ depend on the number of monomers?\n",
"You can refer to the Flory theory of polymers, and assume you are simulating a real polymer in a\n",
"good solvent, with Flory exponent $\\nu \\approx 0.588$."
"The second task is to estimate the polymer diffusion coefficient for different polymer lengths\n",
"using two methods:\n",
"* the center of mass mean squared displacement method (introduced in a previous part of this tutorial)\n",
"* the center of mass velocity autocorrelation method (also known as Green–Kubo method)\n",
"\n",
"For this purpose we can again use the multiple tau correlator."
]
},
{
Expand Down Expand Up @@ -97,8 +98,10 @@
"\n",
"N_MONOMERS = np.array([6, 8, 10])\n",
"\n",
"tau_results = []\n",
"msd_results = []\n",
"com_pos_tau_results = []\n",
"com_pos_msd_results = []\n",
"com_vel_tau_results = []\n",
"com_vel_acf_results = []\n",
"rh_results = []\n",
"rf_results = []\n",
"rg_results = []\n",
Expand Down Expand Up @@ -147,12 +150,19 @@
" system.integrator.run(1000)\n",
" logging.info(\"Warming up the system with LB fluid finished.\")\n",
"\n",
" # configure correlator\n",
" # configure MSD correlator\n",
" com_pos = espressomd.observables.ComPosition(ids=range(N))\n",
" correlator = espressomd.accumulators.Correlator(\n",
" com_pos_cor = espressomd.accumulators.Correlator(\n",
" obs1=com_pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n",
" corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n",
" system.auto_update_accumulators.add(correlator)\n",
" system.auto_update_accumulators.add(com_pos_cor)\n",
"\n",
" # configure Green-Kubo correlator\n",
" com_vel = espressomd.observables.ComVelocity(ids=range(N))\n",
" com_vel_cor = espressomd.accumulators.Correlator(\n",
" obs1=com_vel, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=10,\n",
" corr_operation=\"scalar_product\", compress1=\"discard1\")\n",
" system.auto_update_accumulators.add(com_vel_cor)\n",
"\n",
" logging.info(\"Sampling started.\")\n",
" rhs = np.zeros(LOOPS)\n",
Expand All @@ -172,16 +182,15 @@
" chain_start=0,\n",
" number_of_chains=1,\n",
" chain_length=N)[0]\n",
"\n",
" logging.info(\"Sampling finished.\")\n",
"\n",
" # store results\n",
" correlator.finalize()\n",
" corrdata = correlator.result()\n",
" tau = correlator.lag_times()\n",
" msd = np.sum(corrdata, axis=1)\n",
" tau_results.append(tau)\n",
" msd_results.append(msd)\n",
" com_pos_cor.finalize()\n",
" com_pos_tau_results.append(com_pos_cor.lag_times())\n",
" com_pos_msd_results.append(np.sum(com_pos_cor.result(), axis=1))\n",
" com_vel_cor.finalize()\n",
" com_vel_tau_results.append(com_vel_cor.lag_times())\n",
" com_vel_acf_results.append(com_vel_cor.result())\n",
" rh_results.append(rhs)\n",
" rf_results.append(rfs)\n",
" rg_results.append(rgs)\n",
Expand All @@ -195,8 +204,18 @@
"rh_results = np.array(rh_results)\n",
"rf_results = np.array(rf_results)\n",
"rg_results = np.array(rg_results)\n",
"tau_results = np.array(tau_results)\n",
"msd_results = np.reshape(msd_results, [len(N_MONOMERS),-1])"
"com_pos_tau_results = np.array(com_pos_tau_results)\n",
"com_pos_msd_results = np.reshape(com_pos_msd_results, [len(N_MONOMERS),-1])\n",
"com_vel_tau_results = np.array(com_vel_tau_results)\n",
"com_vel_acf_results = np.reshape(com_vel_acf_results, [len(N_MONOMERS),-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will calculate the means of time series with error bars obtained from\n",
"the correlation-corrected standard error of the mean [8,9]."
]
},
{
Expand Down Expand Up @@ -237,9 +256,10 @@
" '''\n",
" Calculate the mean and the correlation-corrected standard error\n",
" of the mean of time series by integrating the autocorrelation\n",
" function. Due to the short simulation length, it is not possible\n",
" to fit an exponential to the long-time tail. Instead, return a\n",
" percentile.\n",
" function. See Janke 2002 [8] and Weigel, Janke 2010 [9].\n",
"\n",
" Due to the short simulation length, it is not possible to fit an\n",
" exponential to the long-time tail. Instead, return a percentile.\n",
" '''\n",
" summary = []\n",
" fig = plt.figure(figsize=(10, 6))\n",
Expand Down Expand Up @@ -269,6 +289,17 @@
" return np.array(summary)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 5.2.1 Distance-based macromolecular properties\n",
"\n",
"How do $R_h$, $R_g$, $R_F$ and the diffusion coefficient $D$ depend on the number of monomers?\n",
"You can refer to the Flory theory of polymers, and assume you are simulating a real polymer in a\n",
"good solvent, with Flory exponent $\\nu \\approx 0.588$."
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -391,7 +422,15 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate the diffusion coefficient of the polymers."
"#### 5.2.2 Diffusion coefficient using the MSD method\n",
"\n",
"Calculate the diffusion coefficient of the polymers using the mean-squared displacement.\n",
"\n",
"Recalling that for large $t$ the diffusion coefficient can be expressed as:\n",
"\n",
"$$6D = \\lim_{t\\to\\infty} \\frac{\\partial \\operatorname{MSD}(t)}{\\partial t}$$\n",
"\n",
"which is simply the slope of the MSD in the diffusive mode."
]
},
{
Expand All @@ -408,12 +447,12 @@
"plt.figure(figsize=(10, 10))\n",
"plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n",
"plt.ylabel(r'MSD [$\\sigma^2$]')\n",
"for index, (tau, msd) in enumerate(zip(tau_results, msd_results)):\n",
"for index, (tau, msd) in enumerate(zip(com_pos_tau_results, com_pos_msd_results)):\n",
" plt.loglog(tau[1:120], msd[1:120], label=r\"$N=${}\".format(N_MONOMERS[index]))\n",
"plt.loglog(2 * [tau[tau_f_index]], [0, np.max(msd_results)], '-', color='black')\n",
"plt.text(tau[tau_f_index], np.max(msd_results), r'$\\tau_{f}$')\n",
"plt.loglog(2 * [tau[tau_max_index]], [0, np.max(msd_results)], '-', color='black')\n",
"plt.text(tau[tau_max_index], np.max(msd_results), r'$\\tau_{max}$')\n",
"plt.loglog(2 * [tau[tau_f_index]], [0, np.max(com_pos_msd_results)], '-', color='black')\n",
"plt.text(tau[tau_f_index], np.max(com_pos_msd_results), r'$\\tau_{f}$')\n",
"plt.loglog(2 * [tau[tau_max_index]], [0, np.max(com_pos_msd_results)], '-', color='black')\n",
"plt.text(tau[tau_max_index], np.max(com_pos_msd_results), r'$\\tau_{max}$')\n",
"plt.legend()\n",
"plt.show()"
]
Expand All @@ -424,17 +463,17 @@
"metadata": {},
"outputs": [],
"source": [
"diffusion_results = np.zeros(len(N_MONOMERS))\n",
"diffusion_msd = np.zeros(len(N_MONOMERS))\n",
"plt.figure(figsize=(10, 8))\n",
"plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n",
"plt.ylabel(r'MSD [$\\sigma^2$]')\n",
"for index, (tau, msd) in enumerate(zip(tau_results, msd_results)):\n",
"for index, (tau, msd) in enumerate(zip(com_pos_tau_results, com_pos_msd_results)):\n",
" a, b = np.polyfit(tau[tau_f_index:tau_max_index], msd[tau_f_index:tau_max_index], 1)\n",
" x = np.array([tau[1], tau[tau_max_index - 1]])\n",
" p = plt.plot(x, a * x + b, '-')\n",
" plt.plot(tau[1:tau_max_index], msd[1:tau_max_index], 'o', color=p[0].get_color(),\n",
" label=r'$N=${}'.format(N_MONOMERS[index]))\n",
" diffusion_results[index] = a / 6\n",
" diffusion_msd[index] = a / 6\n",
"plt.legend()\n",
"plt.show()"
]
Expand All @@ -443,11 +482,16 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the dependence of the diffusion coefficient on the hydrodynamic radius.\n",
"\n",
"Recalling the formula for the diffusion coefficient of a short polymer in the Kirkwood–Zimm model:\n",
"\n",
"$$D = \\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta} \\left\\langle \\frac{1}{R_h} \\right\\rangle$$\n",
"\n",
"where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid."
"where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid.\n",
"\n",
"Hint: use $D = \\alpha_1 N^{-1} + \\alpha_2 N^{-\\beta}$ with `rh_exponent` for $\\beta$\n",
"and solve for $\\alpha_1, \\alpha_2$."
]
},
{
Expand All @@ -463,25 +507,148 @@
"\n",
"(a, b), _ = scipy.optimize.curve_fit(\n",
" lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n",
" N_MONOMERS, diffusion_results)\n",
" N_MONOMERS, diffusion_msd)\n",
"\n",
"label = f'''\\\n",
"$D^{{\\\\mathrm{{fit}}}} = \\\n",
" \\\\frac{{{a:.2f}}}{{N}} + \\\n",
" \\\\frac{{{b * 6 * np.pi:.2f} }}{{6\\\\pi}} \\\\cdot \\\n",
" \\\\frac{{{b * 6 * np.pi:.3f} }}{{6\\\\pi}} \\\\cdot \\\n",
" \\\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \\\n",
"'''\n",
"\n",
"fig = plt.figure(figsize=(10, 6))\n",
"x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n",
"plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n",
"plt.plot(N_MONOMERS, diffusion_msd, 'o', label='$D^{\\\\mathrm{simulation}}$')\n",
"plt.xlabel('Number of monomers $N$')\n",
"plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 5.2.3 Diffusion coefficient using the Green–Kubo method\n",
"\n",
"Plot the autocorrelation function and check that the decay is roughly exponential.\n",
"\n",
"Hint: use $D = \\alpha e^{-\\beta \\tau}$ and solve for $\\alpha, \\beta$. You can leave out\n",
"the first data point in the ACF if necessary, and limit the fit to the stable region\n",
"in the first 20 data points."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def exponential(x, a, b):\n",
" return a * np.exp(-b * x)\n",
"\n",
"fig = plt.figure(figsize=(10, 8))\n",
"for N, tau, acf in zip(N_MONOMERS, com_vel_tau_results, com_vel_acf_results):\n",
" popt, _ = scipy.optimize.curve_fit(exponential, tau[:20], acf[:20])\n",
" x = np.linspace(tau[0], tau[20-1], 100)\n",
" p = plt.plot(x, exponential(x, *popt), '-')\n",
" plt.plot(tau[:20], acf[:20], 'o',\n",
" color=p[0].get_color(), label=f'$R(\\\\tau)$ for $N = {N}$')\n",
"plt.xlabel('$\\\\tau$')\n",
"plt.ylabel('Autocorrelation function')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Green–Kubo integral for the diffusion coefficient take the following form:\n",
"\n",
"$$D = \\frac{1}{3} \\int_0^{+\\infty} \\left<\\vec{v_c}(\\tau)\\cdot\\vec{v_c}(0)\\right>\\, \\mathrm{d}\\tau$$\n",
"\n",
"Since our simulation is finite in time, we need to integrate up until $\\tau_{\\mathrm{int}}$. To find\n",
"the optimal value of $\\tau_{\\mathrm{int}}$, plot the integral as a function of $\\tau_{\\mathrm{int}}$\n",
"until you see a plateau. This plateau is usually followed by strong oscillations due to low\n",
"statistics in the long time tail of the autocorrelation function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"diffusion_gk = []\n",
"fig = plt.figure(figsize=(10, 6))\n",
"for N, tau, acf in zip(N_MONOMERS, com_vel_tau_results, com_vel_acf_results):\n",
" x = np.arange(2, 28)\n",
" y = [1 / 3 * np.trapz(acf[:j], tau[:j]) for j in x]\n",
" plt.plot(tau[x], y, label=f'$D(\\\\tau_{{\\\\mathrm{{int}}}})$ for $N = {N}$')\n",
" diffusion_gk.append(np.mean(y[10:]))\n",
"plt.xlabel('$\\\\tau_{\\\\mathrm{int}}$')\n",
"plt.ylabel('$\\\\frac{1}{3} \\\\int_{\\\\tau=0}^{\\\\tau_{\\\\mathrm{int}}} \\\\left<\\\\vec{v_c}(\\\\tau)\\\\cdot\\\\vec{v_c}(0)\\\\right>\\\\, \\\\mathrm{d}\\\\tau$')\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(1, 4))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the dependence of the diffusion coefficient on the hydrodynamic radius."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(a, b), _ = scipy.optimize.curve_fit(\n",
" lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n",
" N_MONOMERS, diffusion_gk)\n",
"\n",
"label = f'''\\\n",
"$D^{{\\\\mathrm{{fit}}}} = \\\n",
" \\\\frac{{{a:.2f}}}{{N}} + \\\n",
" \\\\frac{{{b * 6 * np.pi:.3f} }}{{6\\\\pi}} \\\\cdot \\\n",
" \\\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \\\n",
"'''\n",
"\n",
"fig = plt.figure(figsize=(10, 8))\n",
"x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n",
"plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n",
"plt.plot(N_MONOMERS, diffusion_results, 'o', label='$D^{\\\\mathrm{simulation}}$')\n",
"plt.plot(N_MONOMERS, diffusion_gk, 'o', label='$D^{\\\\mathrm{simulation}}$')\n",
"plt.xlabel('Number of monomers $N$')\n",
"plt.ylabel('Diffusion coefficient')\n",
"plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us compare the value of the diffusion coefficients calcualted with the MSD and Green–Kubo methods:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(f'N\\tMSD\\t\\tGK\\t\\tdifference')\n",
"for N, d_msd, d_gk in zip(N_MONOMERS, diffusion_msd, diffusion_gk):\n",
" print(f'{N}\\t{d_msd:.2e}\\t{d_gk:.2e}\\t{np.ceil(np.abs(d_msd-d_gk) * 100 / d_msd):.0f}%')"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -494,7 +661,9 @@
"[4] P. G. de Gennes. *Scaling Concepts in Polymer Physics*. Cornell University Press, Ithaca, NY, 1979. \n",
"[5] M. Doi. *Introduction to Polymer Physics.* Clarendon Press, Oxford, 1996. \n",
"[6] Michael Rubinstein and Ralph H. Colby. *Polymer Physics.* Oxford University Press, Oxford, UK, 2003. \n",
"[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002."
"[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002. \n",
"[8] W. Janke, Statistical analysis of simulations: Data correlations and error estimation, *Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms*, Lecture Notes, J. Grotendorst, D. Marx, A. Muramatsu (Eds.), John von Neumann Institute for Computing, 10:423–445, 2002. https://www.physik.uni-leipzig.de/~janke/Paper/nic10_423_2002.pdf \n",
"[9] M. Weigel, W. Janke, Error estimation and reduction with cross correlations, *Phys. Rev. E*, 81:066701, 2010, doi:[10.1103/PhysRevE.81.066701](https://doi.org/10.1103/PhysRevE.81.066701); Erratum-ibid 81:069902, 2010, doi:[10.1103/PhysRevE.81.069902](https://doi.org/10.1103/PhysRevE.81.069902). "
]
}
],
Expand Down
Loading

0 comments on commit 8841251

Please sign in to comment.