Electric-field induced entropic effects in liquid water

In addition to quantum mechanics, matter’s behavior is ultimately ruled by the subtle interplay between local electric fields.¹ Hence, it is not surprising that the application of external electric fields can be used to alter the thermodynamics of complex molecular systems in much the same way as temperature and pressure have been traditionally exploited.^2–5 This way, a widespread series of chemical and physical phenomena may be triggered by opportunely tuned oriented fields.

At the present time, many laboratories are able to produce controllable static fields on the order of 1–10 V nm⁻¹ at the proximity of high-capability field emitter tips in Scanning Tunneling Microscopy (STM) and Atomic Force Microscopy (AFM) apparatus,^6–8 rendering evident the possibility that tailored static electric fields will routinely be employed as smart catalytic tools in the imminent future.⁹ Indeed, it has been proven that electrostatic potential gradients can selectively promote and catalyze Diels–Alder¹⁰ and other industrially relevant¹¹ chemical reactions. These findings represent some of the first experimental evidences, indicating that electric fields with magnitudes on the order of $∼1−10$ V nm⁻¹ can drive and finely control specific chemical transformations. On the other hand, detailed conclusions on the effects produced by the application of intense electric fields were achieved through traditional quantum-based computational approaches on single-molecule systems^12,13 and—by means of ab initio and advanced metadynamics techniques—on condensed phases of matter.^14,15 Field strengths on the order of $∼5$ V nm⁻¹ have been employed in order to reproduce in silico Miller¹⁶ and Miller-like¹⁷ experiments and to highlight some crucial aspects of the role played by electrostatics on general chemical reactivity.^9,14,18 Those field intensities couple with local electric fields naturally present in condensed matter, whose intensity falls exactly in the same range of magnitude.¹⁹ Indeed, electric fields on the order of 10 V nm⁻¹ are generated by the solvent in aqueous systems.²⁰ In addition, as proven by recent investigations carried out via two-dimensional infrared (2D-IR) spectroscopy, field strengths of 3 V nm⁻¹ are ubiquitous in neat liquid water and are responsible for the solvation state of protons.²¹ Pioneering Car–Parrinello molecular dynamics simulations²² and more recent and sophisticated ab initio molecular dynamics (AIMD) investigations,^23,24 also including subtle Nuclear Quantum Effects (NQEs) via path-integral AIMD (PI-AIMD),²⁵ have partly elucidated the role played by static electric fields when applied on bulk liquid water, shedding some light on its spectroscopic, molecular, and proton conductivity response upon exposure to electric fields. In vibrational Stark spectroscopy, frequency shifts induced by similarly intense electric fields in selected vibrational modes are conveniently exploited as local probes.²⁶ Furthermore, proton transfer reactions triggered by strong electric fields spontaneously present at the interface between water and different solid materials have recently been characterized.^27–29 This way, electrostatic potential gradients are capable not only of reorganizing molecular dipoles in such a way to freeze liquids such as ammonia³⁰ but also of opening pathways in complex reaction networks, shaping the chemistry of aqueous systems.^16,17,31

Notwithstanding some thermodynamic investigations reported on the effects generated by electric fields on the solid–liquid equilibrium,³² electric-field-induced entropy modifications associated with the water hydrogen-bond (H-bond) reorganization have only very recently been investigated at electrified interfaces by means of multiscale simulations.³³ In particular, Estejab et al., by means of classical MD and AIMD simulations, studied the solvation thermodynamics of a series of surface species on the assumption that local electric fields could be able to reorient hydration water molecules.³³ Although a key finding—concerning the fact that the entropy of solvation depends on both the strength and direction of the electric field—was reported,³³ this study was affected by the (inextricable) competition between the field and surface effects at water/electrode interfaces. In addition, many computational efforts were traditionally devoted to the well-known hydrophobic effect,³⁴ whose molecular-scale origins are modified by local charges (and hence fields).^35,36 However, to the best of our knowledge, neither qualitative nor quantitative assessments on the role played by electric fields in shaping the local entropy of bulk water have been reported so far.

The aim of the current study is to exploit force-field MD and AIMD simulations of bulk liquid water under the effect of static and homogeneous electric fields to unveil field-induced changes on the structure and entropy of the system. Moreover, here, we also track the outcomes stemming from local entropy descriptors capable of returning a spatially resolved quantification of the entropy. In particular, we use an MD-based analysis (3D-2PT)^37,38 that spatially resolves the entropy per water molecule and the local number density of bulk water under an electric field and their changes relative to the zero-field regime. As for these local entropy indicators, we investigate the bulk water environment of two different scenarios at diverse field intensities: (i) the ideal system of a water molecule perfectly aligned along the field (named hereinafter as fixed) and (ii) the more realistic system of a water molecule freely rotating around a locked point in the bulk liquid (named hereinafter as pinned). While the ideal system allows us to characterize the complex water environment along specific directions of the perfectly aligned molecule (e.g., bond vectors and dipole moment), the more realistic system is needed to describe the response of the central water molecule to the externally imposed electric field, i.e., the combined effect of field alignment and local structure preferences.

