key: cord-1034860-wnaayk4j
authors: Chaimovich, Mark; Mathematics, Aviel Chaimovich Russian School of; Bethesda, North; Maryland,; University, Drexel; Philadelphia,; Pennsylvania,
title: Relative Resolution: A Computationally Efficient Implementation in LAMMPS
date: 2021-04-20
journal: Journal of chemical theory and computation
DOI: 10.1021/acs.jctc.0c01003
sha: 489278394835f7b45ae7076a19e7816959d3a561
doc_id: 1034860
cord_uid: wnaayk4j

Recently, a novel type of a multiscale simulation, called Relative Resolution (RelRes), was introduced. In a single system, molecules switch their resolution in terms of their relative separation, with near neighbors interacting via fine-grained potentials yet far neighbors interacting via coarse-grained potentials; notably, these two potentials are analytically parameterized by a multipole approximation. This multiscale approach is consequently able to correctly retrieve across state space, the structural and thermal, as well as static and dynamic, behavior of various nonpolar mixtures. Our current work focuses on the practical implementation of RelRes in LAMMPS, specifically for the commonly used Lennard-Jones potential. By examining various correlations and properties of several alkane liquids, including complex solutions of alternate cooligomers and block copolymers, we confirm the validity of this automated LAMMPS algorithm. Most importantly, we demonstrate that this RelRes implementation gains almost an order of magnitude in computational efficiency, as compared with conventional simulations. We thus recommend this novel LAMMPS algorithm for anyone studying systems governed by Lennard-Jones interactions.

Ever since the advent of molecular simulations, a main goal of the scientific community has been enhancing their computational efficiency, and this has been frequently done by modifying the interactions between the molecules. [1] [2] [3] [4] These days, multiscale algorithms have become especially popular for such a task: At the most basic level, these schemes link together detailed fine-grained (FG) molecular models and simplified coarse-grained (CG) molecular models. [5] [6] While there have been various recipes for combining the two models in a single system, [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] formulations which switch between the FG and CG interactions in terms of a spatial variable have been particularly common. [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] In this work, we specifically focus on the multiscale framework which we call Relative Resolution (RelRes). [23] [24] [25] [26] The main signature of RelRes is that the molecular model is a function of the relative separation: Given two arbitrary molecules, they are modeled via FG interactions if near to each other, yet via CG interactions if far from each other. Such an idea was introduced by Izvekov and Voth in molecular simulations of water, 23 whilst Shen and Hu employed a reminiscent approach in biomolecular systems. 24 Most importantly, while other multiscale approaches typically employ a numerical parametrization between the FG and CG models, [27] [28] Chaimovich at al. discovered a crucial benefit of RelRes: It naturally possesses an analytical relation between the FG and CG potentials, stemming in a multipole approximation. [25] [26] In fact, by correspondingly employing tabulated potentials in GROMACS, 29 it was shown that RelRes overcomes the transferability and representability issues that hinder most multiscale protocols [30] [31] : In particular for several nonpolar mixtures, RelRes correctly captures across state space numerous structural correlations and thermal properties of static and dynamic behavior. [25] [26] In principle, employing RelRes for molecular simulations can yield a significant gain in computational efficiency.

Foremost in preparing a multiscale system, the proposed parametrization from the FG model to the CG model can be done analytically on a paper sheet rather than numerically on a computer cluster; this already saves many days of work. Regarding the multiscale simulation itself, we also expect a computational speedup: Although RelRes maintains all degrees of freedom, it can significantly reduce the number of pairwise interactions.

Compared with conventional simulations (which have just FG and no CG interactions), RelRes has a similar number of near neighbors but a modest number of far neighbors. While this computational efficiency is expected for RelRes, it has not yet been unequivocally demonstrated.

In this work, we most notably show that if a molecular simulation is implemented via RelRes, it gains almost an order of magnitude in computational efficiency (as compared with a conventional simulation). We in fact show this computational enhancement by implementing

RelRes specifically for the Lennard-Jones (LJ) potential in the computational package of LAMMPS. 32 The LJ potential is frequently used by the LAMMPS community in various applications, [33] [34] [35] [36] [37] [38] especially in studying elementary alkanes, as well as oligomeric and polymeric hydrocarbons. [39] [40] [41] [42] [43] [44] The main benefit of LAMMPS for RelRes is that it inherently supports multiple neighbor lists (based on various cutting distances), and thus, the computational efficiency for RelRes is successfully achieved. Ultimately, we exemplify our LAMMPS algorithm on several alkane liquids (propane, neopentane, alternate cooligomer, block copolymer, etc.). Our current work thus permits the scientific community for the automated use of RelRes in LAMMPS, for the purpose of efficiently studying any nonpolar system.

As mentioned earlier, the main goal of our work is implementing RelRes in LAMMPS, 32 specifically for the LJ potential. Here, we foremost review the theoretical foundation of the RelRes framework, [25] [26] introducing several practical modifications along the way. We also thoroughly describe all the technical aspects of our implementation in LAMMPS.

