Simulations of water nano-confined between corrugated planes

TL;DR

This study uses computer simulations to explore how water confined between corrugated planes responds structurally and dynamically, revealing phase changes and inhomogeneous behavior influenced by wall oscillations.

Contribution

It introduces a novel periodic confining potential to simulate atomistic wall oscillations, analyzing their effects on water's phase behavior with both empirical and first-principles methods.

Findings

Melting temperature increases with commensurate lattice parameters.

Hexatic phase is replaced by trilayer crystalline phase at high modulation amplitudes.

Structural and dynamical responses are surprisingly insensitive to wall rugosity.

Abstract

Two-dimensionally nanoconfined water between ideal planar walls has been the subject of ample study, aiming at understanding the intrinsic response of water to confinement, avoiding the consideration of the chemistry of actual confining materials. In this work, we study the response of such nanoconfined water under a periodic confining potential by means of computer simulations, both using empirical potentials and from first-principles. We propose a periodic confining potential emulating the atomistic oscillation of the confining walls, which allows varying the lattice parameter and amplitude of the oscillation. We do it for a triangular lattice, with several values of the lattice parameter: one which is ideal for commensuration with layers of Ih ice, and other values that would correspond to more realistic substrates. For the former, an overall rise of the melting temperature is observed. The liquid maintains a bi-layer triangular structure, however, despite the fact that it is not favoured by the external periodicity. The first-principles liquid is significantly affected by the modulation in its layering and stacking even at relatively small amplitudes of the confinement modulation. Beyond some critical modulation amplitude the hexatic phase present in flat confinement is replaced by a trilayer crystalline phase unlike any of the phases encountered for flat confinement. For more realistic lattice parameters, the liquid does not display higher tendency to freeze, but it clearly shows inhomogeneous behaviour as the strength of the rugosity increases. In spite of this expected inhomogeneity, the structural and dynamical response of the liquid is surprisingly insensitive to the external modulation. Although the first-principles calculations give a more triangular liquid than the one observed with empirical potentials, both agree remarkably well for the main conclusions of the study.

Figures

Tables

Table 1. Table 1: Values of various simulation parameters as defined in the text, that depend on the chosen lattice parameter.

$a$ (Å)	${\rho }_{\text{LJ}}$ (Å-3)	${ϵ}_{12\text{-}6}$ (kcal/mol)	$l_{z_1}$ (Å)	$l_{z_2}$ (Å)	$L_x^{\text{FFMD}}$ (Å)	$L_y^{\text{FFMD}}$ (Å)	$L_x^{\text{AIMD}}$ (Å)	$L_y^{\text{AIMD}}$ (Å)
2.50	0.0905	0.0538	-0.446	-2.041	37.500	34.641	22.500	25.980
2.75	0.0680	0.0716	-0.456	-2.245	35.750	33.341	24.750	23.816
3.00	0.0524	0.0929	-0.464	-2.449	36.000	36.373	24.000	25.980
4.78	0.0129	0.3757	-0.487	-3.903	33.460	33.117	23.900	24.838

Table 2. Table 2: Lowest energy modulation E b subscript 𝐸 𝑏 E_{b} (see text) for a given value of a 𝑎 a and α 𝛼 \alpha .

$\alpha$	$E_b^{a=2.50\text{Å}}$ (K)	$E_b^{a=2.75\text{Å}}$ (K)	$E_b^{a=3.00\text{Å}}$ (K)	$E_b^{a=4.78\text{Å}}$ (K)
0.2	3.88	6.93	11.05	44.38
0.4	7.64	13.55	21.5	85.75
0.6	11.32	20.01	31.98	136.61
0.8	14.95	26.41	42.47	207.66
1.0	18.54	32.79	53.16	342.19

Equations

$$\begin{split}U(x,y,z)&=U_{9\text{-}3}(z)+\alpha\big{\{}U_{12\text{-}6}(x,y,z-l_{z_{1}})\\ &\quad-[U_{9\text{-}3}(z)-U_{9\text{-}3}(z-l_{z_{2}})]\big{\}},\end{split}$$

$$U_{9\text{-}3}(z)$$

$$U_{12\text{-}6}(x,y,z)$$

$$\begin{array}{rl}\left|{\mathbf{r}}_i\right|& =[(x-x_i{)}^2\\ & +(y-y_i{)}^2+(z-z_i{)}^2{]}^{1/2}.\end{array}$$

$$\begin{array}{rl}U(x,y,z)& =(1-\alpha )U_{9\text{-}3}(z)\\ & +\alpha [U_{12\text{-}6}(x,y,z-l_{z_1})+U_{9\text{-}3}(z-l_{z_2})].\end{array}$$

$${\sigma }_{12\text{-}6}$$

$${ϵ}_{12\text{-}6}$$

$$l_{z_2}$$

$${\rho }_{\text{LJ}}$$

$$L_z^{\prime }=L_z-\frac{{\sigma }_{\text{9-3}}+{\sigma }_{\text{O}}}{2},$$

$${\Delta }_{\text{xy}}(t)=[r(t)-r(0){]}^2,$$

$$⟨{\Delta }_{\text{xy}}(t)⟩=\frac{1}{N_{\text{cell}}{t}_{\text{steps}}}\sum _{\tau =0}^{t_{\text{steps}}}\sum _{i=1}^{N_{\text{cell}}}[r_i(t+\tau )-r_i(\tau ){]}^2.$$

