Merge pull request #169 from sandialabs/rescale-potential-stiffness

rescale potential stiffness

Author: mrbuche

docs/source/physics/single_chain/fjc/thermodynamics/modified_canonical/example_asymptotic.ipynb

@@ -27,17 +27,11 @@
   },
   {
    "cell_type": "markdown",
-   "id": "8abeab96",
-   "metadata": {},
-   "source": [
-    "## Strong potential"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "8c7dadae",
+   "id": "1592d36b",
    "metadata": {},
    "source": [
+    "## Strong potential\n",
+    "\n",
     "For sufficiently strong potentials, the modified canonical ensemble can be accurately approximated using the reference system (the isometric ensemble) and an asymptotic correction. For example, the nondimensional force $\\eta$ as a function of the nondimensional potential distance $\\gamma$ is approximated as\n",
     "\n",
     "$$\n",
@@ -68,11 +62,10 @@
    "source": [
     "gamma = np.linspace(1e-3, 99e-2, 256)\n",
     "for varpi in [1e0, 1e1, 1e2]:\n",
-    "    w = varpi*fjc.number_of_links**2\n",
-    "    eta = fjc.nondimensional_force(gamma, w)\n",
+    "    eta = fjc.nondimensional_force(gamma, varpi)\n",
     "    line = plt.plot(gamma, eta, label=r'$\\varpi=$' + str(varpi))\n",
     "    eta_asymptotic = \\\n",
-    "        fjc.asymptotic.strong_potential.nondimensional_force(gamma, w)\n",
+    "        fjc.asymptotic.strong_potential.nondimensional_force(gamma, varpi)\n",
     "    plt.plot(gamma, eta_asymptotic, ':', color=line[0].get_color())\n",
     "plt.legend()\n",
     "plt.xlim([0, 1])\n",
@@ -84,17 +77,11 @@
   },
   {
    "cell_type": "markdown",
-   "id": "1d21a64d",
-   "metadata": {},
-   "source": [
-    "## Weak potential"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "dce6f3d2",
+   "id": "8a77ab06",
    "metadata": {},
    "source": [
+    "## Weak potential\n",
+    "\n",
     "For sufficiently distant potentials, the modified canonical ensemble can be accurately approximated using the reference system (the isotensional ensemble) and an asymptotic correction. The potential is considered sufficiently distant when the length of center of the potential well to the end of the chain experiencing it is much larger than the expected end-to-end length of the chain. This disparity in length is only typically possible considering weak potentials. For example, if $\\eta/N_b\\varpi$ is the nondimensional potential distance, the nondimensional end-to-end length per link $\\gamma$ as a function of the effective nondimensional potential force $\\eta$ is approximated as\n",
     "\n",
     "$$\n",
@@ -105,14 +92,9 @@
     "\n",
     "$$\n",
     "  \\gamma(\\eta) = \\frac{1}{N_b}\\frac{\\partial}{\\partial\\eta}\\,\\ln\\left[\\iiint Q_0(\\boldsymbol{\\gamma}') \\, e^{N_b\\boldsymbol{\\eta}\\cdot\\boldsymbol{\\gamma}'} e^{-\\tfrac{\\varpi}{2} N_b^2\\left(\\boldsymbol{\\gamma}'\\right)^2} d^3\\boldsymbol{\\gamma}'\\right].\n",
-    "$$\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "8a77ab06",
-   "metadata": {},
-   "source": [
+    "$$\n",
+    "\n",
+    "\n",
     "This exact relation is plotted below along with the asymptotic relation while varying $\\varpi$, the nondimensional potential stiffness. As $\\varpi$ decreases and/or the nondimensional force $\\eta$ increases (the nondimensional potential distance $\\eta/N_b\\varpi$ increases), the asymptotic approach appears to do increasingly well. Notably, the asymptoic approach appears to succeed for sufficiently distance potentials for any value of $\\varpi$. This is because the freely-jointed chain model has inextensible links, so even a stiff potential that is distant ($\\eta/N_b\\varpi\\gg 1$) will not stretch the chain past $\\gamma=1$. For chain models which have extensible links, the link stiffness will compete with the potential stiffness, such that the potential would need to be weak ($\\varpi\\ll 1$) in addition to distant ($\\eta/N_b\\varpi\\gg 1$) in order for the asymptotic approach to become accurate."
    ]
   },
@@ -125,11 +107,10 @@
    "source": [
     "gamma_in = np.linspace(5e-1, 5, 256)\n",
     "for varpi in [1e0, 1e-1, 1e-2]:\n",
-    "    w = varpi*fjc.number_of_links**2\n",
-    "    gamma = fjc.nondimensional_end_to_end_length_per_link(gamma_in, w)\n",
+    "    gamma = fjc.nondimensional_end_to_end_length_per_link(gamma_in, varpi)\n",
     "    line = plt.plot(gamma, gamma_in, label=r'$\\varpi=$' + str(varpi))\n",
     "    gamma_asymptotic = \\\n",
-    "        fjc.asymptotic.weak_potential.nondimensional_end_to_end_length_per_link(gamma_in, w)\n",
+    "        fjc.asymptotic.weak_potential.nondimensional_end_to_end_length_per_link(gamma_in, varpi)\n",
     "    plt.plot(gamma_asymptotic, gamma_in, ':', color=line[0].get_color())\n",
     "plt.legend()\n",
     "plt.xlim([0, 1])\n",

src/physics/single_chain/fjc/thermodynamics/modified_canonical/asymptotic/strong_potential/mod.rs

@@ -30,8 +30,7 @@ pub struct FJC
 /// The expected force as a function of the applied potential distance, potential stiffness, and temperature, parameterized by the number of links and link length.
 pub fn force(number_of_links: &u8, link_length: &f64, potential_distance: &f64, potential_stiffness: &f64, temperature: &f64) -> f64
 {
-    let contour_length = (*number_of_links as f64)*link_length;
-    BOLTZMANN_CONSTANT*temperature/link_length*nondimensional_force(number_of_links, &(*potential_distance/contour_length), &(potential_stiffness*contour_length.powi(2)/BOLTZMANN_CONSTANT/temperature))
+    BOLTZMANN_CONSTANT*temperature/link_length*nondimensional_force(number_of_links, &(*potential_distance/((*number_of_links as f64)*link_length)), &(potential_stiffness*link_length.powi(2)/BOLTZMANN_CONSTANT/temperature))
 }
 
 /// The expected nondimensional force as a function of the applied nondimensional potential distance and nondimensional potential stiffness, parameterized by the number of links.
@@ -46,21 +45,19 @@ pub fn nondimensional_force(number_of_links: &u8, nondimensional_potential_dista
     let sum_1: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p - 1)).sum();
     let sum_2: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p - 2)).sum();
     let sum_3: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p - 3)).sum();
-    (1.0/nondimensional_potential_distance + (0.5*number_of_links_f64 - 1.0)*sum_1/sum_0)/number_of_links_f64 + 0.5/nondimensional_potential_stiffness/number_of_links_f64*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0*((number_of_links_f64 - 2.0)*(sum_1/sum_0).powi(2) - (number_of_links_f64 - 3.0)*sum_2/sum_0) - (0.5*number_of_links_f64 - 1.5)*((0.5*number_of_links_f64 - 1.0)*sum_1*sum_2/sum_0.powi(2) - (0.5*number_of_links_f64 - 2.0)*sum_3/sum_0)) + 2.0*nondimensional_potential_distance.powi(-3) - 2.0*((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0 + nondimensional_potential_distance.powi(-1))*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*(sum_1/sum_0).powi(2) - (0.5*number_of_links_f64 - 1.5)*sum_2/sum_0) - nondimensional_potential_distance.powi(-2)))
+    (1.0/nondimensional_potential_distance + (0.5*number_of_links_f64 - 1.0)*sum_1/sum_0)/number_of_links_f64 + 0.5/nondimensional_potential_stiffness/number_of_links_f64.powi(3)*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0*((number_of_links_f64 - 2.0)*(sum_1/sum_0).powi(2) - (number_of_links_f64 - 3.0)*sum_2/sum_0) - (0.5*number_of_links_f64 - 1.5)*((0.5*number_of_links_f64 - 1.0)*sum_1*sum_2/sum_0.powi(2) - (0.5*number_of_links_f64 - 2.0)*sum_3/sum_0)) + 2.0*nondimensional_potential_distance.powi(-3) - 2.0*((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0 + nondimensional_potential_distance.powi(-1))*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*(sum_1/sum_0).powi(2) - (0.5*number_of_links_f64 - 1.5)*sum_2/sum_0) - nondimensional_potential_distance.powi(-2)))
 }
 
 /// The Helmholtz free energy as a function of the applied potential distance, potential stiffness, and temperature, parameterized by the number of links, link length, and hinge mass.
 pub fn helmholtz_free_energy(number_of_links: &u8, link_length: &f64, hinge_mass: &f64, potential_distance: &f64, potential_stiffness: &f64, temperature: &f64) -> f64
 {
-    let contour_length = (*number_of_links as f64)*link_length;
-    nondimensional_helmholtz_free_energy(number_of_links, link_length, hinge_mass, &(potential_distance/contour_length), &(potential_stiffness*contour_length.powi(2)/BOLTZMANN_CONSTANT/temperature), temperature)*BOLTZMANN_CONSTANT*temperature
+    nondimensional_helmholtz_free_energy(number_of_links, link_length, hinge_mass, &(potential_distance/((*number_of_links as f64)*link_length)), &(potential_stiffness*link_length.powi(2)/BOLTZMANN_CONSTANT/temperature), temperature)*BOLTZMANN_CONSTANT*temperature
 }
 
 /// The Helmholtz free energy per link as a function of the applied potential distance, potential stiffness, and temperature, parameterized by the number of links, link length, and hinge mass.
 pub fn helmholtz_free_energy_per_link(number_of_links: &u8, link_length: &f64, hinge_mass: &f64, potential_distance: &f64, potential_stiffness: &f64, temperature: &f64) -> f64
 {
-    let contour_length = (*number_of_links as f64)*link_length;
-    nondimensional_helmholtz_free_energy_per_link(number_of_links, link_length, hinge_mass, &(potential_distance/contour_length), &(potential_stiffness*contour_length.powi(2)/BOLTZMANN_CONSTANT/temperature), temperature)*BOLTZMANN_CONSTANT*temperature
+    nondimensional_helmholtz_free_energy_per_link(number_of_links, link_length, hinge_mass, &(potential_distance/((*number_of_links as f64)*link_length)), &(potential_stiffness*link_length.powi(2)/BOLTZMANN_CONSTANT/temperature), temperature)*BOLTZMANN_CONSTANT*temperature
 }
 
 /// The relative Helmholtz free energy as a function of the applied potential distance, potential stiffness, and temperature, parameterized by the number of links and link length.
