A Macromolecule in a Solvent: Adaptive Resolution Molecular Dynamics Simulation

0.27 MB·English

by  Matej Praprotnik

A Macromolecule in a Solvent: Adaptive Resolution Molecular Dynamics Simulation
Matej Praprotnik,∗ Luigi Delle Site, and Kurt Kremer
Max-Planck-Institut fu¨r Polymerforschung, Ackermannweg 10, D-55128 Mainz, Germany
We report adaptive resolution molecular dynamics simulations of a flexible linear polymer in
solution. The solvent, i.e., a liquid of tetrahedral molecules, is represented within a certain radius
fromthepolymer’scenterofmasswithahighlevelofdetail,whilealowercoarse-grainedresolution
7 is used for the more distant solvent. The high resolution sphere moves with the polymer and
0 freely exchanges molecules with the low resolution region through a transition regime. The solvent
0 molecules change their resolution and number of degrees of freedom on-the-fly. We show that
2 our approach correctly reproduces the static and dynamic properties of the polymer chain and
n surroundingsolvent.
a
J PACSnumbers: 02.70.Ns,61.20.Ja,61.25.Em,61.25.Hq

0
1
I. INTRODUCTION of the macromolecule’s structure, dynamics, and func-

] tion. To determine the interactions of a solvent with a
ft The structure of polymers in solution is determined macromolecularsolute chemistry specific interactions on
o theatomiclevelofdetailhavetobeconsidered. However,
by the solvent-polymer interaction. In the case of a
s
theresultingsolvatingphenomenamanifestthemselvesat
. nonpolar polymer in an nonpolar solvent, one typi-
t mesoscopic and macroscopic scales[2] and in the overall
a cally distinguishes three ”types” of solvent, good, Θ or
m structure of the chains. Due to large number of degrees
marginal, and poor. In the case of a good solvent the
offreedom(DOFs)suchsystemsaredifficulttotackleus-
- solvent-solvent, solvent-polymer, and polymer-polymer
d ingall-atomcomputersimulations[3]. Moreover,thevast
interactions effectively result in a situation, where the
n majority of the simulation time is typically spent treat-
chain monomers are preferably surrounded by solvent
o ingthesolventandnotthepolymerorprotein. Astepto
c molecules. As a result the chains are extended and the
[ size scales as hR2 ∝ N2νi with N being the number of bridgethe gapbetweenthe time andlength scalesacces-
monomers and ν ∼= 0.6 in three dimensions. For poor sible to simulations that still retain an atomistic level of
1 detail and the solvating phenomena on longer time and
solventoneobservesjusttheoppositeandthechainscol-
v
9 lapse into a dense globule, hR2 ∝ N2/3i. The Θ regime largerlengthscales,isgivenbyhybridmultiscaleschemes
thatconcurrentlycoupledifferentphysicaldescriptionsof
0 is where these twoeffects compensateandthe chainsbe-
the system (see, e.g., Refs. [4, 5, 6, 7, 8, 9, 10, 11]).
2 haves to a first approximation as a random walk, i.e.,
1 hR2 ∝ Ni. In the limit of N → ∞ the Θ-point is a Recently, we have proposed an adaptive resolution
0 moleculardynamics(MD)scheme(AdResS)thatconcur-
tricriticalpointinthephasediagram. Aslongasthesol-
7 rentlycouplestheatomisticandmesoscopiclengthscales
ventdoesnotinduce speciallocalcorrelationsbeyondan
0 ofagenericsolvent[12,13,14]. Inthefirstapplicationwe
unspecific attraction/repulsion and one is not studying
/
t dynamical properties the collapse of polymers is usually studied a liquid oftetrahedralmolecules where anatom-
a
istic region was separated from the mesoscopic one by a
m studied with an implicit solvent. The complicated local
flat or a spherical boundary. The two regimes with dif-
interactions are accounted for by an effective interaction
-
ferentresolutionsfreelyexchangedmoleculeswhilemain-
d betweenthechainbeads. Studiesofthatkindhavealong
n traditioninpolymerscienceandthebehaviorisnowwell taining the thermodynamical equilibrium in the system.
o understood. Beyond that there are however many situa- The spatial regions of different resolutions, however, re-
c mained constant during the course of the simulations.
tions, where it becomes difficult or even questionable, to
:
v ignorethelocalstructureofthesolvent. Solventcanplay More recently this approach was extended to the study
i an important role in the functional properties of macro- of water[15]. In the present paper, we generalize our
X
approach to the study of a polymer chain in solution.
molecules. For example, dehydration studies of proteins
r The chainis surroundedby solventwith”atomistic”res-
a solvatedinwaterdemonstratedthatatleastamonolayer
olution. When the chain moves around, the sphere of
of water is needed for full protein functionality[1]. The
atomistically resolved solvent molecules moves together
influence of a macromolecule on the structure and dy-
with the center of mass of the chain. In this way the
namics of the surrounding solvent is also an important
chain is free to move around, although the explicit res-
issue. Therefore, a detailed study of interactions of a
olution sphere is much smaller than the overall simula-
macromolecule with a solvent beyond effective coupling
tionvolume. Thisenablesustoefficientlytreatsolvation
parameters is quite often required for an understanding
phenomena, because only the solvent in the vicinity of
a macromolecule is represented with a sufficiently high
level of detail to take the specific interactions between
∗OnleavefromtheNationalInstituteofChemistry,Hajdrihova19, the solvent and the solute into account. Solvent farther
SI-1001Ljubljana,Slovenia. ElectronicMail: [email protected] away from the solute, where the high resolution is no
2
longer required, is represented on a more coarse-grained sive Weeks-Chandler-Andersenpotential
level. In this work, a macromolecule is represented by
a generic flexible polymer chain[16] embedded in a sol- Uatom(r )=
rep iαjβ
ventoftetrahedralmoleculesintroducedinRefs. [12,13].
Thisstudyrepresentsafirstmethodologicalsteptowards 4ε riασjβ 12− riασjβ 6+ 41 ; riαjβ ≤21/6σ (1)
adaptiveresolutionMDsimulationsofsystemsofbiolog- ( (cid:2)(cid:0) (cid:1) (cid:0) (cid:1) 0(cid:3); riαjβ >21/6σ
ical relevance, e.g., a protein in water.
Thepaperisorganizedasfollows: InsectionIIthedual with the cutoff at 21/6σ. σ and ε are the standard
scale model of a polymer chain in a liquid is presented. Lennard Jones units for lengths and energy respectively.
The hybrid numerical scheme and computational details riαjβ isthedistancebetweentheatomiofthemoleculeα
aregiveninsectionIII.Theresultsanddiscussionarere- andtheatomjofthemoleculeβ. Theneighboringatoms
portedinsectionIV,followedbyasummaryandoutlook in a given molecule α are connected by finite extensible
in section V. nonlinear elastic (FENE) bonds
Uatom(r )=
bond iαjα
II. MULTISCALE MODEL −21kR02ln 1− riRα0jα 2 ; riαjα ≤R0 (2)
(cid:26) (cid:2) (cid:0) (cid:1)∞(cid:3); riαjα >R0
Westudyasinglegenericbeadspringpolymersolvated
in a molecular liquid as illustrated in figure 1. Solvent with divergence length R0 = 1.5σ and stiffness k =
30ε/σ2,sothattheaveragebondlengthisapproximately
0.97σ for k T = ε, where T is the temperature of the
B
system and k is Boltzmann’s constant. For the coarse-
B
grained solvent model in the low resolution regime we
use one-site spherical molecules interacting via an effec-
tive pair potential[13], which was derived such that the
statistical properties, i.e., the center of mass radial dis-
tribution function and pressure, of the high resolution
liquid are accurately reproduced. This is also needed for
thepresentstudy,sincethemotionofthehighresolution
sphere should not be linked to strong rearrangements in
the liquid. The high and low resolution freely exchange
molecules through a transition regime containing hybrid
molecules(see figure1), wherethe molecules withno ex-
tra equilibration adapt their resolution and change the
number of DOFs accordingly[12, 13, 14].
The polymer is modeled as a standard bead-spring
polymer chain[16]. It contains N monomers, which rep-
FIG.1: (Color online) Aschematic plot ofasolvated generic
resent chemical repeat units, usually comprising several
bead-springpolymer. Thesolventismodeledondifferentlev-
els of detail: solvent molecules within a certain radius from atoms. The interactions between monomers (beads) are
the polymer’s center of mass are represented with a high defined using Eqs. (1) and (2) with the rescaled val-
u(asteodmfoisrtitch)ermesoorluetdioisntawnhtisloelvaenlotw.eTrhmeehsiogshcorpesicoluretsioonlutsipohneries ukσes2/σσB2 ==310.ε8/σσ,2R,0sBuch=thRa0tσtBhe/σsiz=eof1t.5hσeBp,olaynmderkbBea=d
B B
moveswiththepolymer’scenterofmass. Thepolymerbeads σ is approximately the same as the size of the solvent
B
arerepresentedsmallerthanthesolventmoleculesforpresen- molecule[13]. The average bond length between beads
tation convenience, for details see text.
is rescaled accordingly. The bead mass is also increased
mB = 5m0 to make them behave more like Brownian
moleculeswithin adistance r0 fromthe polymer’scenter particles. Standard Lorentz-Berthelot mixing rules[17]
of mass are modeled with all ’atomistic’ details to prop- are used for the interaction between monomers and the
erly describe the specific polymer-solvent interactions. ’atoms’ of the solvent molecules.
For the description of the solvent farther away, where
the high resolution is not required, we use a lower reso-
lution. The solvent molecules then, depending on their III. NUMERICAL SCHEME AND
distance to the polymer’s center of mass, automatically COMPUTATIONAL DETAILS
adapt their resolution on-the-fly.
Themodelsolventisaliquidofntetrahedralmolecules To smoothly couple the regimes of high and low level
as introduced in Refs. [12, 13]. The solvent molecules in of detail of the description of the solvent molecules, we
the high resolution regime are composed of four equal applytherecentlyintroducedAdResSscheme[12]. There
atoms with mass m0. Their size σ is fixed via the repul- the molecules can freely move between the regimes,they
3
are in equilibrium with eachother with no barrier in be- change. The validity condition for our approachthus re-
tween. Thetransitionisgovernedbyaweightingfunction quiresD ≪D ,whereD andD
polymer solvent polymer solvent
w(r) ∈ [0,1] that interpolates the molecular interaction the corresponding diffusion constants. This condition is
forces between the two regimes, and assigns the identity trivially fulfilled in polymeric solutions and thus also in
of the solvent molecule. We resort here to the weighting our simulations (see the next section).
function defined in Ref. [13]: WeconductedallMDsimulationsusingtheESPResSo
package[22]. We integrated Newtons equations of mo-
1; r0 >r ≥0 tion by a standard velocity Verlet algorithmwith a time
w(r)= 0; r ≥r0+d (3) step ∆t = 0.005τ and coupled the motion of the parti-
cos2[2πd(r−r0)]; r0+d>r ≥r0 cles to a DPD theromstat[21] with the temperature set
to T = ε/k . The DPD friction constant ζ = 0.5τ−1,
where r0 is the radius of the high resolution region and whereτ =(εB/m0σ2)−1/2,andtheDPDcutoffradiuswas
d the interface region width, cf. 1. The radius r0 must
setequaltothecuttoffradiusoftheeffectivepairinterac-
be chosen sufficiently large so that the whole polymer
tion between solventmolecules, i.e., 3.5σ[13]. The width
always stays within the high resolution solvent regime.
of the transition regime is 2.5σ[13]. Periodic boundary
w(r) is defined in such a way that w = 1 corresponds
conditions and the minimum image convention[17] were
to the high resolution, w = 0 to the low resolution, and
employed. Afterequilibration,trajectoriesof5000τ were
values 0 < w < 1 to the transition regime, respectively.
obtained,withconfigurationsstoredevery5τ. Thesepro-
Thisleadstointermolecularforceactingbetweencenters 9
ductionrunswereperformedwitha10 ε/σforcecapping
of mass of solvent molecules α and β:
to prevent possible force singularities that could emerge
due to overlaps with the neighboring molecules when a
F =w(|R −R|)w(|R −R|)Fex
αβ α β αβ givenmoleculeentersthetransitionlayerfromthecoarse-
+[1−w(|Rα−R|)w(|Rβ −R|)]Fcαgβ. (4) grained side [12]. The temperature was calculated using
the fractional analog of the equipartition theorem:
F isthe totalintermolecularforceactingbetweencen-
αβ
ters of mass of the solvent molecules α and β. Fex αk T
αβ hK i= B , (5)
is the sum of all pair ’atom’ interactions between ex- α 2
plicit tetrahedral ’atoms’ of the solvent molecule α and
explicit tetrahedral ’atoms’ of the solvent molecule β, where hKαi is the average kinetic energy per fractional
Fcg is the effective pair force between the two solvent quadratic DOF with the weight w(r) = α[14]. Via
αβ
molecules,andR ,R ,andRarethecentersofmassof Eq. (5) the temperature is also rigorously defined in
α β
the molecules α, β and the polymer, respectively. Note the transition regime in which the vibrational and rota-
that one has to interpolate the forces and not the in- tional DOFs are partially ’switched on/off’. The molec-
teraction potentials in Eq. (4) if the Newton’s Third ular number density of the solvent is ρ = 0.175/σ3,
Law is to be satisfied [14]. To suppress the unphysi- which corresponds to a typical high density Lennard-
cal density and pressure fluctuations emerging as arti- Jones liquid[13]. We considered three different sys-
facts of the scheme given in Eq. (4) within the transi- tem sizes with corresponding cubic box sizes: L =
tion zone we employ aninterface pressurecorrection[13]. 25.0σ,30.6σ,34.2σ. ThereducedLennard-Jonesunits[17]
The latter involves a reparametrization of the effective are used in the remainder of the paper.
potential in the system composed of exclusively hybrid
molecules with w = 1/2. Each time a solvent molecule
IV. RESULTS AND DISCUSSION
crossesaboundarybetweenthedifferentregimesitgains
or looses on-the-fly (depending on whether it leaves or
enters the coarse-grained region) its equilibrated rota- TovalidatetheAdResSapproachforthepresentpoly-
tional and vibrational DOFs while retaining its linear mer solvent system, we carried out the analysis of the
momentum[14, 18, 19, 20]. This change in resolution re- structural and dynamic properties of a polymer chain
quires to supply or remove ”latent heat” and thus must in the hybrid multiscale solvent compared to the corre-
be employedtogether with a thermostat that couples lo- spondingfullyexplicitsystemwhereallsolventmolecules
cally to the particle motion [12, 14]. This is achieved by aremodeled with a highlevel ofdetail, i.e., as a tetrahe-
coupling the particle motion to the Dissipative Particle dral molecules.
Dynamics (DPD) thermostat[21]. This bears the addi-
tional advantage of preserving momentum conservation
and correct reproduction of hydrodynamics in our nVT A. Statics of the Polymer Chain and Solvent
MD simulations. Because of the freely moving polymer
chain and solvent molecules, the above scheme requires First, we focus on the explicit (ex) systems where the
the centerofthe highresolutionspheretomovewiththe solvent is modeled with the high resolution all over the
polymerbut slowlycomparedto the surroundingsolvent simulationbox. Theseresultsareconsideredastherefer-
molecules, so that they at the boundary between differ- encetocheckhowwellAdResSproducesthesamephysics
ent regimes have enough time to adapt to the resolution as the all-atom MD simulation. The reference average
4
thermodynamic properties of the corresponding ex sys- 2.5
tems (polymer+explicitely resolvedsolvent) are listed in exN=10
exN=20
table I. 2 exN=30
N 10 20 30
1.5
m
L 25.0 30.6 34.2 Fc
<p> 2.01±0.04 2.01±0.03 2.02±0.01 D
R 1
<T > 1.0±0.01 1.0±0.01 1.0±0.01
<T > 1.0±0.5 1.0±0.3 1.0±0.2
polymer
<T > 1.0±0.01 1.0±0.01 1.0±0.01 0.5
solvent
TABLEI:Thermodynamicpropertiesofthefullyexplicitsys- 0
0 1 2 3 4 5
tems(ex) of different chain lengths N and box sizes L: aver-
r
age total pressure hpi, average total temperature of the sys-
temhTi,averagetemperatureofthepolymerhT i,and (a)
polymer
average temperature of thesolvent hT i.
solvent
10 exN=10
exN=20
The static properties of the solvent are characterized exN=30
by the solvent radial center of mass distribution (RDF)
functiondepictedinfigures2(a). Thisdistributionfunc-
tion is within the thickness of the lines the same for all (q)
systems studied, including the hybrid ones. This is to 2qS 1
be expected from our previous studies and the fact, that
the polymer fraction of volume is very small compared
to that of the solvent.
Thestatisticalpropertiesofpolymersareconveniently
0.1
described by a number of quantities, namely the radius 0.1 1 10
of gyration q
(b)
1
R2 = (r −R)2 , (6)
G N i
i 1.4 exN=10
(cid:10) (cid:11) X(cid:10) (cid:11) exN=20
exN=30
where r is the position vector of the ith monomer and 1.2
i
R=N−1 r is the polymer’s center ofmass, the end-
i i 1
to-end distance
P ρ0
/ 0.8
RE2 = (rN −r1)2 , (7) ρ
0.6
(cid:10) (cid:11) (cid:10) (cid:11)
and the hydrodynamic radius
0.4
1 1 1
0.2
= , (8)
R N2 r 2 4 6 8 10 12 14
(cid:28) H(cid:29) Xi6=j(cid:28) ij(cid:29) r
(c)
where r =|r −r |[23].
ij i j
2 2
R and R scale as
G E FIG. 2: (a) The solvent center of mass RDF for three dif-
ferentpolymerlengths(N=10,20,30), (b)thestaticstructure
(cid:10) (cid:11) (cid:10) (cid:11)R2 ∝ R2 ∝N2ν (9) factor of the polymer in the Kratky representation, and (c)
G E
the solvent density around the center of mass of the chains,
(cid:10) (cid:11) (cid:10) (cid:11)
with the number of monomers N where ν =0.5 in θ sol- which illustrates theso called correlation hole.
ventandν ≈0.588ingoodsolventconditions[24, 25,26]
withrathersmallfinite sizecorrections,while the hydro-
probestheselfsimilarstructurewithinthescalingregime
dynamic radius is known to show significant deviations
and thus provides an accurate way to determine ν. S(q)
from asymptotic behavior up to very long chains[27].
scales as
The single-chain static structure factor S(q)
S(q)∝q−1/ν →q2S(q)∝q2−1/ν (11)
S(q)= 1 exp(iq·(r −r )) (10) in the regime R−1 ≪ q ≪ b−1, where b is the typical
i j G
N
* + bond length. By fitting a power law to the computed
ij
X
5
q2S(q) plotted in figure 2 (b) we obtained the values for
ν reportedin table II. Table II summarizes the values of 18 r0=7.0
r0=11.0
allquantitiesdefinedabove,whichcharacterizethestatic 16 r0=12.0
properties of the polymer chain. The calculations were 14
performed for N =10,20,30. 12
x
ma 10
N 10 20 30 ∆r 8
L 25.0 30.6 34.2
6
h∆r i 4.2±0.8 5.7±1.0 8.1±1.4
max
RRGE ==˙˙RR−GE221¸¸11−//221 26..77±±02..50 38.8.6±±03.6 5.102±±03.8 420
RH =˙RH ¸ 3.3±0.3 4.0±0.3 4.7±0.4 0 1000 2000 3000 4000 5000
ν 0.63 0.54 0.57 t
2/z 0.59 0.71 0.67
τ =R2/(6D) 152 481 1390 FIG. 3: (Color online) Time evolution of the maximal
G
monomerdistancefromthepolymer’scenterofmass,∆r ,
max
TABLE II: Summary of some polymer (embedded in the for polymers with N = 10 beads and the radius of the high
explicitely resolved ex solvent) properties: number of poly- resolution regime r0 = 7.0 (red line), N = 20 and r0 = 11.0
mer beads N, size of the simulation box L, average maximal (green line), and N =30 and r0=12.0 (blueline).
distance of a monomer from the polymer’s center of mass
h∆r i,radiusofgyrationR ,end-to-enddistanceR ,hy-
max G E
drodynamic radius R , the static exponent ν, the exponent
H
10
2/z, wherez isthedynamicexponent,andthelongest relax-
ation time τ calculated using data from table VI. The error
bar for theexponentsν and 2/z is roughly 10%.
Another property, which directly reveals the fractal (q)
structure of the chains is the correlation hole, which is 2qS 1
showninfigure2(c). Itdirectlyshows,towhichdistance
fromthecenterofmassofthechains,the solventdensity
is perturbed by the chain beads. For the later applica- ex−ecxgNN==3300
tion of the hybrid scheme it is important to define the ex−cgicN=30
0.1
explicitsolventregimelargeenoughinordertocoverthe 0.1 1 10
correlation hole completely. q
Thevaluesofν actuallydifferslightlyfromtheasymp- (a)
totical value for the good solvent due to the finite chain
lengths. Nevertheless, the agreement improves with the
1.4 exN=30
increasing N, as expected. ex−cgN=30
Let us now turn our attention to the hybrid solvent 1.2 ex−cgicN=30
studied by MD simulation using AdResS.
1
To assure that the polymer is surrounded only by
the explicitely resolvedmolecules, we determine first the ρ0
/ 0.8
ρ
maximalmonomerdistance fromthe polymer’s centerof
mass, ∆r , as shown in figure 3 for all chain lengths 0.6
max
studied. As shown ∆r always stays within the high
max
0.4
resolution regime.
Because for N = 30 ∆rmax gets rather close to r0, 0.2
we checked the static polymer properties for that case 2 4 6 8 10 12 14 16
again. In figure 4 we compare the all explicit simulation r
tothetwohybridsimulationschemes(with(ex-cgic)and (b)
without (ex-cg) the pressure correction[13] in the transi-
tionregime)forthechainformfactorandthecorrelation FIG.4: (a)Thethestaticstructurefactorofthepolymerwith
hole. The agreement is excellent, showing that the pro- N =30in theKratkyrepresentation for all threecases stud-
posed scheme should at least be capable of properly re- ied: thefullyexplicite,theAdResSschemewithandwithout
producing the conformationalstatistics of the embedded theinterfacepressurecorrection. (b)Thecorrelation holefor
polymer in solution. thesame systems as in (a).
Thisisfirstcheckedbycomparingtheaveragethermo-
dynamicpropertiesasgivenintableIIIforthehybridex-
cg andex-cg systems(polymer+hybridsolvent). While the temperatures are identical the pressure correctionin
ic
6
the interface layer reduces the pressure slightly, so that N 10 20 30
the hybrid system now also there agrees quite well with L 25.0 30.6 34.2
the all explicit simulation. The agreement with the ref- r0 7.0 11.0 12.0
<∆r > 4.0±0.8 5.9±1.2 8.6±1.5
max
N 10 20 30 RG =˙RG2¸1/2 2.7±0.5 3.9±0.7 5.2±0.8
L 25.0 30.6 34.2 RE =˙RE2¸1/2 6.6±2.1 9.4±3.0 13.3±3.3
<p>ex−cg 2.03±0.02 2.04±0.01 2.04±0.03 RH =˙RH−1¸−1 3.2±0.3 4.0±0.4 4.8±0.4
<p>ex−cgic 2.01±0.01 2.01±0.01 2.01±0.01 ν 0.59 0.55 0.57
<T > 1.0±0.02 1.0±0.01 1.0±0.01 2/z 0.69 0.71 0.77
<Tpolymer > 1.0±0.5 1.0±0.3 1.0±0.2 τ =RG2/(6D) 152 362 1127
<T > 1.0±0.02 1.0±0.01 1.0±0.01
solvent
TABLEV:Same dataas in table IV, butnowfor thehybrid
TABLE III: Thermodynamic properties of systems with the ex-cg , where interface pressure correction is applied.
ic
polymer solvated in the hybrid ex-cg solvent and the hy-
brid ex-cg solvent: average total pressure hpi, average to-
ic
tal temperature of the system hTi, average temperature of
of freedom not only the structure but also the dynam-
the polymer hT i, and average temperature of the sol-
polymer ical properties are altered, however in a way which is
vent hT i. For the temperatures the results cannot be
solvent less understood. It is also not a priori clear, whether an
distinguished.
approach, which produces a precise coarse graining for
structural properties, does this for dynamical properties
erence values from table I is very good. This is in line
aswell. Inarecentstudyofsmalladditivemoleculestoa
with the generalstatic propertiesofthe polymers,which
polymermeltitwasshown,thatwhile the lengthscaling
aregiveninIVandV,andcompareverywelltothedata
is identical, the time scaling can be different[28]. In the
from table II.
present situation the influence of the transition regime
poses additional difficulties.
N 10 20 30
In order to determine the dynamical properties of the
L 25.0 30.6 34.2
solvent and solute we calculated the respective diffusion
r0 7.0 11.0 12.0
coefficients. The diffusion coefficient of a species is com-
h∆r i 4.2±0.8 6.6±1.2 7.7±1.3
max
RG =˙RG2¸1/2 2.7±0.4 4.0±0.6 4.6±0.7 pEuintestdeifnrormelatthioencenter of mass displacements using the
RE =˙RE2¸1/2 6.7±2.0 10.4±2.8 10.8±2.5
RH =˙RH−1¸−1 3.3±0.3 4.0±0.3 4.5±0.4 D = 1 lim h|Ri(t)−Ri(0)|2i = 1 lim h∆R2i, (12)
ν 0.63 0.58 0.54 6t→∞ t 6t→∞ t
2/z 0.56 0.69 0.62
τ =RG2/(6D) 122 381 882 whereRi(t)isthecenter-of-masspositionofthemolecule
i (which can be either a solvent or a solute molecule)
TABLEIV:Summaryofsomepolymer(embeddedinthehy- at time t and averaging is performed over all choices of
brid ex-cg solvent) properties: number of polymer beads N, time origin and, in the case of solvent, over all solvent
size of the simulation box L, radius of the high resolution molecules.
regime r0, average maximal distance of a monomer from the Figure 5 shows this for the solvent molecules’ centers
polymer’s center of mass h∆r i, radius of gyration R ,
max G of mass as a function of time for the different systems
end-to-enddistanceR ,hydrodynamicradiusR ,thestatic
E H indicated. All the curves in figure 5, except the one for
exponent ν, the exponent 2/z, where z is the dynamic expo-
thecoarse-grainedsolventcoincide. Thustheeffectofthe
nent,andthelongest relaxation time τ calculated usingdata
from table VI. The error bar for the exponents ν and 2/z is polymeronthediffusivityofthesolventmoleculesisneg-
roughly 10%. ligible. Inotherwords,thedilutionisstrongenoughthat
the polymereffectonthe solventdynamicsis verysmall.
Thecoarsegrainedsolventmoleculeshowevermovefaster
From the presented results we can conclude that
than the explicit ones. This is a consequence of the re-
AdResS faithfully reproduces the reference statics ob-
ducednumberofDOFs causingatimescaledifferencein
tained from the simulations with a polymer embedded
thedynamicsofthecoarse-grainedsystem[12,13]. While
in the explicitely resolvedsolvent.
this can be very advantageous in some cases, one can
also adjust D by an increasedbackgroundfriction in the
DPD thermostat[16, 29]. The diffusion coefficient D of
B. Dynamics of the Polymer Chain and Solvent
the solvent was obtained by fitting a linear function to
thecurvesdepictedinfigure5(a)andtheobtainedvalues
While the conformationalpropertiesof the polymer in are D = 0.036 and D = 0.057 for the explicit
bulkex bulkcg
solution are well understood and properly described by andcoarse-grainedsolvent,respectively. Thequestionto
the adaptive resolution approach, the situation for the ask here is, to what extend does this have any influence
dynamics is much less clear. By changing the degrees on the dynamics of the embedded polymer.
7
1000 100
exN=10
exN=20
exN=30
bulkex 10
bulkcg
100
>
2∆R MSD 1
<
110 0.1 <<<<<<∆∆∆∆∆∆RRRrrr222222>>>>>>eeeeeexxxxxxNNNNNN======132123000000
0.01
1 10 100 1000
1 10 100 1000
t
t
(a)
(a)
FIG. 5: Log-log plot of the time dependence of the mean 100
square displacement of the solvent molecule’s center of mass
in the. time interval [0,5000]: explicitely resolved (bulk )
ex
and coarse-grained (bulk ) solvents without solvated poly- 10
cg
mer and theexplicitely resolved solvent for the systems with
threedifferentlengthsofthesolvatedpolymer(N=10,20,30). D
S 1
M
iitccasskW,eiwnsithdhiiniicltnhuottiehsaecksconZolouiumwtninmotnttmsohoeodfdehspeyclo[rd2liyrb5om]edfetoyrhrnseaprsmoacltaiyhclmienirengrwteocerfhlaltachaitneniodddnyyswnn,haatimmhche-- 0.00.11 <<<<<<∆∆∆∆∆∆RRRrrr222222>>>>>>eeeeeexxxxxx−−−−−−ccccccggggggNNNNNN======321321000000
polymer diffusion coefficient scales as 1 10 100 1000
t
D ∝N−ν ∝R−1 ∝R−1. (13) (b)
H G
Thelongestrelaxationtimeτ =R2/(6D),i.e.,theZimm 100
G
timeτ ∝R3 =Rz,isthetimethechainneedstomove
Z G G
its own size. z = 3 is the dynamic exponent. Note that
10
the motion of inner monomers within the appropriate
scalingregimeshouldbeindependentofN. Forthemean
D
square displacements of the monomers a scaling analysis S 1
M
immediately yields for the mean square displacement of
a monomh∆erri2,i=h(ri(t)−ri(0))2i∝t2/z =t2/3, (14) 0.1 <<<<<<∆∆∆∆∆∆RRRrrr222222>>>>>>eeeeeexxxxxx−−−−−−ccccccgggggg//////iiiiiiccccccNNNNNN======123321000000
fordistancessignificantlylargerthanthebondlengthand 0.01
smaller than hR2i, i.e., times smaller than τ . For the 1 10 100 1000
Z t
center of mass of the chains a diffusive behavior for the
mean square displacement h∆R2(t)i is always observed. (c)
Althoughthechainsarerelativelyshort,atleastforN =
FIG. 6: Log-log plot of the time dependence of the mean
30 one expects a behavior relatively close to the above
mentionedidealizedscheme[27]. Figure6showsh∆r2(t)i squaredisplacementofasinglemonomer(consideredareonly
monomers near the chain’s center of mass) for polymers and
and h∆R2(t)i for polymer chains with N = 10,20,30
their centers of mass with N = 10,20,30 solvated in the ex
embedded in the different solvent scenarios studied. solvent (a), the hybrid ex-cg solvent (b) and the ex-cg sol-
In all cases the observed exponents for h∆r2(t)i and vent(c) as indicated. ic
h∆R2(t)iarecloseto0.7±0.05and1±0.05respectively.
Also the amplitudes of the displacements of the inner
beads are almost the same for all chain lengths stud-
ied. This is in good agreement with earlier studies on thermostat in our simulations. This suggests that our
different generic polymer models in a explicit solvent as hybrid scheme is also applicable to study dynamic prop-
well as studies of chains in a hybrid lattice Boltzmann erties of a polymer in a solution. Small deviations how-
solvent[6, 27]. This is to be expected since we preserve ever occur in the diffusion constant itself. The diffusion
the hydrodynamic interactions by employing the DPD constantsforthepolymerswithN =10,20,30inthehy-
8
brid ex-cg and ex-cg solvents were obtained by fitting V. SUMMARY AND OUTLOOK
ic
the straight curve to the polymer’s center of mass mean
square displacement presented in figure 6. The fit yields
In this paper we presented a hybrid multiscale MD
data listed in table VI. The corresponding static and
simulation of a generic macromolecule in a solvent us-
dynamic exponents and the longest relaxation times are
ing the recently proposed AdResS method. The solvent
given in tables IV and V.
surroundingthemacromoleculeisrepresentedwithasuf-
ficiently high level of detail so that the specific interac-
N 10 20 30
tions between the solvent and the solute are correctly
D(ex) 0.008 0.005 0.003
taken into account. The solvent farther away from the
D(ex-cg) 0.009 0.006 0.0045
D(ex-cg ) 0.0085 0.006 0.0035 macromolecule, where the high resolution is not needed,
ic
is represented on a coarse-grained level. The high and
lowresolutionregimesfreelyexchangesolventmolecules,
TABLE VI: Diffusion constant of the polymer chain embed-
ded in three different solvents: explicitely resolved ex, hy- which change their resolution accordingly. To correctly
brid ex-cg, and hybrid ex-cg . Though the statistics of the simulate momentum transport through the solvent, we
ic
data is rather poor we can estimate the error bar roughly to use the DPD thermostat. The simulation results show
10−15%. For comparison, the diffusion coefficients of the that AdResS accurately reproduces the thermodynamic
explicit and coarse-grained solvents are Dbulkex = 0.036 and and structural properties of the system. The presented
Dbulkcg = 0.057, respectively with an effect of the polymers methodology is an extension of AdResS to simulations
too small to determine here. Hence D ≪D .
polymer solvent of a solvation cavity, and represents a first step towards
thetreatmentofmorerealisticsystemssuchasbiomacro-
While the ratio of the diffusion constants for different molecules embedded in water. Work along these lines is
chain lengths roughly follow the expected scaling, even already under way[15].
though it cannot hold precisely due to the different box
sizes,wehereobserveatendencytoaweaklyaccelerated
diffusion in the hybrid regime. This is most evident for
thehybridex-cg case. Twodifferentaspectsmightplaya
role here. First the viscosity in the coarse grained outer
ACKNOWLEDGMENTS
regimeissmaller,whichmusthaveaneffectonthediffu-
sion. Second,thesmallpressureanddensityfluctuations
in the transitionregimemight contribute to the effect as We thank Thomas Vettorel, Vagelis Harmandaris,
well. BenedictReynolds,andBurkhardDu¨nwegforusefuldis-
Althoughthisisonlyaveryfirstandincompletetest,it cussions. ThisworkissupportedinpartbytheVolkswa-
showsthatwithintheAdResSschemeessentialaspectsof genfoundation. Oneoftheauthors(M.P.)acknowledges
thedynamicalpropertiesoftheembeddedpolymerchain thefinancialsupportfromthe statebudgetbythe Slove-
are reasonably well reproduced. nian research Agency under grant No. P1-0002.
[1] G.Careri, inHydrationProcessesinBiology: Theoretical Coveney, Phys. Rev.Lett. 97, 134501 (2006).
andExperimentalApproaches,editedbyM.C.Bellisient- [12] M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem.
Funel,pages 143–155, IOS,Amsterdam, 1999. Phys. 123, 224106 (2005).
[2] P.Das,S.Matysiak,andC.Clementi, Proc.Natl.Acad. [13] M.Praprotnik,L.DelleSite,andK.Kremer, Phys.Rev.
Sci. U.S. A.102, 10141 (2005). E 73, 066701 (2006).
[3] E. Villa, A. Balaeff, and K. Schulten, Proc. Natl. Acad. [14] M.Praprotnik,K.Kremer,andL.DelleSite, Phys.Rev.
Sci. U.S.A.102, 6783 (2005). E 75, 017701 (2007).
[4] H.Rafii-Tabar,L.Hua,andM.Cross,J.Phys.: Condens. [15] M. Praprotnik, L. Delle Site, K. Kremer, S. Matysiak,
Matter 10, 2375 (1998). and C. Clementi, cond-mat/0611544 (2006).
[5] J. Q. Broughton, F. F. Abraham, N. Bernstein, and [16] K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057
E. Kaxiras, Phys.Rev. B 60, 2391 (1999). (1990).
[6] P. Ahlrichs and B. Du¨nweg, J. Chem. Phys. 111, 8225 [17] M.P.AllenandD.J.Tildesley, Computer Simulationof
(1999). Liquids, Clarendon Press, Oxford, 1987.
[7] A.MalevanetsandR.Kapral, J.Chem.Phys.112,7260 [18] D. Janeˇziˇc, M. Praprotnik, and F. Merzel, J. Chem.
(2000). Phys. 122, 174101 (2005).
[8] G. Csanyi, T. Albaret, M. C. Payne, and A. DeVita, [19] M. Praprotnik and D. Janeˇziˇc, J. Chem. Phys. 122,
Phys.Rev.Lett. 93, 175503 (2004). 174102 (2005).
[9] C. F. Abrams, J. Chem. Phys. 123, 234101 (2005). [20] M. Praprotnik and D. Janeˇziˇc, J. Chem. Phys. 122,
[10] M. Neri, C. Anselmi, M. Cascella, A. Maritan, and 174103 (2005).
P.Carloni, Phys. Rev.Lett. 95, 218102 (2005). [21] T. Soddemann, B. Du¨nweg, and K. Kremer, Phys. Rev.
[11] G. D. Fabritiis, R. Delgado-Buscalioni, and P. V. E 68, 046702 (2003).
9
[22] http://www.espresso.mpg.de. [25] M. Doi and S. F. Edwards, The Theory of Polymer Dy-
[23] Note that the definition of R given in Eq. (8) is only namics, Clarendon, Oxford, 1986.
H
correctforaninfinitesimulationbox.Forafiniteboxwe [26] A. D. Sokal, in Monte Carlo and Molecular Dynam-
havetoalsoconsiderthehydrodynamicinteractionswith ics Simulation in Polymer Science, edited by K. Binder,
theperiodicimages.Thevalueofinsuchawaycorrected chapter 2, Clarendon, Oxford, 1995.
R deviates from the value obtained using the formula [27] B. Du¨nweg and K. Kremer, J. Chem. Phys. 99, 6983
H
(8)[6, 27]. However, since our aim in the present work is (1993).
onlytocomparethepropertiesofthepolymerchainfrom [28] V. A. Harmandaris, N. P. Adhikari, N. F. A. van der
the hybrid simulation with the corresponding ones from Vegt, K. Kremer, R. Voelkel, H. Weiss, and CheeChin
theall-atomsimulationweheredonottakethefinitesize Liew, Macromolecules, submitted.
correction into account and use the formula (8) for the [29] S. IzvekovandG. A.Voth, J. Chem. Phys. 125, 151101
definition of R . (2006).
H
[24] P.-G. de Gennes, Scaling Concept in Polymer Physics,
Cornell University,Ithaca, 1979.