Achilles: A Tool for Contact Angle Estimation from Molecular Dynamics Simulations

In this work, a tool for estimating the contact angle from the molecular dynamics simulations is developed and presented. The tool (Achilles) can detect water droplet on hydrophobic and hydrophilic surfaces. The tool can reconstruct the droplets broken across the periodic boundaries. Further a neighbor density based accurate filter is used to find the droplet liquid vapor interface and a circle is fitted using it after removing the dense layers of water next to solid surface. This fitted circle is solved for contact angle and results are outputted in the form of graphical images and text. The entire content of the internal computations of the tool is broken down into 4 phases and users can monitor the outcomes at every phase through output images. The tool is tested using sample molecular dynamics results of water droplet on hydrophobic and hydrophilic surfaces. We believe this tool can be a good addition to the molecular dynamics simulation community who work on the interfacial physics, droplet evaporation, super hydrophobic surfaces, and wettability etc.


Introduction
Contact angle measurement plays an important role in determining the surface wettability of nano coatings, surface modifications and natural surfaces. Some of the recent studies include characterization method of intermediate wetting states, which determines the condition of anisotropic sliding on micro-grooves [1], Electromagnetically induced change in interfacial tension and contact angle of oil droplet using dielectric nanofluids [2], effects of hysteresis window on contact angle hysteresis behavior at large Bond number [3], Droplet behaviors on inclined surfaces with dynamic contact angle [4], estimating critical surface tension from droplet spreading area [5], geomaterial wettability and contact area versus three-phase contact line [6], contact angle measurements using sessile drop on six sandstones [7], analysis of the contact angle for impacting droplets using polynomial fitting [8] etc. The molecular dynamics simulations are very useful in studying nanoscale phenomena [9]- [11].
There exist many contact angle studies using molecular dynamics tool at nanoscale level. These studies aim to investigate the contact angle hysteresis of droplets on surfaces, contact angle relaxation, contact angle of soil minerals, graphite surface wettability, static an dynamic contact angle, droplets in carbon nanotubes, drop size influence and surface potential influence on contact angle and temperature dependence of contact angle [12]- [21] etc. Also, the molecular dynamics and multiscale studies are helpful in understanding the fundamental behavior of many systems and phenomena [22]- [40]. Not many dedicated tools are available for estimating the contact angles from the molecular dynamics (MD) simulations [41]- [43]. The estimation of the contact angle from the MD simulation is often challenging due to the presence of dense layers adjacent to solid surface, estimation of liquid vapor interfacial line and movement of droplet across the periodic boundary. For these reasons, most of the studies perform preconditioning of the data (mostly coordinates of the water molecule) after the computer simulations and then estimate the contact angle manually using graphical methods. Some researchers estimated contact angle based on the eigenvector of the covariance matrix of the interfacial molecules [20]. Many other studies use density based image processing to estimate the contact angle [44]- [48]. When we study large number of data sets, automated contact angle estimation will be a helpful approach. With this in mind, we have developed Achilles tool for contact angle detection from molecular dynamics simulations.
The Achilles tool will read the coordinates of the water and the surface plate from a the given xyz file. Based on the conditions and parameters set in the configuration file, the Achilles tool will find the location of effective contact between the water droplet and the solid surface. The band filter based on a threshold neighbor atom density will extract the atoms/beads in the interface of liquid and vapor. These atoms/beads are used to fit a circle which in turn is solved for the contact angle. The tool performs well in situations where the droplet spreads across the periodic boundaries, have dense layers which are not to be included in the contact angle estimation process, etc. The tool also gives useful information like coordinates of the droplet, filtered data, center and radius of the fitted circle to the users at regular intervals if needed through images and in text form.

Theory and Methods
The water droplet can rest on surfaces with wettability ranging from super hydrophilic to super hydrophobic. This means, the contact angle ( ) can vary between 0 deg to 180 deg. Specifically, surfaces are classified based on the wettability as superhydrophobic ( > 150°), hydrophobic ( > 90°), and hydrophilic ( < 90°) [49]. To account for these cases, the software (henceforth in this paper, this refers to our tool -Achilles) will classify the data into two categories -1) hydrophilic and 2) hydrophobic. The first case of hydrophilic surfaces refers to the Figure 1, where the water droplet rests on a flat surface, let us say platinum surface [41], [49]- [51]. Depending on the presence of the strong attractive forces near the surface, highly dense layer(s) can form closer to it [51]. This layer(s) is often disregarded in contact angle estimation [52]. Referring to the Figure 1, is the vertical distance of the platinum or any other solid surface from the origin. The variable indicates the distance from the plate surface at which the effect of the dense layer formation vanishes.
In the case of hydrophilic surfaces, the droplet will form like as shown in the Figure 1. The point of contact of three phases (solid, liquid and vapor) is indicated with coordinates ( , ) , shown in red dot. The center of the fitted circle is indicated as ( , ). The surface of the is shown using blue solid line. The tangent of the circle at the three-phase contact point is shown in green. The contact angle is then shown as and can be estimated from the below equation. For hydrophobic surfaces, the droplet can be considered as shown in the Figure 2. This time, the center of the circle is above the dense layer. The contact angle corresponding to this setup can be estimated as shown in the below equation.