@@ -74,10 +71,12 @@ pub fn relative_helmholtz_free_energy_per_link(number_of_links: &u8, link_length
 {
     helmholtz_free_energy_per_link(number_of_links, link_length, &1.0, potential_distance, potential_stiffness, temperature) - helmholtz_free_energy_per_link(number_of_links, link_length, &1.0, &(ZERO*(*number_of_links as f64)*link_length), potential_stiffness, temperature)
 }
+
 /// The nondimensional Helmholtz free energy as a function of the applied nondimensional potential distance, nondimensional potential stiffness, and temperature, parameterized by the number of links, link length, and hinge mass.
 pub fn nondimensional_helmholtz_free_energy(number_of_links: &u8, link_length: &f64, hinge_mass: &f64, nondimensional_potential_distance: &f64, nondimensional_potential_stiffness: &f64, temperature: &f64) -> f64
 {
     let number_of_links_f64 = *number_of_links as f64;
+    let number_of_links_squared = number_of_links_f64.powi(2);
     let contour_length = number_of_links_f64*link_length;
     let n = *number_of_links as u128;
     let p: i32 = (number_of_links - 2).into();
@@ -86,7 +85,7 @@ pub fn nondimensional_helmholtz_free_energy(number_of_links: &u8, link_length: &
     let sum_0: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p)).sum();
     let sum_1: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p - 1)).sum();
     let sum_2: f64 = (0..=k-1).collect::<Vec::<u128>>().iter().map(|s| (-1.0_f64).powf(*s as f64)*(((1..=n).product::<u128>()/(1..=*s).product::<u128>()/(1..=n-s).product::<u128>()) as f64)*(m - (*s as f64)/number_of_links_f64).powi(p - 2)).sum();
