On phase diagram and the pseudogap state in a linear chiral homopolymer model
A. Sinelnikova,1, ∗ A.J. Niemi,2, 3, 4, † and M. Ulybyshev5, 1, 6, ‡
1
arXiv:1501.05267v1 [cond-mat.soft] 21 Jan 2015
Institute for Theoretical Problems of Microphysics,
Moscow State University, Moscow, 119899 Russia
2
Department of Physics and Astronomy, Uppsala University, P.O. Box 803, S-75108, Uppsala, Sweden
3
Laboratoire de Mathematiques et Physique Theorique CNRS UMR 6083,
F´ed´eration Denis Poisson, Universit´e de Tours, Parc de Grandmont, F37200, Tours, France
4
Department of Physics, Beijing Institute of Technology, Haidian District, Beijing 100081, P. R. China
5
Institute of Theoretical Physics, University of Regensburg,
D-93053 Germany, Regensburg, Universitatsstrasse 31
6
ITEP, B. Cheremushkinskaya str. 25, Moscow, 117218 Russia
The phase structure of a homopolymer chain is investigated in terms of a universal theoretical
model, designed to describe the infrared limit of slow spatial variations. The effects of chirality are
studied and compared with the influence of a short-range attractive interaction between monomers,
at various ambient temperature values. In the high temperature limit the homopolymer chain is in
the self-avoiding random walk phase. But at low temperatures, two different phases are possible:
When short-range attractive interactions dominate over chirality, the chain collapses into a spacefilling conformation. But when the attractive interactions become weaker, there is a low temperature
unfolding transition and the chain becomes like a straight rod. Between the high temperature and
low temperature limits, several intermediate states are observed. For sufficiently high values of
short-range attraction , the conventional θ-regime is observed between the self-avoiding random walk
phase and the space filling collapsed phase. But when chirality increases, there is a trasition from
the θ-regime to a pseudogap state. Moreover, a regime akin the θ-regime is identified between the
pseudogap state, and the low temperature phase where the chain is like a straight rod. Applications
to polymers and proteins, in particular collagen, are suggested.
I.
INTRODUCTION
A linear homopolymer is made of a single type of a
repeat unit. An important example is polyacetylene, an
organic conductive polymer which is the paradigm material for fractional fermion number [1, 2]. Additional examples, among many others, are poly-L-lysine and polyL-glutamic. The former is a food preservative with potential for wider, even pharmaceutical antimicrobial effects [3] while the latter is used for drug delivery against
cancer [4].
From a theoretical point of view, the concept of a homopolymer chain can provide a useful coarse grained approximation, even in the case of a heteropolymer that
exhibits only approximatively repeating patterns: For a
sufficiently long chain the distinct monomers are simply
combined into appropriate subunits, to dispose of the inhomogeneities in the monomer species. For example collagen, which is the most abundant protein in mammals,
displays a repeated glycine-proline-X pattern, where X
is any amino acid other than glycine and proline. DNA,
RNA and the Cα backbone of a protein chain are additional examples where a homopolymer approximation is
occasionally profitably introduced [5–7].
Here the phase structure of linear homopolymers is in-
∗ Electronic
address: [email protected]
address: [email protected]; URL: http://
www.folding-protein.org
‡ Electronic address: [email protected]
† Electronic
vestigated, in terms of a universal energy function [8–
13]. Of particular interest is the effect of chirality on the
phase diagram. A chiral hopolymer is one where parity
is broken; the mirror image of a stable chiral homopolymer conformation is in general not stable. Chiral polymers often have a tendency to form helical structures.
For example in the case of proteins right-handed helical
structures are more common than left-handed ones. Even
though there are also chiral proteins that can form different right-handed and left-handed structures. An important example is collagen where the right-handed polyproline I conformation is more compact than the left-handed
polyproline II.
A general argument which is due to Kadanoff and Wilson [14–18] states, that in the thermodynamical limit the
phase transition properties of a material system are often universal i.e. independent of the atomic level details.
From this perspective, the construction of the phase diagram of a linear and chiral, structureless homopolymer
should be relevant for the understanding of the phase
diagram of more elaborate linear chiral homopolymers,
maybe even that of certain heteropolymers [5–7].
All linear polymers are presumed to have a very similar phase structure, quite independently of their chemical
composition [5–7]. Even though the phase where a particular polymer resides depends on many factors such as
concentration, the quality of solvent, ambient temperature and pressure. Three different, universal phases are
commonly identified. These phases are categorised by
the way how the polymer structure fills the space [5–7]:
Under poor solvent conditions or at low temperatures,
when the attractive interactions between the monomers
2
dominate, a single polymer chain is presumed to collapse
into a configuration which is space filling. On the other
hand, in a good solvent or at high temperatures, when the
repulsive interactions dominate and cause the chain to effectively swell, its geometric structure resembles that of
a self-avoiding random walk (SARW). Between the two,
there is a θ-regime where the attractive and repulsive interactions cancel each other. In the θ-regime the polymer
chain is presumed to have the characteristics of an ordinary random walk (RW). Finally, some polymers such as
collagen for example, are more like straight, riding rods.
Each of the four phases – rigid rod, SARW, RW and
the space filling one – is characterised by the inverse of
their Hausdorff dimension, called the scaling exponent ν
[5–7]. To define this quantity, one introduces the radius of
gyration Rgyr . With ri the coordinates of the individual
monomers, the radius of gyration is [5–7],
s
1 X
Rgyr =
(ri − rj )2
(1)
2N 2 i,j
In the limit where the number of monomers N becomes
very large, the radius of gyration has the asymptotic expansion [19–22]
N large
R02 N 2ν (1 + R1 N −δ1 + ...) ∼ R02 N 2ν + . . .
(2)
Here the length scale R0 is called the Kuhn length. It
is the effective distance between the monomers in the
large-N limit. The Kuhn length is not a universal quantity, its value can in principle be computed from the
atomic level details of the polymer and environment including pressure, temperature and chemical microstructure of the solvent. The dimensionless scaling exponent
ν i.e. the inverse Hausdorff dimension that governs the
large-N asymptotic form of equation (2), is presumed to
be a universal quantity. Its numerical value is independent of the local atomic level structure of the polymer
[5, 19? –21]. The δ1 etc. are critical exponents and
the R1 etc. are the corresponding amplitudes, and together they constitute the finite-size corrections. The δ1
etc. are universal quantities [21], but the R1 etc. are not
universal [21].
The following mean field values are conventionally assigned to ν [5, 7]:
2
Rgyr
−→

1/3


1/2
ν =

 3/5