Software development and architecture
We have chosen C sharp language to develop the tool due to its low production time and relatively better performance in comparison with Python. However, this claim is highly debated in online forums in relation with the optimization, use of precompiled libraries, use of operating system (OS) specific libraries, CPU architecture specific optimizations and so on. A thorough investigation and opinion on this context is not the focus of this paper.
To reduce further overhead in execution of the program, the graphical user interface is dropped and instead we will be using a command line interface. The interaction with the Achilles tool will be through two files 1) config.txt and 2) pathway.txt. A screenshot of the contents of the config.txt is shown in the Figure 3. The contents of this file are case sensitive and follows the sequence "variable = value". The variable plate indicates platinum or any other solid surface, and its value indicates the atom type included in the "xyz" file. The variable water indicates the water droplet, and its value indicates the atom type in the "xyz" file. The variable debug can be set to 1, if the user wants to get detailed information about every phase of contact angle estimation by the tool. Variable walleffect is the value of the dense layer thickness indicated in the Figure 1, in Angstrom units. If the user does not want to exclude the dense layers, then simply set this variable to 0. The remaining variables will be discussed in detail later in this section of the paper. The pathway.txt file contains the list of the xyz files to be processed for estimating the contact angle. A sample file content of this file is shown in the Figure 4. Currently, in this first version of the contact angle tool, only xyz files with a single timestep is supported. This means, if you have a simulation running for 5000 steps, writing output at every 50 steps, and needed contact angle at every this many intervals, then you have to create individual xyz files at those intervals (in this case 5000/50 = 100 files). The format of the xyz file should be according to the standard xyz representation [53]. The Achilles tool will output all the important information to the a_logfile.txt for results and error checking. The sample contents of the logfile and the Achilles command line interface is shown in the Figure 5. This information includes the shifted values of the system, center of mass, density per atom etc. which will be discussed in the following paragraphs. For easier handling of the development process and debugging the tool, the Achilles software is logically divided into four phases as shown in the Figure 6. The contents of the different phases are discussed below.

Phase 1:
This phase starts with reading the config.txt and pathway.txt files. Then based on the contents of the pathway.txt file, the program finds the number of atoms present in it. Based on this number, memory is allocated for storing the x, y, and z coordinates of the atoms. The coordinates of the atoms are then read into it. Since xyz files do not carry information on mass, boundaries etc., they must be estimated or assumed. Achilles tool estimates the boundaries of the system by doing multi-pass reading. This step will yield xlo, xhi, ylo, yhi, zlo, zhi which are the lower and upper bounds along 3 dimensions. Based on the information of config.txt, the atom type of plate will be known and from the xyz file the will be estimated.

Phase 2:
On many occasions, molecular dynamics simulations can lead to droplet crossing across the periodic boundary. This will be automatically addressed by the Achilles tool. First, the center of mass of the data across the periodic boundary system can be estimated based on the method found in literature [54]. Once the center of mass of the atoms are found, then simple translation can reconstruct the droplet in place. To make further computations easy, the droplet is shifted to the center of the system and adjustments are made to make xlo = ylo = zlo = 0.

Phase 3:
In the third phase of the processing, based on the brute force algorithm [55] the number of neighbors of every atom within a cutoff radius of is estimated. The value of this is supplied through config.txt file. Using a quick search, the maximum neighbor value is estimated. Next, atoms are filtered using a threshold range between ℎ ℎ × ℎ and ℎ ℎ ℎ × ℎ . Here, ℎ ℎ and ℎ ℎ ℎ values are specified through config.txt file. The defaults values should work for most cases. However, if the user suspects it, then turn on the debug variable and Achilles will plot the filtered coordinates as a png file and as a xyz file. These xyz files can be opened using any standard molecular visualization software and can inspect its accuracy. If needed, the threshold values can be adjusted until it meets the user requirements.

Phase 4:
The effect of the dense layers of water is removed simply by using a distance cutoff filter. Recall the contact angle estimation equations discussed in the previous (Theory and Methods) section. It needs classification based on hydrophobicity. To achieve it, the filtered coordinate data is used to construct a circle using a fast sphere fit algorithm [50]. This circle's center ( , ) and can determine whether the droplet is hydrophobic or hydrophilic (refer Figure 1and  Using this information and from the equation 1 and 2, the contact angle will be estimated. Additionally, the Achilles tool will output the result in graphical form to a png file. This includes display of the x, y, z coordinate data of atoms, fitted circle, tangent line at contact point, location as a line. Also, title of the file contains the contact angle in degrees, coordinates of center of the circle, contact points etc.

