Machine Learned Coarse Grain Water Models for Evaporation Studies

Evaporation studies of water using classical molecular dynamics simulations are largely limited due to its high computational expense. We aim at addressing the computational issues by developing a coarse grain model for evaporation of water on solid surfaces by combining four water molecules into a single bead. Most commonly used mono atomic pair potentials like Lennard Jones, Morse, Mie and three body potential like Stillinger-Weber are optimized using a combination of Genetic algorithm and Nelder-Mead algorithm. Among them, Stillinger-Weber based model shows excellent agreement of density and Enthalpy of vaporization with experimental results for a wide range of temperatures. Further, the new water model is used to simulate contact angle of water and thin ﬁlm evaporation from surfaces with diﬀerent wettabilities. plate. The interaction parameters between


Introduction
With the advent of super computers and data center based cloud computing, the energy consumption of electronic sector has significantly increased over the past decade and the trend will continue. This demands compact designs with higher energy density and hence calls for highly efficient and faster cooling systems. The cooling system of a data center typically consumes around 40% of the energy input. This forecast to around $150 billion per year cost by the year 2030 1 and hence efficient cooling systems becomes a hot topic to research. Of many conventional cooling techniques, evaporative water cooling stands out as a best candidate due to its high specific heat capacity and high enthalpy of vaporization.
Based on kinetic theory, 2 an evaporative heat flux of 20,000 W/cm 2 can be achieved using water at 1 atm. However, even with the recent developments in Nano-scale fabrication techniques the maximum heat flux is limited to 500 to 1000 W/cm 2 . 3 This points out to our poor understanding of the nanoscale and microscale evaporation of water.
While experimental techniques pose challenges to validate and characterize the nanoscale phenomena, the computational tools like molecular dynamics can give promising results without the need for expensive machines and instrumentation. Studies show that thin film evaporation has high potential to remove heat compared to bulk region. 4,5 To study and extend these studies further, we need computational models with accurate representation of water and reproducibility of experimental results.
Studying water evaporation at nanoscale using experiments is a challenging attempt.
Use of molecular dynamics can become helpful in this scenario, 6-8 however, to date, there is no single water model which can capture or simulate all of its properties. 9 Even the best performing and widely used models like SPCE 10 and TI3P 9 are unsuitable to study evaporation, due to the computational cost. This shifts our focus to computationally faster models called coarse grain molecular dynamics (CGMD) models.
In a typical CGMD model, one or more water molecules are combined into a bead or a super-molecule to represent the bulk properties of the system. Most of the existing CGMD models are developed for bio-molecular studies [11][12][13][14][15][16] and is not tested for heat transfer. The existing CGMD models 13,17 mainly focus on the room temperature behavior of water or even sometimes the behavior below zero degree Celsius. While these models can capture freezing, ice formation and other properties of water, 18 their applicability to high temperature applications is limited or seldom. The field of coarse graining itself is a big research area and hence the readers are advised to refer consolidated reviews found elsewhere. [19][20][21] Among various types of CGMD models, mono bead models are appealing due to its very low computational power consumption. 11,22,23 Specifically the generalized Mie potential model 23 which maps one water molecule to a single bead can accurately represent density, surface tension and enthalpy of vaporization of water at a wide range of temperatures. Compared to this generalized Mie potential model, classical MARTINI model 11,22 combines 4 water molecules into one bead, which makes it computationally efficient.
In this work, we have utilized the combination of meta-heuristic 24 Genetic Algorithm (GA) 24 and Nelder-Mead (NM) method 25 for finding the optimum values of force field 26 parameters. The first step is GA in which the initial sample population is created randomly and optimization is performed by imitating natural selection, cross-over, mutation and reproduction. 27 The process is repeated until the best cost value reaches a threshold or upon reaching certain number of iterations. If the cost is not converged within the conditions set then the NM is initiated with an initial simplex. 25 The simplex is iteratively moved in the multi dimensions in direct search of the optimum parameters and terminates when the standard deviation of the simplex reaches a threshold or upon reaching certain number of iterations. The optimized force field parameters are then tested for reproducibility and its thermodynamic properties in the temperature range of 0°C to 170°C. Our machine learning results show Stillinger-Weber (SW) potential simulates density and Enthalpy of vaporization and agrees with the experimental results accurately.
The machine learned Stillinger Weber (MLSW) model is used to study the contact angle of water droplet on the surface of a platinum plate. The interaction parameters between the platinum and the MLSW beads are tuned to simulate hydrophobic and hydrophilic surfaces. Further, the thin film of water evaporation from the platinum surface is simulated and characterized for various wettability conditions.

