12028
J. Phys. Chem. B 2000, 104, 12028-12034
Liquid Structure via Cavity Size Distributions Pieter J. in ‘t Veld,† Matthew T. Stone,‡ Thomas M. Truskett, and Isaac C. Sanchez*,§ Department of Chemical Engineering, UniVersity of Texas at Austin, Austin, Texas 78712 ReceiVed: May 26, 2000; In Final Form: August 23, 2000
A new algorithm has been developed for determining the cavity size distribution in liquids. The computational method is based on energetic rather than geometric considerations. It is applicable to any liquid structure, including polymers, that can be simulated by molecular dynamic or Monte Carlo methods. To illustrate its utility, it has been applied to hard sphere (HS) and Lennard-Jones (LJ) fluids, SPC/E water, as well as to the two isomeric polyimides. To verify its ability to locate cavities and to determine cavity size frequency, it has also been applied with success to crystalline systems. Significant differences in liquid structure can be detected among HS, LJ, and SPC/E liquids, particularly at low densities where LJ and SPC/E liquids show evidence of “clustering.” The supercooled LJ liquid clearly reveals the development of a solid phase by the preferential appearance of tetrahedral and octahedral cavities. Cavitation in water under negative pressures can be detected. Differences in cavity size distributions are seen in two isomeric polyimides of identical chemical composition. Molecular dynamic simulations of carbon dioxide diffusion in these polyimides confirms that diffusion is faster in the isomer with the larger average cavity size.
I. Introduction The most successful experimental efforts to determine cavity size distributions in liquids and glasses have been undertaken in the areas of positron annihilation lifetime experiments1-3 and photochromic labeling.4-9 Although tremendous strides have been made, the cavity size distribution of a material is very difficult to access experimentally. On the modeling side, numerous Monte Carlo studies have been carried out on polymeric systems.10-16 Scaled particle theory (SPT), solvation theory, and the characterization of pore structures are some of the many applications in which cavity size distributions play a pivoting role. SPT predicts the reversible work of cavity formation and allows calculation of the chemical potential.17 The same reversible work of cavity formation is also believed to govern the process of solvation18,19 through the dependence on cavity size and shape.20 Cavity size distributions are also important in understanding mixing processes through the study of random porous microstructures.21 The pore structure characterization helps clarify the effect of random porous microstructures on the phase behavior of immiscible liquid mixtures. Furthermore, cavity size distributions are linked to pressure derivatives of solvation and association free energy and the resulting partial molar solvation and reaction volumes, which are a function of solute size and solvent density.22 Other more specific applications of cavity size distributions are found in the study of watermembrane systems, water, and other molecular liquids.23-26 The motivation for these studies is 2-fold. First, characterization of unoccupied volume is essential to understanding the microscopic structure of materials. As computer simulations and experimental techniques become more sophisticated, these investigations promise a fine-tuned, conceptual picture of †
E-mail:
[email protected]. E-mail:
[email protected]. § E-mail:
[email protected]. Fax: (512) 471 7060. Phone: (512) 471-1020 or 5014. ‡
molecular geometry. Second, it is believed that a detailed description of the cavity size distribution will aid in predicting bulk, thermodynamic and transport properties, in both simple and polymeric liquids. Cavity size distributions have already been used to interpret gas permeation, molecular mobility, diffusion, viscosity, crystallization, solubility, glass transition, and aging. Although the total unoccupied volume offers a quantitative description, it is the distribution, shape, and connectivity of the unoccupied space that control the underlying physics. Sastry et al.27,28 and Lee and Mattice29 introduce cavity sizing methods related to the one discussed in this paper. Sastry et al. present a refined geometrical approach using Delaunay and Voronoi tesselation (discussed earlier by Arizzi et al.30). They uniquely and exactly subdivide the entire three-dimensional space into tetrahedra. The vertexes of the fundamental space elements coincide with atom centers, creating a descriptor of the unoccupied volume composed by the cavities. Delaunay and Voronoi tessellations have been employed to describe polymeric and other microstructures.31-37 The cavity sizing procedure introduced by Lee and Mattice involves growing a bubble in the empty space between hard spheres. A cavity is defined as a bubble touching four neighboring hard spheres; overlap with the adjacent hard spheres is not allowed, forcing continuous center repositioning during bubble growth. This paper presents a new algorithm for obtaining cavity size information, using energetic rather than geometric means for cavity definition. The algorithm was tested and verified by application to crystalline structures with well-characterized cavities. Its use is demonstrated in four applications: (i) the evolution of the cavity size distribution in a supercooled Lennard-Jones liquid, (ii) a comparison between hard sphere, Lennard-Jones, and SPC/E water models, (iii) cavitation in SPC/E water under conditions of negative pressure, and (iv) the connection between the cavity size distribution and carbon dioxide diffusion in two isomeric polyimides.
10.1021/jp001934c CCC: $19.00 © 2000 American Chemical Society Published on Web 11/28/2000
Liquid Structure via Cavity Size Distributions
J. Phys. Chem. B, Vol. 104, No. 50, 2000 12029
II. Algorithm
by
A. Cavity Centering. We define a cavity as a spherical volume with a well-defined center, and that center as a local minimum in a repulsive particle energy field, which means that no force is exerted. A local minimum is found as follows: After a liquid structure is created by molecular dynamic or Monte Carlo methods, each particle (or united atom) is stripped of all its interactions with other particles and replaced by a LJ repulsive interaction characterized by the usual parameters σ and . Assumptions are made if LJ parameters are not available. For example, consider centering for a hard-sphere (HS) fluid: σi is set equal to σHS,i; i is arbitrarily chosen and identical for all particles. To find the local minimum, the artificially created repulsive force field is probed with a randomly introduced test particle with parameters σt and t. The choice of σt does not greatly influence the position of the local minimum, provided the system’s σ’s are within 50% of each other. t is a scale factor that does not affect the minimum’s position. The interaction energy of the test particle ψt with the other N particles is given by N
ψt ) 4
ti ∑ i)1
() σti
12
rti ) |b rt - b r i|
rti
(1)
where b ri is the position of particle i, and b rt is the position of the test particle. The parameters σti and ti are calculated by
1 σti ) (σt + σi) 2
ti ) xti
(2)
∇ψt
)-
ψt
48 ψt
∑ i)1
()
b r ti σti
N
ti
12
rti2 rti
(3)
The problem of finding the minimum is numerically solved by the multivariable Newton method combined with the steepest descent method. The Newton method requires linearization of eq 3, given by a Jacobian operating on ln ψt:
( )
∂2 ∂x2 2 ∇∇ ln ψt ) ∂ ∂y∂x ∂2 ∂z∂x
∂2 ∂x∂y ∂2 ∂y2 ∂2 ∂z∂y
∂2 ∂x∂z ∂2 ln ψ t ∂y∂z 2 ∂ ∂z2
(4)
Application of this Jacobian in a Newton root finding scheme results in an iterative process:
r t - fdamp,n(∇∇ ln ψt)-1 ∇ ln ψt b r t,new ) b
(5)
where the new position with a lower energy b rt,new is predicted from the current position b rt. The Newton damping factor fdamp,n, which typical value is between 0.5 and 1.0, ensures that the next guess will be inside the cavity and is adjusted during execution. However, if the determinant or a diagonal term of the Jacobian is negative, this scheme could yield a local maximum instead of a local minimum. In this case the method of steepest descent is applied. Cruder than the linear Newton method, it uses the gradient directly to find a local minimum
(6)
Typical values of the steepest descent damping factor fdamp,s are between 0.0005 and 0.001. A switch back to the Newton method occurs when the diagonal terms and the determinant of the Jacobian are positive. A convergence rate measure ξ is used to moderate either fdamp,n or fdamp,s (fdamp below):
ξ)
ψt - ψt,new ψt
(7)
The damping factor fdamp need not be adjusted when ξ g 0.25. Slow convergence needs an increase of fdamp. A negative ξ indicates a nonvalid move up the gradient and away from the cavity center; fdamp is decreased. When the condition ξ < 0 occurs five times in a row, a random move on a sphere with a 0.01σ radius is made to avoid gridlock. Adjustment of fdamp results in ψt,new and hence ξ recalculation. In summary, the rules used to control fdamp are
fdamp,new ) 0.75fdamp fdamp,new ) 1.25fdamp no change
ξ
