tag:blogger.com,1999:blog-28890822945976253392021-08-08T23:02:23.629-07:00This Blog is Untitled呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.comBlogger26125tag:blogger.com,1999:blog-2889082294597625339.post-80676486216169634412021-04-07T19:46:00.021-07:002021-04-08T23:21:37.048-07:00 FEniCS: Read mesh from msh (XDMF) with subdomainsThe 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,呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-10881328332938167202020-11-06T07:04:00.004-08:002021-01-27T07:35:10.090-08:00Bending energy and persistence lengthPersistence 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, 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-32343016408884965102020-10-20T20:26:00.002-07:002020-10-20T20:27:06.040-07:00Steady 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) \\
&=呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-40377780693160961182020-07-06T09:04:00.008-07:002020-09-06T19:56:02.997-07:00Saddle point method and end-to-end distribution of polymer modelsIntroduction
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^\呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com2tag:blogger.com,1999:blog-2889082294597625339.post-38565550898600638682020-05-23T11:15:00.001-07:002020-06-15T08:19:15.985-07:00Fancy way to calculate polymer chain model characteristics via NumPyNumPy 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 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-46648175800123009222020-01-18T09:30:00.003-08:002020-09-08T08:38:17.813-07:00Not 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 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-53167191723153570432020-01-11T05:34:00.003-08:002020-03-03T01:18:06.905-08:00Free energy calculation: umbrella integrationFormulae 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^呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-45238998033019827392019-12-11T04:56:00.001-08:002020-01-05T03:27:26.609-08:00Di-block copolymer analysis: GPCGel Permeation Chromatography is a well known method of measuring molecular weight distribution of polymers. For di-block copolymers, a proper calibration curve can be constructed by using Mark-Houwink parameters of homopolymers (see here). Usually, di-block copolymers are obtained by synthesis one block first then "grow" the other block from the end of preceding block. The molecular weight 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com1tag:blogger.com,1999:blog-2889082294597625339.post-74802870313984519312019-10-28T10:11:00.003-07:002019-10-28T10:11:41.314-07:00A note on Cython ParallelismInstead of using thread-local ndarrays, a very convenient method is to create an array shaped in (num_of_threads, <whatever the structure>) to make the changes thread-safe. For example:
cdef np.ndarray[np.double_t, ndim=2] ret
num_threads = multiprocessing.cpu_count()
ret = np.zeros((num_threads, nbins), dtype=np.float)and in the prange loop:
thread_num = openmp.omp_get_thread_num()
ret[呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-56785788294329691682019-10-18T00:40:00.001-07:002019-10-19T01:14:55.368-07:00Notes on scipy.optimize.minimizeProblem: evenly distributed points on ellipsoid surfaces.
Solution: numerically minimizing energy of charges constrained on the ellipsoid surface.
Assuming the ellipsoid satisfies equation:
$$\left(\frac{x}{a}\right)^2+\left(\frac{y}{b}\right)^2+\left(\frac{z}{c}\right)^2=1$$
the ratios of 3 axes of the ellipsoid surface are $(a,b,c)^T$. The minimization process is:
let
$$\mathbf{x}^\prime:=\frac呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-88296602880235361922019-10-01T11:31:00.000-07:002019-10-03T08:23:32.992-07:00Dealing with aggregates: I don't like tedious workI used to study nanoparticles (NPs) self-assembly in polymer matrices before 2017 by coarse-grained simulation method. Analyzing morphologies of self-assemblies is a tedious work: one needs to repeatedly watch frames from simulation trajectories, finding appropriate view points, counting number of aggregates, measure size of aggregates, etc. Obtaining integrate aggregates under periodic boundary 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-79875905166999091942019-10-01T08:42:00.002-07:002020-09-17T07:36:16.750-07:00Eigenvalues of circulant matricesA circulant matrix is defined as
$$C=\begin{bmatrix}c_{0}&c_{n-1}&\dots &c_{2}&c_{1}\\c_{1}&c_{0}&c_{n-1}&&c_{2}\\\vdots &c_{1}&c_{0}&\ddots &\vdots \\c_{n-2}&&\ddots &\ddots &c_{n-1}\\c_{n-1}&c_{n-2}&\dots &c_{1}&c_{0}\\\end{bmatrix}$$
where $C_{j, k}=c_{j-k \mod n}$, the $k$-th eigenvalue $\lambda_k$ and eigenvector $x_k$ satisfy $C\cdot x_k=\lambda_k x_k$, or, equivalently, $n$ equations as:
$呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-64179144642758473882019-09-22T07:40:00.001-07:002019-10-03T08:26:29.508-07:00Toeplitz, Circulant matrix and discrete convolution theoremA topelitz matrix is a matrix with form of
$$A=\begin{bmatrix} a & b & c & d & e \\ f & a & b & c & d \\ g & f & a & b & c \\ h & g & f & a & b \\ i & h & g & f & a\end{bmatrix}$$
Generally, the elements are denoted as $a_{i-j}$: $A_{ij}=A_{i+1, j+1}=a_{i-j}$. A circulant matrix is a special kind of Toeplitz matrix where each row vector is rotated one element:
$$C=\begin{bmatrix}c_{0}&c_{n-1}&\呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-37351219437135798092019-07-24T08:16:00.004-07:002019-10-03T08:24:14.708-07:00Distribution of segments on Gaussian chainThe probability distribution function of $P_i(\mathbf{r}_i-\mathbf{r}_\mathrm{cm})$ of ideal chain is evaluated, for $P_i$ represents the probability distribution fucntion of $i$th segment with respect to its centre of mass on the ideal chain. An ideal chain is modeled as multidimensional random walk, where the steps are independent, and the distribution of a step at mean length $b$ is given by $呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-5212322507823494562019-07-22T09:39:00.001-07:002019-10-28T09:44:02.513-07:00Simple Examples of Parallel Computing of CythonIt is more convenient to calculate some properties during a molecular simulation process by accessing data from API of the molecular simulation program than calculating after the whole simulation progress and dumping all the coordinates. Especially, for sharply fluctuated properties such as Virial, RDF, etc., require massive frames to calculate. It is expensive to dump densely. For some 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-64493045956060412212019-07-13T17:47:00.001-07:002019-08-05T10:10:49.506-07:00Anisotropy of ideal chainThe Gaussian chain is isotropic when averaged over conformation space and all orientations. Therefore, a Gaussian chain is dealt as sphere with radius of $R_g$, it’s radius of gyration. In the 2nd chapter of Rubinstein’s Polymer Physics, an exercise shows that the $R_g^2$ is asymmetric if one set the coordinate frame on its end-to-end vector, $\mathbf{R}_{ee}$. i.e., the $\mathbf{R}_{ee}$ vector 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-82186876887772075952019-07-12T11:13:00.000-07:002019-07-12T12:19:10.324-07:00Flory characteristic ratio of hindered rotation chainThis is an exercise on Rubinstein's textbook of Polymer Physics.Definition of Flory characteristic ratio:
$$C_\infty:=\lim_{n\to\infty}\frac{\langle \mathbf{R}^2 \rangle}{nb^2}$$
where $\langle \mathbf{R}^2 \rangle$ is the mean square end-to-end vector, $b$ is Kuhn length of the polymer and $n$ is corresponding chain length.
Calculation of $\langle \mathbf{R}^2\rangle$:
Let $\mathbf{r}_i$ be the 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-11571204557423691202019-07-01T22:19:00.002-07:002019-07-24T08:02:35.571-07:00RTFM, pls RTFM!!!I was trying to write a parallel version C-extension with Cython several days ago, to calculate Coulomb energy of $n$ points. I wrote the code as follows:
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
def u(double[:,:] x):
cdef long i
cdef long j,k
cdef double ret=0,tmp
with nogil, parallel():
for i in prange(x.shape[0]-1, schedule='static'):呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-7640248193839159772019-06-16T23:51:00.000-07:002019-06-17T00:30:38.886-07:00Find intersections numerically -- Using NumPyIntersects of 2 functionsAssuming function $f$ and $g$, scipy.optimize.fsolve would give an elegant solution:from scipy.optimize import fsolve
def f(xy):
x, y = xy
return (y - x ** 2, y - np.sin(x))
fsolve(f, [1.0, 1.0])
gives array([0.87672622, 0.76864886]) is the intersect of $y=x^2$ and $y=\sin(x)$.Intersects of 2 arraysThere is a function called numpy.intersect1d that returns the 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-71508677482179469042019-04-22T10:31:00.001-07:002019-04-22T12:08:39.756-07:00CAPVT II -- notes on Wheelock's LatinFirst declension nouns and adjectivesEndings
Singular Plural
Nom. -a -ae
Gen. -ae -ārum
Dat. -ae -īs
Acc. -am -ās
Abl. -ā -īs
Voc. -a -ae
Most 1st declension nouns are feminine, a few nouns denoting individuals engaged in what were among the Romans traditionall male occupations are masculine, e.g., poeta, poet; nauta, sailor; agricola, farmer; auriga, charioteer; incola, 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-41974981821947654652018-04-25T22:50:00.001-07:002019-04-22T12:08:39.694-07:00Veritas in lucem emergit.In contrast to put endings on verbs, aka conjugation, declension is to put endings to nouns and adjectives. The 3rd declension of nouns are usually learnt after the 1st and 2nd declension.
Singular Plural
Nominative / -es
Genitive -is -um
Dative -i -ibus
Accusative -em -es
Ablative -e -ibus
Nominative is derived from Latin word nomen, which means name, this is the form we name 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-61025234093868925472018-04-02T09:06:00.000-07:002019-04-22T12:08:39.730-07:00Esse an non esse?In this lesson, irregular verb sum and possum are learnt.
sum, esse, fui, futurum: be
I am, to be, have been, going to be
And pot- (can, be able).
possum, posse, potui: be able, can
pot- add some be yield can/be able to, the rules are simple: pot- or pos- if the be verb stats with "s".
Present tense active indicative mood: sum, es, est, sumus, estis, sunt
Present subjunctive 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-87804018995171200882018-04-02T07:54:00.000-07:002019-04-22T12:08:39.780-07:00Miser Catulle, desinas ineptire!In this lesson, subjunctive mood was induced. Changing a verb into subjunctive mood from indicative mood is to add vowel "a" and personal endings.
Indicative is to tell the truth, to point out something. Subjunctive mood is to exhort or tendence. Imperative mood is to command someone to do something (always command the 2nd person).
Using subjunctive mood for exhorting: ut/ne, so that/so that 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-56082429462336378812018-04-02T06:12:00.000-07:002019-04-22T12:08:39.718-07:00Veni, vidi, viciLatin verbs can exhibit their number, person, tense, voice and mood. In this lesson, only present tense indicative mood is considered.
Personal endings of words
-o/-m first person singular
-mus first person plural
-s second person singular
-tis second person plural
-t third person singular
-nt third person plural
The present tense: ago
The infinite form: agere
The perfect tense: egi
The 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0tag:blogger.com,1999:blog-2889082294597625339.post-65612154228779539562018-04-02T02:15:00.000-07:002019-04-22T12:08:39.743-07:00In principioGenesis 1:1-3
In principio creavit Deus caelum et terram.
Terra autem erat inanis et vacua.
Et tenebrae super faciem abyssi et spiritus Dei ferebatur super aquas.
Dixitque Deus fiat lux et facta est lux.
Vowels:
A: Short, "uh" in about, long; Long, "ah" in father. Latin A never pronounced as cat or kate in English.
E: Short, "eh" in get, bed; Long, "ay" in cake
I: Short, "ih" in kin, tin; 呵呵http://www.blogger.com/profile/10034675933395730939noreply@blogger.com0