Consider a molecular system. Naturally, each molecule is composed of several atoms; in this article, we will commonly use the term atom for the concept of a "united" atom (e.g., a principal carbon with its adjacent hydrogens). In the many molecular simulations, each atom is represented by a FG site. It may be convenient to break up each molecule into several groups of atoms. In such cases, each group is represented by a CG site. The distinct aspect of the various multiscale strategies is that the FG and CG models are combined in a single molecular simulation. [17] [18] [19] [20] [21] [22] Fig. 1 [23] [24] [25] [26] We will now define a distance parameter at which interactions switch between the FG potential and the CG potential, and we will specifically call it the "switching distance" ;

realize that in principle, may be distinct for each molecular pair , yet we omit these indices for clarity. In the most recent formulation of RelRes, [25] [26] its FG potential is defined as and its CG potential is defined as In essence, these are just Heaviside functions, with the constants and respectively providing continuity of the FG and CG potentials at . The total energy of the system can be expressed by combining Eqs. (1) and (2), while performing the appropriate summation over all sites

Notice that with , we obtain a pure FG system, and with , we obtain a pure CG system. We will call the former the "reference" system, as this is frequently the benchmark used in molecular simulations (i.e., it supposedly describes the correct behavior of a liquid of interest).

In order to ensure a correct representation of the system behavior, multiscale strategies require appropriate parametrization between the FG and CG models. To perform such a parametrization, different numerical route can be employed. [27] [28] [45] [46] In such established procedures, the parametrization is notably performed across the entire domain of the pairwise distance. While the initial efforts in RelRes essentially used such strategies, 23-24 the most recent RelRes framework found a natural way of connecting the FG and CG potentials. [25] [26] Above all,

RelRes only cares about replacing the FG potential with the CG potential beyond . Thus, RelRes proceeds by equating the FG and CG energies at the "infinite limit" of the relative separation between an isolated molecular pair (i.e., ); this yields the following analytical parametrization 25 With this parametrization, it is obvious that the RelRes potential defined by Eq. (3) will approximately capture the behavior of the reference system if for all sites (i.e., the dimension of each group of atoms is relatively small, in comparison with the switching distance); this is in fact why we often partition a single molecule into several groups. Note that the formal derivation of this parametrization proceeds via a Taylor series. 26 Here, we present just the zero-order term of this multipole approximation, which is sufficient for nonpolar molecules; for polar scenarios, the zero-order term may vanish, and thus, the first-order or second-order terms become essential.

In Refs. 25-26, the CG site was placed at the center of mass given by the FG sites, which required the construction of a virtual site; such mapping is depicted in the top portion of Fig.1 .

However, we are not interested in exact trajectories in molecular simulations, and thus, the placement at the center of mass is not an absolute requirement. In this current work, we place the CG site on one of its FG sites, which logically should be the site closest to the center of mass;

such mapping is depicted in the bottom portion of Fig. 1 

Naturally, molecular simulations deal with interactions of a finite range, and thus, we must introduce a cutting distance in the RelRes potential of Eq. (6) (i.e., for ). Note that similar with , can be in principle different for different groups, and of course, .

Importantly, most packages, including LAMMPS, are based on Newtonian motion, and therefore, from a practical perspective, Eq. (6) must be modified, so that there is no singularity in the force (i.e., ) at , as well as at . For this purpose, we define : the FG potential is only present for , and the CG potential is only present for . In the interval (i.e, the smoothing zone), the potential is governed by the smoothing function.

In an analogous manner, a smoothing function is also applied at the cutting distance, with its inner and outer values defined by . Note that we often omit the indices and (e.g., ), if our discussion is relevant for both smoothing zones.

In principle, the smoothing zone can have any functionality; it is just constrained by the appropriate application of the boundary conditions at, , which at the most basic level, just ensure the continuity of the energy , together with the force, . For enhanced smoothness, we also apply an additional continuity requirement for the ensuing derivative,

. We will specifically cast each smoothing function in terms of the relative location within each smoothing zone, , which in turn means that .

In practice, we have three smoothing functions. Reminiscent of the familiar situation, we have two smoothing functions which shut down interactions: turns off the FG force at , smoothing it to zero at , and turns off the CG force at , smoothing it to zero at . Notably, we also have one smoothing function which starts up an interaction: turns on the CG force at , smoothing it from zero at . Note that in our notation, a fourth smoothing function in principle exists: does nothing for the FG force at (it is formally zero across its entire domain).

For all these functions, we have five boundary conditions. Four of the conditions take care of the first and second derivatives of the potential at both boundaries, and . The fifth condition takes care of the continuity of the potential at one of the boundaries. Note that no sixth condition is required since the shifting constants of Eq. (6) take care of the continuity of the potential at the other boundary.

The boundary conditions for and are very similar, while the boundary conditions for are slightly different. In the following definitions of the boundary conditions, we omit their FG and CG superscripts for clarity (as is often done throughout the manuscript, we also omit the and of the subscripts). For and , we have:

and for , we have:

In turn, the shifting constants guaranteeing the continuity of the potential on the other boundary of the smoothing zone will be Note that while in principle exists, it is formally identical with . Now we are ready to formulate our entire RelRes potential for practical use in molecular simulations:

