The breakdown of the direct relation between the density scaling exponent and the intermolecular interaction potential for molecular systems with purely repulsive intermolecular forces

In this work, we question the generally accepted statement that the character of intermolecular interactions can be directly determined from the scaling exponent. Based on detailed studies of polyatomic molecular systems with precisely defined and purely repulsive intermolecular potential, we show that the value of the density scaling exponent evidently differs from the one predicted by the intermolecular virial-potential-energy correlation. Since the latter value directly results from the intermolecular potential, information on the interactions between molecules within the system cannot be immediately gained from the density scaling exponent value. Moreover, we suggest that the recently proposed"molecular force"method also returns the value that varies from the one scaling the dynamics. Finally, basing on our results, it might be deduced that the intramolecular interactions influence the density scaling value for real liquids.


MANUSCRIPT
Many liquids have a glass-forming ability that often makes them capable of becoming supercooled, i.e., they may remain a liquid even when their temperature is lowered below the melting point.The supercooling state can be achieved by isobaric cooling, where a constant pressure ensures that the melting point is fixed, or the liquid can be isothermally compressed.
In the second scenario, transformation from normal liquid to metastable state results from an increase in pressure, which causes the molecules to be more densely packed.In turn, lowering the temperature is equivalent to reducing the average kinetic energy of the molecules while increasing the pressure forces the molecules to get closer to each other, which decreases the volume in which the molecules can move.As a result of both processes, the liquid dynamics slowdown proceeds in the Arrhenius manner: the dependence of dynamic quantities on temperature (or pressure) is exponential. 1In the supercooling regime, however, the dynamics slowdown trend changes and becomes non-exponential. 2,3The critical factor determining the non-trivial nature of liquid dynamics is not fully defined 4,5 , but it is recognized that the temperature and volume effects are responsible.Which of them plays a more important role has been the subject of many works and widely discussed. 6It was found that the temperature and volume relative influence differs among the substances.Moreover, it was suggested 6 that instead of defining the factor with the dominant role, one should rather consider some function of these properties to represent the dynamics properly.A breakthrough that contributed to simplifying the description of liquid dynamics is the density scaling concept. 7,8cording to this concept, dynamics under different thermodynamic conditions can be mapped into a single relationship.The dynamical quantities can be expressed using the function of a single variable,  !, the product of temperature and volume raised to the power of the scaling exponent, from which inferences about the relative influence of temperature/volume on molecular dynamics can be made. 6,9Experimental studies 10 performed for more than 100 real glass-forming liquids showed that the value of the scaling exponent varies for different substances, indicating that  is a material parameter.A problem about what determines the value of the scaling exponent of molecular systems arose.Based on the theoretical research 11,12 and the results of the computer simulations 13,14 of the system described by the inverse power-law interaction potential, the physical interpretation of the scaling exponent was suggested.For the simple monoatomic system whose interactions were described using a potential in the form of an inverse power law 11  "#$ () =  %& +  ( -potential parameters,  -attractive background constant,  -distance), the scaling exponent turned out to be related to the potential as follows:  = /3.By analogy, the explanation of the scaling exponent value of the monoatomic system governed by the Lennard-Jones interaction potential  $' () = 4 1 , ( -the potential well depth,  -the distance at which the potential energy is zero) was made. 15At short distances where the potential is purely repulsive, up to the characteristic radius where it has a minimum,  $' ( ≤ 2 */, ) was fitted by a single effective  "#$ .The fit takes into account that the potential's attractive term modifies the curvature of the repulsive term.
Consequently, the effective  "#$ potential describes the repulsion with exponent  .//≈ 18, which explains  $' ≈ 6.However, at this point, the shortcomings of this consideration are worth mentioning.Only the part corresponding to the purely repulsive potential is taken into account, and this implies the misleading conclusion that repulsive forces exclusively determine the dynamics.If indeed the dynamics are governed only by repulsive forces, then both the considered  system and the corresponding system with  potential truncated at the minimum and shifted (i.e., Weeks-Chandler-Andersen, ) should have the same value of scaling exponent.However, it was shown that both the temperature dependences of the relaxation times and the scaling exponents of the systems are different from each other [16][17][18] , which emphasizes the importance of the role of attraction, which cannot be ignored.The second issue is that the effective potential method works only for the systems characterized by a potential distribution with spherical symmetry.In the case of real systems, it should be considered that molecules commonly consist of many atoms bound into different structures.
Additionally, but not necessarily, these atoms can be of various types.Due to these reasons, the potential around the real molecules is non-uniformly distributed in different directions.
Therefore, anisotropy is an inherent feature of real molecular systems that must be taken into account.In connection with this, it is highly challenging to correctly represent even the simplest anisotropic (i.e., diatomic) system by a point particle.The task requires specifying the reference point where the potential should have its source and considering that the value it takes depends not only on the distance but also on the orientation relative to some specific directions.It is also important that when considering the interaction of two molecules, each atom of the first one affects each atom of the second one according to a given potential, due to which it exerts forces in specific directions.Thus, the forces resulting from the impact of the first molecule on individual atoms of the second one depend on the mutual arrangement of both molecules.Therefore, the effective potential must correctly reflect the net forces between molecules.
A promising and systematically developed approach that seems to overcome the aforementioned problems is the isomorph theory [19][20][21][22][23] , according to which the scaling occurrence results from a strong correlation between potential energy  and the configurational contribution to pressure virial .In the case of the monoatomic  system, the simplicity of the potential implies that the  −  correlation is perfect, so its relationship is linear, and its slope coefficient is equal to /3.Considering more complex potentials whose curvature due to the attractive term varies with distance, such as  , the correlation is weakened.However, the theory states that it is sufficient for the correlation coefficient to be greater than or equal to 0.9, and then for such highly correlated systems (Roskilde-simple liquids), the scaling exponent can be determined as the slope of the  −  dependence.[26] Considering the rigid dumbbells example, it can be concluded that the correlation slope allows the determination of a scaling exponent for anisotropic systems.On the other hand, the results obtained for the monoatomic Gay-Berne system of ellipsoidal-like molecules contradict this statement. 27Similarly, the study of quasi-realistic systems 28 , i.e., rhombus-like  molecules with flexible bonds, , demonstrated that the scaling exponent cannot be determined from the dependency of the total  and  quantities because the contribution from their intramolecular interactions' terms significantly impairs the correlation. 29Therefore, the focus was on the virial and potential energy terms derived from intermolecular interactions, for which high correlation is present.Nevertheless, the slope of their relationship does not provide a scaling exponent.It was suspected that the potential curvature (an exponent of the effective  potential that could match the  potential within narrow intervals at given distances) changing with distance was responsible for this.Hence, in order to eliminate this factor, the consequent studies focused on the simplified system  "#$ , which differs from  in the way that all intermolecular interactions were described by the purely repulsive potential,  "#$ with exponent  = 12.The dynamical quantities of  "#$ were successfully scaled with an exponent equal to 4, i.e., with the value of the intermolecular virial and potential energy correlation slope.In the same work, a tetrahedral-shaped  system ( "#$ ) was also investigated.Surprisingly, the density scaling exponent for the  "#$ system was also found to be equal to 4, which suggested that in the case of  molecules, the value of the scaling exponent may not depend on the molecule's structure.This finding initiated further studies 30 of the  molecular systems with different potentials' parameters, as a result of which, an exact definition of the scaling exponent for polyatomic  molecules was obtained:  is the weighted average of the exponents of potentials describing the particular intermolecular interactions, where the weights are the potential energies corresponding to these interactions.
In this work, we challenge the generally accepted statement that the scaling exponent for molecular systems can be identified with the exponent of the potential describing intermolecular interactions.We conduct comprehensive research on the molecular systems that lack intermolecular attractive forces, investigating whether the scaling exponent is related to the  −  dependency.Moreover, we check whether, for  systems with arbitrary molecular structure, the scaling exponent takes the value of the exponent describing the  potential.Additionally, in the final stage of research, we apply the "molecular force" method, a recently developed tool based on isomorph theory, to trace the invariant dynamics and structure's phase-diagram curves to check if we can use it to determine the scaling exponent.
The systems we study in this work are constructed on the basis of  "#$ and  "#$ 31- 34 , with one additional atom attached to their structures, from now on referred to as  −  "#$ and  −  "#$ .In the case of  −  "#$ , the additional atom is attached to the atom located at the end of the longer diagonal of the rhombus, and for  −  "#$ , the extra atom is bonded to one of the side atoms.The schemes of the molecular structures can be seen in the insets of FIG. 1. &56 ,  = 0.355 , and their nonbonded, intermolecular interaction is described by the  "#$ potential with exponent  = 12, where  = 4 & .The interaction potential  "#$ is truncated at a distance  >?@ = 1.065 (which corresponds to a triple value of the  parameter) and shifted by the constant  = − "#$ ( >?@ ) to ensure it is equal to zero at the cut-off distance.

