Generic Adaptive Resolution Method for Reverse Mapping of

Sep 29, 2016 - Department of Computer Science, KU Leuven, Celestijnenlaan 200A, ... with a degree of coarse-graining ranging from two to ten heavy ato...
0 downloads 0 Views 4MB Size
Subscriber access provided by HKU Libraries

Article

A generic adaptive resolution method for reverse mapping of polymers from coarse-grained to atomistic descriptions Jakub Krajniak, Sudharsan Pandiyan, Erik Nies, and Giovanni Samaey J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00595 • Publication Date (Web): 29 Sep 2016 Downloaded from http://pubs.acs.org on September 30, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A generic adaptive resolution method for reverse mapping of polymers from coarse-grained to atomistic descriptions Jakub Krajniak,† Sudharsan Pandiyan,‡ Eric Nies,∗,‡ and Giovanni Samaey∗,† †KU Leuven Department of Computer Science, Celestijnenlaan 200A, 3001 Leuven, Belgium ‡KU Leuven Division of Polymer Chemistry and Materials, Department of Chemistry, Celestijnenlaan 200F, 3001 Leuven, Belgium E-mail: [email protected]; [email protected] Abstract In this paper, we propose a new generic approach for reverse mapping from coarsegrained to atomistic scale based on the adaptive resolution scheme (AdResS). In AdResS simulation, two spatial domains, modelled at two different scales, are brought together in a concurrent simulation by defining a hybrid region where particles can switch representation from one model to another. We use AdResS as a central part of a reverse mapping algorithm from a different perspective by treating the whole simulation box as a hybrid region and changing the resolution as a function of time during the course of a molecular dynamics simulation. The proposed method depends only on a single parameter that controls the reverse mapping process and it is independent of atomistic and coarse-grained force-fields. We use our method to reconstruct a representative atomistic configuration from a coarse-grained representation of dodecane and polyethylene liquids with degree of coarse-graining ranging from two to ten heavy

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

atoms. The conformational and dynamical properties of the reconstructed systems are in excellent agreement with the reference atomistic simulation.

1

Introduction

Polymers play a crucial role in our day-to-day life due to their variety of applications. Atomistic simulations have been applied successfully to many polymer systems and the results are generally consistent with experimental data. However, the number of particles used to describe complex polymer melts grows rapidly, hence the computationally feasible time and length scales are limited. One of the solutions to overcome this problem is to reduce the number of particles. This is done via a coarse-grained (CG) approach, which not only leads to a smaller number of particles (so-called CG beads), but also allows using a larger time-step because the fast degrees of freedom (e.g., bond length vibrations) are integrated out in CG effective potentials. 1–3 There are multiple techniques to obtain CG effective potentials that preserve conformational properties, 4 match the multi-body potentials of the mean force 5 or reproduce the thermodynamic properties. 6 In the examples for this manuscript, we choose most commonly used the Iterative Boltzmann inversion (IBI) method that preserve the radial distribution function, although the reverse mapping procedure is independent of the CG force-field. The idea of using a coarse-grained approach is already utilized in a wide range of computational experiments, e.g. CG models of polymers 7,8 or more sophisticated examples of polyelectrolytes. 9,10 Some polymer properties can already be analysed from CG simulations. 11–13 However, many other properties strictly depend on the presence of fine-scale degrees of freedom. 2,14 For instance, the dynamics of the molecules are affected by the soft coarse-grained potentials hence the calculated diffusion coefficients does not corresponds to the real values and scaling factors are required to match with the experimental values. 15 The stress vs. strain behaviour of polymeric materials are studied at the CG level 16–19 however,

2

ACS Paragon Plus Environment

Page 2 of 42

Page 3 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

