Calculation of isotropic-nematic phase transition

Question

In this paper, the theory behind the isotropic-nematic phase transition is discussed. Furthermore, an algorithm is given to calculate some properties of this phase transition.

I have written a program in C which gives me the free energy, pressure, chemical potential and nematic order parameter for a given concentration $c$ of the particles in the system. Now I want to find the concentrations $c_I$ and $c_N$, which define the interval where the system phase separates, i.e., it is partly nematic and partly isotropic in order to lower its free energy.

As the article explains, one needs to solve the system of equations $$ p_{\text{iso}}(c_I) = p(c_N) \\ \mu_{\text{iso}}(c_I) = \mu(c_N), $$ for $c_{I,N}$. The expressions $p_{\text{iso}}$ and $\mu_{\text{iso}}$ are analytically known and I can calculate $p$ and $\mu$ for certain discrete values of of the concentration.

My question is: how to solve these two equations in C? I know that one can use Newton-Rhapson to solve zeros of a function numerically. However, as I only know the functions $p$ and $\mu$ for certain values of $c$, I do not really know how to tackle this problem.

Can someone help me find a way to do this? (Preferably an algorithm of some sorts.) Disclaimer: I am starting to learn C, but I am not extremely good at it.

Answer 1

I don't have enough reputation to write a comment, so this one goes down as an answer.

I take it you've written a code to solve the angular distribution $\psi$ numerically for a given concentration $c$, using the self-consistent procedure outlined in the paper. Now that you can write the free energy as a function of phi, and the pressure is the $c$-derivative of the free energy, you can, in principle, get the pressure for any value of $c$.

I am slightly confused by the fact that it seems, in the derivation for pressure, as if the $\psi$s do not depend on $c$. And they obviously do. So I believe there is an implicit approximation here. If you keep to the same approximation, you can write down the derivative of the pressure by hand. You have then a way to evaluate both the pressure and its derivative at any given $c$ (this is why I am not sure what you mean when you say that you only know them for certain values of $c$). Now that you have these functions (and similarly for the chemical potentials), writing down the Newton-Raphson scheme should pose no problem (you have a two variable, two function setting, so you'll have to differentiate both functions with respect to each variable, and then invert the resulting Jacobian matrix).

EDIT: So I decided to write the code myself, in "Python". I say "Python" instead of Python, for I used it as a quick drop in toolbox. I uses ipython with the -pylab switch so it'll load parts of scipy, numpy for computations.

Ntheta  = 128
ks = arange(Ntheta+1)[1:]

thetas = pi/2.*ks/(Ntheta+1)

deltas = zeros_like(thetas)
deltas[0]    = 1 - (cos(thetas[0]) + cos(thetas[1]))/2.
deltas[1:-1] = (cos(thetas[:-2]) - cos(thetas[2:]))/2.
deltas[-1]   = (cos(thetas[-1]) + cos(thetas[-2]))/2.

Nphi = 256
js = arange(Nphi+1)[1:]
phis = 2.*pi*js/(Nphi+1)
Ks = zeros((Ntheta, Ntheta))
for k in range(Ntheta):
  for l in range(Ntheta):
    g = sqrt(1. - (cos(thetas[k])*cos(thetas[l]) + sin(thetas[k])*sin(thetas[l])*cos(phis))**2)
    Ks[k, l] = 2.*pi/(Nphi+1) * (3./2*g[0] + sum(g[1:-1]) + 3./2*g[-1])

def thefunc(c):
  taus = (c/pi)**2 * exp(-2*c**2*thetas**2/pi)

  for i in range(128):
    As = 16.*c/pi * dot(Ks, deltas*taus)
    Z  = 4.*pi*sum(deltas*exp(-As))
    psis = 1./Z * exp(-As)
    taus = psis.copy()

  S = 4.*pi*sum(deltas*psis*(3.*cos(thetas)**2-1.)/2.)
  theint = dot(deltas*psis, dot(Ks, deltas*psis))
  p = c + 32.*c**2*theint
  pder = 1. + 64.*c*theint
  f  = log(c) - 1. + 4.*pi*sum(deltas*psis*log(psis)) + 32.*c*theint
  mu = f + p/c
  muder = pder/c
  return p, pder, mu, muder

def isofunc(c):
  p = c + c**2
  pder = 1. + 2.*c
  mu = log(c/(4.*pi)) + 2.*c
  muder = 1./c + 2.
  return p, pder, mu, muder

def totalfunc(ct, c):
   pt, pdert, mut, mudert = thefunc(ct)
   p, pder, mu, muder     = isofunc(c)
   F = -array([pt - p, mut - mu])
   J = array([[pdert, -pder], [mudert, -muder]])
   dc = solve(J, F)
   return ct+dc[0], c+dc[1]

ct = 4.
c  = 3.
for i in range(100):
  ct, c = totalfunc(ct, c)

Here thefunc computes the self consistent solution for $\psi$ returning the pressure, the chemical potential and their derivatives. isofunc does the same for the isotropic phase. totalfunc does a Newton-Raphson iteration. Sorry for the names of the functions. The variables should be self explanatory, though.