Coarse Grain Models
The concept of coarse-grained (CG) modeling originates back from the 1970s. 28 Coarsegrained models assume various levels of reduced representation of the molecules, most of the time combining a bunch of atoms into a bead (as shown in Fig. 1 ) . These beads may be interconnected using bonds representing bio molecular chains or can be solo particles representing a fluid medium like water. The field of coarse grain molecular dynamics is broad and discussion of the methods is beyond the scope of this article. A few review papers discussing the recent innovations in the development of CG models, which broadly include structure-based, knowledge-based, and dynamics-based approaches can be found in the literature. 29,30 The most widely used CG model is called MARTINI model, 11 which originally developed to simulate lipid molecules. 22 The MARTINI force field maps four heavy atoms to one CG interaction site and is parameterized to reproduce the thermodynamic properties. Simulating water heat transfer or evaporation using original MARTINI model as it is will produce incorrect results, but it is re-parameterized to the temperatures near 100°C. 31 A number of coarse grain models for water have emerged recently with the developments in the machine learning and artificial intelligence field. 18,32,33 Among them, we have selected the most widely used and classical potential models for mono atomic beads. They are (i)

Morse Potential
The Morse potential, is a convenient inter atomic interaction model for a diatomic molecule.
It is a better approximation for the vibrational structure of the molecule than the quantum harmonic oscillator 38 because it explicitly includes the effects of bond breaking, such as the existence of unbound states 34 . Though it was originally developed for diatomic molecules, the Morse potential can also be used to model other interactions such as the interaction between an atom and a surface.
The Morse potential energy function takes the form as shown in the Eqn. 1. Here r is the distance between the atoms, r 0 is the equilibrium distance, D 0 is the potential well depth, and α controls the spread of the potential function For numerical simulations, often a shift and a linear term is added to the potential so that both, potential energy and force, go to zero at the cut-off radius (Eqn. 2). Generally a cutoff radius is used to speed up the calculation using computer, so that atom pairs which distances are greater than the cutoff have an interaction energy of zero.

Mie Potential
Mie potential was proposed by Gustav Mie in 1903. 35 The functional form of Mie potential is given as in Eqn. 3. When γ R = 12 and γ A = 6 then Mie potential becomes a standard Lennard Jones potential. The Mie potential is investigated to represent CG model of water to match the vapor pressure and liquid-vapor co-existence. 23 Here, σ is the value of r at U M ie (r) = 0, is the well depth (energy) and the constant C is given by the Eqn. 4. Also, γ R and γ A are the constants controlling the repulsive and attractive strength of the potential.

Lennard Jones Potential
Lennard Jones (LJ) potential 36 developed by John Lennard-Jones in 1924, may be the most used inter atomic potential across all force fields in computational chemistry. The mathematical form of the potential is simple as given in the Eqn. 5 and has only two parameters to fit the thermodynamic properties of the fluid with experimental data.
A variation of LJ potential for van der Waals interactions used in computer simulations of organic compounds, molecular liquids and divalent metal cations, is the COMPASS force field. 39 The potential takes the functional form as given in the Eqn. 6.
Here, in both cases, σ is the value of r at U (r) = 0, is the energy well depth.