FIG. 1 Reduced relaxation times for mass centers of molecules as a function of the density scaling variable with an exponent determined from the correlation of the virial and potential energy corresponding to intermolecular interactions, see inset in
7][38] Systems consisting of 4000 molecules were examined in the pressure range from 0 to 100 MPa, customary for standard experiments.The procedure included cooling under isochoric conditions of  = 0.12, 0.13, 0.14 and  = 0.14, 0.15, 0.16 0][41] The  simulations for a single thermodynamic point consisted of two parts, the first of which was devoted to the equilibration of the system lasting at least ten times longer than the relaxation time  determined from the incoherent intermediate scattering function () 42 calculated for the molecular centers of mass using the wave vector corresponding to the maximum of the structure factor.In the second part, a simulation was performed during which data was collected.The relaxation time  was estimated as the time after the  function takes the value 1/.Subsequently, because the isomorph theory assumes that scaling can be observed when the quantities of the system are considered in reduced units, the relaxation times were reduced according to the formula , where  0 is the Boltzmann constant and  &56 is the mass of the molecule (the mass of the atom  <@ = 12.011 ). 22en, due to the lack of correlation between total virial and potential energy, we calculated the ,  contributions resulting from intermolecular interactions. 29,30,43Because of the simplicity of the  potential, the relationship is perfectly linear and has a slope coefficient equal to 4, see inset of FIG. 1a).The main parts of FIG. 1 show the relationship between the reduced relaxation times and the density scaling variable with exponent γ = 4, which does not collapse well onto one curve.Hence, the slope of the correlation does not give a scaling exponent.To check whether the relaxation times of the system can be scaled into one curve, we used the density scaling criterion that ( != ) B * G>59H@ and so  = − 1 . Based on the Vogel-Fulcher-Tammann [44][45][46] () function fit to the temperature dependence of the reduced relaxation times, we determined the temperature and volume conditions in which  * is constant ( * = 3, 3.5, 4), see insets of FIG. 2. As a result, we obtained linear volume-temperature relationships with the slopes equal to: 4.87 ± 0.03 for  −  "#$ , and 5.07 ± 0.22 for  −  "#$ .Consequently, in the main parts of FIG.
2, we confirmed that the reduced relaxation times can be scaled into a single curve with the exponent equal to the slope obtained from the density scaling criterion.Unless a slight deviation for long relaxation times may be observed for  −  "#$ .The above result confirms that the scaling exponent for  molecular systems with only one type of inter-molecular interactions may have a different value than the exponent of the  potential.The discrepancy between the scaling exponent and the  potential exponent confirms that the intermolecular interaction potential is not the only factor determining the scaling exponent for molecular systems.A possible explanation is that the scaling exponent is structure-dependent.Therefore, based on the scaling exponent, we cannot directly infer the characteristics of the intermolecular interaction potential, even for simplified molecular systems, i.e., with purely repulsive intermolecular interactions.