to obtain experimentally comparative results the CG potentials have to be re-optimized. This is because the intrinsic conformational properties of polymers contribute to the overall mechanical response. 20 Another example is water transport through materials, which sensitively depends on the presence of hydrogen bonds. 21 Therefore, it would be beneficial to have a generic tool that will transform the CG representation back into the atomistic one. Once a particular coarse-grained technique is chosen, the mapping of the atomistic particle positions to the CG representative positions has a unique solution and it is mostly nothing more than calculating centres of masses (COMs). The tools to achieve that are included into existing software packages that are used to determine the CG effective potentials e.g. VOTCA 22 or OCTA. 23 Conversely, the reverse mapping from CG to a fully atomistic scale has no unique solution and therefore various methods for reconstructing atomistic details from a CG structure have been proposed in the literature. 24–30 We highlight two of them that should be applicable to a wide range of molecules. Rzepiela et al. 29 proposed a method where the atomistic particles belonging to a CG bead are constrained around the centre of mass (COM) by harmonic bonds that are released during a simulated annealing procedure. The method uses the modified GROMACS 31 (v3.3) code; the procedure requires several parameters to be set up, such as an initial value for the maximal force, the rate at which the restraints are released, the restraint force constant and the initial temperature and time for the annealing procedure. The authors showed correctness of going back from the CG scale (described with the MARTINI force-fields 6 ) to the atomistic scale for various systems by analysing the potential energy changes and the root mean squared displacement of resulting structure compared to the target one. Another approach, proposed by Wassenaar et al. 30 uses several short cycles of energy minimization and molecular dynamics (MD) runs with different time steps to obtain an equilibrated atomistic configuration. Also in this case, the MARTINI force-field was used at the CG scale and the initial position of the atomistic particles was put randomly around the CG COM. Although both methods could, in principle, be applied

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to any system, they require from users to set several parameters that control the reverse mapping procedure. Moreover, to support other atomistic and CG force-fields, the available code has to be extended. Recent years have brought a few hybrid schemes that allow performing concurrent multiscale simulation. 32–37 The common idea behind these methods is to define a spatial region where CG or atomistic potentials are used to describe interactions. The CG region is larger than the atomistic one hence the total number of pair interaction is reduced with respect to the atomistic simulation. This is because in the CG region only pair interaction between centres of mass are calculated. Between those regions, there is a hybrid part where the force 32 or potential interpolation 34,38 takes place. In this sense, a molecule that travels from the CG to the fine-scale part is effectively reverse mapped in the hybrid region where both CG and fine-scale potential are presented. In this paper, we propose a generic reverse mapping technique that requires a minimum set of parameters. We extend the idea of a hybrid region to the whole simulation box and propose that the representation will change in time and not in space. Recently, B¨ockemann et al. 39 use a very similar approach to perform the transition between ab initio (QM) and force-field (MM) description. In their approach, the scaling parameter was an additional fictitious particle that evolves in time (as in the Nos´e-Hoover thermostat 40 ). The aim of this work is to construct a generic reverse mapping method that is independent of the specific CG and atomistic force-fields used and moreover requires only a single parameter to control the process of mapping the representation from the CG to fine-scale level. We use the existing implementation of the force based adaptive resolution method (AdResS) 32 as a starting point. The remainder of this article is organized as follows: in section 2, we describe the method. Section 3 contains a number of simulations to illustrate the performance and versatility of the presented method. In Section 4, we use this method to obtain representative configurations of atomistic dodecane and polyethylene (from two different CG mapping schemes) liquids. We conclude in Section 5. The software package

4

ACS Paragon Plus Environment

Page 4 of 42

Page 5 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

created during this work is publicly available 41 and attached to the supportive information.

2

Method

This section describes the reverse mapping method that reconstructs atomistic details (AT) from the underlying coarse-grained (CG) configuration and is the main contribution of the present paper. The procedure requires that the input coarse-grained snapshot represents an equilibrated state and that all CG potentials are available. Moreover, we assume that the mapping between the CG bead and the corresponding set of atoms is available. We use the term molecule to describe a list of beads that together form a single physical entity (e.g., a polymer chain). The whole procedure is divided into two parts (see figure 1). Section 2.1 describes the first part, which has the aim of preparing a hybrid AT-CG configuration that contains both coarse-grained and atomistic particles and serves as an initial condition for the second part. The second part, described in Section 2.2, is responsible for reintroducing correct positions of the atomistic particles during a dedicated molecular dynamics (MD) simulation. This step is required because the replicated fragments from the initial conditions are in a random orientation with respect to each other; therefore, the molecule that those fragments represent is not in an equilibrated state. At the end of the work flow, we get atomistic coordinates that later on can be used in a standard molecular simulation at the atomistic level.

2.1

Preparation of hybrid coordinates

Starting from a pure CG system, we introduce a hybrid system that is composed of the position coordinates of both CG beads and atomistic particles. In this section, we describe the method for only one type of CG beads but the software package supports multiple CG bead types. The CG system is described as a list of M beads divided into J molecules and each of the

5

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 6 of 42

Page 7 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Algorithm 1 Preparing of the hybrid coordinates Input: the representative atomistic molecule K the list of M coarse-grained coordinates {Cij } the mapping function Q ⊲ returns corresponding atomistic fragment for CG bead i Output: output hybrid configuration as a list of the CG and the atomistic coordinates {O} 1: 2: 3: 4: 5: 6: 7: 8: 9:

{O} ← empty output list of coordinates for i ← 1 to M do {c} ← Q(i) O.insert(R[i]) for all p ∈ {c} do P ← p + R[i] O.insert(P ) end for end for

2.2

Reintroducing of atomistic details

In this section, we describe the basis of the reverse mapping method. We build the method on top of AdResS, which combines the CG and the atomistic representation by dividing the simulation box into three different spatial regions. Particles in the first part are simulated with CG potential and particles in the second part are simulated with an atomistic potential. The exchange of particles between those two region is done in a third hybrid part where the mixing of two potentials is governed by a position dependent switching function. The details of the method can be found in Praprotnik et al. 32 In our method, we treat the whole simulation box as a hybrid region and we use a time-dependent switching function, defined as

λ(t) = λ0 + α(t − t0 ), {λ ∈ R; 0 ≥ λ ≥ 1}

(1)

where λ0 is the initial resolution, α is a switching rate (s−1 ) and t0 is the simulation time at 7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 8 of 42

Page 9 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