1
collapsed
SARW
RW
rod
(3)
Under poor solvent conditions or at low temperatures,
the polymer collapses into a space-filling conformation
[23–25] with the mean field exponent ν = 1/3. Folded
proteins are commonly found in this phase. For an ordinary random-walk (RW) the mean field value is ν = 1/2.
This corresponds to the θ regime, that separates the
collapsed phase from the high-temperature self-avoiding
random walk phase for which the Flory value ν = 3/5
is found. Finally, when ν = 1, the polymer loses its inherently fractal structure and behaves like a straight rod
[25].
The transition between the collapsed phase and the
SARW phase can involve the RW phase as a tri-critical
θ-point. More generally, the RW phase can appear as
a transitional θ-regime, separating the high temperature
SARW and low temperature collapsed phases. The transitions between the rigid rod phase and the other three,
is less studied. However, there is a physically and biologically very important scenario where such a transition
could have a rˆole; that of cold denaturation of a protein chain. The presence of all four phases (3) opens the
possibility of a 4-critical point, under proper conditions
[26].
The θ-point value ν = 1/2 is exact for a polymer with
no long range interactions [5]. For a space-filling structure the value ν = 1/3 is also exact, and similarly ν = 1 is
exact for a straight, linear rod-like structure. But in the
case of SARW the mean field value is corrected by fluctuations. Using the argument that the SARW is in the
same universality class with the O(n) symmetric scalar
field theory with a quartic self-interaction, the following
modified value was obtained [19, 20]
ν = 0.5880 ± 0.0015
A subsequent numerical Monte Carlo evaluation, computed directly by using the self-avoiding random-walk
model on a square lattice, gave a very similar value [21]
ν = 0.5877 ± 0.0006
In the present article, the phase structure of a homopolymer is studied using numerical simulations, in the
context of a universal off-lattice energy function. Different Monte Carlo algorithms are introduced and analysed,
to identify an optimal one for describing a homopolymer. The phase digram that corresponds to the geometric characteristics (3) is then constructed, by inspecting
the radius of gyration in terms of ambient temperature
and the various parameters that appear in the Hamiltonian, in particular the one that characterises the chirality
and the one that characterises the attractive (hydrophobic) forces.
II.
THE MODEL
A.
Geometry
Consider a homopolymer chain with i = 1, ..., N
monomers and let ri be the three dimensional space coordinates of the monomers. The unit tangent vectors along
the lines that connect two consecutive monomers are
ti =
ri+1 − ri
|ri+1 − ri |
(4)
3
τi → τi + ∆i−1 − ∆i
The unit binormal vectors are defined by
bi =
ti−1 × ti
|ti−1 × ti |
(5)
The unit normal vectors are defined by
ni = bi × ti
(6)
The three vectors (ni , bi , ti ) determine an orthonormal
frame at the monomer position ri . The discrete bond
angles are
κi ≡ κi+1,i = arccos (ti+1 · ti )
(7)
and the discrete torsion angles are
(14)
The range of τi is [−π, π) mod(2π). The equations (13)
and (14) may be used to extend the range of the bond
angle from [0, π) to κi into [−π, π) mod(2π). The extension is compensated for by the following discrete Z2
symmetry
κk → − κk
τi →
τi − π
for all k ≥ i
(15)
that leaves the chain intact.
In the numerical simulations presented here, all the
distances between nearest neighbour monomers are fixed
to the uniform constant value
τi ≡ τi+1,i = sgn[(bi−1 ×bi )·ti ]×arccos (bi+1 · bi ) (8)
|ri+1 − ri | = δ = 3.8 (˚
A)
Conversely, when the angles (κi , τi ) are known the discrete Frenet equation [27]

 

 
ni+1
cos κ cos τ cos κ sin τ − sin κ
ni
bi+1  =  − sin τ
bi 
cos τ
0 
ti+1
sin κ cos τ sin κ sin τ cos κ i+1,i ti
(9)
determines the frames iteratively, by computing the
frame at the position of the (i + 1)th monomer from
the frame at the position of the ith monomer. Once all
the frames have been constructed, the entire chain is obtained as follows,
This equals the average distance between two consecutive Cα atoms along a protein backbone, measured in
˚
Angstr¨om’s.
A polymer is subject to steric constraints, due to overlapping electron clouds and various short range Born
repulsions. Accordingly the following forbidden volume
constraint is introduced,
rk =
k−1
X
|ri+1 − ri | · ti
|ri − rk | ≥ δ ≡ 3.8 (˚
A)
for |i − k| ≥ 2
(16)
(17)
This is in line with the minimum distance observed between any two Cα atoms, in folded protein structures.
The numerical values (16) and (17) can both be independently modified, with no effect to conclusions.
(10)
i=0
With no loss of generality one can set r0 = 0, and orient
t0 to point into the direction of the positive z-axis.
A framing is necessary for the construction of the chain
from the bond and torsion angles. But the equation (10)
does not involve the vectors ni and bi . Thus any linear
combination of these two vectors could be chosen to define a framing, to construct the chain from the angles.
For this, consider a local SO(2) transformation that rotates the frame (ni , bi ) by an angle ∆i leaving ti intact,
 
  
 
