Equations of motion for polymer chains in a thermostat

A constrained dynamics method suitable for molecular dynamics simulations is considered to study the long-time dynamics of polymer chains. The method is initially discussed on the basis of the Lagrangian and Hamiltonian formalisms for isolated polymer chains with fixed monomer–monomer links. Subsequently, the corresponding equations of motion are obtained for describing the dynamics of such polymer chains in the presence of a thermostat. The approach is applied to a few typical cases to illustrate how the formalism is implemented numerically and to elucidate its convergence properties when studying such systems in equilibrium. As an example, we consider the problem of reconstructing the backbone structure (chain of Cα atoms) of protein Rubredoxin from its contact matrix. It is shown that the target structure is succesfully reached after a long transient regime (typically in the range from 106 to 108 integration steps). A particular attractive extension of the algorithm presented here is to environment-dependent couplings, which could allow the study of the long-time polymer dynamics in realistic environments. The present method is thus expected to have useful applications in the modelling of the complex dynamics of bio-polymers such as proteins, and also in the context of nanoscale polymer materials.


General remarks
Numerical simulations of the dynamics of polymer chains have been addressed by studying transport equations either in the form of a Fokker-Planck equation (often referred to as Smoluchowski equation also) or of Langevin type [1], by means of Monte Carlo (MC) methods [2] or molecular dynamics (MD) techniques [1]- [8]. In all of the above cases, one of the major challenges is a suitable description of the links between subsequent monomers in a polymer chain. One approach is to introduce a relatively stiff potential between monomers [1,2]. In actual simulations, it leads to rather small integration time steps and, as a result, one is forced to spend most of the time in simulating the polymer links rather than the relatively slow polymer chain dynamics. One obvious solution is to introduce holonomic constraints [4]- [6]; another is to define generalized coordinates [1,7,8] and thus deal only with the real degrees of freedom.
A further requirement in the case of polymers in solutions is to have a suitable description of the solvent and its interaction with the polymer. That can be done either by an explicit treatment of the solvent molecules [2,4] (which requires at least an order of magnitude or more molecules than the number of monomers) or by the introduction of a phenomenological friction and fluctuating forces as in the Langevin equation. The introduction of constraints and a phenomenological description of the solvent molecules allows one to concentrate one's attention on the most interesting feature, the polymer dynamics.
In this paper, we discuss the main theoretical aspects involved in the study of polymer chains (with either fixed bond lengths, or fixed bond lengths and fixed angles between successive links) coupled to a thermostat in the framework of a MD approach, in which hydrodynamic interactions are neglected. A very similar Lagrange and Hamiltonian dynamics for systems with holonomic constraints was developed independently in [9]. The Hamiltonian formalism developed here, however, is somewhat more transparent and the manner in which the holonomic constraints are introduced is more direct and in clear relationship with Dirac's formalism [10]. The main difference between our approach and that of [9], for example, is in the way we couple the polymer chain to the thermostat. The coupling to the thermostat described here is more flexible in several aspects, it allows a controlled exchange of the energy flow between the system and the thermostat both in the time and in the spatial domains. In particular, it is very easy in our approach to allow the thermostat to couple to the system only at its surface, for example, which in many instances is clearly the manner by which the energy is exchanged between the system and the thermostat.
The paper is organized as follows: in section 2 we shall briefly review the Lagrangian dynamics of a polymer chain. In section 3 we discuss the Hamiltonian dynamics of a polymer chain, using the so-called Dirac-Poisson brackets [10], a method to treat holonomic and nonholonomic constraints, developed for the needs of quantum field theory. In section 4 we show how one can couple a polymer chain to a thermostat and in this way have a phenomenological description of the solvent and of the polymer-solvent interaction, which obviates an explicit treatment of the solvent molecules with very little effort. We shall briefly exemplify our methods with numerical simulations of the polymer dynamics in section 5 and present our main conclusions in section 6. A number of rather technical aspects of our approach are relegated to appendices A-E.