$${\Delta }_{\text{xy}}(t)=4Dt.$$

$$⟨C_{\text{O-O}}(t)⟩=\frac{1}{N_{\text{cell}}{t}_{\text{steps}}}\sum _{\tau =0}^{t_{\text{steps}}}\sum _{i=1}^{N_{\text{cell}}}\frac{n_i(t+\tau )}{n_i(\tau )}.$$

$$\begin{array}{rl}{⟨U(x,y,z)⟩}_{xy}& =(1-\alpha )U_{9\text{-}3}(z)\\ & +\alpha \left[U_{10\text{-}4}(z-l_{z_1})+U_{9\text{-}3}(z-l_{z_2})\right].\end{array}$$

$$\dfrac{dU_{9\text{-}3}(z)}{dz}\Big{|}_{z_{0}}=0,$$

$$\dfrac{d\left[U_{10\text{-}4}(z-l_{z_{1}})+U_{9\text{-}3}(z-l_{z_{2}})\right]}{dz}\Big{|}_{z_{0}}=0.$$

Full text

I Introduction

Although water has been one of the most studied substances in history, it still maintains a great scientific interest due to its complex behavior still defying understanding Ball (2015); Mishima and Stanley (1998). The study of water under strong confinement has attracted special attention during the last years Bertrand et al. (2013); Kaneko et al. (2014); Chen et al. (2016); Corsetti et al. (2016a); Bai and Zeng (2012); Bai et al. (2003); Koga et al. (1997); Zubeltzu et al. (2016); Corsetti et al. (2016b); Zangi and Mark (2003); Zhao et al. (2014); Zhu et al. (2015); Giovambattista et al. (2009); Zangi (2004); Chen et al. (2017); Zubeltzu et al. (2016); Han et al. (2010). Its properties are strongly altered under these conditions, their understanding being important in biological, geological, and technological contexts Brovchenko and Oleinikova (2008); Ball (2003).

Computer simulations have been used to study the intrinsic behavior of nanoconfined water by means of simple confining potentials that depend on few parameters, avoiding possible effects that may be produced by particularities of the chosen confining substrate. Many new dynamical and structural properties of water have been observed, such as new crystalline phases Kaneko et al. (2014); Chen et al. (2016); Corsetti et al. (2016a); Bai and Zeng (2012); Bai et al. (2003); Koga et al. (1997); Zubeltzu et al. (2016); Corsetti et al. (2016b); Zangi and Mark (2003); Zhao et al. (2014); Zhu et al. (2015); Giovambattista et al. (2009); Zangi (2004); Chen et al. (2017), or unusual phase transitions Zubeltzu et al. (2016); Han et al. (2010).

Han et al. Han et al. (2010) observed by classical molecular dynamics simulations using the TIP5P force-field model Mahoney and Jorgensen (2000) that bilayer water under planar confinement undergoes unusual freezing, which is continuous or discontinuous depending on where in the phase diagram. Bai and Zeng Bai and Zeng (2012) instead observe amorphous ice formation at similar conditions. Recent studies employing the TIP4P/2005 water model Abascal and Vega (2005) combined with ab initio molecular dynamics simulations Corsetti et al. (2016a); Zubeltzu et al. (2016) explain the origin of these observations proposing a scenario in which a continuous melting occurs, with the appearance of an intermediate hexatic phase following the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory for two-dimensional melting Kosterlitz and Thouless (1973); Young (1979); Halperin and Nelson (1978); Nelson and Halperin (1979), and in which the hydrogen atoms are delocalized throughout the phase transition, while the oxygens freeze into a triangular ice. This phase, characterized by a large configurational entropy due to the large hydrogen disorder, can be quenched into amorphous ices, and undergoes a further discontinuous phase transition when the H atoms freeze, thereby accounting for the observations.

In order to obtain a deeper understanding of nanoconfined water, it is necessary to take into account the influence of the confining substrate. Computer simulations using atomistic surfaces Giovambattista et al. (2012); Fitzner et al. (2015); Al-Hamdani et al. (2017); Kiselev et al. (2017); Sosso et al. (2016) and realistic protein-like surfaces Reddy et al. (2010); Khurana et al. (2006); Cheng and Rossky (1998) have shown that the behavior of water changes drastically compared to the one observed using smooth surfaces. It has been observed in simulations with realistic surfaces that factors such as surface geometry, topography, and chemical heterogeneity play crucial roles Giovambattista et al. (2012); Fitzner et al. (2015); Mittal and Hummer (2010). It is therefore important to study the effect of different aspects of the confining system.

In this work we smoothly introduce a modulation to the planar potential employed in Han et al. (2010); Zubeltzu et al. (2016). This has been done before with large modulation periodicities Mittal and Hummer (2010) and in a different context (to study the ice nucleation on mineral surfaces of relevance for cloud formation Fitzner et al. (2015); Cox et al. (2015a, b)). We study the effects produced by the topography of the confining substrate into the liquid, hexatic and solid phases described in Zubeltzu et al. (2016). We do it by computational simulations using both, empirical potentials and first-principles calculations. For the confining walls, we propose a model based on a referential Lennard-Jones 9-3 potential plus Lennard-Jones explicit particles on its surface, allowing us to smoothly add a roughness into the confining potential and control the lattice parameter and energetic amplitude of the oscillation. We consider a triangular lattice: firstly, with a lattice parameter that commensurates with the honeycomb ice previously observed Zubeltzu et al. (2016), and secondly, three more realistic values of the lattice parameter. For the former, although the honeycomb ice stabilizes as expected, the triangular structure disfavoured by the external modulation is maintained in the liquid. The activation of the external modulation has substantial effects on the first-principles liquid, affecting the stacking and the well-defined layering observed in the planar confinement. At high temperatures and densities, the frustration produced on the hexatic phase by the external modulation stabilizes a different tri-layer ice beyond a critical amplitude of the modulation. For more realistic lattice parameters, a clear dynamical and structural inhomogeneous behavior is observed in the liquid density, although other properties are much less affected.

