Calculation of interfacial properties using molecular simulation with the reaction field method: Results for different water models

Second inflection point of supercooled water surface tension induced by hydrogen bonds: A molecular-dynamics study

Surface tension of supercooled water is a fundamental property in various scientific processes. In this study, we perform molecular dynamics simulations with the TIP4P-2005 model to investigate the surface tension of supercooled water down to 220 K. Our results show a second inflection point (SIP) in the surface tension at temperature TSIP ≈ 267.5 ± 2.3 K. Using an extended IAPWS-E functional fit for the water surface tension, we calculate the surface excess internal-energy and entropy terms of the excess Helmholtz free energy. Similar to prior studies [Wang et al., Phys. Chem. Chem. Phys. 21, 3360 (2019); Gorfer et al., J. Chem. Phys. 158, 054503 (2023)], our results show that the surface tension is governed by two driving forces: a surface excess entropy change above the SIP and a surface excess internal-energy change below it. We study hydrogen-bonding near the SIP because it is the main cause of water's anomalous properties. With decreasing temperature, our results show that the entropy contribution to the surface tension reaches a maximum slightly below the SIP and then decreases. This is because the number of hydrogen bonds increases more slowly below the SIP. Moreover, the strengths and lifetimes of the hydrogen bonds also rise dramatically below the SIP, causing the internal-energy term to dominate the excess surface free energy. Thus, the SIP in the surface tension of supercooled TIP4P-2005 water is associated with an increase in the strengths and lifetimes of hydrogen bonds, along with a decrease in the formation rate (#/K) of new hydrogen bonds.

Effects of molecular size and orientation on the interfacial properties and wetting behavior of water/n-alkane systems: a molecular-dynamics study

Molecular dynamics simulations (MD) are performed to study the interfacial structure/tension and wetting behavior of water/n-alkane systems (water/nC5 to water/nC16 where nCx = CxH(2x + 2)). In particular, we study complete-to-partial wetting transitions by changing the n-alkane chain length (N_C) at a constant temperature, T = 295 K. Simulations are carried out with a united-atom TraPPE model for n-alkanes and the TIP4P-2005 model of water. Simulation results are in excellent agreement with the initial spreading coefficients and contact angles calculated using experimental values of the surface and interfacial tensions. In addition, it has been determined that water/(nC5-nC7) and water/(nC8-nC16), respectively, exhibit complete and partial initial wetting modes. Simulations show that the interfacial structures of water/(nC5-nC7) are different from water/(nC8-nC16) systems. In the latter, water preferentially orients near the interface to increase the number of hydrogen bonds and the charge and mass densities. Moreover, the orientation of n-alkane molecules at water/(nC8-nC16) interfaces has a long-range persistence, resulting in layered structures that increase with N_C. In addition, simulation results of the orientational order parameter S_z show alignment behavior of the n-alkane molecules with respect to the interfaces. Simulations predict that the central segments of n-alkane are strongly packed in the interfaces while the end segments (methyl groups) form smaller peaks in the outer edge of the layer. This observation confirms the "horseshoe" or "C-shaped" structure of n-alkane molecules in the water/n-alkane interfaces. At constant temperature, the interface widths of both water and the n-alkanes decrease with increasing n-alkane molecular length. These results suggest that increasing the n-alkane chain length affects the water/n-alkane interfacial properties in a manner similar to that of cooling.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

Molecular dynamics of liquid-liquid equilibrium and interfacial properties of aqueous solutions of methyl esters

In this work, the liquid-liquid phase equilibria and interfacial properties of methyl ester + water binary mixtures are determined at atmospheric pressure and from 278 to 358 K combining the direct coexistence technique and molecular dynamics simulations. Methyl esters are modelled using new parametrizations based on the united atom TraPPE model force field proposed recently by us [E. Feria, J. Algaba, J. M. Míguez, A. Mejía, P. Gómez-Álvarez and F. J. Blas, Phys. Chem. Chem. Phys., 2019, 22, 4974-4983] that are able to predict the vapour-liquid interfacial properties of pure methyl esters with high accuracy. In the case of water, we consider the well-known TIP4P/2005 model, the most popular rigid and non-polarizable model to describe the interfacial properties of pure water. The simulations are performed using the direct coexistence technique in the isothermal-isobaric or NP_zT ensemble in combination with molecular dynamics. We obtain density profiles, temperature-densities and temperature-composition projections of the phase diagrams, and interfacial tensions. The liquid-liquid interfacial tension is calculated from the normal and tangential components of the pressure tensor according to the mechanical virial route. We pay attention particularly to the ability of the molecular models in predicting the experimental behavior of the systems. Simulation results are able to account for the liquid-liquid phase equilibria of these binary mixtures, in good agreement with the experimental data taken from the literature. Unfortunately, experimental values for interfacial tensions are substantially overestimated by predictions from computer simulations in all cases. To our knowledge, this is the first time that the liquid-liquid phase equilibrium and interfacial properties of methyl ester + water mixtures have been predicted from computer simulations.