The liquid water sample contained 128 H₂O molecules (i.e., 384 atoms) arranged in a cubic cell with 15.82 Å edges so as to reproduce a density of 0.97 g cm⁻³. As usual, structures were replicated in space by employing periodic boundary conditions. The intensity of the electric field was gradually increased with a step increment of 0.5 V nm⁻¹ from zero up to a maximum of 2.0 V nm⁻¹. Albeit the latter value is slightly beneath the computational^22,25 and experimental^7,8,45 water dissociation threshold, from our prolonged AIMD trajectories, a non-zero probability of finding field-induced ionizations at 2.0 V nm⁻¹ was recorded, a result in accordance with recent simulations.²⁵ Thus, to make a one-to-one comparison with the entropy estimates stemming from classical force-field simulations, we decided to consider for the analysis of the entropy only the trajectories where proton transfers do not occur. On the other hand, for structural analyses, also data stemming from simulation at the highest field intensity (2.0 V nm⁻¹) were post-processed. In the zero-field case, we performed dynamics of 50 ps, whereas for each other value of the field intensity, we ran dynamics of almost 250 ps with a time step of 0.5 fs.

Wavefunctions of the atomic species have been expanded in the TZVP basis set with Goedecker–Teter–Hutter (GTH) pseudopotentials using the Gaussian plane waves (GPW) method.⁴⁶ A plane wave cutoff of 400 Ry has been imposed. Exchange and correlation (XC) effects were treated with the gradient-corrected Becke–Lee–Yang–Parr (BLYP)^47,48 density functional. Moreover, in order to take into account dispersion interactions, we employed the dispersion-corrected version of BLYP [i.e., BLYP+D3(BJ)].^49,50 A nominal temperature slightly higher than the standard one has been simulated in order to better reproduce the water structure (i.e., T = 350 K). The dynamics of ions was simulated classically within a constant number, volume, and temperature (NVT) ensemble using the Verlet algorithm, whereas the canonical sampling has been executed by employing a canonical-sampling-through-velocity-rescaling thermostat (CSVR)⁵¹ set with a time constant equal to 10 fs.

In the zero-field regime, comparisons of some structural properties of water were performed by means of the more recent and accurate non-empirical strongly constrained and appropriately normed (SCAN) functional⁵² (plus the rVV10 dispersion correction⁵³). In this case, a cubic box containing 97 H₂O molecules (i.e., 291 atoms) arranged in a cubic cell with edges equal to 14.30 Å was simulated for $∼50$ ps. Wavefunctions have been expanded in the TZVP basis set with GTH pseudopotentials specifically optimized for SCAN, while a large plane wave cutoff of 1200 Ry has been imposed after several tests. Moreover, optimized parameters for SCAN were used for mimicking dispersion interactions within the rVV10 method (i.e., b = 11.40 and c = 0.93). A temperature of 300 K was imposed in the SCAN+rVV10 case via a CSVR thermostat⁵¹ set with a time constant equal to 50 fs.

Classical molecular dynamics (MD) simulations of liquid water under static and homogeneous electric fields were carried out with the Gromacs 2018.1 software package.⁵⁴ The TIP4P/2005 model⁵⁵ was used to describe water–water interactions. Although more sophisticated potentials have been developed in recent decades, the reliability of TIP4P/2005 has been demonstrated by recent interesting results, also concerning the application of external electric fields.^5,56

To analyze the bulk entropy of TIP4P/2005 water for comparison with AIMD, we performed one set of simulations with 1024 water molecules in the canonical ensemble with a fixed-volume cubic box with 31.62 Å edges. A second set of simulations for the same system was performed in the isobaric–isothermal ensemble at a constant pressure of 1 bar to study the impact of field-induced density changes on the water entropy, which were not captured in constant-volume AIMD simulations. To spatially resolve local water entropies in the 3D environment of a central water molecule, we used third and fourth sets of simulations with 2147 water molecules in a cubic box with initial edges of 42 Å, which were equilibrated at a pressure of 1 bar for each applied electric field strength prior to sampling at constant volume.

As described in the Introduction, these last two sets of simulations described two different scenarios, which we refer to as fixed and pinned. In the fixed set of simulations, the dipole moment of a central water molecule was aligned with the Z axis along which electric fields were applied. Cartesian coordinates of this water molecule were kept fixed using harmonic position restraints with a large force constant [500 kJ/(mol Ų)]. In the pinned set of simulations, the harmonic position restraint was only applied to the oxygen atom of the central water molecule, thus allowing the latter to rotate freely. These two scenarios allow for different high-resolution analyses of water properties in the environment of the restrained water molecule, as described in detail in Subsection II B.

Each of the four sets of simulations with the TIP4P/2005 model was simulated in the zero-field regime and in the presence of electric fields E of varying magnitude between 0.5 and 2.0 V nm⁻¹. The fields were implemented using an external force along the Z axis that acts on all atoms proportional to their partial atomic charge. In all simulations, long-range electrostatics were treated with the particle mesh Ewald algorithm on a 1.2 Å spatial grid with a fourth order interpolation and tin-foil boundary conditions.⁵⁷ Short-range electrostatics and Lennard-Jones forces were shifted to zero between 9.0 and 10.0 Å. Long-range dispersion corrections were applied for the energy and pressure. Neighbor lists were updated every 10 fs (equilibration) or 8 fs (production) with a 13.0 Å distance cutoff. The time step for integration was set to 1 fs (for both equilibration and production). The SETTLE algorithm⁵⁸ was used to constrain the intramolecular degrees of freedom of water molecules.

