FEniCS: Read mesh from msh (XDMF) with subdomains

The overall procedure is generating Mesh() , MeshValueCollection objects first, and reading information into them. Then generating MeshFunctionSizet accordingly to define dx , ds ... for subdomains. Step 0. Open corresponding XDMF files; e.g.,  with XDMFFile(MPI.comm_world, "mesh.xdmf") as xdmf_infile: Step 1. Create Mesh object, using mesh = dolfin.cpp.mesh.Mesh() Step 2. Read Mesh from XDMF file, e.g., Step 3. Create MeshValueCollection object as container of subdomain information; mvc_subdomain = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim()) Step 4. Read subdomain information,, "name-to-read") The "name-to-read" here is the name given in XDMF file generated by meshio . Step 5. Using step 3, 4 to read boundaries (open another file if necessary): mvc_boundaries = dolfin.MeshValueCollection("size_t", mesh, mesh.topology().dim()-1) xdmf_infile.rea

Bending energy and persistence length

Persistence length $L_p$ is a basic mechanical property quantifying the bending stiffness of a polymer and is defined as the relaxation of $\langle\cos(\theta)\rangle$ of bond angles of the polymer chain: $$ \langle\cos(\theta(s))\rangle = \exp(-sl/L_p) $$ The problem of calculating $L_p$ becomes calculating $\langle\cos(\theta(s))\rangle$, if we choose the $\theta$ as the angle of adjacent bonds, i.e., $s=1$, we have: $$ \langle\cos(\theta)\rangle = \exp(-l/L_p) $$ in simulations, the stiffness constant is given by a bending energy $U_b$, when $U_b$ is large, where excluded volume effect is negligible, we can calculate $\langle\cos(\theta)\rangle$ as $$\langle\cos(\theta)\rangle=\frac{\int_0^\pi \cos(\theta)\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}{\int_0^\pi\sin(\theta)\exp(-\beta U_b)\mathrm{d}\theta}$$ the $\sin(\theta)$ is a geometric weight that when 2 bonds are at an angle of $\theta$, then in 3-D space, the number of bonds is proportion to $\sin(\theta)$. Here is an example

Steady compliance (Linear viscoelasty)

Boltzmann superposition: assume that there is no stress when $t\le 0$, so the integral starts from $0$ and $\sigma(0^{-})=0$: $$\begin{equation} \gamma(t)=\int_0^t J(t-t^\prime)\dot{\sigma}(t^\prime)\mathrm{d}t^\prime \end{equation}$$  Performing Laplace Transform yields:  $$\begin{align} \hat{\gamma}(s)&=\hat{J}(s)\hat{\dot{\sigma}}(s)\\ &=\hat{J}(s)\left(s\hat{\sigma}(s)-\sigma(0^{-})\right) \\ &=\hat{J}(s)\left(s\hat{G}(s)\hat{\dot{\gamma}}(s)-\sigma(0^{-})\right)\\ &=\hat{J}(s)\left(s\hat{G}(s)(s\hat{\gamma}(s)-\gamma(0^-))-\sigma(0^{-})\right)\\ &=s^2\hat{J}(s)\hat{\gamma}(s)\hat{G}(s)  \end{align}$$ Here we let $\hat{f}(s):=\mathcal{L}\lbrace f(t)\rbrace(s)$ be the Laplace transformation and since the stress starts at time $0$, $\gamma(0^-)$ and $\sigma(0^-)$ are simply $0$. The 2nd to 3rd step is derived from the convolution relation $\sigma(t)=\int_0^t G(t-t^\prime)\dot{\gamma}(t^\prime)\mathrm{d}t^\prime$. Cancelling out the $\hat{\gamma}(s)$ and rearrangin

Saddle point method and end-to-end distribution of polymer models