The basis for most molecular simulations is the LJ potential. In general, it is perfectly sufficient for nonpolar molecules, and in conjunction with the Coulomb potential, it can adequately describe polar molecules as well. Indeed, the LAMMPS package supports multiple variations of the LJ potential. Therefore, we consider here the utilization of the LJ potential with

RelRes.

For the interaction between FG sites and , the LJ potential can be expressed as with and being its respective length and energy parameters. Substituting this expression in Eq. (4), we obtain that the interaction between CG sites and is governed by a similar expression where Eq. (13) gives us a simple analytical relation between FG and CG parameters.

Now we need to define the CG interactions for all site combinations in a similar manner as done in Eq. (5) . Specifically for the LJ potential it is fulfilled by setting

Having these relations in mind, while respectively replacing and in Eq. (10) with Eqs. (11) and (12) we attain the expression for RelRes with the LJ potential.

The parameters between the same FG sites, and , are often known. In turn, the mixed parameters between arbitrary sites can be determined via a geometric mean: and . We can then simplify Eq. (13) for CG site comprised from FG sites to obtain and subsequently, the geometric mean holds here as well, and

. To summarize, for the LJ potential Eqs. (13) and (15) are very convenient: They determine the CG parameters analytically, with no need for an actual molecular simulation; only knowledge of the FG parameters of the reference system is required. Eq. (13) can be always employed, but if geometric mixing is applicable in the LJ potential, Eq. (15) is the relevant one for use in molecular simulations.

Let us now discuss our smoothing functions. Recall that RelRes employs three such functions (i.e., , , and ). In our arguments, we may omit superscripts for clarity. In the spirit of LAMMPS, we choose a polynomial for the smoothing function where the five coefficients can ensure the continuity of the potential with its first and second derivatives, and again, is the relative location within the smoothing zone. Applying the boundary conditions of Eqs. (7) and (8), we obtain a system of liner equations for each smoothing function, and in turn we get the following and Besides, with this notation, the shifting constants of Eq. (9) can be expressed as

The introduction of RelRes into the LAMMPS package requires computing pairwise interactions according to the RelRes potential of Eq. (10). To specifically use this equation with the LJ potential, we need to substitute Eqs. (11) and (12) in place of the general of Eq. (10).

Besides, we substitute the polynomial smoothing from Eq. (16) in place of the general . In turn, we obtain the following piecewise potential:

where the coefficients of the smoothing functions are defined by Eqs. (17) and (18) The blue curves represent our ideal model for propane (i.e., a mapping ratio of , with a switching distance of ), and the red curves represent our ideal model for neopentane (i.e., a mapping ratio of , with a switching distance of ). Each dashed vertical line corresponds with the location of the switching distance, while its shaded region represents the smoothing zone. The mapping for these alkanes is described in Fig. 3 , and their LJ parameters are specifically given in Table 1 .

LJ potential). Note here that the interaction parameters must be defined for all atom types, and thus, it generally requires multiple pair_coeff commands. Besides, the pair_coeff command provides an option to override the global values defined by the pair_style command (in the case of lj/smooth, these arguments are the specific inner and outer cutting distances for the specified pair of atom types).

Analogous with Eq. (20) , lj/relres requires eight parameters: Besides the four LJ energy and length parameters of the FG and CG potentials (i.e., ), there are the four distance parameters of the smoothing zones (i.e., ). In our implementation, the distance parameters are defined globally in the pair_style command; there is an option of overriding them for a particular atom type with the pair_coeff command. On the other hand, since the LJ energy and length parameters are specific for each atom type, they are defined in the pair_coeff command. Here is the formal syntax of these commands in the input script:

where and are atom types, and brackets denote optional arguments. The default value is geometric, yet if desired, the other typical mixing rules of LAMMPS are also available.

As an illustration, let us define a hypothetical system which has a group with two atom types: The first type corresponds to hybrid sites, and the second type corresponds to ordinary sites. Here is the input script relevant for such a supposed system:

pair_style lj/relres 0.68e-9 0.72e-9 1.20e-9 1.40e-9 #Line1 pair_coeff 1 1 0.80e-21 0.37e-9 14.9e-21 0.39e-9 #Line2 pair_coeff 2 2 1.00e-21 0.40e-9 0.00 0.39e-9 #Line3 pair_modify shift yes #Line4

Everything is presented here in SI units, with and . For convenience, we numbered the command lines; as usual, a text following a pound sign is processed as a comment. pair_coeff 1 1 0.80e-21 0.37e-9 14.9e-21 0.39e-9 & 0.58e-9 0.62e-9 1.20e-9 1.40e-9 #Line2 pair_coeff 2 2 1.00e-21 0.40e-9 0.00 0.39e-9 & 0.58e-9 0.62e-9 1.20e-9 1.40e-9 #Line3

Note that the ampersand sign indicates the continuation of a command onto the next line.

In addition, the flexibility of LAMMPS allows employing the lj/relres style as one of the subtypes within the hybrid/overlay (also hybrid) pair style, in turn, enabling its use with multiple other potentials in a single molecular simulation. Most notably, it can be combined with the Coulomb potential in studying polar systems. Here is an illustration of a possible combination:

