Multidimensional integral to compute the end-to-end distribution of a FENE ideal chain

How do I evaluate the mean end-to-end squared distance of a FENE ideal chain at fixed inverse temperature $\beta$ in the canonical ensemble?

This quantity is defined as the mean value $$\left<\left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)\right>$$

Under the canonical distribution $$\exp\left(-\beta \sum_{i=1}^NH(\vec r_i)\right)\text{ where } H(\vec r_i) = - \ln\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)\,.$$

Therefore, it would be $$\frac{ \int \prod_{i=1}^{N-1}d\vec r_i\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)^\beta \left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)} {\int \prod_{i=1}^{N-1}d\vec r_i\left(1 - \left(\frac{|\vec r_i| - r_0}{\Delta}\right)^2\right)^\beta}$$ where the integral is to be performed over all the values of $\vec r_i$ that give a positive argument of the logarithm. Edit: I realised that the integral that I wanted to solve was not the one written below, but the one written here above. The following is the question in its first (wrong) version, which is the one to which Yarchik responded.

How do I evaluate the quantity $\left.\frac{d}{d\alpha} \ln Z_\alpha\right|_{\alpha=0}$ for N=30 in Mathematica, when

$$Z_\alpha = V\int e^{-\beta [H(\vec r_1)+H(\vec r_2-\vec r_1)+...+H(\vec r_{N-1}-\vec r_{N-2})]+\alpha|r_{N-1}|^2}d\vec r_1d\vec r_2\cdots d\vec r_{N-1}$$ and $$H(\vec r) =- \ln\left(1-16(|\vec r|-1)^2\right) ?$$ The integration is extended to all the domain in which the integral is defined (that is, the argument of the function $H$ must be a vector with length in [3/4,5/4]).

I know how to declare the function H, but how can I nest 30 integrate functions without writing them explicitly? I don't want to rewrite my code from scratch if I change the value of N I wish to consider. So far I have the following (H is expressed in spherical coordinates, so that r is just the distance between the particles):

a = 1; d = 0.25; r0 = 1
H[r, r0, d] = -a*Log[1 - (r - r0)^2/d^2]
U = FullSimplify[-D[Log[Z[r0, d, b]], b]]


but I do not know how I can write Z.

For the record, I need this to compute the end to end distance distribution of an ideal chain of 30 atoms linked by FENE springs. $H$ is the hamiltonian associated with each of the 29 springs, and Z is the partition function of my system. The derivative I want to compute is a way of computing the mean of the parenthesis that multiplies $\alpha$.

Edit: thanks to Yarchik for noting that the domain of H was wrong. It should be right now, but in case there's an error I want to integrate over all the vectors for which H is defined.

Edit: I also realised that I was writing the end to end distance in a wrong way. As Yarchik pointed out, the way I was writing it I was computing the position of the center of mass, while now I'm actually computing the average end to end squared distance.

• Rather than trying to solve 30 nested integrals at once, perhaps you can do the first one, gain a little experience with that, and then try the second. Trying to solve the whole thing at once seems a rather poor choice. – bill s Nov 28 '14 at 0:19
• What do you mean by $\alpha\left(\sum_{i=1}^{N-1}\vec{r}_i\right)^2$? Do you need a norm in there somewhere? – Kellen Myers Nov 28 '14 at 4:53
• @KellenMyers That is meant to be a norm. The sum of vectors yields a vector, of which then I take the square norm (e.g., the scalar product of the vector by itself), denoted by the $^2$. Is it clearer now? – Ferdinando Randisi Nov 28 '14 at 9:50
• @bills The problem is that, once I do have a solution for a trivial $N$, I'd love it to be some kind of loop/iteration of the Integrate command that allows me to change $N$ quickly. My question would therefore be quickly solved if there were such a language functionality. – Ferdinando Randisi Nov 28 '14 at 9:51
• If you don't have a solution for trivial $N$, then there is no way you'll ever get a loop working. There are lots of different iteration functionalities... – bill s Nov 28 '14 at 14:19

I believe it is not possible to solve this problem with Mathematica in the presented form. Partly because there inconsistencies in the formulation:

Plot[-Log[1 - 16 (Abs[x] - 1)^2], {x, 1 - 1/4, 1 + 1/4}]


We see that the potential is, in fact, defined on a larger domain [3/4,5/4]

An approximate solution can be numerically obtained using a Monte Carlo integration. Notice, the expectation value of the center-of-mass position has to be computed directly, rather then through the partition function. That is to answer the question on how to formulate it for Mathematica in $N$-invariant way. Below I would like to consider a special analytically solvable case as is also suggested by others. We take $N=2$, i.e. a chain consisting of a single particle at infinite temperature, i.e. $\beta=0$. In this case we deal with a 3d integral that in view of the spherical symmetry can be reduced to 1d.

  Z[a_] := 4 \[Pi] Integrate[Exp[a r^2] r^2, {r, 1-1/4, 1 + 1/4}]


 Limit[D[Log[Z[a]], a], a -> 0]