II Methods

We carry out computational simulations based on molecular dynamics with force-fields (FFMD), and ab initio molecular dynamics (AIMD). For the former, we use the LAMMPS code Plimpton (1995) and the TIP4P/2005 force field to model the interaction among the water molecules, where the cutoff of the Lennard-Jones interaction is set to 12 Å and the particle-particle particle-mesh (PPPM) method Hockney and Eastwood (1988) is used to compute the long-range Coulombic interactions. The rigidity of the molecules is treated using the SHAKE algorithm Ryckaert et al. (1977). The timestep is set to 1 fs. Starting from a randomly set configuration of positions and velocities, during the first 60 ns of each FFMD run we adopt the constant particle number, volume, and temperature ensemble (NVT) and use the Nose-Hoover thermostat in order to control the temperature of the system. Then, 2 ns of dynamics are collected for data analysis. The dimensions of the rectangular cell ($L_{x}^{\text{FFMD}}\times L_{y}^{\text{FFMD}}$) are different for each lattice parameter in order to keep the periodicity of the cell (see Table 1). For system size testing, see Corsetti et al. (2016a).

For the AIMD calculations based on density functional theory, we use the SIESTA code Soler et al. (2002) with the vdW-DF functional with a fully non local correlation Dion et al. (2004), devised to describe the van der Waals interactions, and PBE Perdew et al. (1996) exchange, especially tested for water Wang et al. (2011) (as in refs Corsetti et al. (2016a); Zubeltzu et al. (2016)). Previous studies with this functional have shown noticeable improvements of the calculated radial distribution functions of water due to a better description of H bonds and reproduce the maximum of the diffusion with respect to the density Corsetti et al. (2013).

The final configuration obtained from the FFMD calculations is annealed for 5 ps and then the NVE ensemble is used for at least 10 ps in AIMD while data are collected. The timestep is set to 0.5 fs. Due to the larger computational cost of such calculations, we reduce the size of the cell. Table 1 shows the dimensions of the cell ($L_{x}^{\text{AIMD}}\times L_{y}^{\text{AIMD}}$) for each lattice parameter.

The model we propose to confine water to a thin film is a linear combination of a surface made of Lennard-Jones particles and a Lennard-Jones 9-3 potential. The latter is obtained by integrating semi-infinitely a piece of matter made of Lennard-Jones particles considering a density $\rho_{\text{LJ}}$ Han et al. (2010); Zubeltzu et al. (2016). The proposed confining potential is:

[TABLE]

where:

[TABLE]

On the right hand side of Eq. (1) we add to the Lennard-Jones 9-3 potential a layer of Lennard Jones particles located at $z=l_{z_{1}}$, which introduces roughness into the potential, and remove the integrated piece of matter that represents this layer of Lennard-Jones particles \big{[}U_{9\text{-}3}(z)-U_{9\text{-}3}(z-l_{z_{2}})\big{]}. $l_{z_{2}}$ is the distance in $z$ between neighbor layers in the piece of matter made of Lennard-Jones particles. The second term on the right is multiplied by a parameter $\alpha$ which controls the strength of the roughness in the potential. Eq. (1) can be rewritten simply:

[TABLE]

In our case, the $U_{9\text{-}3}$ potential only interacts with the oxygen atoms, and its parameters are chosen to mimic the interaction of water with solid paraffin as proposed in Lee et al. (1984) and used in Han et al. (2010); Zubeltzu et al. (2016): $\epsilon_{9\text{-}3}=0.2985$ kcal/mol and $\sigma_{9\text{-}3}=2.5$ Å. The position of the confining potential is set so that the origin of the $U_{9\text{-}3}$ potential is at $z$ = 0 Å, and its symmetric potential at $z$ = 8 Å, as in Han et al. (2010); Zubeltzu et al. (2016). The two layers of Lennard-Jones on each confining side are stacked in such a way that each atom in one layer has another one directly above (below) it, a disposition we will call AA stacking. From the condition that the piece of matter made of Lennard-Jones particles of density $\rho_{\text{LJ}}$ integrates into the Lennard-Jones 9-3 potential, we can obtain the relationship between the parameters of each potential:

[TABLE]

From Eq. (6) we obtain $\sigma_{12\text{-}6}=3.498$ Å, which is independent of the lattice parameter. The values of $\epsilon_{12\text{-}6}$ for each lattice parameter are shown in Table 1.

During the integration, the Lennard-Jones particles are considered to be arranged into a fcc crystallographic configuration, where the surface is oriented in the $(111)$ direction conforming a triangular lattice. Given a triangular lattice parameter $a$, we deduce:

[TABLE]

