acta mechanica solida sinica 30 (2017) 630–637

Available online at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/CAMSS

The application of nonlocal theory method in the

coarse-grained molecular dynamics simulations of
long-chain polylactic acidR

Xiongjun Li a,b, Tan Xiao a,∗, Neng Xiao c

a Architecture and Civil Engineering Institute, Guangdong University of Petrochemical Technology, Maoming 525000,
China
b School of Engineering, Sun Yat-sen University, Guangzhou 510006, China
c School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907, USA

a r t i c l e i n f o a b s t r a c t

Article history: The micro-capsules used for drug delivery are fabricated using polylactic acid (PLA), which
Received 28 March 2017 is a biomedical material approved by the FDA. A coarse-grained model of long-chain PLA
Revised 10 October 2017 was built, and molecular dynamics (MD) simulations of the model were performed using a
Accepted 13 October 2017 MARTINI force field. Based on the nonlocal theory, the formula for the initial elastic modulus
Available online 25 October 2017 of polymers considering the nonlocal effect was derived, and the scaling law of internal
characteristic length of polymers was proposed, which was used to adjust the cut-off radius
Keywords: in the MD simulations of PLA. The results show that the elastic modulus should be computed
Polylactic acid using nonlinear regression. The nonlocal effect has a certain influence on the simulation
Molecular dynamics simulation results of PLA. According to the scaling law, the cut-off radius was determined and applied
Nonlocal theory to the MD simulations, the results of which reflect the influence of the molecular weight
Scaling law change on the elastic moduli of PLA, and are in agreement with the experimental outcome.
Cut-off radius © 2017 Published by Elsevier Ltd on behalf of Chinese Society of Theoretical and Applied
Mechanics.

wall thickness, significantly increasing the ratio of the wall

1. Introduction thickness to the capsule diameter. Due to the differences in
hydrophilicity, the structure and characteristics of the micro-
Used for amplifying signals and having strong biocompat- capsules are affected when the carried drug is being dis-
ibility, conventional ultrasound contrast agents have basic tributed along the inner aqueous layer or within the capsule
structures of hollow bubbles or micro-capsules, and can be wall [1]. Driven by the factors mentioned, the structural com-
modified into drug carriers [1, 2]. In consideration of the cap- position of the ultrasonic contrast agent developed from the
sule permeability to tissues such as capillary walls, the cap- traditional gas bubble or thin walled capsule into a thick-
sule diameter transcends from microscale to nanoscale, while walled multi-layer capsule, whose functions also evolved from
the functionalization of the capsule surface leads to greater basic ultrasonic signal amplification to the targeted release of

R
Project supported by the National Natural Science Foundation of China (no. 11272360), the Natural Science Foundation of Guangdong
Province (no. 2014A030313793), the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the
second phase), and National Supercomputer Center in Guangzhou.
∗
Corresponding author.
E-mail address: [email protected] (T. Xiao).

https://doi.org/10.1016/j.camss.2017.10.003
0894-9166/© 2017 Published by Elsevier Ltd on behalf of Chinese Society of Theoretical and Applied Mechanics.
 acta mechanica solida sinica 30 (2017) 630–637 631

