A finite element method of the self-consistent field theory on general curved surfaces

TL;DR

This paper introduces a finite element method for accurately simulating self-assembled block copolymer structures on arbitrary curved surfaces using self-consistent field theory, expanding the scope of numerical analysis in soft matter physics.

Contribution

It develops a rigorous surface finite element method and an adaptive surface optimization approach for studying block copolymer phases on general curved surfaces.

Findings

Efficient numerical method for complex curved surfaces.

Ordered structures consistent with known results.

Capable of analyzing phase behavior on arbitrary surfaces.

Abstract

Block copolymers provide a wonderful platform in studying the soft condensed matter systems. Many fascinating ordered structures have been discovered in bulk and confined systems. Among various theories, the self-consistent field theory (SCFT) has been proven to be a powerful tool for studying the equilibrium ordered structures. Many numerical methods have been developed to solve the SCFT model. However, most of these focus on the bulk systems, and little work on the confined systems, especially on general curved surfaces. In this work, we developed a linear surface finite element method, which has a rigorous mathematical theory to guarantee numerical precsion, to study the self-assembled phases of block copolymers on general curved surfaces based on the SCFT. Furthermore, to capture the consistent surface for a given self-assembled pattern, an adaptive approach to optimize the size of the general curved surface has been proposed. To demonstrate the power of this approach, we investigate the self-assembled patterns of diblock copolymers on several distinct curved surfaces, including five closed surfaces and an unclosed surface. Numerical results illustrate the efficiency of the proposed method. The obtained ordered structures are consistent with the previous results on standard surfaces, such as sphere and torus. Certainly, the proposed numerical framework has the capability of studying the phase behaviors on general surfaces precisely.

Tables

Table 1. Table 1: The calculated information of SCFT simulations on a series of subdivision meshes, where “–” means divergence.

# of Node	# of Elem.	$\Gamma$	$\Gamma$-diff.	Hamiltonian	Ham. diff.
162	320	–	–	–
642	1280	3.707	8.07e-1	-3.173
2562	5120	3.579	1.29e-1	-3.161	1.22e-2
10242	20480	3.566	1.13e-2	-3.165	4.89e-3
40962	81920	3.563	2.40e-3	-3.167	1.42e-3

Table 2. Table 2: The iteration details of the SCFT simulation when computing the 12-microdomain spotted phase on the sphere surface used 10242 nodes and 20480 elements.

Iter. of SCFT	Iters. of FSP	$\Gamma$	$\Gamma$-Diff.	Hamiltonian	Ham. Diff.
1	169	2.90	-	-3.10648180	-
2	361	3.372922	4.729e-1	-3.16079157	5.431e-2
3	73	3.656954	2.843e-1	-3.16437826	3.587e-3
4	50	3.530368	1.266e-1	-3.16530810	9.298e-4
5	46	3.580692	5.032e-2	-3.16538583	7.773e-5
6	44	3.559392	2.130e-2	-3.16543307	4.723e-5
7	14	3.568215	8.823e-3	-3.16529309	1.399e-4
8	17	3.564173	4.043e-3	-3.16546705	1.739e-4
9	15	3.566281	2.018e-3	-3.16540431	6.273e-5
10	15	3.565232	1.049e-3	-3.16543986	3.554e-5
11	13	3.565744	5.120e-4	-3.16542586	1.399e-5
12	8	3.565524	2.212e-4	-3.16542759	1.737e-6
13	2	3.565564	4.027e-5	-3.16543001	2.422e-6
14	2	3.565581	1.061e-6	-3.16543102	1.008e-6
15	2	3.565591	1.000e-6	-3.16543136	3.377e-7

Table 3. Table 3: The initial scale factor Γ 0 subscript Γ 0 \Gamma_{0} and optimal Γ o ​ p ​ t subscript Γ 𝑜 𝑝 𝑡 \Gamma_{opt} of the sphere surface for different ordered pattern as shown in Fig.  7 .

	Fig. 7 (a)	Fig. 7 (b)	Fig. 7 (c)	Fig. 7 (d)	Fig. 7 (e)
${\Gamma }_0$	2.9	11.65	8	7.273	8
${\Gamma }_{opt}$	3.565591	10.529	7.273	5.468	6.408

Table 4. Table 4: The relationship between the number of microdomains of a spotted pattern and the optimal scale factors of sphere when χ ​ N = 25.0 𝜒 𝑁 25.0 \chi N=25.0 and f = 0.20 𝑓 0.20 f=0.20 .

${\Gamma }_{opt}$	$3.430$	$4.593$	$5.136$	$5.253$	$7.411$	$8.384$	$9.271$	$9.841$	$10.53$
# of Spots	$12$	$22$	$28$	$29$	$56$	$72$	$89$	$101$	$116$

Table 5. Table 5: The initial R 0 subscript 𝑅 0 R_{0} , r 0 subscript 𝑟 0 r_{0} and optimal scale factors Γ o ​ p ​ t subscript Γ 𝑜 𝑝 𝑡 \Gamma_{opt} for different ordered structure on the torus surface.

	Fig. 8 (a)	Fig. 8 (b)	Fig. 8 (c)	Fig. 8(d)	Fig. 8(e)	Fig. 8(f)	Fig. 8(g)
$\gamma$	8	2	1.389	2	2	2	2
$R_0$	8	4	3.333	4	4	4	4
$r_0$	1	2	2.4	2	2	2	2
${\Gamma }_{opt}$	1.066	1.044	1.045	0.969	1.023	1.119	1.223

Table 6. Table 6: The initial Γ 0 subscript Γ 0 \Gamma_{0} and optimal scale factors Γ o ​ p ​ t subscript Γ 𝑜 𝑝 𝑡 \Gamma_{opt} of double torus, heart and orthocircle surfaces used in the surface simulations for diblock copolymers.

	Double torus		Heart		Orthocircle
	Fig. 9 (a)	Fig. 9 (d)	Fig. 9 (b)	Fig. 9 (e)	Fig. 9 (c)	Fig. 9 (f)
${\Gamma }_0$	12	12	5	5	5	5
${\Gamma }_{opt}$	12.308	14.161	5.022	5.403	5.143	4.939

Table 7. Table 7: The initial scale factor Γ 0 subscript Γ 0 \Gamma_{0} and optimal Γ o ​ p ​ t subscript Γ 𝑜 𝑝 𝑡 \Gamma_{opt} , in computing ordered patterns of diblock copolymers on parabolic surface.

	Fig. 10 (a)	Fig. 10 (b)
