Proceedings of the 13th International Conference on Nanochannels, Microchannels and Minichannels
ICNMM15
July 6-9, 2015, San Francisco, CA, USA

Preprint submitted to ASME

InterPACKICNMM2015-48569
A NEW ALGORITHM FOR CONTACT ANGLE ESTIMATION IN MOLECULAR
DYNAMICS SIMULATIONS
Sumith YD
Department of Mechanical and Aerospace Engg.,
Syracuse University, Syracuse, NY, USA

ABSTRACT
It is important to study contact angle of a liquid on a solid
surface to understand its wetting properties, capillarity and
surface interaction energy. While performing transient
molecular dynamics (MD) simulations it requires calculating
the time evolution of contact angle. This is a tedious effort to
do manually or with image processing algorithms. In this work
we propose a new algorithm to estimate contact angle from MD
simulations directly and in a computationally efficient way.
This algorithm segregates the droplet molecules from the vapor
molecules using Mahalanobis distance (MND) technique. Then
the density is smeared onto a 2D grid using 4th order B-spline
interpolation function. The vapor liquid interface data is
estimated from the grid using density filtering. With the
interface data a circle is fitted using Landau method. The
equation of this circle is solved for obtaining the contact angle.
This procedure is repeated by rotating the droplet about the
vertical axis. We have applied this algorithm to a number of
studies (different potentials and thermostat methods) which
involves the MD simulation of water.

INTRODUCTION
It is quite common in the literature to use Berendsen[1] or
Nose-Hoover[2] thermostat for liquid molecular dynamic
simulations to study different thermodynamic properties of
water. In this paper alongside with the new algorithm, we have
investigated the effect of such thermostats on the contact angle
of water with respect to wall heating algorithm [3].
In the past Bo Shi et al [4] have simulated the contact angle
of water on top of Platinum surface with FCC 111. Their work
was with simulating columbic potential with P3M method [5]
and keeping the temperature constant using Berendsen
thermostat [1]. Maruyama and et al [6] have also done contact
angle studies of water on platinum with truncated potential.

Shalabh C. Maroo*
Department of Mechanical and Aerospace Engg.,
Syracuse University, Syracuse, NY, USA
Email: scmaroo@syr.edu

Both researchers have used ZP potential for water platinum
interaction. There are previous works on contact angle
estimation from MD simulations. Erik et al [7] have come up
with an idea of smearing the molecules into grid with Nearest
Grid Point (NGP) scheme which we will discuss the drawback
later in this paper. In another work Malani et al [8] have shown
a new method to study contact angle from MD. Sergi et al [9]
have calculated contact angle of water from MD using local
averaging and fitting methods. Also there exists a wide range of
studies done on contact angle estimation using MD simulations
[10-13]. Most of these researchers have focused on finding
contact angle and surface energy through different ways.
However almost everyone uses NGP scheme which makes the
prediction of liquid-vapor interface skeptical. Again a
discussion about transient monitoring of contact angle is not
found in the literature.
With the enormous amount of result data estimation of
contact angle using image recognition algorithms or manual
methods becomes almost impossible. There are not many
algorithms which suit well for the direct post processing of
molecular simulations. The paper is divided into two sections,
first explaining the algorithm used for estimating contact angle;
second the details of MD Simulation and models used.

ALGORITHM DETAILS
Often MD simulation results of water and other liquids will
have the information of co-existence of vapor and liquid.
Segregating the liquid droplet or liquid from the vapor will
avoid false detection of Solid-Liquid boundary (LVB) when
using density filtering. Though this step is optional, it would
give more appropriate data for estimating contact angle. For
convenience the whole process is divided into 5 steps.

1

Step 1: Sample preparation
The unwanted vapor molecules (outliers) in the data can be
removed effectively using Mahalanobis Distance [17]
technique.

Fig. 3. Data after removing the vapor molecules using MND

Fig. 1. Equilibrated water droplet on platinum