the potential energy and consequently the Hamiltonian is invariant under the translation. However, the time-dependent λ function in the total potential energy violates the energy conservation. This is remedied by the work done of thermostat that removes the excess of heat produced during the reverse mapping procedure. Second, since there are no regions in space at different resolutions, density imbalances also cannot occur. 42,43 For λ = 0, the positions of the CG and atomistic particles are driven by forces deAT rived from the coarse-grained terms and VIII . In this way, we preserve the conformation of

the molecules on the coarse-grained level and also we keep valid position of the atomistic fragments inside the coarse-grained beads. During the transition time, t > t0 , the forces obtained from the coarse-grained terms are gradually switched off (both bonded and nonbonded terms). At the same time the forces acting between atomistic particles as well as the atomistic non-bonded terms are gradually switched on. At λ = 1 the system is reverse mapped to the fine scale representation. Evidently, the reverse mapping procedure can be also used to switch from the atomistic to the CG representation of the system and the current implementation allows switching back and forth during the simulation from coarse-grained to atomistic representation.

3

Simulation details

We performed the reverse mapping of three different systems, a simple molecule (dodecane), a polymer chain (polyethylene) and a ring molecule (trimethylol melamine). In case of dodecane, the CG beads were composed of two united-atom CH2 segments. For the polyethylene systems, three different coarse-grained representations were simulated. In the first case, the same scheme as for dodecane was used, in the second case four united-atom CH2 segments were placed in the CG bead. For the last case, ten carbon atoms with corresponding hydrogens were placed in the CG bead and for this case the complete all-atom description was used. The 27 atoms forming the ring molecule were decomposed into three CG beads. We

9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 42

used GROMOS 53a6, 44 OPLS/UA 45 and OPLS/AA 46 for building atomistic models (see table 1). We prepared a single molecule atomistic representation for all the test cases using the molecular drawing tool Avogadro 47 and the initial atomistic configuration was then used to prepare the bulk samples in Packmol. 48 The topology files for GROMACS were prepared using MKTOP 49 and verified manually. All the simulations were carried out with GROMACS (5.0.4). 50 The CG potentials were prepared using IBI method available in the VOTCA package (1.3). 22 A modified version of ESPResSo++ 51 was used to run the reverse mapping procedure. The simulation data were stored in the H5MD format 52 and were analysed by our in-house code. 53 The data analysis work-flow was integrated within IPython 54 environment and the plots were made in matplotlib. 55 Table 1: Details about simulated systems where N is the number of particles, M is the number of molecules and T is the temperature. In square bracket we report the level of coarse-graining where 2:1 means that two atomistic particles are put in a single coarsegrained bead. System Dodecane [2:1] Polyethylene [2:1] Polyethylene [4:1] Polyethylene [10:1] Trimethylol melamine [9:1]

N 6000 7500 7500 22650 13500

M 500 75 75 75 500

T (K) 298 423 423 423 800

Density (kg/m3 ) 743.4 814.7 816.1 733.04 643.8

Force field GROMOS 53a6 OPLS/UA OPLS/UA OPLS/AA OPLS/AA

L (nm) 5.75 6.02 6.00 6.2 6.53

The atomistic systems were simulated with periodic boundary conditions in a cubic box. The non-bonded interactions were cut-off at 1.4 nm for GROMOS 53a6 and at 1.2 nm for OPLS/UA and OPLS/AA force-fields. The Verlet list algorithm 56 was used for updating neighbour lists every time-step with a skin parameter set to 0.16 nm. A Nos´e-Hoover thermostat 40 with coupling constant of 0.5 ps was used and the temperature set according to the table 1. The Berendsen 57 barostat with coupling constant 0.5 ps and the pressure of 1 bar was used. The time step was set to ∆t = 0.001 ps. The two last examples (polyethylene: case 3 and trimethylol melamine) possess partial charges hence the Coulombic interactions 10

ACS Paragon Plus Environment

Page 11 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

were computed by the reaction-field method 58 with the cut-off set to 0.9 nm. The initial input configuration was equilibrated by energy minimisation followed by 10 ns of a NPT simulation. Afterwards 10 ns of NVT simulation was performed and the coordinates were collected every 0.5 ps. We used the Inverse Boltzmann method for bonded potentials and the IBI method to obtain non-bonded potentials. The details of the final densities and the box sizes L are shown in table 1. The CG simulation for IBI method were governed by a Langevin thermostat, 59,60 with coupling constant γ = 5.0 ps and the cut-off distance for non-bonded interactions was set to 1.4 nm. The remaining settings were as in the atomistic simulations. The CG force-field was obtained after several iterations and validated by comparing the radial distribution functions between the reference atomistic system and the coarse-grained system. The CG force-field as well as the parameters of the atomistic force field are attached to the supportive information. After obtaining the CG force-field, a simulation at the CG level was carried out. The last frame of the CG simulation was used to prepare the hybrid configuration using the procedure given in section 2. The reverse mapping simulation ran for ∆t/α ps and at the end of the procedure we obtained the atomistic representation. This was followed by 10 ns of NVT simulation with an atomistic resolution (λ = 1) that was used to calculate properties after the reverse mapping and to compare with the reference atomistic simulations. Figure 3 shows the complete workflow.