${\Gamma }_0$	10	10
${\Gamma }_{opt}$	10.310	10.558

Equations

$$\mathcal{S}=\{\mathbf{x}\in U|d(\mathbf{x})=0\},$$

$$\mathbf{n}=\frac{\nabla d(\mathbf{x})}{|\nabla d(\mathbf{x})|},$$

$$\mathcal{P}(\mathbf{x})=\mathbf{x}-d(\mathbf{x})\mathbf{n},$$

$${\nabla }_{S}v=\nabla v-(\nabla v\cdot \mathbf{n})\mathbf{n}=\mathbf{P}(\nabla v),$$

$${\nabla }_S\cdot \mathbf{v}(\mathbf{x})=\nabla \cdot \mathbf{v}-{\mathbf{n}}^t\nabla \mathbf{v}\mathbf{n},$$

$${\Delta }_{\mathcal{S}}v={\nabla }_{\mathcal{S}}\cdot ({\nabla }_{\mathcal{S}}v)=\Delta v-(\nabla v\cdot \mathbf{n})(\nabla \cdot \mathbf{n})-{\mathbf{n}}^t{\nabla }^{2}v\mathbf{n},$$

$$H^1(\mathcal{S})=\{v\in L^2(\mathcal{S})|{\nabla }_{S}v\in (L^2(\mathcal{S}){)}^3\}.$$