If 𝑋𝑐 is the 𝑛 × 2 column centered vector consisting of
(𝑥 − 𝑥̅ , 𝑦 − 𝑦̅)data of n points then variance-covariance matrix
𝐶𝑥 is defined as
𝟏

𝑪𝒙 = (𝒏−𝟏) (𝑿𝒄 )𝑻 (𝑿𝒄 )

(1)

Then the Mahalanobis Distance (MND) is calculated as
𝑴𝑵𝑫𝒊 = √𝑿𝒊 𝑪𝒙 −𝟏 𝑿𝒊 𝑻

(2)

Where 𝑋𝑖 is the mean centered data of ith data point. From
this list of MND we can neglect those data points with
considerably high MND values.
The effectiveness of MND technique can be visualized
from Fig. 2 and Fig. 3. The vapor molecules which would
have been a hindrance to estimate the contact angle accurately
was removed easily with MND method.

Step 2: Density grid using B-Splines
In this step a 2D mesh in XZ plane (see figure 2) is
generated and the entire droplet data is projected and smeared
into the grids using 4th order B-spline function. This makes it
easier to find a smooth transition between the liquid core and
the empty space.
The Nearest Grid Point (NGP) method is a traditional first
order method normally researchers use to assign data to the grid
points. Inspired from Hockney and Eastwood [18] the density
of the molecules at every grid point using NGP scheme can be
calculated by the below equation.
𝝆𝒊𝒋 = ∑𝑵
𝒑=𝟏 𝑾 (

|𝒙𝒊 −𝒙𝒑 |
𝒅𝒙

) ∗ 𝑾(

|𝒚𝒊 −𝒚𝒑 |
𝒅𝒚

)

(3)

Where the weight function is defined as
𝟏, 𝒙 < 𝟏𝟐
𝑾(𝒙) = {
𝟎, 𝒐𝒕𝒉𝒆𝒓𝒘𝒊𝒔𝒆

(4)

The B-spline method is inspired from Smooth Particle
Mesh Ewald [19]. Reproducing the definition gives, for any
real number u, let 𝑀2 (𝑢) denote the linear hat function given
by 𝑀2 (𝑢) = 1 − |𝑢 − 1| for 0 ≤ 𝑢 ≤ 2 and 𝑀2 (𝑢) = 0 for
𝑢 < 0 or 𝑢 > 2. For n greater than 2, define 𝑀𝑛 (𝑢) by the
recursion
𝑴𝒏 (𝒖) =

𝒖
𝒏−𝟏

𝑴𝒏−𝟏 (𝒖) +

𝒏−𝒖
𝒏−𝟏

𝑴𝒏−𝟏 (𝒖 − 𝟏)

(5)

For our case n = 4 and u is in fractional coordinates. From Fig.
4 we can see that the B-spline scheme gives an upper hand in
predicting the LVB.

Fig. 2. Center of mass data with vapor molecules
Fig. 4. a) NGP method, b) B-Spline method

2

Step 3: Density based filtering
In this step the exact LVB is calculated based on threshold
density bands. Mathematically this can be explained using
Equations 6 and 7. This will remove the grid points with high
densities which resemble the liquid region and low densities
that resemble the vapor (noise).
𝐓𝐡𝐫𝐞𝐬𝐡𝐨𝐥𝐝𝐌𝐀𝐗 = 𝐰𝟏 ∗ 𝐌𝐞𝐚𝐧

(6)

𝐓𝐡𝐫𝐞𝐬𝐡𝐨𝐥𝐝𝐌𝐈𝐍 = 𝐰𝟐 ∗ 𝐌𝐞𝐚𝐧

(7)

𝑤1 and w2 are weights which can be fine-tuned according
to the data. For the present studies we have taken them as 2 and
0.5 respectively. Now, all the data which meets the below
criteria will be used for further analysis.
𝐓𝐡𝐫𝐞𝐬𝐡𝐨𝐥𝐝𝐌𝐈𝐍 ≤ 𝐆𝐫𝐢𝐝. 𝐯𝐚𝐥𝐮𝐞 ≤ 𝐓𝐡𝐫𝐞𝐬𝐡𝐨𝐥𝐝𝐌𝐀𝐗