The bulk water systems with 1024 molecules were simulated with a time step of 0.5 fs, and coordinates were stored every 1 fs to allow for direct comparisons with AIMD simulations using identical analysis protocols. The first set of simulations at constant volume was equilibrated for each applied field strength for 100 ps using the Berendsen weak-coupling thermostat with a reference temperature of 300 K and a time constant of 1.0 ps.⁵⁹ This was followed by 250 ps of production simulations with the Nose–Hoover thermostat and a time constant of 2.0 ps.⁶⁰ The second set of simulations of these systems at constant pressure used the same protocol with added barostats at a reference pressure of 1 bar. The Berendsen weak-coupling barostat⁵⁹ and the Parrinello–Rahman barostat⁶¹ were used during equilibration and production, respectively, with a time constant of 1.0 ps.

Simulations of fixed and pinned systems with 2147 water molecules were equilibrated for each electric field strength at constant pressure with the same basic protocol as described above. However, production simulations in the isothermal–isobaric ensemble were extended to 80 ns in this case to sample 800 independent microstates in the NPT ensemble, i.e., coordinates and velocities every 100 ps. Each of these microstates was then used to initiate microcanonical (NVE) simulations of 100 ps length each with 1 fs time steps, which were used for the subsequent analysis. During these 800 production runs per system, coordinates and velocities were saved every 8 fs. Compared to the previous sets of simulations that were directly compared to AIMD simulations, we increased the integration time step to 1 fs for the equilibration of the fixed and pinned systems in consideration of the lack of fast intramolecular vibrations in the rigid TIP4P/2005 model. In NPT simulations, we increased the time step further to 2 fs to increase the sampling speed. For the following simulations in the microcanonical ensemble used for the analysis, we decreased the time step again to 1 fs to ensure strict energy conservation and accurate dynamics in accord with previously established protocols for entropy estimations using 3D-2PT.^37,38,62

The main feature of the 2PT framework is that the VDoS of translational and rotational degrees of freedom obtained from MD simulations are decomposed into separate contributions that are treated by simple analytical model systems. Specifically, the translational degrees of freedom are described as a superposition of a hard sphere fluid and a set of harmonic oscillators that cover the range of all frequencies sampled in simulation. By construction, the hard sphere fluid describes all diffusive components of the VDoS that are sampled at zero frequency.

A similar approach is used for rotational degrees of freedom of simulated molecules. Here, the hard sphere fluid is replaced by the rigid rotor model to describe zero-frequency components of the rotational VDoS. Intramolecular vibrations of water molecules lack diffusive processes at zero frequency and are thus simply treated as a spectrum of harmonic oscillators.

Standard expressions in statistical mechanics for thermodynamic properties of hard sphere fluids, rigid rotors, and harmonic oscillators are then used to describe thermodynamic properties of the simulated liquid based on the decomposition of its VDoS described above. While the accuracy of the approach is intrinsically limited by simplifications introduced through analytical models, absolute molar entropies obtained with 2PT from molecular simulations have been shown to correlate well with experimental data.⁶⁴ In addition, applications within the 3D-2PT framework have shown that changes in entropy, e.g., in response to the presence of a solute, benefit from the cancellation of systematic errors and can be obtained within the accuracy limits of the applied simulation model, e.g., empirical force fields used in simulation.^37,62

Average molar entropies without spatial resolution were obtained for AIMD simulations and TIP4P/2005 simulations of smaller systems with 1024 water molecules. To compute velocity time correlation functions from these simulations, atomic velocities were first reconstructed using finite differences and atomic coordinates saved at 1 fs intervals. For each molecule, atomic velocities were then converted into velocities of the molecular center of mass (translations), angular velocities around the rotational axes with distinct moments of inertia (rotations), and, in the case of AIMD simulations, intramolecular vibrations. Time correlation functions of distinct velocity components were then computed for a maximum correlation time of 1.0 ps, averaged for all water molecules, Fourier transformed into the corresponding VDoS, and used to estimate the entropy with the 2PT formalism described briefly above and in detail elsewhere.^37,63,64

For the fixed and pinned systems simulated with the TIP4P/2005 water model, we used atomic velocities evaluated during simulation and stored in the trajectory together with atomic coordinates at time intervals of 8 fs. To analyze molecular entropies with spatial resolution, we defined a three-dimensional (3D) grid around the restrained central water molecule. The grid contained 80 × 80 × 80 volume elements (voxels) with a grid constant of 0.5 Å. The orientation of the analysis grid was aligned with the axes of the coordinate system and was kept unchanged during the analyses for both the fixed and pinned systems. The analysis of the translational and rotational VDoS described above (the TIP4P/2005 water model lacks intramolecular vibrations) was then performed for each voxel of the analysis grid using a standard protocol adopted from previous 3D-2PT studies.^37,62 Water molecules analyzed for a given time frame of the trajectory were assigned to a voxel based on their center of mass position. Auto-correlation functions of water translational and rotational velocities were computed for correlation times of up to 1.6 ps and Fourier transformed into the corresponding VDoS. The resulting local VDoS for the translational and rotational degrees of freedom for each voxel were averaged over all simulation trajectories for each system and used to estimate the local water entropy with the 2PT formalism.^37,63,64 In the pinned set, at zero field, we expect isotropic physical properties in the surrounding of the freely rotating central water molecule, while at non-zero field, any sign of anisotropy with the increasing field intensity is due to the alignment of the molecules of the system (including the central water molecule). This is not valid for the fixed set: due to the set orientation of the central water molecule, the results are anisotropic even in the zero-field regime. Notably, since both the analysis grid and the central water molecule in the fixed system do not change their orientation with respect to the simulation box and, therefore, do not change their reciprocal orientation either, an advantage of the fixed system is that we can characterize properties of water in the environment along specific directions (i.e., analysis along the O–H bond vectors or along the water dipole moment). On the other hand, in the pinned set of simulations, the rotation of the central water molecule in the fixed analysis grid results in a spherically symmetric environment at zero field strength and cylindrical symmetry along the Z axis with the applied electric fields, as schematically shown in Fig. 1.

