Scaling Properties of Liquid Dynamics Predicted from a Single Configuration: Pseudoisomorphs for Harmonic-Bonded Molecules

Isomorphs are curves of constant excess entropy in the thermodynamic phase diagram along which structure and dynamics are invariant to a good approximation Gnan et al. (2009). We remind that the excess entropy ${S}_{\rm ex}$ is the entropy minus that of an ideal gas at the same density and temperature Hansen and McDonald (2013). Systems with isomorphs, termed R-simple Malins et al. (2013); Flenner et al. (2014); Prasad and Chakravarty (2014); Schrøder and Dyre (2014); Heyes et al. (2015); Khrapak et al. (2016); Kaskosz et al. (2023), are characterized by a strong correlation between the canonical-ensemble equilibrium fluctuations of potential energy $U$ and virial $W$ as quantified by the Pearson correlation coefficient (sharp brackets denote canonical averages and $\Delta$ the deviation from the mean),

The criterion for being R-simple is $R>0.9$ Gnan et al. (2009), but even systems with somewhat lower R values may have good isomorphs Ingebrigtsen et al. (2012).

Isomorph invariance of structure and dynamics refers to the use of units where the energy unit is $e_{0}\equiv k_{B}T$, the length unit is $l_{0}\equiv\rho^{-1/3}$ in which $\rho$ is the particle number density, i.e., the density of atoms, and the time unit is $t_{0}\equiv\rho^{-1/3}\sqrt{m/k_{B}T}$ in which $m$ is a particle mass. Quantities made dimensionless by reference to this unit system are termed “reduced” and marked by a tilde; for instance the reduced particle position $\mathbb{r}$ is given by ${\tilde{\mathbb{r}}}\equiv\mathbb{r}/l_{0}=\rho^{1/3}\mathbb{r}$.

The existence of isomorphs for a given system means that the thermodynamic phase diagram is essentially one-dimensional in regard to structure and dynamics. Thus if one imagines filming how the molecules move, two state points on the same isomorph would give rise to (almost) the same movie – except for an overall scaling of space and time. This provides an significant simplification of the physics, basically saying that the change of dynamics induced by a pressure increase may be counteracted by increasing at the same time the temperature. Many systems, however, do not have isomorphs in the original sense of the word as lines of constant excess entropy Gnan et al. (2009), which raises the question: Can some systems still have lines of (approximately) invariant structure and dynamics? This question motives the below presented investigation of two molecular models.

Because ${S}_{\rm ex}$ is the part of entropy that refers to particle positions, isomorphs are configurational adiabats. These can be traced out in the thermodynamic phase diagram by using the generally valid statistical-mechanical relation Gnan et al. (2009):

Evaluating the right-hand side by equilibrium $NVT$ simulations, an isomorph is traced out by solving the differential equation Eq. (2) numerically using, e.g., the Euler or Runge-Kutta methods Attia et al. (2021).

Isomorphs have been identified and validated for both atomic Gnan et al. (2009); Bøhling et al. (2012); Bacher et al. (2018); Heyes et al. (2019); Attia et al. (2021); Castello et al. (2021) and molecular systems Ingebrigtsen et al. (2012); Veldhorst et al. (2014), and isomorph-theory predictions have also been verified in experiments on glass-forming van der Waals molecular liquids Xiao et al. (2015); Hansen et al. (2018). In molecular systems, isomorphs are found when bonds are modeled as constraints Ingebrigtsen et al. (2012); Veldhorst et al. (2014), but not when the bonds are flexible Veldhorst et al. (2015); Olsen et al. (2016); Kaskosz et al. (2023). As an example, Fig. 1 shows scatter plots of virial versus potential energy for the asymmetric dumbbell (ASD) and 10-bead Lennard-Jones chain (LJC) models, both with harmonic springs as commonly used in simulations Koperwas and Paluch (2022); Kaskosz et al. (2023). Neither model is R-simple; the virial potential-energy correlation coefficients are $0.579$ and $0.284$, respectively. As expected, these models do not have isomorphs (data not shown), i.e., structure and dynamics are not invariant along the curves of constant ${S}_{\rm ex}$. Nevertheless, it has been found that both models have curves in the phase diagram of invariant structure and dynamics; termed “pseudoisomorphs” Veldhorst et al. (2014); Kaskosz et al. (2023).