4

Results

We investigate the role of the reverse mapping parameter α (eq. 1) on the conformational properties such as the radial distribution function (RDF), the bond length, the bend angle, torsional angle probability distribution function (PDFs) and the average volumetric mass density. The change in total potential energy is reported in the reduced unit per segment Ep⋆ = Ep /N RT where N is number of segments, R is the gas constant and T is the absolute

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 12 of 42

Page 13 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

temperature. The RDF is calculated for the systems between non-bonded CH2 units and in the case of all-atom examples (polyethylene case 3 and trimethylol melamine), between the atomistic particles. The standard error for the volumetric mass density is calculated by the block averaging method. 61 We compare probability distribution functions of the radius of gyration, where Rg is defined by

Rg =

*

K 1 X mk krk − rcom k2 M k

! 12 +

where M is the total mass of the chain K, rcom is the centre of mass of the polymer chain and rk is the position of atom k with mass mk . The end-to-end distance REE is defined as REE = kr100 − r1 k where r are the positions of atoms. The characteristic ratio 62 is defined as Cs = hRs2 i/(sb20 )

(3)

where b0 is an average bond length and < Rs2 > is an average squared distance between atoms separated by s bonds along the polymer chain. The characteristic ratio is sensitive to abnormal conformational changes in the polymer chain and it is successfully used to validate equilibrating methods of long polymer chains. 63–65 We calculated it only for polyethylene examples because the dodecane molecules are too short to get any insight from this measurement. We assess the difference between PDFs distributions in two ways, by the square difference test function, which is often use to compare distributions during the iterative methods and by the Kolmogorov-Smirnov two-sample (KS) test. 66,67 The squared difference test function f is defined by: f=

N X

2 (hαi − hAT i )

(4)

i

where hα and hAT are the histograms to compare, the later one is the reference histogram and N is the number of bins. In the KS test, the difference between two empirical distribution

13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 42

functions (edf) is calculated by: d(eAT , eα ) = max |eAT (x) − eα (x)| x

(5)

where eAT is an empirical cumulative distribution of a quantity from the reference simulation and eα is from the simulation after reverse mapping with given α parameter. We computed dynamical properties of the polymer in terms of the mean squared displacement (MSD) for all segments in the chain (with subtraction of the system COM movement) and the orientational autocorrelation function (OACF) 68 of end-to-end vectors between the first and the last segment in the molecules. The MSD has been calculated by:

M SD(τd ) =

TX −τd X N 1 (Ri (t) − Ri (t + τd ))2 N (T − τd ) t=0 i=1

(6)

where T is the number of time frames, τd is the time frame gap, N is the number of particles and R represents the position of the polymer segment. The OACF is defined similarly by: TX −τd X J 1 EE REE CEE (τd ) = j (t) · Rj (t + τd ) J(T − τd ) t=0 j=0

(7)

where T is the number of time frames, τd is a time frame gap, J is the number of polymer chains and REE represents the end-to-end vector of j-th chain. We use a linear least squares j method 69 of 6Dt = hM SD(t)i

(8)

to get diffusion coefficients (D) from the diffusive regime and the OACF to obtain the relaxation time τR by fitting CEE (t) to a linear form of exponential function:

CEE (t) = C0 e−t/τR

where C0 is a constant. 14

ACS Paragon Plus Environment

(9)

Page 15 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 16 of 42

Page 17 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

squares method 69 of equation 8 to get D: 0.37 ± 0.02 nm2 /ns (AT), 0.37 ± 0.01 nm2 /ns (1 ps), 0.37 ± 0.01 nm2 /ns (2 ps), 0.36 ± 0.01 nm2 /ns (20 ps) where numbers in brackets denotes ∆t/α. We use the same method to obtain the relaxation time τR from a linear form of equation 9. The values of τR are: 0.29 ± 0.01 ns (AT), 0.30 ± 0.01 ns (1 ps), 0.29 ± 0.01 ns (2 ps), 0.29 ± 0.01 ns (20 ps). The values, both for diffusions coefficients and relaxation times, are in a very good agreement with respect to values from the atomistic simulation. In addition to conformational and dynamical properties, we also ran 10 ns of NPT simulation to compare the average volumetric mass densities of the simulated systems. The density of our reference atomistic system was 744.05 ± 0.16 kg/m3 and the densities for the different values of α (20 ps, 2 ps, 1 ps), are: 744.19 ± 0.15 kg/m3 , 744.24 ± 0.15 kg/m3 , 744.03 ± 0.16 kg/m3 . Here, we also do not observe any influence of the reverse mapping on the average densities.

4.2

Polyethylene: case 1

For the next two examples, we study a polyethylene (PE) melt that consists of 75 PE chains, each contains 100 CH2 units. We discuss three CG levels; namely with two and four polymer segments grouped into a single CG bead. As in the case of dodecane, we compare conformational and dynamical properties to assess the correctness of the method with respect to α: 10−5 , 5 · 10−5 , 10−4 , 5 · 10−4 , 10−3 that correspond to the reverse mapping times: 100 ps, 20 ps, 10 ps, 2 ps and 1 ps (we refer later to these values).