The position of the layer of Lennard-Jones particles $z=l_{z_{1}}$, is calculated so that the position of the minimum in $z$ of the mean confining potential $\left<U\right>_{xy}(z)$ is independent of $\alpha$ (see Appendix A). There here proposed confining potential requires a redefinition of the confinement width $L_{z}$, which is needed to calculate some physical magnitudes below. We use the one proposed in Kumar et al. (2005):

[TABLE]

where $\sigma_{\text{O}}$ = 3.1589  Å is the $\sigma$ parameter set for the Lennard-Jones interaction between the oxygen atoms in the TIP4P/2005 water model Abascal and Vega (2005). This gives $L^{\prime}_{z}=0.515$ nm. Table 1 shows the values of the parameters used in the simulations for each lattice parameter. The value of $\alpha$ can be associated to an energy scale of relevance to the problem. Table 2 shows the minima at on-top $xy$ positions (on top of Lennard-Jones particles) and the absolute minima of the potential for each lattice parameter and value of $\alpha$.

To obtain dynamical information of the simulations, we calculate the mean square displacement (MSD), defined as:

[TABLE]

where, $r(t)=[x(t)^{2}+y(t)^{2}]^{1/2}$ is the position of the particle at time $t$ in the $xy$ plane. We then average this function over all the oxygens and different initial time steps:

[TABLE]

The diffusivity $D$ is obtained from the Einstein’s relation for a two-dimensional system:

[TABLE]

The oxygen first-neighbor correlation function $C_{\text{O-O}}$ gives the proportion of initial in-plane nearest-neighbors of any particle that remain after time $t$. The nearest-neighborhood is defined by a circle of radius $r_{0}$ obtained from calculating the distance at which the radial distribution function shows a minimum between the first and second-neighbor peaks. Considering the function $n(t)$ that gives the number of nearest-neighbors that remain after time $t$, the averaged first-neighbor correlation function is defined as:

[TABLE]

The lateral pressure $P_{\parallel}$ and perpendicular pressure $P_{\bot}$ are calculated using the virial expression for the $x$ and $y$ directions, and calculating the force produced by oxygens on the wall, respectively, as in Kumar et al. (2005). In order to obtain the density at which $P_{\parallel}=0$ keeping the lattice parameter constant (and hence, the dimensions of the cell), we carry out FFMD simulations on the $NVT$ ensemble with all the possible number of particles within a range such that the $P_{\parallel}=0$ point is known to be crossed. We then choose the $N$ that gives the closest value to $P_{\parallel}=0$.

III Results and Discussion

We first discuss the results obtained with a triangular lattice parameter that is commensurated with honeycomb two-dimensional ice Zubeltzu et al. (2016), and then we consider several values of $a$ in the range of what would be found with realistic confining materials, such as (111) faces of various metals.

III.1 Ideal commensuration

Taking into account that the oxygen-oxygen radial distribution function (RDF) of honeycomb ice shows an oxygen-oxygen $xy$ first-neighbor distance $r_{\text{O-O}}$ = 2.76 Å, we choose a triangular lattice with a lattice parameter $a$ = $\sqrt{3}~{}r_{\text{O-O}}$ = 4.78 Å. This produces an energy surface where the minima are positioned in a honeycomb structure with the distance between the nearest energy minima equal to $r_{\text{O-O}}$. Therefore, it should be ideal a priori for the establishment of a honeycomb ice monolayer, but it would disfavour the formation of triangular ice (one third of the oxygens would have to sit on maxima of the modulation potential). We study the structural and dynamical properties of water at different densities, temperatures, and values of $\alpha$.

III.1.1 Phase diagram

Fig. 1 shows the temperature-density phase diagrams for the system for $\alpha$ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0, obtained from the FFMD calculations. The phases appearing in the diagram for $\alpha$ = 0.0 have been previously described in Zubeltzu et al. (2016). In order to determine the phase at each calculated point, we have used four different indicators: radial distribution function, oxygen diffusion, and the oxygen $xy$ density histogram, the three of them in the xy plane, and the density profile of oxygen atoms along the confining direction, $z$. The $xy$ density histogram is obtained by averaging the $xy$ coordinates of the particles over 1 ns.

At $\rho\leq$ 1.27 g cm*-3*, the effect of increasing $\alpha$ is the expected one: the corrugation of the potential favours the structure of honeycomb ice and the melting point rapidly increases with $\alpha$. The observation of phase coexistence of honeycomb ice and liquid at different values of $\alpha$ indicates that the phase transition that connects these two phases keeps being first-order independently from the amplitude of the modulation (see Appendix B).

When $\alpha\leq 0.6$, at $\rho$ = 1.37 g cm*-3* water remains being liquid. At $\rho$ = 1.47 g cm*-3* and low temperatures, the square-tubes ice observed for flat confinement Zubeltzu et al. (2016) stays stable for $\alpha\leq 0.6$. In this solid, the molecules arrange into tubes and the hydrogen bonds tend to point towards the oxygens from the same tube (see Appendix C). The hexatic phase observed for planar confinement at $\rho$ = 1.47 g cm*-3* and high temperatures demands a deeper analysis to verify whether the orientational long-range order and translational short-range order are kept under the different values of the amplitude of the modulation.

III.1.2 Hexatic phase