pair_style hybrid/overlay lj/relres r s-r s+ r c-r c+ coul/cut r c pair_coeff α β lj/relres ϵ α

Here, a name of a pair style appears as an argument in the pair_coeff command, allowing for the definition of multiple potentials for the same pair of atom types. Note that these other potentials maintain their FG nature, and thus, no CG treatment is required for these interactions.

Calculation of pairwise interactions is the main consumption factor in molecular simulations: Usually, these contribute to more than 90% of the CPU time. Due to a reduced number of pairwise interactions, CG models have a computational advantage over FG models, and, therefore, we should expect a similar effect in a multiscale system like RelRes.

LAMMPS, like most other computational packages, uses a neighbor list for calculating pairwise energetics. In general, only sites within a specific distance (slightly beyond the cutting distance) are included in the neighbor list. Importantly, the interactions are computed merely with these relevant neighbors, and therefore, the CPU time is proportional to the characteristic size of the neighbor list, .

In the RelRes algorithm, the list size depends on the category of the site (an ordinary site or a hybrid site). In either case, the size of the list is significantly reduced as compared with a reference system. Foremost, an ordinary site has less neighbors as its list size is governed by the fairly short rather than the fairly long . On the other hand, in consideration of the mapping ratio , a hybrid site has less neighbors because beyond the switching distance, it interacts only with other hybrid sites yet not with any ordinary sites. Therefore, the decrease of and the increase of are the main factors that can enhance the computational efficiency of RelRes. where is an empirical parameter that depends on the ratio ; when this ratio approaches unity, approaches unity as well, yet is always less than one.

We now shift our attention to testing the RelRes algorithm for the LJ potential in LAMMPS.

Can it correctly describe the structural and thermal behavior of various liquids, and can it do so with an enhanced computational efficiency? The ultimate goal is to show that our LAMMPS implementation can be effectively used for complex mixtures (solutions with cooligomers, copolymers, etc.), and for achieving this aim, we systematically evolve here simple systems of "monomeric" and "dimeric" liquids. We succeed in this piecemeal task by employing the appropriate switching distance (between the FG and CG potentials), together with the appropriate parametrization of Eq. (15) (for the and parameters). We specifically examine the following alkanes:

 Two distinct monomers, C 3 H 8 (propane) and C 5 H 12 (neopentane), which form the basis for the remainder of our investigation.

 Two basic dimers, C 3 H 7 -C 3 H 7 (hexane) and C 5 H 11 -C 5 H 11 (bineopentyl), both of which are symmetric extensions of our monomers.

 A dimer C 3 H 7 -C 5 H 11 (2,2-dimethylhexane) which is an asymmetric construction of our monomers.

 A dimer C 4 H 9 -C 4 H 9 (biisobutyl) which is a symmetric construction from groups, with a mapping ratio of , that have not been investigated as monomers.

 An isobutane solution which contains four cooligomeric chains C 3 H 7 -(C 5 H 10 -C 3 H 6 ) 4 -C 5 H 11 .

 An isobutane solution which contains two copolymeric chains C 3 H 7 -(C 3 H 6 ) 9 -(C 5 H 10 ) 9 -C 5 H 11 .

Before utilizing RelRes, the mapping of atoms into groups must be defined. The previous studies suggested that to keep manageable (i.e., around the pairwise distance associated with nearest neighbors), the CG site should be roughly within one bond length from the FG sites. [25] [26] Taking this into consideration, let us consider an arbitrary hydrocarbon molecule. Each "united" atom can be bonded to one, two, three, or four other "united" atoms. A couple of these scenarios are depicted in the top portion of Fig. 3 ; the "two-case" is shown on the left, and the "four-case"

is shown on the right. In our protocol for this article, the central "united" atom will always be the hybrid site, while the peripheral "united" atoms will always be the ordinary sites. The mapping ratio of each group is equal to the total number of "united" atoms it contains.