carrier drugs. Therefore, the complexity has increased con- system. A standard MARTINI force field would include four
siderably in terms of the structural design and fabrication of types of basic particles: polar (P), nonpolar (N), apolar (C), and
the carrier capsules; wherein if the potential designs are fab- charged (Q). These four types of particles are further divided
ricated one by one for testing, there would be an impractical into eighteen basic particles of P1–P5, C1–C5, Na, Nd, Nda, N0,
amount of time and fiscal investment required to select a fi- and so on. Normally all the diameters of the particles are fixed
nal design. In comparison, computational simulation meth- at 0.47 nm. The movement of particles is controlled by the
ods have the advantages of low cost, high computation speed, strength of interaction of 10 energy levels.
and the ease of manipulating structure and material param- In molecular simulations, the nonlocal effect caused by
eters, which makes it a widely adapted means in biological long-range forces should be taken into consideration. By ap-
membrane research [3]. Computational molecular simulation plying the cut-off radius, the computation time can be re-
is capable of being introduced as a tool for preliminary re- duced on the premise of certain computational accuracy be-
search in the process of identifying a feasible material and ing assured for simulations of rather large molecular systems.
its available structures for the micro-capsule design. It is also Based on nonlocal theory, an attempt has been made to deter-
worth noticing that the elastic moduli of micro-capsule ma- mine the value of cut-off radius required for simulations of the
terials usually need to be obtained after comparing the data same polymer at different molecular weights, and a formula
from the theoretical model with those from the multi-body of elastic modulus was derived for the nonlinear analysis of
experiments [4, 5], because they are difficult to be measured simulation results. The nonlocal theory [17], which was de-
directly through indentation tests on an individual capsule. veloped by Eringen from the early nonlocal integration theory
Fortunately, the elastic moduli as well as phase transition in 1980s, is a set of stress gradient methods simply containing
temperatures of micro-capsule materials can also be evalu- one scaling parameter. The essence of this theory is that at the
ated using computational simulations. microscale, the stress at a point in material will be affected by
To fabricate the micro-capsules, the commonly selected the stress state of other neighboring points. In 1990s, Aifan-
material was initially albumin, but was later replaced by li- tis rewrote the scaling parameter in the linear elastic form
posome which is more bio-compatible. However, the limited of Laplacian operator with a strain gradient [18], and then in
stability of liposome makes the capsules prone to bursting, so 2010 Askes and Gitman exhibited the possibility of combining
an increasing number of bio-compatible polymers are applied the two above-mentioned methods [19]. The historical devel-
to the fabrication of micro-capsules. Among these polymers, opment of these methods were reviewed in reference [20]. As
polylactic acid (PLA) is non-toxic, degradable with products a continuum theory concerning scaling effects, the nonlocal
that have no adverse effects on the human body, and is a fre- theory has been widely used to study problems such as wave
quently used material for drug carriers and in tissue engineer- propagation, dislocation and crack singularity. It has been in-
ing [6]. It is possible to modify PLA by attaching PEG and chain troduced into micro and nano structural mechanics in 2003
segments of the targeting molecule, or other biological macro- and used to solve the problems of beams [21], bars [22], plates
molecules in various proportions to it, which would adjust the [23], cylindrical shells (for example, carbon nanotubes) [24],
physical properties of the material [7, 8] for the fabrication of and spherical shells [25] at microscale.
medical scaffold or drug carriers. There exists considerable In this work, two macromolecular systems of PLA at the
experimental research on the physical properties of PLA [9], molecular weights of 50 k and 100 k were first built using the
and concurrently, some full-atomic molecular dynamics (MD) coarse-grained method. The PLA system was almost fully an-
simulations have been performed for a qualitative analysis on nealed and equilibrated through cooling simulation and then
the physical properties of small molecule PLA-PEG [10]. Since the curves of self-diffusion coefficient and density of material
polymers with greater molecular weight tend to exhibit higher versus temperature were acquired. After that, the two PLA sys-
strength, greater elastic modulus and better stability, the poly- tems at different molecular weights were separately stretched
mer materials of larger molecular weight are often adopted by MD simulations, and the resulting stress–strain curves were
for the design and fabrication of the carrier capsules. How- acquired and linearly regressed. Meanwhile, the scaling law
ever, the computational cost is too immense when full-atomic for describing the change of the internal characteristic length
MD simulation is performed on polymers with large molecular and the formula for computing the polymer elastic modulus
weight, thus a coarse-grained MD simulation method should were derived based on the nonlocal theory. Finally, the cut-off
be used to simplify the computation. radius was adjusted according to the change of internal char-
The MARTINI force field [11], proposed in 2004 by Marrink acteristic length and the simulations were rerun. The stress–
et al., is a coarse-grained force field initially used to study the strain curves were nonlinearly fitted and compared with the
system behavior of surface-active agents. Given the advantage experimental results.
of having open source modules, many researchers were able to
use and gradually improve upon this force field based on their
specific research topics, which greatly expanded the applica- 2. Force field parameterization and
tion of the force field to biomolecules such as proteins and car- pretreatment
bohydrates [12–15], and furthermore to polyester, polyolefin
and other artificially synthesized polymers [16]. The funda- An isotactic PLA molecule of 10 repetitive units was built in
mental concept of the MARTINI force field is to treat four or Materials Studio (MS). The single PLA molecule was simu-
more non-hydrogen atoms as a coarse-grained particle, and lated through molecular mechanics, and then the centroid
then redefine the non-bonding parameters for the interac- distance between two PLA units, 0.434 nm, was adopted as
tion between the particles to simplify the calculation of the the bond length of coarse-grained PLA particles. The elastic
632 acta mechanica solida sinica 30 (2017) 630–637

Fig. 2 – The change of self-diffusion coefficient of PLA

