Effect of attractive interactions on the water-like anomalies of a core-softened model potential Shashank Pant, Tarun Gera, and Niharendu Choudhury Citation: The Journal of Chemical Physics 139, 244505 (2013); doi: 10.1063/1.4851478 View online: http://dx.doi.org/10.1063/1.4851478 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/24?ver=pdfcov Published by the AIP Publishing
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
THE JOURNAL OF CHEMICAL PHYSICS 139, 244505 (2013)
Effect of attractive interactions on the water-like anomalies of a core-softened model potential Shashank Pant,1 Tarun Gera,2 and Niharendu Choudhury3,a) 1
Department of Chemical Sciences, Indian Institute of Science Education and Research-Kolkata, Mohanpur-741252, India 2 Department of Chemistry, Indian Institute of Technology-Delhi, New Delhi, 110016, India 3 Theoretical Chemistry Section, Bhabha Atomic Research Centre, Mumbai-400 085, India
(Received 17 September 2013; accepted 4 December 2013; published online 31 December 2013) It is now well established that water-like anomalies can be reproduced by a spherically symmetric potential with two length scales, popularly known as core-softened potential. In the present study we aim to investigate the effect of attractive interactions among the particles in a model fluid interacting with core-softened potential on the existence and location of various water-like anomalies in the temperature-pressure plane. We employ extensive molecular dynamic simulations to study anomalous nature of various order parameters and properties under isothermal compression. Order map analyses have also been done for all the potentials. We observe that all the systems with varying depth of attractive wells show structural, dynamic, and thermodynamic anomalies. As many of the previous studies involving model water and a class of core softened potentials have concluded that the structural anomaly region encloses the diffusion anomaly region, which in turn, encloses the density anomaly region, the same pattern has also been observed in the present study for the systems with less depth of attractive well. For the systems with deeper attractive well, we observe that the diffusion anomaly region shifts toward higher densities and is not always enclosed by the structural anomaly region. Also, density anomaly region is not completely enclosed by diffusion anomaly region in this case. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4851478] I. INTRODUCTION
Although water is one of the most essential and abundant substances on earth, the origin of its numerous anomalous properties is not yet fully understood.1–3 It has been a challenge for scientists over the decades to correlate and understand the dynamic and thermodynamic anomalies of water with its microscopic local structure. One of the simplest and well-known anomalies of water is the density anomaly. The phase diagram of water has regions in the temperaturepressure (T-P) parameter space where the thermal expansion coefficient of water is negative, i.e., its specific volume increases when it is cooled. Water has many other anomalous properties2, 3 as well; notable among these are diffusion anomaly, compressibility anomaly, viscosity anomaly, etc. In the anomalous diffusivity region of the T-P (or T-ρ) plane, unlike a simple fluid, diffusion coefficient D of water increases on isothermal compression and outside this region diffusivity again follows the normal behavior of decreasing mobility on isothermal compression. It has been proposed by Stanley and coworkers4 that on increasing the pressure or density, a water molecule comes to close proximity in the interstitial region between the two hydrogen bonded water molecules. As a result, the hydrogen bonded structure of water becomes weak and consequently, it becomes easier for the hydrogen bonded water molecule to move before being connected to another molecule through a a) Author to whom correspondence should be addressed. Electronic
addresses: [email protected] and [email protected]. Also in Homi Bhabha National Institute. 0021-9606/2013/139(24)/244505/9/$30.00
new hydrogen bond. In order to understand these anomalies better, we need to find out their relationship with the local structure. In a pioneering work, Errington and Debenedetti5 have demonstrated that a well-defined correlation exists between structural orders and the anomalies of water. They have studied6 structural and orientational orders of model water using two order parameters, namely, translational order parameter t and orientational order parameter q assuming local geometry of water to be tetrahedral. A more general orientational order parameter Q6 applicable to any geometry has been formulated by Steinhardt et al.7 by considering only the nearest neighbors of a molecule. Using these order parameters for an atomistic model of water, Errington and Debenedetti5 have shown that a region in the P-T phase diagram exists where apart from density and diffusivity, structural and orientational parameters also show anomalous trend on isothermal compression. They have also found that these structural, thermodynamic, and dynamic anomalies constitute a cascade: thermodynamic (density) anomaly region is enclosed by dynamic (diffusion) anomaly region, which in turn is enclosed within the structural anomaly (anomalies related to order parameters t and q) region. Many investigations8–10 involving atomistic description of water have dealt with the anomalies of water and their locations in the T-P (or T-ρ) plane. Although it has been presumed initially that water-like anomalies can be observed only in fluids with hydrogenbonded network structure, many experimental and simulation investigations have indicated the presence of such anomalies even in many non-hydrogen-bonded liquids such as Te,11–18 Ga, S,12 Bi, Ge15 Te85 ,13 graphite,14 silica,15, 16 silicon,17 and
139, 244505-1
© 2013 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-2
Pant, Gera, and Choudhury
BeF2 .18 Extensive efforts are being made in recent years to understand the origin of such anomalies encountered in such a diverse set of systems encompassing both hydrogen-bonded and non-hydrogen-bonded fluids. In order to encompass all these anomalous liquids, it is essential to develop a more general, isotropic model capable of reproducing these anomalies or a subset of it. An isotropic model is convenient for understanding the physics behind these anomalies as it is easy to implement and interpret. In this context, Robinson and co-workers19 have observed that water can be viewed as a mixture of two components interacting with two different lengthscales. It is thought that an isotropic, simplified pair interaction with two lengthscales may be a suitable candidate for reproducing waterlike anomalies. In fact, long back in 1970, Hemmer and Stell20 successfully designed an interaction potential with two lengthscales (henceforth referred to as core-softened (CS) potentials) by conjugating a Lennard-Jones potential with a Gaussian peak. This potential has a repulsive core and a softening region as a shoulder or ramp. Stanley and coworkers21 and others22 used this potential in 2-dimensions to obtain various water-like anomalies. Since then, many attempts23–27 have been made to use a core-softened potential for three-dimensional liquid. Jagla28–31 successfully developed a discontinuous version of the two-scale potential, which reproduced32 water-like anomalies in 3-dimensions. A continuous version of the core-softened potential of the form given by Hemmer and Stell20 was then successfully parameterized by Barbosa and co-workers33 to obtain water-like anomalies in 3-dimensions. The cascade of anomalies involving structural, thermodynamic, and dynamic quantities of this core-softened potential is found to be of the same pattern as that obtained by Errington et al.5 for the SPC/E model34 of water. Both the continuous and discontinuous versions of the CS potential have no attractive part and therefore cannot reproduce liquid–gas or liquid–liquid critical points. All real molecules including inert gases have attractive interaction among themselves due to dispersion or electrostatic interactions and, therefore, presence of attractive region in the CS interaction will make it a more realistic potential. Recently, Barraz et al.35 introduced attractive interaction into the CS potential by considering a combination of four different CS potentials and studied existence of different anomalies by varying depth of the repulsive shoulder without changing the attractive part of the potential. In a new study, Silva et al.36 designed another potential form in which an additional Gaussian term is added into the usual Hemmer-Stell CS potential. They considered the effect of variation of attractive well depth without changing the shoulder part of the potential on the appearance and location of different water-like anomalies. They have observed that the anomalous region disappears if the attractive interaction is strong enough to reduce considerably the flux of particles between the two scales. Apart from the Hemmer and Stell form of the potential, many other forms of two-scale potential also show water-like anomalies. For example, Franzese and co-workers have come up37–40 with a new form of the two-scale potential using exponential and r−n terms and it is shown to possess all the water-like anomalies.
J. Chem. Phys. 139, 244505 (2013)
Similarly, Fomin et al.41–45 have put forward a different form of the two-scales potential using r−n term with tan hyperbolic function and have shown that this potential, depending on the parameters used, can reproduce,41 both water-like and silicalike anomalous behaviour. It is observed42, 43 that with the addition of attractive well in the potential, hierarchy of anomalies changes from water-like to silica-like. It is, thus, observed that the results for the anomalous properties obtained from different two-scale potentials are quite different and this difference arises probably because of the different functional forms used in these studies. In the present work, we intend to use the original form of the CS potential as given by Hemmer and Stell20 and incorporate attractive region in the CS potential by adjusting the width parameter of the Gaussian part of the CS potential. In this case both attractive and shoulder regions of the CS potential change simultaneously. Our aim is to investigate the effect of this change in the potential on different anomalous properties.
II. MODELS AND METHODS
The model of the fluid considered here consists of a system of N particles interacting with each other with a coresoftened potential of the form used by Hemmer and Stell,20 expressed as U (r)∗ = U (r)/ε
1 r − r0 2 σ 12 σ 6 + a exp − 2 =4 − , r r c σ (1)
where σ and ε are the size and energy parameters, respectively. In the above equation, a and c, respectively, are dimensionless parameters associated with the height and width of the Gaussian part of the potential. By tuning these parameters, we can obtain different potentials with varying depth and width of the well corresponding to the second lengthscale. All the N particles are initially placed in a fcc lattice in a cubic box of volume V corresponding to a number density ρ = N/V. We have constituted different systems with potentials of different depths of the attractive well by suitably choosing various potential parameters of Eq. (1). The value of a is chosen to be 5.0 (in reduced unit) in all the cases stud∗ ied here. We have considered r0 = 1.13 (in reduced unit, σ ); thus, the potential experiences a change in the slope around that region due to the weakening of the repulsive forces acting between the particles. We vary the 1/c2 (= C) parameter (see Eq. (1)) to make the attractive well deeper and in the present study we have considered C = 11, 12, and 20 and these three different cases will be referred hereinafter as systems A, B, and C, respectively. In Fig. 1 we have shown plots of the three potentials used in this study. As seen in this figure, with the increase in the value of the parameter C, both depth and width of the attractive well of the potential increase. We have employed molecular dynamics simulations46 in NVE ensembles with 1372 particles placed inside a cubic box with periodic boundary conditions in all three directions.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
3.0
*
U(r)
Repulsive shoulder (First lengthscale)
System A
c
1.2 1.0 1.6
Attractive well (Second lengthscale)
0.5
b
1.4
1.5 1.0
a
t
2.5 2.0
1.6
C = 11 C = 12 C = 20
System B
1.4
t
244505-3
1.2 1.0
0.0 -0.5 1.0
1.5
r*
2.0
0.2
2.5
0.3 0.4
0.5 0.6
1.6
t
1.4
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
FIG. 1. Three different core-softened potentials used in the present study [cf. ∗ Eq. (1) with r0 = 1.13, a = 1.0, and C = 11 (black line, system A), 12 (red line, system B), and 20 (green line, system C)].
A cutoff distance of 5σ has been used for potential truncation. All the parameters and quantities are expressed in reduced units as r∗ = r/σ , U∗ = U/ε, T∗ = kB T/ε, P∗ = Pσ 3 /ε, ρ ∗ = ρσ 3 , time step dt∗ = dt (ε/mσ 2 )1/2 . In the present study, we have used time step dt∗ = 0.005. In order to fix the temperature to a desired value, we have used velocity scaling method for around 50 000 steps and then another 50 000 steps without velocity scaling have been discarded for equilibration. After equilibration, 5 × 106 steps are used (without velocity scaling) as production run in which trajectories are stored in every 200 steps.
III. RESULTS AND DISCUSSION A. Translational order parameter
In order to calculate translational order in a system, we calculate translational order parameter of the system defined as5
ξ c |g (ξ ) − 1| dξ,
t=
(2)
0
where ξ = rρ 1/3 is the distance r in units of mean interparticle separation ρ −1/3 . Here, ρ is the average density, ρ = N/V, where N is the number of particles and V is the volume of the system. The quantity g(ξ ) in the above equation is the radial distribution function. ξ c is the cutoff distance, which in the present study is chosen to be ξ c = ρ 1/3 L/2, L being the side of the cubic simulation box with V = L3 . It is evident from the above equation that the translational order parameter t vanishes for a completely uncorrelated system (ideal gas) since g(ξ ) is unity and on the other hand, for a crystalline system g(ξ ) = 0 irrespective of r, and thus t assumes a large value. In Fig. 2, we have shown translational order parameter, t, as a function of density along different isotherms for three different systems considered in this work. These three systems differ from each other in the depths of the attractive wells of the potentials. In a normal fluid, translational order parameter along an isotherm generally increases upon compression or on increasing the density. In this case also at very high temperatures (T∗ > 0.45, not shown here), translational orders follow the normal behavior. But at lower temperatures below
1.2 System C
1.0 0.2 0.3 0.4 0.5 0.6 0.7 ρ*
FIG. 2. Translational order parameter t (cf. Eq. (2)) as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (top to bottom) for three different systems A, B, and C. The whole density range has been divided (shown by horizontal arrows in the top panel) into three zones a, b, and c, where zone b is the anomalous region and zones a and c are the normal regions before and after the anomalous region. The breaks in the curves at lower temperatures just before zone b is due to amorphous or glassy behaviour (see text).
T∗ = 0.45, it is found that in all the systems, t initially rises with increasing density showing a normal behavior but above a certain density ρ t-max (the density at which t is maximum), t starts to decrease with increasing ρ and thus ρ t-max defines the onset of an anomalous region, where translational order decreases on increasing density or compression. This anomalous region ends at a density ρ t-min (the density at which t is minimum), above which t again increases on compression showing normal behavior. It is important to notice that along some of the isotherms at lower temperatures, there are breaks in the curves. In this region, short-ranged liquid-like structure of water vanishes and long-ranged structural correlation sets in and thus calculation of t becomes difficult. As seen from Fig. 2, this (broken) region in the temperature-density plane increases with the increase of interparticle attraction in the potential from systems A to C. In order to look into the local orders in these regions, we have plotted radial distribution function (RDF) g(r) for different (T∗ , ρ ∗ ) points corresponding to these regions for the three different systems and are shown in Fig. 3. It is clearly observed that order is longranged with many peaks appearing throughout the system. It may be due to spontaneous crystallization as observed in an earlier study by Fomin et al.43 or the system in this region may be in amorphous or glassy states. Further, from the calculation of mean squared displacement (MSD) of the particles corresponding to these systems it is found that MSD does not increase with time (see inset of Fig. 3(a)) indicating solid or glass like environment in these cases. The anomalous behavior of t can be understood by analyzing the radial distribution functions as shown in Fig. 4 in the three different density regions as designated in Fig. 2 as a, b, and c. Systems in the regions a and c follow normal liquid behavior in the sense that translational order increases with increasing density, whereas in the region b, unlike normal liquid, translational order parameter t decreases with increasing
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-4
Pant, Gera, and Choudhury
T*=0.17, ρ*=0.35
0.05
3
System A
0.00
2
0
100 200 300
Time
1
g(r)
lous trend in t in this region can be explained as follows: an increase in peak height at r∗ = 1.0 from almost zero (in the region a) (see upper panels of Fig. 4) makes |g(r)−1| to decrease, and at the same time, decrease in peak height at around r∗ = 1.5 also causes t to decrease. In the lower panels of Fig. 4 we find that peaks at r∗ = 1 increase sharply (above 1) with density, making t value to increase. This region corresponds to the region-c of the t plots. The effect of increasing attraction causes the anomalous t-region to shift towards higher densities (compare density range in the 2nd panel of System C with that in System A of Fig. 4). As the density increases along an isotherm, particles move from the attractive well to the repulsive region of the potential. For the potential with deeper attractive well, more pressure (density) is required to push the particles from the second lengthscale (r∗ = 1.5) to the first lengthscale (r∗ = 1.0).
0.10
MSD
4
g(r)
J. Chem. Phys. 139, 244505 (2013)
0 4
T*=17, ρ*=0.35
3
System B
2 1
g(r)
0 4 3
T*=0.22 ρ*=0.40 ρ*=0.45 ρ*=0.50
System C
2 1 0 0
1
2
3
r*
4
5
6
7
B. Orientational order parameter
FIG. 3. Radial distribution functions, g(r) for the phase-space points corresponding to the broken regions (as shown in Fig. 2) in the translational order parameter vs. ρ ∗ curves of the three systems A, B, and C.
density. In Fig. 4, we have shown g(r) in these three different density regions, namely, a, b, and c for the three different systems A, B, and C considered here. In all the cases same general trend is observed. In the density range corresponding to region a of the translational order parameter plots (see Fig. 2), the RDF shows first major peak at around r∗ = 1.5 (see Fig. 4) corresponding to the attractive well of the interparticle potential. In the density range corresponding to region b, another peak starts developing at around r∗ = 1.0 corresponding to the repulsive part of the potential. The anoma-
In order to measure orientational order in the system, we calculate orientational order parameter as defined by Steinhardt et al.7 and implemented earlier.9, 47, 48 First, we need to consider k vectors rij that connects a central particle i to its k nearest neighbors j. These are also called bonds. The azimuthal θ and polar φ angles with respect to an arbitrary reference axis are computed for all the k bonds and then the average of the spherical harmonics Y lm (θ, φ) is calculated. The orientational order parameter pertaining to each particle, i, is then calculated from Qli =
System-B
System-A
g(r)
System-C ρ*=0.20 ρ*=0.25 ρ*=0.30 ρ*=0.35
T*=0.35 ρ*=0.20 ρ*=0.25 ρ*=0.30
2
m=l 4π 2 Ylm 2l + 1 m=−l
1/2 .
T*=0.30 ρ*=0.20 ρ*=0.25 ρ*=0.30 ρ*=0.35
1 Region-a
g(r)
0
Region-a
ρ*=0.40 ρ*=0.45 ρ*=0.50
2
Region-a
ρ*=0.40 ρ*=0.45 ρ*=0.50
ρ*=0.45 ρ*=0.50 ρ*=0.55
g(r)
1 Region-b
Region-b
Region-b
ρ*=0.50 ρ*=0.55 ρ*=0.60 ρ*=0.70
ρ*=0.55 ρ*=0.60 ρ*=0.65
ρ*=0.60 ρ*=0.65 ρ*=0.70
0 3 2 1
Region-c 0
1
2
r*
3
Region-c 1
2
r*
Region-c 3
1
2
r*
3
4
FIG. 4. Radial distribution functions, g(r), for the phase-space points corresponding to the three zones a, b, and c (as shown in Fig. 2) in the translational order parameter vs. ρ ∗ curves of the three systems A, B, and C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-5
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
0.25 System-A
0.36
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
0.28
Q6
0.36 0.32
0.15 0.25
T*
Q6
0.32
0.20
0.2
0.3
0.4
0.20 0.15
System B 0.28
System-B
0.3
0.5
0.6
0.2 System-C
0.36
0.1 0.32 0.40 0.48 0.56 0.64
Q6
ρ*
0.32 System C 0.2
0.3
0.4
ρ*
0.5
0.6
FIG. 6. Temperature of maximum density as a function of ρ ∗ for the three systems A, B, and C.
0.7
FIG. 5. Orientational order parameter Q6 (cf. Eq. (3)) as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (top to bottom) for three different systems A, B, and C. The breaks in the curves at lower temperatures are due to amorphous or glassy behaviour.
Finally, the average orientational order of the system of N particles is calculated for l = 6 by averaging over all the particles, viz., Q6 =
N 1 Q6i . N 1
(3)
The orientational parameter Q6 assumes a large value for a crystal and decreases with the decreasing order in the system. For a completely uncorrelated system like ideal gas, Q6 = √1k , and for a crystal, Q6 value depends on specific crystal structure and the number of nearest neighbors considered for the calculation. Like translational order parameters, value of Q6 also generally increases under isothermal compression in case of a normal liquid. In the present case of CS potential, however, the Q6 parameter as calculated from Eq. (3) increases (see Fig. 5) with increasing density up to ρ Q-max and then starts decreasing on further increase in density following normal behavior. Also, ρ Qmax in all the cases follow the trend ρ t-max ≤ ρ Q-max < ρ t-min , which implies that both Q6 and t behave anomalously in the density range between ρ Q-max and ρ t-min and this region is generally called structurally anomalous region. C. Density anomaly
Anomalous density region is generally observed in waterlike fluids. The temperature of maximum density (TMD) at ∂ρ )P = 0, corresponds to the minconstant pressure, i.e., ( ∂T ) = 0. It can be imum pressure along an isochore, ( ∂P ∂T ρ ) = easily seen by analyzing thermodynamic relation ( ∂V ∂T P ∂V ) ( ) that a minimum pressure with respect to tem−( ∂P ∂T V ∂P T perature represents a maximum in the density with respect to ) = 0. The TMD is the boundary of the thermoT, i.e., ( ∂V ∂T P dynamic anomaly region, where either decrease in tempera-
ture along an isobar or increase in pressure (density) along an isotherm results in decrease in density. Here, we calculate the temperature of minimum pressure at a constant density to find out the TMD and the results are shown in Fig. 6 for three different systems. It is found that increasing attraction causes not only to shift the anomalous region to higher densities but also at higher temperatures as well. For the systems A and B, in which depths of the attractive well are not much different, the anomalous regions (shown in Fig. 6 by black and red colors) are almost coinciding, but at a much higher value of the well depth (system C), this region (green line) shifts at higher densities and temperatures. Anomalous density trend can be observed in systems in which particles can stay at two preferential distances represented by the first and the second scales of the CS potential. For a normal liquid, percentage of particle at the closet distance decreases with increase in temperature. However, for the anomalous liquid, there is a region in the T-P plane where, as the temperature is increased, the percentage of particles at the closest distance increases. This is possible if particles initially residing at the second length scale moves to the first scale corresponding to the closest distance. The shifting of the anomalous region to higher densities and temperatures when attractive well is deeper can be explained in the following way: when attractive well is deeper, at a fixed density, higher temperature (energy) (as compared to the same for a potential with shallower attractive well) is required for the particles in the second length scale to migrate over to the first length scale and produce density anomaly. Similarly, at a fixed temperature, higher pressure or density is required for the system with deeper potential well (as compared to the same for the potential with not so deep potential well) to shift the particle from the second length scale to the first length scale. D. Diffusion anomaly
Water-like fluids generally have another anomaly, which is related to the dynamics of the system, and is known as diffusion anomaly. Self-diffusivity (D) of a fluid can be easily calculated from the slope of the Mean Square Displacements
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
0.25 0.20 0.15 0.10 0.05 0.00 0.20
J. Chem. Phys. 139, 244505 (2013)
System A
-1.2 -1.8
s2
System A
-2.4 -3.0
System B
-0.6 -1.2 -1.8 -2.4 -3.0
s2
0.15 0.10 0.05
T*=
0.00 0.2
0.3
0.4
0.15 D*
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
-0.6
ρ*
0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
0.5
System C
0.10 0.05 0.00 0.2
0.3
0.4
0.5
ρ*
0.6
0.7
System B 0.2
0.3
0.4
0.5
0.6
-0.6 -1.2 -1.8
s2
D*
D*
244505-6
-2.4 -3.0
System C 0.3
0.8
0.4
ρ*
0.5
0.6
0.7
FIG. 7. Diffusion coefficient D∗ as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (bottom to top) for three different systems A, B, and C.
FIG. 8. Excess entropy s2 as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (bottom to top) for three different systems A, B, and C.
(MSDs) of the fluid particles. The MSD is a good measure of translational mobility of a fluid in a given environment and, therefore, gives us information about the diffusion processes in the medium. At long time, MSD is related to the selfdiffusion coefficient (D) of the fluid through the well-known Einstein relation46, 49
tion box size.38 In Fig. 8, we have shown the variation of s2 as a function of density along various isotherms for the three different systems A, B, and C. In all the cases, we have found non-monotonic variations of s2 with density along different isotherms. As mentioned earlier, break in some of the curves indicates that in those regions, long ranged order has been established due to glassy or solid like states (see g(r) in Fig. 3). It is also observed that for the system C, which has higher attractive interactions, the non-monotonic behavior of s2 occurs in the higher temperature-density region. relation with various The density derivative of s2 has anomalies. In Fig. 9, we have shown ex = (∂s ex /∂ ln ρ)T
1 r 2 |r(t) − r(0)|2 1 lim = lim , 6 t→∞ t 6 t→∞
t
(4)
where r(t) is the position vector at time t. For normal liquids, diffusivity at a constant temperature decreases with increasing density. In Fig. 7, we have shown diffusivity as a function of density along various isotherms and it is found that in a certain density-temperature region, the diffusivity of the fluid increases with increasing density showing anomalous trend. From the comparison of the plots for the three different systems, it is found that for the system with deeper attractive potential well, diffusion anomaly region shifts towards higher densities.
4
0
-2 4
sex ≈ s2 = −2πρ ∫[g(r)lng(r) − g(r) + 1]r dr. 2
(5)
Here, g(r) is the radial distribution function. In the present work, we calculate s2 by numerical integration of the above equation using a cutoff distance of half of the simula-
Σ2
E. Excess entropy
System B
2 0
-2 4
0.2
0.3
0.4
ρ*
0.5
0.6 T*= 0.25 0.30 0.35 0.40 0.45
System C 2
Σ2
The excess entropy sex is the difference between the entropy of a real liquid and of an ideal gas at a similar condition of temperature and density. To determine the excess entropy one has to focus on all available configurations for the given system and relate it to ideal system. It is quite tedious to do. Hence, Errington et al.50 approximated the excess entropy by configurational entropy involving only two-body contributions, s2 , viz.,
T*= 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
2
Σ2
D=
0
-2 0.3
0.4
0.5
ρ*
0.6
0.7
FIG. 9. Derivative ( 2 ) of theexcess entropy (s2 ) as a function ofdensity ρ ∗ along differentisotherms for three different systems A, B, and C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
as a function of density along different isotherms for all the three systems A, B, and C. Errington et al.50 also proposed that the density dependence of s2 has relations with density, diffusion, and structural anomalies. It is proposed that the condition
for the density anomaly to appear follows the relation ex = (∂s ex /∂ ln ρ)T > 1 (green horizontal line in Fig. 9). Rosenfeld51 has shown that the logarithm of diffusion coefficient is directly proportional to the excess entropy and it implies that in order to this approximation to be valid, a non-monotonic behavior in excess entropy would result in a non-monotonic behavior in diffusivity. Using the empiri50 have also cal Rosenfeld’s parameterization, Errington et al. found the condition for diffusion anomaly to be ex > 0.42 (red horizontal line in Fig. 9). They also argued that structural anomaly region is determined by ex > 0 (black horizontal line in Fig. 9). It is, however, important to note that in some core-softened systems, Rosenfeld’s approximation relating diffusivity and excess entropy breaks down.45 However, in the present study we have observed that the regions bounded by these criteria as obtained from Fig. 9 for different anomalies approximately matches with the same obtained directly from the simulations (see density ranges corresponding to the anomalous regions in Figs. 2, 6, and 7). Therefore, in the present set of systems, no breakdown of the Rosenfeld’s approximation has been observed.
F. Order map
An interesting way of investigating structural order in liquids is to map state points into t-Q6 plane. Originally developed by Torquato and co-workers52 for hard sphere packing, Errington and Debenedetti5 used it to investigate structural order in water by mapping state points into t-q plane, where q is the tetrahedral order parameter. This representation is also
(a) System-A
1.50
t
1.35 1.20 Inaccessible region
1.05
T*=
0.90 1.50
(b) System-B
t
1.35 1.20 1.05
Dmin
0.45
Dmax
0.40
tmax tmin
0.35
Qmax
0.25
TMD min S2
0.20
S2
0.30
max
0.15 0.30
0.35
0.40
ρ*
0.45
0.50
FIG. 11. The relationship between different anomalies presented in the present study for the model system A. Translational order parameter t is anomalous in the region between the tmax and tmin lines (red), whereas orientational order parameter Qmax line (blue) is the starting point of the orientational anomaly region; above the density at a particular temperature corresponding to this curve, change in Q6 with density is anomalous. In the region between Qmax and tmin both translational and orientational order parameters are anomalous. This region is known as structural anomaly region. Change in diffusivity with density is anomalous in the region enclosed by Dmin and Dmax (dark yellow) lines and this region is called dynamic anomaly region. In the region enclosed by green s2 min and s2 max lines excess entropy change is anomalous. The region enclosed by the black TMD line (with symbols) is the density (thermodynamic) anomaly region.
used47 in case of core-softened model fluid, where the t-Q6 order map is used instead of t-q order map used for water. It has been found that state points accessible to liquid states in the order map fall into a two-dimensional region signifying that t and Q6 are, in general, independent. The order maps for all the model potentials considered here are shown in Fig. 10. In all the cases, t and Q6 values of the state points fall into a two-dimensional region; i.e., t and Q6 can be changed independently. As in the case of silica and water, for all the three systems studied here, we find an inaccessible region where no liquid state points can be found. The non-monotonic behavior of Q6 gives rise to loops along isotherms and the loop region indicates that state points with structural anomalies are lying on a narrow strip of the order map adjacent to the boundary between accessible and inaccessible regions.
System B
T*
(c) System-C
t
0.90 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0
Inaccessible region
0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
T*
244505-7
Dmax
0.45
Dmin
0.40
tmax tmin
0.35
Qmax
0.30 Inaccessible region
0.30
0.32
0.34
0.36
0.38
Q6 FIG. 10. The order map representations of the state points corresponding to different isotherms for three different systems. Upper panel: System A; Middle panel: System B, and Lower panel: System C. The slanted line in each plot is drawn to delineate inaccessible region from the accessible region.
0.25
TMD S2
0.20
S2
min
max
0.15 0.30
0.35
0.40
ρ*
0.45
0.50
FIG. 12. Same as that in Fig. 10(a) but for the system B.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-8
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
System C Dmin
0.45
Dmax
T*
0.40
tmax
0.35
tmin
0.30
Qmax
0.25
TMD min S2
0.20
max
S2
0.15 0.40 0.45 0.50 0.55 0.60 0.65 0.70
ρ*
FIG. 13. Same as that in Fig. 10(a) but for the system C.
G. Relative locations of different anomalies in the T-ρ plane
In this section we discuss relative positions of various anomalies discussed so far in the plane for all the three models studied here. In water, the region in the T-P or T-ρ plane corresponding to density anomaly is surrounded by diffusion anomaly region, which in turn is bounded by structural anomaly region, whereas in silica, the diffusion anomaly region contains the structural anomaly region, which in turn envelops the density anomaly region. The hierarchy of anomalies for three different systems A, B, and C is shown in Figs. 11–13, respectively. For system A, in which attraction among the particles is the minimum of the three systems, the density anomaly region is enclosed by the diffusion anomaly region, which in turn, is enclosed by structural anomaly region. All these anomalous regions are finally bounded by entropy anomaly region. It is important to mention that at lower density regions of the density anomaly region, systems are in amorphous or solid-like states at lower temperatures. As we increase the attractive well depth of the potential from systems A to B to C, it is found that the anomalous region as a whole shifts towards higher densities (compare Figs. 11 and 13). Not only that; in case of system C, the density anomaly region is not completely surrounded by diffusion anomaly region. Lower boundary of the diffusion anomaly region in this case cuts the TMD line. Thus, there is a part of the density anomaly region where there is no diffusion anomaly. The diffusion anomaly region in this case (system C, see Fig. 13) is shifted more (in comparison to boundaries of other anomalous regions) towards higher densities and it is coinciding with upper density boundary of the entropy anomaly region.
within the density anomaly region, where the system exists in an amorphous or a solid-like phase. As the depth of the attractive well of the potential is increased (from system A to B and C) by changing the parameter c of the potential (cf. Eq. (1)), the slope between the two length scales also changes and as a result the anomalous regions shift towards higher densities. The density anomaly region is shifted towards higher densities and temperatures. In a recent work,36 depth of the attractive well has been changed by adding an additional Gaussian to the potential presented in Eq. (1) and it has been shown that systems with sufficient attraction do not show any anomalous behavior. In the present case, however, we have found that the anomalous regions shift towards higher densities and temperatures. When potential well is deep, more temperature or pressure (density) is required for the particles to be shifted to the repulsive shoulder region from the well region. It is not very surprising, as in a recent study42 considering a different functional form of the potential, it is shown that with increasing attraction among the particles, relative positions of various anomalies change and water-like anomaly, in which structural anomaly region encloses diffusion anomaly region, changes to silica-like anomaly, in which relative position of the two above mentioned anomalies is interchanged. Thus, the present investigation demonstrates that the presence or absence of various anomalies and their positions (if anomalies are present at all) in the T-ρ or T-P plane depend on the minute details of the potential and its specific form. Although it is generally believed that the anomalous behavior of a fluid is a consequence of two lengthscales in the intermolecular potential, recently it has been demonstrated53 that anomalous behaviors can be observed even in fluids with only one local structure. It will be interesting to study the effect of attractive inter-particle interaction on the manifestation of various anomalies exhibited by such a model system. ACKNOWLEDGMENTS
We dedicate this paper to Dr. S. K. Ghosh on the occasion of his 64th birthday. It is a pleasure to thank Dr. S. K. Ghosh for his support and encouragement. One of the authors (S.P.) would like to acknowledge HBCSE, Mumbai for NIUS fellowship and thank Dr. Ashwani Kumar Tiwari of IISER, Kolkata and Savita Ladage of HBCSE, Mumbai for help and support and another author (T.G.) would like to thank the Indian Academy of Sciences for the Summer Research Fellowship. Thanks are due to the Computer Division, BARC for providing ANUPAM Supercomputing Facility. 1 F.
IV. CONCLUSIONS
In this work, we have considered one of the simplest forms of the continuous version of the core-softened potential and varied the width and the depth of the attractive well by varying the width of the Gaussian. For moderate depth of the attractive potential well (system A and B) we found waterlike cascade of anomalies: structural anomaly region contains diffusion anomaly region, which in turn, contains density anomaly region. It is also observed that there is a region
Franks, Water. A Comprehensive Treatise (Plenum, New York, 1973). Ludwig, “Water: From clusters to the bulk,” Angew. Chem., Int. Ed. 40, 1808–1827 (2001). 3 G. S. Kell, J. Chem. Eng. Data 20, 97–105 (1975). 4 P. A. Netz, F. W. Starr, M. C. Barbosa, and H. E. Stanley, Physica A 314, 470 (2002). 5 J. R. Errington and P. G. Debenedetti, Nature (London) 409, 318 (2001). 6 J. R. Errington, P. G. Debenedetti, and S. Torquato, J. Chem. Phys. 118, 2256 (2003). 7 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). 8 P. G. Debenedetti, Nature Phys. 9, 7–8 (2012); M. El Mekki Azouzi, C. Ramboz, J.-F. Lenain, and F. Caupin, ibid. 9, 38–41 (2012). 2 R.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-9 9 Z.
Pant, Gera, and Choudhury
Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. E 76, 051201 (2007). 10 A. Scala, F. W. Starr, E. La Nave, F. Sciortino, and H. E. Stanley, Nature (London) 406, 166–169 (2000); P. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, ibid. 360, 324–328 (1992). 11 H. Thurn and J. Ruska, J. Non-Cryst. Solids 22, 331 (1976). 12 G. E. Sauer and L. B. Borst, Science 158, 1567 (1967). 13 T. Tsuchiya, J. Phys. Soc. Jpn. 60, 227 (1991). 14 M. Togaya, Phys. Rev. Lett. 79, 2474 (1997). 15 M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Phys. Rev. E 66, 011202 (2002). 16 R. Sharma, S. N. Chakraborty, and C. Chakravarty, J. Chem. Phys. 125, 204501 (2006). 17 S. Sastry and C. A. Angell, Nature Matter 2, 739 (2003); T. Morishita, Phys. Rev. E 72, 021201 (2005). 18 C. A. Angell, R. D. Bressel, M. Hemmatti, E. J. Sare, and J. C. Tucker, Phys. Chem. Chem. Phys. 2, 1559 (2000). 19 C. H. Cho, S. Singh, and G. W. Robinson, Faraday Discuss. 103, 19–27 (1996). 20 P. C. Hemmer and G. Stell, Phys. Rev. Lett. 24, 1284 (1970). 21 M. R. Sadr-Lahijany, A. Scala, S. V. Buldyrev, and H. E. Stanley, Phys. Rev. Lett. 81, 4895 (1998). 22 P. J. Camp, Phys. Rev. E 71, 031507 (2005) 23 N. Choudhury and S. K. Ghosh, Phys. Rev. E 66, 021206 (2002). 24 G. Stell and P. C. Hemmer, J. Chem. Phys. 56, 4274–4286 (1972). 25 A. Scala, M. R. Sadr-Lahijany, N. Giovambattista, S. V. Buldyrev, and H. E. Stanley, J. Stat. Phys. 100, 97–106 (2000). L. Xu, S. V. Buldyrev, N. Giovambattista, and H. E. Stanley, Int. J. Mol. Sci. 11, 5184–5200 (2010). 26 P. Mausbach and H.-O. May, Proc. Appl. Math. Mech. 2, 531–532 (2003). 27 H. M. Gibson and N. B. Wilding, Phys. Rev. E 73, 061507 (2006). 28 E. A. Jagla, Phys. Rev. E 58, 1478 (1998). 29 E. A. Jagla, J. Chem. Phys. 111, 8980 (1999). 30 E. A. Jagla, Phys. Rev. E 63, 061501 (2001). 31 E. A. Jagla, Phys. Rev. E 63, 061509 (2001). 32 L. Xu, S. V. Buldyrev, N. Giovambattista, C. A. Angell, and H. E. Stanley, J. Chem. Phys. 130, 054505 (2009).
J. Chem. Phys. 139, 244505 (2013) 33 A.
B. de Oliveira, P. A. Netz, T. Colla, and M. C. Barbosa, J. Chem. Phys. 125, 124503 (2006); E. Salcedo, A. B. de Oliveira, N. M. Barraz, C. Chakravarty, and M. C. Barbosa, ibid. 135, 044517 (2011). 34 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269–6271 (1987). 35 N. M. Barraz, Jr., E. Salcedo, and M. C. Barbosa, J. Chem. Phys. 131, 094504 (2009). 36 J. N. da Silva, E. Salcedo, A. B. de Oliverira, and M. C. Barbosa, J. Chem. Phys. 133, 244506 (2010). 37 G. Franzese, J. Mol. Liq. 136, 267 (2007). 38 P. Vilaseca and G. Franzese, J. Chem. Phys. 133, 084507 (2010). 39 A. B. de Oliveira, G. Franzese, P. A. Netz, and M. C. Barbosad, J. Chem. Phys. 128, 064901 (2008). 40 P. Vilaseca and G. Franzese, J. Non-Cryst. Solids 357, 419–426 (2011). 41 Yu. D. Fomin and V. N. Ryzhov, Phys. Lett. A 375, 2181–2184 (2011). 42 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, Eur. Phys. J. Special Topics 216, 165–173 (2013). 43 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, J. Chem. Phys. 135, 234502 (2011) 44 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, Phys. Rev. E 87, 042122 (2013). 45 Yu. D. Fomin and V. N. Ryzhov, Phys. Rev. E 81, 061201 (2010). 46 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University, New York, 2004). 47 Z. Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. Lett. 95, 130604 (2005). 48 D. Bandyopadhyay, S. Mohan, S. K. Ghosh, and N. Choudhury, J. Phys. Chem. B 117, 8831–8843 (2013). 49 J.-P. Hansen and I. McDonald, Theory of Simple Liquids, 3rd ed. (Elsevier, New York, 2006). 50 J. R. Errington, T. M. Truskett, and J. Mittal, J. Chem. Phys. 125, 244502 (2006). 51 Y. Rosenfeld, J. Phys.: Condens. Matter 11, 5415 (1999). 52 S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000). 53 S. Prestipino, F. Saija, and G. Malescio, J. Chem. Phys. 133, 144504 (2010).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
THE JOURNAL OF CHEMICAL PHYSICS 139, 244505 (2013)
Effect of attractive interactions on the water-like anomalies of a core-softened model potential Shashank Pant,1 Tarun Gera,2 and Niharendu Choudhury3,a) 1
Department of Chemical Sciences, Indian Institute of Science Education and Research-Kolkata, Mohanpur-741252, India 2 Department of Chemistry, Indian Institute of Technology-Delhi, New Delhi, 110016, India 3 Theoretical Chemistry Section, Bhabha Atomic Research Centre, Mumbai-400 085, India
(Received 17 September 2013; accepted 4 December 2013; published online 31 December 2013) It is now well established that water-like anomalies can be reproduced by a spherically symmetric potential with two length scales, popularly known as core-softened potential. In the present study we aim to investigate the effect of attractive interactions among the particles in a model fluid interacting with core-softened potential on the existence and location of various water-like anomalies in the temperature-pressure plane. We employ extensive molecular dynamic simulations to study anomalous nature of various order parameters and properties under isothermal compression. Order map analyses have also been done for all the potentials. We observe that all the systems with varying depth of attractive wells show structural, dynamic, and thermodynamic anomalies. As many of the previous studies involving model water and a class of core softened potentials have concluded that the structural anomaly region encloses the diffusion anomaly region, which in turn, encloses the density anomaly region, the same pattern has also been observed in the present study for the systems with less depth of attractive well. For the systems with deeper attractive well, we observe that the diffusion anomaly region shifts toward higher densities and is not always enclosed by the structural anomaly region. Also, density anomaly region is not completely enclosed by diffusion anomaly region in this case. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4851478] I. INTRODUCTION
Although water is one of the most essential and abundant substances on earth, the origin of its numerous anomalous properties is not yet fully understood.1–3 It has been a challenge for scientists over the decades to correlate and understand the dynamic and thermodynamic anomalies of water with its microscopic local structure. One of the simplest and well-known anomalies of water is the density anomaly. The phase diagram of water has regions in the temperaturepressure (T-P) parameter space where the thermal expansion coefficient of water is negative, i.e., its specific volume increases when it is cooled. Water has many other anomalous properties2, 3 as well; notable among these are diffusion anomaly, compressibility anomaly, viscosity anomaly, etc. In the anomalous diffusivity region of the T-P (or T-ρ) plane, unlike a simple fluid, diffusion coefficient D of water increases on isothermal compression and outside this region diffusivity again follows the normal behavior of decreasing mobility on isothermal compression. It has been proposed by Stanley and coworkers4 that on increasing the pressure or density, a water molecule comes to close proximity in the interstitial region between the two hydrogen bonded water molecules. As a result, the hydrogen bonded structure of water becomes weak and consequently, it becomes easier for the hydrogen bonded water molecule to move before being connected to another molecule through a a) Author to whom correspondence should be addressed. Electronic
addresses: [email protected] and [email protected]. Also in Homi Bhabha National Institute. 0021-9606/2013/139(24)/244505/9/$30.00
new hydrogen bond. In order to understand these anomalies better, we need to find out their relationship with the local structure. In a pioneering work, Errington and Debenedetti5 have demonstrated that a well-defined correlation exists between structural orders and the anomalies of water. They have studied6 structural and orientational orders of model water using two order parameters, namely, translational order parameter t and orientational order parameter q assuming local geometry of water to be tetrahedral. A more general orientational order parameter Q6 applicable to any geometry has been formulated by Steinhardt et al.7 by considering only the nearest neighbors of a molecule. Using these order parameters for an atomistic model of water, Errington and Debenedetti5 have shown that a region in the P-T phase diagram exists where apart from density and diffusivity, structural and orientational parameters also show anomalous trend on isothermal compression. They have also found that these structural, thermodynamic, and dynamic anomalies constitute a cascade: thermodynamic (density) anomaly region is enclosed by dynamic (diffusion) anomaly region, which in turn is enclosed within the structural anomaly (anomalies related to order parameters t and q) region. Many investigations8–10 involving atomistic description of water have dealt with the anomalies of water and their locations in the T-P (or T-ρ) plane. Although it has been presumed initially that water-like anomalies can be observed only in fluids with hydrogenbonded network structure, many experimental and simulation investigations have indicated the presence of such anomalies even in many non-hydrogen-bonded liquids such as Te,11–18 Ga, S,12 Bi, Ge15 Te85 ,13 graphite,14 silica,15, 16 silicon,17 and
139, 244505-1
© 2013 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-2
Pant, Gera, and Choudhury
BeF2 .18 Extensive efforts are being made in recent years to understand the origin of such anomalies encountered in such a diverse set of systems encompassing both hydrogen-bonded and non-hydrogen-bonded fluids. In order to encompass all these anomalous liquids, it is essential to develop a more general, isotropic model capable of reproducing these anomalies or a subset of it. An isotropic model is convenient for understanding the physics behind these anomalies as it is easy to implement and interpret. In this context, Robinson and co-workers19 have observed that water can be viewed as a mixture of two components interacting with two different lengthscales. It is thought that an isotropic, simplified pair interaction with two lengthscales may be a suitable candidate for reproducing waterlike anomalies. In fact, long back in 1970, Hemmer and Stell20 successfully designed an interaction potential with two lengthscales (henceforth referred to as core-softened (CS) potentials) by conjugating a Lennard-Jones potential with a Gaussian peak. This potential has a repulsive core and a softening region as a shoulder or ramp. Stanley and coworkers21 and others22 used this potential in 2-dimensions to obtain various water-like anomalies. Since then, many attempts23–27 have been made to use a core-softened potential for three-dimensional liquid. Jagla28–31 successfully developed a discontinuous version of the two-scale potential, which reproduced32 water-like anomalies in 3-dimensions. A continuous version of the core-softened potential of the form given by Hemmer and Stell20 was then successfully parameterized by Barbosa and co-workers33 to obtain water-like anomalies in 3-dimensions. The cascade of anomalies involving structural, thermodynamic, and dynamic quantities of this core-softened potential is found to be of the same pattern as that obtained by Errington et al.5 for the SPC/E model34 of water. Both the continuous and discontinuous versions of the CS potential have no attractive part and therefore cannot reproduce liquid–gas or liquid–liquid critical points. All real molecules including inert gases have attractive interaction among themselves due to dispersion or electrostatic interactions and, therefore, presence of attractive region in the CS interaction will make it a more realistic potential. Recently, Barraz et al.35 introduced attractive interaction into the CS potential by considering a combination of four different CS potentials and studied existence of different anomalies by varying depth of the repulsive shoulder without changing the attractive part of the potential. In a new study, Silva et al.36 designed another potential form in which an additional Gaussian term is added into the usual Hemmer-Stell CS potential. They considered the effect of variation of attractive well depth without changing the shoulder part of the potential on the appearance and location of different water-like anomalies. They have observed that the anomalous region disappears if the attractive interaction is strong enough to reduce considerably the flux of particles between the two scales. Apart from the Hemmer and Stell form of the potential, many other forms of two-scale potential also show water-like anomalies. For example, Franzese and co-workers have come up37–40 with a new form of the two-scale potential using exponential and r−n terms and it is shown to possess all the water-like anomalies.
J. Chem. Phys. 139, 244505 (2013)
Similarly, Fomin et al.41–45 have put forward a different form of the two-scales potential using r−n term with tan hyperbolic function and have shown that this potential, depending on the parameters used, can reproduce,41 both water-like and silicalike anomalous behaviour. It is observed42, 43 that with the addition of attractive well in the potential, hierarchy of anomalies changes from water-like to silica-like. It is, thus, observed that the results for the anomalous properties obtained from different two-scale potentials are quite different and this difference arises probably because of the different functional forms used in these studies. In the present work, we intend to use the original form of the CS potential as given by Hemmer and Stell20 and incorporate attractive region in the CS potential by adjusting the width parameter of the Gaussian part of the CS potential. In this case both attractive and shoulder regions of the CS potential change simultaneously. Our aim is to investigate the effect of this change in the potential on different anomalous properties.
II. MODELS AND METHODS
The model of the fluid considered here consists of a system of N particles interacting with each other with a coresoftened potential of the form used by Hemmer and Stell,20 expressed as U (r)∗ = U (r)/ε
1 r − r0 2 σ 12 σ 6 + a exp − 2 =4 − , r r c σ (1)
where σ and ε are the size and energy parameters, respectively. In the above equation, a and c, respectively, are dimensionless parameters associated with the height and width of the Gaussian part of the potential. By tuning these parameters, we can obtain different potentials with varying depth and width of the well corresponding to the second lengthscale. All the N particles are initially placed in a fcc lattice in a cubic box of volume V corresponding to a number density ρ = N/V. We have constituted different systems with potentials of different depths of the attractive well by suitably choosing various potential parameters of Eq. (1). The value of a is chosen to be 5.0 (in reduced unit) in all the cases stud∗ ied here. We have considered r0 = 1.13 (in reduced unit, σ ); thus, the potential experiences a change in the slope around that region due to the weakening of the repulsive forces acting between the particles. We vary the 1/c2 (= C) parameter (see Eq. (1)) to make the attractive well deeper and in the present study we have considered C = 11, 12, and 20 and these three different cases will be referred hereinafter as systems A, B, and C, respectively. In Fig. 1 we have shown plots of the three potentials used in this study. As seen in this figure, with the increase in the value of the parameter C, both depth and width of the attractive well of the potential increase. We have employed molecular dynamics simulations46 in NVE ensembles with 1372 particles placed inside a cubic box with periodic boundary conditions in all three directions.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
3.0
*
U(r)
Repulsive shoulder (First lengthscale)
System A
c
1.2 1.0 1.6
Attractive well (Second lengthscale)
0.5
b
1.4
1.5 1.0
a
t
2.5 2.0
1.6
C = 11 C = 12 C = 20
System B
1.4
t
244505-3
1.2 1.0
0.0 -0.5 1.0
1.5
r*
2.0
0.2
2.5
0.3 0.4
0.5 0.6
1.6
t
1.4
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
FIG. 1. Three different core-softened potentials used in the present study [cf. ∗ Eq. (1) with r0 = 1.13, a = 1.0, and C = 11 (black line, system A), 12 (red line, system B), and 20 (green line, system C)].
A cutoff distance of 5σ has been used for potential truncation. All the parameters and quantities are expressed in reduced units as r∗ = r/σ , U∗ = U/ε, T∗ = kB T/ε, P∗ = Pσ 3 /ε, ρ ∗ = ρσ 3 , time step dt∗ = dt (ε/mσ 2 )1/2 . In the present study, we have used time step dt∗ = 0.005. In order to fix the temperature to a desired value, we have used velocity scaling method for around 50 000 steps and then another 50 000 steps without velocity scaling have been discarded for equilibration. After equilibration, 5 × 106 steps are used (without velocity scaling) as production run in which trajectories are stored in every 200 steps.
III. RESULTS AND DISCUSSION A. Translational order parameter
In order to calculate translational order in a system, we calculate translational order parameter of the system defined as5
ξ c |g (ξ ) − 1| dξ,
t=
(2)
0
where ξ = rρ 1/3 is the distance r in units of mean interparticle separation ρ −1/3 . Here, ρ is the average density, ρ = N/V, where N is the number of particles and V is the volume of the system. The quantity g(ξ ) in the above equation is the radial distribution function. ξ c is the cutoff distance, which in the present study is chosen to be ξ c = ρ 1/3 L/2, L being the side of the cubic simulation box with V = L3 . It is evident from the above equation that the translational order parameter t vanishes for a completely uncorrelated system (ideal gas) since g(ξ ) is unity and on the other hand, for a crystalline system g(ξ ) = 0 irrespective of r, and thus t assumes a large value. In Fig. 2, we have shown translational order parameter, t, as a function of density along different isotherms for three different systems considered in this work. These three systems differ from each other in the depths of the attractive wells of the potentials. In a normal fluid, translational order parameter along an isotherm generally increases upon compression or on increasing the density. In this case also at very high temperatures (T∗ > 0.45, not shown here), translational orders follow the normal behavior. But at lower temperatures below
1.2 System C
1.0 0.2 0.3 0.4 0.5 0.6 0.7 ρ*
FIG. 2. Translational order parameter t (cf. Eq. (2)) as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (top to bottom) for three different systems A, B, and C. The whole density range has been divided (shown by horizontal arrows in the top panel) into three zones a, b, and c, where zone b is the anomalous region and zones a and c are the normal regions before and after the anomalous region. The breaks in the curves at lower temperatures just before zone b is due to amorphous or glassy behaviour (see text).
T∗ = 0.45, it is found that in all the systems, t initially rises with increasing density showing a normal behavior but above a certain density ρ t-max (the density at which t is maximum), t starts to decrease with increasing ρ and thus ρ t-max defines the onset of an anomalous region, where translational order decreases on increasing density or compression. This anomalous region ends at a density ρ t-min (the density at which t is minimum), above which t again increases on compression showing normal behavior. It is important to notice that along some of the isotherms at lower temperatures, there are breaks in the curves. In this region, short-ranged liquid-like structure of water vanishes and long-ranged structural correlation sets in and thus calculation of t becomes difficult. As seen from Fig. 2, this (broken) region in the temperature-density plane increases with the increase of interparticle attraction in the potential from systems A to C. In order to look into the local orders in these regions, we have plotted radial distribution function (RDF) g(r) for different (T∗ , ρ ∗ ) points corresponding to these regions for the three different systems and are shown in Fig. 3. It is clearly observed that order is longranged with many peaks appearing throughout the system. It may be due to spontaneous crystallization as observed in an earlier study by Fomin et al.43 or the system in this region may be in amorphous or glassy states. Further, from the calculation of mean squared displacement (MSD) of the particles corresponding to these systems it is found that MSD does not increase with time (see inset of Fig. 3(a)) indicating solid or glass like environment in these cases. The anomalous behavior of t can be understood by analyzing the radial distribution functions as shown in Fig. 4 in the three different density regions as designated in Fig. 2 as a, b, and c. Systems in the regions a and c follow normal liquid behavior in the sense that translational order increases with increasing density, whereas in the region b, unlike normal liquid, translational order parameter t decreases with increasing
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-4
Pant, Gera, and Choudhury
T*=0.17, ρ*=0.35
0.05
3
System A
0.00
2
0
100 200 300
Time
1
g(r)
lous trend in t in this region can be explained as follows: an increase in peak height at r∗ = 1.0 from almost zero (in the region a) (see upper panels of Fig. 4) makes |g(r)−1| to decrease, and at the same time, decrease in peak height at around r∗ = 1.5 also causes t to decrease. In the lower panels of Fig. 4 we find that peaks at r∗ = 1 increase sharply (above 1) with density, making t value to increase. This region corresponds to the region-c of the t plots. The effect of increasing attraction causes the anomalous t-region to shift towards higher densities (compare density range in the 2nd panel of System C with that in System A of Fig. 4). As the density increases along an isotherm, particles move from the attractive well to the repulsive region of the potential. For the potential with deeper attractive well, more pressure (density) is required to push the particles from the second lengthscale (r∗ = 1.5) to the first lengthscale (r∗ = 1.0).
0.10
MSD
4
g(r)
J. Chem. Phys. 139, 244505 (2013)
0 4
T*=17, ρ*=0.35
3
System B
2 1
g(r)
0 4 3
T*=0.22 ρ*=0.40 ρ*=0.45 ρ*=0.50
System C
2 1 0 0
1
2
3
r*
4
5
6
7
B. Orientational order parameter
FIG. 3. Radial distribution functions, g(r) for the phase-space points corresponding to the broken regions (as shown in Fig. 2) in the translational order parameter vs. ρ ∗ curves of the three systems A, B, and C.
density. In Fig. 4, we have shown g(r) in these three different density regions, namely, a, b, and c for the three different systems A, B, and C considered here. In all the cases same general trend is observed. In the density range corresponding to region a of the translational order parameter plots (see Fig. 2), the RDF shows first major peak at around r∗ = 1.5 (see Fig. 4) corresponding to the attractive well of the interparticle potential. In the density range corresponding to region b, another peak starts developing at around r∗ = 1.0 corresponding to the repulsive part of the potential. The anoma-
In order to measure orientational order in the system, we calculate orientational order parameter as defined by Steinhardt et al.7 and implemented earlier.9, 47, 48 First, we need to consider k vectors rij that connects a central particle i to its k nearest neighbors j. These are also called bonds. The azimuthal θ and polar φ angles with respect to an arbitrary reference axis are computed for all the k bonds and then the average of the spherical harmonics Y lm (θ, φ) is calculated. The orientational order parameter pertaining to each particle, i, is then calculated from Qli =
System-B
System-A
g(r)
System-C ρ*=0.20 ρ*=0.25 ρ*=0.30 ρ*=0.35
T*=0.35 ρ*=0.20 ρ*=0.25 ρ*=0.30
2
m=l 4π 2 Ylm 2l + 1 m=−l
1/2 .
T*=0.30 ρ*=0.20 ρ*=0.25 ρ*=0.30 ρ*=0.35
1 Region-a
g(r)
0
Region-a
ρ*=0.40 ρ*=0.45 ρ*=0.50
2
Region-a
ρ*=0.40 ρ*=0.45 ρ*=0.50
ρ*=0.45 ρ*=0.50 ρ*=0.55
g(r)
1 Region-b
Region-b
Region-b
ρ*=0.50 ρ*=0.55 ρ*=0.60 ρ*=0.70
ρ*=0.55 ρ*=0.60 ρ*=0.65
ρ*=0.60 ρ*=0.65 ρ*=0.70
0 3 2 1
Region-c 0
1
2
r*
3
Region-c 1
2
r*
Region-c 3
1
2
r*
3
4
FIG. 4. Radial distribution functions, g(r), for the phase-space points corresponding to the three zones a, b, and c (as shown in Fig. 2) in the translational order parameter vs. ρ ∗ curves of the three systems A, B, and C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-5
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
0.25 System-A
0.36
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
0.28
Q6
0.36 0.32
0.15 0.25
T*
Q6
0.32
0.20
0.2
0.3
0.4
0.20 0.15
System B 0.28
System-B
0.3
0.5
0.6
0.2 System-C
0.36
0.1 0.32 0.40 0.48 0.56 0.64
Q6
ρ*
0.32 System C 0.2
0.3
0.4
ρ*
0.5
0.6
FIG. 6. Temperature of maximum density as a function of ρ ∗ for the three systems A, B, and C.
0.7
FIG. 5. Orientational order parameter Q6 (cf. Eq. (3)) as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (top to bottom) for three different systems A, B, and C. The breaks in the curves at lower temperatures are due to amorphous or glassy behaviour.
Finally, the average orientational order of the system of N particles is calculated for l = 6 by averaging over all the particles, viz., Q6 =
N 1 Q6i . N 1
(3)
The orientational parameter Q6 assumes a large value for a crystal and decreases with the decreasing order in the system. For a completely uncorrelated system like ideal gas, Q6 = √1k , and for a crystal, Q6 value depends on specific crystal structure and the number of nearest neighbors considered for the calculation. Like translational order parameters, value of Q6 also generally increases under isothermal compression in case of a normal liquid. In the present case of CS potential, however, the Q6 parameter as calculated from Eq. (3) increases (see Fig. 5) with increasing density up to ρ Q-max and then starts decreasing on further increase in density following normal behavior. Also, ρ Qmax in all the cases follow the trend ρ t-max ≤ ρ Q-max < ρ t-min , which implies that both Q6 and t behave anomalously in the density range between ρ Q-max and ρ t-min and this region is generally called structurally anomalous region. C. Density anomaly
Anomalous density region is generally observed in waterlike fluids. The temperature of maximum density (TMD) at ∂ρ )P = 0, corresponds to the minconstant pressure, i.e., ( ∂T ) = 0. It can be imum pressure along an isochore, ( ∂P ∂T ρ ) = easily seen by analyzing thermodynamic relation ( ∂V ∂T P ∂V ) ( ) that a minimum pressure with respect to tem−( ∂P ∂T V ∂P T perature represents a maximum in the density with respect to ) = 0. The TMD is the boundary of the thermoT, i.e., ( ∂V ∂T P dynamic anomaly region, where either decrease in tempera-
ture along an isobar or increase in pressure (density) along an isotherm results in decrease in density. Here, we calculate the temperature of minimum pressure at a constant density to find out the TMD and the results are shown in Fig. 6 for three different systems. It is found that increasing attraction causes not only to shift the anomalous region to higher densities but also at higher temperatures as well. For the systems A and B, in which depths of the attractive well are not much different, the anomalous regions (shown in Fig. 6 by black and red colors) are almost coinciding, but at a much higher value of the well depth (system C), this region (green line) shifts at higher densities and temperatures. Anomalous density trend can be observed in systems in which particles can stay at two preferential distances represented by the first and the second scales of the CS potential. For a normal liquid, percentage of particle at the closet distance decreases with increase in temperature. However, for the anomalous liquid, there is a region in the T-P plane where, as the temperature is increased, the percentage of particles at the closest distance increases. This is possible if particles initially residing at the second length scale moves to the first scale corresponding to the closest distance. The shifting of the anomalous region to higher densities and temperatures when attractive well is deeper can be explained in the following way: when attractive well is deeper, at a fixed density, higher temperature (energy) (as compared to the same for a potential with shallower attractive well) is required for the particles in the second length scale to migrate over to the first length scale and produce density anomaly. Similarly, at a fixed temperature, higher pressure or density is required for the system with deeper potential well (as compared to the same for the potential with not so deep potential well) to shift the particle from the second length scale to the first length scale. D. Diffusion anomaly
Water-like fluids generally have another anomaly, which is related to the dynamics of the system, and is known as diffusion anomaly. Self-diffusivity (D) of a fluid can be easily calculated from the slope of the Mean Square Displacements
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
0.25 0.20 0.15 0.10 0.05 0.00 0.20
J. Chem. Phys. 139, 244505 (2013)
System A
-1.2 -1.8
s2
System A
-2.4 -3.0
System B
-0.6 -1.2 -1.8 -2.4 -3.0
s2
0.15 0.10 0.05
T*=
0.00 0.2
0.3
0.4
0.15 D*
T*= 0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
-0.6
ρ*
0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
0.5
System C
0.10 0.05 0.00 0.2
0.3
0.4
0.5
ρ*
0.6
0.7
System B 0.2
0.3
0.4
0.5
0.6
-0.6 -1.2 -1.8
s2
D*
D*
244505-6
-2.4 -3.0
System C 0.3
0.8
0.4
ρ*
0.5
0.6
0.7
FIG. 7. Diffusion coefficient D∗ as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (bottom to top) for three different systems A, B, and C.
FIG. 8. Excess entropy s2 as a function of density ρ ∗ along different isotherms ranging from T∗ = 0.15 to T∗ = 0.45 (bottom to top) for three different systems A, B, and C.
(MSDs) of the fluid particles. The MSD is a good measure of translational mobility of a fluid in a given environment and, therefore, gives us information about the diffusion processes in the medium. At long time, MSD is related to the selfdiffusion coefficient (D) of the fluid through the well-known Einstein relation46, 49
tion box size.38 In Fig. 8, we have shown the variation of s2 as a function of density along various isotherms for the three different systems A, B, and C. In all the cases, we have found non-monotonic variations of s2 with density along different isotherms. As mentioned earlier, break in some of the curves indicates that in those regions, long ranged order has been established due to glassy or solid like states (see g(r) in Fig. 3). It is also observed that for the system C, which has higher attractive interactions, the non-monotonic behavior of s2 occurs in the higher temperature-density region. relation with various The density derivative of s2 has anomalies. In Fig. 9, we have shown ex = (∂s ex /∂ ln ρ)T
1 r 2 |r(t) − r(0)|2 1 lim = lim , 6 t→∞ t 6 t→∞
t
(4)
where r(t) is the position vector at time t. For normal liquids, diffusivity at a constant temperature decreases with increasing density. In Fig. 7, we have shown diffusivity as a function of density along various isotherms and it is found that in a certain density-temperature region, the diffusivity of the fluid increases with increasing density showing anomalous trend. From the comparison of the plots for the three different systems, it is found that for the system with deeper attractive potential well, diffusion anomaly region shifts towards higher densities.
4
0
-2 4
sex ≈ s2 = −2πρ ∫[g(r)lng(r) − g(r) + 1]r dr. 2
(5)
Here, g(r) is the radial distribution function. In the present work, we calculate s2 by numerical integration of the above equation using a cutoff distance of half of the simula-
Σ2
E. Excess entropy
System B
2 0
-2 4
0.2
0.3
0.4
ρ*
0.5
0.6 T*= 0.25 0.30 0.35 0.40 0.45
System C 2
Σ2
The excess entropy sex is the difference between the entropy of a real liquid and of an ideal gas at a similar condition of temperature and density. To determine the excess entropy one has to focus on all available configurations for the given system and relate it to ideal system. It is quite tedious to do. Hence, Errington et al.50 approximated the excess entropy by configurational entropy involving only two-body contributions, s2 , viz.,
T*= 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
2
Σ2
D=
0
-2 0.3
0.4
0.5
ρ*
0.6
0.7
FIG. 9. Derivative ( 2 ) of theexcess entropy (s2 ) as a function ofdensity ρ ∗ along differentisotherms for three different systems A, B, and C.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
as a function of density along different isotherms for all the three systems A, B, and C. Errington et al.50 also proposed that the density dependence of s2 has relations with density, diffusion, and structural anomalies. It is proposed that the condition
for the density anomaly to appear follows the relation ex = (∂s ex /∂ ln ρ)T > 1 (green horizontal line in Fig. 9). Rosenfeld51 has shown that the logarithm of diffusion coefficient is directly proportional to the excess entropy and it implies that in order to this approximation to be valid, a non-monotonic behavior in excess entropy would result in a non-monotonic behavior in diffusivity. Using the empiri50 have also cal Rosenfeld’s parameterization, Errington et al. found the condition for diffusion anomaly to be ex > 0.42 (red horizontal line in Fig. 9). They also argued that structural anomaly region is determined by ex > 0 (black horizontal line in Fig. 9). It is, however, important to note that in some core-softened systems, Rosenfeld’s approximation relating diffusivity and excess entropy breaks down.45 However, in the present study we have observed that the regions bounded by these criteria as obtained from Fig. 9 for different anomalies approximately matches with the same obtained directly from the simulations (see density ranges corresponding to the anomalous regions in Figs. 2, 6, and 7). Therefore, in the present set of systems, no breakdown of the Rosenfeld’s approximation has been observed.
F. Order map
An interesting way of investigating structural order in liquids is to map state points into t-Q6 plane. Originally developed by Torquato and co-workers52 for hard sphere packing, Errington and Debenedetti5 used it to investigate structural order in water by mapping state points into t-q plane, where q is the tetrahedral order parameter. This representation is also
(a) System-A
1.50
t
1.35 1.20 Inaccessible region
1.05
T*=
0.90 1.50
(b) System-B
t
1.35 1.20 1.05
Dmin
0.45
Dmax
0.40
tmax tmin
0.35
Qmax
0.25
TMD min S2
0.20
S2
0.30
max
0.15 0.30
0.35
0.40
ρ*
0.45
0.50
FIG. 11. The relationship between different anomalies presented in the present study for the model system A. Translational order parameter t is anomalous in the region between the tmax and tmin lines (red), whereas orientational order parameter Qmax line (blue) is the starting point of the orientational anomaly region; above the density at a particular temperature corresponding to this curve, change in Q6 with density is anomalous. In the region between Qmax and tmin both translational and orientational order parameters are anomalous. This region is known as structural anomaly region. Change in diffusivity with density is anomalous in the region enclosed by Dmin and Dmax (dark yellow) lines and this region is called dynamic anomaly region. In the region enclosed by green s2 min and s2 max lines excess entropy change is anomalous. The region enclosed by the black TMD line (with symbols) is the density (thermodynamic) anomaly region.
used47 in case of core-softened model fluid, where the t-Q6 order map is used instead of t-q order map used for water. It has been found that state points accessible to liquid states in the order map fall into a two-dimensional region signifying that t and Q6 are, in general, independent. The order maps for all the model potentials considered here are shown in Fig. 10. In all the cases, t and Q6 values of the state points fall into a two-dimensional region; i.e., t and Q6 can be changed independently. As in the case of silica and water, for all the three systems studied here, we find an inaccessible region where no liquid state points can be found. The non-monotonic behavior of Q6 gives rise to loops along isotherms and the loop region indicates that state points with structural anomalies are lying on a narrow strip of the order map adjacent to the boundary between accessible and inaccessible regions.
System B
T*
(c) System-C
t
0.90 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0
Inaccessible region
0.15 0.17 0.20 0.22 0.25 0.30 0.35 0.40 0.45
System A
T*
244505-7
Dmax
0.45
Dmin
0.40
tmax tmin
0.35
Qmax
0.30 Inaccessible region
0.30
0.32
0.34
0.36
0.38
Q6 FIG. 10. The order map representations of the state points corresponding to different isotherms for three different systems. Upper panel: System A; Middle panel: System B, and Lower panel: System C. The slanted line in each plot is drawn to delineate inaccessible region from the accessible region.
0.25
TMD S2
0.20
S2
min
max
0.15 0.30
0.35
0.40
ρ*
0.45
0.50
FIG. 12. Same as that in Fig. 10(a) but for the system B.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-8
Pant, Gera, and Choudhury
J. Chem. Phys. 139, 244505 (2013)
System C Dmin
0.45
Dmax
T*
0.40
tmax
0.35
tmin
0.30
Qmax
0.25
TMD min S2
0.20
max
S2
0.15 0.40 0.45 0.50 0.55 0.60 0.65 0.70
ρ*
FIG. 13. Same as that in Fig. 10(a) but for the system C.
G. Relative locations of different anomalies in the T-ρ plane
In this section we discuss relative positions of various anomalies discussed so far in the plane for all the three models studied here. In water, the region in the T-P or T-ρ plane corresponding to density anomaly is surrounded by diffusion anomaly region, which in turn is bounded by structural anomaly region, whereas in silica, the diffusion anomaly region contains the structural anomaly region, which in turn envelops the density anomaly region. The hierarchy of anomalies for three different systems A, B, and C is shown in Figs. 11–13, respectively. For system A, in which attraction among the particles is the minimum of the three systems, the density anomaly region is enclosed by the diffusion anomaly region, which in turn, is enclosed by structural anomaly region. All these anomalous regions are finally bounded by entropy anomaly region. It is important to mention that at lower density regions of the density anomaly region, systems are in amorphous or solid-like states at lower temperatures. As we increase the attractive well depth of the potential from systems A to B to C, it is found that the anomalous region as a whole shifts towards higher densities (compare Figs. 11 and 13). Not only that; in case of system C, the density anomaly region is not completely surrounded by diffusion anomaly region. Lower boundary of the diffusion anomaly region in this case cuts the TMD line. Thus, there is a part of the density anomaly region where there is no diffusion anomaly. The diffusion anomaly region in this case (system C, see Fig. 13) is shifted more (in comparison to boundaries of other anomalous regions) towards higher densities and it is coinciding with upper density boundary of the entropy anomaly region.
within the density anomaly region, where the system exists in an amorphous or a solid-like phase. As the depth of the attractive well of the potential is increased (from system A to B and C) by changing the parameter c of the potential (cf. Eq. (1)), the slope between the two length scales also changes and as a result the anomalous regions shift towards higher densities. The density anomaly region is shifted towards higher densities and temperatures. In a recent work,36 depth of the attractive well has been changed by adding an additional Gaussian to the potential presented in Eq. (1) and it has been shown that systems with sufficient attraction do not show any anomalous behavior. In the present case, however, we have found that the anomalous regions shift towards higher densities and temperatures. When potential well is deep, more temperature or pressure (density) is required for the particles to be shifted to the repulsive shoulder region from the well region. It is not very surprising, as in a recent study42 considering a different functional form of the potential, it is shown that with increasing attraction among the particles, relative positions of various anomalies change and water-like anomaly, in which structural anomaly region encloses diffusion anomaly region, changes to silica-like anomaly, in which relative position of the two above mentioned anomalies is interchanged. Thus, the present investigation demonstrates that the presence or absence of various anomalies and their positions (if anomalies are present at all) in the T-ρ or T-P plane depend on the minute details of the potential and its specific form. Although it is generally believed that the anomalous behavior of a fluid is a consequence of two lengthscales in the intermolecular potential, recently it has been demonstrated53 that anomalous behaviors can be observed even in fluids with only one local structure. It will be interesting to study the effect of attractive inter-particle interaction on the manifestation of various anomalies exhibited by such a model system. ACKNOWLEDGMENTS
We dedicate this paper to Dr. S. K. Ghosh on the occasion of his 64th birthday. It is a pleasure to thank Dr. S. K. Ghosh for his support and encouragement. One of the authors (S.P.) would like to acknowledge HBCSE, Mumbai for NIUS fellowship and thank Dr. Ashwani Kumar Tiwari of IISER, Kolkata and Savita Ladage of HBCSE, Mumbai for help and support and another author (T.G.) would like to thank the Indian Academy of Sciences for the Summer Research Fellowship. Thanks are due to the Computer Division, BARC for providing ANUPAM Supercomputing Facility. 1 F.
IV. CONCLUSIONS
In this work, we have considered one of the simplest forms of the continuous version of the core-softened potential and varied the width and the depth of the attractive well by varying the width of the Gaussian. For moderate depth of the attractive potential well (system A and B) we found waterlike cascade of anomalies: structural anomaly region contains diffusion anomaly region, which in turn, contains density anomaly region. It is also observed that there is a region
Franks, Water. A Comprehensive Treatise (Plenum, New York, 1973). Ludwig, “Water: From clusters to the bulk,” Angew. Chem., Int. Ed. 40, 1808–1827 (2001). 3 G. S. Kell, J. Chem. Eng. Data 20, 97–105 (1975). 4 P. A. Netz, F. W. Starr, M. C. Barbosa, and H. E. Stanley, Physica A 314, 470 (2002). 5 J. R. Errington and P. G. Debenedetti, Nature (London) 409, 318 (2001). 6 J. R. Errington, P. G. Debenedetti, and S. Torquato, J. Chem. Phys. 118, 2256 (2003). 7 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). 8 P. G. Debenedetti, Nature Phys. 9, 7–8 (2012); M. El Mekki Azouzi, C. Ramboz, J.-F. Lenain, and F. Caupin, ibid. 9, 38–41 (2012). 2 R.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
244505-9 9 Z.
Pant, Gera, and Choudhury
Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. E 76, 051201 (2007). 10 A. Scala, F. W. Starr, E. La Nave, F. Sciortino, and H. E. Stanley, Nature (London) 406, 166–169 (2000); P. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, ibid. 360, 324–328 (1992). 11 H. Thurn and J. Ruska, J. Non-Cryst. Solids 22, 331 (1976). 12 G. E. Sauer and L. B. Borst, Science 158, 1567 (1967). 13 T. Tsuchiya, J. Phys. Soc. Jpn. 60, 227 (1991). 14 M. Togaya, Phys. Rev. Lett. 79, 2474 (1997). 15 M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Phys. Rev. E 66, 011202 (2002). 16 R. Sharma, S. N. Chakraborty, and C. Chakravarty, J. Chem. Phys. 125, 204501 (2006). 17 S. Sastry and C. A. Angell, Nature Matter 2, 739 (2003); T. Morishita, Phys. Rev. E 72, 021201 (2005). 18 C. A. Angell, R. D. Bressel, M. Hemmatti, E. J. Sare, and J. C. Tucker, Phys. Chem. Chem. Phys. 2, 1559 (2000). 19 C. H. Cho, S. Singh, and G. W. Robinson, Faraday Discuss. 103, 19–27 (1996). 20 P. C. Hemmer and G. Stell, Phys. Rev. Lett. 24, 1284 (1970). 21 M. R. Sadr-Lahijany, A. Scala, S. V. Buldyrev, and H. E. Stanley, Phys. Rev. Lett. 81, 4895 (1998). 22 P. J. Camp, Phys. Rev. E 71, 031507 (2005) 23 N. Choudhury and S. K. Ghosh, Phys. Rev. E 66, 021206 (2002). 24 G. Stell and P. C. Hemmer, J. Chem. Phys. 56, 4274–4286 (1972). 25 A. Scala, M. R. Sadr-Lahijany, N. Giovambattista, S. V. Buldyrev, and H. E. Stanley, J. Stat. Phys. 100, 97–106 (2000). L. Xu, S. V. Buldyrev, N. Giovambattista, and H. E. Stanley, Int. J. Mol. Sci. 11, 5184–5200 (2010). 26 P. Mausbach and H.-O. May, Proc. Appl. Math. Mech. 2, 531–532 (2003). 27 H. M. Gibson and N. B. Wilding, Phys. Rev. E 73, 061507 (2006). 28 E. A. Jagla, Phys. Rev. E 58, 1478 (1998). 29 E. A. Jagla, J. Chem. Phys. 111, 8980 (1999). 30 E. A. Jagla, Phys. Rev. E 63, 061501 (2001). 31 E. A. Jagla, Phys. Rev. E 63, 061509 (2001). 32 L. Xu, S. V. Buldyrev, N. Giovambattista, C. A. Angell, and H. E. Stanley, J. Chem. Phys. 130, 054505 (2009).
J. Chem. Phys. 139, 244505 (2013) 33 A.
B. de Oliveira, P. A. Netz, T. Colla, and M. C. Barbosa, J. Chem. Phys. 125, 124503 (2006); E. Salcedo, A. B. de Oliveira, N. M. Barraz, C. Chakravarty, and M. C. Barbosa, ibid. 135, 044517 (2011). 34 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269–6271 (1987). 35 N. M. Barraz, Jr., E. Salcedo, and M. C. Barbosa, J. Chem. Phys. 131, 094504 (2009). 36 J. N. da Silva, E. Salcedo, A. B. de Oliverira, and M. C. Barbosa, J. Chem. Phys. 133, 244506 (2010). 37 G. Franzese, J. Mol. Liq. 136, 267 (2007). 38 P. Vilaseca and G. Franzese, J. Chem. Phys. 133, 084507 (2010). 39 A. B. de Oliveira, G. Franzese, P. A. Netz, and M. C. Barbosad, J. Chem. Phys. 128, 064901 (2008). 40 P. Vilaseca and G. Franzese, J. Non-Cryst. Solids 357, 419–426 (2011). 41 Yu. D. Fomin and V. N. Ryzhov, Phys. Lett. A 375, 2181–2184 (2011). 42 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, Eur. Phys. J. Special Topics 216, 165–173 (2013). 43 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, J. Chem. Phys. 135, 234502 (2011) 44 Yu. D. Fomin, E. N. Tsiok, and V. N. Ryzhov, Phys. Rev. E 87, 042122 (2013). 45 Yu. D. Fomin and V. N. Ryzhov, Phys. Rev. E 81, 061201 (2010). 46 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University, New York, 2004). 47 Z. Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. Lett. 95, 130604 (2005). 48 D. Bandyopadhyay, S. Mohan, S. K. Ghosh, and N. Choudhury, J. Phys. Chem. B 117, 8831–8843 (2013). 49 J.-P. Hansen and I. McDonald, Theory of Simple Liquids, 3rd ed. (Elsevier, New York, 2006). 50 J. R. Errington, T. M. Truskett, and J. Mittal, J. Chem. Phys. 125, 244502 (2006). 51 Y. Rosenfeld, J. Phys.: Condens. Matter 11, 5415 (1999). 52 S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000). 53 S. Prestipino, F. Saija, and G. Malescio, J. Chem. Phys. 133, 144504 (2010).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 93.180.53.211 On: Wed, 15 Jan 2014 09:41:46