Before determining a spatially resolved molecular picture of structural and thermodynamic properties of the bulk water environment, it is worth to preliminarily test the accuracy—in the absence of the field—of the commonly employed TIP4P/2005⁵⁵ water parameterization with respect to two different dispersion-corrected Density Functional Theory (DFT) exchange–correlation (XC) functionals belonging to diverse classes, i.e., BLYP+D3^47–50 and SCAN+rVV10.^52,53 Whereas the former represents by far one of the most employed Generalized Gradient Approximation (GGA) density functionals in ab initio MD (AIMD) simulations of bulk liquid water, the latter is considered the forefront of meta-GGA DFT functionals for reproducing the behavior of aqueous systems. In particular, the SCAN functional corrected by the rVV10 dispersion parameterization is capable of overcoming some of the serious drawbacks characterizing the previous generation of DFT XC functionals belonging to the GGA class. As shown in Fig. 2(a), SCAN+rVV10 (blue dotted curve) is capable of faithfully reproducing the oxygen–oxygen radial distribution function g_OO(r) emerging from neutron scattering experiments⁶⁶ (profile on top of the gray shaded area) at ambient conditions. Interestingly, while Grimme’s D3^49,50 dispersion corrections are not capable of adequately softening intermolecular H₂O interactions, which result to be sizably overestimated by the BLYP functional (red dashed curve) even at a temperature higher than the ambient one, the TIP4P/2005 force-field reproduces spatial correlations between the water molecules in bulk liquid water at room temperature remarkably well (black solid curve). Thus, the hydrogen-bond (H-bond) network organization emerging from g_OO(r) appears to be adequately described by the computationally efficient and rigid TIP4P/2005 water model, similarly to cutting-edge and demanding dispersion-corrected first-principles meta-GGA functionals.

The over-structuring induced by the employment of GGA XC functionals visibly emerges from the evaluation of structural order parameters, such as the well-known orientational tetrahedral order parameter q.⁶⁷ As shown in Fig. 2(b), indeed, the distribution P(q) in the absence of the external electric field shows that the population of tetrahedrally ordered water molecules is overestimated in the BLYP+D3 framework with respect to the TIP4P/2005 force-field case. Incidentally, the orientational order observed with this latter classical potential is adherent to results stemming from the predictive meta-GGA SCAN+rVV10 functional [Fig. 2(b)]. Albeit meta-GGA DFT functionals are currently computationally too demanding to allow for the generation of prolonged trajectories in the presence of an electric field, the adherence of the structural data stemming from the rigid TIP4P/2005 water model gives us confidence with regard to the quality of its description of structure-related thermodynamic properties, such as the entropy. On the other hand, before considering local entropy density estimates for the TIP4P/2005 model only, we compared total entropy changes induced by the application of an external electric field in simulations with the TIP4P/2005 force-field and BLYP+D3, one of the most employed gradient- and dispersion-corrected DFT functionals.

Since the structural ordering emerging from the analysis of the data shown in Fig. 2 is tightly related to the estimate of the entropy, it is not surprising that the total entropy estimated for BLYP+D3 water is lower than that determined for TIP4P/2005 (i.e., 44.50 ± 0.26 J mol⁻¹ K⁻¹ vis-à-vis 54.17 ± 0.06 J mol⁻¹ K⁻¹, respectively), as shown in Fig. 3 (black dashed lines). Here, statistical uncertainties of the computed molar entropies are reported as the standard error of the mean obtained by averaging over non-overlapping 10 ps segments of the trajectory, as shown in Fig. 3 (the first 10 ps segment in each trajectory is treated as additional equilibration time and omitted from the analysis). Due to the already discussed over-structuring of water in GGA XC functionals,⁶⁸ the recorded value is significantly lower than the experimentally measured one at room temperature (i.e., 69.95 J mol⁻¹ K⁻¹).⁶⁹ In fact, our AIMD results for zero field are comparable to an earlier 2PT calculation reported by Galli’s group for AIMD simulations of water with the Perdew–Burke–Ernzerhof (PBE) functional (i.e., 42.74 J mol⁻¹ K⁻¹).⁷⁰ However, we note that a study by Kühne and co-workers observed a higher value for the entropy of PBE water (i.e., 51.32 J mol⁻¹ K⁻¹).⁶⁵ While our trajectory for the meta-GGA SCAN+rVV10 functional is relatively short (i.e., $∼50$ ps), we obtained an average water entropy of 52.17 ± 0.97 J mol⁻¹ K⁻¹, which is still significantly lower than the experimental value but only marginally lower than the entropy obtained from classical simulations with the TIP4P/2005 model. 2PT provides quantitative estimates of the molar entropy of liquids, which can contain systematic errors in addition to errors in the simulation model. For most water models, 2PT tends to underestimate the absolute molar entropy of water⁶⁵ as observed here for both the TIP4P/2005 model and SCAN+rVV10 simulations in comparison to the experimental reference value. Furthermore, by construction, the formalism is sensitive to the diffusion coefficient and related finite-size effects.⁷¹ However, such systematic errors cancel out for entropy differences,⁶² for example, between simulations of water with the same model but in the presence of distinct external fields.