Since pseudoisomorphs are not configurational adiabats, Eq. (2) cannot be used to identify such curves in the phase diagram. In 2016 Olsen et al. presented a method for tracing out pseudoisomorphs involving the following steps Olsen et al. (2016): 1) An equilibrium configuration is quenched to the nearest local minimum in the high-dimensional potential energy landscape, the so-called inherent state Stillinger and Weber (1983); 2) the Hessian matrix is diagonalized to find the vibrational spectrum of the inherent state; 3) the high-frequency part of the spectrum (related to the springs) is identified; 4) the scaling properties upon a density change of the remaining part of the spectrum is used to identify the pseudoisomorph. This method is very computationally demanding, but works well. Olsen et al. demonstrated that the it traces out pseudoisomorphs for the ASD and the flexible LJC models demonstrating good, though not perfect invariance of the dynamics probed via the incoherent intermediate scattering function. For comparison with later results, Fig. 2 shows the dynamics probed via the mean-square displacement, incoherent intermediate scattering function, and orientational time-autocorrelation function, for the state points identified in Olsen et al. (2016).

The present paper proposes simpler methods for generating pseudoisomorphs based on the scaling properties of the forces of a single configuration Schrøder (2022). This idea has been shown to work well not only for atomic systems like the Kob-Andersen binary Lennard-Jones mixture, but also for molecular models like the ASD model and the Lewis-Wahnstrom OTP model with rigid bonds Ingebrigtsen et al. (2012); Koperwas et al. (2020); Schrøder (2022); Sheydaafar et al. (2023).

The ASD model is a toy model of toluene Vrabec et al. (2001); Galbraith and Hall (2007); Schrøder et al. (2009a); Chopra et al. (2010a, b); Ingebrigtsen et al. (2012); Fragiadakis and Roland (2017); Santiago (2018); Dombrowski and Klotsa (2020). We simulated a system of 5000 ASD molecules defined as two different-sized spheres, a large (A) and a small (B) one Vrabec et al. (2001); Schrøder et al. (2009b); Milinkovic et al. (2013). The spheres interact via Lennard-Jones (LJ) potentials, and in the units defined by the large sphere ($\sigma_{AA}\equiv 1$, $\epsilon_{AA}\equiv 1$, and $m_{A}\equiv 1$) the other LJ-parameters are $\sigma_{AB}=0.894$, $\sigma_{BB}=0.788$, $\epsilon_{AB}=0.342$, $\epsilon_{BB}=0.117$, $m_{B}=0.195$. Bonds are modeled as harmonic springs with equilibrium length $0.584$ and spring constant $k=3000$.

The LJC model is a generic coarse-grained polymer model Bennemann et al. (1998); Aichele et al. (2003); Puosi and Leporini (2011); Shavit et al. (2013). We simulated 1000 10-bead LJC molecules. Non-bonded particles interact via a standard LJ potential, cutting and shifting the forces at $2.5\sigma$. All particles are of same type and the potential parameters are set to unity: $\sigma=1,\epsilon=1$. Bonds are modeled as harmonic springs with equilibrium length $1$ and spring constant $k=3000$.