As previously mentioned, the triangular modulation with the lattice parameter $a$ = 4.78 Å is expected to disfavour highly triangular structured phases, and therefore, the hexatic phase. In addition, one expects that a modulation inserting a lattice of energy maxima and minima along the plane should hinder continuous shear, and thus constrain the appearance and diffusion of the dislocations, responsible for the lack of long-range translational order of the hexatic phase Kosterlitz and Thouless (1973); Young (1979); Halperin and Nelson (1978); Nelson and Halperin (1979); von Grünberg et al. (2007). These dislocation appear when the density of the ideal triangular ice (expected to be at $\rho$ = 1.76 g cm*-3* Zubeltzu et al. (2016)) is decreased.

For the characterization of the hexatic phase, we also calculate the oxygen-oxygen first-neighbor correlation function (C${}_{\text{OO}}$), as was already used in Zubeltzu et al. (2016). From C${}_{\text{OO}}$, we calculate the characteristic time $\tau$ of the C${}_{\text{OO}}$ curve, which is defined by C${}_{\text{OO}}(\tau)=$ 0.5. For $\alpha\leq 0.4$, the oxygens arrange themselves into a triangular lattice showing similar structural and dynamical features as the ones observed in the hexatic phase at $\alpha$ = 0 Zubeltzu et al. (2016): similar RDFs (see Appendix D), diffusivities ($D=$ 6 10*-7* cm2s*-1* at $T$ = 300 and $\alpha$ = 0.4 versus $D=$ 9.5 10*-7* cm2s*-1* for $\alpha$ = 0.0 Zubeltzu et al. (2016)) and oxygen-oxygen first-neighbor correlation functions [Fig. 2 (a)]. We also observe the existence of shear motion along the main directions of the triangular lattice in the oxygen $xy$ density histogram [Fig. 2 (b)]. All these indicators combined suggest that the hexatic phase remains stable for $\alpha\leq 0.4$. At $\alpha$ = 0.6, however, although there is still a clear structuring of the oxygens into a triangular lattice, we do not observe shear motion of oxygens [Fig. 2 (c)]. The increase in the diffusivity ($D$ = 3 10*-6* cm2s*-1*), decrease of the characteristic time of the C${}_{\text{OO}}$ function in Fig. 2 (a), and the less structured RDF (see Appendix D) suggest that the hexatic phase is sufficiently frustrated to be no longer stable under these conditions. Although the triangular structuring is still present in the $xy$ density histogram of the oxygens, water seems to behave as a dense triangular liquid. Therefore, the indicators employed in this section suggest that the continuous melting reported in Zubeltzu et al. (2016) washes out for $\alpha\geq 0.6$.

III.1.3 Intercalated honeycomb ice

The structural and dynamical properties of the honeycomb ice, square tubes ice, and the triangular phases appearing in Fig. 1 have been already well described in previous works Corsetti et al. (2016a); Zubeltzu et al. (2016) and $\alpha$ does not affect them significantly (see Appendices B and E). However, when $\alpha$ = 0.8 and 1.0 a different form of ice stabilizes at $\rho$=1.47 g cm*-3*, which we call intercalated honeycomb ice. Fig. 3 shows the $xy$ density histogram of oxygens and the $z$ density profile characteristic of this solid. The oxygen atoms arrange themselves into three layers: the two exterior layers conform a AA stacked honeycomb structure favoured by the energy landscape produced by the Lennard-Jones confining particles. The oxygens from the central layer fill the centre of the honeycomb hexagons coinciding with the $xy$ coordinates of the confining Lennard-Jones particles. Between the optimal densities of the honeycomb ice ($\rho$=1.17 g cm*-3*), and the intercalated honeycomb ice ($\rho$=1.47 g cm*-3*), we observe regions of coexistence of both phases independently of the temperature (see Appendix B). This is a signature of having a first-order phase transition in a $NVT$ ensemble.

III.1.4 Liquid

The effect of the corrugation on the structure of the liquid is different depending on its density. For densities ranging between $\rho$ = 0.97$-$1.17 g cm*-3*, the liquid rapidly freezes at the highest temperatures considered as $\alpha$ increases, as shown in Fig. 1. Before it freezes, however, the liquid does not show appreciable changes on its structure (see Appendix F). At higher densities, the liquid resists to freeze, specially at $\rho$ = 1.37 g cm*-3*, which only freezes after $\alpha$ = 0.8. Fig. 4 shows the RDFs, and $z$ density profiles of liquid water at $T$ = 300 K, and $\rho$ = 1.37 g cm*-3* for different values of $\alpha$. The vertical dashed lines on Fig. 4(a) highlight the position of the second and third neighbor peaks in triangular ice $r_{h}$ = $\sqrt{3}r_{\text{O-O}}$ = 4.78 Å and $r_{t}$ = 2$r_{\text{O-O}}$ = 5.52 Å, respectively, the latter being absent in honeycomb ice Zubeltzu et al. (2016). With no corrugation, the liquid at this density shows structural features on the RDF and density profile characteristic of triangular ice Zubeltzu et al. (2016): high inter-planar correlation (high value on $g_{\text{O-O}}(0)$), a pronounced peak on the RDF at $r_{t}$, and two pronounced peaks on the density profile with little density and flux of molecules between them. When $\alpha$ is applied, the RDFs barely increase their value at $r_{h}$, which shows the clear resistance of the liquid to maintain the triangular structure despite the fact that it is disfavoured by the external modulation. The liquid, instead of becoming more honeycomb-like, shows a destructuring effect under the modulation (Fig. 4): as $\alpha$ increases, the RDFs show smaller peaks and decrease their value at the origin that measures the inter-layer correlation. The two peaks on the density profile decrease while the flux of molecules between them increases. The destructuring tendency with the external modulation is also supported by the observation of the disappearance of the solid-liquid coexistence point observed at $T$ = 220 K, $\rho$ = 1.37 g cm*-3*, and $\alpha$ = 0, and 0.2 in the phase diagram (Fig. 1).