When the external electric field is switched on at 0.5 V nm⁻¹, the initial decrement of the total entropy appears to be very modest, especially in the TIP4P/2005 series of simulations [Fig. 3(a)]. When electrons are treated as quantum entities in BLYP+D3 simulations, the observed field-induced entropy changes are also moderate, accounting to a reduction—in the first-principles case—of less than 2 J mol⁻¹ K⁻¹ [Fig. 3(b)]. However, when we compare Figs. 3(a) and 3(b), it is apparent that the entropy of TIP4P/2005 water is significantly less sensitive to polarization than BLYP+D3 water. In fact, whereas a maximum entropy decrement of $∼3%$ is observed at the highest field regime in the classical simulations (i.e., from 54.17 ± 0.06 to 52.33 ± 0.08 J mol⁻¹ K⁻¹), a reduction of $∼14%$ (i.e., from 44.50 ± 0.26 to 38.49 ± 0.30 J mol⁻¹ K⁻¹) emerges for BLYP+D3 AIMD simulations at 1.5 V nm⁻¹. This is likely due to the fact that the system treated at the BLYP+D3 level is able to respond to the external electric field up to higher field strengths owing to an explicit treatment of the electronic polarization. In fact, a more evident conversion toward more ordered molecular configurations, associated with a smaller total entropy, is recorded in the BLYP+D3 case, as displayed in Fig. 3(b). Beyond this threshold, AIMD simulations executed at the BLYP+D3 level have shown only a marginal further decrease of the total entropy. In addition, for field intensities beyond 1.5 V nm⁻¹, not negligible field-induced molecular dissociations and consequent proton transfers are recorded,²⁵ rendering the entropy analysis not comparable to the classical—non-polarizable and non-dissociable—force-field case. In further contrast to the rigid TIP4P/2005 water model, vibrational degrees of freedom are taken into account in the AIMD case, in addition to translational and rotational contributions to the total entropy. In a nutshell, the spectral features emerging upon field exposure in the VDoS reveal an increased formation of H-bonds, in agreement with infrared (IR) and Raman spectroscopic analyses conducted on equivalently treated ab initio water.²³ 

AIMD and TIP4P/2005 simulations discussed above were performed at constant volume and thus ignore field-induced effects in the structure of water that affect its density. To estimate the importance of density effects for the reported molecular entropy, we performed a second set of simulations for the TIP4P/2005 water model at constant pressure as described in Sec. II. At zero field, TIP4P/2005 water simulation at constant pressure increased its density to 1.00 g cm⁻³, which reduced its entropy to 51.91 ± 0.16 J mol⁻¹ K⁻¹ relative to 54.17 ± 0.06 J mol⁻¹ K⁻¹ observed at the constant volume with a density of 0.97 g cm⁻³. With the increasing field strength, the density increased gradually to 1.04 g cm⁻³ for the highest field strength of 1.5 V nm⁻¹. However, the response of the molecular entropy remained very similar to constant volume simulations. For the largest field strength of 1.5 V nm⁻¹, the molar entropy decreased relative to the zero-field regime by $∼2%$ to 50.69 ± 0.11 J mol⁻¹ K⁻¹ compared to an $∼3%$ decrease observed in constant volume simulations. Hence, compared to the $∼14%$ decrease of BLYP+D3 water, TIP4P/2005 water simulations remained equally insensitive to the field strength under constant volume and constant pressure conditions.

In light of relatively recent path integral AIMD (PI-AIMD) simulations executed under the effect of externally applied electric fields,²⁵ it is tempting to consider the impact that the inclusion of Nuclear Quantum Effects (NQEs) would have on the entropy estimate of liquid water. However, the 2PT formalism here adopted is based on the evaluation of velocity time correlation functions, which is, by definition, a dynamical property that cannot be reliably measured from path integral simulations.⁶¹ 

On the other hand, some effort has been done in the literature to determine the impact of NQEs on some structural properties of liquid water. It is well-known, indeed, that zero-point motion and proton tunneling lead to faster reorientational dynamics and to a slight de-structuring of radial distribution functions,⁷² presumably producing larger entropy estimates. However, it is worth stressing the fact that the inclusion of NQEs in classical force-fields molecular dynamics simulations and in the first-principles ones may generate entropic contributions of completely different physical nature. In fact, PI-AIMD trajectories simulating the quantum nature of nuclei and electrons exhibit a non-negligible amount of autoprotolysis events even in the absence of external electric fields.^25,73 Such a circumstance, which cannot be caught and reproduced within non-dissociable water force-fields, prevents reliable one-to-one comparisons of the genuine impact of NQEs on the entropy of classical and ab initio simulations of liquid water.

Although it has been shown that the application of increasingly stronger fields produces a progressive reduction of the total entropy of water, a corresponding order-maker action is not clearly reflected by field-induced changes of the tetrahedral order parameter of water. As shown in Fig. 4(a), the tetrahedral order actually decreases with the increasing field strength in water simulations with the TIP4P/2005 force field. For ab initio simulations with BLYP+D3, the tetrahedral order also initially decreases for a field strength of 1.0 V nm⁻¹, but then increases when the strength of the polarizing field is raised to 2.0 V nm⁻¹. Hence, the decrease in water entropy with the increasing field strength observed in Fig. 3 is not associated with an increase in tetrahedral order. Instead, the decrease of the water entropy is primarily associated with the increase in the orientational polarization of the water molecules. This can be observed in Fig. 4(b), where we show the distributions of the instantaneous angles θ between molecular dipoles of water with the Z axis of the coordinate system. Small values of the angle θ correspond to molecular dipole moments that are aligned with the electric fields applied along the Z axis. Albeit in the absence of the field all dipole orientations are equivalent, the imposition of the external electrostatic potential gradient induces a prompt reorganization of the water dipole vectors. As has been shown for other computational models of water,⁷⁴ the gradual saturation observed at the strongest fields resembles the behavior of free dipoles described by the Langevin–Debye equation for simulations with TIP4P/2005 force field and ab initio BLYP+D3 simulations.⁷⁵ This observation indicates a weak resistance of molecular dipoles to polarization by very intense fields. Such a field-induced dipole locking leads, among other things, to a progressive reduction of the tetrahedral translational order parameter q_k,⁶⁷ as listed in Table I.