In figure 7 the change in the reduced total potential energy Ep⋆ is

plotted. As for dodecane, the reverse mapping procedure starts at t0 = 20 ps. Here we do not observe a rapid increase of the energy after the procedure started and only for the highest value of α there is a peak above the reference value (purple line in figure 7). However, the energy barrier between the CG and atomistic state is larger than in the case of dodecane 2.3. The values of Ep⋆ at the end of the reverse mapping procedure (at t > t0 + ∆t/α) with respect to the α parameters are: 3.4 ± 0.2 (1 ps), 3.4 ± 0.2 (2 ps), 3.4 ± 0.2 (20 ps), 3.4 ± 0.2 18

ACS Paragon Plus Environment

Page 18 of 42

Page 19 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 20 of 42

Page 21 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 22 of 42

Page 23 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 24 of 42

Page 25 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.4

Polyethylene: case 3

With the last example of polyethylene, we show that even for bigger CG beads, the reverse mapping procedure works very well and the hybrid configuration is shown in figure 15. The

Figure 15: The visualization of the hybrid configuration for the third case of polyethylene at the end of the reverse map procedure. The whole simulation box contains of 75 chains, here we only show a single chain. coarse-grain force-field was obtained with the IBI method as in all previous examples. We follow the same protocol that is depicted in figure 3 and described in section 3. In figure 16 the change in the reduced total potential energy Ep⋆ is plotted. Here, as in all previous examples, we start from the coarse-grained configuration after 10 ns of simulation and then use the reverse mapping method to go back to fine grain state and the first 5 ns simulation with the coarse-grained resolution is only for the comparison. The reduced total potential energy Ep⋆ after the reverse mapping (at t > t0 + ∆t/α) is 0.12 ± 0.03 and the reference value is 0.12 ± 0.02. We study in detail the conformational properties of polyethylene. Here, because of all atom description, we compare all available conformational distributions. Following the results depicted in figure 17, we do not observe any significant difference between the distributions from the reference atomistic simulation and the one obtained after reverse mapping. Next, we compare the dynamical properties of the system, the mean squared displacement and the rotational autocorrelation function. The results are shown in figure 18. The values of the diffusion coefficient D are 0.11 ± 0.01 nm2 /ns for the reference simulation and 0.11 ± 0.01 nm2 /ns for the reversed mapped system. The 26

ACS Paragon Plus Environment

Page 26 of 42

Page 27 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 28 of 42

Page 29 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

relaxation times τR of OACF are: 12.10 ± 0.04 ns for reference system and 12.00 ± 0.03 ns for the reversed mapped system. Both, the configurational properties and the dynamical are not affected by the reverse mapping method as well as the average value of radius of gyration are: 1.93 ± 0.01 nm (AT) and 1.96 ± 0.01 nm (100 ps), the values of end-to-end distance are 3.01 ± 0.09 nm (AT) and 2.98 ± 0.09 nm (100 ps). The KS and the squared difference tests did not indicate any significant differences between obtained distributions. To conclude, this example demonstrates that the reverse mapping method presents in this work is capable to work with high degree of coarse-graining where ten heavy atoms are grouped into single CG bead. The analysis of conformational as well as dynamical measurements do not show any significant difference between the reference simulation and the one obtained after reverse mapping.

4.5

Trimethylol melamine

In this example, we reverse mapped a system composed of 500 ring molecules. The coarsegrain system was built by mapping 27 atoms into three CG beads (see figure 19). For evaluating the procedure, we followed the same protocol that is shown in figure 3 and described in section 3. The change in the total reduced potential energy is shown in figure 20. The average value of Ep⋆ after the reverse mapping (at t > t0 + ∆t/α) is 12.72 ± 0.11 whereas the average value of the reference simulation is 12.66±0.13. We examined the conformational properties of the molecules by comparing RDFs between atom pairs C − C, C − N and O − H. The results are shown in figure 21. All of the distributions are recovered properly, the visual inspection of the molecules after the reverse mapping confirms that the ring structure is correctly put back. The KS and squared difference tests do not show any significant difference between the distributions obtained after reverse mapping and the reference simulation. In the last two calculations, we check the MSD and the average number of hydrogen bonds 29

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 30 of 42

Page 31 of 42

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The dynamical and configurational properties are very well reproduced. The average number of hydrogen bonds that can be formed in this system are in very good agreement with the reference simulation. With this last example we demonstrate that the method is capable to reverse map the ring molecule.

5

Conclusions and future work