$${\phi }_i({\mathbf{x}}_j)=\{\begin{array}{ll}1,& \text{ if }i=j,\\ 0,& \text{ if }i\ne j.\end{array}$$

$${\nabla }_{S_h}{v}_h{|}_{{\tau }_h}:=\nabla v_h-(\nabla v_h\cdot {\mathbf{n}}_h){\mathbf{n}}_h=(\mathbf{I}-{\mathbf{n}}_h{\mathbf{n}}_h^t)\nabla v_h={\mathbf{P}}_h\nabla v_h.$$

$$x(\theta ,\varphi )$$

$$y(\theta ,\varphi )$$

$$z(\theta ,\varphi )$$

$$x^2(x^2-1)[x^2(x^2-1)+2y^2]+y^4+z^2-0.04=0.$$

$$(x-z^2{)}^2+y^2+z^2-1=0.$$

$$[(x^2+y^2-1{)}^2+z^2][(y^2+z^2-1{)}^2+x^2][(z^2+x^2-1{)}^2+y^2]-0.075^2[1+3(x^2+y^2+z^2)]=0.$$

$$H[w_{+},w_{-}]=$$

$$\frac{\delta H}{\delta w_{+}(\mathbf{x})}={\varphi }_A(\mathbf{x})+{\varphi }_B(\mathbf{x})-1=0,$$

$$\frac{\delta H}{\delta w_{-}(\mathbf{x})}=\frac{2w_{-}(\mathbf{x})}{\chi N}-[{\varphi }_A(\mathbf{x})-{\varphi }_B(\mathbf{x})]=0,$$

$$Q=\frac{1}{|\mathcal{S}|}\int \mathrm{d}\mathbf{x}q(\mathbf{x},1)=\frac{1}{|\mathcal{S}|}\int \mathrm{d}\mathbf{x}q(\mathbf{x},s)q^{†}(\mathbf{x},s),\forall s\in [0,1],$$

$${\varphi }_A(\mathbf{x})=\frac{1}{Q}{\int }_0^f\mathrm{d}sq(\mathbf{x},s)q^{†}(\mathbf{x},s),$$

$${\varphi }_B(\mathbf{x})=\frac{1}{Q}{\int }_f^1\mathrm{d}sq(\mathbf{x},s)q^{†}(\mathbf{x},s).$$

$$\frac{\partial }{\partial s}q(\mathbf{x},s)$$

$$q(\mathbf{x},0)$$

$$w(\mathbf{x},s)$$

$$\frac{\partial }{\partial s}{q}^{†}(\mathbf{x},s)$$

$$q^{†}(\mathbf{x},1)$$

$$w^{†}(\mathbf{x},s)$$

$$\frac{\partial }{\partial t}{w}_{+}(\mathbf{x},t)$$

$$\frac{\partial }{\partial t}{w}_{-}(\mathbf{x},t)$$

$$\underset{\mathcal{S}}{\min }\underset{w_{+}}{\max }\underset{w_{-}}{\min }H[w_{+}(\mathbf{x}),w_{-}(\mathbf{x}),\mathcal{S}].$$

$${\mathcal{S}}_{\Gamma }=\{\Gamma \cdot \mathbf{x}:\mathbf{x}\in {\mathcal{S}}_0\},$$

$$\underset{\Gamma }{\min }\underset{w_{+}}{\max }\underset{w_{-}}{\min }H[w_{+}(\mathbf{x}),w_{-}(\mathbf{x}),\Gamma ].$$

$${\left(\frac{\partial }{\partial s}q,v\right)}_{\mathcal{S}}=-({\nabla }_{\mathcal{S}}q,{\nabla }_{\mathcal{S}}v{)}_{\mathcal{S}}-(wq,v{)}_{\mathcal{S}},\text{ for all }v\in H^1(\mathcal{S}),$$

$${\left(\frac{\partial }{\partial s}{q}_h,v_h\right)}_{{\mathcal{S}}_h}=-({\nabla }_{{\mathcal{S}}_h}{q}_h,{\nabla }_{{\mathcal{S}}_h}{v}_h{)}_{{\mathcal{S}}_h}-(w_h{q}_h,v_h{)}_{{\mathcal{S}}_h},\text{ for all }{v}_h\in {\mathcal{V}}_h,$$

$$M\frac{\partial }{\partial s}\mathbf{q}(s)=-(A+F)\mathbf{q}(s),$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Geometry plays a key role in many scientific fields [1], including Hamiltonian mechanics with constraints, general relativity, quantum mechanics, quantum field theory, the arrangement of electrons on a sphere (Thomson problem), defect motion on a curved surface, polymer field theory, generic reaction-diffusion model occurring in chemistry and biology. Recent years, the microphase separation of block copolymer under various types of geometrical confinements both in the bulk and on surfaces has been attracted tremendous attention. Many novel patterns of block copolymer emerge from rearrangements of traditional ordered phases when adapting to these restrictions [2, 3, 4, 5]. The self-assembly of block copolymers under geometrical constraint provides an efficient means on the nanometer scale for such applications as the construction of high-capacity data storage devices, waveguides, quantum dot arrays, dielectric mirrors, nanoporous membranes, nanowires, and interference lithography [5, 6]. The curvature effects on the self-assembled behavior of block copolymers confined to a curved surface, however, are still far from being fully understood [7, 8, 9]. The main obstacle may be that it still lacks a general theoretical framework which can properly describe the microphase separation of block copolymers on general curved surfaces.

Theoretically, the self-consistent field theory (SCFT) [10] proposes a successful framework in the study of the self-assembling of block copolymers. It has predicted a variety of ordered bulk phases that have been observed in experiments, e.g., lamella, hexagonal cylinder, sphere and gyroid phases of diblock copolymers [11], and many intricate ordered phases of triblock copolymers [12]. Mathematically, the SCFT model is a very complicated variational problem, possessing many unsatisfactory features, such as saddle-point, nonlinearity, multi-solutions, and multi-parameters. Analytically solving this model goes beyond current technologies. An alternating approach is the numerical technique.

The numerical methods for solving SCFT model mainly consists of four components, namely screening initial values [13, 14, 15], solving time-dependent partial differential equations (PDEs), evaluating (monomer) density operators, and finding saddle-points via iteration methods [7, 16, 17]. In this paper, we focus on the spatial discretization of the time-dependent PDEs in the SCFT model. In flat spaces, many numerical approaches have been developed for spatial variables which can generally be divided into two classes of discretization schemes. The first type is the subspace discretization method which discretizes equations in a constrained space based on specific properties. For example, according to the crystallographic symmetry of a given pattern in microphase-separated block copolymers, the PDEs in SCFT model can be expanded in terms of a set of symmetric basis functions. One of the typical representatives is the symmetry-based Fourier method proposed by Masten and Schick [11]. The second type is numerically approximating PDEs in the whole space. These methods can be carried out either in the real space [18] or in the Fourier-space [19]. It has also been demonstrated that the whole-space discretization methods are able to discover new patterns [18, 19]. In recent years, an efficient Fourier pseudospectral method has been introduced to solve PDEs in SCFT model [20, 21]. It fully takes advantage of the best performance of real space and Fourier-space and reduces the computational complexity to $O(M\log M)$, with the number of degrees $M$, based on the Fast Fourier Transformation (FFT).

Besides space discretization, many time discretization schemes, such as the operator-splitting method, the Crank-Nicokson scheme, and linear multi-step approaches, have been used in the time discretization [10]. The total computational complexity in solving time-dependent PDE should consider both the time and spatial discretization schemes.

From flat spaces to curved surfaces, the Laplace operator in the PDEs becomes the Laplace-Beltrami operator which is the heat kernel on the manifold. Numerical methods should be developed to discretize the differential operator on a curved surface. There are two approaches to represent the curved surface. The first one is viewing the curved surface as a level set through embedding it into a high dimensional flat space. In another viewpoint, the curved surface is represented by the inner coordinates of the manifold. From the first viewpoint, the mask method, which uses a large enough domain in high dimensional flat space to cover the curved surface, can be chosen to attack curved surface problems [22]. These methods developed in the flat space can be used to solve the curved surface problem. Then these results on a manifold can be approximated by restricting the results of high dimensional flat space on the curved surface. From another viewpoint, one shall directly discretize the function defined on a curved surface. For some special curved surfaces, such as the sphere surface, one can expand the function defined on the sphere by the eigenfunctions, i.e., the spherical harmonic functions, of the Laplace-Beltrami operator [7]. For a general case, however, it is impossible to find out the eigenfunctions of the Laplace-Beltrami operator. It causes inconvenience in solving PDEs in the SCFT model on general curved surfaces. In 2014, Li el al. [8] mimicked the process of finite difference method and developed an extended spherical ADI finite difference [9] to solve SCFT model on a general curved surface. However, the precision of this method can not be guaranteed.

It is well-known that as an effective numerical method for solving PDE, the finite element method plays an important role in modern scientific and engineering computing. Based on variational principle, finite element method subdivides the definition domain into small and simple regions, such as triangle, quadrilateral, tetrahedron, or hexahedron, etc., and then simple algebraic equations on each small region can be created. At last these simple equations are assembled into a large sparse algebraic system to approximate the PDE problem. The advantage of the finite element method is that it can handle complex geometrical domains with complex boundaries. Based on the similar principal, surface finite element method has been developed to solve PDEs on a general manifold [24, 25, 26, 27, 28, 29, 30]. The approach was firstly proposed by Dziuk to solve the Laplace-Beltrami equation on arbitrary surfaces using linear finite element [24]. Then Dziuk and Elliott [25] applied the linear surface finite element to parabolic equations, and gave error bounds in the case of semi-discretization in space for a fourth order linear problem. Our previous work [26] generalized the superconvergence results and several gradient recovery approaches of linear finite element methods from flat spaces to general curved surfaces for the Laplace-Beltrami equation with mildly structured triangular meshes. In this work, we will continue to extend the surface finite element method to study the microphase separation of block copolymer on general curved surface based on the SCFT model.

On the other hand, the surface size also affects self-assembling patterns. For example, the self-assembled periodic structures are affected by the computational domain in flat space [11]. On the sphere surface, the radius can affect the number of microdomains of a cylindrical phase as well as its energy value, and phase behavior of lamellar phases. Therefore, for a given self-assembled phase, there exists the most appropriate size of curved surface. To capture the optimal size, in this work, we will propose an adaptive method to optimize the surface, and further obtain accurate energy value for a given ordered structure.

The remaining portion is arranged as follows. In Sec. 2, we introduce the preliminary knowledge of surface finite element method and curved surfaces to be used in this paper. In Sec. 3 we derive the SCFT model on the general curved surface. The numerical algorithm of linear surface finite element for solving SCFT model is given in Sec. 4. The efficiency of our method and the self-assembled structures of diblock copolymers on several curved surfaces are presented in Sec. 5. In the Sec. 6, we will draw the conclusion and outline the further work.

2 Preliminaries

In this section, we firstly introduce some background knowledge about the surface finite element method, and refer to  [30] for details. Then we present the surfaces used in the following calculations and the approach of how to generate optimal meshes on these surfaces.

2.1 Surface finite element method

Let $\mathcal{S}$ be a two-dimensional, compact, and $\mathcal{C}^{2}$-hypersurface in $\mathbb{R}^{3}$. Here we assume there exists a smooth level set function $d(\mathbf{x})$ that can describe the surface $\mathcal{S}$ as

[TABLE]

where $U$ is open subset of $\mathbb{R}^{3}$ in which $\nabla d(\mathbf{x})\not=0$, and $d\in\mathcal{C}^{2}(U)$. We define the unit normal vector on $\mathcal{S}$ by

[TABLE]

which means that the orientation of $\mathcal{S}$ is fixed through normal vector $\bm{n}$. For any $\mathbf{x}\in U$, one can define a projection

[TABLE]

such that $\mathcal{P}(\mathbf{x})\in\mathcal{S}$.

For $v(\mathbf{x})\in\mathcal{C}^{1}(\mathcal{S})$, since $\mathcal{S}$ is $\mathcal{C}^{2}$, we can extent $v$ to $\mathcal{C}^{1}(U)$ and still denote the extension by $v$. The tangential gradient of $v$ on $\mathcal{S}$ can be written as

[TABLE]

where $\bm{P}(\mathbf{x})=(I-\bm{n}\bm{n}^{t})(\mathbf{x})$ is the projection operator to the tangent plane at a point $\mathbf{x}\in\mathcal{S}$. Notice that we use the extension of $v$ to define the surface gradient. However, it can be shown that $\nabla_{S}v$ depends only on the value of $v$ on $\mathcal{S}$ but not on the extension [29]. Namely, $\nabla_{S}$ is an intrinsic operator.

Similarly, for a vector field $\bm{v}\in(\mathcal{C}^{1}(\mathcal{S}))^{3}$, we can also extend it to $(\mathcal{C}^{1}(U))^{3}$ and define the tangential divergence operator of $\bm{v}$ on $\mathcal{S}$ as

[TABLE]

The Laplace-Beltrami operator on $\mathcal{S}$ reads as follows:

[TABLE]

provided $v(\mathbf{x})\in\mathcal{C}^{2}(\mathcal{S})$ and $\nabla^{2}v$ is the Hessian matrix of $v$ (suitably extended as a $\mathcal{C}^{2}(U)$ matrix function). The Sobolev spaces on surface $\mathcal{S}$ can be defined as:

[TABLE]

Let $\mathcal{S}$ be approximated by a triangular mesh $\mathcal{S}_{h}$ with node set $\mathcal{N}_{h}=\{\mathbf{x}_{i}\}_{i=1}^{N}$ and triangle set $\mathcal{T}_{h}=\{\tau_{h}\}$. We assume that these triangles are shape-regular and quasi-uniform of a diameter $h$ and their vertices lie on $\mathcal{S}$. For any $\tau_{h}\in\mathcal{T}_{h}$, let $\bm{n}_{h}$ be the unit outward normal vector of $S_{h}$ on $\tau_{h}$. Let $V_{h}$ be the continuous piecewise linear finite element space on $S_{h}$, with linear Lagrangian basis functions $\{\varphi_{i}(\mathbf{x})\}_{i=1}^{N}$ defined on $S_{h}$. $\varphi_{i}(\mathbf{x})$ is the piecewise linear function on each triangle face $\tau_{h}$,

[TABLE]

For $v_{h}\in\mathcal{C}(S_{h})$ and $v_{h}|_{\tau_{h}}\in\mathcal{C}^{1}(\tau_{h})$, we have

[TABLE]

2.2 Surface representation and mesh generation

In this subsection, we will discuss how to generate high-quality finite element meshes on surfaces. A high-quality mesh on a surface means: the triangle elements are all almost equilateral and more mesh nodes located at the place with large curvature [26]. In this work, five closed curved surfaces including sphere, torus, double torus, heart, and orthocircle surfaces, and an unclosed parabolic surface are considered. In general, there are three ways to represent a continuous surface in $\mathbb{R}^{3}$: explicitness, implicitness, and parametrization.

Sphere surface

Any point on the unit sphere surface satisfies the implicit equation $\sqrt{x^{2}+y^{2}+z^{2}}-1=0$. A good triangular mesh on sphere can be easily generated by uniformly refining an icosahedron and projecting the new mesh nodes, which are the old mesh edge center, onto the sphere. The refinement procedure is illustrated in Fig. 1.

Torus surface

The second example is the torus which can be defined parametrically by

[TABLE]

where $\theta,\phi\in[0,2\pi]$, $R$ is the major radius which is the distance from the center of the tube to the center of the torus, $r$ is the minor radius which is the radius of the tube. The ratio of $r/R$ is the so-called ”aspect ratio”. With the above parameteric representation, the triangular mesh on torus can be generated by mapping a triangular mesh on the square domain $[0,2\pi]^{2}$, see Fig. 2.

The implicit expressions of double torus, heart, and orthocircle are given as follows:

Double torus surface

[TABLE]

Heart surface

[TABLE]

Orthocircle surface

[TABLE]

For these complicated surfaces, it is difficult to generate a high quality surface triangluar mesh as done for sphere or torus. To get the high quality mesh, we firstly use the 3D Surface Mesh Generation package in CGAL [31] to generate an initial triangluar mesh (see Fig. 3), then further optimized these initial mesh through the Centroidal Voronoi Tessellation (CVT) technique [32] (see Fig. 4).

Parabolic surface

Finally, for the unclosed parabolic surface, we firstly generate a good quality triangluar mesh on the unit disk, then lift the mesh nodes onto the parabolic surface, as shown in Fig. 5.

3 SCFT on general curved surfaces

In this section, we will introduce a standard field theory model for an incompressible $AB$ diblock copolymer melt on general curved surfaces. We consider $n$ $AB$ diblock copolymers with the polymerization of $N$ confined on surface $\mathcal{S}$ with a total surface area $|\mathcal{S}|$. The volume fraction of block $A$ is $f$ and that of block $B$ is $1-f$. The two blocks disfavor with each other characterized by the Flory-Huggins parameter $\chi N$. We assume that the statistical segment length and volume of two blocks are equal, i.e., $b_{A}=b_{B}=b$ and $v_{A}=v_{B}=v_{0}$. A characteristic length of the copolymer chain can be defined by the unperturbed radius of gyration, which is used as the unit of length so that all spatial lengths are presented in units of $R_{g}=b\sqrt{N/6}$. With the incompressible melt assumption, the average segment density is uniform in space and given by $v_{0}=1/\rho_{0}=V/(nN)$, $\rho_{0}$ is the number density. Here we use a continuous Gaussian chain to model the block copolymer [10, 33]. Based on the statistical mechanics, and using Hubbard-Stratonovich transformation and infinite dimensional Fourier transformation [10, 34, 35], one can derive the field-based model to describe the phase behavior of the block copolymer system. However, directly solving the field-based model is nontrivial due to the very complicated nature of the functional integral exhibited in the partition function [10]. In order to simplify the field model, we often use an analytic approximation technique called the mean-field (saddle-point) approximation which ignores the field fluctuations and assumes that the functional integral is dominated by a single field configuration. This approximation method become accurate when the dimensionless chain concentration $C=\rho_{0}R_{g}^{3}/N$ goes to infinity. In the case of high molecular weight block copolymer melts, where $C$ can be very large, the approximation will make sense. From the viewpoint of asymptotic theory, the functional integral, i.e., the partition function, can be approximated by the integrand function at the critical point (actually the saddle point) when the parameter $C$ becomes very large. Then one can obtain the mean-field model. Specifically, within the mean-field approximation, the effective Hamiltonian of this system is

[TABLE]

In the above equation, $w_{+}(\mathbf{x})$ and $w_{-}(\mathbf{x})$ are the fluctuating pressure and exchange chemical potential fields, respectively. The pressure field enforces the local incompressibility, while the exchange chemical potential is conjugate to the density operators. $Q$ is the single chain partition functional subjected to external fields $w_{+}$ and $w_{-}$. $\mathrm{d}\mathbf{x}$ is the area element of the surface. First-order variations of the free energy functional with respect to the fields lead to the following SCFT equations,

[TABLE]

where $\phi_{A}(\mathbf{x})$ and $\phi_{B}(\mathbf{x})$ are the monomer densities of block $A$ and $B$, respectively. Then using the continuous Gaussian chain model,

[TABLE]

The forward propagator, $q(\mathbf{x},s)$, represents the probability weight that the chain of contour length $s$ has its end at surface position $\mathbf{x}$, where the variable $s$ is used to parameterize each copolymer chain. $s=0$ represents the tail of block $A$ and $s=f$ is the junction between block $A$ and $B$. From the continuous Gaussian chain model, $q(\mathbf{x},s)$ satisfies the PDE

[TABLE]

The backward propagator, $q^{{\dagger}}(\mathbf{x},s)$, satisfies an almost identical PDE of Eqn. (11), except that it represents the probability weight from $s=1$ to $s=0$. In particular, it satisfies

[TABLE]

4 Methodology

The equilibrium point, as well as the saddle point, of the SCFT model corresponds to the ordered phase. And finding the saddle points requires iteration methods. The previous study has been shown that the effective Hamiltonian (5) can reach its local minima along the exchange chemical field $w_{-}$ and maxima along the pressure field $w_{+}$ [34, 35]. Although it still lacks a rigorous theoretical guarantee, the numerical behavior has demonstrated its effectiveness. Based on this observation, we can use an alternating direction scheme to obtain the saddle points. Specifically, we introduce a fictitious time variable $t$, and at each time step, the saddle-point search method is given as

[TABLE]

Clearly, Eqns. (6) and (7) are satisfied when Eqns. (13) and (14) are stationary.

Within the standard framework of SCFT, finding saddle point (FSP) of $w_{\pm}$ can be obtained by the following iteration scheme:

Step 1

Given initial estimations of fields $w_{\pm}(\mathbf{x},0)$.

Step 2

Compute forward and backward propagator operators $q(\mathbf{x},s)$ and $q^{\dagger}(\mathbf{x},s)$ on a general curved surface (see Sec. 4.1).

Step 3

Obtain $Q$, $\phi_{A}(\mathbf{x})$ and $\phi_{B}(\mathbf{x})$ by integral equations (see Sec. 4.4 and Sec 4.5), and evaluate the value of effective Hamiltonian $H$.

Step 4

Update fields $w_{+}(\mathbf{x},t)$ and $w_{-}(\mathbf{x},t)$ by Eqns. (13) and (14) using iteration methods (see Sec. 4.3).

Step 5

Repeat steps 2-4 until a convergence criterion is met.

In the polymeric systems, each self-assembled microphase possesses its characteristic scale due to the statistical length of the polymer chain. In the flat space, the size and shape of the computational domain can affect the energy value of self-assembled structures [11, 15] and the properties of the polymer materials, such as the mechanical property [36, 37]. As in the flat domain cases, the ordered pattern formed on the curved surface has its characteristic scale. In contrast, the surface size also affects the self-assembled patterns. For example, the radius of sphere affects the self-assembled phases of diblock copolymers and their energy values [7]. However, it still lacks a systematical method to select an optimal size for a self-assembled structure on a general curved surface. Here, we will propose a theoretic tool to efficiently handle this problem, i.e., an adaptive approach of optimizing the size of a general curved surface for a given ordered structure.

In order to obtain the optimal surface size, the effective Hamiltonian of SCFT shall be viewed as a function of the surface size, as well as a functional of the field functions. Therefore, the completed optimization problem of solving SCFT equations becomes

[TABLE]

In practice, we parameterize the surface $\mathcal{S}$ by

[TABLE]

where $\mathcal{S}_{0}$ is the initial surface, $\Gamma>0$ is a scale factor which describes the surface size. For example, the parameter $\Gamma$ of a sphere is its radius. Then we can optimize the surface size by adjusting the scale factor $\Gamma$.

Therefore the optimization problem (15) becomes

[TABLE]

It can be solved by the following iteration procedure, here we name SCFT iteration, including finding the saddle point of SCFT and adaptively optimizing the scale factor $\Gamma$ of a given surface,

Step 1

Give proper parameters $\chi N$, $f$, a surface $\mathcal{S}$, and reasonable initial distributions of $w_{\pm}$.

Step 2

Fix $\mathcal{S}$, finding the saddle point of SCFT through FSP procedure and obtaining the effective Hamiltonian.

Step 3

Fix $w_{\pm}$, optimizing the size of $\mathcal{S}$ by the method described in Sec. 4.2, and evaluating the value of effective Hamiltonian.

Step 4

Repeat steps 2-3 until the effective Hamiltonian discrepancy is lower than a given convergence criterion.

4.1 Surface finite element method for PDE

In the iteration scheme outlined above, the most time consuming step is solving the PDEs of (11) and (12). For a special surface, such as a sphere, the global basis of spherical harmonics can be used to expand the spatial functions. Therefore the (pseudo) spherical harmonic spectral method can be used to solve PDEs on sphere. For a general curved surface, however, the global basis might not exist. Therefore, the methods based on local basic functions shall be developed to solve SCFT model on general curved surface, such as the finite volume or element method. In this work, we focus our attention on the development of the finite element method. We discretize the Laplace-Beltrami operator in (11) and (12) by the linear surface finite element method.

Here we just present the finite element discretization for Eqn. (11). For Eqn. (12) of backward propagator $q^{{\dagger}}$, one can obtain the discretization scheme similarly. Firstly, we rewrite (11) into the variational formulation: find $q$ in $H^{1}(\mathcal{S})$, which satisfies

[TABLE]

where $(\cdot,\cdot)_{S}$ represents the $L^{2}$ inner product on $\mathcal{S}$.

Replacing the infinite dimensional space $H^{1}(\mathcal{S})$ by the finite dimension space $\mathcal{V}_{h}$, we obtain the linear finite element discretization: find $q_{h}=\sum_{i=1}^{N}q_{i}(s)\varphi_{i}(\mathbf{x})$ in $\mathcal{V}_{h}$, which satisfies

[TABLE]

where $(\cdot,\cdot)_{\mathcal{S}_{h}}$ represents the $L^{2}$ inner product on the triangular mesh $\mathcal{S}_{h}$, and $w_{h}(\mathbf{x},s)=\sum_{i=1}^{N}w(\mathbf{x}_{i},s)\varphi_{i}(\mathbf{x})$ is the linear interpolation of $w(\mathbf{x},s)$. Replacing $v_{h}$ by $\varphi_{j}$, $j=1,2,\cdots,N$, we can have the matrix form of (19),

[TABLE]

where

[TABLE]

and

[TABLE]

For the discretization of $s$-direction, we use the Crank-Nicolson method,

[TABLE]

where $\Delta s$ is the time step size. It should be noted that any stable time discretization scheme can be used here. Combining with the above two steps, we obtain the full discretization scheme

[TABLE]

If the linear algebraic system (22) is relatively small, we solve it with the direct method. Otherwise the algebraic multigrid method can be applied [38]. Furthermore, the order of the above discretizations for the approximation to the solution of the modified diffusion equations is $O(\Delta s^{2}+h^{2})$, where $h$ is triangular surface mesh size [28].

4.2 Adaptively optimizing the curved surface

When fixed fields $w_{\pm}(\mathbf{x})=w_{\pm}^{f}(\mathbf{x})$, the optimization problem of (17) becomes

[TABLE]

Any appropriate optimization method can be chosen to solve this problem, such as the simplest steepest descend method, i.e., $d\Gamma/dt=-dH/d\Gamma$. However, due to the complexity of SCFT, it is hardly written down the derivative of the SCFT’s Hamiltonian with respect to the scale factor of a curved surface $dH/d\Gamma$ analytically. Alternatively, we can approximate $dH/d\Gamma$ numerically through the finite difference method. To improve the effectiveness of algorithm, in practical implementation, we adapt the nonlinear conjugate-gradient (CG) method to minimize objective function (23) [39], i.e.,

[TABLE]

where $d^{k}$ is the conjugate gradient direction in $k$th step. $\alpha_{k}$ is the step size obtained by the linear search approach. Meanwhile, the restart technology has been used in nonlinear CG method to avoid the oscillation and improve the convergent rate.

4.3 Iteration method of finding saddle point

The iteration methods to update the fields are dependent on the mathematical structure of SCFT. An important fact is that the effective Hamiltonian (5) of diblock copolymers can reach its local minima along the exchange chemical field $w_{-}(\mathbf{x})$, and reach the maxima along the pressure field $w_{+}(\mathbf{x})$ [10].

Then a series of accelerate convergence methods, such as the semi-implicit method, can be designed to find the saddle point. More general results for multicomponent polymeric systems with more than two chemically distinct blocks, including a simple analysis of SCFT model and a class of efficient iteration methods of finding saddle points, can be found in Ref. [34]. It should be noted that the semi-implicit scheme is obtained through the asymptotic expansion and Fourier transformation. Based on similar techniques, it can be extended to sphere surface problems using spherical harmonic transformation [7]. However, it is hardly to find a representation of spatial functions defined on a general curved surface using global basis.

Therefore we choose the alternative direction explicit Euler method to update the fields $w_{\pm}(\mathbf{x})$. Specifically, this approach is expressed as

[TABLE]

4.4 Integral formula along $s$-direction

A modified numerical integration formula for a closed interval is chosen to evaluate integral equations (9)-(10) that can guarantee fourth-order precision in $s$-direction whether the number of discretization points is even or odd [40].

[TABLE]

where $n$ is the number of discretization nodes along $s$-direction, $\Delta s$ is a constant step, and $s_{n}=n\Delta s$. $f_{k}=f(k\Delta s)$, $k=0,1,\dots,n$.

4.5 Surface integral

Given a linear finite element function $v_{h}=\sum_{i=1}^{N}v_{i}\varphi_{i}(\mathbf{x})$ on $\mathcal{S}_{h}$, we can calculate the integration of $v_{h}$ on $\mathcal{S}_{h}$ from Eqn. (1),

[TABLE]

where $\varphi_{i}(\mathbf{x})$ is a linear function about $\mathbf{x}$ on each triangle face, $v_{i}$, $v_{j}$ and $v_{k}$ are the function values on the three vertices of a triangle $\tau_{h}=(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k})$, respectively, and $|\tau_{h}|$ is the area of $\tau_{h}$. Notice that, for linear function on plain triangle, the above single point numerical integration formula is accurate. Given a function $v\in H^{1}(\mathcal{S})$ and its interpolation $v_{I}$ on $\mathcal{S}_{h}$, the error estimate can be obtained by Cauchy-Schwarz inequality and interpolation error estimate [24]:

[TABLE]

where $\bar{v}(\mathbf{x})=v(\mathcal{P}(\mathbf{x}))$ and $|\mathcal{S}_{h}|$ is the area of $\mathcal{S}_{h}$.

5 Numerical Results

In this section, we will use five closed curved surfaces including sphere, torus, double torus, heart, orthocircle, and an unclosed parabolic surface to demonstrate the efficiency of our proposed numerical method. To ensure the accuracy, we set the substep of contour length $s$ to $0.005$ in the following simulations. The convergence criterion of FSP is

[TABLE]

and that of SCFT iteration is the change of effective Hamiltonian smaller than $10^{-6}$.

5.1 Efficiency of the numerical method

Since it is hard to obtain a nontrivial analytical solution of SCFT on a general curved surface, we have to verify the convergence of the numerical solution of our proposed method by running the simulations on a series of subdivision meshes. Here we take the 12-microdomain spots structure on the sphere with model parameters $\chi N=30.0$, $f=0.20$, $\Gamma=2.9$ as an example to show the convergence of our method, see Tab. 1.

In these series of simulations, the initial fields on each node $(x,y,z)$ is $w_{-}=\chi N\sin(5\theta)$, $w_{+}=0$, where $\theta=\arctan(y/x)\in[0,2\pi]$. The initial surface is a unit sphere with scale factor $\Gamma=2.9$. We uniformly refine the initial mesh and the number of nodes and triangular elements can be found in the first and second columns of Tab. 1. From the table, one can find that, along with the refinement of the mesh, the value of effective Hamiltonian is indeed convergent, and furthermore, the scale factor $\Gamma$ is adaptively optimized and convergent. The finally convergent morphology can be found in Fig 7 (a). In order to further demonstrate the iteration procedure, we give the iteration details of the SCFT iteration when using $10242$ nodes and $20480$ elements as shown in the fifth row of Tab. 1. The information is given in Tab. 2. The effective Hamiltonian and its discrepancy with respect to iteration are given in Fig. 6. It costs $15$ SCFT iterations, and $831$ FSP iterations. The scale factor $\Gamma$ is optimized from $2.9$ to $3.5656$. Correspondingly, the discrepancy of the effective Hamiltonian values is $5.9\times 10^{-2}$ (from $-3.106$ to $-3.165$). The change amount of effective Hamiltonian is enough to determine the thermodynamic stability, as well as the phase boundary in the phase diagram, of self-assembled structures in block copolymers [7, 10, 11, 15]. The results are consistent with the previous theoretical results [7]. However, it should be noted that, in [7], the authors obtained the accurate effective Hamiltonian values through manually changing the scale factor $\Gamma$ for each time. Here, $\Gamma$ is optimized adaptively.