Introduction Saddle point methods are widely used in estimation of integrals with form $$ I = \int \exp\left(-f(x)\right)\mathrm{d}x $$ where function $f$ can be approximated by first 2 terms of its Taylor series around some $x_0$, i.e. $$ f(x)\approx f(x_0) + f^\prime (x_0)(x-x_0) + \frac{1}{2}f^{\prime\prime}(x_0)(x-x_0)^2 $$ The integral is thus approximated by its saddle point, where $f^\prime (x_0)=0$ and $f^{\prime\prime}(x_0)>0$: $$ \begin{align} I&\approx \int \exp\left(-f(x_0) - \frac{1}{2}f^{\prime\prime}(x_0)(x-x_0)^2\right) \mathrm{d}x\\ &=\exp(-f(x_0))\sqrt{\frac{2\pi}{f^{\prime\prime}(x_0)}} \end{align}$$ Examples Stirling's formula: With the knowledge of $\Gamma$ function we know that $$N!=\int_0^\infty \exp(-x)x^N\mathrm{d}x$$ let $f(x):=x-N\ln(x)$, with large $N$, the negative part is negligible, solving $f^\prime (x) = 0$, we have: $$N!\approx\exp(-N+N\ln(N))\sqrt{2\pi{}N}=\sqrt{2\pi{}N}\left(\frac{N}{e}\right)^N$$ Partition function:

Fancy way to calculate polymer chain model characteristics via NumPy

NumPy is the fundamental package for scientific computing with Python. It contains among other things: a powerful N-dimensional array object sophisticated (broadcasting) functions tools for integrating C/C++ and Fortran code useful linear algebra, Fourier transform, and random number capabilities Besides its obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide variety of databases. The ndarray is also the basic data structure of JIT tools like Numba, which boosts python dramatically. Recently, I wrote a program about calculation of the gyration tensor, center of mass and end to end vectors of polymer chain models under periodic boundary condition, and I found that using NumPy and Numba would finish the task in very few lines of codes, and also applicable to various situations. Here is the code: def pbc(r, d): return

Not SQ again...

Structure factor (SQ) is a characterization quantity which is frequently calculated in molecular simulations. It is easily to code the calculation program according to the definition. In most molecular simulations, periodic boundary condition is adopted therefore $$S(\mathbf{q}):=\mathcal{FT}\lbrace\langle\rho(\mathbf{r})\rho(\mathbf{0})\rangle\rbrace(\mathbf{q})$$ is actually a circular convolution, where FFT can give a dramatic boost in contrast to calculate directly. The steps using FFT are: Calculate $\rho(\mathbf{r})$, which is a summation of Dirac-Delta function that can be estimated as a hisotgram; Calculate $\hat{\rho}(\mathbf{q})$ by FFT; $S(\mathbf{q})=|\hat{\rho}(\mathbf{q})|^2$; Calculate mean over modulus: $S(q)=\int S(\mathbf{q})\delta(|\mathbf{q}|-q)\mathrm{d}\mathbf{q}/\int \delta(|\mathbf{q}|-q)\mathrm{d}\mathbf{q}$ Efficiency If the simulation box is divided into $N$ (in 3D systems, $N=N_xN_yN_z$ for example) bins, the FFT gives $O(N\log(N))$ complexity and step

Free energy calculation: umbrella integration

Formulae of umbrella sampling method can be found on wikipedia. In umbrella sampling, a reaction coordinate $\xi$ is pre-defined from the atomic coordinates, a bias potential is added to the atoms of interest to keep $\xi$ of the system at a specific window $\xi_w$. The bias form is usually a harmonic potential: $$u^b_w(\xi)=0.5k_w(\xi-\xi_w)^2$$ Therefore, the energy of biased system $A^b_w = A^{ub}_w + u^b_w(\xi)$. The superscript $ub$ is short for "un-biased". In the simulations, we can sample the reaction coordinate in each window and evaluate their distribution $P^b(\xi)$, since the free energy $A=-k_BT\ln(P)$, we have: $$A_w^{ub}(\xi) = -k_BT\ln(P^b_w(\xi))-u^b_w(\xi)-F_w$$ with $F_w$ is a reference free energy of each window and remains an unknown constant. One method to derive $F_w$ is WHAM, in year 2005, Kästner et al.  (The Journal of Chemical Physics 123, no. 14 (October 8, 2005): 144104.) have proposed a new method whose idea is to take derivative of $A^u_w$ w