With this work, we proposed a generic reverse mapping algorithm that can be controlled only by a single parameter and that is independent of fine-scale and coarse-grained force-fields. We have used this method on two different representative polymer chains with different levels of coarse-graining and on a ring molecule. We investigate the effect of reverse mapping on the final structure by examined the conformational and dynamical properties and we have not observed any significant difference between the results from the fine-scale simulations and from the simulations after reverse mapping. Moreover, we have not encountered any systematic influence of the α parameter on the final structures. As a result of this work, we deliver a software package with the examples used to prepare this manuscript as well as an extension to ESPResSo++ that allows easily perform the reverse mapping. While the method we propose here is very generic and require little tuning of parameters, future work will be necessary to apply the method to more complex situations. For instance, the rising importance of multiscale simulation of highly cross-linked polymer networks brings new challenges for reverse mapping. The main difficulties arise from the variable connectivity between the CG beads, which has to be reflected by the different atomistic representations. Because of its modular design, the method presented in this work can still be used for such systems by modifying only the first (preparation) step. This new application is the focus of our current investigation.

32

ACS Paragon Plus Environment

Page 32 of 42

Page 33 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Acknowledgement This work was supported by the “Strategic Initiative Materials” in Flanders (SIM) under the InterPoCo program. The computational resources used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Hercules Foundation and the Flemish Government – department EWI. The authors also wish to thank the anonymous reviewers for their helpful suggestions.

Supporting Information Available The following files are available: • krajniak generic reverse mapping si.pdf: The details of the force-fields parameters as well as the protocol with the list of commands to run the code. • bakery-v1.3.zip: The software package to perform the reverse mapping with the examples. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Li, C.; Strachan, A. Molecular scale simulations on thermoset polymers: A review. J. Polym. Sci., Part B: Polym. Phys. 2014, 103 – 122. (2) Noid, W. G. Perspective: Coarse-grained models for biomolecular systems. J. Chem. Phys. 2013, 139, 090901 – 090926. (3) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J. Molecular dynamics algorithm for multiple time scales: Systems with long range forces. J. Chem. Phys. 1991, 94, 6811 – 6815. 33

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 42

(4) Reith, D.; P¨ utz, M.; M¨ uller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. J. Comput. Chem. 2003, 24, 1624 – 1636. (5) Lyubartsev, A.; Laaksonen, A. Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach. Phys. Rev. E 1995, 52, 3730 – 3737. (6) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; Vries, A. H. D. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812 – 7824. (7) Woloszczuk, S.; Mineart, K. P.; Spontak, R. J.; Banaszak, M. Dual modes of selfassembly in superstrongly segregated bicomponent triblock copolymer melts. Phys. Rev. E 2015, 91, 010601 – 010607. (8) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A moleculardynamics simulation. J. Chem. Phys. 1990, 92, 5057 – 5086. (9) V¨ogele,

M.;

electrolyte

Holm,

complexes:

C.;

Smiatek, MARTINI

J.

Coarse-grained

models

for

simulations

poly(styrene

of

sulfonate)

polyand

poly(diallyldimethylammonium). J. Chem. Phys. 2015, 143, 243151 – 243159. (10) Uusitalo, J. J.; Ingolfsson, H. I.; Akhshi, P.; Tieleman, D. P.; Marrink, S. J. Martini Coarse-Grained Force Field: Extension to DNA. J. Chem. Theory Comput. 2015, 11, 3932 – 3945. (11) Li, T.; Jiang, Z.; Yan, D.; Nies, E. A polyethylene chain investigated with replica exchange molecular dynamics simulation: Equilibrium lamellar thickness and melting point, ordering and free energy. Polymer 2010, 51, 5612 – 5622. (12) Salerno, K. M.; Agrawal, A.; Perahia, D.; Grest, G. S. Resolving Dynamic Properties

34

ACS Paragon Plus Environment

Page 35 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

of Polymers through Coarse-Grained Computational Studies. Phys. Rev. Lett. 2016, 116, 058302 – 058307. (13) Wu, R.; Li, T.; Nies, E. Langevin Dynamics Simulation of Chain Crosslinking into Polymer Networks. Macromol. Theory Simul. 2012, 21, 250 – 265. (14) Dunn, N. J. H.; Noid, W. G. Bottom-up coarse-grained models that accurately describe the structure, pressure, and compressibility of molecular liquids. J. Chem. Phys. 2015, 143, 243148 – 243164. (15) Pandiyan, S.; Parandekar, P. V.; Prakash, O.; Tsotsis, T. K.; Basu, S. Systematic Coarse Graining of a High-Performance Polyimide. Macromol. Theory Simul. 2015, 24, 513 – 520. (16) Jang, C.; Abrams, C. Thermal and mechanical properties of thermosetting polymers using coarse-grained simulation. Eur. Phys. J.: Spec. Top. 2016, 1 – 9. (17) Tsige, M.; Lorenz, C. D.; ; Stevens, M. J. Role of Network Connectivity on the Mechanical Properties of Highly Cross-Linked Polymers. Macromolecules 2004, 37, 8466 – 8472. (18) Sliozberg, Y. R.; Mrozek, R. A.; Schieber, J. D.; Kr¨oger, M.; Lenhart, J. L.; Andzelm, J. W. Effect of polymer solvent on the mechanical properties of entangled polymer gels: Coarse-grained molecular simulation. Polymer 2013, 54, 2555 – 2564. (19) Venkatesan, S.; Basu, S. Investigations into crazing in glassy amorphous polymers through molecular dynamics simulations. J. Mech. Phys. Solids 2015, 77, 123 – 145. (20) Li, C.; Strachan, A. Molecular simulations of crosslinking process of thermosetting polymers. Polymer 2010, 51, 6058 – 6070. (21) Pandiyan, S.; Krajniak, J.; Samaey, G.; Roose, D.; Nies, E. A molecular dynamics