Stillinger Weber Potential
The Stillinger-Weber (SW) potential 37 has a two-body and three-body terms of the standard form as shown in the Eqn. 7.
The two body interaction term is φ 2 and is given in Eqn. 8.
The three body interaction term is φ 3 and is given by Eqn. 9.
The SW potential is originally created for simulating Silicon metal interaction, however later has been extended to simulate other elements and also water. 18,40 The readers are advised to refer the original Stillinger Weber research paper 37 for the nomenclature and details of the parameters shown in Eqns 7,8,9. Genetic Algorithm The potential models described in the previous sections are used to map four water molecules into a CG bead, just like the standard MARTINI model. 11 This means, a CG bead represents 12 atoms which can significantly reduce the computational cost. To have a consistent comparison among different potentials, 4 water molecules are mapped to a bead for Mie, Morse, LJ 12-6, LJ 9-6 and SW potential models. For every model, the initial population (parameters) are created based on the original atomic or CG versions with a 10% spread using a uniformly distributed random noise using MATLAB. 41 This random population is used to search for an optimized parameter set using Genetic Algorithm.
Genetic algorithm (GA) is a meta-heuristic inspired by the process of natural selection that belongs to the larger class of evolutionary algorithms (EA). Genetic algorithm is commonly used to optimization and search problems by relying on biologically inspired operators such as mutation, crossover and selection. For our work, the selection, pairing, mating and mutations are performed based on the textbook by Haupt and Ellen. 42 The selection rate is 50% and mutation rate is 20%. The mating step follows the BLX-α method. 43 The population size is selected as 50, and maximum number of generations is set as 200. To save the computational effort, the iterative process is stopped when cost GA value falls below 5% or the standard deviation of 10 consecutive cost values falls less than 1%. The steps involved in the parameter search using GA is shown in the Fig. 2a.
All Coarse Grain Molecular Dynamics (CGMD) simulations are performed using LAMMPS software. 44 The temperature and pressure is controlled using a Nose-Hoover thermostat and barostat 45-47 with a chain length of four. The simulations start with an energy minimization for 1000 steps, followed by equilibration runs of 50,000 steps and production runs of 50,000 steps. The time step of integration is chosen as 5f s which is relatively low compared to a standard CGMD simulation. Number of steps for equilibration and production are selected based on the time required for the equilibration of energy, temperature and pressure and the standard deviations of energy and temperature is within 1%. The cut off radius for all the simulations are taken as 12Å, same as the original MARTINI model.
The density is estimated using a liquid cubic system of 1000 beads using NPT ensemble, which is pre-equilibrated for a density of 1000kg/m 3 . The Enthalpy of vaporization is estimated based on both liquid cubic system of 1000 beads (density of 1000kg/m 3 ) and vapor This comparison is made using a percentage error based cost function as described in the Eqn. 10. Initially, surface tension values of water was also included in the optimization loop.
However, since it lead to a never converging scenario, we excluded the surface tension from the list.
Here, ψ represents the quantities (density and enthalpy) to compare with experimental results, subscripts EXP and CGM D means the values from experiments and simulations respectively. The summation runs from i to D, which represents the number of quantities used to compare (2 in our case). Thus obtained cost values for each CG water models are shown in Fig. 2b. The legends represents the corresponding CG water models with MLSW being Machine Learned Stillinger-Weber model. and the rest (LJ 12-6, LJ 9-6 and Mie) will be unsuitable.
Even though Genetic algorithm is an efficient technique, like in many problems, GA have a tendency to converge towards local optima or even arbitrary points rather than the global optimum of the problem. In such occassions, local optimization algorithms like Nelder-Mead 25 or Levenberg-Marquardt algorithm. 49 The SW models are proven to accommodate a broader temperature range 18 and hence we will further parameterize SW and Morse potential based CG models using NM algorithm as described in the next section. There exists many variations for the numerical implementation of NM method, of which we used the original method 25 to further optimize the parameters of Morse and MLSW potential based CG water models. The algorithm of NM implementation in our CG strategy is shown in Fig. 3. The final parameter set obtained from the GA optimization is used as the initial simplex for the iterative NM method. The integration of NM method with the Coarse graining procedure is shown as in Fig. 3a.

Nelder-Mead Multi-Temperature Optimization
Based on the parameters (simplex) selected, the cost function is estimated using the Eqn. 11. This time, instead of comparing and optimizing the results for a single temperature, thermodynamic properties corresponding to four different temperatures are considered. This multi point (Temperatures of 0°C, 50°C, 100°C, 150°C) optimization will allow a better agreement of the CG model with the experimental results.
Here, cost N M is the cost function value associated with the NM method, s 1 and s 2 are scales, ρ is density and ∆H is the enthalpy of vaporization corresponding to the temperature T . The subscript i runs from 1 to 4, representing T 1 = 0°C to T 4 = 150°C. The scales (for our case, s 1 = 0.25 and s 2 = 0.75) are used for controlling the convergence between two quantities whose percentage error makes big difference in absolute error.
The cost function value associated with the NM optimization is shown in the Fig. 3b.
Both CG models Morse and MLSW didn't show any further convergence after 50 iterations.
The cost value at the 200 th iteration is 5.46% and 1.44% for Morse and MLSW respectively.
The corresponding parameters of the Morse CG force field (Eqn. 2) and MLSW CG force field is shown in Table 2 and Table 3 respectively.  These are the final optimized parameters of the CG water models Morse and MLSW. In the next section, further validation of these models and sensitivity testings will be explained.

Machine Learned Models and Sensitivity Studies
The water CG models based on Mie, Lennard Jones, Morse and MLSW will be used to estimate the density, enthalpy of vaporization, and surface tension. These properties are estimated for the temperature range of 0°C to 150°C with an interval of 10°C.

