András Baranyai

A systematic development of a polarizable potential of water

Based on extensive studies of existing potentials we propose a new molecular model for water. The new model is rigid and contains three Gaussian charges. Contrary to other models, all charges take part in the polarization of the molecule. They are connected by harmonic springs to their gas-phase positions: the negative charge to a prescribed point on the main axis of the molecule; the positive charges to the hydrogens. The mechanical equilibrium between the electrostatic forces and the spring forces determines the polarization of the molecule which is established by iteration at every timestep. The model gives excellent estimates for ambient liquid properties and reasonably good results from high-pressure solids to gas-phase clusters. We present a detailed description of the development of this model and a large number of calculated properties compared to the estimates of the nonpolarizable TIP4P∕2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)], the polarizable GCPM [P. Paricaud, M. Predota, A. A. Chialvo, and P. T. Cummings, J. Chem. Phys. 122, 244511 (2005)], and our earlier BKd3 model [P. T. Kiss and A. Baranyai, J. Chem. Phys. 137, 084506 (2012)]. The best overall performance is shown by the new model.

New algorithm for constrained molecular-dynamics simulation of liquid benzene and naphthalene

Two algorithms are in common use for performing molecular-dynamics simulations of fluids composed of molecules with holonomic bond constraints. The SHAKE algorithm is widely used but is conceptually complex. The algorithm of Edberg et al. (EEM), is simple in concept and structure except in regard to the manner in which numerical drift is handled. We describe a simple way of modifying the EEM method so that numerical drift in the holonomic constraints can be treated in a simple non-iterative fashion. The method is applied to fluids composed of planar molecules: liquid benzene and naphthalene. The results for the former system are in agreement with previous simulations using site-site potential models.

Clusters of classical water models

The properties of clusters can be used as tests of models constructed for molecular simulation of water. We searched for configurations with minimal energies for a small number of molecules. We identified topologically different structures close to the absolute energy minimum of the system by calculating overlap integrals and enumerating hydrogen bonds. Starting from the dimer, we found increasing number of topologically different, low-energy arrangements for the trimer(3), the tetramer(6), the pentamer(6), and the hexamer(9). We studied simple models with polarizable point dipole. These were the BSV model [J. Brodholt et al., Mol. Phys. 86, 149 (1995)], the DC model [L. X. Dang and T. M. Chang, J. Chem. Phys. 106, 8149 (1997)], and the GCP model [P. Paricaud et al., J. Chem. Phys. 122, 244511 (2005)]. As an alternative the SWM4-DP and the SWM4-NDP charge-on-spring models [G. Lamoureux et al., Chem. Phys. Lett. 418, 245 (2006)] were also investigated. To study the impact of polarizability restricted to the plane of the molecule we carried out calculations for the SPC-FQ and TIP4P-FQ models, too [S. W. Rick et al., J. Chem. Phys. 101, 6141 (1994)]. In addition to them, justified by their widespread use even for near critical or surface behavior calculations, we identified clusters for five nonpolarizable models of ambient water, SPC/E [H. J. C. Berendsen et al., J. Phys. Chem. 91, 6269 (1987)], TIP4P [W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983)], TIP4P-EW [H. W. Horn et al., J. Chem. Phys. 120, 9665 (2004)], and TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)]. The fifth was a five-site model named TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)]. To see the impact of the vibrations we studied the flexible SPC model. [K. Toukan and A. Rahman, Phys. Rev. B 31, 2643 (1985)]. We evaluated the results comparing them with experimental data and quantum chemical calculations. The position of the negative charge in the models plays a crucial role. In this respect models with SPC geometry provided structures different from the TIP4P-type potentials, including polarizable ones. The TIP4P variants form configurations similar to one another. Results for TIP4P-EW and for TIP4P/2005 were especially close to each other in every respect. This is also true for the BSV and the DC pair. The charge-on-spring models (SWM4-DP and SWM4-NDP) are also very similar to each other, despite the sign exchange of charges on the spring particle and the oxygen. The spherical polarization of water is crucial. Due to the planar polarization of the SPC-FQ and the TIP4P-FQ models, they prefer planar arrangements contrary to other polarizable models and quantum chemical calculations. The tetrahedral geometry of TIP5P stabilizes additional clusters with peculiar geometries and small O-O distances. Inclusion of vibrations causes only insignificant changes in the characteristic geometries but decreases the internal energy relative to its reference rigid version. Comparing with quantum mechanical calculations the GCP model provided the best overall results.

Steady state simulation of planar elongation flow by nonequilibrium molecular dynamics

We present a novel method for performing steady state nonequilibrium molecular dynamics simulation of planar elongation flow based on the studies of Kraynik and Reinelt [Int. J. Multiphase Flow 18, 1045 (1992)]. These authors identified the orientation of the unit cell which leads to periodic behavior of the square lattice with the minimum period. This way the exponential deformation of the system periodically returns to a state where replacing some of the original particles with their images the initial state boundaries are recovered. We adopted their theoretical results to nonequilibrium molecular dynamics simulations and performed representative calculations for simple fluids. The new method solves the long-standing problem of simulating planar elongation flow in the steady state.

A transferable classical potential for the water molecule.