Fig. 6. Radial rays from centroid to find the extremities

(8)

This will lead to detection of the points as shown in Fig. 5.

Fig. 7. Detected droplet boundary
Fig. 5. Density based filtered data

Step 4: Removing the wall effects
As one may notice, an unwanted lower set of contour is
also detected due to the interface between water and platinum.
This has to be taken care as explained next. Find the centroid of
the data from the Step 3 and set the z coordinate as the top of
the platinum plate. Then we radially move outwards from that
point until we reach the boundary of data. This procedure is
repeated for 0 ≤ 𝜃 ≤ 𝜋. This will help us to identify the most
appropriate data. At the end we will be having a collection of
points which forms the exterior of the droplet.
Fig. 6 shows this procedure graphically and Fig. 7 shows
the detected droplet boundary. This will also eliminate the
unwanted data points inside the droplet due to the occasional
void formation due to the density fluctuations in the MD
simulations. Also using a z-location based cutoff we can
eliminate the unwanted monolayer or interface formed next to
the platinum wall.

If the centroid of the data was above the platinum surface
then we should leave the z –coordinate of centroid as it is and
angular search should be from 0 ≤ 𝜃 ≤ 2𝜋. This will guarantee
the contact angle estimation of droplet with contact angle
greater than 90°. Also this is required if the droplet is floating in
vacuum and not attached to the solid.

Step 5: Solving for contact angle
With the detected LVB step 4, we can fit a circle using
Landau method [20] and get the circle’s equation which
represents the droplet boundary. Landau method is an efficient
method to fit a circle using non iterative geometric fit. This
method relies on minimizing the error of fit between the set of
points and the estimated arc. The main equations for this
method are mentioned below. For a detailed version and
nomenclature we encourage readers to see the original
reference.

𝒙
̅=

3

𝒄𝟏 𝒃𝟐 −𝒄𝟐 𝒃𝟏
𝒂𝟏 𝒃𝟐 −𝒂𝟐 𝒃𝟏

(9)

̅=
𝒚

𝒂𝟏 𝒄𝟐 −𝒂𝟐 𝒄𝟏
𝒂𝟏 𝒃𝟐 −𝒂𝟐 𝒃𝟏
𝟏