5.2 Patterns on several curved surfaces

In this subsection, we apply our numerical method to several different curved surfaces, including sphere, torus, double torus, heart, orthocircle and parabolic surfaces. Since the SCFT is a high nonlinear system and has multi-solutions, initial values play an important role in determining the final morphologies and accelerating convergent speed. Starting from random initial conditions usually does not obtain energetically favorable patterns. In order to obtain these stable solutions, preferred structures are usually chosen as initial values in the simulations.

5.2.1 Sphere

We firstly present the SCFT results on the sphere surface. To capture these self-assembled structures precisely, we use enough discretization nodes in these simulations. The final equilibrium ordered structures are represented in Fig. 7. The initial scale factor $\Gamma_{0}$ and optimal $\Gamma_{opt}$ for each self-assembled pattern are given in Tab. 3.

For the asymmetric composition of $AB$ diblock copolymers with $f=0.20$, spotted patterns are energetically favorable. To ensure the accuracy, $10242$ nodes and $20480$ elements have been used to calculate these asymmetric cases. In contrast to the case in 2D flat space where only hexagonal pattern is globally stable, various spotted patterns appear on the sphere surface. As shown in Fig. 7 (a), a 12-microdomain spotted phase whose spots locate at vertices of a regular icosahedron. The initial fields of this structure have been given in the above subsection. Then we enlarge the scale factor $\Gamma$ to $11.65$ and use random initial fields. An ordered pattern, with 116-microdomain spots, as shown in Fig. 7 (b), appears in the simulation. It has hexagonal lattice structure along with a small number of pentagonal patterns. Besides Fig. 7 (a) and (b), many spotted phases with different number of microdomains have been also discovered in our simulations. Tab. 4 gives the optimal scale factor of sphere and its corresponding number of microdomains of a spotted phase when $\chi N=25.0$ and $f=0.20$. From this table, one can find that, for spotted patterns, the number of microdomains of a spotted structure is linearly dependent on the sphere surface area.