(Mw = 50 k) with temperature.
Fig. 1 – Schematic diagram of repetitive units in a
coarse-grained chain of PLA molecule. The grey, red, white
balls denote carbon, oxygen, hydrogen atoms, respectively.
(For interpretation of the references to colour in this figure
legend, the reader is referred to the web version of this and to attain ideal lower-state conformations at 300 K. The de-
article.) tailed process of the simulation is: all the coarse-grained MD
simulations were performed in the NPT ensemble using the
Berendsen method to control temperature and pressure. The
equilibrium simulation was performed for 80 ns at tempera-
constants of the bond can be calculated from the energy ture 700 K and environmental pressure 1 bar, with a time step
change with the bond length and angle. The tensile stiffness of 8 fs. Then the system was cooled down to 300 K at a speed of
Kbond is 8000 kJ/mol nm2 , consistent with the elastic constant 5 K/ns. During this period the density change of this PLA sys-
of the bond of Na particle reported in reference [16]. The bend- tem was recoded and the self-diffusion coefficient was calcu-
ing stiffness Kangle is 118 kJ/mol rad2 , and the equilibrium an- lated. It can be seen in Fig. 2 that the self-diffusion coefficient
gle formed by three PLA units, as shown in Fig. 1, measured declines with the decrease of temperature, which indicates a
120° in MS. The electric charge constant is not required be- rapid reduction of system fluidity. When temperature falls be-
cause PLA monomer does not exhibit electric property. The low 350 K, the system behaves like solid and the whole cooling
default number of water molecules in a water particle in MAR- process is equivalent to crystallization.
TINI force field is four, and the total molecular weight is 72, It is shown in Fig. 2 that the self-diffusion coefficient of
equivalent to that of a PLA monomer. As shown in Fig. 1, ev- PLA fluctuates at around 400 K, because in this temperature
ery PLA monomer corresponds to a Na particle, and a PLA range the PLA is in a rubber-like state, where chain segments
molecule is equivalent to a coarse-grained chain of Na par- are free to move relative to one another locally. Despite this
ticles. freedom, the whole chains remain sticky, resulting in the un-
The parameters of particle interactions were added into stable relaxation process of PLA conformation in simulations.
MARTINI force field [26, 27]. The input files of topologies and The contingency of PLA conformation has an obvious effect
parameters of PLA system were compiled respectively. Two on the simulation results due to the limited number of PLA
types of PLA systems were built: one simulation box of 21 nm molecules built for simulations. When the simulation tem-
on each edge, containing 27 coarse-grained chains, with each perature is above 440 K, the PLA is in a liquid state; thus the
chain consisting of 695 repetitive units, was built as the PLA whole chain of molecules can move freely and the change of
system at the molecular weight of 50 k (Mw = 50 k); and an- self-diffusion coefficient tends to stabilize.
other simulation box of 40 nm on each edge, containing 27 The glass transition temperature of PLA can also be pre-
coarse-grained chains, with each chain consisting of 1390 dicted by MD simulations based on dilatometry [10, 28]. Equi-
repetitive units, was built as the PLA system at the molec- librium simulations of 400 ns were performed again at 10 K
ular weight of 100 k (Mw = 100 k). The simulations described or 20 K intervals from 500 K to 280 K. As shown in Fig. 3, two
in this paper were performed with the GROMACS simulation straight lines representing the change of PLA density with
package. temperature were determined by the average values of sys-
These systems were optimized through molecular me- tem density calculated from the simulation results at differ-
chanics simulations to eliminate the force between particles ent temperatures. The glass transition temperature of PLA,
caused by improper conformations at the preliminary stage. 328.5 K, was determined from the point of intersection of
After that, the systems underwent annealing treatment to en- these two lines, which is in agreement with 328.8 K, the ex-
sure sufficient migration of macromolecules in the system, perimental result in reference [10].
 acta mechanica solida sinica 30 (2017) 630–637 633

Fig. 3 – The change of PLA (Mw = 50 k) density with Fig. 5 – Stress–strain curves of PLA (Mw = 100 k) in tensile
temperature and the glass transition temperature. simulations when cut-off radius r = 1.2 nm.

simulation process, and then the stress–strain curves were de-