FIG. 2 Reduced relaxation times for mass centers of molecules as a function of the density scaling variable with an exponent determined from the density scaling criterion, see insets. 𝑅𝑀 − 𝑡𝑎𝑖𝑙
The molecular system described by the  potential for which the scaling exponent was obtained from the  −  correlation is a rigid dumbbell. 24In our work, inspired by the above example, we check whether the  −  relationship will allow us to determine the scaling exponent when we consider an analogous system but with a flexible bond.
Appropriately, we designed the  "#$ symmetric dumb-bell system, which are described by the same parameters of inter-molecular interactions as tailed systems.The only intramolecular interaction is a harmonic bond between two identical atoms characterized by the same length, and the same bond stretching constant as the bond for the tetrahedral system.
To consider the analogous pressure range, we performed the simulations for the system consisting of 5000 molecules in constant volume conditions  = 0.06, 0.07, 0.08 &56 .The simulation procedure and the analysis of the obtained data were identical, as in the case of systems with tails.Using the dependency of relaxation times calculated from  , we determined a scaling exponent for  "#$ equal to 3.09, see FIG. 3.For real molecules,  < 4 was observed for extremely polar systems. 47In the case of  "#$ , the system is non-polar, yet the scaling exponent takes a value relatively much smaller than 4.

FIG. 3 Reduced relaxation times for mass centers of molecules as a function of the density scaling variable with an exponent determined from the correlation of the virial and potential energy corresponding to intermolecular interactions on the left and
from the density scaling criterion on the right.The inset of section a) shows the scheme of the  !"# molecular structure, and the inset of section b) shows the density scaling criterion.
In FIG. 3, it can be seen that the relaxation times collapse well to one curve when using the exponent determined from the scaling criterion while the  −  correlation does not provide a scaling exponent.Similar result is presented in the work on the system of molecules of Lennard Jones chains with flexible bonds. 43,48It was suggested that for systems with harmonic bonds the phase diagram curves with invariant structure and dynamics be called pseudoisomorphs.This is because the invariance for molecular systems is approximate.
Namely, it can be found for the dynamics and structure of the centers of mass, but it is highly impaired for the atomic structure and dynamics, especially for systems with flexible bonds.
Due to the inability to use the virial and potential energy relationships to determine the scaling exponent of molecular systems, a search for new ways of tracking invariance began.
One of the recently developed approaches is the so-called "force method" 25,49 , according to which, based on a single configuration obtained from a system simulation, it is possible to find isomorphic thermodynamic points.Briefly speaking, let's consider a starting configuration containing information about the atoms' position in the system.Its isotropic rescaling allows for obtaining a target configuration whose structure is enlarged or shrunk in the same way in all directions.In reduced units, both configurations are structurally identical.
According to the theory, for both states (starting and target) to be isomorphic, the reduced forces exerted by the system on its atoms in both states should also be the same.In consequence, considering a system with the density  L9L@ in the starting configuration characterized by the vector  L9L@ = (  ,   , … ,   ) L9L@ of the  atom's positions, for a given target density  @<)=.@, we can create a target configuration  @<)=.@= g P +,+- P -./01- h */Q  L9L@ and estimate the target temperature from the equation  @<)=.@= RT -./01-VR |( +,+-)| g P +,+- where |()| is a length of the vector of forces exerted by the system on its atoms.As a result, the obtained states of the system at the initial and target conditions should be isomorphic.
Consequently, on their basis, the scaling exponent can be determined from the formula . This method has been successfully applied to atomic systems and systems of rigid molecules. 25,49w, we check whether this approach can be successfully applied to molecules with flexible bonds.The method will be tested for our system  −  "#$ .As the starting point, we choose the thermodynamic point  L9L@ = 0.14 9& " &56 ,  L9L@ = 80, and as the target state  @<)=.@= 0.12 9& " &56 that is about 17% more than the starting density.The starting point used in testing the molecular force method is marked in FIG.2a) with the star symbol.In the case of molecular systems, configuration rescaling cannot be made at the atomic level because it would change the bond lengths.Due to that, the arrangement of individual molecules is left the same as in the starting configuration, and only the positions of the centers of mass are rescaled.To take statistics into account, we analyzed 200 configurations from the equilibrated system simulation.In the molecular systems case, only forces from non-bonded interactions of surrounding molecules contribute because the forces resulting from intra-molecular interactions reduce.Hence, the method does not take into account intra-molecular interactions.Using the formula  @<)=.@= RT -./01-VR |( +,+-)| g P +,+- P -./01- h */Q  L9L@ , we determined temperature in the target state  @<)=.@= 196.58.Consequently, we simulated the system at the target thermodynamic point to verify the system's characteristics directly in reference to the starting state.In the FIG. 4, we compare the dynamics via  for the centers of mass as a function of reduced time for  −  "#$ .= 169.48≈ 170).As a result, the overestimation of the temperature estimated from the molecular force method for our system is at least 26 .
Summarizing, this work focused on verifying the generally accepted claim that the density scaling exponent is directly related to the intermolecular potential.Within it, we carried out detailed studies of molecular systems characterized by harmonic intramolecular interactions and no intermolecular attraction.We checked that the methods using the interaction potential, i.e., virial-potential energy correlation and molecular force, did not enable us to determine the scaling exponent for our systems.Therefore, relying solely on the intermolecular potential is insufficient to determine the scaling exponent of the molecular systems, and that, in turn, highlights the possible importance of their structural anisotropy and intramolecular interactions.In conclusion, it must be emphasized that based on the scaling exponent's value, we cannot directly infer the characteristics of the intermolecular potential of real liquids.
FIG. 1 Reduced relaxation times for mass centers of molecules as a function of the density scaling variable with an exponent determined from the correlation of the virial and potential energy corresponding to intermolecular interactions, see inset ina) section.Schematic representation of the molecular structure of  −  !"# on the left and  −  !"# on the right.
FIG. 2 Reduced relaxation times for mass centers of molecules as a function of the density scaling variable with an exponentdetermined from the density scaling criterion, see insets. −  !"# in the a) section and  −  !"# in the b) section.The point marked with the star symbol is used in testing the molecular force method.

FIG. 4
FIG. 4 Incoherent intermediate scattering function for centers of mass in function of the reduced time calculated for  −  !"# at points: starting, target, and determined from the density scaling criterion.The relaxation time is determined as the time the  function gets the value 1/ .

FIG. 4
FIG.4shows that the functions are different, so the starting and the target points are not characterized by invariant dynamics.The faster decay of the  function calculated for the target point confirms that the temperature determined using the molecular force method is overestimated.To quantify the overestimation, we placed the  function in the figure for one more thermodynamic point, which has invariant dynamics with the starting point and was