For the symmetric diblock copolymers, i.e., $f=0.50$, a two-component alternating flat lamellar phase is energetically favorable. On the sphere surface, however, various striped patterns can be formed dependent on the initial fields and the sphere size. As shown in Fig. 7 (c-e), the standard ring and spiral patterns are formed. In these cases, $2560$ nodes and $5120$ elements have been used to calculate these symmetric cases. Among these patterns, Fig. 7 (c) is composed of ordinary ring-form ribbons with point defects at only two opposite site on the sphere surface. The pattern in Fig. 7 (d) is a single spiral ribbon crawling on the sphere surface with two defects. The initial fields for the classical striped phase Fig. 7 (c) and a single spiral ribbon Fig. 7 (d) is $w_{-}=\chi N\sin(8\theta+8\phi),w_{+}=0$, where $\theta=\arctan(y/x)$, $\phi=\arctan(z/y)\in[0,2\pi]$, but with different initial scale factor $\Gamma_{0}$ given in Tab. 3. It is obvious that the scale factor dominates the morphologies of striped phases. Fig. 7 (e) is a semi-ring striped phase whose initial fields is chosen randomly. More abundant spiral patterns on sphere surface have been found in [9] with a large range of $\Gamma$. The stability of striped patterns when $\chi N=12.5$ and $f=0.50$ has been presented by Chantawansri et al. [7] when $\Gamma=3.1\sim 4.9$. However, due to the complexity of these striped patterns, more systematical studies are still absent.