We compare the effect of the corrugated wall on the liquid with FFMD calculations with the one obtained with AIMD calculations. The red points in Fig. 1 represent the density, temperature, and $\alpha$ values for which AIMD calculations have been done. The direct comparison of the RDFs obtained at $\rho$ = 1.37 g cm*-3* and $T$ = 300 K (Fig. 5) show the good agreement between the structural results obtained by the two different methods. Previous works Corsetti et al. (2016a); Zubeltzu et al. (2016) with planar confinement, together with the ones shown here, indicate that the results obtained with the TIP4P/2005 empirical force-field and the ab inito molecular dynamics with the vdW-DF${}^{\text{PBE}}$ functional agree surprisingly well.

The small differences between the RDFs obtained by AIMD and FFMD calculations agree with what presented in previous reports: the RDFs obtained by AIMD are more structured, showing larger correlation peaks positioned on the O-O distances characteristic of the triangular ice. Both calculation methods show that the liquid is triangularly structured independently of the density. It is particular of the AIMD RDFs however, that the activation of the modulation affects significantly the maximum at $r_{\text{xy}}$ = 0, which measures the interplanar on-top correlation. The peak at the origin of the RDF is characteristically pronounced for a triangular liquid where there is a tendency for AA stacking Zubeltzu et al. (2016). The peak decreases abruptly for AIMD when $\alpha$ = 0.2 for all the sampled densities. In order to understand this behavior, we calculate the interlayer RDF, which takes into account the correlation of an oxygen atom with the oxygens of the other layer.

Fig. 6 shows the interlayer RDFs and density profiles at $T$ = 300 K and $\rho$ = 1.37 g cm*-3* for different modulation amplitudes obtained by AIMD calculations. The $\rho(z)$ density profiles show that once the corrugation is activated the two water layers get noticeably closer to each other. The intensity of the peaks decreases as the amplitude of the modulation increases. As the AIMD liquid is more triangularly structured than the FFMD one, and the external modulation frustrates triangularly structured phases, the effect of the modulation in the AIMD calculation is greater than in the FFMD case. Therefore, the oxygens closer to maxima of the modulation tend to locate closer to the center of the film displacing the peaks of the density profiles. This displacement explains the loss of AA stacking: the distance between the oxygen layers becomes smaller than the optimal distance for a vertical H-bond, and therefore, the vertical bonds get tilted with respect to the $z$ axis. The effect on the stacking is also seen by the slight differences of the peak away from $r_{\text{xy}}=0$ in Fig 6 (a), which is accompanied by a small displacement of the other peaks towards smaller distances, reflect of the change in stacking correlations.

The obtained diffusivities for the liquid and solid phases for all the values of $\alpha$ are very similar to the ones previously calculated in the absence of modulation Zubeltzu et al. (2016) and observed for bulk water Corsetti et al. (2013): At $T$ = 300 K, the liquid has a diffusivity of the order of $D$ $\sim$ $10^{-5}$ cm2s*-1* while the solid phases $D$ $\sim$ $10^{-8}-10^{-9}$ cm2s*-1*. The calculations of the mean-square displacement (MSD) at the different points of the phase diagrams show that for a given $T$ and $\rho$, the diffusivity does not significantly change with $\alpha$, as long as a phase transition does not occur (see Appendix F).

III.2 Realistic lattice parameters

After analyzing a lattice parameter that is ideally commensurated with honeycomb ice, we study the effect of corrugated walls with more realistic lattice parameters. We chose three different triangular lattices with $a$ = 2.50, 2.75, and 3.00 Å that cover a range of typical values for closed-packed metal surfaces, such as Ni (2.49 Å), Cu (2.56 Å), Pt (2.78 Å), Au (2.88 Å), and Ag (2.89 Å). In this section we restrict our simulation to $T$ = 300 K, which corresponds to the highest $T$ in Fig 1.

III.2.1 Structure

We carry out FFMD calculations for the same densities as in the previous section, and $\alpha$ = 0.2, 0.4, 0.6, 0.8, and 1.0. We observe no significant change in the different phases under the effect of the corrugation, staying liquid at low densities and triangularly structured at high densities independently of $\alpha$ [as in Fig. 1(a) at $T$ = 300 K]. The high-density triangular phase shows the features of the hexatic phase previously mentioned: similar RDFs, and a clear triangular lattice in the $xy$ density histograms of the oxygens with the usual shear motions along the main directions of the lattice (see Appendix G). The liquid shows almost no change in the RDFs and in the density profiles, when varying $\alpha$. Fig. 7 shows the RDFs for $a$ = 3.0 Å at $\rho$ = 1.17, 1.27, and 1.37 g cm*-3* (different vertical shifts for each $\rho$), and $\alpha$ = 0.2, 0.6, and 1.0. When $r_{\text{xy}}>$ 1.0 Å, the RDFs with the same density and different values of $\alpha$ are almost indistinguishable from each other. The peak around $r_{\text{xy}}$ = 0 Å, which contains information about the interlayer on-top correlation, shows a small tendency to increase with $\alpha$ (inset in Fig. 7) the opposite to the behavior for the ideal $a$. This means that as the corrugation increases, the liquid tends to structure into better AA stacking. From these results we deduce that for $a$ = 3.0 Å, although the relative distances among oxygen atoms are almost unaffected by the corrugation, the water molecules from one layer tend to be on top of another from the other layer as the corrugation increases. The figures are very similar for $a$ = 2.5, and 2.75 Å and the same conclusions can be drawn for them (see Appendix H).