35

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

study of water transport inside an epoxy polymer matrix. Comput. Mater. Sci. 2015, 106, 29 – 37. (22) R¨ uhle, V.; Junghans, C.; Lukyanov, A.; Kremer, K.; Andrienko, D. Versatile objectoriented toolkit for coarse-graining applications. J. Chem. Theory Comput. 2009, 5, 3211 – 3223. (23) Aoyagi, T.; Sawa, F.; Shoji, T.; Fukunaga, H.; Takimoto, J.-i.; Doi, M. A generalpurpose coarse-grained molecular dynamics program. Comput. Phys. Commun. 2002, 145, 267 – 279. (24) Brocos, P.; Mendoza-Espinosa, P.; Castillo, R.; Mas-Oliva, J.; Pi˜ neiro, A. Multiscale molecular dynamics simulations of micelles: coarse-grain for self-assembly and atomic resolution for finer details. Soft Matter 2012, 8, 9005 – 9014. (25) Peter, C.; Kremer, K. Multiscale simulation of soft matter systems – from the atomistic to the coarse-grained level and back. Soft Matter 2009, 5, 4357 – 4366. (26) Nielsen, S. O.; Ensing, B.; Moore, P. B.; Klein, M. L. Multiscale Simulation Methods for Nanomaterials; John Wiley & Sons, Inc., 2007; pp 73–88. (27) Ghanbari, A.; B¨ohm, M. C.; M¨ uller-Plathe, F. A simple reverse mapping procedure for coarse-grained polymer models with rigid side groups. Macromolecules 2011, 44, 5520 – 5526. (28) Brasiello, A.; Crescitelli, S.; Milano, G. A multiscale approach to triglycerides simulations: from atomistic to coarse-grained models and back. Faraday Discuss. 2012, 158, 479 – 492. (29) Rzepiela, A. J.; Sch¨afer, L. V.; Goga, N.; Risselada, H. J.; De Vries, A. H.; Marrink, S. J. Reconstruction of atomistic details from coarse-grained structures. J. Comput. Chem. 2010, 31, 1333 – 43. 36

ACS Paragon Plus Environment

Page 36 of 42

Page 37 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(30) Wassenaar, T. A.; Pluhackova, K.; B¨ockmann, R. A.; Marrink, S. J.; Tieleman, D. P. Going Backward: A Flexible Geometric Approach to Reverse Transformation from Coarse Grained to Atomistic Models. J. Chem. Theory Comput. 2014, 10, 676 – 690. (31) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701 – 1718. (32) Praprotnik, M.; Delle Site, L.; Kremer, K. Adaptive resolution molecular-dynamics simulation: changing the degrees of freedom on the fly. J. Chem. Phys. 2005, 123, 224106 – 224120. (33) Ensing, B.; Nielsen, S. O.; Moore, P. B.; Klein, M. L.; Parrinello, M. Energy Conservation in Adaptive Hybrid Atomistic/Coarse-Grain Molecular Dynamics. J. Chem. Theory Comput. 2007, 3, 1100 – 1105. (34) Nielsen, S. O.; Moore, P. B.; Ensing, B. Adaptive multiscale molecular dynamics of macromolecular fluids. Phys. Rev. Lett. 2010, 105, 237802 – 237806. (35) Potestio, R.; Espa˜ nol, P.; Delgado-Buscalioni, R.; Everaers, R.; Kremer, K.; Donadio, D. Monte Carlo Adaptive Resolution Simulation of Multicomponent Molecular Liquids. Phys. Rev. Lett. 2013, 111, 060601 – 060606. (36) Shen, L.; Hu, H. Resolution-adapted all-atomic and coarse-grained model for biomolecular simulations. J. Chem. Theory Comput. 2014, 10, 2528 – 2536. (37) Pezeshki, S.; Lin, H. Recent developments in QM/MM methods towards open-boundary multi-scale simulations. Mol. Simul. 2015, 41, 168 – 189. (38) Potestio, R.; Fritsch, S.; Espa˜ nol, P.; Delgado-Buscalioni, R.; Kremer, K.; Everaers, R.; Donadio, D. Hamiltonian Adaptive Resolution Simulation for Molecular Liquids. Phys. Rev. Lett. 2013, 110, 108301 – 108306.

37

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 42

(39) B¨ockmann, M.; Doltsinis, N. L.; Marx, D. Adaptive switching of interaction potentials in the time domain: An extended Lagrangian approach tailored to transmute force field to QM/MM simulations and back. J. Chem. Theory Comput. 2015, 11, 2429 – 2439. (40) Evans, D. J.; Holian, B. L. The Nose–Hoover thermostat. J. Chem. Phys. 1985, 83, 4069 – 4074. (41) Krajniak,