5.2.2 Torus

We next show the SCFT results on the torus surface. For every simulation, we fix the rate of the major and minor radius of torus surface $\gamma=R/r$ and set the initial scale factor $\Gamma_{0}$ as 1. The initial $R_{0}$, $r_{0}$ and optimal scale factor $\Gamma_{opt}$ for each ordered pattern are presented in Tab. 5. The converged self-assembled structures are given in Fig. 8.

For asymmetric $AB$ diblock copolymers, the spotted patterns appear, as shown in Fig 8 (a)-(c) where $\chi N=25.0$ and $f=0.20$. 51200 nodes and 102400 elements are used in these calculations. Using random initial conditions, the spotted phases with different radii can be found as shown in Figs. 8 (a) and (b). To test the influence of initial condition for spotted phases, we use two ways to choose initial condition. For Fig. 8 (c), the first initial condition is the deterministic configuration $w_{-}=0.5(\chi N\cos(3\theta)\cos(3\phi)+\cos(6\phi)),w_{+}=0$, and the second is also the random initial condition. The same spotted patterns surprisingly appear after calculating. It indicates that the morphologies of spotted phases of asymmetric diblock copolymers on torus may be mainly dominated by the geometry of the surface, i.e., radii $R$ and $r$ of torus, and $f$, $\chi N$.