Of course, such groups, if bonded between each other, can be the building blocks for various oligomeric and polymeric chains. Hence, we have the nomenclature for our work: The orange circles delineate an effective boundary of a group; light orange is for , and dark orange is for (the former is always on the left, and the latter is always on the right). Actual atoms are given in green: Hybrid sites are represented by dark green circles, while ordinary sites are represented by light green circles. The top portion focuses on two examples of our building blocks for any alkane chain. Subscripts "a", "b", "c" and "d" represent the number of hydrogens in a united atom. The bottom portion, with its three sections, presents the various molecules in our study. Each monomer consists just of one building block (C 3 H 8 or C 5 H 12 ). Each dimer consists of two building blocks: Some of the dimers are just symmetric combinations of our two monomers (C 3 H 7 -C 3 H 7 or C 5 H 11 -C 5 H 11 ); another dimer is an asymmetric combination of our two monomers (C 3 H 7 -C 5 H 11 ) having a characteristic mapping ratio of , while another dimer is a symmetric combination of a different monomer (C 4 H 9 -C 4 H 9 ) with . There are also two chain molecules with a characteristic mapping ratio of : an alternate cooligomer (C 3 H 7 -(C 5 H 10 -C 3 H 6 ) 4 -C 5 H 11 ) and a block copolymer (C 3 H 7 -(C 3 H 6 ) 9 -(C 5 H 10 ) 9 Block copolymer … … C 3 H 7 -(C 3 H 6 ) 9 -(C 5 H 10 ) 9 -C 5 H 11

From the perspective of LAMMPS, we must make an important clarification regarding the definitions of atom types in the context of such mapping. As an illustration, we present here an excerpt of the data file of LAMMPS for the dimer C 3 H 7 -C 5 H 11 (shown in the middle section of the bottom portion of Fig. 3 We notably focus on the Atoms section of the data file, showing just two molecules. Each line corresponds to a particular site, with the position omitted for clarity. The first column specifies the identifier of the site, while the second column specifies the identifier of its molecule. Most importantly, column three defines the atom type of the site, and the comment following the pound sign uniquely characterizes this atom type: It is given by its category (ordinary or hybrid),

its atom (C, CH 2 , or CH 3 ), and its group (C 3 H 7 or C 5 H 11 ). Following all possible combinations for the dimer, there are six atom types in total, specifically three atom types for each group (two correspond with ordinary sites, and one corresponds with a hybrid site). Do realize that here, CH 2 as an ordinary site and CH 2 as a hybrid site have same FG parameters yet different CG parameters, and therefore, they must be defined as distinct atom types.

On another note, this asymmetric dimer brings further subtleties for RelRes: The mapping ratio is different for its two groups. In this case, it is convenient to introduce the concept of a characteristic mapping ratio for the entire system. We define it as the ratio between the number of all sites and the number of hybrid sites. With this definition, observe that for C 3 H 7 -C 5 H 11 the characteristic mapping ratio is (the building blocks themselves have and ). Note that the oligomeric-polymeric solutions also have 4 as the characteristic mapping ratio.

Without loss of generality, for the non-bonded interactions of the reference models, we use here the set of LJ parameters of OPLS_UA. 47 These parameters are presented in To compare the thermal behavior of the fluids, we calculated intra-molecular and intermolecular potential energy. In addition, we calculated the virial associated with intramolecular and inter-molecular forces. For the oligomeric-polymeric solutions, the energy was calculated separately for the solvent and solute molecules. All relevant data was outputted every timesteps. This data was used for calculating static properties such as arithmetic means (e.g., , , etc.) and standard deviations (e.g., , , etc.). This data was also used for evaluating dynamic properties, such as auto-correlations in terms of time (e.g., , , etc.). Note that we frequently normalize thermal properties, denoted by the symbol "*": This corresponds with a ratio in a property between the RelRes system and the reference system.

We begin our investigation with the two monomeric liquids, each one having its own mapping ratio (i.e., C 3 H 8 with and C 5 H 12 with ; see mapping in the upper section of the bottom portion of Fig. 3 ). Our main goal here is in determining the ideal switching distance for a given mapping ratio. Besides, we also validate here the parametrization of Eq. (15) .

Consequently, the monomers become the basis for all other systems we test.

As discussed earlier, a shorter corresponds with a faster computation, and thus, the ideal is given as the shortest switching distance that allows for the correct retrieval of structural and thermal behavior. In the previous work, a conjecture was made that the proper switching distance is between the first and second coordination shells of a certain liquid. 25 It was in turn shown that this specifically corresponds with the midpoint between the maximum and minimum of (between the centers of mass), which is roughly in the vicinity of the respective inflection point. 26 For the reference liquids, Fig. 4 We also show in Fig. 4 between two ordinary sites, and the middle panels show between a hybrid site and an ordinary site. Table 2 . Signature distances in the reference of all studied liquids. These are specifically the midpoint between the two extrema, as well as its adjacent inflection point. The rightmost column also shows our ultimate suggestions for the switching distances in RelRes. by about . This is why, from a structural perspective, our recommendation for the proper switching distance for the mapping ratio of is (its is off by less than 1%). Table 3 also presents the diffusion coefficient. RelRes adequately represents this facet of the reference system, having less than 10% discrepancy for it, which is within its statistical error.

This good representation stems mostly in the fact that RelRes keeps all degrees of freedom.

Moreover, as the diffusivity is well correlated with many other transport coefficients (e.g., viscosity), we expect that RelRes will yield an adequate representation for these as well.

Besides these structural correlations, we find that RelRes also captures thermal properties satisfactorily well. Table 4 shows the potential and virial associated with the inter-molecular energetics, specifically presenting the values of their corresponding arithmetic means and standard deviations. In all cases, we naturally find the overall trend that the larger yields a smaller discrepancy in both the potential and virial energies. For C 3 H 8 , we find sufficiently good results with : The discrepancy in the mean values is always less than 1%, while the discrepancy in the deviation values is always less than 7%. For C 5 H 12 , we find sufficiently good results only with : The discrepancy in the mean values is always less than 3%, while the discrepancy in the deviation values is always less than 3%. These discrepancies are more or Table 3 . Comparison of structural characteristics for the studied liquids. The coordinates of the minima and maxima for the radial distributions, as well as the diffusion coefficients, are presented for the monomers, dimers, and the solvent of the cooligomer and copolymer systems. For the cooligomeric and copolymeric solutes themselves, the radius of gyration is presented instead. For the systems comprised from different building blocks, characteristic values (marked with *) are given for their mapping ratios, as well as for their switching distances. For C 3 H 7 -C 5 H 11 , results for between different hybrid sites (i.e., CH 2 :C of the middle panel of Fig. 6 ) are provided; the intra-molecular extrema of Fig. 6 less within the range of the statistical errors. Note that we also look at the intra-molecular energetics of the potential and virial (not presented in Table 4 ): In all cases, we observe no noticeable discrepancies, and this is expected, since our multiscale strategy does not alter this aspect of the reference liquid. This concludes our recommendation for the proper switching distance: for the mapping ratio of , and for the mapping ratio of . Overall, this again nicely corresponds with the signature distances we emphasized in Table 2 .

For these monomers, we conducted an additional study to verify that the parametrization from the FG potential to the CG potential, as defined by Eq. (15) , is in fact optimal. For this purpose, we introduce two tuning factors, and , and we vary them while running the RelRes simulations with CG parameters and (note that and yield the ideal parameters of Eq. (15)). In Fig. 5 , we plot several thermal properties as functions of (bottom abscissa, squares, darker solid lines) and (top abscissa, diamonds, lighter dashed lines): Specifically, we plot the inter-molecular potentials in the top panels, and the intramolecular potentials in the bottom panels; the arithmetic mean is in red, and the standard deviation is in green. Again, the left panels are for C 3 H 8 , and the right panels are for C 5 H 12 .

Error bars are not shown: For the energy mean, they are within the size of the markers, and for for such optimization may be helpful, [45] [46] but it will obviously require significant computational resources as compared with the analytical relation of Eq. (15) , and thus, it is beyond the scope of our current work.

Proceeding with RelRes for dimers, we adopt the recommendations from the analysis of the monomers: for (C 3 H 7 ) and for (C 5 H 11 ). In turn, we run simulations with the same protocol as for the monomers. Our main aim here is in determining if the switching distances are mostly governed by the mapping ratio (rather than by any peculiarities in the molecular complexities).

Let us discuss briefly the results for the two symmetric dimers (i.e., C 3 H 7 -C 3 H 7 and C 5 H 11 -C 5 H 11 , see the mapping for these dimers in the middle section of the bottom portion of Fig. 3) .

FIG. 5. The components of the potential energy, normalized by their respective values from the reference system. The top panels are for the inter-molecular energy, and the bottom panels are for the intra-molecular energy; the left panels are for C 3 H 8 , and the right panels are for C 5 H 12 . The data is plotted in terms of (squares connected by solid lines, plotted on the bottom abscissa) and (diamonds connected by dashed lines, plotted on the top abscissa). The red color is for the arithmetic mean of the energies (plotted on the left ordinate), and the green color is for the standard deviation of the energies (plotted on the right ordinate). Tables 3 and 4 :

The data confirms that the recommended switching distances provide a good representation of the arithmetic mean and standard deviation of the energy, the virial, as well as the diffusion coefficient. Fig. 7 shows a comparison of a dynamic correlation of two thermal functions, intramolecular and inter-molecular components of the virial. It also shows a good match between the RelRes and reference models. Graphs are not shown as they are very similar to all other liquids studied.

The study of the dimers confirms the robustness of our multiscale strategy. Foremost, the selection of the switching distance is mainly dependent on the mapping ratio of the group.

Interestingly, with an increase of size of molecules (i.e., from monomers to dimers, etc.), the performance of RelRes even slightly improves with the same switching distance. This gives a good hope that this improvement will be even better for more complex systems.

In the last stage of our verification, we dilute a few cooligomers or copolymers in a solvent, building in turn two complex solutions. For the solvent, we use isobutane (C 4 H 10 ) in both systems. In one system, we immerse an alternate cooligomer C 3 H 7 -(C 5 H 10 -C 3 H 6 ) 4 -C 5 H 11 . This cooligomer has 10 groups: one C 3 H 7 group, four C 3 H 6 groups, four C 5 H 10 groups and one C 5 H 11 group. In the second system, we immerse a block copolymer C 3 H 7 -(C 3 H 6 ) 9 -(C 5 H 10 ) 9 -C 5 H 11 .

This copolymer has 20 groups: one C 3 H 7 group, nine C 3 H 6 groups, nine C 5 H 10 groups and one C 5 H 11 group. The mapping of these two systems is shown in the lower section of the bottom portion of Fig. 3 .

For the complex solutions, the switching distances were based on the mapping ratio of each group: Relevant for the solutes, for the C 3 H 7 and C 3 H 6 groups ( ), and for the C 5 H 11 and C 5 H 10 groups ( ); regarding the solvent, for the C 4 H 10 group ( ). Similar to the dimer C 3 H 7 -C 5 H 11 , the value of represents a characteristic switching distance of the systems.

First, let us discuss the behavior of the solute. In both systems, we are getting a good representation of the properties. From a structural perspective (see Table 3 ), the radius of gyration in both cases is within 1% from the reference system, which is in the range of the statistical error (about 2%). From a thermal perspective, intra-molecular energy is mainly not impacted by RelRes. For inter-molecular energy, the results are also good: Both the mean and deviation in the energy is within the statistical error of roughly 10%. To make these results evident, we produced graphs of probability distribution of components of the potential energy, specifically for the copolymer molecules. They are presented in Fig. 8 . The top panel is for the inter-molecular energy (upper abscissa), and the bottom panel is for the intra-molecular component (lower abscissa). The reference system is shown as the black solid line, and the

RelRes system is shown as the violet dashed line.

Let us now discuss the solvent. With , in both solutions, RelRes perfectly describes structural (see Table 3 ) and thermal (not shown) behavior of the reference model, and the properties of the solvent are very similar in both cases. Besides, curves of dynamic correlation of virial components for the copolymer solution are plotted in Fig. 9 , which is analogous in format to Fig. 7 . It also shows a good match between the models. For cooligomers, graphs are very similar, yet they are not shown here.

At this point, we have shown that RelRes in LAMMPS can capture the behavior of nonpolar liquids very well. However, is it actually worth using in consideration of computational efficiency? Here, we in fact show that this is the case. Focus for now just on the region encompassed by the dashed ellipse. This is the region of all systems with the recommended switching distances, as given by Table 2 . Most importantly, the molecular simulations using RelRes are almost fivefold faster than the corresponding reference systems in calculating the pairwise interactions, and in turn, they are about fourfold faster in terms of the total CPU time. Besides, notice that in this region, all symbols of the same 21) and (22) using .

The dashed ellipse indicates the region in the graph which encompasses the systems with the recommended switching distances.

Complex system color essentially collapse on each other: This means that with the recommended choice for the switching distance, systems with the same mapping ratio have an almost identical computational efficiency, rendering the computational enhancement of RelRes independent of the molecular complexity. Therefore, we expect a similar speedup of almost an order magnitude for other nonpolar mixtures using RelRes.

We now shift our focus on the triangles of Fig. 10 , which represent the two monomers, C 3 H 8 and C 5 H 12 . In the elliptical region, the blue triangle is for with , and the red triangle is for with . Going beyond the elliptical region, the lowest blue triangle is for with , and the highest red triangle is for with (realize that these switching distances are not the recommended ones). Regarding the switching distance, given the same color of triangles, the higher consistently yields the lower computational efficiency. Regarding the mapping ratio, the higher consistently renders the higher computational efficiency. Altogether, this confirms our earlier intuition that decreasing and increasing are the main factors for improving the computational efficiency of RelRes.

Finally, focus on the black dashed lines in the top and bottom panels of Fig. 10 . These respectively confirm our predictions formulated in Eqs. (21) and (22): The computational efficiency is directly proportional to the reduction in the characteristic size of the neighbor list.

Besides, we also perform a linear regression through the origin (not shown in Fig. 10 ): For , we find a slope of 1.0065 with a determination coefficient of 0.94, and for , we find a slope of 0.8248 with a determination coefficient of 0.96. These trends affirm that RelRes can reduce the CPU time of reference systems by almost an order of magnitude.

Here, we have continued the work on RelRes, the multiscale algorithm which describes near neighbors with FG models and far neighbors with CG models. [23] [24] [25] [26] Foremost, we automated RelRes in the computational package of LAMMPS, particularly for the LJ potential. Indeed, we developed a new C++ class, pair_lj_relres, with its pair style, lj/relres, that calculates pairwise interactions according to Eq. (20) . Most importantly in the current work, we demonstrated that RelRes is in fact computationally efficient. In particular, our implementation in LAMMPS reduces the computational time of molecular simulations by a factor of about -.

We have essentially tested this implementation on several alkane liquids, being as complex as solutions of alternate cooligomers and block copolymers (i.e., C 3 H 7 -(C 5 H 10 -C 3 H 6 ) 4 -C 5 H 11 and C 3 H 7 -(C 3 H 6 ) 9 -(C 5 H 10 ) 9 -C 5 H 11 , respectively). Note that just as in earlier studies of RelRes that parameterize between the FG and CG potentials via a multipole approximation, 25 Interestingly, we have found that the proper switching distance for RelRes is generally driven by the mapping ratio. For example, for both propane and hexane (i.e., C 3 H 8 and C 3 H 7 -C 3 H 7 , respectively), with a mapping ratio of , , yet for both neopentane and bineopentyl (i.e., C 5 H 12 and C 5 H 11 -C 5 H 11 , respectively), with a mapping ratio of , .

If one is interested in an even better representation of a reference system, a larger switching distance can be employed, while giving up some computational efficiency (e.g., moving from to costs about 20% of CPU time). Note that the current work restricts the mapping for groups whose ordinary sites are directly bonded to hybrid sites. Removing this limitation can favorably increase the mapping ratio but unfavorably increase the switching distance; exploring the computational efficiency of such alternatives becomes especially of interest.

Although the current implementation only deals with the LJ potential, it can be also used for systems with other potentials as well, just by invoking the hybrid pair style in LAMMPS.

Moreover, there is the possibility to develop RelRes variations of the Coulomb potential (e.g., pair styles coul/relres or lj/coul/relres in the LAMMPS standard). This kind of an algorithm would be very useful for a group of atoms with a net charge (e.g., nitrate, ammonium, etc.); this would just require replacing and in Eq. (4) with the appropriate Coulomb expressions, in a similar manner as was done with Eqs. (11) and (12) for the LJ potential.

However, for a group of atoms with no net charge (e.g., water), a more elaborate approach is necessary, since the zero order term of a multipole approximation (i.e., Eq. (4)) is not sufficient:

Chaimovich et al. establishes the appropriate framework by deriving the first and second terms of the relevant Taylor series. [25] [26] Implementing such functionalities would require employing the dipole in LAMMPS (e.g., pair styles dipole/relres or lj/dipole/relres in the LAMMPS standard). All these routes are very practical directions for future research, since they could all substantially speedup molecular simulations of polar mixtures.

In summary, we recommend our automated version of RelRes in LAMMPS for anyone studying molecular systems that involve LJ interactions. As compared with conventional simulations, this RelRes implementation can yield almost an order of magnitude enhancement in computational efficiency, whilst capturing the behavior of nonpolar mixtures with negligible error.

The equation of state of the classical hard sphere fluid

Intermolecular forces and the nature of the liquid state: Liquids reflect in their bulk properties the attractions and repulsions of their constituent molecules

Equilibrium Structure of Simple Liquids

Role of the first coordination shell in determining the equilibrium structure and dynamics of simple liquids

From quantum to subcellular scales: multi-scale simulation approaches and the SIRAH force field

Molecular systems with open boundaries: Theory and simulation

Coarse-Grained Model of Proteins Incorporating Atomistic Detail of the Active Site

Mixed Atomistic and Coarse-Grained Molecular Dynamics: Simulation of a Membrane-Bound Ion Channel

Hybrid simulations: Combining atomistic and coarse-grained force fields using virtual sites

Further optimization of a hybrid united-atom and coarse-grained force field for folding simulations: Improved backbone hydration and interactions between charged side chains

Hierarchical modeling of polystyrene: From atomistic to coarse-grained simulations

Coarse-grained and reverse-mapped unitedatom simulations of long-chain atactic polystyrene melts: Structure, thermodynamic properties, chain conformation, and entanglements

Resolution Exchange Simulation

Multigraining: An algorithm for simultaneous fine-grained and coarse-grained simulation of molecular systems

Direct Mixing of Atomistic Solutes and Coarse-Grained Water

A different approach to dual-scale models

Adaptive resolution molecular-dynamics simulation: Changing the degrees of freedom on the fly

Concurrent dual-resolution monte carlo simulation of liquid methane

Energy conservation in adaptive hybrid atomistic/coarse-grain molecular dynamics

Hamiltonian Adaptive Resolution Simulation for Molecular Liquids

Adaptive resolution simulation of an atomistic protein in MARTINI water

Adaptive resolution simulation of a biomolecule and its hydration shell: Structural and dynamical properties

Mixed resolution modeling of interactions in condensed-phase systems

Resolution-adapted all-atomic and coarsegrained model for biomolecular simulations

Relative resolution: A hybrid formalism for fluid mixtures

Relative resolution: A multipole approximation at appropriate distances

Perspective: Coarse-grained models for biomolecular systems

Perspective on the Martini model

Algorithms for highly efficient, load-balanced, and scalable molecular simulation

Beware of density dependent pair potentials

On the representability problem and the physical meaning of coarse-grained models

Fast Parallel Algorithms for Short-Range Molecular Dynamics

Computing the crystal growth rate by the interface pinning method

Graphene Interactions: Friction, Superlubricity, and Exfoliation

Finite-Size Effects of Binary Mutual Diffusion Coefficients from Molecular Dynamics

Interfacial Free Energy, and Crystal Nucleation in a Supercooled Lennard-Jones Liquid

Hamiltonian Transformation to Compute Thermo-osmotic Forces

Grand Canonical Inverse Design of Multicomponent Colloidal Crystals

Independence between friction and velocity distribution in fluids subjected to severe shearing and confinement

Volume-based mixing rules for viscosities of methane + n-butane liquid mixtures

Evaluation and Refinement of the General AMBER Force Field for Nineteen Pure Organic Electrolyte Solvents

Heterogeneous nucleation of an n-alkane on graphene-like materials

Stress Relaxation in Highly Oriented Melts of Entangled Polymers

Vapor-Liquid Equilibrium Simulations of Hydrocarbons Using Molecular Dynamics with Long-Range Lennard-Jones Interactions

Relative entropy as a universal metric for multiscale errors

Coarse-graining errors and numerical optimization using a relative entropy framework

Optimized intermolecular potential functions for liquid hydrocarbons

A new force field for molecular mechanical simulation of nucleic acids and proteins

We are grateful for all those who gave us moral and spiritual support in this project during the Coronavirus Pandemic.