-    -(0.125/PI/nondimensional_potential_distance*(n.pow(n as u32) as f64)/((1..=n-2).product::<u128>() as f64)*sum_0/contour_length.powi(3)).ln() - (number_of_links_f64 - 1.0)*(8.0*PI.powi(2)*hinge_mass*link_length.powi(2)*BOLTZMANN_CONSTANT*temperature/PLANCK_CONSTANT.powi(2)).ln() - 1.5*(2.0*PI/nondimensional_potential_stiffness).ln() - 3.0*(contour_length).ln() + 0.5/nondimensional_potential_stiffness*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*(sum_1/sum_0).powi(2) - (0.5*number_of_links_f64 - 1.5)*sum_2/sum_0) - nondimensional_potential_distance.powi(-2) - ((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0 + nondimensional_potential_distance.powi(-1)).powi(2))
+    -(0.125/PI/nondimensional_potential_distance*(n.pow(n as u32) as f64)/((1..=n-2).product::<u128>() as f64)*sum_0/contour_length.powi(3)).ln() - (number_of_links_f64 - 1.0)*(8.0*PI.powi(2)*hinge_mass*link_length.powi(2)*BOLTZMANN_CONSTANT*temperature/PLANCK_CONSTANT.powi(2)).ln() - 1.5*(2.0*PI/nondimensional_potential_stiffness/number_of_links_squared).ln() - 3.0*(contour_length).ln() + 0.5/nondimensional_potential_stiffness/number_of_links_squared*((0.5*number_of_links_f64 - 1.0)*((0.5*number_of_links_f64 - 1.0)*(sum_1/sum_0).powi(2) - (0.5*number_of_links_f64 - 1.5)*sum_2/sum_0) - nondimensional_potential_distance.powi(-2) - ((0.5*number_of_links_f64 - 1.0)*sum_1/sum_0 + nondimensional_potential_distance.powi(-1)).powi(2))
 }
 
 /// The nondimensional Helmholtz free energy per link as a function of the applied nondimensional potential distance, nondimensional potential stiffness, and temperature, parameterized by the number of links, link length, and hinge mass.