Interestingly, the water over-structuring induced by the adoption of the BLYP+D3 functional leads to a more complete and concerted reorientation of the dipoles with respect to the TIP4P/2005 case at the same field strengths, in accordance with the results of Fig. 3 and as displayed in Fig. 4(b). Finally, our finding that the orientational polarization of BLYP+D3 water is already nearly complete for 1.0 V nm⁻¹—as well as the respective translational order parameter—is consistent with the observation that also the entropy reduction reaches a substantial saturation above this threshold [Fig. 3(b)]. By contrast, in the TIP4P/2005 case, the orientational polarization of water molecules at 1.0 V nm⁻¹ is much less pronounced [Fig. 4(b)] and increases further with the increasing field strength—similarly to the reduction of the respective translational order parameter—while the total entropy decreases [Fig. 3(a)].

To analyze molecular entropies of TIP4P/2005 water as a function of the electric field, we used the 3D-2PT approach introduced by Persson et al.^37,38,63,64 to spatially resolve the local water entropy in the environment of a central water molecule. It is noteworthy that, differently from standard calculations of isotropically averaged distribution functions, the spatially resolved approach adopted here allows us to resolve structural anisotropies with a high resolution owing to the 0.5 Å grid constant of the 3D analysis grid. For every volume element of the analysis grid, we have calculated the entropy per water molecule for the two investigated systems: the fixed and the pinned system (see Sec. II and Fig. 1).

In Tables II and III, we report the bulk water entropy, obtained as an average for voxels separated by 10 Å or more from the central reference water molecule in the investigated field regimes for the fixed system (Table II) and the pinned one (Table III). The changes of the entropy per molecule relative to the zero-field regime are visualized in the histograms of Fig. 5. For the total molar entropy of bulk water in the zero-field regime, $StotE0$⁠, we obtain 55.77 ± 0.02 J mol⁻¹ K⁻¹ (fixed system) and 55.75 ± 0.02 J mol⁻¹ K⁻¹ (pinned system). These results deviate slightly from the molar entropy of water obtained for NVT simulation of 1024 water molecules at zero field (54.17 ± 0.02 J mol⁻¹ K⁻¹) used for comparison to AIMD simulations in Sec. III A, which is the result of finite size effects in the latter case and a lower frequency resolution in the analyzed VDoS.

Bulk water entropies in the fixed and pinned system are statistically identical at zero field because the position restraints applied to the central water molecule in both cases have no effect on an isotropic system at rest (no net translations or rotations). Minor differences between both systems are only observed for non-zero fields as the systems become polarized, even though only water molecules with a minimum distance of 10 Å to the central water molecule are considered for the analysis of bulk water properties in both systems. The position restraint applied to the oxygen atom of the central water molecule in the pinned system still has no impact on thermodynamics. However, under polarized conditions, the restrained orientation of the central water molecule in the fixed system modifies interactions between water molecules and thus has thermodynamic consequences, which are still detectable at distances of 10 Å and beyond.

In addition to the total entropy per bulk water molecule, we also report separate contributions from translational and rotational degrees of freedom in Tables II and III and in Fig. 5. These are computed separately from fluctuations of the center-of-mass velocity and angular velocities for rotations around the three distinct rotational axes of water molecules. The largest contribution $(∼84%)$ to the entropy of bulk water molecules in the absence of a polarizing field comes from the translational degrees of freedom, e.g., $StransE0=46.87±0.02$ J mol⁻¹ K⁻¹ in the fixed system. In contrast, the rotational degrees of freedom contribute $∼16%$ to the entropy per water molecule, e.g., $SrotE0=8.90$ J mol⁻¹ K⁻¹. With the increasing field strength, the total entropy per bulk water molecule decreases by $∼4%$ for the highest field strength of 2.0 V nm⁻¹. This total entropy decrease is primarily associated with the translational degrees of freedom due to their dominant contribution. However, as shown in Fig. 5, the relative change of the entropy contributions is significantly larger for the rotational degrees of freedom, e.g., $∼6%$ and $∼7%$ in the fixed and pinned systems, respectively. This corresponds to the expectation of preferentially oriented water molecules in the presence of electric fields. However, given that the rotational entropy in our decomposition is only a minor contributor to the total entropy of water, even a large relative change remains a minor contribution to the field-induced change of the total entropy. Instead, the total entropy change in response to the field remains dominated by the small relative change of the significantly larger translational entropy.

