DOI:
10.1039/C5SM02604K
(Paper)
Soft Matter, 2016,
12, 24162439
Open boundary molecular dynamics of sheared starpolymer melts
Received
20th October 2015
, Accepted 15th January 2016
First published on 15th January 2016
Open boundary molecular dynamics (OBMD) simulations of a sheared star polymer melt under isothermal conditions are performed to study the rheology and molecular structure of the melt under a fixed normal load. Comparison is made with the standard molecular dynamics (MD) in periodic (closed) boxes at a fixed shear rate (using the SLLOD dynamics). The OBMD system exchanges mass and momentum with adjacent reservoirs (buffers) where the external pressure tensor is imposed. Insertion of molecules in the buffers is made feasible by implementing there a low resolution model (blobmolecules with soft effective interactions) and then using the adaptive resolution scheme (AdResS) to connect with the bulk MD. Straining with increasing shear stress induces melt expansion and a significantly different redistribution of pressure compared with the closed case. In the open sample, the shear viscosity is also a bit lowered but more stable against the viscous heating. At a given Weissenberg number, molecular deformations and material properties (recoverable shear strain and normal stress ratio) are found to be similar in both setups. We also study the modelling effect of normal and tangential friction between monomers implemented in a dissipative particle dynamics (DPD) thermostat. Interestingly, the tangential friction substantially enhances the elastic response of the melt due to a reduction of the kinetic stress viscous contribution.
1. Introduction
Die swelling^{1} is a well known phenomenon in polymer melts and most viscoelastic liquids which consist of the sudden expansion of the liquids after exiting a slit or orifice. The most frequent explanation has a microscopic origin: molecules elongate in the stream direction and compress perpendicularly exerting extra elastic pressure in the normal planes. This leads to the socalled normal stress differences which are the landmark of viscoelasticity. Despite this accepted view, the devil is in the details and although considerable effort has been carried out since the middle of the last century (see e.g.ref. 1 for historic details), accurate modelling of polymer melts is a very difficult task. Any of the many constitutive relationships^{1} for continuum models cannot generally predict the rheology of a new molecular polymer design. In turn, the huge span in time scales in any standard polymer melt limits the scope of molecular dynamics (MD) to simple rheological tests with extremely small samples under simple (usually steady) flows. However, a detailed account of bonded and nonbonded interactions in atomistic (AT) simulations (see e.g.ref. 2–4) is able to grasp relevant information, maybe then to feed continuum (fluid dynamics) models. In between, coarse grained (CG) molecular modelling is useful due to many reasons. We give at least a couple of reasons: first, polymer science has some degree of universality which benefits the use of simplified models, quite often able to provide insight and valuable predictions.^{4,5} Second, the theory of coarse graining to extract precise coarse potential interactions pertaining to the atomistic model at hand is now advancing at a relatively fast rate.^{4,6,7} More recently, the community has started to recognize the relevance of the dynamic aspects of coarse graining either based on GENERIC^{8} or the Mori–Zwanzig formalism.^{9} The idea is to perform short atomistic simulations^{8} to extract the quantities determining the reversible and irreversible dynamics of the slow variables, such as the friction kernels.^{9} These friction kernels (between polymer “blobs”) are naturally implemented in the dissipative particle dynamics (DPD) method as, notably, Español et al.^{9} showed (under the Markovian and pairwise approximations) that DPD can be formally derived from the Mori–Zwanzig coarse graining route. Another relatively newer route is the use of hybrid models concurrently combining continuum and molecular simulations. Yasuda et al.^{10} and other (mostly Japanese) groups have exploited a version of these hybrids to model polymer melts. We have also worked in this field using essentially the same technique used in this work for simple fluids.^{11,12} Contrary to Yasuda et al.^{10} and other hybrid schemes, the present method is designed to open up the simulation box so as to consider melt expansions. Here, we apply the idea to study the rheology of much larger, polymeric molecules. At present, we restrict to the open boundary molecular dynamics (OBMD) simulation,^{13} without connection to the continuum side. However, the present method naturally connects with continuum fluid dynamics via hybrid schemes.^{14–20}
Another aspect of this work is the use of DPD as a tool to check the effect of monomer friction on the rheological behaviour of these CG models. In particular, the tangential friction between blobs naturally arises while performing dynamic coarse graining^{9} but is seldom included in these sort of analyses (see for example the DPD study for melts by Fedosov et al.^{21}). To conclude our comments on methodological aspects, we also note that the present simulations in open domains use a useful trick which consists of using an even lower detailed molecular model to feed the molecular reservoir close to the open system boundaries. The idea is taken from the socalled “adaptive resolution scheme” (AdResS)^{22–34} and permits to generalize the use of AdResS in standard periodic (closed) boxes (see e.g.ref. 11, 12 and 35–37).
The main physical question we pose here is what are the rheological consequences of imposing a fixed pressure load to a sheared sample of polymer melt, compared to the case of shearing at a fixed volume. This question was also raised by a bunch of groups spread over the last two decades^{38–43} with contradictory results. It is indeed a particularly relevant question for molecular simulations because the vast majority of numerical studies on melts have considered closed (usually periodic) systems, while many rheological experiments are carried out under a normal load: melts across slabs in ref. 44 (cited experiments therein) or a cone and plate rheometry under a fixed load.^{1}
From the fundamental side, this question connects with the already mentioned die swelling phenomenon whose details are still not completely understood. It also connects with another interesting question which is, what consequences do boundary constraints have on flowing (farfromequilibrium) polymer melts (see ref. 45 for a recent study). The OBMD is a flexible tool for these questions on boundary constraints because it can be tuned to permit different ensembles^{46} such as the grand canonical, isoenthalpic, isothermal, constant stress (Neumannlike) or constant shear (Dirichletlike).^{14} In particular, it could be useful to validate theories for nonequilibrium thermodynamics (such as extended thermodynamics^{47,48}), or the far from trivial fate of fluctuations of mass (related to sound modes) and momentum in sheared complex fluids, which often lead to undesirable instabilities in sheared or extruded melts, like the shark skin.^{1,49}
Even in unentangled melts, the influence of boundary or global constraints on the density expansion of sheared melts is still poorly understood with the studies presenting contradicting results. At least there is a consensus on the fact that for a given shear, the shear viscosity η is larger in the isochoric (NVT) constraint than under either a constant pressure^{39,40} or a constant load.^{43,44} A clear manifestation of this effect was presented in a numerical study of Thompson et al.,^{44} in which a slab of lubricating liquid (20mers) flowing between two solid walls presented a shear thinning exponent β_{η} ≃ 2/3 under a constant normal load, while just 0.5 under a constant volume (here η ∼ ^{−βη} with the shear rate). This effect has been ascribed to the shear dilatancy manifesting in a larger hydrostatic pressure under a constant volume.^{39,40} Indeed, the viscosity of polymer melts often increases with the pressure^{49} as also observed in ref. 44. In Section 6, we offer a more precise analysis showing that the spring stress is reduced in the open system due to the smaller intermolecular friction in the expanded melt.
There are several kinds of “isobaric” conditions: a constant hydrostatic pressure (P_{iso}) is usually termed isobaric, while a constant normal load (here P_{22}) is closer to industrial processes like slit extrusion^{1} or lubricants.^{44} Both can show substantially different rheological behaviours when compared with constant volume studies (see e.g.ref. 44). Experimental studies of the Couette flow under a constant load are also scarce and indicate shear expansion and a measurable increase of the melt viscosity with increasing external pressure.^{50,51}
The most striking differences in the computational literature are found in the density variation with shear. Dlugogorski^{38} (for a FENE dumbbell) and Daivis and Evans (modelling decane)^{39} seem to be the first to perform a molecular simulation showing the density decrease under shear (using isobaric conditions). Daivis and Evans use the term “shear dilatancy”, following the term used by Reynolds on the same phenomenon (see the quotation in ref. 39). Note that “shear dilatancy” has been later often used instead of “shear thinning”, but here it is not. By contrast, Xu et al.^{40} (attractive linear chains up to 50mers) reported just the opposite result (compression under shear). For (purely repulsive) branched chains (under NVT), they also reported a reduction in the hydrostatic pressure with shear. A subsequent study by Matin^{43} for linear chains at a constant load (and chain lengths up to 50 monomers) found shear dilatancy and also a nonmonotonous trend for the hydrostatic pressure (as Xu et al.^{40} and others^{2,3} did). Shear dilatancy was later also found by Bosko et al.^{41} while analysing dendrimer melts under isobaric conditions. Consistently, they found a pressure increase under a constant volume. Shear dilatancy is the trend we also observe but, under our imposed constant normal load (which is closer to the experimental setup^{1,39,44}) we find that the density (and viscosity) is controlled by the load and not by the hydrostatic pressure. However, there is a lack of studies and the question remains about universality of shear dilatancy of polymer melts under a Couette flow (consistent with the die swelling phenomenon under a Poisson flow^{1}).
To conclude the introduction some words should be said about the technological relevance of star molecules (see discussion in ref. 52–54) which certainly arises from their unique dynamic features. Star molecules present a broad range of relaxation times associated with different molecules pulsations (rotation, elastic deformation, and arm entanglements) analysed in ref. 53. Each relaxation time triggers a change in the rheological regime once the external perturbation (shear) exceeds the corresponding threshold rate. The present results also illustrate this phenomenon. Moreover, star molecules bridge the gap between linear polymers and colloids^{54,55} and can present interesting (colloidallike) ordering effects, sometime enhanced due to their softer character.^{53,56–59} In this context, a suggesting observation in these simulations is the onset of some ordering in the neutral direction at large shear rates whose origin (hydrodynamic or entropic) remains to be established.
In Section 2, we briefly present the OBMD method which was otherwise more thoroughly explained in ref. 13. The star molecule melt model presented in Sections 3 and 4 shows that the OBMD correctly reproduces thermodynamic equilibrium according to the grand canonical ensemble. In Section 5, we present the results for sheared melts in the absence of tangential friction between monomers and analyze the results in Section 6 according to the pressure balance. This serves to enlighten the discussion on the effect of tangential friction in Section 8. Finally, Section 9 discusses some interesting results obtained for melts presenting severe viscous heating that depends on the characteristics of the DPD thermostats and their friction kernels. Comparison with the previous results is made in Section 10, while conclusions and future outlook are given in Section 11.
2 Open boundary molecular dynamics
We begin by briefly explaining the OBMD method which combines features of the open MD^{46} and adaptive resolution.^{11} The reader is referred to a review on open MD^{14} and ref. 13 for a more detailed presentation of the present OBMD implementation to star polymers. The OBMD simulation is carried out in an open rectangular box which, in the present setup, permits the MD domain to exchange mass and momentum through two of its boundaries (along the x_{2} direction) with a reservoir (called buffer) which is maintained at some desired thermomechanic state. The two domains of finite extent of the buffer embed the central part of the box (the MD domain). They allow for molecular insertion or deletion so as to keep their average molecular density fixed (typically to a fraction between 0.5 and 0.7 of the bulk density). The OBMD is therefore not periodic in the coupling direction. Molecules are free to enter or leave the buffer from or to the MD domain, but in doing so they cross another layer where they gradually change their atomistic resolution (73 monomers for the star molecule considered hereby) to a reduced CG model, comprising one only spherical “blob” per molecule. Obviously the CG layer is placed inside the buffer domain (which here also contains a smaller atomistic part). This strategy permits to perform an otherwise impossible task: the insertion of new polymeric molecules into the melt. New molecules are inserted into the low resolution layer of the buffer, where soft CG interactions govern the dynamics of the blobmodel polymers. Soft effective interactions can be obtained from the Boltzmann iteration procedure,^{60,61} although we shall see that in principle, the consistency of OBMD (in terms of pressure balance across the layers) does not depend on the CG potential chosen. The insertion of these blob molecules is carried out by the USHER scheme^{62} and the change from CG to monomer molecular resolution (usually termed atomistic resolution, AT) is carried out by the AdResS.^{22–34}
The dynamics of the monomers can be described by the following equations of motion:

 (1) 

 (2) 
Here
r_{i} denotes the position of the
ith particle,
v_{i} its velocity, and
m_{i} its mass. The total force acting on this particle has three contributions: the external force
F^{ext}_{i} acting only on the particles at the buffer (to impose the desired momentum flux); the adaptive resolution force
F^{ad}_{i}, which accounts for all types of particle–particle interactions, and the thermostat
F^{th}_{i} contribution (here, it is applied to the whole system).
The adaptive resolution force F^{ad}_{i} is constructed to allow for a momentum conserving “alchemic” transformation of the molecules, which takes place gradually along the transition layer (where 0 < w < 1, see below). The transition is achieved by the following linear combination of the AT and CG forces (F^{AT}_{αβ} and F^{CG}_{αβ}, respectively),^{22}

F^{ad}_{αβ} = w(x_{α})w(x_{β})F^{AT}_{αβ} + (1 − w(x_{α})w(x_{β}))F^{CG}_{αβ}.  (3) 
Both expressions are correspondingly weighted by a position dependent function
w(
X), whose value equals 0 in the CG region and 1 in the AT one and gradually changes in between (transition layer). The adaptive resolution force provided by the AdResS is not derived from a Hamiltonian and does not conserve the energy.
^{23} It is however constructed to obey Newton's third law, which ensures the conservation of total linear momentum of the system and can thus be used to study fluid flows. This fact is unimportant for the present study which targets isothermal sheared systems (notconserving energy anyhow). The OBMD model might be generalized by using the recent Hamiltonian AdResS (HAdResS)
^{63,64} which is based on a global Hamiltonian (
i.e. also allows Monte Carlo simulations). In such a hypothetical case, extra care should be taken with the momentum conservation because of the presence of driftforces in HAdResS coming from the free energy difference between the AT and CG models.
^{65}
An essential function of the buffer region is the imposition of boundary conditions to the open MD box. This is done by adding an extra “external” force at the buffer regions, F^{ext}, calculated from eqn (4) (see e.g.ref. 66).

 (4) 