Lagrangian formalism
Let us assume that a polymer chain with N + 1 monomers with coordinates r k , where k = 0, . . . , N, is described by the Lagrangian L where U is the potential energy. (The formalism is identical for the case of unequal monomer masses.) If the links of the chain are rigid, one can derive the following Lagrange equations of motion for the chain where we have introduced the link vectors ρ k and the constraints X k (for k = 1, . . . , N) In the above equations, λ k s are Lagrange multipliers, corresponding to the fixed bond lengths. By using the second time derivative of the constraintṡ after some straightforward manipulations of the equations of motion, one can derive the following set of equations for the Lagrange multipliers λ k : As a last remark, once a chain has been specified, the initial velocities of the monomers cannot be specified in an arbitrary way, but have to be consistent with the constraints, namelẏ for k = 1, . . . , N. In this way, the length of each particular link will be conserved and also the relative velocity of two consecutive monomers will always be orthogonal to the instantaneous link vector ρ k = r k − r k−1 .

Hamiltonian formalism
In order to derive the Hamiltonian equations of motion for a polymer chain with fixed links, we shall use Dirac's formalism [10]. The reader who is not interested in the formalism and the derivation of the equations of motion can simply skip the whole section and refer only to the last equations in this section, equations (33) and (34). According to Dirac's formalism, besides the so-called primary constraints X k = 0, (see (6)), one has to introduce also the so-called secondary constraints Y k = 0 (for k = 1, . . . , N) where p k are canonical conjugate momenta to coordinates r k with respect to the usual Poisson bracket between two arbitrary functions A and B on the phase space r k , p k Following [10], let us gather together both primary and secondary constraints and denote them as Z µ , with µ = 1, . . . , 2N, where the first N quantities are the primary constraints X k and the following N quantities are the secondary constraints Y k . The Hamiltonian dynamics in the presence of constraints is now completely defined by the so-called Dirac-Poisson brackets It is convenient to introduce the following antisymmetric N × N matrices: As we will show below, see equations (24)-(27), the matrices U and T have a rather simple structure. Then the 2N × 2N matrices C and D can be easily shown to have the following structure: whereŨ stands for the transpose of U. The Dirac-Poisson bracket has the property that for an arbitrary Hamiltonian Since the dynamical evolution of any quantity A is now defined as [10] dA where all the primary and secondary constraints are automatically satisfied along the trajectory of the system, if they were satisfied initially. In other words, both primary and secondary constraints are integrals of motion of the Dirac-Poisson dynamics. Above we have introduced the quantities e µ as The set of quantities e µ is thus entirely defined in terms of the phase space variables r k , p k and the Hamiltonian of the system. Alternatively, the quantities e µ can be determined as the solution of the following set of linear equations: which follows from equations (13) and (22). It is straightforward to show that only the following Poisson brackets among the constraints are nonvanishing: and as a result the previous linear system of equations for e µ has essentially a tridiagonal form. After defining a k = e k and b k = e N+k , equations (23) can be rewritten as follows (with the convention that a 0 = a N+1 where In equation (29) one can introduce the quantities b k under the Poisson bracket because Y k = 0 along the physical trajectory. The equations of motion for the phase space variables r k , p k can be written as follows: A further simplification of the equations of motion occurs if one remembers that along the physical trajectory Y k = 0. In such a case, b k = 0, see equation (33), and the Hamiltonian equations for the phase space variables r k , p k becomė These are the Hamiltonian equations of motion we shall be using later. The quantities ρ k have been introduced in equation (4), Y l in equation (9) and V l in equations (14) and (15), respectively. It is not obvious that the system of linear equations (28) is always degenerate and the only solution is b k = 0. In particular, for a perfectly straight polymer chain this system of linear equations has a nontrivial solution, even though the right-hand side is vanishing identically. We shall neglect here such refinements of the theory, which can however be incorporated into the formalism [10]. It is somewhat surprising to find out that along the physical trajectory one has the usual relationship between velocities and momenta, see equations (33), as in the absence of constraints. For example, it is easy to convince oneself that r k is not canonically conjugate to p k with respect to the Dirac-Poisson bracket (12), see appendix A. Because of this somewhat unexpected coincidence, the right hand side of equation (34), has exactly the same form as in the Lagrangian formalism.

Coupling the polymer chain to a thermostat
Once the Hamiltonian dynamics in the presence of constraints is constructed, the coupling of the polymer dynamics to a thermostat can be performed as described in [11], which represents a generalization of the improved Nosé-Hoover dynamics [12] developed in [13]. In [11] two different types of couplings were introduced, so-called projected and twisted. We shall use here the twisted coupling of the polymer to the thermostat, since it is much easier to implement and it is also extremely efficient. One way of introducing the coupling to a thermostat can be described by adding to the Hamiltonian equations of motion couplings through the so-called pseudo-friction coefficients according to the following rule: where ζ is a pseudo-friction coefficient, τ is a dimensionless constant of order one, characterizing the strength of the coupling to the thermostat, and t 0 is a time constant of the order of the shortest period of the polymer chain. The reason for having the factor √ 2N in the above relations will be explained below. The form of this coupling can vary, for example, one can consider a different pseudo-friction coefficient for each cartesian direction and/or consider different other functions of p k and ζ, e.g., τζ 3 p 3 k . Even though we shall use more complicated couplings in the simulations, see appendix C, in order not to clutter the paper with complicated relations, that can be derived in a straightforward manner and have been discussed previously [11,13,14], we shall exemplify the method in this section only with the above coupling. Rather often ergodicity is not achieved if only one pseudo-friction coefficient is introduced and for that reason a coupling with several pseudo-friction coefficients is preferable [11,13,14] and such a scheme will be used in actual calculations.
In the presence of a thermostat the equations of motion becomė Note that the pseudo-friction coefficient ζ becomes a dynamical variable in this approach. T is the temperature of the thermostat (here in energy units, Boltzmann constant k B = 1). In the equation for p k , we have taken into account that ρ l · π l = 0 for all l = 1, . . . , N. (2N + 3) represents the number of degrees of freedom of a polymer chain with N + 1 monomers and N fixed bond lengths and thus NT is the canonical average of the total internal kinetic energy. In this way the average of the term N k=0 p 2 k /mT is equal to 2N if the internal motion is ergodic. If the dynamics is ergodic, the quantity in parentheses in equation (38) is a sum of independent random variables, with zero mean and unit variance. Thus the mean amplitude of this quantity is √ 2N and in this way the mean value ofζ is ≈τ/t 0 . Note that this type of coupling to the thermostat conserves the total linear momentum and also the direction of the total angular momentum of the chain (see appendix C).

The dynamics of a monomeric chain
We illustrate the method in the case of a chain of N = 50 identical monomers of unit mass M = 1, connected by rigid bonds of unit length. The interaction potential between monomers is taken as where R 0 = 1 and ε = 0.3 yields the strength of the attractive interaction. The minimum of V occurs at r min = √ 2/ε ≈ 2.582, yielding V min = −ε 2 /4 = −0.0225. The temperature is taken as T = 0.025 and the friction ζ obeys the equation of motion (see appendix C): with similar expressions for the y and z components. Here we use t 0 = 1 and τ = 2.42. Typical chain configurations attained in the early stages of the simulation are shown in figure 1. After about 10 7 integration steps, corresponding to a time t = 150, the chain attains a Each plot has been obtained after every 2 × 10 6 iterations, using the Adams-Bashforth-Moulton scheme [15] with an integration step τ = 1.5 × 10 −5 , corresponding to time intervals t = 30. The different chain configurations have been shifted from each other for clarity. rather compact structure. The resulting radius of gyration at equilibrium is R G ≈ 2.65 in units of bond length. Applications of this method to study protein structure will be considered in the next section.
The temporal behaviour of friction variables ζ x , ζ y and ζ z are shown in figure 2 for the simulation trajectory considered here. The probability distribution functions (PDFs) for friction and momenta have been calculated by recording all the corresponding values for times t > 150. While the standard deviation for each friction component is unity (cf equation (38) and text), the standard deviation of a monomer momentum component is given by, a value which is used in figure 3 in order to display the numerically obtained PDF. As one can see from figure 3, the PDF of the momentum components have reached their asymptotic form, while the PDF for the frictions still display departures from Gaussian behaviour for times t = 510.

Reconstruction of the backbone structure of a protein
In the following, we apply the present method to the reconstruction of the backbone structure of a protein, represented by the linear chain connecting its C α atoms. The structure of the backbone uniquely determines the contact matrix of the protein, from which the positions of all atoms in  the native structure can be obtained (see e.g., [16,17]). We consider as target structure a (wildtype) Rubredoxin protein, choosing in particular the 1e8j structure [18]. From the all-atoms' coordinates [19], we extract the positions of the N = 51 C α atoms, which define our target backbone chain (see figure 4). For convenience we take R u as the unit length, the C α -C α mean distance in 1e8j which turns out to be R u = 3.8 Å.
We define two C α atoms, i and j, to be in contact when their distance in space is smaller than the cutoff R c = 1.9 (in units of R u ), and their positions along the linear chain obey |i − j| 3.

Native structure
Reconstructed structure Figure 4. The native backbone structure of C α atoms in Rubredoxin 1e8j. In the second panel, we display a typical reconstructed structure obtained using the present algorithm.
In this way, the target structure is characterized by having 93 contacts. For each pair i, j that is in contact with the target structure, we assume that the following interaction potential V(r) is acting between them: where R a = 0.8, R r = 0.5 and ε = 0.7, independently of the type of amino acid, and r is the actual distance between monomers i, j at any instant of time. For two amino acids that are not in contact with the target structure, the interaction potential is given by equation (42) with ε = 0, i.e., we assume their interaction to be purely repulsive. With this simple scheme, we aim to reconstruct the target structure, using the thermostat rules described in the previous section, by starting from an elongated chain configuration. Here we use a mass M = 1 for all monomers, the temperature T = 0.005, the thermostat parameters τ = 26 and t 0 = 1, and the integration step τ = 10 −6 . We show in figure 5 a typical result of the reconstruction simulation by starting from an elongated chain configuration. As indicated in the figure, both the number of native and wrong contacts grow initially as a function of time. The latter remains bounded below about 30, while the true contacts are continuously formed up to a time of about 200, when a plateau is reached. This plateau is found in all simulations performed and seems to be typical of the dynamical process. After the search period is over, at t ≈ 350, the number of true contacts increases rapidly again and enters the (small) corridor delimited by the dashed lines in figure 5. These lines represent characteristic bounds for the fluctuations observed in separate simulations in which the initial configuration coincides with the target structure itself. This is done to control the stability of the relaxed structure and permits to determine the model parameters appropriately. One observes that the fluctuations in both the number of true and wrong contacts are consistent with the ones found for a relaxed structure. A typical reconstructed structure obtained using the present algorithm is displayed in figure 4.
The actual contacts taking place in a relaxation and a reconstruction simulation are shown in figure 6, where one can distinguish more clearly the different local environments displayed by the different structures. The examples shown in the figure report two cases having the same number of true and wrong contacts.

Final remarks
We have discussed a new efficient algorithm aimed at studying the long-time dynamics of polymer chains in the presence of a thermostat, by keeping the monomer-monomer link distances fixed over the whole dynamical trajectory of the system. A particular attractive extension of the algorithm presented here is to environmental-dependent couplings, which could allow the study of the long-time polymer dynamics in realistic environments. The well-known Nosé-Hoover thermostat couples all the particles of a system at all times with the same strength to a fictitious thermostat. In this manner all the parts of a system exchange energy with the thermostat at all times and at the same average rate. A polymer however, depending on its instantaneous structure, can have either its whole length or only a small fraction of its monomers exposed to the surrounding medium. As a result, the polymer dynamics is expected to depend dramatically on what fraction of its monomers can exchange energy with the surrounding medium in a given time interval. This allows for the investigation of the long timescale dynamics of polymer chains, which can not be otherwise accessible with traditional methods. The present results should also be useful in the context of protein structure predictions and in the modelling of nanoscale polymer materials.
In the following formulae, the subscripts for r k and p k take values in the interval k = 0, . . . , N, while for the quantities X k , Y k only in the interval k = 1, . . . , N and for D k,l the indices take values in the interval 1, . . . , 2N. If a subscript is outside these intervals the corresponding quantity is identically vanishing. Only the following Poisson brackets between r k , p k and the constraints X k , Y k are nonvanishing: The explicit form of the Dirac-Poisson brackets demonstrates that the phase space variables r k and p k are not canonically conjugate with respect to the Dirac symplectic structure: Here I is the 3 × 3 identity matrix and ⊗ stands for the direct product of two 3-vectors, i.e. a symmetric second-rank tensor.

Appendix B. An alternative derivation of the Lagrangian and Hamiltonian dynamics
Let us define the following 3(N + 1)-dimensional vectors: The primary and secondary constraints can then be written as follows: and the trivial scalar product in this 3(N + 1) space is defined as expected A new Lagrangian, in which the constraints are taken into account automatically can be defined as follows:Ṙ which is equivalent to considering only those coordinate variations that leave the links fixed. It is a relatively straightforward exercise to show that in this way one recovers the Lagrange equations of motion derived in the main text. The momenta are defined as The last equality holds because only such motions where δR 0 |R k = 0 are considered. It is obvious also from the above relations that and consequently momenta are not at all independent and they have to satisfy the secondary constraints. The set of secondary constraints has thus been obtained in a somewhat unusual fashion, by not allowing explicitly any dynamics in the Lagrangian L 0 in those directions, which correspond to the primary constraints. The new Hamiltonian of the system now reads as The last two equalities in equation (B.18) hold only when the secondary constraints are fulfilled, i.e., P 0 |R k = 0 for all k.

Appendix C. Coupling of a polymer chain to a thermostat
If one desires to couple the polymer with a thermostat and study its dynamics, the knowledge of the Hamiltonian dynamics becomes imperative. A way to avoid the study of the Hamiltonian dynamics in the presence of constraints is to use generalized coordinates [1,7,8]. One aspect that is often overlooked in this case is the fact that the introduction of generalized or internal coordinates makes the description very complicated. The main reason is the appearance of nonabelian gauge fields, a fact which has been brought to light only relatively recently in molecular physics [20]. The presence of nonabelian gauge fields means that in the equations of motion one has strong singularities (poles), which requires quite a sophisticated treatment (use of several different coordinate atlases), not always simple to implement and to deal with. In this respect the situation is similar to the treatment of rigid body rotations, a much simpler problem, which requires the introduction of quaternions [21] instead of the well-known Euler angles, i.e., again another form of constraints that can be dealt with in a similar way which we have chosen in the present work [22]. The treatment of holonomic constraints in a Langevin approach is a challenging problem, and accurate numerical solutions of such stochastic equations are not a trivial problem. On one hand, the presence of white noise requires in principle an infinitesimally small integration time step. On the other hand, very high frequencies should not be relevant for the physics one is usually interested in. In the following we discuss the treatment of a thermostat within the Hamiltonian dynamics by making use of the Dirac-Poisson brackets formalism. Since some of the ensuing equations are rather long if one uses standard notations, we shall introduce here a somewhat more compact one.
Let us denote the phase space variables r k , p k with k = 0, . . . , N by z α , where α = 1, . . . , 6(N + 1). The Dirac-Poisson bracket between two arbitrary functions, see equation (12), can then be written as and the dependence on z stands for dependence on all phase space variables z α . The Hamiltonian equation of motion reads aṡ When coupled to a thermostat through one pseudo-friction coefficient only the equations of motion becomeż where g(ζ) is a function of the pseudo-friction coefficient ζ, A β (z) is a vector function of the phase space variables z and κ is a constant, all of these are to be determined later. This type of coupling to a thermostat was called twisted in [11]. We are interested in generating a phase space density distribution f(z, ζ, t) satisfying the generalized Liouville equation [11] ∂f ∂t + where the symbol D/Dz α is a generalized partial derivative defined as The angle brackets stand for the scalar product of the corresponding vectors (here, the gradients of the corresponding functions) in the 6(N + 1) phase space. By definition, this generalized partial derivative satisfies the following relations: Df Dz By taking the equilibrium density distribution to be a solution of the generalized Liouville equation, we shall derive an equation of motion for the pseudo-friction coefficient ζ. We shall also require that [11] where σ is a constant, the choice of which we shall discuss later. After inserting the equilibrium density distribution into the generalized Liouville equation and using the above relations one readily obtains The second and the third term vanish identically, which can be checked by considering the generalized Liouville equation for an isolated system (in the absence of a thermostat) with an arbitrary Hamiltonian (in particular all Hamiltonians of the form H = z α ). Thus one obtainṡ (C.14) As one can see, in the final analysis for the twisted type of coupling to the thermostat all generalized partial derivatives are equivalent to ordinary partial derivatives.
We shall present now an explicit type of coupling to a thermostat, which, in spite of its simplicity, in most cases seems to lead to a very efficient thermalization of such systems. Even though we have derived the above equations by assuming one pseudo-friction coefficient only, in practice we use several pseudo-friction coefficients in order to ensure ergodic properties of this system of equations of motion and a more efficient coverage of the phase space. Quite often the simpler Nosé-Hoover coupling does not display an ergodic behaviour and the introduction of several pseudo-friction coefficients is the only cure of this defficiency (see [11]- [14]). In particular, the coupling we have described in the main text conserves the total linear momentum of the system and also the direction of the total angular momentum. Since we are interested in the internal dynamics of a polymer chain, the conservation of the total linear momentum is irrelevant, however the fact that at least the direction of the total angular momentum is conserved, might lead to nonergodic behaviour in the internal degrees of freedom, since overall rotational motion does not decouple exactly from the internal motion. For that reason, the coupling we shall consider amounts to modification of the usual Hamiltonian equations of motion given by (C.3) and (C.4) in the following form: where the last term is a 3-vector. Thus we have introduced a specific vector function A β (z).
In particular, one has to use this prescription when one determines the constants a k appearing in equations (29), (32) and (34), so as to have the constraints satisfied in the presence of the thermostat as well. Here p xk stands for the x-component of the linear momentum p k of the kth monomer. One can easily check that the total linear momentum of the chain is an integral of motion even in the presence of the thermostat, however the direction of the total angular momentum is not conserved anymore for this type of coupling, as opposed to the coupling described in the main text. We have assumed here that G(ζ x ) = ζ 2 x /2, etc. For the three pseudo-friction coefficients ζ x , ζ y , ζ z , one thus derives the following equations of motioṅ {x k , p xk } D (C. 16) and similar equations for ζ y and ζ z . An explicit expression for the Dirac-Poisson bracket {x k , p xk } D was given in appendix A, see equation (A.8). An appropriate value for the coupling constant κ is where τ is a constant (generally of order O(1)). If chosen in this way, the phase space variables r k , p k and the pseudo-friction coefficients ζ x , ζ y , ζ z will have comparable time derivatives and in this way the energy exchange between the thermostat and the polymer chain is the most effective.
One can of course consider other situations as well and it certainly makes sense to do so, however, the simulations will have to be much longer in such a case. These different types of coupling, which differ by the rate at which the energy is exchanged with the thermostat can be simulated by an appropriate choice of the coupling constants σ and κ. For the time being we shall use σ = 1 and relate further refinements to the future. In order to speed up the simulations and without a significant loss of accuracy, one can replace the last term N k=0 {x k , p xk } D by its average value 2N/3 (when the contribution of the overall translation of the chain is excluded).

Appendix D. Environment-dependent coupling of a polymer chain to a thermostat
There is a certain freedom in the coupling of a system to a thermostat, which has not been explored until now. In particular, the role of a thermostat for a polymer chain is played by solvent molecules. One might reasonably assume that the coupling of a given monomer to the thermostat should depend on how many other monomers of the chain are in the immediate vicinity of the given monomer. If the polymer chain is unfolded, all parts of the chain will be exposed to the solvent. However, after the chain has folded, either partially or fully, the monomers inside interact only very weakly with the solvent, their thermalization being thus mediated through the neighbouring monomers. One can simulate such a situation by changing the form of coupling we have described in appendix C as follows: where 0 < g(r k ) 1 is a function that depends on the environment, representing the degree of exposure of monomer k (with coordinates r k ) to the thermostat. (Implicitly g(r k ) depends on all monomer coordinates.) Indeed, we can imagine two opposite situations. On one hand, if the monomer k has only one (if it is at either end of the polymer chain) or two (otherwise) close monomer neighbours then g(r k ) ≈ 1. On the other hand, if the given monomer is surrounded by many other monomers of the chain then g(r k ) ≈ 0. Except for the appearance of the function g(r k ), the equations of motion for the phase space variable remain the same as those in the main text or those in appendix C. The equation for the pseudo-friction coefficient ζ x for example becomes noẇ Some care should be taken now in choosing appropriate values for the coupling constants κ and σ, since the average magnitude of the quantity in square brackets in the above equation for ζ x is now environment dependent. The time dependence of this quantity is controlled now by two factors: (i) thermal fluctuations, which are relatively fast changes, and (ii) environmental changes, which are relatively slow. Previously, in appendix C, we have chosen the value for the coupling constant κ such that the characteristic frequency for the pseudo-friction coefficient ζ x was comparable with the fastest mode of the polymer chain. In this case the coupling to the thermostat is most effective. It is not difficult to realize that the formalism we have described, does not actually require the coupling constants κ and σ to be constants, and can be chosen in particular as time dependent. This is what we propose to do here, namely to choose where the effective (time dependent) number of x-degrees of freedom coupled to the thermostat is and the overline stands for some type of time average over the chain configuration, consequently over the slowest mode of the polymer chain. There is no reason to compute this average exactly, only to have a fair guess of its value. In general, the described formalism should work even if the coupling constants are arbitrary (as long as the dynamics is ergodic). Only the efficiency of the simulation is affected by the actual values of these couplings, which can be varied freely, within an order of magnitude or so from their 'optimal' values, without affecting significantly the quality of a given simulation. Since we have at our disposal two 'coupling constants' κ and σ, by varying κ or τ for example, while maintaining though equation (D.3), one can still control the strength of the coupling of the polymer chain to the thermostat, i.e., the rate at which energy is exchanged between the polymer chain and the solvent. Note that the results quoted in appendix C are recovered when g(r k ) = 1, and equation (D.4) yields g eff (t) = N. The reader should have realized by now that the class of couplings to a thermostat, we have described, allow for significant freedom in modelling various situations in a very effective manner. Moreover, one can easily consider other schemes as well, e.g., coupling only part(s) of the chain to a thermostat, or even coupling each monomer to its 'own' thermostat, or increasing the number of pseudo-friction coefficients and considering other types of coupling as well, see e.g., [11,13,14], and so on.