picted in Fig. 4. Similarly, three PLA systems at the molecular
weight of 100 k were built and simulated when 1.2 nm was also
adopted for the cut-off radius r. The stress–strain curves were
depicted in Fig. 5.
Typically, the elastic modulus was computed by the lin-
ear regression of the stress–strain data acquired from simu-
lations [16]. For example, three curves of groups A, B, and C in
Fig. 4 were linearly fitted and their elastic moduli are 0.55 GPa,
0.57 GPa, and 0.50 GPa, respectively. As control subjects, linear
regression of the curves of groups A, B, and C in Fig. 5 gave the
elastic moduli of 0.52 GPa, 0.50 GPa, and 0.52 GPa, respectively.
However, the experimental result reported in reference [9] are:
1.2 GPa for the elastic modulus of l-PLA (Mw = 50 k), 2.7 GPa for
l-PLA (Mw = 100 k), and 1.9 GPa for d,l-PLA (Mw = 107 k). Here
two problems arise: the elastic moduli given by linear regres-
Fig. 4 – Stress–strain curves of PLA (Mw = 50 k) in tensile
sion are significantly different from the experimental results
simulations when cut-off radius r = 1.2 nm.
in reference [9]; and the simulation results of elastic moduli
do not show the trend of increase when the molecular weight
of PLA is increased from 50 k to100 k. Therefore, an attempt is
made to introduce the nonlocal theory and use the nonlinear
3. Tensile simulations and linear regression form of stress–strain relationship to analyze the simulation
of simulation results results.

Tensile simulations were performed by nonequilibrium kinet-

ics, for understanding the change of stress when the PLA sys- 4. Derivation of elastic modulus based on
tem deforms along a certain single direction and computa- nonlocal theory
tion of elastic modulus. Five PLA systems at the molecular
weight of 50 k were built and simulated to examine the re- The basic formula of nonlocal stress for isotropic elasticity is
peatability of MD simulations, as shown in Fig. 4, and to meet [25]:
the data requirement for the computation of elastic moduli

based on the nonlinear regressions at the late stage. At first,  
T (x ) = α(x − x, γ )σ (x )d(x ) (1)
the PLA systems which have been fully equilibrated in the pre- 
vious section underwent MD equilibrium simulation again for
24 ns. The tensile speed (strain rate) was set at 10−4 /ns and the
σ =C:ε (2)
cut-off radius r was 1.2 nm. The simulation box was gradually
elongated in the z direction until a 2% strain was achieved. The where T(x) and σ (x) are respectively the nonlocal stress and
average stress (the average value of pressure Pzz at both sides the classical local stress at a point x in a body, C is the elastic
of the box in the z direction in the corresponding time domain) coefficient, and ε the local strain. The function α(|x − x|, γ ) de-
and strain data were recorded at intervals of 8 ns during the scribes the nonlocal modulus, in which |x − x| is the distance
634 acta mechanica solida sinica 30 (2017) 630–637

Table 1 – Simulation results of PLA (Mw = 50 k) when

r = 1.2 nm.

Group C1 C3 Groups κ E
−5
A 1.28 1556 A+B 4.96 × 10 1.20
B 1.12 −1661 A+C 4.67 × 10−5 1.21
C 1.21 44 A+D 1.15 × 10−4 1.10
D 1.39 2567 A+E 7.56 × 10−5 1.16
E 1.14 −264 B+C 5.22 × 10−5 1.21
The average value of E is B+D 6.53 × 10−5 1.23
1.18 GPa with the standard B+E 1.57 × 10−5 1.14
deviation of 0.04 C+D 7.42 × 10−5 1.20
Fig. 6 – Micro unit cell is stretched along the z-axis to the C+E 2.18 × 10−4 1.20
range of cut-off radius r. D+E 8.98 × 10−5 1.16

Note: The units of C1 , C3 , and E are GPa, and κ is a dimensionless