Here,
P_{out} and
P_{in} represent the total linear momenta of the particles that were removed and inserted into the simulation in the last time step of integration
δt.
is the momentum flux tensor, while
n is the outward normal vector of an openend plane of the box.
^{11} In general, the pressure tensor contains normal and tangential contributions,
i.e.. The external force is designed to exactly conserve the linear momentum over the whole particle system (buffer + MD) and it is distributed among the buffer particles according to
F^{ext}_{i} =
G(
x_{i})
F^{ext}. To allow a different distribution of the normal and tangential forces, the distribution function
G is chosen to be a tensor defined by
eqn (5),
^{13} 
 (5) 
with
g_{∥} determining the spatial distribution of the normal force and
g_{⊥} the distribution of shear stress. Both functions are depicted in
Fig. 2.
Many OBMD applications (as those illustrated hereby) involve transfer of momentum (the pressure tensor) from outside the MD domain. This requires a momentum conserving thermostat. In production runs, we used the DPD thermostats (while a strong damping Langevin for equilibration purposes). Our choice for the DPD thermostat is not only to conserve momentum (in principle one could use the LoweAndersen^{67}) but also because of modelling purposes. The Mori–Zwanzig formalism, under Markovian conditions leads to coarsegraining dynamics with DPDlike equations of motion.^{9} The solid theoretical results show that the friction kernels introduce an important modelling aspect. The friction kernels in a CG model of some real melt, should ideally be measured from force–force correlations of the detailed allatom model.^{9} Here, we adopt a simpler but yet useful route, which is to study how friction affects the rheology of the model melt. The generic form of the DPD thermostat used is

 (6) 
where
_{ij} is the fluctuating force constructed to satisfy the fluctuationdissipation under equilibrium conditions. We refer the reader to
ref. 68–70 for details on the DPD implementation. As in
ref. 9, 69 and 70 the friction kernel has normal and tangential components,
Γ = γ_{∥}n_{ij}n_{ij} + γ_{⊥}t_{ij}t_{ij}, 
where
n_{ij} is the vector joining two monomers
i and
j and
t_{ij} determines the directions in the perpendicular plane. The kernels
γ_{∥} and
γ_{⊥} are distance dependent. Its shape (for a particular allatom model) can be obtained from dynamic CG,
^{9} here we will use Heaviside functions with a certain cutoff distance (see
Table 1 for details). Most thermostats require or introduce some form of friction, albeit, the great majority of simulations of polymer melts do not consider thermostating (and its added friction) as part of the molecular model, but just a way to remove the heat dissipated under shear. Also most DPD simulations, such as the relatively recent work on sheared melts
^{21} do not introduce the tangential friction between blobs, but rather take the form of the most standard DPD kernel (normal friction alone). However, as also pointed out by Padding and Briels,
^{5} friction should be considered as a part of the CG model. Indeed, as shown by Hijon
^{9}et al. (CG of a star molecule as the unit blob), tangential and normal frictions can be quite different from each other. Here, we start to explore how tangential and normal frictions affect the rheology of a star molecule under farfromequilibrium conditions.
Table 1 Thermostats used in simulations. Standard means a standard DPD thermostat and the adaptive DPD is explained in eqn (7). The transverse DPD thermostat from ref. 70 is denoted as “tdpd”
Label 
Kernel cutoff (R^{DPD}_{cut}) 
γ
_{∥}

γ
_{⊥}

Sdpdshort 
2^{1/6}σ 
[1.0–20.0] 
0 
Sdpdlong 
1.5 × 2^{1/6}σ 
1.0, 5.0 
0 
Adpd 
2^{1/6}σ 
1.0, 5.0 
0 
Tdpd 
1.5 × 2^{1/6}σ 
1.0 
1.0 
In the following sections, we present and analyse the results obtained for the DPD thermostat with no tangential friction between monomers (blobs). This analysis is then used to understand the effect of the tangential friction, considered in Section 8. In the case of the normal friction alone, the substantial heat dissipated by the sheared melt requires from us to implement a slight modification of the DPD thermostat to maintain a constant (kinetic) temperature at the largest shear rates. Details of this modified DPD thermostat (we call it “adaptive”) and friction kernels are given in Section 3.
Importantly, simulations in closed isothermal (periodic) boxes (NTV) are used as a reference to investigate the effect of open boundaries. Shear flow in the closed system has been simulated using Lees–Edwards boundary conditions^{71} and the SLLOD algorithm.^{72,73}
3 Setup and melt models
The simulation setup is illustrated in Fig. 1 and 2. The polymer melt is exposed to a Couette flow in the x_{1} direction being sheared along the x_{2} direction (gradient direction). The vorticity or neutral direction is x_{3}. In the closed box, we use the SLLOD^{72,73} dynamics to impose the desired shear rate in a closed periodic box (with a constant particle number and volume). In the open setup, the surfaces located at x_{2} = ±L_{2}/2 are subjected to equal normal pressures p^{ext}_{22} and (opposite sign) tangential stresses ±p^{ext}_{12}, in such a way that the rotational part of the shear flow turns counterclockwise in the flow–gradient (x_{1} − x_{2}) plane. No constraint is imposed to the remaining component of the pressure tensor, resulting in 〈p_{13}〉 = 0 and a self determined 〈p_{33}〉. The box is periodic in the other two directions x_{1} and x_{3} so this setup corresponds to an slice of polymer melt with a fixed load at two of its boundaries (at x_{2} = ±L_{2}/2).

 Fig. 1 Schematic representation of the open system at equilibrium along the longitudinal direction. All monomers and polymers in the system are of the same kind. Different colors are used for the sake of clarity of the picture. In the region of interest polymers are represented in the highest, i.e. atomistic (AT) resolution. Buffer regions (see text) are heterogeneous, i.e. containing regions of different resolutions. The change of resolution from AT (dots) to coarsegrained, i.e. CG, (spheres) occurs in the hybrid region (HY) of each buffer, carried out by the AdResS. New molecules are inserted into the CG part of buffers. The system is open at both ends of the box. The upper part of the figure depicts schematic representations of a molecule inside regions of different resolutions.  

 Fig. 2 Buffer distribution function. The force F^{ext}, that acts on buffer regions in order to impose the boundary conditions at open ends of the box to the region of interest, is distributed among molecules inside buffers via the depicted functions, forming the distribution function tensor G, given by eqn (5). The force acting on each molecule equals F^{ext}_{i} = G(x_{i})F^{ext}.  