We note that the polarization-induced change in the total entropy per water molecule remains small compared to structural changes associated with phase transitions. ΔS_tot associated with the freezing of liquid water at 273 K is $∼22$ J mol⁻¹ K⁻¹. Thus, it is straightforward to conclude that electric fields below the ionization threshold are not capable to cause orientational polarization-induced electrofreezing in bulk water (no confinement) at room temperature. This finding is in agreement with—and extends the range of applicability of—recent experiments that demonstrate that the ice nucleation temperature of water is independent of the application of external electric fields of up to 0.1 V nm⁻¹.⁷⁶ This is in contrast with recent AIMD and PI-AIMD simulations that have shown the possibility to electrofreeze another important H-bonded system, such as ammonia, at temperatures and pressures at which it is in the liquid state.³⁰ Of course, determination of electric-field-induced effects on enthalpy and free energy, in addition to those on the entropy, would further clarify on this intriguing phenomenon of enormous importance. However, the estimate of the free energy and of the enthalpy in the presence of an external electric field is not straightforward. In fact, the value of the electric field strength to be included in the thermodynamic equations is, in general, quite different from the applied external field due to the intrinsic additional field generated by the polarized surface of the cavity, which depends on the specific geometry of the system and on the anisotropy of the dielectric constant.⁵⁶ 

In Fig. 6, the entropy density defined as the product of the water entropy per molecule and the water number density for each voxel is shown for the fixed and pinned systems in the XZ-plane. The color code describes changes in the local entropy density in the environment of the molecule placed at the origin of the analysis grid relative to bulk water at zero field. Due to the fully restrained position of the central water molecule in the fixed system, the results are already anisotropic in the absence of an externally applied electric field. Interactions with the restrained central water molecule force water molecules to adapt their position and orientation. In the pinned system, the environment of the central water molecule is isotropic at zero field, while anisotropy is gradually induced when an electric field is applied. In the fixed system, a local increase in the entropy density relative to bulk water with the increasing field strength is particularly visible for sites localized in the second coordination layer (4–5 Å) as shown in Fig. 6 (top) despite the overall decrease in the entropy of bulk water. In the pinned system (Fig. 6, bottom), whose thermodynamic properties are equivalent to a bulk water system, the applied fields restructure the isotropic entropy density into a cylindrically symmetric density around the Z axis. A decrease of the local entropy density is observed in the first coordination layer (3 Å) in the XY-plane and along the Z axis, which is followed by a local increase in the second coordination layer (4–5 Å). The trend is opposite at ±45° angles with the Z axis, where an increase in the entropy density is observed in the first coordination layer, while a decrease is found for the second coordination layer. Notably, the local water entropy densities in the environment of the central water molecule observed in Fig. 6 are strongly affected by the local water number density in both systems, which effectively describes the 3D water–water pair correlation function. Thus, we analyze in the following the spatially resolved water number density and its response to the applied fields in more detail.

We analyzed the average local water number density in the environment of the central water molecule on the same analysis grid as used for the analysis of local entropy densities above. In Fig. 7, we show the contour plots of the water number density, ρ(E), normalized to the bulk value at zero field, ρ_bulk(E0), in the environment of the fixed molecule placed at the origin of the analysis grid. Data are shown for different field intensities (i.e., 0.0, 1.0, and 2.0 V nm⁻¹) in the XZ-plane (top), the YZ-plane (middle), and the XY-plane (bottom), while—as stated above—the field is applied along the Z-direction. The average number density of bulk water at zero field has been obtained from voxels separated by 10 Å or more from the central water molecule. As seen earlier in Fig. 6, interactions with the restrained central water molecule already result in an anisotropic environment at zero field (left panels in Fig. 7) in accord with other reports of 3D pair correlation functions of water in a molecular reference frame.⁷⁷ Increased local densities of water are observed in the first coordination layer, specifically in H-bond acceptor sites in the XZ-plane (top panel) and the H-bond donor sites in the YZ-plane (middle panel) at 3 Å distance from the center. H-bond acceptor and donor sites are surrounded by a reduced number density nearby due to repulsions with water molecules in the first shell. Additional enhanced densities can be observed at ∼4 Å, which correspond to water molecules in the second coordination layer of the central water molecule. The anisotropic structure of water changes notably as a function of the applied electric field. In particular, Fig. 7 shows that locally increased water densities in the second hydration shell increase further with electric field strength in both the XZ- and YZ-planes. This is accompanied by a less structured first hydration shell that is best observed in modified peak intensities of oxygen–oxygen radial distribution functions (see the discussion of Fig. 8 below). The first peak in the radial distribution function decreases slightly with the increasing field strength, while the second peak increases significantly.

Our analysis of the local environment of the fixed water molecule thus shows that the application of a static electric field induces an enhanced layering of the first two hydration shells. The second solvation shell is affected by the field to a larger extent with respect to the first shell, in fairly good agreement with previous force-field molecular dynamics simulations by Shafiei et al.⁷⁴ and classical AIMD and path-integral AIMD simulations by some of us.²⁵ This observation can be explained by comparatively weaker correlations in the second hydration shell, which are easier to perturb by the applied electric field compared to the strong correlations in the first hydration shell that are the result of direct H-bonds with the central reference water molecule.

In Fig. 8, we analyzed changes in the local water structure around the central water molecule in the pinned system using contour plots of the local water number density in the XZ- and XY-planes. Due to the unhindered rotation of the central water molecule, its environment is isotropic in the absence of an electric field (left panels) and retains cylindrical symmetry around the Z axis once the field is applied (i.e., the XZ- and YZ-planes are symmetrically equivalent). The coordination layers, in particular the first one identifiable by a red ring at $∼3$ Å, can be easily located in field-free contour plot maps of Fig. 8 (left column). By applying an increasingly intense electric field, we can observe that the progressive dipolar alignment of the entire system—and in particular of the central water molecule—along the field direction results in manifest structural anisotropies in the XZ- and YZ-planes (Fig. 8, top). In other words, by increasing the field strength, the coordination shells show an enhanced layering toward specific directions due to the preferential orientation of all water molecules in the system and the resulting changes of intermolecular interactions.