We developed a new model for the water molecule which contains only three Gaussian charges. Using the gas-phase geometry the dipole moment of the molecule matches, the quadrupole moment closely approximates the experimental values. The negative charge is connected by a harmonic spring to its gas-phase position. The polarized state is identified by the equality of the intermolecular electrostatic force and the spring force acting on the negative charge. In each timestep the instantaneous position of the massless negative charge is determined by iteration. Using the technique of Ewald summation, we derived expressions for the potential energy, the forces, and the pressure for Gaussian charges. The only properties to be fitted are the half-width values of the Gaussian charge distributions and the parameters of the nonelectrostatic repulsion-attraction potential. We determined the properties of gas-phase clusters up to six molecules, the internal energy and density of ambient water and hexagonal ice. We calculated the equilibrium density of ice VII as a function of pressure. As an additional test, we calculated the pair-correlation function, the isotherm compressibility, the heat capacity, and the self-diffusion coefficients for ambient water. As far as we know, this is the first classical model of water which is able to estimate both ends of the phase diagram, the high pressure ice VII, and the gas clusters of water with excellent accuracy.

A new polarizable force field for alkali and halide ions

We developed transferable potentials for alkali and halide ions which are consistent with our recent model of water [P. T. Kiss and A. Baranyai, J. Chem. Phys. 138, 204507 (2013)]. Following the approach used for the water potential, we applied Gaussian charge distributions, exponential repulsion, and r(-6) attraction. One of the two charges of the ions is fixed to the center of the particle, while the other is connected to this charge by a harmonic spring to express polarization. Polarizability is taken from quantum chemical calculations. The repulsion between different species is expressed by the combining rule of Kong [J. Chem. Phys. 59, 2464 (1972)]. Our primary target was the hydration free energy of ions which is correct within the error of calculations. We calculated water-ion clusters up to 6 water molecules, and, as a crosscheck, we determined the density and internal energy of alkali-halide crystals at ambient conditions with acceptable accuracy. The structure of hydrated ions was also discussed.

Statistical geometry of molten alkali halides

The generalized Dirichlet–Voronoi polyhedra were used to characterize the adjacent neighborhood of the ions in computer simulated configuration set of molten RbCl and LiI as examples of binary Coulomb liquids with like and unlike size ions, respectively. The statistical distribution of various features of the DV polyhedra has been studied with the following main results: (i) radial and angular distributions of adjacent neighbors indicate traces of face‐centered‐cubic symmetry in molten RbCl and wurtzite‐like structure in LiI; (ii) the radial internal point distribution function on the DV polyhedra reveal more compact coordination spheres for Cl−, Rb+, and I− ions than for Li+ ions in the corresponding melts; (iii) the volume fluctuation of the DV polyhedra makes it possible to calculate isothermal compressibility coefficients even from (N,V,T) simulation data; (iv) identification of vacancies by seeking for the most protruding vertices of the DV polyhedra provides information on vacancy–vacancy pair correlation functions.

Grand canonical Monte Carlo simulation of liquid argon

A grand canonical Monte Carlo procedure with fixed values of the chemical potential μ, volume V, and temperature T, is described which is suitable to simulate simple fluids with only a minor increase in computer time in comparison with canonical (N,V,T) simulations and considerably faster than (N,p,T) ones. The method is rapidly convergent for rather dense systems with a reduced density of about ρσ3=0.88. The rapid convergence is attained by decreasing the vain attempts in the regime when new particles are added. The chance to find a place for an additional particle is increased by locating the cavities suitable to house a particle with the aid of the Dirichlet–Voronoi polyhedra. As an example, liquid argon is simulated with Lennard‐Jones potentials at T=86.3 K and μ=−73.4 J/mol. The simulated density has been found to be 1.468 g/cm3 which is to be compared with the experimental value of 1.425 g/cm3. The same density was obtained by starting the procedure with both 216 and 250 particles in the simulation ...

Monte Carlo simulation of the complete set of molten alkali halides

The results of Monte Carlo simulations of the complete family of molten alkali halides with a consistent set of inter-ionic potentials are reported. It is found that the structures of the melts with a common cation are isomorphous, the density scaling approximately with the sum of the ionic radii. Within groups with a common anion the local coordination varies from tetrahedral for Li+ salts to octahedral for Cs+ salts. General agreement is found with the available experimental results and with the results of other simulations using different potentials.

NONEQUILIBRIUM MOLECULAR DYNAMICS STUDY OF SHEAR AND SHEAR-FREE FLOWS IN SIMPLE FLUIDS

Nonequilibrium molecular dynamics simulations have been performed in order to compare the characteristics of planar Couette, planar elongation, uniaxial stretching, and biaxial stretching flows in simple fluids at different strain rates. After deriving the periodic boundary conditions for general flow fields and introducing some methodological improvements for elongation flow calculations we simulated the combination of shear and shear‐free flows as well. We found that even at high strain rates where simple fluids exhibit strong non‐Newtonian behavior (shear‐thinning) it is a reasonable approximation to consider the two planar flows to be rotationally equivalent. This is because in planar Couette flow the in‐plane normal stress difference of simple fluids is approximately zero even far from equilibrium. Similarly to planar Couette flow, the trace of the pressure tensor and the internal energy vary approximately as function of the 3/2 power of the strain rate in shear free flows. However, the individual di...