Sensitivity studies
The initial optimization of the parameters are based on the cutoff radius 1.2 nm, number of water molecules mapped into a bead as 4, and a time step of integration of 5 f s. Additional studies are performed to check how these simulation parameters affect the force field's ability to reproduce thermodynamic properties like density and Enthalpy. For MLSW, the cutoff radius is part of the potential function and is optimized as the part of the CG algorithm.
For Morse, the cutoff radius sensitivity is tested by changing it to 1.2 nm , 1.8 nm and 2.2 nm. The resulting effect in the density and the Enthalpy is shown in the Fig. 5a and 5b respectively. The density and ∆H is very sensitive to the cutoff radius. Both quantities increase with increasing cutoff radius and hence won't be able to predict the correct thermodynamic properties. To study the computational speed of various CG models, we varied the system size with 100 beads (12 nm 3 ) to 20,000 beads (1800 nm 3 ) at 400 K for 5,000 steps. The simulations are performed using LAMMPS 44 and the time consumed in seconds using a single processor is shown in Fig. 6a. The results show LJ 12-6 potential is fastest among the three. Despite of being a 3-body potential, MLSW is 7 to 22 % faster than the Morse depending on the system size.
The time step of integration is chosen in a way that it will have minimum impact in the total energy fluctuation. It was varied from 1f s to 100f s. A temperature of 500K is used for testing the sensitivity for higher temperatures, with a system size of 1000 beads, equilibrated for 50000 steps using NVT ensemble and production runs of 50000 using NVE ensemble. The standard deviations of total energy for Morse and MLSW cases are shown in the Fig. 6b. The system became unstable for Morse potential after 70f s and for MLSW after 40f s. The ideal choice of time step for Morse potential can be somewhere between 30f s to 40f s. On the other hand, the time step is sensitive with MLSW potential. Since our validation and optimization studies are performed at 5f s, we will be using 5f s as the time step for the future studies using MLSW potential.
Based on our studies and analysis so far, MLSW predicts density and ∆H better than Morse potential. One might often be interested in simulating water using Morse potential due its simplicity in implementing. Hence, we used GA algorithm again to search for the optimum parameters with different cutoff radii, number of water molecules per bead, time step of integration etc. The maximum number of generations are set as 100, and the cases shown in the Table 4 is optimized for the temperature 100°C using GA. The corresponding results of different Morse models (A, B, C, D) are shown in Table 4, with the optimized cost in the last column.

Results and Discussion
In comparison with other machine learned coarse grained models, MLSW was able to reproduce the density, enthalpy of vaporization and surface tension better. Hence we chose MLSW model for studying the contact angle and thin film evaporation properties of water at various types of surfaces.

Evaporation of thin water film
Using the contact angle study, we could identify the wall-water interaction parameters for creating super-hydrophobic, hydrophobic and hydrophilic surfaces. We extend this system to further study evaporation of a thin water film and its effect on surface wettability 54 characteristics. Two parallel platinum (FCC 100) plates at 15 nm apart are taken with an MLSW water film of 6 nm thickness is equlibrated next to the lower plate as shown in the Fig. 8 a, b. The dimensions of the system are shown in the Fig. 8 b, and Fig. 8 a represents multiple periodic images to visualize the system. A width of 3.6 nm (3 times cutoff radius) is chosen to avoid any self-interaction of beads across periodic boundaries. The interaction of MLSW beads with each other is using MLSW potential (Table 3), and with the platinum wall is using LJ 12-6 potential. A time step of 5 fs is used for the studies with an equilibration of films for 300,000 steps (1.5 ns) at 300 K followed by a 300,000 steps (1.5 ns) of production runs at 400 K. To simulate the evaporation and condensation of the thin films, the aforementioned system is defined with heating and cooling regions as shown in the Fig. 8 c. The temperature of lower region is kept at 400 K and that of the upper region is maintained at 300 K using a simple velocity re-scaling algorithm. This temperature gradient led to a sudden jumping of half of the film from the lower region to upper plate. Eventually the remaining liquid film is evaporated and deposited on the top surface as shown in the Fig. 8  due to strong wall-water interactions. The results show that hydrophobic surfaces are more favorable towards evaporation than the hydrophilic ones. The quantification of the evaporation data using the trajectory evolution is done using a self written Python computer code and shown in Fig. 9. The Fig. 8

Conclusion
In this work we used the most commonly used molecular force fields to develop the best coarse grain model to perform evaporation studies of water at nano to micro scale. We combined 4 water molecules into a coarse grained bead and using machine learning algorithms we examined Lennard-Jones potential, Morse potential, Mie potential and Stillinger-Weber (SW) potential and their suitability to reproduce the experimental properties of water. We found that Morse and SW potentials performs better than the rest with SW potential performing well across a wide range of temperatures. Using the Machine Learned Stillinger Weber (MLSW) water model, we performed the contact angle and thin film evaporation studies at surfaces with different wettabilities. The contact angle studies were used to determine the wall-water interaction parameters to simulate hydrophobic and hydrophilic surfaces. The evaporation study simulated the thin film of water evaporating from a surface and condensing to another platinum surface. The surface wettability characteristics of this platinum surface was tuned and evaporation was performed at hydrophobic and hydrophilic surfaces.
The results show a favorable evaporating condition with the hydrophobic surfaces. We believe this study will serve as the first coarse grain water evaporation study and in the future multiscale and mesoscale studies of water at complex geometries can be performed with minimal computational expenditure.