constant.
that point x moves before and after deformation, and γ the
material constant defined by γ = e0 a/l. Here e0 is a material-
dependent constant, while a and l are the internal character- deformation of ε = ε max , namely, zmax = (1 + ε max )z0 = r, then
istic length (e.g., lattice scale, effective molecular radius, etc.) becomes the variable z = (1 + ε)z0 = (1 + ε) r/(1 + εmax ). Sub-
and external characteristic length (e.g., crack length, etc.), re- stituting it into Eq. (4), the elastic modulus at ε → 0 is derived,
spectively.  
d  e a 2 2
0 2d T
Chains interpenetrate and entangle copiously in an amor- E = lim T− (1 + εmax ) (5)
ε→0 dε r dε2
phous polymer. The segment between two adjacent entangle-
ment points on a chain is basically considered ideal and can where T is the stress Pzz at the cut-off radius (i.e., z = r) .
move freely [29, 30]. The elastic properties of polymer chains Under comparatively small deformation conditions, the
primarily depend on the distribution of entanglement points nonlocal stress can be approximated below after the high-
[30], and the entanglement behavior is mainly controlled by order small quantities are neglected
chain tortuosity [31]. Meanwhile, the internal forces in the
polymer transfer primarily through the entanglement points T = C0 + C1 ε + C2 ε 2 + C3 ε 3 (6)
between chains. The longer the chain is, the more the entan-
glement points there are, and the higher the efficiency of force Substituting Eq. (6) into Eq. (5), the initial elastic modulus
transference is. Graessley and Edwards [32] noted that melt is obtained
elasticity must involve two independent length scales: Kuhn
length and contour length per unit volume. Therefore, based E = C1 − κC3 (7)
on the nonlocal theory, the internal characteristic length is de-
e a
fined here in proportion to the root mean square of end-to-end
√ where the coefficient κ = 6 ( 0r )2 , reflecting the influence of
distance of ideal polymer chains, a ∝ < R2 >. nonlocal effect on the elastic modulus of polymer material.
Then, Eq. (1) can be transformed into an equivalent differ-
ential form
  5. Nonlinear regression of simulation results
2
1 − ( e0 a ) ∇ 2 T = σ (3) and discussions

When the unit cell system is subjected to uniaxial tension There are five sets of data, from group A to E, in Fig. 4. Each set
in the z direction, as shown in Fig. 6, Eq. (3) can be simplified was nonlinearly fitted by a cubic polynomial of Eq. (6) to calcu-
after substituting Hooke’s law σ = Eε: late coefficients C1 and C3 . Five sets of C1 and C3 were acquired,
  from which each two sets were selected to form a system of
