adtzlr / matadi

Material Definition with Automatic Differentiation
GNU General Public License v3.0
21 stars 2 forks source link

Bug in Calculation of free energy of single chain #37

Closed adtzlr closed 2 years ago

adtzlr commented 2 years ago

The driving force of the 1d micro-sphere model has to be integrated w.r.t. the stretch, in order to obtain the 1d strain energy function. The Padé approximation of the inverse Langevin function will be integrated w.r.t. the stretch. Here is the result of WolframAlpha

integrate (3*N - λ^2)/(N-λ^2) dλ

which gives:

λ + 2 sqrt(N) tanh^(-1)(λ/sqrt(N))

This is transformed into a Python function:

def langevin(stretch, mu, N):
    """Langevin model (Padé approximation) given by the free energy 
    of a single chain as a function of the stretch."""

    return mu * (stretch + 2 * sqrt(N) * atanh(stretch / sqrt(N)))

Originally posted by @adtzlr in https://github.com/adtzlr/matadi/issues/28#issuecomment-963997189

adtzlr commented 2 years ago

It has to be multiplied by λ

integrate (3*N - λ^2)/(N-λ^2) * λ dλ

which gives

λ^2/2 - N log(-N + λ^2)

Again, as Python function:

def langevin(stretch, mu, N):
    """Langevin model (Padé approximation) given by the free energy 
    of a single chain as a function of the stretch."""

    return mu * (stretch ** 2 / 2 - N * log(stretch ** 2 - N))

Note that this strain energy function assumes a complex logarithm and only its derivative gives a real result for N > λ^2!