The melt is made of the star polymer model already presented in ref. 9 and 74. Each polymer consists of 73 monomers, i.e. 12 arms of 6 monomers attached to the central monomer. In what follows, we use m_{0}, σ_{0} and ε_{0} for mass, length and energy units and we will arbitrarily set these units to m_{0} = 1, σ_{0} = 1, and ε_{0} = 1, respectively. The resulting time unit is τ_{0} = σ_{0}(m_{0}/ε_{0})^{1/2} = 1. Excluded volume interactions of monomers are modelled by the repulsive Weeks–Chandler–Anderson interaction with diameter parameter σ = 2.415σ_{0} and energy parameter ε = 1. The interactions between two adjacent bonded monomers are harmonic with a linear spring of stiffness constant K = 20.0ε/σ_{0}^{2}. The equilibrium distance between noncentral monomers is r^{eq}_{ij} = 2.77σ_{0}, while the equilibrium distance between the central monomer and the first monomer of an arm is r^{eq}_{ij} = 3.9σ_{0}. The size of the simulation box is 390 × 117 × 117 (in units of σ_{0}).
Simulations were performed at a fixed constant (monomer's kinetic) temperature (T = 4.00 ± 0.01). The results presented in the following sections correspond to simulations obtained using a modified DPD thermostat which we label as “adpd”. This adpd thermostat has no tangential friction γ_{⊥} = 0 while γ_{∥} = 1 (some results also with γ_{∥} = 5). The cutoff of the Heaviside friction kernel is R_{dpd} = 2^{1/6}σ. We refer to Table 1 for thermostat and kernel details. We recall that the effect of the tangential friction is analysed in Section 8. In Section 9, we illustrate the heat dissipation and temperature increase observed upon using standard thermostats DPD with normal friction.
Thermostating sheared polymer melts is a delicate issue due to the large amount of heat they dissipate. At a large enough shear rate or large shear stress, the temperature of the melts increases. The same phenomenon is also observed in experiments and industrial processes (extrusion) at high shear rates (typically above 500 Hz).^{47,75,76} Phenomenological temperature “corrections” for the melt viscosity are often used in industry and experiments.^{77} Although this problem goes beyond the present manuscript, the temperature of a system under a nonequilibrium steady state is also a fundamental problem because equipartition is lost and the different temperature definitions present slight variations (kinetic versus configurational temperature^{78}). In the literature, few studies report on the problem of viscous heating in molecular simulations (see the exceptions in ref. 79–81) and many published materials elude reporting on possible temperature variation in their sheared thermostated systems. After the present experience, we believe that some of the data presented in previous papers might be somehow biased by temperature. We will show an indication later in Section 9.
We considered four different thermostats: the “standard” (sdpd)^{69} and “transverse” (tdpd)^{70} DPD thermostats and also a modified DPD thermostat. The latter is able to extract larger amounts of heat by selfadapting its temperature parameter T_{DPD} which controls its random force term. This “adaptive DPD thermostat” (adpd) as we call it, dynamically adjusts T_{DPD} according to a sort of coupled heat equation,

 (7) 
where
T_{tgt} = 4 is the temperature of the target system and
T_{MD} is the kinetic temperature obtained from the variance of the peculiar velocities of the monomer
r_{i},
u_{i} =
v_{i} −
v_{f}(
r_{i}). Here,
v_{f}(
r_{i}) is the flow velocity at the position of the monomer
r_{i} evaluated onthefly from (time averaged) the binned
x_{2} coordinate.
Eqn (7) resembles the characteristic equation of the Berendsen thermostat,
^{82} where the linear differential equation in time is solved for the current temperature of the system. In our case, on the other hand,
T_{DPD} is just the temperature incorporated in the equations for the DPD thermostat and not the actual temperature of the system. The basic idea is simple: if the system under nonequilibrium sheared state is producing substantial heat due to the friction, the noise term (or in physical terms, the hypothetical reservoir temperature
T_{DPD}) of the thermostat should be made colder to rapidly extract heat. The adpd thermostat does not alter the equilibrium state (indeed we first checked this fact), but reduces the noise term under substantial shear.
The thermostat time τ_{DPD} (acting like a coefficient of heat transfer to the “reservoir”) was set to τ_{DPD} = 100δt. The thermostat nominal temperature T_{DPD} was updated using a simple explicit Euler scheme for eqn (7) with a time step of 100δt. All the thermostats used are applied to the relative velocities of monomers. The friction kernels of the DPD thermostats (damping and noise terms are constructed using the same kernel^{70}) are chosen to be Heaviside functions with a cutoff distance R^{DPD}_{cut}, i.e. γ(R) = 1 for R ≤ R^{DPD}_{cut} and zero otherwise.
All the results presented here correspond to a normal friction kernel γ_{∥} = 1 (some results for the adpd thermostat were also carried out for γ_{∥} = 5 to test sensitivity). To test the effect of tangential friction we also run simulations with γ_{∥} = γ_{⊥} = 1 in the tdpd thermostat. The thermostat details and labels used are given in Table 1.
The integration step ranges from δt = 0.01τ_{0} to δt = 0.005τ_{0} for the highest shear rates. Note that τ_{0} = 1 is smaller than the standard LennardJones time (monomer–monomer interaction), τ = σ(m/ε)^{1/2} = 2.415τ_{0} where m = 1 is the monomer mass.
The equilibration of the lump of melt in the OBMD simulation is conducted by a modified version of the AdResS, whereby the weighting function w in eqn (3) is gradually increased in time starting from w = 0 (CG model).^{83} The weighting function is therefore switched from a positiondependent to a timedependent one. The resolution is thus gradually sharpened from CG to AT. The procedure is described in detail in ref. 13. After equilibration, each simulation is run for 10000τ_{0}.
4 Melt at equilibrium
4.1 Characteristic times
In view of the close relationship between the structure and dynamics that takes place in sheared melts (and complex fluids in general), it is interesting to present the range of physical times of the melts before analysing their structural transformation with shear. In star polymers, one can observe three types of relaxation phenomena.^{84} First is elastic deformation of the overall shape of polymers, second relaxation occurs via the rotational diffusion, and the third one regards disentanglement of arms of every star polymer. Each of the relaxation processes can be estimated from the integral of the corresponding normalized autocorrelation function (ACF), via, and these are given by eqn (8)–(10), respectively. 
 (8) 

 (9) 

 (10) 
R_{i} represents the centerend vector of arm i, R_{i} its length, t time, and f the number of arms of each polymer. And i and j are indices of different arms within the same polymer. Each autocorrelation function decays with its characteristic time of the relaxation process^{84,85} and these are given in Table 2 for the different DPD thermostat friction kernels considered. We checked that under an equilibrium state all thermostats produce consistent results, in terms of pressure and density, while the correlations and the characteristic times differ, as they should. The adpd and the sdpdshort correspond to the same friction kernel and produce similar relaxation times at equilibrium (within statistical uncertainty) indicating that the adpd modification does not essentially alter the dynamics. We observe that the disentanglement of the arms occurs more rapidly than rotational diffusion and that is due to the short length of the arms, as each contains only 6 monomers. The longest time is the diffusion time for the center of mass (CoM) of molecules τ_{dif} = R_{g}^{2}/D, where R_{g} = 7.6σ_{0} is the radius of gyration of polymers and D the diffusion constant of CoM, which are different for every friction kernel. The differences in the relaxation times between the sdpdshort and sdpdlong thermostats are due to the larger thermostat cutoff distance R^{DPD}_{cut} leading to a larger interparticle friction in the sdpdlong case.^{70,86} Interestingly, the ratio between relaxation times and viscosity is similar for all the thermostats in Table 1 (values coincide within error bars), regardless of the kernel and R^{DPD}_{cut}.^{87} This indicates that the translational and orientational dynamics are affected in a similar way by the thermostats.
Table 2 Polymer time scales obtained at equilibrium for different setups and thermostats (see text). All times are in units of τ_{0} = τ/2.415, where τ = σ(m/ε)^{1/2} is the standard LennardJones time scale for monomers. In all standard and adaptive DPD cases, the thermostat damping constant is γ_{DPD} = 1m_{0}/τ_{0}, while for the tdpd thermostat γ = 1m_{0}/τ_{0} and γ_{⊥} = 1m_{0}/τ_{0}. τ_{elas}, τ_{rot}, and τ_{arm} are defined as characteristic decay times of the ACF given by eqn (8)–(10), respectively. Diffusion time is defined as τ_{dif} =R_{g}^{2}/D, with R_{g} = 7.6σ_{0} the average radius of gyration in equilibrium
Simulation 
τ
_{arm}

τ
_{elas}

τ
_{dif}

τ
_{rot}

Open sdpdshort 
33 ± 2 
6 ± 1 
800 ± 100 
55 ± 5 
Closed sdpdshort 
33 ± 2 
6 ± 1 
800 ± 100 
55 ± 5 

Open sdpdlong 
55 ± 2 
10 ± 1 
1200 ± 100 
100 ± 5 
Closed sdpdlong 
55 ± 2 
10 ± 1 
1100 ± 100 
100 ± 5 

Open adpd 
33 ± 2 
3 ± 1 
700 ± 100 
59 ± 5 
Closed adpd 
33 ± 2 
3 ± 1 
800 ± 100 
59 ± 5 

Open tdpd 
195 ± 5 
35 ± 5 
4800 ± 400 
332 ± 5 
Using the adpd thermostat (normal friction alone) we get τ_{elas} = 3 ± 1, τ_{arm} = 33 ± 2, τ_{rot} = 59 ± 5 and τ_{dif} = 700 ± 100, which illustrates the wide range in time scales involved in these sort of simulations. These times are similar in the open and closed systems (in equilibrium) and are compared with those obtained with other thermostats and kernels in Table 2. However, as shown below the rheology of the melt is not determined by the molecular diffusion, but rather by the molecular relaxation times. For the star molecule under study with f = 12 arms of length l_{a} = 6σ (i.e. 6 monomers per arm), molecular rotation is slower than the elastic relaxation of arms, the ratio being τ_{rot}/τ_{arm} > 10 in all the cases considered (see Table 2). This ratio determines the type of rheological behaviour of the melt according to a theoretical approach based on solving the Fokker–Planck equation for the bond distribution.^{88,89} We will come back to this issue in Section 5.3.
4.2 Equation of state
As explained in ref. 13, the OBMD simulations of the melt at equilibrium provide the correct average thermodynamic variables. The pressure equation of state obtained in OBMD agrees with that obtained from a standard NVT simulation for all the volume fractions studied.^{13} More precisely, the average equilibrium density 〈ρ〉 (p^{ext}_{22}) obtained with OBMD at a fixed normal pressure p^{ext}_{22} is consistent with the equilibrium pressure p calculated in a closed MD simulation at a fixed density ρ = 〈ρ〉 (p^{ext}_{22}). In passing, it should be highlighted that the same CG potential was used for all the melt densities considered, indicating that the pressure consistency at the MD domain is independent of the coarsegrained potential used in the buffer.
In the following, we work with the polymer volume fraction, defined as Φ = ρ_{N}πσ^{3}/6, where ρ_{N} = N/V is the number density (N and V respectively represent the number of monomers and the volume of the region of interest and the mass density is ρ = mρ_{N} with monomer mass m = 1, thus ρ = 0.136Φ). For the NVT ensemble, the study conducted here corresponds to p = (0.093 ± 0.001) and Φ = 0.20, fixed for any shear rate. In the open box, we fix p^{ext}_{22} = 0.093 and for zero shear rate (equilibrium) we get 〈Φ〉 = (0.20 ± 0.01). Fig. 3 shows the normalised density profile of polymers (ρ^{slab}_{M}/ρ^{avg}_{M}). Here ρ^{slab}_{M} denotes mass density of polymers in each slab, where it is measured, and ρ^{avg}_{M} its average value. The latter is constant in closed simulations, i.e. ρ^{avg}_{M} = 0.0271, while its calculated value in open cases equals ρ^{avg}_{M} = 0.0271 ± 0.0001. We observe that the obtained density profile is flat in the region of interest with some minor artefacts, which are due to the lack of statistics. Along the buffer zones the density gradually decreases as a consequence of the application of the external pressure p^{ext}_{22}. As it has been explained in previous related studies on hybrid atomisticcontinuum schemes,^{14,66} this inhomogeneity does not affect the transfer of pressure and stress from the exterior, provided that the density profile is flat around the hybrid interface. This is indeed the case, as it can be seen in Fig. 3 (see the interface between the “region of interest” and “buffer”). The radial distribution function (RDF) of CoMs of molecules is in perfect agreement with NVT simulations.

 Fig. 3 Normalised density profile (NDP – top panel) and comparison of RDFs (bottom). The NDP is depicted in the direction, in which the box is open. It is normalised to the desired value of the density, which corresponds to the occupational factor Φ = 0.2. RDFs are calculated from molecules inside the region of interest in open simulation (red) and closed systems (green).  
4.3 Mass fluctuations
In the grand canonical ensemble, mass fluctuations are related to the integral of the RDF,^{90} so the excellent agreement between the RDFs obtained from NVT and open boxes suggest that the fluctuations in the number of molecules inside the domain of interest should be thermodynamically consistent with the grand canonical prescription: namely, the relative fluctuation in polymer mass M = ρ_{M}V should be Std[M]/M = [k_{B}T/(Mc_{T}^{2})]^{1/2}, where c_{T} is the isothermal sound velocity c_{T}^{2} = (∂p/∂ρ)_{T}. c_{T}, which is evaluated using the pressure equation of state from ref. 13, is shown in Fig. 4. The agreement between OBMD and the NVT ensemble is excellent. In the same figure, we plot the results for the relative mass fluctuations at different pressures (average densities) and compare with the grand canonical prediction. We find that the agreement is very good, notably because of the tiny relative mass fluctuations in the (not small V = 3504384σ_{0}^{3} = 248805σ^{3}) volume considered for our open box: which range between 3.6% and 0.26% at the largest density considered (Φ = 0.2). In terms of mass density variance σ_{ρ}^{2} = ρk_{B}T/(Vc_{T}^{2}), for the state we considered hereafter under shear (p = 0.093ε_{0}/σ_{0}^{3}, Φ = 0.2 and ρ = 0.0271m_{0}/σ_{0}^{3}) the OBMD result is σ_{ρ}^{2} = (5.0 ± 0.3) × 10^{−9}m_{0}/σ_{0}^{3} in excellent agreement with the variance predicted for the μVT ensemble 4.79 × 10^{−9}. The conclusion of this study strongly supports our claim that the OBMD equilibrium simulation samples the grand canonical ensemble without any (or negligible) bias. It has to be said that the value of the external chemical potential μ^{ext} = μ(p^{ext}_{22},T) however cannot be imposed in OBMD, although it could be reconstructed following the standard Gibbs–Duhem route with varying external pressure. New implementations of the AdResS^{29,65} might also be used to evaluate μ.

 Fig. 4 Top panel: The relative fluctuation of the polymer mass (standard deviation over average polymer mass M) in the interest MDdomain of the open setup (see Fig. 2). Comparison is made with the grand canonical theoretical result (green dashed line) Std[M]/M = [k_{B}T/(Mc_{T}^{2})]^{1/2} and, also included, the ideal gas limit (red dashed) obtained with the isothermal sound velocity c^{id}_{T} = (k_{B}T/M_{m})^{1/2} with M_{m} = 73m_{0} the molecular mass. The volume of the interest domain is V = 156 × 117 × 117σ_{0}^{3}. Recall that the monomer LJ diameter is σ = 2.415σ_{0} and its mass is m_{0} = 1. The right ordinate axis shows the values of the isothermal sound velocity c_{T} = (∂p/∂ρ)_{T} (from the pressure equation of state) that is compared with ω/k_{eff} obtained from the oscillation frequency ω of the total mass autocorrelation function in the MD domain (see bottom panel). The effective wavenumber is k_{eff} = (π/370)σ_{0}^{−1} and the total open box is L_{2} = 390σ_{0}. Bottom panel: the time autocorrelation function of the mass in the MD domain at equilibrium with imposed external pressure p^{ext}_{22} = 0.001ε_{0}/σ_{0}^{3} (dashed line is the fit to extract ω, see text); the inset illustrates the evolution of the total mass in the MD domain for this case.  
A typical outcome for the time evolution of the total mass of polymers in the MD domain is shown in Fig. 4. It presents oscillations, suggesting that it might contain some information about the sound velocity of the system.
Sound propagation can indeed be studied in our simulations because we work with momentum preserving (DPD) dynamics (by contrast, sound is damped in the standard Langevin dynamics). It is noted that an “ideal” open system should be transparent to all waves, meaning that all waves, either created by inner mass fluctuations of any wavelength or by external waves travelling across, should leave the system and do not reflect back. This implies, in particular, that in the absence of external longitudinal forces, fluctuations of the total mass should have no memory, being a white noise (or at least a broadband signal).^{91} The presence of correlations in the total mass of our “region of interest” is in fact due to partial reflections of waves at the rarefied buffers, where density fluctuations are reduced (recall we fix the average total mass of these reservoirs, see e.g.ref. 66). In the autocorrelation (ACF) of the total mass fluctuation, the dominant wavelength should obviously be the largest possible one compatible with this condition and the total simulation box length. A density mode with wavenumber k = 2π/λ decays like 〈ρ(k,t)ρ(k,0)〉 = 〈δρ^{2}〉exp[−ν_{L}k^{2}t]cos[ωt] with ν_{L} being the sound attenuation coefficient^{92} and ω = c_{T}k the oscillation frequency. Inspection of Fig. 3 indicates that the density profile of the whole system (MD + buffer) roughly conforms to Dirichlet boundary conditions with a fixed density at the end of the buffer domains ρ(x_{2} = ±L_{2}/2)<ρ_{0}. In such a case, the longest wavelength available to mass fluctuations of the system would be λ_{eff} ≃ 2L_{2}. We fitted the time autocorrelation function of the MD mass to extract ω and compare it with the ansatz ω = c_{T}k_{eff} using ω = c_{T}k_{eff} with k_{eff} = 2π/λ_{eff}. The best fit to the simulation results corresponds to λ_{eff} = 760σ_{0}, while 2L_{2} = 780σ_{0} and it is shown in Fig. 4 in terms of ω/k_{eff} and compared with the sound velocity of melt c_{T}. The excellent match confirms that the mass in the MD domain has memory induced by reflections of sound waves against the low density domains (buffers). Such reflections could be reduced by coupling the MD with a continuum hydrodynamic field outside,^{66} or by imposing a nonreflecting boundary condition^{91} (still to be generalized to MD, see also ref. 92). However, in the present scenario we find that these results are quite interesting because they suggest the possibility of measuring the sound velocity c_{T} from the fluctuations of the total MD mass. In particular, it might allow to measure how the sound velocity is modified in a sheared melt, c_{T}(Wi). Although a detailed study would be certainly required, just as an indication, we find that c_{T}(Wi) (estimated in this way, i.e. from mass fluctuations) decreases.
5 Sheared melt with normal friction
This section presents results for the star molecule model with zero tangential friction between monomers. Kernel and thermostat details (adpd) are given in Sections 3 and 4. We decided to first focus on this model to avoid embarrassing the discussions with details of different cases (open versus closed, normal versus tangential friction) and also because this model presents a richer dynamical behaviour, whose analysis will be useful to understand the effect of tangential friction in Section 8.
5.1 Weissenberg number
The Weissenberg number Wi is a useful number to compare and dissect different regimes in polymer rheology. It is defined as 
Wi = τ_{rel},  (11) 
where τ_{rel} is the longest molecular relaxation time and ^{−1} is the “shear flow time” needed to affinely deform a square box of sheared fluids into a parallelepiped with an angle of 45° between adjacent planes. As stated, in eqn (9) and Table 2, the longest relaxation of our star polymer is related to the molecular rotation so τ_{rel} = τ_{rot} in eqn (11). In fact, the diffusion of the CoM of the molecules is much slower (see Table 2). The CoM diffusion does not directly sample the structure of molecules, whose modification under flow is related to the nonNewtonian character of polymers. The Peclet number determines the ratio between CoM diffusion and flow advection Pe = τ_{dif} and for our setup it is about 10 times larger than Wi. In colloidal suspensions, the shear thinning typically starts for Pe > 1 due to shear banding. Interestingly, star polymers constitute a sort of bridge between the open structures of linear polymers and the compactness of colloids. Thus, one might elucubrate that the onset of shear thinning in compact stars could well be due to shear banding (collective molecular ordering in lanes), rather than by (individual) polymer elongations. We shall see later on that both (collective and individual) effects take place in our sheared system, although we advance that the transition to the nonNewtonian regime takes place at Wi = τ_{rot} > 1. Hence at least for the (not so compact) star polymer studied here, shear thinning is not determined by collective ordering at straining rates faster than the CoM diffusion. The hybrid character of star molecules (between colloids and polymers) is the subject of current research.^{56–58}
5.2 Density and hydrostatic pressure
Denoting the symmetric pressure tensor (exerted by the melt) as P_{α,β}, the hydrostatic pressure is defined as 
 (12) 
We find that our model for a star polymer melt expands when sheared under the normal load (see Fig. 17 which collects the results from different cases). This behaviour is consistent with the die swelling phenomenon of polymer melts^{1} and the increase in the normal pressure P_{22} observed in all cases. By contrast, the hydrostatic pressure decreases with Wi in the open domain, while in the NVT box presents a nonmonotonous trend (with not large variations). As explained in Section 6 P_{iso} is influenced by several molecular mechanisms, whose relevance changes with along with the molecular structure. The density (in the open domain) and pressures (P_{22} in NVT and also P_{iso} in both cases) are quite sensible to temperature changes, as shown in Section 9.
5.3 Rheology
We start by determining the zero shear viscosity, η_{0} = η( = 0). A standard way is the Green–Kubo equilibrium route (via the integral of the time autocorrelation of shear stress).^{93} To maintain consistency with the nonequilibrium route we follow the experimental approach which fits the trend for η() for a range of values of the shear rate with some suitable expression such as the popular Carreau fit.^{50,94} We also checked that Green–Kubo viscosities^{95} are similar to the Carreaufitted ones within statistical uncertainty (about 10%). The Carreau fit has the following expression, 
η = η_{0}[1 + (τ_{η})^{2}]^{−βη/2},  (13) 
shown in Fig. 5 along with simulation data. Fig. 5 shows the shear viscosity η obtained under open and closed setups for several models with different thermostats and friction kernels (see Table 1). η() was calculated from the offdiagonal pressure tensor component P_{12} = −η, which was measured in simulations. From eqn (13), the viscosity shear thinning exponent β_{η}, is such that,
η → ^{−βη} for large 
and for polymer melts β_{η} ranges between 0.4 and 1.^{1,88} This fit of eqn (13) also provides an estimation of the zeroshear viscosity η_{0} and a characteristic time τ_{η} related to the onset of the shear thinning regime. η_{0} differs for different friction kernels of the DPD thermostat (see Table 1). The adpd and sdpdshort thermostats (identical friction kernels) consistently provide the same zero shear rate viscosity η_{0} = 0.5 (that does not increase largely for γ_{∥} ≤ 10). The sdpdlong model with an increased kernel cutoff R^{DPD}_{cut} = 1.5 × 2^{1/6}σ presents η_{0} = (0.60 ± 0.1) + 0.29γ_{∥} an increase consistent with the increasing relaxation time with friction† as deduced in the analysis of Kindt and Briels.^{96} Tangential friction (tdpd) leads to η_{0} = 2.6 (γ_{∥} = γ_{⊥} = 1). The meaning of τ_{η} becomes clear when it is compared with the estimated star rotational relaxation time τ_{rot}. Here are their values for different thermostats: τ_{η} = 61 and τ_{rot} = 59 for adpd; τ_{η} = 58 and τ_{rot} = 55 for sdpdshort; τ_{η} = 125 and τ_{rot} = 100 for sdpdlong; and τ_{η} = 287, τ_{rot} = 332 for tdpd. The error bar of the given values is approximately ± 5. We observe τ_{η} ≃ τ_{rot} indicating that the onset of shear thinning, which takes place at τ_{rot} > 1, coincides with the molecular deformation altering the equilibrium rotational diffusion.

 Fig. 5 Normalized shear viscosity obtained for several models under open and closed conditions. Details of the model (varying in thermostats and kernels) are given in Table 1 and 2. Lines are best Carreau fits.  
We now focus on the adpd model and defer the discussion on the tangential friction to Section 8. The system temperature is fixed to T = 4.0 ± 0.01. In this case, we get a zero shear viscosity η_{0} = 0.50 ± 0.05. At larger Wi, the shear thinning exponents β obtained from Carreau fits (eqn (13)) are found to be slightly steeper in the open domain β_{η} = 0.41(2) than in the NVT box β_{η} = 0.35(4). Thus at a fixed shear rate, the open system is slightly less viscous than the closed sample. This is in agreement with previous studies for linear and branched melts carried out at a constant and unconstrained density (see Section 10).
The first and second normal stress coefficients Ψ_{1} = N_{1}/^{2} and Ψ_{2} = N_{2}/^{2} (for the adpd model, discussed in this section) are shown in Fig. 6. For Wi < 20 we find a decrease in Ψ_{1} consistent with the Carreau standard behaviour (eqn (13)), providing Ψ_{1}( = 0) = 21 ± 1 and an exponent β_{Ψ1} ≃ 1.30 ± 0.04 (i.e. N_{1} ∼ ^{0.7}) quite similar for both ensembles. The relaxation time for Ψ_{1} obtained from the fit is also consistent with τ_{η} = τ_{rot} within error bars. At Wi > 20, we find a measurable decrease of Ψ_{1} with respect the Carreau trend (see Fig. 6), which takes place at slightly smaller Wi in the open case. This corresponds to a loss of the elastic component of the melt at large shear rates. The second normal stress coefficient Ψ_{2} = N_{2}/^{2}, shown in Fig. 6 is also quite similar in both environments and at large scales like Ψ_{2} ∼ −^{−1}. The similar behaviour for N_{2} under open and closed boxes might be due to the fact that we just fix the normal load (in the x_{2} direction) and not the hydrostatic pressure as in some other studies,^{39,40,42} presenting different trends for N_{2} in NVT and NP_{iso}T constraints.

 Fig. 6 First and second normal stress coefficients of the star polymer melt under open and closed conditions and the adpd thermostat for γ_{∥} = 1 (blue) and γ_{∥} = 5 (purple). Dashed line for Ψ_{1} is the Carreau fit (see text).  
The normal stress ratio VR ≡ −Ψ_{2}/Ψ_{1} and the recoverable shear strain^{87}

SR ≡ Ψ_{1}/(2η) = (P_{22} − P_{11})/(2P_{12})  (14) 
are standard indicators of viscoelasticity
^{87} (
e.g. SR vanishes for a Newtonian fluid). As shown in
Fig. 7 in our melt model, SR increases with the shear rate, as expected; however, both indicators (SR and VR) clearly show that a change in the elastic component of the model takes place at Wi > 20. Notably, SR decreases, so the melt becomes less compliant with shear strain and stores less elastic energy across the flow–gradient plane, with a jump in the VR. We advance that the amount of the elastic energy loss at large Wi depends on friction forces (notably on the presence of tangential friction as shown in Section 8). This is a clear indication that calibration of friction from detailed allatom models
^{9} is crucial to represent or simulate some particular real melt.

 Fig. 7 (top pannel) Normal stress ratio and (bottom pannel) recoverable shear strain index for some model melts under open and closed environments. The results for the case with tangential friction (tdpd) are also included.  
5.4 Gyration tensor
The thermodynamic and rheological properties of any polymeric system are intimately related to the deformation of the molecules induced by flow. It is therefore convenient to start by presenting the results for the average gyration tensor of our star polymer model under shear, whose components are shown in Fig. 8. At Wi > 1 the flow induces the alignment of the star molecules in the flow direction and at the same time a compression in the gradient direction. The width of the stars in the neutral direction slightly increases up to Wi ∼ 20 and AT larger shear rates, they also start to contract in this direction. It is also instructive to observe the shape of molecules in their principal deformation axes, obtained from the diagonalization of the gyration tensor. The molecular axes are s_{1} = cos(θ_{G})x_{1} − sin(θ_{G})x_{2}, s_{2} = −sin(θ_{G})x_{1} + cos(θ_{G})x_{2}, and s_{3} = x_{3} (since G_{13} = G_{23} = 0). The ith eigenvalue is noted as _{i}. The angle θ_{G} is the average molecular tilt over the flow direction that satisfies: 
 (15) 

 Fig. 8 Eigenvalues of the gyration tensor scaled with the monomer diameter σ^{2}. The results include cases with the adpd friction kernel γ_{∥} = 1 (no tangential friction) and the tdpd model γ_{∥} = γ_{⊥} = 1 discussed in Section 8. Estimated error bars of the data are approximately 5%.  
The results for the eigenvalues of gyration tensor are shown in Fig. 8 and reveal the strong contraction of the molecules in the s_{2} direction. In particular, at Wi ∼ 50 we observe that _{2}^{1/2} ∼ 1.12σ so the chain width in the contracted axis reaches the monomer (or blob) diameter. In a real polymer at this shear rate, the flow starts altering the structure of the chain “blobs” smaller than σ. But this is precisely our modelling unit length so we will take this as the maximum Wi to be explored here (roughly Wi < 100). Above this (size dependent) shear rate, one should also expect a shear rate dependence of the blob–blob friction kernel (see ref. 97 for a related study on this issue). It is also interesting to note that the total volume of the molecules (evaluated from the product of the eigenvalues of the gyration tensor) first increases (molecules expand) up to Wi ∼ 10 and due to contraction in the neutral direction and then tends to contract (in the form of a highly extended ellipsoid) at larger Wi. In this sense, the behaviour of the arms in the neutral direction of the star have some similarities with what is found in tethered polymers,^{98,99} which also present a maximum volume at intermediate Wi.
6 Analysis of the stress–structure and dynamics relationship
The present analysis provides the relationship between stress components and the molecular structure and dynamics. It focuses on the results presented in the previous section, which corresponds to the adpd star model, with zero tangential friction. The conclusions will be useful for the comparison done in Section 8.
6.1 Pressure balance in closed system
Pressure is the leading mechanical variable in rheology as it directly connects the microscopic world with material properties. We now elaborate on the pressure balance in an attempt to connect the molecular structure with the observed rheology.
The pressure exerted by the melt‡ is calculated using the IrvingKirkwood method and the method of planes,^{100} which is particularly suitable for the open setup as it was designed for non periodic boundaries. The pressure can be first decomposed in virial and kinetic parts. The virial pressure includes contributions from spring forces and intramolecular forces (both acting amongst pairs of monomers of the same molecule) and intermolecular forces (between monomers of different molecules). The pressure balance allows to analyze the different molecular contributions,

P_{α,β} = P^{kin}_{α,β} + P^{spring}_{α,β} + P^{intra}_{α,β} + P^{inter}_{α,β}.  (16) 
The pressure tensor is symmetric for the type of molecules considered hereby (in fact for most polymers).
^{1} The “average” pressure on the system is given by the hydrostatic pressure which is just the third part of the trace of
P (see
eqn (12)). The hydrostatic pressure is not involved in (although indirectly affects) the main rheological quantities such as the shear stress
P_{12} and the first and second normal stress differences, defined as

N_{1} = P_{22} − P_{11},  (17) 

N_{2} = P_{33} − P_{22}.  (18) 
It is customary to introduce the traceless pressure tensor
≡ P − P_{iso}I, 
with I = δ_{α,β} the unit second rank tensor. Indeed, the hydrostatic pressure can be decomposed in different molecular contributions with A = kinetic, springs, etc. Also, due to the linearity of the trace operator we can decompose the traceless stress tensor P in a sum of contributions of traceless tensors, i.e., with ^{A} = ^{A} − (1/3)Tr[P^{A}]I. This is useful to simplify the analysis. In fact, before discussing the differences between open and closed ensembles, we inspect some generalities of the pressure balance focusing on the results obtained for the closed (NVT) ensemble. To that end, we present in Fig. 9 the contributions to the pressure components of the melt under shear to the hydrostatic pressure and the different traceless stress tensors (componentwise: gradient _{22}, flow _{11}, neutral _{33} and shear P_{12} obviously P_{13} = P_{23} = 0).

 Fig. 9 Components of the different contributions to the traceless tensor (, with A= kinetic, springs, etc.) of the traceless stress tensors ^{A} = P^{A} − (1/3)Tr[P^{A}]. The results correspond to the closed setup using the adpd thermostat with γ_{DPD} = 5. P_{0} is the (isotropic) pressure at = 0, the rest of terms are explained in the text. The error bars of the data are around 0.002–0.005ε/σ_{0}^{3}.  
6.1.1 Spring stress.
Bonded interactions are hereby modelled by linear springs with a non zero equilibrium distance. The contribution of the springs of star molecules to the traceless stress tensor is illustrated in Fig. 9 (again, for the case of closed systems). As expected, as the shear rate is increased, molecules tend to be stretched in the flow direction and compressed in the gradient direction. Springs tend to restore the equilibrium (spherical) shape of the molecules by producing a negative stress (compressing) in the x_{1} direction and a positive (expanding) in the gradient direction. As in any material with elastic properties, both effects contribute to enlarge the first normal stress difference N_{1} = P_{22} − P_{11}. In turn, we find that elastic contribution to the second normal stress difference N_{2} = P_{33} − P_{22} is negative, as it happens in real polymeric fluids.^{1}Fig. 9 shows that the offdiagonal component of springs' contribution to the pressure tensor P^{spring}_{12} is the most important contribution to the shear stress, although as discussed below, the kinetic pressure becomes significant at large Wi (above Wi > 20 for the adpd model discussed in this section).
6.1.2 Kinetic pressure.
The kinetic pressure tensor was obtained from the peculiar velocities δv_{i} = v_{i} − u(r_{i}), where u(r) = u_{1}(x_{2})e_{1} is the average velocity field obtained by binning the velocity gradient direction x_{2}. Indeed, the ideal part of the hydrostatic pressure is one third of the trace of the kinetic pressure tensor, (1/3)Tr[P^{kin}] = ρT. However, its traceless part, ^{kin} ≡ P^{kin} − ρTI, contains relevant information about the correlations in velocity fluctuations over different axes. In particular, velocity fluctuations of monomers in the x_{1} and x_{2} directions become correlated at > 0 and contribute to the shear stress and viscosity. In our setup, this correlation becomes negative (owing to the imposed counterclockwise shear) and leads to a negative contribution to P^{kin}_{12} with a nonnegligible viscous contribution to the shear viscosity η = −P_{12}/ which tends to reduce the shear thinning exponent (see the comparison with the case with tangential friction in Section 8). Although seldom pointed out in the literature (see ref. 4), at high enough shear rates the kinetic pressure of monomers becomes significant in polymer melts. This is in fact what we clearly observe in Fig. 9 for the star model with normal friction (adpd). The anisotropy in the kinetic pressure observed in the covariance of peculiar velocities of monomers is certainly much larger than what would be observed in a simple fluid.^{47,78} It is also of a different nature because in the melt, the monomer motions are constrained and they fluctuate and rotate being tethered to the center of their star molecule. This leads to different monomer peculiar velocities in the flow, gradient, and neutral directions which are reflected in the kinetic pressure. In Fig. 9, we show ^{kin}versus Wi (recall that P^{kin} = ^{kin} + ρTI and diagonal components of P^{kin} are positive). At Wi > 1, the kinetic pressure contributes to the first and second normal stress N^{kin}_{1} and N^{kin}_{2}. As soon as the molecules become stretched in the flow direction (for Wi > 1), velocity fluctuations in the flow direction are enhanced with respect to those along the gradient and neutral directions: ^{kin}_{11} > 0 and ^{kin}_{22} < 0, ^{kin}_{33} < 0. We find that N^{kin}_{2} < 0 (see Fig. 10 below), so the kinetic stress acts in the same way as the elastic (springs) components. However, the kinetic contribution to the first normal stress differences is negative N^{kin}_{1}<0 so it is just opposite of the elastic contribution of the chain. Our conclusion is that the kinetic pressure of monomers tends to reduce the elasticity of the melt in the flow–gradient direction at large shear rates.

 Fig. 10 Contributions of different mechanisms to the summing up the total normal stress differences, shear stress, and hydrostatic pressure. The results for the adpd thermostat under open and closed conditions. The error bars of the data in graphs are around 0.002–0.005ε/σ_{0}^{3}.  
6.1.3 Intramolecular pressure.
Intramolecular pressure gives an indication of excluded volume effects and molecular collisions within one molecule. For our moderatesize molecules, it has a minor contribution to the total stress. Intramolecular collisions induce a viscous (notrestoring) stress in the melt which slightly contributes to the shear stress and tends to counterbalance the elastic first normal stress difference N^{intra}_{1} < 0 (just like the kinetic pressure – see Fig. 10). It has the same effect in the gradient–neutral plane (N^{intra}_{2} > 0) where due to the molecule stretching, monomers tend to collide less in the gradient direction. However, above Wi ∼ 20 a shallow maximum is observed in P^{intra}_{33} (and a minimum for P^{intra}_{22} – see Fig. 9) suggesting that the reduction of the arm extension in the neutral direction has the consequence of increasing monomer collisions in the gradient direction (thus N^{intra}_{2} diminishes – see Fig. 10). Evidence of arm retraction in the neutral direction is shown in Fig. 8 as a clear the reduction of _{3} = G_{33} for Wi > 20 (also observed in “sdpdshort” and “sdpdlong” thermostats). A somewhat similar rebound in the intramolecular pressure at a large shear rate has been observed in simulations of entangled and disentangled linear polymer melts (see ref. 3 and reference therein). However, in these studies the monomer–monomer interaction is attractive, and this fact modifies (reverses) the intramolecular contributions (polymer compression then reduces the intramolecular pressure, unlike for our purely repulsive potential).
6.1.4 Intermolecular pressure.
According to Fig. 9, the contribution of the intermolecular pressure to the melt is minor. The seemingly irrelevant contribution of intermolecular stress is assumed in many theoretical models^{87} to explain viscoelasticity (stressoptic rule). However, this is a simplistic view because intermolecular forces are responsible for spreading the external momentum introduced through the system boundaries (note that, in this respect, the open boundary setup behaves like a real experiment). In fact, internal forces (between monomers of the same molecule) sum up to zero so they cannot modify the CoM velocity of the molecule. For the present star molecule (with relatively short arms), momentum in the flow direction is transferred and maintained across x_{2} by intermolecular collisions.^{9} These friction forces gradually build up the elastic stress in the molecules, before finally collapsing to a stationary value in the steady state. The central role of intermolecular forces can be also seen by considering the alternative molecular formulation of the pressure tensor^{101} (based on molecules (μ) CoM and virial pressure proportional to ). This “molecular pressure” formulation is however less informative because the effects of all internal forces, like the springs, are hidden in the spatial CoM distribution g({R}). Nevertheless, it serves to illustrate the central role of intermolecular forces and shed some light on the apparently striking similarity between the intermolecular pressure and the total pressure dependence with Wi, which can be seen by comparison between the corresponding (inter and total) panels of Fig. 9 (note the difference in values). This similarity between intermolecular (IM) and global magnitudes (such as IM energy and hydrostatic pressure, see Fig. 9) has also been reported as “striking” in previous sheared melt simulations.^{2,3} It is interesting to note that the dominant velocity gradient component of the IM force P^{inter}_{22} reaches a plateau around Wi ∼ 20 and slightly decreases at larger Wi (as the total P_{22} does – see Fig. 9). This is indicative of a change in the dynamics of polymers, which according to the concomitant increase observed in P^{inter}_{11}, most probably start to rotate and collide more often in the flow direction.
6.1.5 Hydrostatic pressure.
The hydrostatic pressure is key in shear induced polymeric phenomena, such as shear induced crystallization^{87} or separation of blends.^{47} Its dependence on the shear rate is not well understood. P_{iso} depends on the molecule size^{40,43} and their architecture^{3,40,42,88}. However, there have been reports of both increase and decrease in different (sometimes contradicting) studies. Fig. 9 indicates that P_{iso} can present a nonmonotonous trend with Wi due to changes in the molecular structure under shear. Nonmonotonous trends for P_{iso}() have also been observed for different polymers in ref. 2, 3, 40 and 43. The present analysis (see Fig. 9) reveals in fact that P_{iso} depends on a competition of several mechanisms. Below Wi ≃ 10 the hydrostatic pressure varies little but presents a slight descend, probably due to the chain expansion (revealed by the analysis of the gyration tensor) which decreases the intramolecular pressure. At larger shear rates Wi > 10, two opposite mechanisms enter into the play: First, an increase in P_{11} due to the strong increase of kinetic pressure (and to a lesser extent to intermolecular collisions). This is probably due to molecular rotations similar to the tankthread motion reported for a star molecule solution^{102} and also pointed out in ref. 2 and 3. And second, a decrease in P_{33} due to the contraction of the stars in the neutral direction (see Fig. 8) and consequent reduction of the kinetic pressure P^{kin}_{33}. Both effects nearly counterbalance each other (see Fig. 9) leaving small variations in P_{iso}. The relevance of these effects depends on the boundary conditions (the open case is analysed in next section) and type of friction (Section 8). More generally variations of P_{iso} depend on the presence of attractive monomer interactions and the molecular structure.
It is interesting to note that the decrease of neutral kinetic pressure in our adpd star model starts to take place around Wi ≃ 20, which is precisely the ratio τ_{rot}/τ_{elas} ≃ 20. Above this shear rate the flow strains faster than the elastic relaxation of the molecules thus reducing the fluctuations of the arms in the neutral direction. This effect finally induces a net decrease in P_{iso} above Wi > 30. We shall see (Section 10) that this effect is absent (or at least delayed) while introducing the tangential friction, again indicating that friction should be an essential part of any CG model.^{9}
6.1.6 Intermolecular forces and the Hookean limit.
To summarize, at low shear rates the elastic energy stored by the melt grows in response to intramolecular (nonbonded) interactions (here mainly friction forces). At large enough shear rates viscous (Newtonian) effects coming from the kinetic and intra molecular pressures, tend to modify (normally reduce) the elastic response of the melt (notably first normal stress differences). Our findings are in agreement with the conclusions of Kroger et al.^{4} clearly and succinctly summarized in their book (page 144) in relationship with the breakdown of the linear stressopticresponse (SOR) due to the Newtonian viscous transport at large Wi. In our simulations we find a linear relationship between G_{22} − G_{11} and P_{22} − P_{11} (essentially similar in open and closed setups), indicating the validity of the linear SOR up to Wi ≃ 10. In the approach by Kroger et al. (applied to linear multibead FENE chains), both, kinetic and intramolecular pressures, are collected in what they call the “simple fluid” stress contribution.^{4,89} We find here that both contributions can have different roles, which more generally should probably depend on the molecular shape (kinetic pressure) and monomer (intramolecular) interactions (attraction/repulsion).
To highlight the relevancy of the intermolecular stress in the linear SOR regime, we plot in Fig. 11 (top panel) the total and elastic contributions (springs) to the normal and shear stress against the intermolecular counterparts. At low shear rates, about Wi < 10, the same linear relation is found for the first normal stress difference and shear stress, N_{1} ≃ 8N^{inter}_{1} and P_{12} ≃ 8P^{inter}_{12}, while we found N_{2} ≃ 4N^{inter}_{1}. Approximately the same linear relation holds for the elastic stress and also for the normal molecular strains evaluated with the gyration tensor (G_{11} − G_{22} and G_{22} − G_{33}, scaled in Fig. 11(top panel)). This provides the Hookean limit of the melt N_{1} ≃ C(G_{11} − G_{22}) and N_{2} ≃ (C/2)(G_{33} − G_{22}) with C = 62.5 (we find that P_{12} ∼ CG_{12} holds only for smaller Wi < 5). In the open (adpd) setup, the linear regime for N_{1} and P_{12} perfectly agrees with closed simulations; however we found deviation from linearity in the case of N_{2}, an issue which deserves further investigation. The bottom panel of Fig. 11 presents the results for stars with the tangential friction, analysed in Section 10.

 Fig. 11 Total stress, elastic stress components, and molecular strain differences in normal planes measured from the average gyration tensor (G_{11} − G_{22} and G_{33} − G_{22}) plotted against the intermolecular stress component corresponding to each case (i.e., N^{intra}_{1}, N^{intra}_{2}, etc.). The results for adpd (top panel) and tdpd models under open conditions (bottom panel). Estimated error bars of the depicted quantities are approximately 5%.  
7 Effect of open environment
7.1 Density and hydrostatic pressure
In our open domain, we fix the load of the melt in the gradient direction (P^{ext}_{22} = P_{22}) and this produces a redistribution of the pressure tensor, reducing its component in the flow and neutral directions and also, indirectly, its shear stress. This is deduced from Fig. 12 where we compare the traceless stress tensor and the hydrostatic pressure in the open and closed setups at the fixed T = 4 temperature. Note that P^{ext}_{22} = P_{22} = _{22} + (1/3)P_{iso} is constant in the open setup (within statistical uncertainties). In the open domain, the sheared melt expands in the gradient direction; a phenomenon similar to the die swelling observed in polymer extrusion at the pipe orifice and related to other viscoelastic phenomena.^{103} In the open domain, this corresponds to a decreasing melt density (at faster shear rates) and brings about a smaller hydrostatic pressure than in the closed environment at similar Wi (see Fig. 12). However, the relative decrease of P_{iso} is larger than the density jump. This fact is due to several related effects we now analyze from the inspection of Fig. 10. Indeed, at the fixed temperature, a lower density leads to a lower kinetic pressure P^{kin}_{iso} = ρT found in the open domain (this trend also applies here to the intramolecular pressure P^{intra}_{iso} because our model considers purely repulsive nonbonded interactions, the opposite effect could arise for attracting chains^{2,3}). However, an even larger reduction in P_{iso} with respect to the closed box is observed due to the smaller intermolecular pressure in the open box (see Fig. 10 for Wi > 10). Indeed, at a high shear, a less dense melt presents less molecular collisions, less intermolecular friction and thus less elastic load. As stated, at the fixed temperature, the elastic strain is essentially activated by intermolecular friction in the melt. Notably, at Wi ∼ 20, these intertwined effects induce a reduction of about 25% of hydrostatic pressure of the open domain mainly arising from the decrease of elastic stress. At this Wi, density has only decreased about 10% (see Fig. 17). In agreement with this comment, we note that this depressurizing effect is doubled when the tangential friction is added as commented in Section 10.

 Fig. 12 Hydrostatic pressure and components of the traceless stress obtained in closed and open domains at fixed T = 4 and under the adpd thermostat (normal friction). Error bars of the measured quantities are approximately 0.005ε/σ_{0}^{3}.  
7.2 Normal stress differences and shear stress
Fig. 10 compares the different contributions to the normal stress differences and the shear rate in the open and closed systems. Much of what has been already said in the previous section applies here. The results for the shear stress nicely corroborate what we pointed out before about the close relationship between intermolecular and elastic stresses. Fig. 10 (right bottom) clearly shows that the kinetic and intra molecular shear stresses are essentially equal in the open and closed domains. The decrease in elastic shear stress found in the open case is due to the reduction in the intermolecular friction at a lower density (although seemingly paradoxical, the “nominal” contribution to P^{intra}_{12} in the monomer pressure balance is minor). In all instances, at low enough shear (Wi < 10) the elastic stress is close to the total stress (and proportional to the intramolecular stress). This is the regime of validity of the stress optic rule which is broken at larger shear due to viscous (and compressible) effects related to the “simple fluid” of monomers. Remarkably, for the present model of the star polymer melt, the kinetic normal stress becomes the dominant “viscous” contribution and at large Wi > 20 it even induces a decrease in the first normal stress difference N_{1}.
7.3 Molecular ordering under shear
Fig. 13 presents the angle of the largest eigenvector of different contributions to the pressure tensor with the flow direction (measured according to eqn (15)). We also include the molecular orientation, measured from the angle associated with the gyration tensor (eqn (15)). This plot condenses what has been already mentioned in previous sections. Recall that a spherical molecular structure provides θ_{G} = 45° and similarly from eqn (15) a Newtonian fluid without elastic component, necessarily presents θ = 45°. This is what is observed for the total pressure tensor angle θ_{P} and the molecular orientation θ_{G} at Wi → 0 in Fig. 13. As the straining rate is made faster both angles decrease in a similar fashion, however, at Wi > 10, the pressure tensor angle θ_{P} presents a minimum and starts increasing towards 45°. By contrast, the molecular orientation θ_{G} keeps aligning with the flow direction. This in an indication of the loss of the Hookean behaviour of the melt which here is mainly due to the kinetic pressure (see its principal direction in Fig. 13). The stress direction of springs also aligns with the flow, although their angle is larger than the molecular orientation (a similar outcome was observed in ref. 3). Finally, note the close match between the direction of intermolecular forces and the total pressure tensor. As stated, intermolecular forces are the driving mechanism of transformation between the viscous flow and elastic energy. Lastly, as observed by other authors,^{39,42,43} the open boundary does not modify the molecular structure or orientation with Wi, when compared with the closed case. Here, Fig. 13 shows another remarkable result: the significant redistribution of pressure in the open case (see Fig. 12) does not alter the orientation of the different contributions of the pressure tensor at the increasing shear rate. The orientation of pressure eigenvectors is a function of their eigenvalues which in turn determine the material properties of the polymer (its viscoelasticity^{87}). A relevant example is the recoverable shear strain (SR) given by eqn (14). Material properties should not depend on the constraints used to perturb the polymer and this is precisely what our analyses provide.

 Fig. 13 Top panel: angle between the principal eigenvectors of the contributions to the stress tensor and the principal eigenvector of the average gyration tensor. Bottom panel: comparison with the model including transversal friction (tdpd model). Angles are with the same color coding.  
To observe the collective order of the star molecular, we calculated the CoM pair distribution function g(R_{ij}). Fig. 14 illustrates the marginal distributions for different planes and at increasing shear rates. It is illustrative to draw the directions of the principal components of the pressure and the gyration tensor to observe the departure from the Hookean (linear SOR) regime. Above Wi > 20 molecules start to orient in lanes in the flow–gradient plane, as indicated by the elongated shape of the CoM distribution in Fig. 14. In linear polymers at large shear rates, this effect creates order and even crystallization (ref. 87). This collective order is also clearly visible in one snapshot of the system, in Fig. 15. Molecules with different relative velocities (closeby in the gradient direction) slide over creating a stress which is larger along a direction differing from the orientation of individual molecules. This direction of maximal stress is correlated with the CoM distribution (see Fig. 14) which shows bright spots where molecules slide over and depleted regions in the “wake” of each molecule. At the largest shear rates considered, we also observe some collective ordering in the neutral (vorticity) direction. In sheared colloids, lanes of particles in the neutral direction appear due to hydrodynamic interactions.^{104} This could be a plausible hypothesis, in view of the hybrid colloid–polymer nature of the star molecule direction.

 Fig. 14 The pair distribution of molecules, g(R) providing the probability of finding the CoM of a molecule at a distance (vector) R from the target molecule CoM. Left: Marginal probability in the flow–gradient and (right) gradient–neutral planes. The results for increasing Wi in a closed domain (open simulations at similar Wi are visually indistinguishable). Green line denotes the direction of the deviatoric stress and the blue line the molecular orientation obtained from the gyration tensor eqn (15).  

 Fig. 15 Snapshot of a closed box simulation under shear rate (Wi = 42). For the sake clarity only 10% of randomly chosen molecules are shown. Colors indicate the position in the neutral direction. The snapshot clearly shows the formation of lanes tilted in the flow–gradient plane created due to sliding over rows of oriented molecules. Compared to Fig. 14, the snapshot is rotated for 180° with respect to the gradient axis.  
8 Effect of tangential friction
To investigate how the tangential friction between monomer blobs alters the rheology of the melt model, we use a standard (i.e. not adaptive) DPD thermostat^{70} with γ_{∥} = 1 and γ_{⊥} = 1 (see Section 3, 4, and 5.3 for details). As stated the friction kernels are Heaviside functions, in this case with a cutoff distance R_{dpd} = 1.5 × 2^{1/6}σ. In the following, we label this tangential friction thermostat as “tdpd”. This tdpd thermostat was found to be strong enough to keep the relative increase of the system temperature smaller than 5% at the largest shear rates considered.
In a real system friction acts by reducing the relative velocities of interacting monomers, generally in the normal and also in the tangential directions. Under Markovian and pairwise interaction assumptions, this form of friction leads to the DPD equations as shown by the Mori–Zwanzing dynamic coarse graining applied to the microscopic Liovillian dynamics.^{9} The same effect is properly captured by the tdpd thermostat, although here at a qualitatively level. Reducing the monomer relative velocities immediately leads to a reduction in kinetic pressure which has large consequences on the system rheology. In particular, the behaviour of the melt is essentially ruled by its elastic component, activated by the more effective intermolecular friction. Just to illustrate this point, we plot in Fig. 16 the contributions of the first normal stress difference (a direct measure of viscoelasticity) for the tdpd case. Compared with Fig. 10 (for the γ_{⊥} = 0 adpd case) the tdpd model has a much smaller kinetic pressure and N_{1} is essentially determined by the elastic stress (particularly as Wi increases). The same conclusion applies to P_{12} in Fig. 16.

 Fig. 16 Balances for the first normal stress difference and the shear stress for the tdpd model including tangential friction. Comparison with Fig. 10 for the zero tangential friction case. Error bars of the measured quantities are approximately 0.002–0.005ε/σ_{0}^{3}.  
An interesting difference related to the presence of tangential friction concerns shear dilatancy. Fig. 17(a) presents the relative density expansion δρ/ρ_{0} = 1 − ρ/ρ_{0} for different models. Let us now focus on the adpd and tdpd models which are kept isothermal (nonisothermal cases are discussed in Section 9). The density expansion of the adpd model (without tangential friction) scales like δρ/ρ_{0} ∼ Wi, while tangential friction (tdpd) leads to a much softer trend δρ/ρ_{0} ∼ Wi^{0.5} (although it expands relatively more at moderate shear rates). Under a constant normal load, shear dilatancy is a consequence of the growth of pressure in the velocity–gradient direction. In the case of small kinetic effects (tdpd) this growth is controlled by the expansion force arising from the compressed springs. This elastic pressure appears as soon as molecules start to align with the flow and compress in the gradient direction. Under enough tangential friction, the elastic stress is dominant and also controls the hydrostatic pressure, which in the absence of kinetic pressure effects, presents a faster decrease at large Wi compared to the adpd case (see Fig. 18). Of course, this decrease is also related to the fact that the tpdp simulations were done in the open system; notably for tdpd we get about 50% reduction in hydrostatic pressure for less than 10% reduction in density (see Fig. 17).

 Fig. 17 Top panel: Relative density variation for the different cases studied (open and closed, and different models and thermostats, see Table 1). Bottom panel: relative temperature increase observed in simulations with standard DPD thermostats with vanishing tangential friction. The relative variation of T for the adpd thermostat is smaller than 0.01 and less than 0.045 for the tdpd case (both of them not shown in the graph). Details of the thermostats are given in Table 1. The lines correspond to eqn (21), with the characteristic constant A defined there. Error bar estimations of the data are approximately 0.005 and 0.01 for the upper and lower graph, respectively.  

 Fig. 18 Hydrostatic pressure and components of the traceless stress for the adpd model (normal friction between blobs) compared with the tdpd model (including tangential friction). Both cases in open boxes at fixed T = 4. Error bars of the measured quantities are approximately 0.002–0.005ε/σ_{0}^{3}.  
The eigenvalues of the gyration tensor shown in Fig. 8 also indicate that adding tangential friction makes star molecules “stiffer”, in the sense that one needs larger values of the Weissemberg number to deform them. This observation is however somewhat misleading because for a fixed Wi, the real (physical) shear rate = Wi/τ_{rot} is now smaller due to the increase in τ_{rot} with the friction. In any case, the tangential friction is expected to alter the stress–strain relations in the Hookean regime (related to the linear stressoptic rule coefficient). This is (indirectly) seen in Fig. 11 where we plot the normal strain differences G_{11} − G_{22} and G_{22} − G_{33} (also G_{12}) against the corresponding intermolecular stress differences (against P^{intra}_{12}). We choose this plot to illustrate two facts: first, if the kinetic pressure is minor, the intermolecular friction is the leading mechanism driving the molecular deformation (and its elastic response). Second, molecular strains (and elastic stresses, not shown in the figure) in the flow–gradient and gradient–neutral planes (corresponding to the first and second normal stress differences) present a quite similar linear scaling (the shear deformation as well) with the intermolecular stress. This is to be contrasted with the top pannel of Fig. 11 (model in the absence of tangential friction), where the second normal stress (and G_{22} − G_{33}) is half of the first N_{1} counterpart (also G_{22} − G_{11}). A perfect alignment between intermolecular forces, elastic stress and molecular strain was also found by Kroger et al.^{89} in linear FENE chains. It thus seems that tangential friction (tdpd) helps to reduce the second normal stress in such a way that N_{1} and N_{2} present similar scaling laws N_{1} ∼ ^{0.68±0.02} with N_{2} ∼ −0.3N_{1}. This is to be compared with the adpd case in Fig. 6 (N_{1} ∼ ^{0.70} and N_{2} ∼ ^{1.0}). As shown in Fig. 7 the tdpd model yields −N_{2}/N_{1} ≃ 0.3 at any Wi < 100. This value is characteristic of disentangled melts (being −N_{2}/N_{1} = 2/7 the theoretical prediction for small shear rates^{1,88}).
The monotonous increase of elastic storage with found in the tdpd model is reflected in the recoverable shear strain (SR) shown in Fig. 7. Somewhat paradoxically, adding tangential friction increases the melt elasticity. Albeit, this reinforces the conclusions in Section 6: the intermolecular friction is the principal mechanism loading elastic stress into an disentangled melt. In passing, we note that in the tpdp model the orientational resistance parameter m_{G} = Witan(2θ_{G}) grows like m_{G} = 3.7Wi^{0.65} (at least for Wi < 100), a scaling which agrees with that reported for stars in solution^{105} (the prefactor being however about twice larger in our melt).
Not unexpectedly, the zero shear viscosity for the tdpd star model is larger η_{0} = 2.6 than the adpd case η_{0} = 0.5. The relaxation time is also larger τ_{η} = 287 (compared with 60). However, the tdpd shear viscosity thinning is faster; we find β_{η} ≃ 0.5 compared with 0.4 for the adpd case (recall η ∼ ^{−βη}). Again, this is also a consequence of a much smaller contribution of the viscous stress coming from kinetic effects. If the tangential friction is absent, the kinetic (and intramolecular) contribution increases the shear stress and the viscosity at any Wi leading to a softer shear thinning exponent.
9 Thermostats and heat dissipation
9.1 Density and temperature
We now briefly analyze the effect of the temperature increase due to heat dissipation in the sample. In all cases, the temperature reaches a steady state, but as shown in Fig. 17 (bottom panel) plotting δT/T_{0} = T/T_{0} − 1, we face severe viscous heating while performing the first row of simulations with the sdpdshort and (stronger) sdpdlong thermostats (see Table 1). Heating is observed at Wi > 10 and irrespective of the damping parameter (we tried up to γ_{DPD} = 50m_{0}τ_{0}^{−1}). Tangential friction drastically reduces heating (δT/T_{0} < 0.045 for Wi < 70) however the adpd thermostat enabled us to simulate zero tangential friction at a fixed temperature, providing δT/T_{0} < 0.01.
Let us focus on the “heated” runs at the increasing shear rate to illustrate the effect of uncontrolled temperature. Fig. 17 shows that heating introduces further melt expansion under shear and this kinetic energy induces larger hydrostatic pressure (which it is seen to increase with shear in the sdpdshort case). Pressurization due to viscous heating can also alter the rheology response. This is seen in Fig. 5 where the sdpdshort case presents shear thickening at Wi > 10, but only under closed conditions. The shear thickening reported in some of the published studies on polymer melts (closed box simulations) might be due to viscous heating (see e.g.ref. 39). More interesting than this elucubration is the results of Fig. 5 for the sdpdshortopen. The viscosity obtained for the same thermostat in an open environment is not affected by the temperature increase with . In fact it presents the very same trend as the adpd case. This observation indicates that shear viscosity is dominated by the normal load. In fact, the same outcome was also observed for the sdpdlong model (with shear exponent β_{η} = 0.39) and for all γ_{∥} ≤ 10 considered. Thus, this insensitivity of viscosity to temperature found only under normal load is not probably due to a cancellation of effects sometimes observed in experiments^{50} (viscosity decreasing with T and increasing with P_{22}). Rather it should be due to the viscosity dominated by pressure as happens in highly pressurized melts.^{50} Here, monomers interact via purely repulsive forces (WCA) which might map a sample under large pressure, in fact, adding attractive interactions would probably trigger temperature effects on viscosity.
9.2 Viscous heating
The rate of heat production per unit volume due to viscous dissipation is _{η} = η^{2} leading to a larger steady temperature whose value depends on the heat extraction rate. The onset of temperature increase is usually determined by a nondimensional parameter which depends on Wi and on the rate of cooling _{c} (see the recent computational study in ref. 10). Although in this work we shall not focus on heat and entropy production, we believe it is interesting to share our observations on this phenomenon, partly because of the relatively small simulation literature accurately reporting heating effects in sheared, thermostatted melts. A simple equation for the heat produced in the sheared melt includes frictional gain and cooling, 
 (19) 
where c_{X} is the specific heat capacity (molar) at a constant pressure (X = p) or volume (X = V) and dQ_{c} = ρ_{n}c_{X}dT (here ρ_{n} is the monomer number density). The DPD thermostat extracts (kinetic) energy upon pair collisions, at a rate which is proportional to the temperature difference T − T_{0}, where T is the temperature of the system (kinetic) and T_{0} the nominal temperature of the thermostat.§ The value of cooling rate α (which has units of the number of colliding pairs divided by volume and time) should scale as 
 (20) 
where we note that w(r) is the DPD kernel, which is simply a Heaviside function in our thermostat. The heat cooling rate of the thermostat increases with its damping coefficient γ_{DPD}, the kernel cutoff R^{DPD}_{cut}, and with the square of the monomer local density, which determines the number of thermalizing monomer collisions. In the steady state = 0 we thus have, 
 (21) 
Using the viscosity trend η = η(Wi) obtained from simulations (see Section 5.3 above), we plot prediction eqn (21) in Fig. 17 (lines) for different cases considered. The agreement is quite reasonable, indicating that the temperature increase in the melt can be forecast using a simple thermodynamic argument. Best fits to eqn (21) provide A ≡ αc_{P} = 0.011 in the opensdpdshort thermostat, while somewhat larger A = 0.017 for the closedsdpdshort thermostat. For the closed sdpdlong thermostat, which has about twice as much colliding partners within the dpd kernel, we consistently get A = 0.03. A preliminary calculation of c_{X} from the variance of the system energy in a closed (NVT) equilibrium simulation c_{V} = 〈δU^{2}〉_{NVT}/(NT^{2}) provides c_{V} ≃ 1.7. From eqn (21) we get the same order of magnitude (A ≃ 0.05 for sdpdshort and A = 0.09 for sdpdlong). A better agreement is found while comparing with the tdpd (sdpdlong) thermostat whose best fit provides A = 0.06(5) against the prediction 0.09. This indicates that eqn (20) should also depend on the number of degrees of freedom the thermostat acts upon (3 in the case of tdpd, 1 otherwise). A more refined calculation would also require including the dependence of c_{X} and g(r) on the shear rate.
10 Comparison with previous studies
It is interesting to compare our results with previous rheological studies, some of them were carried out at isobaric (constant P_{iso}) or constant load (P_{22}) constraints. As stated in the Introduction, the number of studies of flowing melts under a constant pressure (either hydrostatic pressure or normal load) is not large. However, they present significant discrepancies in the density and pressure variations with shear. For instance, ref. 42 presents the results for dendrimer melts under the isobaric conditions (constant hydrostatic pressure) revealing a decrease in the melt density under shear. For linear chains, ref. 40 presents just the opposite effect (contraction under shear), while ref. 43 (constant load) reports shear expansion (density increase).
The shear thinning exponents found here for a star polymer melt are consistent with those found in other simulations for somewhat similar systems such as hyperbranched and dendrimer polymers.^{41,42} Recall that the shear thinning exponent of any quantity Φ is β_{Φ} with Φ ∼ ^{−βΦ} at large . For viscosity we find β_{η} ≃ 0.4 (adpd) and β_{η} ≃ 0.5 (tdpd), for the first normal stress coefficient β_{Ψ1} ∼ 1.30 (adpd) and 1.31 (tdpd), while for the second one the values are β_{Ψ2} ∼ 1.0 (adpd) and 1.31 (tdpd). For dendrimers, Bosko et al.^{41} reported shear thinning exponents increasing with M under NVT (β_{η} ∈ [0.28–0.36]) while, for NP_{iso}T the were roughly independent of M (β_{η} ∈ [0.37–0.39]). The same work reported β_{Ψ1} ≃ 1.27 and β_{Ψ2} ≃ 1.23 under NP_{iso}T, while β_{Ψ1} ≃ 1.1 β_{Ψ2} ≃ 1.0 under NVT. Closed simulations for hyperbranched polymer melts^{42} predicted β_{η} ≃ 0.3 (slightly increasing with molecular mass M), while β_{Ψ2} ≃ 0.95 and β_{Ψ1} ≃ 1 (both roughly independent on M). The scaling of (first and second) normal coefficients are probably sensitive to the type of external constraint (either isobaric P_{iso} or constant load P_{22}). As an indication, a numerical study at constant load for linear chains^{43} reported the same exponents found in this work (β_{Ψ1} = 1.35 and β_{Ψ2} ≃ 1 respectively). Concerning shear viscosity, if the density is allowed to relax (under isobaric or constant load conditions) the shear thinning exponents reported in the literature show much less variation with the molecular weight than under a fixed density. This has been shown for linear chains under a constant load^{43} and also for dendrimers.^{42} In both cases, the shear thinning exponent in the variable density case was found to be close to 0.4, while it ranged from 0.25 to 0.4 in the closed system (increasing with the molecular weight in both cases). This larger insensitivity of shear thinning exponents under fixed load conditions agrees with our observations in simulations presenting viscous heating mentioned in Section 9.
More recently, very similar papers^{106,107} studied the rheology and dynamics^{108} of star polymers with different functionalities in solution. They reported that the contribution of the stars to the sheared solution viscosity scales like Wi^{−0.4} which is quite close to the shear thinning exponents we find here for the melt. According to these simulations,^{106,107} the first normal stress coefficient scales like Ψ_{1} ∼ Wi^{−1} in solution, although the authors claim an exponent of −4/3 at very high shear rates Wi > 100. In the melt, we observe the −4/3 power law at smaller shear rates. The second normal coefficient seems also to scale slightly differently in solution Ψ_{2} ∼ Wi^{−4/3} than in the melt Ψ_{2} ∼ Wi^{−1}. The strong similarities in solution and melt indicate that rheological properties are mainly ruled by conformational changes of the chains and maybe also that the hydrodynamic effects in the melt are somewhat similar to those in a liquid solvent. Finally, it is interesting to note that the range of values for β_{η} and β_{Ψ1} for stars, dendrimers and hyperbranched molecules are consistent with the theoretical calculation of Kroger et al.^{88} based on the Fokker–Planck equation for the bond vector distribution of multibead linear chains having slow rotational diffusion (compared with the entanglement relaxation) and a finite deformation energy.^{88} These two features are indeed consistent with the nature of shortarmed stars, dendrimers and hyperbranched molecules with internal excluded volume interactions.
The number of experimental studies on rheology of sheared star polymers is not abundant either, but it is growing fast due to the foreseen technological applications.^{52,53} Most of these studies consider long armed stars which present entanglements. However, the shear thinning exponent we found for the tangential friction case β_{η} = 0.5 is quite close to the star case of the experiments by Tezel et al.^{109} (4 arms stars with M_{a} = 140 kg mol^{−1}; the lowest molecular weight considered in these experiments). Although entanglements in star polymers are not yet fully understood in the nonNewtonian regime,^{52} they should be responsible for the substantially increase in shear thinning (exponents close to 0.8) observed in the recent experiments with star molecule melts of Snijkers et al.^{52} Shear thinning exponents found in melts of linear chains are around β_{η} = 1 (predicted by reptation theory). An open question is why star molecules significantly reduce the shear thinning.^{52}
The analyses presented in Section 6 based on the (exact) balance of the pressure coming from the different contributions (kinetic, intramolecular, intermolecular and bonds) provide some clues which can be useful in a broader study. These types of analysis were performed by Baig et al.,^{2,3} Matin^{43} and Tschopp.^{45} However these balances were based on the energy budget which has less direct rheological consequences than the stress. Also, surprisingly, the kinetic contribution was not considered (leading to unbalanced analyses). The sole exception we found is the work of Kroger et al.^{89} for linear FENE chains which, like ours, is based on the exact stress budget. In elongational flow of linear chain melts, Kroger et al.^{89} found that viscous effects came mainly from intramolecular collisions (they called it “simple fluid” stress offset). Here, we find that in star molecules the kinetic pressure can be the dominant viscous contribution to the melt rheology at large shear rates (a possibility in fact recognized by Kroger et al. in their book^{4}). The kinetic pressure is however reduced with the tangential friction and this warns about the need of dynamic coarsegraining^{9} to represent a realistic model. While linear chains tumble in shear flow,^{110} star molecules perform a quite different motion called tankthreading^{102,105} (whereby arms rotate around the molecule center). One can speculate that the reduction of shear thinning observed in star melts is due to an enhanced kinetic pressure related to the tankthreading. Above Wi > 10 we observe the tankthreading motion in our melt with monomer angular momentum growing much faster (ω ∼ ^{0.6} from the preliminary results) than it does in stars in solution (where ω is constant^{102,111}). Concerning intermolecular interactions, quite often neglected in the literature (see e.g.ref. 87), we find that these are the key to transfer the externally imposed stress into molecular strain and elastic stress. This observation agrees with that of Kroger et al.^{89} Intermolecular friction is the sole possible mechanism if entanglements are not relevant (like in our short arm star molecules).
11 Conclusions
We have conducted OBMD simulations of the melt of star polymers (73 monomers, 12 arms of 6 monomers per arm) using a relatively new modelling technique which combines the adaptive resolution and openMD (respectively, introduced in ref. 22 and 46 and used in many other studies). From a technical point of view, this work substantially enlarges the size of molecules exchanged between the open system and the reservoir.^{13} In this work, we however focus on what happens if a melt is sheared under a constant normal load instead of a constant volume. The constant load is in fact the condition in many experiments (see ref. 40 and 44) although the simulation literature on the subject is not abundant. We have also presented some conclusions on how the tangential and normal monomer–monomer frictions affect these CG molecular models.
The OBMD permits to perform several new features in a molecular simulation. At equilibrium, the OBMD correctly represents the mass fluctuations of the grand canonical ensemble and interestingly, it could permit to study how fluctuations and sound velocity are altered under nonequilibrium (e.g. sheared) states. OBMD also allows to impose the external shear stress P^{ext}_{12} as required for (flux based) hybrid continuumMD simulations.^{15,66} It could also enable the validation of theories like extended thermodynamics^{47} predicting different outcomes for the “conjugate” constant stress and constant shear rate nonequilibrium ensembles.
Concerning the present study, we observe that under a constant normal pressure, the melt expands when sheared (shear dilatancy) leading to substantial depressurization and slightly decreasing the shear viscosity. This behaviour was observed in most previous simulations on sheared melts, but surprisingly it is still unclear if it is the general (universal) trend (e.g. see ref. 40). This study provides new information on the rheology of sheared melts: notably, we see that the type of monomer friction is a key aspect for the system rheology. From a theoretical standpoint,^{9} friction between monomers (or rather “polymer blobs”) arises as soon as one consider a CG view of the complex molecule (which is the standard case in polymer science). The pressure balance analysis reveals that in the absence of tangential friction, the monomer kinetic stress becomes significant at large shear rates, increasing the system viscosity (reducing the shear thinning exponent) and diminishing its elastic response (e.g. normal stress difference and stress recovery). We also observed viscous heating in some simulations (e.g. using the sdpdshort and sdpdlong thermostats) revealing a viscosity jump (shear thickening) in closed systems. By contrast, the viscosity of the heated (and less dense) open samples did not change in trend (shear thinning) under a constant load. This indicates that the melt viscosity is controlled by the normal pressure, at least for the present type of molecules with purely repulsive interactions and no significant entanglements. We expect the OBMD to become useful in other studies, such as adding the energy transfer^{63} or extending the incompressible coupling in ref. 10 to perform hybrid simulations of compressible melts including the transfer of dissipated heat through/across system boundaries. Thus, the properties related to heat conduction could be investigated.^{112}
Acknowledgements
We thank J. J. Freire for initial configurations and K. Kremer for useful discussions. J. S. and M. P. acknowledge financial support through grants P10002 and J14134 from the Slovenian Research Agency. J. S. acknowledges financial support from Slovene Human Resources Development and Scholarship Fund (186. JR). R. DB acknowledges support from the Spanish government under national MINECO project FIS201347350C51R. Partial support from the COST Action MP1305 is kindly acknowledged.
References

R. Larson and H. Brenner, Constitutive Equations for Polymer Melts and Solutions: Butterworths Series in Chemical Engineering, Elsevier Science, 2013 Search PubMed.
 C. Baig and V. G. Mavrantzas, J. Chem. Phys., 2010, 132, 014904 CrossRef PubMed.
 C. Baig, V. G. Mavrantzas and M. Kröger, Macromolecules, 2010, 43, 6886–6902 CrossRef CAS.

M. Kröger, Models for Polymeric and Anisotropic Liquids, Springer, Berlin Heidelberg, 2005 Search PubMed.
 J. Padding, E. van Ruymbeke, D. Vlassopoulos and W. Briels, Rheol. Acta, 2010, 49, 473–484 CrossRef CAS.
 V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719 CrossRef CAS.
 C. F. Abrams, J. Chem. Phys., 2005, 123, 234101 CrossRef PubMed.
 P. Ilg, H. C. Öttinger and M. Kröger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 011802 CrossRef PubMed.
 C. Hijon, P. Espanol, E. VandenEijnden and R. DelgadoBuscalioni, Faraday Discuss., 2010, 144, 301–322 RSC.
 S. Yasuda and R. Yamamoto, Phys. Rev. X, 2014, 4, 041011 Search PubMed.
 R. DelgadoBuscalioni, K. Kremer and M. Praprotnik, J. Chem. Phys., 2008, 128, 114110 CrossRef PubMed.
 R. DelgadoBuscalioni, K. Kremer and M. Praprotnik, J. Chem. Phys., 2009, 131, 244107 CrossRef PubMed.
 R. DelgadoBuscalioni, J. Sablić and M. Praprotnik, Eur. Phys. J.: Spec. Top., 2015, 224, 2331–2349 CrossRef CAS.
 R. DelgadoBuscalioni, Numerical Analysis of Multiscale Computations, 2012, 145–166 Search PubMed.
 G. De Fabritiis, R. DelgadoBuscalioni and P. V. Coveney, Phys. Rev. Lett., 2006, 97, 134501 CrossRef CAS PubMed.
 K. Mohamed and A. Mohamad, Microfluid. Nanofluid., 2010, 8, 283–302 CrossRef.
 J. H. Walther, M. Praprotnik, E. M. Kotsalis and P. Koumoutsakos, J. Comput. Phys., 2012, 231, 2677–2681 CrossRef CAS.
 S. Yasuda and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 036308 CrossRef PubMed.
 W. E, B. Engquist, X. Li, W. Ren and E. VandenEijnden, Commun. Comput. Phys., 2007, 2, 367–450 Search PubMed.
 N. D. Petsev, L. G. Leal and M. S. Shell, J. Chem. Phys., 2015, 142, 044101 CrossRef PubMed.
 D. A. Fedosov and G. E. K. B. Caswell, J. Chem. Phys., 2010, 132, 144103 Search PubMed.
 M. Praprotnik, L. Delle Site and K. Kremer, J. Chem. Phys., 2005, 123, 224106 CrossRef PubMed.
 M. Praprotnik, L. Delle Site and K. Kremer, Annu. Rev. Phys. Chem., 2008, 59, 545–571 CrossRef CAS PubMed.
 M. Praprotnik, S. Poblete and K. Kremer, J. Stat. Phys., 2011, 145, 946–966 CrossRef.
 S. Bevc, C. Junghans, K. Kremer and M. Praprotnik, New J. Phys., 2013, 15, 105007 CrossRef.
 H. Wang, C. Hartmann, C. Schütte and L. Delle Site, Phys. Rev. X, 2013, 3, 011018 Search PubMed.
 J. Zavadlav, M. N. Melo, A. V. Cunha, A. H. de Vries, S. J. Marrink and M. Praprotnik, J. Chem. Theory Comput., 2014, 10, 2591–2598 CrossRef CAS PubMed.
 J. Zavadlav, M. N. Melo, S. J. Marrink and M. Praprotnik, J. Chem. Phys., 2014, 140, 054114 CrossRef PubMed.
 A. Agarwal, H. Wang, C. Schütte and L. Delle Site, J. Chem. Phys., 2014, 141, 034102 CrossRef PubMed.
 J. Zavadlav, M. N. Melo, S. J. Marrink and M. Praprotnik, J. Chem. Phys., 2015, 142, 244118 CrossRef PubMed.
 J. Zavadlav, R. Podgornik and M. Praprotnik, J. Chem. Theory Comput., 2015, 11, 5035–5044 CrossRef CAS PubMed.
 H. Wang and A. Agarwal, Eur. Phys. J.: Spec. Top., 2015, 224, 2269–2287 CrossRef CAS.
 A. Agarwal, J. Zhu, C. Hartmann, H. Wang and L. D. Site, New J. Phys., 2015, 17, 083042 CrossRef.
 K. Kreis, A. C. Fogarty, K. Kremer and R. Potestio, Eur. Phys. J.: Spec. Top., 2015, 224, 2289–2304 CrossRef CAS.
 D. Mukherji and K. Kremer, Macromolecules, 2013, 46, 9158–9163 CrossRef CAS.
 P. Szymczak and M. Cieplak, J. Chem. Phys., 2007, 127, 155106 CrossRef CAS PubMed.
 P. Szymczak and M. Cieplak, J. Phys.: Condens. Matter, 2011, 23, 033102 CrossRef PubMed.
 B. Z. Dlugogorski, M. Grmela and P. J. Carreau, J. NonNewtonian Fluid Mech., 1993, 48, 303–335 CrossRef CAS.
 P. J. Daivis and D. J. Evans, J. Chem. Phys., 1994, 100, 541–547 CrossRef CAS.
 Z. Xu, J. J. de Pablo and S. Kim, J. Chem. Phys., 1995, 102, 5836–5844 CrossRef CAS.
 J. T. Bosko, B. D. Todd and R. J. Sadus, J. Chem. Phys., 2005, 123, 034905 CrossRef PubMed.
 T. C. Le, B. D. Todd, P. J. Daivis and A. Uhlherr, J. Chem. Phys., 2009, 131, 044902 CrossRef CAS PubMed.

M. L. Matin, PhD thesis, Royal Melbourne Institute of Technology, Melbourne, 2001.
 P. A. Thompson, G. S. Grest and M. O. Robbins, Phys. Rev. Lett., 1992, 68, 3448–3451 CrossRef PubMed.

M. Tschopp, J. L. Bouvard, D. Ward, D. Bammann and M. F. Horstemeyer, 2013, arXiv:1310.0728, 1–28.
 E. G. Flekkoy, R. DelgadoBuscalioni and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 026703 CrossRef PubMed.

D. Jou, J. CasasVázquez and M. CriadoSancho, Thermodynamics of Fluids Under Flow, Springer, Netherlands, 2010 Search PubMed.
 P. Daivis, J. NonNewtonian Fluid Mech., 2008, 152, 120–128 CrossRef CAS.

J. Dealy and J. Wang, Melt Rheology and its Applications in the Plastics Industry, Springer, Netherlands, 2013 Search PubMed.

J. Aho, Rheological Characterization of Polymer Melts in Shear and Extension: Measurement Reliability and Data for Practical Processing, Tampere University of Technology, 2011 Search PubMed.
 W. Friesenbichler, I. Duretek, J. Rajganesh and S. R. Kumar, Polimery, 2011, 56(1), 58–62 Search PubMed.
 F. Snijkers, K. Ratkanthwar, D. Vlassopoulos and N. Hadjichristidis, Macromolecules, 2013, 46, 5702–5713 CrossRef CAS.
 D. Vlassopoulos, G. Fytas, T. Pakula and J. Roovers, J. Phys.: Condens. Matter, 2001, 13, R855–R876 CrossRef CAS.
 C. N. Likos, Phys. Rep., 2001, 348, 267–439 CrossRef CAS.
 C. N. Likos, H. Löwen, M. Watzlawek, B. Abbas, O. Jucknischke, J. Allgaier and D. Richter, Phys. Rev. Lett., 1998, 80, 4450–4453 CrossRef CAS.

S. M. T. Divoux, M. A. Fardin and S. Lerouge, 2015, arXiv:1503.04130, 1–16.
 S.Q. Wang, S. Ravindranath and P. E. Boukany, Macromolecules, 2011, 44, 183–190 CrossRef CAS.
 S. T. Milner, T. C. B. McLeish and A. E. Likhtman, J. Rheol., 2001, 45, 539–563 CrossRef CAS.
 C. Koch, A. Z. Panagiotopoulos, F. Lo Verso and C. N. Likos, Soft Matter, 2013, 9, 7424–7436 RSC.
 D. Reith, M. Pütz and F. MüllerPlathe, J. Comput. Chem., 2003, 24, 1624–1636 CrossRef CAS PubMed.
 S. Bevc, C. Junghans and M. Praprotnik, J. Comput. Chem., 2015, 36, 467–477 CrossRef CAS PubMed.
 R. DelgadoBuscalioni and P. V. Coveney, J. Chem. Phys., 2003, 119, 978–987 CrossRef CAS.
 R. Potestio, S. Fritsch, P. Español, R. DelgadoBuscalioni, K. Kremer, R. Everaers and D. Donadio, Phys. Rev. Lett., 2013, 110, 108301 CrossRef PubMed.
 R. Potestio, P. Español, R. DelgadoBuscalioni, R. Everaers, K. Kremer and D. Donadio, Phys. Rev. Lett., 2013, 111, 060601 CrossRef PubMed.
 P. Español, R. DelgadoBuscalioni, R. Everaers, R. Potestio, D. Donadio and K. Kremer, J. Chem. Phys., 2015, 142, 064115 CrossRef PubMed.
 R. DelgadoBuscalioni and G. De Fabritiis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 036709 CrossRef CAS PubMed.
 E. A. Koopman and C. P. Lowe, J. Chem. Phys., 2006, 124, 204103 CrossRef CAS PubMed.
 T. Soddemann, B. Dünweg and K. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 046702 CrossRef PubMed.
 P. Espanol and P. Warren, Europhys. Lett., 1995, 30, 191 CrossRef CAS.
 C. Junghans, M. Praprotnik and K. Kremer, Soft Matter, 2008, 4, 156–161 RSC.
 A. W. Lees and S. F. Edwards, J. Phys. C: Solid State Phys., 1972, 5, 1921 CrossRef.
 D. J. Evans and G. P. Morriss, Phys. Rev. A: At., Mol., Opt. Phys., 1984, 30, 1528–1530 CrossRef CAS.
 A. J. Ladd, Mol. Phys., 1984, 53, 459–463 CrossRef CAS.

J. Freire, in Branched Polymers II, ed. J. Roovers, Springer, Berlin Heidelberg, 1999, vol. 143, pp. 35–112 Search PubMed.
 M. SJ, KoreaAust. Rheol. J., 2008, 20, 117–125 Search PubMed.
 F. Legrand and J.M. Piau, J. NonNewtonian Fluid Mech., 1998, 77, 123–150 CrossRef CAS.
 R. A. Mendelson, Trans. Soc. Rheol., 1965, 9, 53–63 CrossRef CAS.
 J. CasasVázquez and D. Jou, Rep. Prog. Phys., 2003, 66, 1937 CrossRef.
 S. Bernardi, B. D. Todd and D. J. Searles, J. Chem. Phys., 2010, 132, 244706 CrossRef PubMed.
 P. Padilla and S. Toxvaerd, J. Chem. Phys., 1996, 104, 5956–5963 CrossRef CAS.
 R. Khare, J. de Pablo and A. Yethiraj, J. Chem. Phys., 1997, 107, 2589–2596 CrossRef CAS.
 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690 CrossRef CAS.
 M. Christen and W. F. van Gunsteren, J. Chem. Phys., 2006, 124, 154106 CrossRef PubMed.
 G. S. Grest, K. Kremer, S. T. Milner and T. A. Witten, Macromolecules, 1989, 22, 1904–1910 CrossRef CAS.
 G. S. Grest, K. Kremer and T. A. Witten, Macromolecules, 1987, 20, 1376–1383 CrossRef CAS.
 A. Eriksson, M. N. Jacobi, J. Nyström and K. Tunstrom, J. Chem. Phys., 2008, 129, 024106 CrossRef PubMed.

G. Strobl, The Physics of Polymers, SpringerVerlag, NY, 1997 Search PubMed.
 M. Kröger, W. Loose and S. Hess, J. Rheol., 1993, 37, 1057–1079 CrossRef.
 M. Kröger, C. Luap and R. Muller, Macromolecules, 1997, 30, 526–539 CrossRef.

J.P. Hansen and I. R. McDonald, in Theory of Simple Liquids, ed. J.P. Hansen and I. R. McDonald, Academic Press, Burlington, 3rd edn, 2006, pp. 11–45 Search PubMed.
 R. DelgadoBuscalioni and A. Dejoan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 046708 CrossRef CAS PubMed.
 R. DelgadoBuscalioni and P. V. Coveney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 046704 CrossRef CAS PubMed.

P. Allen and D. Tildesley, Computer simulation of liquids, Clarendon Press, 1987 Search PubMed.
 K. Yasuda, J. Text. Sci. Eng., 2006, 52, 171–173 CrossRef.
 J.P. Ryckaert, A. Bellemans, G. Ciccotti and G. V. Paolini, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 259–267 CrossRef.
 P. Kindt and W. J. Briels, J. Chem. Phys., 2005, 123, 224903 CrossRef CAS PubMed.
 C. Baig and V. A. Harmandaris, Macromolecules, 2010, 43, 3156–3160 CrossRef CAS.
 S. Barsky, R. DelgadoBuscalioni and P. V. Coveney, J. Chem. Phys., 2004, 121, 2403–2411 CrossRef CAS PubMed.

R. DelgadoBuscalioni, AIP Conf. Proceedings, 2007, 913, 114.
 B. D. Todd, D. J. Evans and P. J. Daivis, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 1627–1638 CrossRef CAS.
 R. L. C. Akkermans and G. Ciccotti, J. Phys. Chem. B, 2004, 108, 6866 CrossRef CAS.
 M. Ripoll, R. G. Winkler and G. Gompper, Phys. Rev. Lett., 2006, 96, 188302 CrossRef CAS PubMed.
 J. Brandao, E. Spieth and C. Lekakou, Polym. Eng. Sci., 1996, 36, 49–55 CAS.
 X. Cheng, X. Xu, S. A. Rice, A. R. Dinner and I. Cohen, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 63–67 CrossRef CAS PubMed.
 M. Ripoll, R. Winkler and G. Gompper, Eur. Phys. J. E: Soft Matter Biol. Phys., 2007, 23, 349–354 CrossRef CAS PubMed.
 D. A. Fedosov, S. P. Singh, A. Chatterji, R. G. Winkler and G. Gompper, Soft Matter, 2012, 8, 4109–4120 RSC.
 S. P. Singh, A. Chatterji, G. Gompper and R. G. Winkler, Macromolecules, 2013, 46, 8026–8036 CrossRef CAS.
 S. P. Singh, C.C. Huang, E. Westphal, G. Gompper and R. G. Winkler, J. Chem. Phys., 2014, 141, 084901 CrossRef PubMed.
 A. K. Tezel, J. P. Oberhauser, R. S. Graham, K. Jagannathan, T. C. B. McLeish and L. G. Leal, J. Rheol., 2009, 53, 1193–1214 CrossRef CAS.
 R. DelgadoBuscalioni, Phys. Rev. Lett., 2006, 96, 088303 CrossRef PubMed.
 S. P. Singh, D. A. Fedosov, A. Chatterji, R. G. Winkler and G. Gompper, J. Phys.: Condens. Matter, 2012, 24, 464103 CrossRef PubMed.
 J. Li and J. Lee, Acta Mech., 2014, 225, 1223–1233 CrossRef.
Footnotes 
† An interesting remark is that the recoverable shear compliance^{87} obtained from J_{e} = τ_{rot}/η_{0} results to be J_{e} ≃ 110 independent of the friction kernel. 
‡ Traditionally, following the experiment standpoint, rheological magnitudes are expressed in terms of the pressure on the system (here given by −P). For instance, the normal stress difference is then N_{1} = −(P_{11} − P_{22}) and this explains the opposite sign used in our (and other) simulation studies. 
§ This can be proved from the equation for the time dependence of the covariance 〈dv_{i}/dt(t)dv_{i}/dt(t′)〉 of the DPD Langevin's equation and can be easily checked (and the cooling rate α measured) upon observation of an exponential convergence of T towards T_{0} after an instantaneous change (increase or decrease) of T_{0}. 

This journal is © The Royal Society of Chemistry 2016 