2 d2 equations as per Eq. (7), for the solution of elastic modulus E
1 − ( e0 a ) T = Eε (4)
d z2 and coefficient κ. Ten results were obtained and tabulated in
Table 1. The average value of the elastic modulus of PLA at
The physical essence of nonlocal theory is taking the in- molecular weight of 50 k is 1.18 GPa with standard deviation
fluence from the neighboring material into account, which is of 0.04, which matches well with the experimental result of
in correspondence with the concept of cut-off radius in sim- 1.20 GPa reported in reference [9]. It shows a good reproducibil-
ulations. For a simulated system, the internal characteristic ity of the five simulations. The initial elastic modulus of PLA
length should rest within the range of cut-off radius r adopted is mainly determined by the coefficient C1 and modified by C3 ,
for simulation computation, beyond which the noncovalent which describes the influence of nonlocal effect.
forces acting on molecules will be isolated and the nonlo- Similarly, the three sets of data in Fig. 5 were nonlinearly
cal effect no longer exists. As shown in Fig. 6, a micro unit fitted and calculated. The results were tabulated in Table 2,
cell of length z0 is stretched to the cut-off radius r after the which are different from the experimental results (with the
 acta mechanica solida sinica 30 (2017) 630–637 635

Table 2 – Simulation results of PLA (Mw = 100 k) when Table 3 – Simulation results of PLA (Mw = 100 k) when
r = 1.2 nm. r = 1.7 nm.

Group C1 C3 Groups κ E
Group C1 C3 Groups κ E
−5
A 1.67 2363 A+B 9.76 × 10 1.44
B 1.35 −981 B+C 5.79 × 10−4 1.92 A 2.28 1319 A+B −5
8.35 × 10 2.17
C 0.14 −3061 C+A 2.82 × 10−4 1.01 B 2.05 −1417 A+C 2.64 × 10−4 1.93
C 2.08 553 A+D 5.16 × 10−4 1.60
Note: The units of C1 , C3 , and E are GPa, and κ is a dimensionless D 3.01 2746 A+E 1.64 × 10−4 2.06
constant. E 2.05 −54 B+C 1.33 × 10−4 2.07
The average value of E is B+D 2.32 × 10−4 2.38
2.02 GPa with the standard B+E 2.49 × 10−6 2.05
deviation of 0.21 C+D 4.28 × 10−4 1.84
C+E 3.76 × 10−5 2.05
D+E 3.43 × 10−4 2.07

Note: The units of C1 , C3 , and E are GPa, and κ is a dimensionless

constant.

Table 4 – Simulation results of PLA (Mw = 100 k) when

r = 2.4 nm.

Group C1 C3 Groups κ E
−4
A 2.158 −1464 A+B 5.93 × 10 3.03
B 1.639 −2339 B+C 2.12 × 10−4 2.13
C 2.647 2421 C+A 1.26 × 10−4 2.34

Note: The units of C1 , C3 , and E are GPa, and κ is a dimensionless

constant.
Fig. 7 – Stress–strain curves of PLA (MW = 100 k) in tensile
simulations when cutoff radius r = 1.7 nm.

proportion, hence the values of elastic moduli in Table 3 are

close to the lower limit of the experimental results (with the
elastic moduli of l-PLA (Mw = 100 k) being 2.7 GPa, and d,l- elastic modulus of d,l-PLA (Mw = 107 k) being 1.9 GPa) [9].
PLA (Mw = 107 k) being 1.9 GPa) in reference [9]. This difference For a more complete understanding on the influence of
could be attributed to the increase of internal characteristic the value selection of cut-off radius on the PLA simulation re-
length a of PLA when the molecular weight of PLA grows from sults, the cut-off radius r was further increased from 1.7 nm
50 k to 100 k. The simulations with unchanged cut-off radius to 2.4 nm, and then three PLA systems at a molecular weight
of r = 1.2 nm cannot reflect the change in the mechanical prop- of 100 k underwent simulations again. The simulation results
erties of PLA with its increase of molecular weight. were tabulated in Table 4. The comparison of the data in Tables
In the section above, the internal characteristic length of 2–4 shows a general trend that the elastic moduli increase
√
polymer a ∝ < R2 >. According to Flory’s theory, the mean with the growth of cut-off radius, but at a diminishing rate.
square end-to-end distance of ideal polymer chains is propor- Although sufficient long PLA molecules were selected and en-
tional to the degree of polymerization, < R2 > ∝N. Thus, the in- tanglement behavior was mimicked in the simulations, the
√
ternal characteristic length should be 2 times of the original elastic modulus of group A in Table 4 is much greater than the
value when the molecular weight of PLA doubles. Therefore, other two. The distribution of molecular conformation of the
the cut-off radius r increases from 1.2 nm to 1.7 nm for the PLA system was checked but no sign of predominant proportion
system at the molecular weight of 100 k before the tensile sim- of l-PLA was found. One possible explanation might be that
ulations resume. Five PLA systems at the molecular weight of the PLA system consisting of 27 molecular chains, whose size
100 k were built and simulated, and the stress–strain curves is limited by the computation time, is not sizable enough to
were depicted in Fig. 7. suppress the influence of conformation contingency of some
The five sets of data in Fig. 7 were nonlinearly fitted and cal- molecules on the statistical regularity of macroscopic proper-
culated, and the results were tabulated in Table 3. The average ties of materials. In short, using the scaling law to determine
value of the elastic modulus of PLA at the molecular weight the cut-off radius in molecular simulations seems to be a typ-
of 100 k, after the proportional adjustment of cut-off radius, is ical solution for the trade-off between computation time and
2.02 GPa with standard deviation of 0.21, which falls within the accuracy.
scope of the experimental results of 1.9–2.7 GPa [9]. The PLA By comparing Table 1 with Table 3, it seems that the elastic
molecular system was built in the form of amorphous chains moduli of PLA increase with the growth of molecular weight.
before the simulation commenced. When the simulation was The increase of molecular weight has strengthened the entan-
completed, it was found that two molecular conformations glement between polymer molecules, and therefore the con-
coexist in the system, in which d,l-PLA takes predominant nection of molecules in all directions. The chain entanglement
636 acta mechanica solida sinica 30 (2017) 630–637