For the symmetric diblock copolymer systems with $f=0.50$, the equilibrium striped phases can be obtained, as shown in Fig. 8 (d)-(g) where $\chi N=16.0$. In these calculations, $64000$ nodes and $128000$ elements are used. From simulations, we find that striped phases are sensitive to initial values which is different from spotted ones. In order to obtain different striped structures, the initial configuration is chosen as

[TABLE]

with different $n_{\phi}$ and $n_{\theta}$, where $\phi,\theta\in[0,2\pi]$ are the parameter coordinates of torus surface. If we choose the initial condition with $n_{\phi}=0$, but with different $n_{\theta}$, the unconnected stripes are obtained. For example, when $n_{\theta}=8$, the stripes of Fig. 8 (d) emerge. If the initial condition is chosen such that $n_{\phi}\neq 0$, we can find equilibrium connected striped patterns with periods $n_{\theta}=n_{\phi}=6,7,8$, as shown in the Fig. 8 (e)-(g). The results are consistent with previous calculations [23].

5.2.3 Three general closed surfaces

To further demonstrate the power of our proposed method, we next show the results of the SCFT simulations on the double torus, heart, and orthocircle surfaces in turn. The parameter setting is $\chi N=25.0$, $f=0.20$ for asymmetric diblock copolymer system, and $\chi N=15.0$, $f=0.50$ for symmetric case $f=0.50$. The number of nodes (elements) used on double torus, heart and orthocircle surfaces is $20642$ ($41288$), $172482$ ($344960$), and $21964$ ($43952$), respectively.

As can be seen in Fig. 9, we have obtained spotted and striped phases of diblock copolymers on these closed surfaces. The initial $\Gamma_{0}$ and optimal scale factor $\Gamma_{opt}$ of the three closed surfaces are given in Tab. 6. The expressions of initial surface $\mathcal{S}_{0}$ have been given by Eqns. (2)-(4). Obviously, the scale factors of these three surfaces have been successfully optimized dependent on the specific ordered structures. For spotted structures (Fig. 9(a)-(c)), the initial condition of the field functions is chosen randomly. For symmetric systems, the striped phases (Fig. 9(d)-(f)) are sensitive to the initial values. On double torus and orthocircle surfaces, the morphologies of Fig. 9 (d) and (f) are obtained by setting the initial value of $w_{+}$ and $w_{-}$ on each mesh node with coordinates $(x,y,z)$ as

[TABLE]

For striped phase on heart surface as shown in Fig. 9 (e), we use the initial value as follows

[TABLE]

where $\theta=\arctan(y/x)$, $\phi=\arctan(z/y)$, and $\theta,\phi\in[0,2\pi]$.

For the asymmetric diblock copolymers, the authors in [8] have found that the elongated spots of a spotted pattern usually locate at the saddle points of a general curved surface. In our simulations, this phenomenon does not appear. The discrepancy may be attributed to the discretization precision and the interaction strength. In [8], only $1912$ nodes were used in spotted phase simulations. However, our calculation uses, at least, $21560$ nodes (10782 elements) to discretize the surface which is enough to ensure the discretization precision. Furthermore, the interaction parameter used here is $\chi N=25.0$ which is stronger than that in [8] where $\chi N=13.0$. The large interaction parameter means strong segregation which contributes to the microscopic phase separation.

5.2.4 Parabolic surface

The above numerical experiments are all applied on closed surfaces. In this subsection we apply the surface finite element method to the SCFT model on a unclosed surface. In particular, we solve SCFT equations on the parabolic surface. For the unclosed surface, the boundary conditions should be considered. In the finite element framework, Dirichlet, Neumann, mixed, or other boundary conditions are all easily supported. Here, we use the homogeneous Neumann boundary condition as an example to test our method. Other boundary conditions can be also used in our proposed numerical framework. The initial values of $w_{+}$ and $w_{-}$ on coordinates $(x,y,z)$ are given by $w_{-}=\chi N\sin(kz)$, $w_{+}=0$ both for asymmetric and symmetric systems. $8515$ nodes and $16720$ elements are used in the following simulations. When $k=5$, we obtain spotted ($\chi N=25.0$, $f=0.20$) and striped ($\chi N=15.0$, $f=0.50$) phases, as shown in the Fig. 10. The initial scale factor $\Gamma_{0}$ and optimal $\Gamma_{opt}$ are given in Tab. 7. Results show that our linear surface finite element method can efficiently capture the standard self-assembled structures of block copolymers both in asymmetric and symmetric cases on parabolic surface. Meanwhile, the size of the parabolic surface can be adaptively optimized during the SCFT iteration.

6 Conclusion and Outlook

In this paper, we proposed a linear surface finite element method to solve the SCFT model and study the self-assembly behaviors of block copolymers on general curved surfaces. At the same time, the surface size indeed has been optimized during iteration procedure which can capture the characteristic size of a given self-assembled structure, and give a more accurate value of effective Hamiltonian. To demonstrate the ability of this algorithm, we applied it to diblock copolymer systems on several distinct curved surfaces, including 5 closed surfaces and an unclosed surface. Numerical results illustrate that our method can successfully obtain the ordered structures of block copolymers on these general curved surfaces.

Block copolymers provide a perfect platform to study the self-assembly behaviors in related physical, chemical and biology systems. To study the phase behaviors of block copolymers precisely, it still requires a high accurate numerical method. To ensure enough precision, we have to use many nodes and elements to describe the ordered structures which costs much CPU time. In the future, we will further improve this computational approach in many directions. The high order finite element method, non-uniform mesh, and adaptive technique will be chosen to improve the precision and reduce computational cost. Another interesting problem is to extend this computational framework to more real-life polymeric systems, for example rigid polymers.

Acknowledgements

This work is supported by National Science Foundation of China (11771368, 11871413, 91430213, 91530321), and Project of Scientific Research Fund of Hunan Provincial Science and Technology Department (2018WK4006). H. Wei is partially supported by Hunan Provincial Civil-Military Integration Industrial Development Project. K. Jiang is partially supported by the research grand from Hunan Science Foundation of China (2018JJ2376), and Youth Project Hunan Provincial Education Department of China (Grant No. 16B257).