In particular, in the XY-plane and along the Z axis, we observe a decrease in the local water number density in the first hydration shell (∼3 Å) and an increase in the second hydration shell (∼4–5 Å). However, the opposite trend is observed at an angle of 45° with respect to the Z axis, especially at the most intense field strength explored here (i.e., 2.0 V nm⁻¹), as shown in Fig. 8 (top right panel). This latter evidence suggests a field-induced strengthening of the H-bonds of the central reference molecule with the first-neighboring molecules along specific directions, a finding in agreement with the computational spectroscopic indication of the tendency of liquid water toward the formation of more “ice-like” structures upon field exposure at these intensities.²³ 

The increased local density of water in the second hydration shell at a distance of ∼4–5 Å from the center in the XY-plane and along the Z axis is accompanied by a depletion of water density in the first hydration shell. Depleted zones are observed along the Z axis at about ∼3–4 Å (Fig. 8, top right panel) as a consequence of the steric exclusion imposed by the nearby water molecules, which—under the field action—come closer to and establish stronger H-bonds with the central water molecule. On the other hand, a different response to the field can be visualized in the XY-plane (Fig. 8, bottom panel) where, at the highest field intensity, the peak of the first coordination layer essentially disappears, while that of the second coordination layer is increased. Such a finding strongly indicates that whereas first neighbors approach each other in the field direction (Z axis), they are pushed farther in the plane perpendicular to the field (XY-plane). This peculiar scenario is also caught and somehow exaggerated within the simpler fixed scenario, as shown in Fig. 7 (bottom panel, right).

In this work, we analyze the electric-field-induced changes in the total and local entropy of bulk liquid water by means of classical and first-principles simulations. In the absence of an external field, a comparison between the rigid force-field TIP4P/2005, GGA BLYP+D3, and meta-GGA SCAN+rVV10 ab initio functionals has shown that key structural properties of water appear to be adequately described by the computationally efficient TIP4P/2005 water model, which is in close agreement with the structure of water obtained from AIMD simulations with the cutting-edge SCAN+rVV10 functional. Related to this, the molecular entropy of bulk water that we estimated for the TIP4P/2005 water model using the 2PT formalism (54.17 J mol⁻¹ K⁻¹) is in close agreement with the corresponding result for AIMD simulations with the meta-GGA SCAN+rVV10 functional (52.17 J mol⁻¹ K⁻¹). That both estimates are significantly lower than the experimental reference value of 69.95 J mol⁻¹ K⁻¹ can be attributed to a systematic underestimation of the water entropy using the 2PT formalism, which is based on the vibrational density of states. However, both the TIP4P/2005 force field and the meta-GGA SCAN+rVV10 functional result in significantly higher estimates for the molar entropy of water than GGA DFT functionals, such as BLYP+D3, which results in significant overstructuring of the liquid even at elevated temperatures and yields a molar entropy of only 44.50 J mol⁻¹ K⁻¹. However, the molar entropy of TIP4P/2005 water seems to be significantly less sensitive to polarization-induced changes compared to BLYP+D3 water due to the inclusion of the electronic polarization in the latter.

In addition to these general results, we have presented a molecular-dynamics-based analysis (3D-2PT) that provides spatially resolved information on the local entropy and number density of water in the environment of a reference molecule. This method allowed us to generate 3D maps of changes in the water structure that are induced by electric fields. By means of this approach, a field-induced strengthening of the H-bonds of the reference molecule with the first-neighboring molecules along specific directions has been observed. This resulted to be in agreement with the vibrational Stark effect, indicating the tendency of liquid water evolving toward more “ice-like” structures upon field exposure at intensities of $∼1.0$ V nm⁻¹.²³ However, the largest field-induced entropy changes are too feeble (i.e., about 4% decrease for TIP4P/2005 water and 14% in the BLYP+D3 case) for realizing the so-called electrofreezing phenomenon. The fact that this finding holds up to strengths just beneath the molecular dissociation threshold (2.0 V nm⁻¹) extends the range of applicability of recent experiments, showing the independence of the water nucleation temperature from the application of external electric fields up to intensities of 0.1 V nm⁻¹.⁷⁶ 

Owing to the 3D-2PT formalism, we have shown that the field-induced local entropy density changes are strongly modulated by field-induced water number density changes, allowing us to depict a microscopic picture of the molecular entropic effects of the field. It turns out that the main effect produced by a static electric field on the water number density is that of inducing an enhanced layering, i.e., an over-coordination in correspondence of the second shell, though with anisotropic structural modifications. Under intense fields, water molecules undergo a progressive dipolar alignment and adopt preferential orientations parallel to the field lines, forming H-bond structures preferably along a direction that is diagonal with respect to the field axis.

G.C. acknowledges the financial support from ICSC, Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union, NextGenerationEU, PNRR, Missione 4 Componente 2 Investimento 1.4.

The authors have no conflicts to disclose.

Valeria Conti Nibali: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Sthitadhi Maiti: Data curation (supporting); Formal analysis (equal); Investigation (supporting); Methodology (supporting). Franz Saija: Conceptualization (equal); Investigation (equal); Supervision (equal). Matthias Heyden: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Writing – review & editing (equal). Giuseppe Cassone: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.