n
n
cos
∆
sin
∆
0
n
i
i
3
b → e∆i T b = − sin ∆i cos ∆i 0b (11)
t i
t i
0
0
1
t i
where (T a )bc = abc are the SO(3) generators
[T a , T b ] = abc T c
When the vectors ni and bi are combined into a complex
vector, (11) can be recast as follows [27],
ni + ibi → e
i∆i
(ni + ibi ) ≡
e1i
+
ie2i
(12)
The frame rotation (11) introduces transformations in
the equation (9), that can be summarised as follows:
3
κi T 2 → e∆i T (κi T 2 ) e−∆i T
3
(13)
B.
Hamiltonian
The bond and torsion angles constitute a complete set
of geometric variables in the case of protein Cα backbones [28]. Furthermore, according to (2) (3) the structural phase diagram of a homopolymer is determined by
the three dimensional chain geometry. Thus the bond
and torsion angles are a complete set of order parameters,
in the sense of Kadanoff and Wilson: The geometrically
defined, structural phase diagram of a homopolymer can
be fully determined by a thermodynamical free energy
which is constructed from these order parameters only.
Consider a homopolymer chain in thermal equilibrium
and let F be the thermodynamical Helmholtz free energy. Thus, the minimum energy configuration of F describes the chain, under thermodynamical equilibrium
conditions. The free energy is the sum of the internal
energy U and the entropy S, at temperature T
F = U − TS
(18)
It is a function of all the inter-atomic distances
F = F (rαβ ) ;
rαβ = |rα − rβ |
(19)
where the indices α, β, ... extend over all the atoms in
the homopolymer system, including those of the solvent
4
environment. Consider the infrared, long distance limit
where the characteristic length scales of spatial deformations along the homopolymer chain around its thermal
equilibrium configuration are large in comparison to the
distance (16) between neighboring monomers. This is
synonymous to an assumption that there are no abrupt
wrenches and buckles along the chain, that there are only
gradual long wavelength bends and twists (which is the
limit of adiabatic deformations). The completeness of the
bond and torsion angles implies that, in order to determine the thermodynamical phase state of the homopolymer chain, it is sufficient to consider the response of all
the distances between all the atoms to the variations in
the bond and torsion angles only,
rαβ = rαβ (κ, τ )
Here, and in the sequel, (κ, τ ) denotes collectively all the
variables κi and τi .
Suppose that at a local extremum of the free energy,
the bond and torsion angles along the homopolymer
chain have the values
The expansion (21) may be re-arranged in terms of of the
differences in the angles, as follows:
X
F (ρ) =
{ Vk (ρk ; ρ0k ) +
k
+ Zk (ρk ; ρ0k )(ρk ρk+1 + ρk ρk−1 ) + . . . }
(22)
Here ρ0k denotes a combination of the various parameters
(κi0 , τi0 ) along the chain. But Vk (ρk ; ρ0k ), Zk (ρk ; ρ0k )
and so forth depend on the variable ρk only on the site
k; these functions are ultralocal. The terms that are
not shown explicitly, consist of higher order differences
ρk ρk+i with i ≥ 2, and higher powers of the differences.
The local terms Vk (ρk ) constitute the effective potential
X
Vef f =
Vk (ρk )
(23)
k
Start by Taylor expanding the infrared limit Helmholtz
free energy (18) around the extremum,
The structure of the effective potential is commonly used
to conclude whether a spontaneous symmetry breaking
takes place.
The transition from (21) to (22) involves, a priori, an
infinite re-arrangement of the terms in the Taylor expansion (21). In particular, the expansion (22) has been
designed so that in the continuum limit where distance
between neighboring monomers vanishes i.e. δ → 0 in
(16), it becomes, at least naively, an expansion of the free
energy in powers of momentum about the point where
momentum vanishes: For a single scalar variable ρk with
continuum limit
F [rαβ = rαβ (κi , τi )] ≡ F (κ, τ ) = F (κ0 , τ0 )
ρk → φ(x)
X ∂F
∂F
∆κk +
∆τk
+
∂κk |0
∂τk |0
the corresponding continuum limit of (22) is the derivative expansion [29]
Z 1
F (φ) =
V (φ) + (∂µ φ)2 Z(φ) + . . .
2
(κi , τi ) = (κi0 , τi0 )
Consider a conformation where the (κi , τi ) deviate from
these extremum values. The deviations are
∆κi = κi − κi0
∆τi = τi − τi0
(20)
k
+
X 1 ∂2F
∂2F
∆κk ∆κl +
∆κk ∆τl +
2 ∂κk ∂κl |0
∂κk τl |0
The free energy must remain invariant under the local
frame rotations (13), (14); the physical properties of a
chain can not depend on the framing. Accordingly, it
has been concluded [8–13] that to the leading non-trivial
order the free energy must have the form
k,l
1 ∂2F
+
∆τk ∆τl
2 ∂τk ∂τl |0
+ O(∆3 )
H=−
The first term in the expansion evaluates the free energy
at the extremum. Since (κi0 , τi0 ) correspond to the extremum, the second term vanishes. Denote in the sequel
(κi , τi ) collectively, as the variable ρi . Then,
= F (ρ0 ) +
1 X ∂2F
∆ρk ∆ρl + O(∆3 )
2
∂ρk ∂ρl |0
k,l
2κi+1 κi
i=1
+
F (κ, τ ) ≡ F (ρ)
N
−1
X
i=N
X
2
2κi + q(κ2i − m2 )2 + cdκ2i τi2
i=1
(21)
+
i=N
X
i=1
2
cτi − aτi − bκ2i τi
(24)
5
In (24) q, m, a, b, c, d depend on the atomic level physical properties and the chemical microstructure of the
homopolymer chain and its environment. In principle,
these parameters can be computed from this knowledge.
(Note the combination cd in the last term of second sum;
this choice is made for later convenience.)
The free energy (24) is the most general, universal and
gauge i.e. frame rotation (11) invariant free energy that
can be introduced to model the homopolymer in the limit
where the characteristic length scales of spatial deformations around the minimum energy configuration become
large in comparison to the distance (16) between consecutive monomers. It is essentially the naively discretized
continuum Abelian Higgs Model with one complex scalar
field, when expressed in the unitary gauge and the U(1)
gauge transformation identified with the frame rotation
(12); the term with parameter a is the Chern-Simons
term which is commonly introduced in gauge theories to
break parity. Thus, like the Abelian Higgs Model in general, it is unique in the sense of Kadanoff and Wilson.
Implicit in (24) is the assumption that there are no
abrupt wrenches and buckles along the polymer chain.
Only small, gradual bends and twists are present in deviations around the energy minimum configuration, which
is obtained by minimising the free energy (24). This is
the limit of adiabatic deformations.
The free energy (24) displays the discrete symmetry
κi → −κi
This symmetry may become spontaneously broken, by
the ground state, in line with the conventional Abelian
Higgs Model scenario; for the latter, see [30] in the
present context. Moreover, the free energy (24) is not
invariant under the local Z2 gauge symmetry (15): The
free energy (24) is the leading non-trivial contribution to
an expansion around a fixed background. To recuperate
the Z2 symmetry one should interpret e.g.
τ =
1
sin 2τ + O(τ 2 )
2
etc.
The free energy (24) can alternatively be interpreted
as a deformation of the standard energy function of the
discrete nonlinear Schr¨
odinger equation (DNLS) [31, 32].
The first two sums coincide the energy of the standard
DNLS equation, in terms of the discretized Hasimoto
variable [10]. The first (c) term in the third sum is the
Proca mass that has a claim of gauge invariance; here it
has been added as a ”regulator” for reasons explained in
[13]. The second (a) term is the helicity, and the last (b)
term is the conserved momentum. The last two terms
break the Z2 parity symmetry, they are responsible for
helicity of the homopolymer chain.
In Table I the parameter values that are kept fixed during the simulations, have been listed. These parameter
values have been chosen in conformity with those, that
are commonly encountered in the case of proteins: For a
protein, the torsion angles are much more flexible than
TABLE I: The parameters in (24) that are kept fixed during
our simulations.
q
3.5
m
1.5
b
0
d
10−4
bond angles, and the value of m in Table I corresponds
to an α-helical structure.
The parameter a, which is not fixed, is of particular
interest in the sequel. This is the parameter that breaks
chirality. Note that the momentum of the DNLS hierarchy is not considered here i.e. b = 0, for simplicity.
This term lacks a direct interpretation in the context of
the Abelian Higgs Model. It turns out that the effects of
this term are largely accounted for by the a dependent
helicity, in any case.
C.
Pseudogap
In [10] the following combination of the bond and torsion angles is considered
ψi = σi eiϑi ≡ tan
κi iϑi
e
2
(25)
where the phase is the anti-symmetrized combination
!
i
N
X
1 X
ϑi =
τk −
τk
2
k=1
k=i+1
In [10] it has been proposed to interpret (25) as the discretised version of the Hasimoto variable, in the Frenet
frames. Accordingly ψi is the complex variable that converts (24) into a generalised version of the discrete nonlinear Schr¨odinger equation [31, 32].
When a symmetry becomes broken the effective potential (23) acquired a minimum where the modulus of
the complex order parameter develops a non-vanishing
expectation value,
< |ψ| > 6= 0
(26)
while in the unbroken, symmetric phase this quantity
vanishes, at the level of the effective potential. The pseudogap state is a refinement of this phase structure [33–
37]. The system is in the pseudogap state when
<ψ>= 0
(27)
even though the modulus of ψ retains its non-vanishing
symmetry breaking value (26), at the level of the effective
potential. Such a situation occurs in the presence of a
strong phase decoherence,
< eiϑ > = 0
(28)
6
Consequently (27) can be observed even though the gap
(26) is there. But the presence of a pseudogap can be difficult to detect, in terms of the effective potential alone,
as the effective potential is commonly insensitive to the
phase. A dynamical computation that engages fluctuations is needed.
The pseudogap state is akin a symmetric phase precursor state in the broken symmetry phase. In particular the
transition between the broken symmetry phase and the
pseudogap state is commonly not a phase transition but
a cross-over prelude to the fully symmetric state that enters when the lowest energy state of the effective potential
is one where the modulus vanishes.
It should be kept in mind, that (26)-(26) are commonly
deduced by inspection of the effective potential. In the
full theory, there are corrections due to fluctuations. In
particular, in the full theory, at finite temperature, the
modulus (26) never vanishes identically. The modulus is
a positive definite quantity and thus, due to fluctuations,
it always acquires a non-vanishing value in the full theory,
as also shown in the simulations presented here.
can annihilate each other, thus any pair of kinks can be
removed by continuous deformations of the chain. A single kink can be translated, so that it becomes removed
through the ends of the chain. However, it should be
noted that on a discrete lattice the translation invariance is commonly broken, by the Peierls-Nabarro barrier
[40–43]. Thus, in general it costs (thermal) energy to
translate a kink along the chain.
The present scenario is different from the one that appears in the case of kinks in folded proteins [38, 39].
There, the parameter values in (24) are different for different super-secondary structures i.e. helix-loop-helix
motifs. A folded protein is described by a heteropolymer generalisation of (24), and the ground state is not a
straight helix such as (31). A short analysis of a simple
heteropolymer is presented in the sequel, in sub-Section
IV E.
E.
D.
Attractive interaction
Zero temperature ground state
The zero temperature ground state of the Hamiltonian
(24) is a solution to the equations of motion,
τi =
a
1
2
2c dκi + 1
κi+1 = 2κi − κi−1 + 2q(κ2i − m2 )κi + cdτi2 κi
(29)
(30)
Thus, with the present choice of parameters, the minimum energy ground state of (24) is a helix with
q
cd 2
τi ≈ ±m
κi = ± m2 − 2q
(31)
a
1
a
τi =
2c dm2 +1 ≈ 2c
The parameter m has the fixed value given in Table I
throughout. As a consequence the pitch of the ground
state helix is controlled by the ratio a/c. In the sequel
the effect of helicity on the phase structure is investigated
by varying the ratio a/c within the range a/c ∈ [0, 4π].
For some parameter values, there can be a conflict between the values of (κi , τi ) given by (31) and the steric
constraint (17). In that case the lowest energy ground
state configuration is obtained by minimising (24), subject to the constraint (17)
The range of the bond angle has been extended to negative values, by the Z2 symmetry (15). The two ground
states κi = ±m have the same energy. In addition of
these two ground states, there can also be local minima
of (24) that have the profile of a kink [38, 39] i.e. a
domain wall that interpolates between the two ground
states κi = ±m. The energy of a kink is higher than
the energy of the ground state helix (31). Two kinks
In the presence of an ambient temperature the bond
and torsion angles fluctuate. This gives rise to a nonvanishing amplitude for finite energy barrier penetration.
Thus, thermal fluctuations eventually remove any kink
there is, and the system reaches a thermal equilibrium
which minimises the Helmhotz free energy. At very low
temperature values, this is a configuration which is very
close to the ground state of (24), (17). At very low temperatures, the chain is then essentially a straight, linear
helical rod. On the other hand, in the limit of very high
temperatures the constraint (17) dominates. The chain
enters the self-avoiding random walk phase, characterised
by the mean field compactness index value ν = 3/5.
Between the low temperature rod and high temperature SARW limits there is a priori a transition region.
Recall that the conventional θ-regime is located between
the collapsed phase and the self-avoiding random walk
phase [5]. Thus this transition regime between a rod and
SARW, if it indeed exists, may be called the η-regime.
Among the goals of the present investigation is to search
for such an η-regime, and the relation between the θ- and
η-regimes and a pseudogap state.
It should be noted that the constraint (17) is purely
repulsive. There are no attractive interactions present
in the model. For a collapsed phase, an attractive interaction between the monomers must be introduced [5].
In particular, the θ-regime is a consequence of a long
distance balance between attractive interactions and repulsive interactions, to the effect that no long-range interactions prevail [5].
Thus, for a collapsed phase, the model (17), (24) is
amended with an attractive interaction. Accordingly,
(17) is generalised into an interaction potential with the
7
following profile,
U (r) = +∞
0<r<δ
U (r) = U0 {tanh(r − R0 ) − 1} δ < r < +∞
(32)
Here δ is the radius of the self-avoiding condition (17).
Hence for r < 0 the steric constraint (17) persists. But
for r > δ there is a short range attractive interaction with
strength determined by the parameter U0 . In the sequel
this parameter will be varied, jointly with the ratio a/c
that characterises helicity.
The parameter R0 which determines the range of the
attractive interaction, has the following value
R0 = 5.0 (˚
A)
throughout. It is reminded, that in all-atom molecular
dynamics simulations any long range interaction between
atoms is commonly cut off, sharply, beyond distances
around 10 ˚
Angstr¨
om or so.
The attractive interaction has a physical interpretation, in terms of ”hydrophobic” forces: In the case of e.g.
a protein chain under physiological conditions, there is
an effective attractive interaction between those amino
acids which are considered ”hydrophobic” [5–7]. Thus
in the presence of the attractive interaction (32) the energy function (24) models a chain made of ”hydrophobic” residues, with ”hydrophobicity” that depends on the
value of U0 .
It is noted that qualitatively, the present results have
been found to be quite insensitive to the details of the
profile of the potential U (r).
III.
SIMULATION DETAILS
Extensive simulations have been performed, to test
three different algorithms [44]. The steric constraint (17)
has been used in lieu of the more general (32). The aim is
to identify a Monte Carlo algorithm that has the fastest
rate of convergence towards a thermal equilibrium state.
The three Markovian algorithms that have been tested
are
1) Heat Bath algorithm
2) Metropolis algorithm
3) Mixed algorithm
The curvature and torsion angles are updated according to a probability distribution, that satisfies the detailed balance condition
of a canonical ensemble is obtained in the limit of an
infinite number of updates. Each update consists of a
”walk” through the entire chain with a provisional revision of each value (κi , τi ). The three algorithms differ
from each other only in the manner how the new values
, τinew ) are generated.
(κnew
i
1.
Heat Bath algorithm
In the Heat Bath algorithm, new values (κnew
, τinew )
i
are generated randomly, according to probability distributions
1
P (κnew
) =
exp {−βHi,κ (κnew
)}
(33)
i
i
Zi,κ
and
P (τinew ) =
1
exp {−βHi,τ (τinew )}
Zi,τ
(34)
Here Hi,κ and Hi,τ are the sum of all those terms in the
Hamiltonian (24) that contain κi resp. τi with the given
index i and the Zi,κ and Zi,τ are normalisation factors.
The updated values of κnew
and τinew do not depend on
i
the previous values of κi and τi .
The probability density for κnew
has the form
i
P (κi ) ∼ exp{−c1 κ4i − c2 κ2i − c3 κi }
(35)
where
c1 = βq
c
c2 = β(2 − 2qm2 + dτi2 − abτi )
2
c3 = β(−2(κi+1 + κi−1 ))
(36)
Thus (35) is non-Gaußian. On the other hand, the probability density P (τinew ) has the Gaußian profile
c
P (τi ) ∼ exp{−β( [dκ2i + 1]τi2 − aτi )}
2
(37)
Rejection sampling has been used to generate random
numbers according to these probability distribution: After generating κnew
and τinew the self-avoidance condii
tion (17) is checked, and the update is rejected when the
condition is violated.
2.
Metropolis algorithm
Here β is the inverse Monte Carlo temperature. The
equilibrium distribution
New values of κ and τ are generated according to
Gaußian probability distributions, which is centered at
the old values. The dispersion of each Gaußian can be
adjusted, to enhance the convergence of the algorithm.
The new values of κ and τ are accepted or rejected, in
the same manner as in the conventional Metropolis algorithm. For example, in the case of τi the probability of
acceptance of a new value is
exp{−βH(κ, τ )}
P (τi ) = min{1, exp(−β∆H)}
P ({κnew , τnew }, {κold , τold }) exp(−βH({κold , τold }))
= P ({κold , τold }, {κnew , τnew }) exp(−βH({κnew , τnew }))
(38)
8
where ∆H is the difference of the energy between the
new and the old configurations. In addition, the selfavoidance condition (17) is also verified at each step.
Mixed algorithm
The values of κnew
are generated in the same mani
ner as in the Heat Bath algorithm, while for τinew the
Metropolis algorithm is used. The convergence of the algorithm can be adjusted, by changing the dispersion of
the Gaußian distribution in the τ update.
80
Metropolis
Mixed, ∆τ = 0.02
Mixed, ∆τ = 0.20
Mixed, ∆τ = 0.50
Mixed, ∆τ = 1.00
Mixed, ∆τ = 2.00
Heat Bath
70
60
Rgyr
3.
50
40
30
A.
Comparison of algorithms
20
10
A series of comparative simulations have been performed, using the three algorithms. The following parameter values have been used in the hamiltonian (24),
(39)
These values have been chosen to reproduce a
monotonous α-helical structure of a protein, as the lowest
energy conformation; the choice is not unique. The simulations have been performed with varying chain lengths,
from N=100 to N=900. In each case, simulated annealing has been used and the initial configuration is always
a linear straight rod with
κi = τi = 0
-10
-8
-6
-4
log10T
-2
0
2
4
FIG. 1: Dependence of Rgyr on temperature for identical homopolymer chains simulated with different
algorithms. Here ∆τ is dispersion of Gaußian distribution for generation of τnew in mixed algorithm. In
conventional Metropolis algorithm the dispersions are
∆κ = ∆τ = 0.01
45
n = 100
n = 200
n = 300
n = 400
40
35
Rgyr
a = −1.0 × 10−4
c = 1.0 × 10−4
d = −5.8 × 10−4
-12
30
25
20
The initial configuration is first heated to very high temperature values (up to Tmax = 100), where the structure
is fully thermally randomised. This is followed by a slow
cooling period, to the target temperature. The cooling
takes place with small temperature steps, with each step
equal to ∆T = 0.05 in logarithmic scale. After each step,
103 − 5 × 104 Monte Carlo updates are performed along
the whole chain, to ensure that it becomes thermalised to
the ambient temperature; the longer the chain, the longer
the thermalisation. The final phase diagrams have been
calculated using 5000 Monte Carlo updates per one step
of cooling.
In Figure 1 the value of the radius of gyration (1) is
compared, as a function of Monte Carlo temperature for a
chain with N=100 monomers. The results obtained using
the Metropolis algorithm are found to be very different
from those obtained using the Heat Bath and Mixed algorithms; for small dispersion ∆τ the Mixed algorithm
coincides with the Metropolis algorithm. But when ∆τ
increases, the Mixed algorithm approaches the Heat Bath
algorithm as shown in the Figure.
There is no a priori reason why the results for the
three algorithms should be different: The stationary distribution is the same, in each of the three algorithms.
But apparently one of the three converges very slowly
15
10
10
100
1000
10000
Termalization length
100000
FIG. 2: Dependence of Rgyr on the thermalisation
length in Metropolis algorithm. The number of updates,
per one step in simulated annealing process, is shown
along the horizontal axis.
towards the equilibrium distribution. The Markovian
processes have been modelled by simulations, which are
finite in time. Thus, in order to conclude which of the
algorithms is slow to converge, the dependence of the result on the length of simulation needs to be analysed.
In the case of the Metropolis algorithm, the results are
shown in the Figure 2, for T = 10−13 which is the lowest Monte Carlo temperature value that has been used
in the present simulations. As shown in the Figure, Rgyr
continues to increase with increasing length of thermalisation. Thus, it is concluded that the results obtained
with the Metropolis algorithm apparently approach those
of the Heat Bath algorithm very slowly as a function of
the simulation time. Accordingly, it is concluded that the
9
200
35
log T = 2
log T = -2
log T = -10
160
n = 1000
n = 2000
n = 5000
n = 10000
30
25
Rgyr
Rgyr
120
80
20
15
ν2 = 0.593 ± 0.022
ν-2 = 0.622 ± 0.014
ν-10 = 0.917 ± 0.073
40
0
200
400
600
800
10
5
1000
0
0.5
Metropolis algorithm is the least optimal for the present
purposes. Either the Heat Bath algorithm, or alternatively the Mixed algorithm with sufficiently large ∆τ ,
should be used, to try and describe the thermal equilibrium configurations.
In addition, the compactness index ν i.e. the inverse
of the Hausdorff dimension has been inspected using the
three different algorithms. The results for the Heat Bath
algorithm are shown in Figure 3.
Between T = 100 and T = 1000 only very small variations are detected in the value of the radius of gyration
Rgyr ; the value of ν is essentially temperature independent, and apparently corresponds to the SARW phase.
At very low temperatures a transition to the rigid rod
state, with ν ≈ 1, is observed. This result is in line
with the general arguments that have been presented in
sub-section II E, on the expected phase structure of the
homopolymer model (24), (17).
IV.
COLLAPSED PHASE
A.
The algorithm
Extensive Monte Carlo simulations have been performed with the model (24) when subject to the attractive interaction (32). On the basis of the test simulations,
described in sub-section III A, the following algorithm is
proposed: The probability distribution is
P =
X
1
exp{−β[H +
U (~ri − ~rj )]}
Z
i<j
Here H is the Hamiltonian (24) and U is the potential
(32). The Metropolis algorithm is used for acceptance,
but with a proposal distribution that coincides with the
Heat Bath algorithm. That is, new values of κi and τi are
generated using the distributions (33) and (34). The ensuing homopolymer configuration is then accepted, provided it satisfies both the self-avoidance condition (17)
1.5
2
2.5
log10T
N
FIG. 3: Compactness index ν at different temperatures
for a homopolymer chain. Simulation performed with
Heat Bath algorithm
1
FIG. 4: Dependence of Rgyr on temperature T for a
chain with 100 monomers, in the truncated model. Here
n is the number of Monte Carlo updates per one simulated annealing step in temperature.
and the Metropolis accept-reject condition that utilises
the residual energy
X
EU =
U (~ri − ~rj )
i<j
The acceptance criterion is
exp(−β∆EU ) > λ
where λ is a random number which is uniformly distributed between 0 and 1, and ∆EU is the change in
EU under the update of κi and τi .
B.
Truncated model
The algorithm has been tested in a truncated model,
where the Hamiltonian H is removed and only the contribution from the potential U is retained. The radius
of gyration Rgyr has been evaluated, as a function of
the Monte Carlo temperature. The results are shown in
Figure 4 for several thermalization lengths. The value
U0 = 15 is used in (32) for all the runs; the exact numerical value is not relevant in these simulations, it can be
adjusted by adjusting the scale of MC temperature T A
smooth transition is observed in the radius of gyration,
from larger values at high temperatures to smaller values
at low temperatures. The Figure 5 shows how Rgyr depends on the length of polymer chain N , and fitted to the
leading order contribution in (2). It is found that Rgyr
decreases as the thermalisation length increases. The
value of the compactness index ν also decreases, and converges towards the mean field value 1/3 of the collapsed
phase, as shown in Figure 6. For high temperatures a
ν ≈ 3/5 is found, corresponding to the SARW phase.
It is concluded that the truncated model displays two
phases: At high temperatures there is the self-avoiding
random walk phase (SARW) and at low temperatures
there is the collapsed phase.
10
20
18
16
14
Rgyr
mer enters a regime of decreasing Rgyr . Finally, there is
the low temperature regime where the radius of gyration
Rgyr has a small value; see Figure 9. The compactness
index ν shows that when the thermalisation increases the
value of ν converges towards ν ≈ 1/3. This is shown in
Figure 10.
n = 1000
n = 2000
n = 5000
n = 10000
12
10
ν1000 = 0.431 ± 0.012
ν2000 = 0.3997 ± 0.0079
ν5000 = 0.3668 ± 0.0011
ν10000 = 0.3679 ± 0.0036
8
6
4
100
200
300
400
500
2.
600
N
FIG. 5: Dependence of Rgyr on the length of polymer
chain N for different thermalisation length with T = 1,
in the truncated model.
According to (31) the classical ground state profile of τi
remains intact when the parameters a and c are changed
in such a manner that the ratio a/c is constant. To study
the effect of such a change in a and c at finite temperature, in particular how it conspires with the parameter
U0 , simulations have first been performed with
a = c = 10−2 > U0 = 10−4
0.46
0.44
ν
0.42
0.4
0.38
0.36
2000
4000
6000
8000
Length of thermalization
10000
FIG. 6: Dependence of compactness index ν on thermalisation length when T = 1, in the truncated model.
C.
Parameters a and c
Full model and parameter effects
The variant of the Heat Bath algorithm described in
sub-section IV A has been used to study the effect of
various parameters to the phase structure of the full homopolymer model with attractive interaction (24), (32).
These results are summarised in Figures 7-11.
1.
Parameter U0
Figures 7-9 describe the properties of the radius of gyration for three representative values of the strength parameter U0 in (32),

 10−4
10−2
U0 =
 10−1
For the other parameters the values (39) are used. For
each value of U0 , three different characteristic regimes are
observed. In the high temperature limit the homopolymer is found in the SARW phase; this is confirmed in Figure 8. When the temperature decreases, the homopoly-
(40)
with the values (39) for the remaining parameters. Thus,
unlike in the previous simulations now the characteristic
scale of the attractive interaction is smaller than that of
the torsion angle dependent terms in the Hamiltonian.
The results for the radius of gyration are presented in
the Figure 11. At high temperatures the chain is again
in the SARW phase. Then, as temperature decreases,
there is a transition to a regime akin the intermediate
regime shown in Figure 7. Finally, there is a low temperature regime where the chain fluctuates around the
classical solution (31). The scale of transition to the low
temperature regime is controlled by parameters a and c
in Hamiltonian.
It is concluded that when the scale U0 of the shortrange attractive interaction is smaller than the scale of
the parameters a and c, the low temperature limit is described by helical structures.
D.
Analysis of different phase regimes
The bond and torsion angles form the complete set of
local order parameters to probe the phase structure (3),
in the case of the present homopolymer model. These
order parameters have the following characteristics, in
the different regimes that have been analysed in Figures
7-11; the results are summarised in Figures 12-15.
In the very high temperature SARW phase, both the
bond angle and the torsion angle are subject to large
fluctuations; the simulation results are shown in Figure
12. For the bond angles, the values are distributed in the
range
0 ≤ <|κ|> ≤ κmax ∼ 2.2
where the upper limit is due to the steric constraint (17).
As temperature increases, the values of κ become increasingly evenly distributed over this range so that in the
T → ∞ limit the distribution is fully uniform; the Figure 12 shows the bond and torsion angle distribution at a
11
24
U = 0.0001
U = 0.01
U = 0.1
22
20
50
45
40
16
Rgyr
Rgyr
18
14
35
12
30
10
25
8
6
-8
-6
-4
-2
0
20
2
-8
-6
-4
-2
0
2
log10T
log10T
FIG. 7: Dependence of Rgyr on temperature for a chain
with 50 monomers. Thermalisation length is n = 104
updates per one step of simulated annealing
FIG. 11: Dependence of Rgyr on temperature. The
chain has 100 monomers, and parameters (a, c, U0 ) are
given in (40).
90
70
Rgyr
60
50
40
30
20
ν = 0.597 ± 0.019
10
0
100
200
300
N
400
500
600
n = 1000
n = 2000
n = 5000
n = 10000
n = 50000
24
Rgyr
20
16
8
4
100
200
300
N
400
500
8
4
0
-2.5 -1.5 -0.5
0.5
1.5
12
8
4
0
2.5
-7
-4.2 -1.4
1.4
4.2
7
τ
FIG. 12: Distribution of (κi , τi ) in the high temperature
SARW phase.
generic but high temperature value. It is apparent from
this Figure that both angles are disordered.
It is concluded that the SARW phase is a disordered
phase.
A priori, the intermediate transition regime should be
the canonical θ-regime; the transition is found to take
place between temperature values within the range
10−3 < T < 101
ν1000 = 0.4719 ± 0.0045
ν2000 = 0.442 ± 0.012
ν5000 = 0.4327 ± 0.0058
ν10000 = 0.434 ± 0.015
ν20000 = 0.3932 ± 0.0013
12
12
16
κ
FIG. 8: Dependence of Rgyr on the length of chain for
high temperature T = 101.7 and with U0 = 0.01.
28
16
Probability distribution
Probability distribution
80
600
as can be seen in Figures 7 and 11. In this putative θregime the values of the bond angle are found to become
ordered. This is shown in Figure (13): The values of κi
are thermally fluctuating at around κ ≈ 1.5. The values
of the torsion angle remain largely disordered. But the
0.48
0.46
ν
0.44
0.42
0.5
κ
1.5
2.5
16
12
8
4
0
-7
-4.2 -1.4
1.4
4.2
τ
FIG. 13: Distribution of (κi , τi ) in the intermediate transition regime.
0.4
0.38
48
40
32
24
16
8
0
-2.5 -1.5 -0.5
Probability distribution
Probability distribution
FIG. 9: Dependence of Rgyr on the length of chain for
low temperature T = 10−6 and various thermalisation
lengths, and with U0 = 0.01.
1000
10000
Length of thermalization
100000
FIG. 10: Dependence of compactness index ν on thermalisation length. Temperature is equal to T = 10−6
and U0 = 0.01.
7
12
60
40
20
0
-2.5 -1.5 -0.5
0.5
1.5
25
40
15
35
10
5
30
25
0
2.5
n = 1000
n = 2000
n = 5000
n = 10000
45
20
Rgyr
Probability distribution
Probability distribution
50
80
-7
-4.2 -1.4
κ
1.4
4.2
7
20
τ
15
FIG. 14: Distribution of (κi , τi ) in the collapsed phase.
10
-8
-6
-4
-2
0
2
80
Probability distribution
Probability distribution
log10T
60
40
20
0
-2.5 -1.5 -0.5
0.5
κ
1.5
2.5
80
FIG. 16: Dependence of Rgyr on temperature for a heteropolymer chain with 150 monomers Thermalization
length n is number of updates per one step of simulated
annealing
60
40
20
16
0
-2
-1
0
1
2
15
τ
14
FIG. 15: Distribution of (κi , τi ) in the helical rod-like phase.
Rgyr
13
emergence of a structure which is indicative of a double
well effective potential is observed to be present in the
τ -distribution; see Figure 13.
The intermediate region is identified as a pseudogap
state, as described in sub-Section II C.
Finally, there are two different low temperature phases:
The collapsed phase shown in Figure 7 where the attractive short-distance interaction dominates and the helical
rod-like phase shown in Figure 11 where the attractive
short-distance interaction becomes weak.
The distribution of bond and torsion angles in the collapsed phase are displayed in Figure 14. The bond angle
is highly ordered around the classical value (31) but the
torsion angle remains disordered. However, there is an
apparent spontaneous symmetry breaking that has taken
place; the double well structure seen in the τ distribution
of Figure 13 has been removed, in a way that resembles
the familiar spontaneous symmetry breaking in a Z2 symmetric potential well.
In the helical rod-like phase, both the bond and torsion angles become peaked around the classical values
(31). The configurations are akin straight helical rods;
a little like e.g. collagen when biologically active. The
κ distribution reflects the discrete Z2 gauge symmetry.
But the Z2 symmetry in the τ distribution observed in
Figure 13 is fully broken.
The transition between the collapsed phase and the helical rod-like phase entails a transitition where the torsion
angles become ordered. Due to the very low temperature
values involved, the fluctuations are strongly suppressed
and a simulation becomes tedious. It is conjectured that,
when the temperature is kept in the low temperature
regime, an initial helical rod-like structure but with parameter values corresponding to the collapsed phase, is in
12
11
10
ν = 0.269 ± 0.039
9
8
50
100
150
200 250
N
300
350
400
FIG. 17: Dependence of Rgyr on the length of polymer
chain for low temperature T = 10−7.5 and various thermalisation lengths.
a glassy phase. Vice versa, an initial collapsed configuration with parameter values in the helical rod-like phase,
will eventually become subject to cold denaturation.
E.
heteropolymer
The effect of perturbations that break the homogeneity
of homopolymer have been investigated as follows: A
collapse has been found to take place when the parameter
U0 that characterises the strength of self-interaction is
larger than the parameters a and c that characterise the
torsion angle dependent terms. Thus, a short segment
is introduced along the chain, where U0 is less than a
and c. Accordingly, a simulation is performed where a
heteropolymer is constructed so that for a short sub-chain
of 12 monomers, the values of the parameters a and c is
increased from a = c = 10−6 to a = c = 10−2 , while
U0 = 10−4 along the entire chain. Results of simulations
are presented in Figures 16 and 17, for a chain with 150
monomers. The dependence of Rgyr on temperature is
shown in figure 16; the phase diagram is very similar to
7. For example, at low temperatures it is found that
13
-5.4
-5.7
-6
-2
log10U0
-3
-4
-5
-6
FIG. 18: (left) Homopolymer does not display helices
in the ground state. (right) Addition of a subchain that
breaks homogeneity, gives rise to a helix in the subchain
region.
4
3
a/(πc)
2
1
0
-7
-5
-3
-1
1
log10T
FIG. 19: The phase diagram on the (U0 , T, a) space.
the compactness index is close to the mean field value
ν = 1/3 of the collapsed phase. However, the geometry
of a configuration in the collapsed phase is different: As
shown in Figure 18 a helical structure appears only in
that sub-chain where the parameter values a and c have
been increased.
V.
PHASE DIAGRAM
The homopolymer phase is found to depend on three
relevant scales.
– There is the extrinsic temperature scale where the
values of κi become ordered. This scale can also be controlled intrinsically, by the parameter q in (24), but the
details have not been addressed here.
– There is temperature scale where the values of τi become ordered. This scale can be controlled intrinsically,
by the parameter ratio a/c in (24).
– Finally, the effects of the scale U0 for the short-range
attractive interactions have been investigated. This
parameter determines an intrinsic scale that controls
the transition temperature alternatively to the collapsed
phase, or to the helical rod-like phase.
Thus, by changing the relations between the three
scales the phase diagram of the homopolymer can be
identified; the phase diagram is constructed here in terms
of (T, a, U0 ). All the remaining parameters are fixed, and
given by the values in Table I.
The Figure 19 shows the three-dimensional phase diagram in the (T, a, U0 ) space. The Figures 20-23 show
various cross-sections, taken at selected values of U0 i.e.
these Figures show the phase diagram in Figure 19 on
the (a, T ) plane, with different values of U0 .
Figure 20 shows the phase diagram, for a quite large
value of U0 (strong coupling) At high temperatures there
is the SARW phase. When temperature decreases there
is the pseudogap state, that becomes the collapsed phase
at low temperatures. Between the pseudogap state and
the low temperature collapsed phase there is a θ-regime,
or rather a θ-point as it is observed only over a very
narrow temperature range.
• Figure 21 displays the phase diagram, when the value
of U0 is lowered but still relatively large (intermediate but
not weak coupling). At high temperatures there is again
the SARW phase, followed by the pseudoogap state as the
temperature decreases. At low temperatures the pseudogap state becomes converted either to the collapsed phase
or to the straight rod phase, depending on the value of
the helicity parameter a. In addition, there is a range of
values of a, when a intermediate similar to the θ-regime
is observed between the pseudogap state and the straight
rod phase. This is the η-regime. It is notable, that there
is a possibility of a 4-critical point involving the pseudogap state, the η-regime, and the collapsed and straight
rod phases.
• Figure 22 shows the phase diagram, as the value of U0
becomes further decreased (intermediate but not strong
coupling). The collapsed state has entire disappeared,
and replaced by the straight rod phase at very low temperatures. The η-regime displays a periodic structure
in the parameter a. There appears to be a tri-critical
point involving the pseudogap state, the η-regime and
the straight rod phase.
• Finally, in Figure 23 the weak coupling U0 phase
diagram is displayed. The overall topology of the phase
diagram is similar to the one in Figure 20, but with the
straight rod phase as the low temperature phase instead
of the collapsed phase.
The Figure 24 shows the phase diagram on the (U0 , T )
plane, with helicity fixed. It is notable that there might
be a 5-criticality involving the η- and θ-regimes, the pseudogap state and the collapsed and straight rod phases.
However, the detailed investigation of this region of the
phase diagram is beyond the capacity of the computer
power which is presently available to us.
In summary, it is remarkable how the phase diagram
14
2
SARW
-3
pseudogap
-2
2.5
-2.5
log10U = -3
0
2
θ-regime
-4
collapsed
-6
log10U0
log10T
-3.5
e
re
θ-
-4
-4.5
1.5
gim
collapsed
1
pseudogap
0.5
SARW
-5
-8
0
0.5
1
1.5
2
a/(πc)
2.5
3
3.5
0
straight
-5.5
rod
-6
-6
η-regime
-5
-4
FIG. 20: A cross-section of the phase diagram in 14 at
U0 = 10−3 .
-0.5
-3
1
2
-1
log10T
FIG. 24: The phase diagram on the (U0 , T ) plane, with
a = 10−4 .
2
SARW
log10U = -4.8
0
pseudogap
log10T
-2
collapsed
-4
η-regime
straight
rod
-6
-8
0
0.5
1
1.5
2
a/(πc)
2.5
3
3.5
FIG. 21: A cross-section of the phase diagram in 14 at
U0 = 10−4.8 ..
2
SARW
log10U = -5.1
0
pseudogap
log10T
-2
η-regime
-4
-6
straight rod
-8
0
0.5
1
1.5
2
a/(πc)
2.5
3
3.5
(19)-(23) is periodic in the helicity parameter a (when c
is fixed). Moreover, there is a rapid transition between
the helical phase where the compactness index ν = 1
and the collapsed phase where ν ≈ 1/3 at low temperatures, as shown in Figure (21). There is a pseudogap
state that appears as a transition regime, akin the conventional θ-regime, between the collapsed and the SARW
phases. It is conjectured that in the thermodynamical limit where the number of monomers N → ∞ the
value of the compactness index in this pseudogap state
becomes ν = 1/2. For the transition regime between
SARW and helical phases shown in Figure (22) there is
a pseudogap state that should essentially coincide in its
properties with the θ-regime pseudogap state. But due
to computational limitations the present analysis is not
sufficient to confirm this. There is an apparent 4-critical,
even 5-critical point as shown in Figures 21 and 24, but
the detailed analysis of this region in the phase diagram
needs to be performed using more extensive simulations,
which is postponed to a future project. However, it is
observed that the potential presence of a 4-critical point
in a theoretical context very similar to the present one
(i.e. Abelian Higgs model) has been reported in [26]
FIG. 22: A cross-section of the phase diagram in 14 at
U0 = 10−5.1 .
VI.
2
SARW
log10U = -6
0
pseudogap
log10T
-2
-4
η-regime
-6
straight rod
-8
0
0.5
1
1.5
2
a/(πc)
2.5
3
3.5
FIG. 23: A cross-section of the phase diagram in 14 at
U0 = 10−6
SUMMARY
The phase structure of chiral homopolymers under
varying ambient temperature values have been investigated theoretically, in terms of a universal infrared limit
energy function in combination with steric constraints
(self-avoidance) and a short-range attractive interaction
between residues. As such, the model should provide a
realistic description even in the case of heteropolymers
that display an approximatively repeating monomer pattern, provided the scale of the repeat can be considered
small in comparison to the chain length. A biologically
important example is given by collagen, the most prevalent protein in our body, in which case the short range
attractive interaction models weak hydrophobicity of the
amino acids.
15
It is found that the phase diagram displays a high level
of complexity, in terms of the parameters that control the
helicity, and the strength of the attractive interaction.
In particular, the low energy phase is either like a linear
one dimensional straight rod, or a space filling collapsed
configuration. At intermediate temperatures, there is a
state which can be identified as an example of the pseudogap state, and there are also intermediates that are more
like the conventional θ-regime. It is possible that these
regimes merge, in the thermodynamical limit, at least for
some range of parameter values. However, the possibility
of the existence of 4-critical, even 5-critical points in the
phase diagram is also proposed, but can not be confirmed
with presently available computer power.
The extension of the present approach to investigate
properties of proteins and other heteropolymers remains
a challenge to future research.
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
A.J. Heeger, Rev. Mod. Phys. 73 681 (2001)
A.J. Niemi,G.W. Semenoff, Phys. Repts. 135 99 (1986)
J. Hiraki, Journ. Antib. Antif. Ag. 23 349 (1995)
C. Li, D.F. Yu, A. Newman, F. Cabral, C. Stephens, N.R.
Hunter, L. Milas, S. Wallace Cancer Res. 58 2404 (1998)
P.G. DeGennes, Scaling Concepts in Polymer Physics
(Cornell University Press, Ithaca, 1979)
A.Yu. Grosberg, A.R. Khokhlov Statistical Physics of
Macromolecules (AIP Series in Polymers and Complex
Materials, Woodbury, 1994)
L. Sch¨
afer, Excluded Volume Effects in Polymer Solutions, as Explained by the Renormalization Group
(Springer Verlag, Berlin, 1999)
A.J. Niemi, Phys. Rev. D67 106004 (2003)
U.H. Danielsson M. Lundgren, A.J. Niemi, Phys. Rev.
E82 021910 (2010)
S. Hu, Y. Jiang, A.J. Niemi, Phys. Rev. D87 105011
(2013)
T. Ioannidou, Y. Jiang, A.J. Niemi, Phys. Rev. D90
025012 (2014)
A.J. Niemi, Theor. Math. Phys. 181 1235 (2014)
A.J. Niemi, arXiv:1412.8321 [cond-mat.soft]
L.P. Kadanoff, Physics 2 263 (1966)
K.G. Wilson, Phys. Rev. B4 3174 (1971)
K.G. Wilson, J. Kogut, Phys. Repts. 12 75 (1974)
B. Widom, J. Chem. Phys. 43 3892 (1965)
M.E. Fisher, Rev. Mod. Phys. 46 597 (1974)
P.G. De Gennes, Phys. Lett. 38A 339 (1972)
J.C. LeGuillou, J. Zinn-Justin, Phys. Rev. B21 3976
(1980)
B. Li, N. Madras, A. Sokal, Journ. Stat. Phys. 80 661
(1995)
A. Krokhotin, A. Liwo, A.J. Niemi, H.A. Scheraga Journ.
Chem. Phys. 137 035101 (2012)
M.L. Huggins, Journ. Chem. Phys. 9 440 (1941)
P.J. Flory, Journ. Chem. Phys. 9 660 (1941)
VII.
ACKNOWLEDGEMENTS
AJN acknowledges support from Region Centre Recherche d0 Initiative Academique grant, Vetenskapsr˚
adet,
Carl Trygger’s Stiftelse f¨or vetenskaplig forskning, and
Qian Ren Grant at BIT. The work of MU was supported
by the DFG grant SFB/TR-55 and by Grant RFBR-1402-01261-a. Numerical calculations were performed at
the ITEP computer systems “Graphyn” and “Stakan”
and at the supercomputer center of Moscow State University. AJN thanks A. Sieradzan and Nevena Ilieva for
discussions. We dedicate our research to the memory of
M. Polikarpov. A close colleague, mentor and teacher
without whom our collaboration would not be. Thank
you, Misha.
[25] P.J. Flory, Principles of Polymer Chemistry (Cornell
University Press, Ithaca, 1953)
[26] M. Bock, M.N. Chernodub, E.-M. Ilgenfritz, A. Schiller
Phys.Rev. B76 (2007) 184502
[27] S. Hu, M. Lundgren, A.J. Niemi, Phys. Rev. E83 061908
(2011)
[28] K. Hinsen, S. Hu, G.R. Kneller, A.J. Niemi Journ. Chem.
Phys. 139 124115 (2013)
[29] S. Coleman, E. Weinberg Phys. Rev. D7 1888 (1973)
[30] M.N. Chernodub, L.D. Faddeev, A.J. Niemi, JHEP 12
014 (2008)
[31] L.D. Faddeev, L.A. Takhtajan, Hamiltonian methods in
the theory of solitons (Springer Verlag, Berlin, 1987)
[32] M.J. Ablowitz, B. Prinari, A.D. Trubatch, Discrete
and Continuous Nonlinear Schr¨
odinger Systems (London
Mat. Soc. Lect. Note Series 302, London, 2003)
[33] V.J. Emery, S.A. Kivelson, Nature 374 434 (1995)
[34] J. Corson, R. Mallozzi, J. Orenstein, J.N. Eckstein, I.
Bozovic, Nature 398 221 (1999)
[35] N. Mannella, et.al. Nature 438 474 (2005)
[36] H. Kleinert, E. Babaev, Phys. Lett. B438 311 (1998)
[37] K. Zarembo, JETP Letters 75 59 (2002)
[38] M. Chernodub, S. Hu, A.J. Niemi, Phys. Rev. E82
011916 (2010)
[39] N. Molkenthin, S. Hu, A.J. Niemi, Phys. Rev. Lett. 106
078102 (2011)
[40] R. Peierls, Proc. Phys. Soc. 52 34 (1940)
[41] F. Nabarro, Proc. Phys. Soc. 59 256 (1947)
[42] F. Nabarro, Mat. Sci. Eng. A234 67 (1997)
[43] A.K. Sieradzan, A.J. Niemi, X.Peng, Phys. Rev. E90
062717 (2014)
[44] B.A. Berg, Markov Chain Monte Carlo Simulations and
Their Statistical Analysis (World Scientific, Singapore,
2004)