̅ + 𝑵𝒙
̅𝟐 + ∑ 𝒚𝟐 − 𝟐 ∑ 𝒚
̅ + 𝑵𝒚
̅𝟐 ]
𝑹𝟐 = [∑ 𝒙𝟐 − 𝟐 ∑ 𝒙
𝑵

(10)

(11)

Once we get the equation of circle we can solve it for the
angle made by the tangent at the platinum-water interface. This
step is graphically shown in Fig. 8.
During a fresh simulation it is advisable to do a test case
and fine tune the parameters to make this algorithm well suited
for the selected data. These parameters are including, but not
limited to the weight functions of step 3, z – cutoff of step 4.

platinum plate as shown in Fig. 9. The lateral boundaries are
applied with periodic boundary condition (PBC) and vertical
boundary is closed with platinum wall on bottom and mirror
boundary on top. The platinum wall is of square shape with 20
nm sides and mirror boundary is kept at a distance of 12.5 nm
from platinum. The water-platinum interaction is simulated
using ZP-potential. A self-written and validated C++ code is
used to perform the analysis. There are two sub cases of
simulations based on this model. They are Berendsen and
surface heating algorithm. More details of this simulation are
described in the following paragraphs.

Fig. 9. Simulation model for contact angle estimation

Fig. 8. Final boundary of the droplet and fitted circle

CONTACT ANGLE ESTIMATION STUDIES

Case 1: Water-platinum interaction using LJ
potential
A cube shaped water droplet with 6 nm sides (7221
molecules) is kept on top of platinum wall interacting using LJpotential. The sides are under periodic boundary condition and
top is covered with another platinum surface to prevent the
escape of molecules. The simulations are performed using
Gromacs software. The platinum surfaces are modeled with
FCC 111 structure and are square shaped with sides of 25 nm.
The platinum surfaces are kept apart at a distance of 14 nm.
There are three sub cases for this model based on the thermostat
used. They are Berendsen, Nose-Hoover and Velocity re-scale
thermostats. The simulation is performed from 0 ps to 1000 ps.

Case 2: Water-platinum interaction using ZPpotential
In this case a pre-equilibrated cubic shaped water droplet
with sides of 4 nm (2133 molecules) is kept on top of FCC 111

The equations of motions are solved using velocity
verlet scheme. A shifted scheme suggested by Stoddard and
Ford [14] with 7σ was used instead of Ewald summation
methods or P3M. The feature of this potential is that both
potential and force goes smoothly to zero at the cut off radius.
12
6

  
  
U  r   4      
r 
 r 





   12   6   r 2   12   6 

 6   - 3      - 7    4   
 rc    rc 
 rc 
 rc  
  rc 


(12)

The water molecules are modeled using SPCE model
introduced by Berendsen et al [14]. The intra molecular bonds
are kept rigid throughout the simulation using the RATTLE
algorithm [15].
Zhu-Philpott [16] model is used to model the interaction
between water and platinum molecules. The advantage of such
model is to obtain reasonable droplet shape formation during
the simulation [4].

ETotal  EConduction  EIsotropic  EAnisotropic

EConduction 


real ,image

4

(13)

qreal qimage
2r

(14)

Cw Pt w10 Pt
rwi10
i 1

N Pt

EIsotropic  4 w Pt 

6
3

2
2
 



w

Pt
w

Pt
 4 w Pt  
 

2
2 
    2  z 2  
i 1  

z


wi 
wi 
 wi
 wi


(15)

N Pt

E Anisotropic

(16)

The parameters are taken from Maruyama’s paper [6].
For simulations with thermostat based heating we have
used Berendsen thermostat [1]. For those with wall heating we
have used the modified Maroo-Chung model [3]. The details of
this surface heating algorithm will be discussed in a separate
paper which is yet to be published.
Since we are not modeling a contiguous array of droplets
we can ignore the long range effects of the single water droplet.
Simulation of droplet and films with SPME method will be
studied extensively and will be published in another work.
Hence we limit our studies with truncated potentials for the
present work.
For the sub cases we have equilibrated the system at 300K
for 200 ps. From 200 ps to 1000 ps the Berendsen thermostat or
surface heating is applied for the two sub cases. For the last sub
case we left the system uncoupled from thermostat from 200 ps
to 1000 ps.
At the end of simulation the trajectory files are processed
and contact angles are determined using the algorithm
mentioned in the next section. The time evolution of contact
angle is shown in Fig. 11 and Fig. 10.

Fig. 10. Contact angle evolution for case 1

Fig. 11. Contact angle evolution for case 2

Since the main focus of this paper is not the effect and
causes of the contact angle we have not included those
discussions here. However we have made such a detailed study
and currently that work is under review elsewhere. A working
version of MATLAB code for this algorithm can be obtained
from authors upon request.
CONCLUSION
A new algorithm for contact angle estimation from
transient data of MD simulation was discussed. The efficiency
and accuracy of this algorithm was demonstrated using
different types of simulations. The MD simulations were done
for water platinum interaction with ZP potential and LJ
potential. Further, we used different thermostatting schemes for
the droplet and show how it influences the contact angle. The
results also showed the influence of interaction potential and
also the thermostatting method used. The algorithm can be used
to study the transient behavior of droplets on different surfaces.

REFERENCES
[1] Berendsen, H. J., Postma, J. P. M., van Gunsteren, W. F.,
DiNola, A., and Haak, J., 1984, "Molecular dynamics with
coupling to an external bath," The Journal of chemical
physics, 81(8), pp. 3684-3690.
[2] Hoover, W. G., 1985, "Canonical dynamics: equilibrium
phase-space distributions," Physical Review A, 31(3), p.
1695.
[3] Maroo, S. C., and Chung, J., 2010, "A novel fluid–wall heat
transfer model for molecular dynamics simulations," Journal
of Nanoparticle Research, 12(5), pp. 1913-1924.
[4] Shi, B., Sinha, S., and Dhir, V. K., "Molecular simulation of
the contact angle of water droplet on a platinum surface,"
Proc. ASME 2005 International Mechanical Engineering
Congress and Exposition, American Society of Mechanical
Engineers, pp. 93-97.
[5] Eastwood, J., Hockney, R., and Lawrence, D., 1980,
"P3M3DP—The three-dimensional periodic particle-

5

particle/particle-mesh program,"
Computer
Physics
Communications, 19(2), pp. 215-261.
[6] Kimura, T., and Maruyama, S., 2002, "Molecular Dynamics
Simulation of Water Droplet in contact with a Platinum
Surface," Heat Transfer, 1, pp. 537-542.
[7] Santiso, E. E., Herdes, C., and Müller, E. A., 2013, "On the
Calculation of Solid-Fluid Contact Angles from Molecular
Dynamics," Entropy, 15(9), pp. 3734-3745.
[8] Malani, A., Raghavanpillai, A., Wysong, E. B., and
Rutledge, G. C., 2012, "Can Dynamic Contact Angle Be
Measured Using Molecular Modeling?," Phys. Rev. Lett.,
109(18), p. 5.
[9] Sergi, D., Scocchi, G., and Ortona, A., 2012, "Molecular
dynamics simulations of the contact angle between water
droplets and graphite surfaces," Fluid Phase Equilibria, 332,
pp. 173-177.
[10] Ohler, B., and Langel, W., "Molecular Dynamics
Simulations on the Interface between Titanium Dioxide and
Water Droplets: A New Model for the Contact Angle,"
Journal of physical chemistry. C, 113(23), pp. 10189-10197.
[11] Blake, T. D., Clarke, A., DeConinck, J., and deRuijter, M.
J., 1997, "Contact angle relaxation during droplet spreading:
Comparison between molecular kinetic theory and
molecular dynamics," Langmuir, 13(7), pp. 2164-2166.
[12] Nijmeijer, M. J. P., Bruin, C., Bakker, A. F., and
Vanleeuwen, J. M. J., 1989, "A VISUAL MEASUREMENT
OF CONTACT ANGLES IN A MOLECULARDYNAMICS SIMULATION," Physica A, 160(2), pp. 166180.
[13] Hong, S., Ha, M. Y., and Balachandar, S., 2009, "Static and
dynamic contact angles of water droplet on a solid surface
using molecular dynamics simulation," Journal of colloid
and interface science, 339(1), pp. 187-195.
[14] Berendsen, H., Grigera, J., and Straatsma, T., 1987, "The
missing term in effective pair potentials," Journal of
Physical Chemistry, 91(24), pp. 6269-6271.
[15] Andersen, H. C., 1983, "RATTLE: A “Velocity” version of
the SHAKE algorithm for molecular dynamics
calculations," Journal of Computational Physics, 52(1), pp.
24-34.
[16] Zhu, S.-B., and Philpott, M. R., 1994, "Interaction of water
with metal surfaces," DTIC Document.
[17] De Maesschalck, R., Jouan-Rimbaud, D., and Massart, D.
L., 2000, "The mahalanobis distance," Chemometrics and
intelligent laboratory systems, 50(1), pp. 1-18.
[18] Hockney, R. W., and Eastwood, J. W., 2010, Computer
simulation using particles, CRC Press.
[19] Essmann, U., Perera, L., Berkowitz, M. L., Darden, T.,
Lee, H., and Pedersen, L. G., 1995, "A smooth particle mesh
Ewald method," The Journal of chemical physics, 103(19),
pp. 8577-8593.
[20] Thomas, S. M., and Chan, Y., 1989, "A simple approach
for the estimation of circular arc center and its radius,"
Computer Vision, Graphics, and Image Processing, 45(3),
pp. 362-370.

6