(*   4323/3920  *)


Analogically, the case of finite temperatures can be considered and analytically integrated for integer values of $\beta$:

 r = Table[{b,
Limit[
D[Log[
Integrate[
Exp[a r^2] r^2 (1 - 16 (r - 1)^2)^b, {r, 1 - 1/4, 1 + 1/4}]], a],
a-> 0]}, {b, 0, 10}]
fr = Interpolation[r]
Plot[fr[b], {b, 0, 10}]


Finally, consider another extreme $N\rightarrow\infty$. Just for physical clarity it makes sense to redefine alpha such that the whole Hamiltonian has a form of interacting points with an extra potential acting on the center-of-mass of the chain. Notice, that in the case each particle is subject to this (quadratic) potential individually the problem is solvable.

Edit

Concerning the modified question on how to determine the distance to the end point $d$. An exact solution is probably not feasible, but a very good approximation can be suggested. First, modify the inter-particle potential $V(x)=-\ln(1-(x-1)^2/b^2)$. In the present setup $b=1/4$. What happens when $b\rightarrow 0$? The system is essentially a chain with fixed link-length ($=1$). In this case the end-to end distance is exactly $\sqrt{N}$. This follows from the fact that the average value $\langle \cos \theta\rangle=0$, where $\cos\theta$ is the angle between two neighbouring links. In the case when $1>b>0$ a very accurate result can be obtained using calculations in hands. Second figure shows the average link-length squared $\langle r^2\rangle$ as a function of the inverse temperature $\beta$. The end-to-end distance is therefore:

$d(\beta)\approx\sqrt{N\langle r^2\rangle(\beta)}.$

When temperature is zero, i.e., $\beta=\infty$ another exact result reads $d(\infty)=\sqrt{N}.$

• Thanks! I know of all these approximations, since they both have the effect that the distance between beads is fixed, but I want an exact result for my model. I did manage to reduce the problem to a much simpler one, but mathematica cannot compute the relevant integral even though it looks simple. I'm not sure of why. I'm starting another question for that integral. – Ferdinando Randisi Dec 4 '14 at 22:22

First, let's rewrite the mean value as the sum of the mean values of the scalar products of the single particle relative positions:

\begin{align*}r^2_{e2e} & = \left<\left(\sum_{i=1}^{N-1} \vec r_i\right)\cdot\left(\sum_{j=1}^{N-1} \vec r_j\right)\right>\\ &=\sum_{i=1}^{N-1} \sum_{j=1}^{N-1} \left< \vec r_i\cdot \vec r_j\right>\\ & = (N-1) \left<|\vec r_i|^2\right>+\sum_{i<j}\left<\vec r_i\cdot \vec r_j\right>\\ & = (N-1) \left<|\vec r_i|^2\right>+\sum_{i<j}\left<\vec r_i\right>\cdot \left<\vec r_j\right>\\\end{align*} Where the last passage can be done since the different positions are uncoupled, i.e. evolve indipendently. The interesting thing is that both the sincgle vector means in the last passage are both zeroes, since each component yelds an odd function integrated on a symmetric domain. Therefore, we can focus on computing the first average in the last passage. Averaging over all the degrees of freedom different from $i$ yields a factor of one, therefore now the problem is simply to find the average distance between two particles connected by the FENE potential (that is, by using a single coordinate).

As is customary in statiscal mechanics, let us define the one oscillator partition sum (that is,the normalisation constant of the system with only two particles), and solve it in polar coordinates. Since $H$ is radial, we get: $$Z = \int d\vec r_1\left(1 - \left(\frac{|\vec r_1| - r_0}{\Delta}\right)^2\right)^\beta = 4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^2\,.$$

Similarly we could write the numerator of the quantity that we wish to calculate, finally yielding:

$$r^2_{e2e} = (N-1)\frac{4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^4}{4\pi\int^{r_0+\Delta}_{r_0-\Delta} dr \left(1 - \left(\frac{r - r_0}{\Delta}\right)^2\right)^\beta r^2}$$

This can be computed with the Mathematica function Nintegrate (for some reason, using Integrate instead will give problems for some values of $\beta$, like $\beta = 10$.

(* by this point, b, N, r0 and D must have already been defined*)
I[a_,b_]:=NIntegrate[(1 - (r - r0)^2/D^2)^b r^a, {r,r0-D,r0+D}]
( N - 1 ) I[4,b]/I[2,b]