√
a ∝ < R2 >, was proposed based on the concept that force
Table 5 – Elastic moduli E versus self-diffusion coeffi-
cients D at temperature 300 K. transfers with entanglement points, and a method for esti-
mating the value of cut-off radius for molecular simulations
PLA (Mw = 50 k, r = 1.2 nm) PLA (Mw = 100 k, r = 1.7 nm) is given accordingly. The research results suggest that: (1) the
Group E D Group E D elastic modulus should be acquired through nonlinear regres-
sion. (2) The nonlocal effect has a certain influence on the
A 1.20 1.48 A 2.17 1.9
simulation results of PLA. (3) According to the scaling law, the
B 1.21 1.47 B 2.07 2.01
C 1.21 1.47 C 1.93 1.98 cut-off radius is determined and applied to the MD sim-
ulations, of which the results reflect the influence of the
Note: The unit of E is GPa, and the unit of D is × 10 − 11 m2 /s. molecular weight change on the elastic moduli of PLA with
considerable accuracy, and are in good agreement with the ex-
perimental outcome.
impedes the relative slippage between polymer molecules
during the process of tensile simulations and raises the value
of elastic modulus. Conversely, the decrease of chain entan- references
glement with the decrease of molecular weight weakens the
intermolecular force, and the elastic modulus declines corre-
spondingly. When the molecular weight decreases to a certain
[1] S. Tinkov, R. Bekeredjian, G. Winter, C. Coester, Microbubbles
critical value, the polymer no longer has tensile strength. as ultrasound triggered drug carriers, J. Pharm. Sci. 98 (2009)
The elastic moduli and related self-diffusion coefficients 1935–1961.
at 300 K of PLA for two different molecular weights were tab- [2] S.P. Qin, C.F. Caskey, K.W. Ferrara, Ultrasound contrast
ulated in Table 5, respectively. It shows that elastic moduli as microbubbles in imaging and therapy: physical principles
well as self-diffusion coefficients have a tendency to increase and engineering, Phys. Med. Biol. 54 (2009) R27–R57.
[3] M. Stepniewski, A. Bunker, M. Pasenkiewicz-Gierula,
when PLA molecular weight grows. Larger self-diffusion co-
M. Karttunen, T. Rog, Effects of the lipid bilayer phase state
efficient implies an increase in mobility and a decrease in
on the water membrane interface, J. Phys. Chem. B 114 (2010)
frictional interaction between PLA chain segments, which is 11784–11792.
unfavorable to load transfer. But as mentioned earlier, the [4] L. Hoff, Oscillations of polymeric microbubbles: effect of the
elastic performance of polymers is in fact predominately de- encapsulating shell, J. Acoust. Soc. Am. 107 (2000) 2272–2280.
termined by the chain entanglement of molecules. PLA with [5] A.N. Damir, B. Khismatullin, Radial oscillations of
larger molecular weight has longer chains and a greater num- encapsulated microbubbles in viscoelastic liquids, Phys.
Fluids 14 (2002) 3534–3557.
ber of entanglement points, and thus greater elastic modulus.
[6] D. Garlotta, A literature review of poly(lactic acid), J. Polym.
The fluctuation of curves in Fig. 7 is quite different from Environ. 9 (2001) 63–84.
that in Fig. 4 because the intermolecular entanglement and [7] A.P. Mathew, K. Oksman, M. Sain, Mechanical properties of
the conformation of molecules are more stable with the in- biodegradable composites from poly lactic acid (PLA) and
creases of molecular weight and chain length, leading to microcrystalline cellulose (MCC), J. Appl. Polym. Sci. 97 (2005)
lower probability of slippage, plastic deformation, and sudden 2014–2025.
change of polymer conformation. In actual polymer systems, [8] C.C. Chen, J.Y. Chueh, H. Tseng, H.M. Huang, S.Y. Lee,
Preparation and characterization of biodegradable PLA
the degree of crystallization increases with larger molecular
polymeric blends, Biomaterials 24 (2003) 1167–1173.
weight and better stability of conformation, so more elasticity [9] I. Engelberg, J. Kohn, Physico-mechanical properties of
and brittleness is observed. On the contrary, more plasticity is degradable polymers used in medical applications: a
observed in the polymer system while the degree of crystal- comparative study, Biomaterials 12 (1991) 292–304.
lization decreases. [10] H.Z. Xiang, W.H. Yang, C.B. Hu, Y.H. Mu, J.F. Li, Study of the
mechanical and thermal properties of poly(lactic acid) and
poly(ethylene glycol) block copolymer with molecular
dynamics, J. Polym. Anal. Charact. 15 (2010) 235–244.
6. Conclusions [11] S.J. Marrink, A.H. de Vries, A.E. Mark, Coarse grained model
for semiquantitative lipid simulations, J. Phys. Chem. B 108
MD simulation is widely used in the study of polymer prop- (2004) 750–760.
erties, through which the intrinsic characteristics of certain [12] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de
polymers in its stable conformation can be roughly deter- Vries, The MARTINI force field: coarse grained model for
biomolecular simulations, J. Phys. Chem. B 111 (2007)
mined. In this work, the coarse-grained models of long-chain
7812–7824.
PLA were built and simulated through the MD method, us- [13] L. Monticelli, S.K. Kandasamy, X. Periole, R.G. Larson,
ing the MARTINI force field. The change of self-diffusion D.P. Tieleman, S. Marrink, The MARTINI coarse-grained force
coefficient and density with temperature, and the glass transi- field: extension to proteins, J. Chem. Theory Comput. 4 (2008)
tion temperature were acquired from the simulations. Using 819–834.
a third-order polynomial to approximate the nonlocal stress [14] C.A. Lopez, A.J. Rzepiela, A.H. de Vries, L. Dijkhuizen,
P.H. Hunenberger, S.J. Marrink, Martini coarse-grained force
of polymer under the condition of comparatively small strain,
field: extension to carbohydrates, J. Chem. Theory Comput. 5
the formula for computing the initial elastic modulus consid-
(2009) 3195–3210.
ering the nonlocal effect based on the nonlocal theory was de- [15] J.J. Uusitalo, H.I. Ingolfsson, P. Akhshi, D.P. Tieleman,
rived. By combing polymer physics with the nonlocal theory, S.J. Marrink, Martini coarse-grained force field: extension to
the scaling law of internal characteristic length of polymers, DNA, J. Chem. Theory Comput. 11 (2015) 3932–3945.
 acta mechanica solida sinica 30 (2017) 630–637 637