Preferential Orientations and Anomalous Interfacial Tensions in Aqueous Solutions of Alcohols

Literature studies on interfacial tension versus temperature between normal alcohols and water show that it increases with temperature and exhibits a maximum value at a given temperature depending on the molecular weight of the alcohol. This very unusual behavior is supposedly accompanied by the formation of monolayers of alcohol molecules oriented preferentially at the interface, a structural issue not confirmed until now. We use molecular-based models for water and alcohols in combination with molecular dynamics simulations to provide physical insights, from a molecular perspective, into the structural and thermodynamic behavior at the liquid-liquid interfaces of aqueous solutions of alcohols.

Calculation of the interfacial tension of the graphene-water interaction by molecular simulations

We report the calculation of the solid-liquid interface tension of the graphene-water interaction by using molecular simulations. Local profiles of the interfacial tension are given through the mechanical and thermodynamic definitions. The dependence of the interfacial tension on the graphene area is investigated by applying both reaction field and Ewald summation techniques. The structure of the interfacial region close to the graphene sheet is analyzed through the profiles of the density and hydrogen bond number and the orientation of the water molecules. We complete this study by plotting the profiles of the components of the pressure tensor calculated by the Ewald summation and reaction field methods. We also investigate the case of a reaction field version consisting in applying a damped shifted force in the case of the calculation of the pressure components.

An excess electron at polyethylene/vacuum interfaces using a reaction-field technique

An excess electron in amorphous and crystalline polyethylene-vacuum interfaces. A precise reaction-field method is used to compute the surface states.

The liquid-vapor equilibria of TIP4P/2005 and BLYPSP-4F water models determined through direct simulations of the liquid-vapor interface

In this paper, the surface tension and critical properties for the TIP4P/2005 and BLYPSP-4F models are reported. A clear dependence of surface tension on the van der Waals cutoff radius (rvdw) is shown when van der Waals interactions are modeled with a simple cutoff scheme. A linear extrapolation formula is proposed that can be used to determine the infinite rvdw surface tension through a few simulations with finite rvdw. A procedure for determining liquid and vapor densities is proposed that does not require fitting to a profile function. Although the critical temperature of water is also found to depend on the choice of rvdw, the dependence is weaker. We argue that a rvdw of 1.75 nm is a good compromise for water simulations when long-range van der Waals correction is not applied. Since the majority of computational programs do not support rigorous treatment of long-range dispersion, the establishment of a minimal acceptable rvdw is important for the simulation of a variety of inhomogeneous systems, such as water bubbles, and water in confined environments. The BLYPSP-4F model predicts room temperature surface tension marginally better than TIP4P/2005 but overestimates the critical temperature. This is expected since only liquid configurations were fit during the development of the BLYPSP-4F potential. The potential is expected to underestimate the stability of vapor and thus overestimate the region of stability for the liquid.

Computer modelling of the surface tension of the gas–liquid and liquid–liquid interface

This review presents the state of the art in molecular simulations of interfacial systems and of the calculation of the surface tension from the underlying intermolecular potential.

The zero-multipole summation method for estimating electrostatic interactions in molecular dynamics: analysis of the accuracy and application to liquid systems

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

Simple and accurate scheme to compute electrostatic interaction: zero-dipole summation technique for molecular system and application to bulk water

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 Å for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.

Non-Ewald methods: theory and applications to molecular systems

Several non-Ewald methods for calculating electrostatic interactions have recently been developed, such as the Wolf method, the reaction field method, the pre-averaging method, and the zero-dipole summation method, for molecular dynamics simulations of various physical systems, including biomolecular systems. We review the theories of these approaches and their potential applications to molecular simulations, and discuss their relationships.

Effect of atom- and group-based truncations on biomolecules simulated with reaction-field electrostatics

The performance of the reaction-field method of electrostatics is tested in molecular dynamics simulations of protein human interleukin-4 and a short DNA fragment in explicit solvent. Two truncation schemes are considered: one based on the position of atomic charges in water molecules and the other on the position of groups of charges. The group-based truncation leads to the melting of the DNA double helix. In contrast, the atom-based truncation maintains the helical structure intact. Similarly for the protein, the group-based truncation leads to an unfolding at pH 2 while the atom-based truncation produces stable trajectories at low and normal pH, in agreement with experiment. Artificial repulsion between charged residues associated with the group-based truncation is identified as the microscopic reason behind unfolding of the protein. Implications of different truncation schemes in reaction-field simulations of biomolecules are discussed.