src/physics/single_chain/fjc/thermodynamics/modified_canonical/asymptotic/strong_potential/test.jl

@@ -93,7 +93,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         force = model.force(potential_distance, potential_stiffness, temperature)
@@ -128,7 +128,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         helmholtz_free_energy = model.helmholtz_free_energy(
@@ -169,7 +169,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         helmholtz_free_energy_per_link = model.helmholtz_free_energy_per_link(
@@ -209,7 +209,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         relative_helmholtz_free_energy = model.relative_helmholtz_free_energy(
@@ -249,7 +249,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         relative_helmholtz_free_energy_per_link =
@@ -285,7 +285,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         helmholtz_free_energy = model.helmholtz_free_energy(
@@ -325,7 +325,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         relative_helmholtz_free_energy = model.relative_helmholtz_free_energy(
@@ -438,7 +438,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         helmholtz_free_energy = model.helmholtz_free_energy(
@@ -482,7 +482,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         helmholtz_free_energy_per_link = model.helmholtz_free_energy_per_link(
@@ -605,7 +605,7 @@ end
         temperature =
             parameters.temperature_reference + parameters.temperature_scale * (0.5 - rand())
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         relative_helmholtz_free_energy_0 = model.relative_helmholtz_free_energy(
@@ -633,7 +633,7 @@ end
         temperature =
             parameters.temperature_reference + parameters.temperature_scale * (0.5 - rand())
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         relative_helmholtz_free_energy_per_link_0 =
@@ -707,7 +707,7 @@ end
         potential_distance =
             nondimensional_potential_distance * number_of_links * link_length
         potential_stiffness =
-            nondimensional_potential_stiffness / (number_of_links * link_length)^2 *
+            nondimensional_potential_stiffness / link_length^2 *
             BOLTZMANN_CONSTANT *
             temperature
         force = model.force(potential_distance, potential_stiffness, temperature)