Although the RDFs and density profiles of oxygens suggest that the structure of the liquid is barely affected by the corrugation, the $xy$ density histogram of the oxygens shows clear effects on the structure of the liquid. The insets in Fig. 8 show the oxygen $xy$ density histogram for $a$ = 3.0 Å, at $\rho$ = 1.27 g cm*-3*, and $\alpha$ = 0.2, 0.6, 1.0. When the corrugation is activated, the oxygens avoid being close to the Lennard-Jones confining particles due to the effective repulsive force they represent, and they are structured anisotropically, resulting in the heterogeneous images shown in the insets of Fig. 8. The maximum value of the density oscillations in the insets are $138\%$, $186\%$, and $237\%$ of the mean value of the density for $\alpha$ = 0.2, 0.6, and 1.0, respectively. To quantify the effect of this structuring of the liquid, we calculate the dipole distribution function (DDF): we average the distribution of the polar angle (projected in the $xy$ plane) of the molecular dipoles. Fig. 8 shows the DDFs for $a$ = 3.0 Å, at $\rho$ = 1.27 g cm*-3*, and $\alpha$ = 0.2, 0.6, 1.0. We can clearly observe that as $\alpha$ increases, the DDFs show six pronounced peaks in the multiples of 60*∘*. This anisotropic structuring effect is observed for all sampled points of the liquid phase, and it is greater as $\rho$, and $a$ are increased. These results show that the (not surprising) anisotropy in the molecular density of the liquid has a surprisingly small effect in the liquid structure (correlation) and dynamics (diffusivity), in spite of the orientational effect.

We now compare the results obtained by FFMD and AIMD calculations at $T$ = 300 K, and $P_{\parallel}$ = 0, for different values of $\alpha$ and $a$. Fig. 9 shows the RDFs obtained by FFMD and AIMD calculations at $a$ = 2.50, 2.75, and 3.00 Å, and $\alpha$ = 0.2, 0.6, and 1.0. The effect of the corrugation in the RDFs is as negligible for AIMD as it was for FFMD. From the direct comparison of the FFMD and AIMD RDFs we obtain the same conclusions to the ones in the previous section: both methods give similar structural features of the liquid, but with the difference that the AIMD RDFs tend to be more triangularly structured showing a larger peak at $r_{\text{t}}$ = 2$r_{\text{O-O}}$. This is also supported by the comparison of the density profiles $\rho(z)$ obtained by both calculation methods: although they are very similar, the ones obtained by AIMD show more pronounced peaks. These results agree with the results obtained for $a$ = 4.78 Å and previous works where the same empirical force-field and DFT functional Corsetti et al. (2016a); Zubeltzu et al. (2016) were used. The differences between both methods are minor.

III.2.2 Diffusivity

As previously mentioned, with no corrugation ($\alpha$ = 0) the sampled liquid points of the phase diagram show diffusivities of the order of $D\sim 10^{-5}$ cm2s*-1*, while for the hexatic phase $D\sim 10^{-7}$ cm2s*-1* Zubeltzu et al. (2016). When the corrugation is applied with these three lattice parameters, and different values of $\alpha$, the diffusivities barely change in the liquid and in the triangular phase (see Appendix H). The small changes with $\alpha$ in the mean square displacements curves are within the noise of the signal and do not display any clear trend; therefore, we conclude that within our accuracy, the diffusivities are barely affected by the corrugation formed with these three lattice parameters.

IV Conclusions

In this work we simulate two-dimensionally confined water between two corrugated walls using empirical potentials (TIP4P/2005) and first principles. We propose a periodic confining potential based on Lennard-Jones particles that play the role of a confining substrate. We control the lattice type, lattice parameter and amplitude of the corrugation.

For $a$ = 4.78 Å, ideal for the stabilization of two layers of honeycomb ice, at low-mid densities, the honeycomb ice rapidly stabilizes within the phase diagram as the modulation amplitude increases. However, before the liquid freezes, its triangular structure, although disfavoured by the external modulation, remains almost unchanged, independently of its density. This is supported by AIMD results. At high densities and low temperatures, the square-tubes ice keeps being stable for $\alpha\leq$ 0.6. The hexatic phase appearing for planar confinement is destroyed for $\alpha\geq$ 0.6, which seems to liquify. When $\alpha\geq$ 0.8, a new solid phase appears, intercalated honeycomb ice. This solid is composed of two external honeycomb layers and a third middle layer with the molecules located in the centre of the honeycomb hexagons. The AIMD liquid tends to be slightly closer to triangular structure as previously reported for $\alpha$ = 0.0. The imposed modulation disfavours triangular-like structures, and therefore affects more significantly the AIMD liquid, which loses the AA stacking and well-defined layering shown in the planar confinement.