[16] G. Rossi, I. Giannakopoulos, L. Monticelli, N.K.J. Rostedt, [24] S.A. Fazelzadeh, E. Ghavanloo, Nonlocal anisotropic elastic
S.R. Puisto, C. Lowe, A.C. Taylor, I. Vattulainen, T. Ala-Nissila, shell model for vibrations of single-walled carbon nanotubes
A MARTINI coarse-grained model of a thermoset polyester with arbitrary chirality, Compos. Struct. 94 (2012) 1016–1022.
coating, Macromolecules 44 (2011) 6198–6208. [25] E. Ghavanloo, S.A. Fazelzadeh, Nonlocal elasticity theory for
[17] A.C. Eringen, On differential equations of nonlocal elasticity radial vibration of nanoscale spherical shell, Eur. J. Mech. A
and solutions of screw dislocation and surface waves, J. Appl. Solid 41 (2013) 37–42.
Phys. 54 (1983) 4703–4710. [26] L. Monticelli, S. Kandasamy, X. Periole, R. Larson,
[18] E.C. Aifantis, On the role of gradients in the localization of D.P. Tieleman, S.J. Marrink, The MARTINI coarse grained force
deformation and fracture, Int. J. Eng. Sci. 30 (1992) 1279– field: extension to proteins, J. Chem. Theory Comput. 4 (2008)
1299. 819–834.
[19] H. Askes, I.M. Gitman, Review and critique of the stress [27] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de
gradient elasticity theories of Eringen and Aifantis, Adv. Vries, The MARTINI force field: coarse grained model for
Mech. Math. 21 (2010) 203–210. biomolecular simulations, J. Phys. Chem. B 111 (2007)
[20] H. Askes, E.C. Aifantis, Gradient elasticity in statics and 7812–7824.
dynamics: an overview of formulations, length scale [28] C. Wei, Thermal expansion and diffusion coefficients of
identification procedures, finite element implementations carbon nanotube-polymer composites, Nano Lett. 2 (2002)
and new results, Int. J. Solids Struct. 48 (2011) 1962–1990. 647–650.
[21] J. Fernandez-Saez, R. Zaera, J.A. Loya, J.N. Reddy, Bending of [29] P.-G. de Gennes, Scaling Concepts in Polymer Physics, Cornell
Euler–Bernoulli beams using Eringen’s integral formulation: University Press, New York, 1985.
a paradox resolved, Int. J. Eng. Sci. 99 (2016) 107–116. [30] P.J. Flory, Molecular theory of rubber elasticity, Polym. J. 17
[22] W. Sumelka, R. Zaera, J. Fernandez-Saez, A theoretical (1985) 1–12.
analysis of the free axial vibration of non-local rods with [31] S.H. Wu, Chain structure and Entanglement, J. Polym. Sci.
fractional continuum mechanics, Meccanica 50 (2015) Part B Polym. Phys. 27 (1989) 723–741.
2309–2323. [32] W.W. Graessley, S.F. Edwards, Entanglement interactions in
[23] P. Lu, P.Q. Zhang, H.P. Lee, C.M. Wang, J.N. Reddy, Non-local polymers and the chain contour concentration, Polymer 22
elastic plate theories, Proc. R. Soc. A 463 (2007) 3225–3240. (1981) 1329–1334.