All Molecular Dynamics simulations were performed in the $NVT$ ensemble with a Nose-Hoover thermostat in periodic boundary conditions, using RUMD that is an open-source molecular dynamics package optimized for GPU computing Bailey et al. (2017) (http://rumd.org).

Single-configuration force-based methods for generating isomorphs were introduced recently Schrøder (2022); Sheydaafar et al. (2023). The idea is the following. Given a configuration $\mathbf{R}_{1}$ at state point $(\rho_{1},T_{1})$, a uniform scaling to density $\rho_{2}$ is performed leading to $\mathbf{R}_{2}=\left(\rho_{1}/\rho_{2}\right)^{1/3}\mathbf{R}_{1}$. For molecules, two variants of this scaling can be applied: “center-of-mass scaling” where the molecular center-of-masses are scaled while all orientations and internal degrees of freedom are kept fixed, or “atomic scaling” where a uniform scaling is applied to all atoms thus modifying also the intramolecular bond lengths. After scaling a configuration, the forces associated with the two configurations, $\mathbf{F}(\mathbf{R}_{1})$ and $\mathbf{F}(\mathbf{R}_{2})$, are compared ($\mathbf{F}$ is the long vector of all forces). The temperature $T_{2}$ at the density $\rho_{2}$ is identified from the condition of invariant reduced forces. Specifically, $|\mathbf{\tilde{F}}(\mathbf{R}_{1})|=|\mathbf{\tilde{F}}(\mathbf{R}_{2})|$ implies Schrøder (2022)

If the two force vectors are parallel, this leads to the forces being identical in reduced units, $\mathbf{\tilde{F}}(\mathbf{R}_{1})=\mathbf{\tilde{F}}(\mathbf{R}_{2})\,$. Assuming this is representative of all relevant configurations, it follows that structure and dynamics are invariant in reduced units because the same reduced-unit equation of motion applies at the two state points Schrøder and Dyre (2014); Dyre (2018a). We do not expect the force vectors before and scaling to be perfectly parallel, however, but Eq. (3) can still be used. Reference Schrøder (2022) proposed that for force-based method to work well, both the Pearson and Spearman correlation coefficients of the force components should be larger than 0.95 (the latter is defined as the Pearson correlation coefficient of the rank of the data).

Different variants of the method are arrived at by different interpretations of what exactly $\mathbf{F}(\mathbf{R})$ represents, e.g., the forces on all the atoms or just on the center-of-mass forces. We consider below also a variant based on invariance of torques in reduced units, $\boldsymbol{\tilde{\tau}_{1}}=\boldsymbol{\tilde{\tau}_{2}}$, which leads to

The three methods are applied to the ASD model in Fig. 3, using center-of-mass scaling, i.e., $\mathbb{R}_{2,cm}=(\rho_{1}/\rho_{2})^{1/3}\mathbb{R}_{1,cm}$. The state point $(\rho_{1},T_{1})=(0.785,0.174)$ is used as reference and $\rho_{2}=0.856$, i.e., a 9% density increase is considered. Figure 3(a) shows a scatter plot of the atomic-force components before and after scaling for a single configuration; here $T_{2}=0.197$ is found by applying atomic forces in Eq. (3). Figure 3(b) shows a similar plot based on the “molecular” forces, i.e., the center-of-mass forces defined as the sum of all forces on the atoms of a given molecule. In this case the quite different $T_{2}=0.299$ is arrived at and a significantly better correlation is obtained. Finally, Fig. 3(c) shows the torque correlations of molecules of unscaled and scaled configurations. The correlation is here not far from that of the molecular forces; using Eq. (4) gives $T_{2}=0.310$.

For each of the three methods, Fig. 4 compares results for the (same) dynamics at the state points, which the different methods propose to give identical dynamics. While the methods in principle require only a single configuration, for carefully comparing the methods, $T_{2}$-values were obtained by averaging over 195 independent configurations. We find that the molecular-force method gives almost invariant dynamics, except at the lowest density (black curve) where negative pressure and phase separation is observed. The torque method gives similar, though slightly inferior results, while the atomic force method does not work well.

Next we turn to the 10-bead harmonic-spring LJC model for which one can use not only the three above methods, but also a fourth one based on invariant reduced segmental forces. These are defined as

in which $\mathbf{F}_{Seg,j}$ is the total force on particle $j$, i.e., including also the non-bonded interaction contributions, and $d_{j}$ and $d_{j+1}$ are the number of bonds that particles $j$ and $j+1$ are involved in. Clearly

The results are quite different from the ASD results. Thus the dynamics are to a good approximation invariant for the atomic- and segmental-force methods (Fig. 5 (g), (h), (i)), except at the lowest density from reference point $(\rho_{1},T_{1})=(1.00,0.700)$ at which the virial becomes negative. On the other hand, neither the molecular force method (d, e, f) nor the torque method (j, k, l) produce good pseudoisomorphs.

In Fig. 4 we applied the three force-based single-configuration methods to the ASD model and found best invariance of the dynamics using the molecular-force and, to a slightly less degree, torque methods. Figure 6 shows the results of applying the three methods at higher densities. The invariance of the dynamics is not as good as at lower densities (Fig. 4). The molecular-force method works best, but no method works as well as the method of Ref. Olsen et al., 2016, compare Fig. 2.

Figure 7(a) shows the distributions of bond lengths for simulations at the state points generated by the molecular-force method. The bonds are compressed when the density is increased, an effect that is not seen at the lower densities of Fig. 4. This means that when a configuration from the reference state point is scaled to a higher density, the scaled configuration is not representative of equilibrium configurations at the new state point because the bonds are too long, an issue that becomes more serious at higher densities.

To eliminate the effects of the harmonic bonds, we introduce the following procedure: For a given configuration fix the center-of-mass and orientation of each molecule. For this “constrained” system add a length $l$ to all bond lengths and minimize the potential energy with respect to $l$. We refer to this procedure as “quenching”, but note that the minimization only involves a single variable, the length $l$ added to all bonds. When applied to an unscaled configuration, the bond length distribution remains largely unchanged (black and orange curves in Fig. 7 (b)). Now, consider the scaling of the original configuration to higher density, with the aim to achieve a bond length distribution close to that of the equilibrium at the new state point (green curve). Using center-of-mass scaling of the molecules leaves the bond lengths unchanged (black curve), whereas atomic scaling lead to too small bond lengths (red curve). Applying the quenching procedure after center-of-mass scaling shifts the distribution to shorter bond lengths (from the black line to the blue line), and the average bond length approaches that of the equilibrium distribution at the higher density (green line). This effectively eliminates the non-scaling degrees of freedom thought to cause the poor invariance observed in Fig. 6.

Based on pairs of quenched configurations, one may again apply the atomic-force, molecular-force, and torque methods to generate state points with, hopefully, same dynamics. Specifically, we proceeded as follows. A single configuration is selected from an equilibrium simulation at the reference state point (density $\rho_{1}$). This configuration is scaled uniformly to the density of interest, $\rho_{2}$. Both scaled and unscaled configurations were then quenched as described above in order to eliminate the bond vibrational degrees of freedom. After this the relevant forces / torques were evaluated and $T_{2}$ determined from Eq. (3) and Eq. (4), respectively. Figure 8(a) shows the forces of a single scaled configuration versus those of the unscaled configuration before quenching, while (b) demonstrates better correlation after quenching; (c) and (d) show the correlations between the center-of-mass forces before and after quenching. The quench method leads to significantly better correlation, with correlation coefficient increasing from 0.850 to 0.934 in the atomic-force case and from 0.975 to 0.990 in the molecular (center-of-mass) force case.

The effect of the quenching procedure is tested in Fig. 9, which is analogous to Fig. 6 except that the temperatures $T_{2}$ are calculated based on quenched configurations. The best results are obtained with the molecular-force method, which was also the one that worked best in Fig. 4. For this method we find excellent collapse of the reduced center-of-mass mean-square displacement as a function of time, as well as of the center-of-mass incoherent intermediate scattering function. The directional time-autocorrelation function shows a slightly worse collapse, but is nevertheless significantly better than without quenching. Comparing the results of the torque method with and without quenching shows that quenching also here significantly improves the invariances. – Using atomic scaling gives the same results for the molecular force and torque methods after the system has been quenched (Table 1), while the atomic-force method gives results different from Fig. 9.

Genuine isomorphs do not exist in systems without strong virial potential-energy correlations, e.g., systems of molecules with harmonic bonds. For the ASD model and the flexible LJC model with harmonic bonds curves nevertheless exist along which the structure and dynamics are invariant to a good approximation, so-called pseudoisomorphs Olsen et al. (2016). These curves do not have invariant excess entropy, but otherwise behave much like isomorphs by having approximately invariant structure and dynamics in reduced units. It would be interesting to investigate whether some coarse-grained description corresponds to an excess entropy that is invariant along the pseudoisomorphs in this paper, but we have not attempted this and instead prioritized to directly search for methods that identify pseudoisomorphs.

A previously discussed method for tracing out pseudoisomorphs in the thermodynamic phase diagram works well, but is quite complicated to apply in practiceOlsen et al. (2016). The present paper explored the possibility of tracing out pseudoisomorphs based on the simple requirement of invariant length of the reduced force vector of a single configuration after scaling. Although the focus above was on the dynamics, we note that for all three methods discussed the structure is also invariant to a good approximation, both with and without quenching (Fig. 10).

Several closely related but different concepts have been introduced the last couple of decades in order to rationalize important findings of invariant structure and/or dynamics in the thermodynamic phase diagram. The present case of pseudoisomorphs for two realistic molecular models presents a good occasion for summarizing these concepts. One speaks about

• •

  Excess-entropy scaling when the (reduced-unit) transport coefficients are invariant along the lines of constant excess entropy Rosenfeld (1977); Dyre (2018b). Here DC signals the zero-frequency limit of, e.g., the diffusion coefficient or viscosity;

• •

  Density scaling when invariant transport coefficients are observed along specific lines in the phase diagram, often described as $\rho^{\gamma}/T=$Const. Alba-Simionesco et al. (2004); Roland et al. (2005);

• •

  Isochronal superposition when entire linear-response relaxation functions are invariant along specific lines in the phase diagram Roland et al. (2003);

• •

  Isomorph invariance when structure and dynamics are invariant along the lines of constant excess entropy then termed isomorphs Gnan et al. (2009); Schrøder and Dyre (2014); Dyre (2018a);

• •

  Pseudoisomorphs when structure and dynamics are invariant along lines that are not of constant excess entropy Veldhorst et al. (2014);

• •

  Isodynes when dynamics, but not structure, is invariant along specific lines Knudsen et al. (2021, 2024).

The most recent additions – pseudoisomorphs and isodynes – were introduced to describe liquids of more complex molecules, often with internal degrees of freedom. The vision is that some modification of isomorph theory, e.g., arrived at via a suitable coarse-graining focusing on the most important degrees of freedom, may apply also here. At this point in time, however, it is clear that more work is needed in terms of gathering data for many different models to arrive at a coherent picture. Even the relation between pseudoisomorphs and isodynes is not clear. For instance, it is an open question whether large molecules with many internal degrees of freedom are pseudoisomorphs or just isodynes, e.g., without structural invariants. This will depend on the relevant coarse-graining; in fact, systems previously classified as isodynes do have important parts of their structure factor being invariant along the same lines as the dynamics Knudsen et al. (2021, 2024).

For the harmonic-bond ASD model the best method for tracing out pseudoisomorphs is the molecular-force method with center-of-mass scaling and quenching of the harmonic bonds (Fig. 9). At low densities, the quenching can be skipped (Fig. 4). We conjecture that this method will work on all small molecules with harmonic bonds, assuming of course that the systems in question have pseudoisomorphs. The fact that quenching improves the tracing out of pseudoisomorphs is significant by indicating that the internal degrees of freedom of medium sized molecular systems can largely be scaled out.

The molecular-force method does not work well for the 10-bead harmonic-bond LJ-chain model for which the best method is the atomic-force method with atomic scaling (Fig. 5). That the atomic force method works best for long chains again is intuitively reasonable as a bead will not be aware that it is interacting with one from another chain or one on the same chain (apart from its nearest neighbors).

From a practical point of view one would like to have a single method that works for all molecular models with pseudoisomorphs, including the two models investigated here. The differences found can be interpreted as reflecting differences regarding which forces are important for the motion of the molecules. Despite the harmonic bonds, the ASD is relatively stiff and its motion is largely governed by the center-of-mass “molecular” force and the torque. If the scaling were perfect, both the molecular-force and the torque methods would lead to invariant dynamics. The fact that we find the molecular-force method to work better than the torque method can be interpreted as the molecular forces being more important than the torques for the motion of the molecules. It would be interesting to investigate whether molecular models exist for which the torque method works better than the molecular-force method. The molecular-force method does not work well for the LJC model. We interpret this as reflecting the fact that the center-of-mass motion does not determine the motion of this complex, flexible molecule, which should rather be be thought of in terms of the motion of its beads.

In summary, simple single-configuration force-based methods exist for tracing out pseudoisomorphs, which work as well as the complicated method of diagonalizing the Hessian before and after scaling.Olsen et al. (2016) No single-configuration method of general validity has been identified, however. More work is needed to investigate these and other methods. In this regard, we find it encouraging that the introduction of quenches after scaling improves all the methods studied.