For more realistic lattice parameters ($a$ = 2.50, 2.75, and 3.00 Å) we do not observe significant changes in the phase behavior, staying liquid at low-mid densities, and triangular with similar features to the hexatic phase at high densities. Although the $xy$ density histograms of the oxygens and DDFs show that the molecules from the liquid display a significant anisotropy, there is hardly noticeable change in the structural (RDF) and dynamical (diffusivity) behavior. Finally, we have found that the differences between AIMD and TIP4P/2005 which we set out to explore, have shown to be significantly smaller than anticipated: after all, these structures are quite removed from the bulk liquid the forces were fitted to.

V Acknowledgements

This work was partly funded by grants FIS2012-37549-C05 and FIS2015-64886-C5-1-P of the Spanish Ministerio de Econom a, Industria y Competitividad, and Exp. 97/14 (Wet Nanoscopy) from the Programa Red Guipuzcoana de Ciencia, Tecnología e Innovación, Diputación Foral de Gipuzkoa. We thank José M. Soler and Fabiano Corsetti for useful discussions. The calculations were performed on the Arina HPC cluster (Universidad del País Vasco/Euskal Herriko Unibertsitatea, Spain) and MareNostrum (Barcelona Supercomputing Center). SGIker (UPV/EHU, MICINN, GV/EJ, ERDF and ESF) support is gratefully acknowledged.

Appendix A Position of Lennard-Jones particles

We calculate the distance between the origin (divergence) of the confining wall and the Lennard-Jones confining particles $l_{z_{1}}$ by making the position of the minimum in $z$ of the mean total confining potential $\left<U(x,y,z)\right>_{xy}$ independent from $\alpha$. Taking into account that a layer made of Lennard-Jones particles located at $z$ = $l_{z_{1}}$ integrates into the Lennard-Jones 10-4 potential, $\left<U_{12\text{-}6}(x,y,z-l_{z_{1}})\right>_{xy}=U_{10\text{-}4}(z-l_{z_{1}})$, the total mean confining potential is:

[TABLE]

We then find the coordinate $z$ = $z_{0}$ at which the first term on the right has its minimum,

[TABLE]

and we obtain $l_{z_{1}}$ by imposing the second term on the right to have its minimum at the same coordinate $z_{0}$:

[TABLE]

Following this procedure, the minimum of the total mean potential is located at $z_{0}$ and is independent from $\alpha$.

Appendix B Phase coexistence for $a$ = 4.78 Å

For $a$ = 4.78 Å we observe phase coexistence in some regions of the phase diagrams in Fig. 1. Fig 10 shows the $xy$ density histogram of the oxygen atoms where the coexistence of the liquid-honeycomb and honeycomb-intercalated honeycomb phases are observed.

Appendix C Square tubes ice for $a$ = 4.78 Å

For $a$ = 4.78 Å, and $0.2\leq\alpha\leq 0.6$, we observe the similar high density phases observed in Zubeltzu et al. (2016). At low temperatures, we observe the formation of the square tubes ice, where both the hydrogen and oxygen atoms are fixed, with the hydrogen atoms tending to point towards the molecules from the same tube. The phase is observed to stay stable for $\alpha\leq$ 0.6. Fig. 11 (a) and (b) show the $xy$ density histogram of the oxygen and hydrogen atoms at $T$ = 220 K and $\rho$ = 1.47 g cm*-3* for a flat confining potential ($\alpha$ = 0) and for $\alpha$ = 0.6, respectively.

Appendix D Hexatic phase for $a$ = 4.78 Å

For $a$ = 4.78 Å at $\rho$ = 1.47 g cm*-3* and $T$ = 300 K corresponding to the reported hexatic phase for $\alpha$ = 0.0 Zubeltzu et al. (2016) similar structural features are observed in the RDFs up to $\alpha\leq 0.4$ (Fig. 12) while it clearly changes at $\alpha$ = 0.6.

Appendix E Honeycomb ice for $a$ = 4.78 Å

We observe that for $a$ = 4.78 Å the honeycomb ice stabilizes with $\alpha$ (Fig. 13). When $\alpha$ = 1.0 it is stable at $\rho\leq$ 1.17 g cm*-3* for all the sampled temperatures (see Fig. 1).

Appendix F Liquid for $a$ = 4.78 Å

The liquid shows similar structural and dynamical features before it freezes for $a$ = 4.78 Å (Fig. 14). The RDFs and MSDs at the different points of the phase diagrams show that for a given $T$ and $\rho$, the diffusivity and structure do not significantly change with $\alpha$, as long as a phase transition does not occur.

Appendix G Hexatic phase for realistic lattice parameters

The hexatic phase shows similar features under the modulation for realistic lattice parameters (Fig. 15). The $xy$ density histograms show that the oxygen atoms are arranged into a triangular lattice with the usual shear motion along the main directions of the lattice. The MSDs show similar diffusivities for different realistic lattice parameters ($a$ = 2.50, 2.75, and 3.00 Å).

Appendix H Liquid for realistic lattice parameters

The liquid shows similar structural and dynamical features before it freezes for realistic lattice parameters (Fig. 16). For the different lattice parameters ($a$ = 2.50, 2.75, and 3.00 Å) the RDFs are almost indistinguishable at a given $\rho$ and the values of the diffusivity of the liquid are around $10^{-5}$ cm2s*-1* as for a planar confinement Zubeltzu et al. (2016).