Results and Discussion
To test the Achilles tool, three different cases are considered -1) super hydrophobic surface, 2) hydrophobic surface and 3) hydrophilic surface. An equilibrated molecular system with 1000 coarse grained molecular dynamics (CGMD) water beads is used on the top of a platinum surface with an FCC 111 structure with an edge length of 3.912 Ang [56]. The beads follow Morse potential [57] and represents 4000 water molecules based on the CGMD water bead force fields [49], [58]. The platinum plate is modeled with dimensions of 15 nm by 15 nm sides and a single layer. The screenshot of the equilibrated droplet on platinum is shown in the Figure 7. For our 3 cases, the time step of integration is kept as 10 fs, equilibration time of 50,000 steps, followed by a production time of another 50,000 steps. The temperature is maintained at 298 K using a Nose Hoover thermostat [59], [60]. The Morse potential takes the form as below.
Here, is the potential energy in kcal/mol, 0 is the potential well depth, 0 is the equilibrium distance at which the pair potential becomes its minimum, is the spread controlling factor of the potential. Also, r is the distance between atoms or beads and for practical computational purposes, the inter atomic or inter CGMD bead interactions are limited by a cut off distance of .
The values of the potential parameters are given by the Table 1. The Platinum plate is modeled as a single flat layer and kept rigid throughout the simulations. The interaction between the water beads and platinum atoms are done using a Lennard Jones potential ( ) as given below.
Here, is the inter bead distance at which the potential becomes zero, is the inter bead/atom distance in Angstrom and is the strength parameter in kcal/mol. For our studies, is considered as 4.7 nm from the standard MARTINI model [61]. The strength parameter can be varied to change the wettability of the surface. We conducted a series of simulations with epsilon parameter ranging from 0.04 to 1.0. From these simulations, three selected cases are selected and its contact angle estimation results from Achilles tool is shown in the coming sections.

Case 1: Super Hydrophobic Surface ( = . )
For this case, value is 0.04 and the surface interaction with the water beads became superhydrophobic. The debug option is turned on while running the Achilles tool and results at the end of various phases are given in the Figure 8 (phase 1), Figure 9 (phase 2), Figure 10     In Figure 12, the contact point is shown outside the circle, tangent line and flat plate becoming colinear. This is because Achilles tool could not find any solution for contact angle yet want to give a visual representation of the data and fitted circle. The remaining details of the estimation process are shown as the title of the Figure 12. This includes contact angle in degrees, coordinates of the circle, and coordinates of the contact point. The units will be same as the xyz file which is Angstrom.

Case 2: Hydrophobic Surface ( = . )
For this case, value is 0.2 and the surface interaction with the water beads became hydrophobic. The debug option is turned on while running the Achilles tool and results at the end of various phases are given in the Figure 13 (phase 1), Figure 14 (phase 2), Figure 15            These 3 case studies, along with the visual results shown in here, the Achilles tool does a good job in finding the accurate contact angle.

Achilles Heel
The tool that we developed serves the purpose of estimating the contact angle. However, the author believe it can be further improvised and below are some of the ideas.
• Processing of the xyz files with multiple timesteps. • Development of the documentation and examples for getting started with beginners.
• Graphical user interface (GUI) for better user experience. • Migration to C++ language for faster execution and for using widely available standard libraries. • Multithreaded processing of the data for faster handling of large data. • Pipeline processing of the data while reading very large files. • Extending to multiple droplet detection.
• Extending the software for bubble detection and boiling. • Extend the tool for cross platform deployment (Mac, Linux, and Windows)

Conclusion
A computational tool for estimating the contact angle of water droplets from the molecular dynamics simulations is developed and presented in this work. The tool is named as Achilles and can read the single frame information of xyz files and estimates the contact angle. The periodic boundary wrapping of the droplet is carried out and reconstructs the accurate droplet from it. Further a neighbor density based accurate filter is used to find the droplet liquid vapor interface and a circle is fitted using it after removing the dense layers of water next to solid surface. This circle is solved with the contact point of the surface to get the contact angle of the droplet. Sample molecular dynamics simulations are performed for a system with water droplet on a platinum plate and Achilles tool is used to estimated 3 different cases of wettability. The results are shown graphically and through a logfile text which shows accurate estimation of contact angle. The tool is highly customizable by the user and can monitor the data processing at 4 different phases through screenshots of the data. We believe this tool will benefit a larger growing community of interfacial studies using molecular simulations and is freely available to download and use.

Appendix 1: Hosting of the software
The software is available as free to use and is hosted permanently at below locations. In case of any issues in installing and using the software, please send an email to the author. Currently the tool supports only windows-based computers.