J.

bakery:

the

reverse

mapping

tool

v1.3.

2016;

http://dx.doi.org/10.5281/zenodo.61233. (42) Potestio, R.; Peter, C.; Kremer, K. Computer Simulations of Soft Matter: Linking the Scales. Entropy 2014, 16, 4199. (43) Delle Site, L. Some fundamental problems for an energy-conserving adaptive-resolution molecular dynamics scheme. Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 2007, 76, 047701. (44) Daura, X.; Mark, A. E.; Van Gunsteren, W. F. Parametrization of aliphatic CHn united atoms of GROMOS96 force field. J. Comput. Chem. 1998, 19, 535 – 547. (45) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638 – 6646. (46) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225 – 11236. (47) Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. Cheminf. 2012, 4, 1 – 17. (48) Mart´ınez, L.; Andrade, R.; Birgin, E. G.; Mart´ınez, J. M. PACKMOL: A package for

38

ACS Paragon Plus Environment

Page 39 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157 – 2164. (49) Ribeiro, A.; Horta, B. A. C.; Alencastro, R. B. d. MKTOP: a program for automatic construction of molecular topologies. J. Braz. Chem. Soc. 2008, 19, 1433 – 1435. (50) Pronk, S.; P´all, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; Van Der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845 – 854. (51) Halverson, J. D.; Brandes, T.; Lenz, O.; Arnold, A.; Bevc, S.; Starchenko, V.; Kremer, K.; Stuehn, T.; Reith, D. ESPResSo++: A modern multiscale simulation package for soft matter systems. Comput. Phys. Commun. 2013, 184, 1129 – 1149. (52) de Buyl, P.; Colberg, P. H.; H¨ofling, F. H5MD: A structured, efficient, and portable file format for molecular data. Comput. Phys. Commun. 2014, 185, 1546 – 1553. (53) Krajniak, J. lab-tools: Set of Python scripts for analysing molecular dynamics simulations of polymers: release v1.3. 2016, http://dx.doi.org/10.5281/zenodo.59804. (54) P´erez, F.; Granger, B. E. IPython: a System for Interactive Scientific Computing. Comput. Sci. Eng. 2007, 9, 21 – 29. (55) Hunter, J. D. Matplotlib: A 2D graphics environment. Computing In Science & Engineering 2007, 9, 90 – 95. (56) Loup, V. Computer ”Experiments” on Classical Fluis. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98 – 103. (57) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684 – 3690. 39

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(58) Barker, J.; Watts, R. Monte Carlo studies of the dielectric properties of water-like models. Mol. Phys. 1973, 26, 789 – 792. (59) Goga, N.; Rzepiela, A. J.; de Vries, A. H.; Marrink, S. J.; Berendsen, H. J. C. Efficient Algorithms for Langevin and DPD Dynamics. J. Chem. Theory Comput. 2012, 8, 3637 – 3649. (60) Espa˜ nol, P. In Novel Methods in Soft Matter Simulations; Karttunen, M., Lukkarinen, A., Vattulainen, I., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2004; Chapter Statistica, pp 69 – 115. (61) Flyvbjerg, H.; Petersen, H. G. Error estimates on averages of correlated data. J. Chem. Phys. 1989, 91, 461 – 466. (62) Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G.; Theodorou, D. N. Atomistic Monte Carlo simulation of strictly monodisperse long polyethylene melts through a generalized chain bridging algorithm. J. Chem. Phys. 2002, 117, 5465 – 5479. (63) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. Equilibration of long chain polymer melts in computer simulations. J. Chem. Phys. 2003, 119, 12718 – 12728. (64) Moreira, L. A.; Zhang, G.; M¨ uller, F.; Stuehn, T.; Kremer, K. Direct Equilibration and Characterization of Polymer Melts for Computer Simulations. Macromol. Theory Simul. 2015, 24, 419 – 431. (65) Harmandaris, V.; Reith, D.; van der Vegt, N. F. a.; Kremer, K. Comparison Between Coarse-Graining Models for Polymer Systems: Two Mapping Schemes for Polystyrene. Macromol. Chem. Phys. 2007, 208, 2109 – 2120. (66) Stephens, M. a. EDF Statistics for Goodness of Fit and Some Comparisons. J. Am. Stat. Assoc. 1974, 69, 730 – 737. 40

ACS Paragon Plus Environment

Page 40 of 42

Page 41 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(67) Cao, Y.; Petzold, L. Accuracy limitations and the measurement of errors in the stochastic simulation of chemically reacting systems. J. Comput. Phys. 2006, 212, 6 – 24. (68) Allen, M.; Tildesley, D. Computer Simulation of Liquids; Oxford Science Publications; Clarendon Press, 1989. (69) Weisstein, E. W. Least squares fitting. From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/LeastSquaresFitting.html. (70) Fetters, L. J.; Graessley, W. W.; Krishnamoorti, R.; Lohse, D. J. Melt chain dimensions of poly(ethylene-1-butene) copolymers via small angle neutron scattering. Macromolecules 1997, 30, 4973 – 4977.

41

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 42 of 42