Pytim: A python package for the interfacial analysis of molecular simulations.

Pytim is a versatile python framework for the analysis of interfacial properties in molecular simulations. The code implements several algorithms for the identification of instantaneous interfaces of arbitrary shape, and analysis tools written specifically for the study of interfacial properties, such as intrinsic profiles. The code is written in the python language, and makes use of the numpy and scipy packages to deliver high computational performances. Pytim relies on the MDAnalysis library to analyze the trajectory file formats of popular simulation packages such as gromacs, charmm, namd, lammps or Amber, and can be used to steer OpenMM simulations. Pytim can write information about surfaces and surface atomic layers to vtk, cube, and pdb files for easy visualization. The classes of Pytim can be easily customized and extended to include new interfacial algorithms or analysis tools. The code is available as open source and is free of charge.

Introduction

The interface between two fluids is not a well‐defined concept as soon as one leaves the continuum description and reaches the molecular detail level.1 The determination of atoms belonging to an interface is intrinsically dependent on the scale at which the observer probes the interface itself. In addition, thermally activated capillary waves, by corrugating the interfaces, are adding further complexity to the determination of the physicochemical properties of fluids in these boundary regions.

Historically, several algorithms were devoted to the calculation of the accessible surface (to the solvent or other small molecules) of macromolecules,2 and several approaches were developed to compute what is now known as the solvent accsssible surface area (SASA)3, 4 or the connolly surface.5, 6 More recently, the field has been stirred up by the necessity of investigating the structure of fluid/fluid interfaces from an intrinsic point of view, namely, by removing the smearing effect of thermal capillary waves,7 and several methods focused on liquid/liquid interfaces were developed, such as Willard and Chandler instantaneous liquid interfaces method,8 the circular variance method of Mezei,9 the method of Škvor and colleagues10 or the method of Yesylevskyy and Ramseyer for the calculation of the local curvature.11

The Pytim package includes a series of algorithms to identify phases and surfaces or surface atoms, and perform various types of analyses related to interfaces. The main algorithms implemented in Pytim are: (1) the identification of truly interfacial molecules12 (itim) for macroscopically planar systems; (2) the generalized itim 13 (gitim), for arbitrarily shaped systems; (3) the algorithm of Willard and Chandler,8 to compute continuous surfaces based on a smoothed estimate of the atomic volumetric density; (4) the Lee‐Richards algorithm3 to compute the SASA and the associated solvent‐exposed molecules; (5) an efficient cutoff‐based clustering algorithm to distinguish liquid and vapor phases,14 and (6) an improved, density based clustering algorithm for highly miscible systems.15

All algorithms that identify surface atoms can be used not only to compute the interfacial layer, but also the successive layers beneath it. The algorithms and the peculiarities related to their implementation in Pytim are briefly reviewed in the next section.

The underlying philosophy of Pytim is to provide a flexible, extendable, and easily scriptable system for the calculation of interfaces and interfacial properties in molecular systems that free the user from the burden of writing ad hoc code for different simulation packages. Pytim is built on top, and extends, the MDAnalysis library,16 which provides the backend for reading different trajectory formats. All algorithms are implemented by making use of the efficient numerical libraries wrapped by numpy and scipy that take advantage of the SIMD architecture of modern processors. Multiprocess parallelization has been implemented, when possible and Cython has been used to speed up inner loops. The resulting code has a speed comparable to that of a compiled one. The user interface has been designed to be intuitive and straightforward, in the spirit of the python language, and allows to perform the fundamental analyses using just a couple of scripting lines. This task is also facilitated by providing meaningful default values for all options. However, the user has complete control over the details of the algorithms and can quickly perform complex tasks by combining the basic building blocks of Pytim. In particular, the structure of the interface is meant to encourage the user to experiment with the possibilities offered by combining Pytim with other libraries and create their own analysis tools. MDAnalysis has been chosen for its high level of abstraction, which allowed to reduce development time and simplicity in implementing other backends. In particular, the pytim code can be used in a transparent way with MDTraj,17 and compute interfacial properties online during OpenMM18 simulation runs.

Methods

Algorithms to compute interfaces and interfacial atoms

ITIM

This algorithm12 can be loosely described as a molecular version of the popular pinscreen toy. It consists in determining the atoms in contact with a virtual surface made of probe spheres bound to move perpendicularly to the surface plane along test lines. The present implementation follows that of Jorge et al.,19 which is based on sorting atoms according to their distance from the center of mass of the slab. Then, starting from the furthermost atom, one finds which test lines allow a sphere of radius α, which moves along the line, to touch the atom with a radius . Here, σ corresponds typically to the Lennard‐Jones parameter.20 The sorting of N atoms can be performed in steps, while finding which of the test lines are within a given radius requires steps using a kd‐tree.21 As the test lines are usually 10 times more than the surface atoms (for the algorithm to work accurately), one has typically , and in the worst case . Therefore, the algorithm scales globally like . By assuming convergence at a large number of test lines, this method is left with one free parameter, the probe sphere radius α, which sets the scale at which the interface is probed.

Figure 1 reports an example of how to use the itim algorithm in Pytim.

Figure 1

Usage of the class to identify the first four molecular layers in a water/vapor interface. The image has been generated in a jupyter notebook using nglview. [Color figure can be viewed at http://wileyonlinelibrary.com]

GITIM

In this generalization of the itim algorithm, no macroscopic orientation is assumed. The touching spheres, instead of “raining down” along probe lines, can be thought of as being “inflated” at all points in space, up to their maximum radius α (in which case the surrounding atoms are considered to be interfacial ones) or until they start touching neighboring atoms. This procedure is realized in gitim using a modification of the alpha‐shapes algorithm22 that takes into account the excluded volume of the atoms,13 and implements in its essence a filter of the Delaunay23 triangulation of atomic positions. In other words, all the triangular faces belonging to those tetrahedra, whose inner touching sphere has a radius larger than α, are considered to be interfacial ones, and so are considered the atoms at the vertices of those triangular faces. The present implementation of gitim can calculate the Delaunay triangulation using the scipy library, which uses, in turn, a wrapper to the Quickhull24 implementation Qhull. In addition, it can use the faster pytetgen python wrapper to the tetgen software.25 Both algorithms scale on average as and in the worst case as , setting, therefore, the scaling of gitim. Also in this case the only free parameter is the probe sphere radius α.

Figure 2 reports an example of how to use the gitim algorithm in Pytim.

Figure 2

Usage of the GITIM class to identify three solvation shells of a glucose molecule. The image is a section cut of the system, as visualized using nglview. [Color figure can be viewed at http://wileyonlinelibrary.com]

Willard and Chandler's instantaneous liquid interface

This algorithm8 defines a continuous interface as the isodensity surface of a Gaussian kernel density estimate.26, 27 A continuous Gaussian density function, with given height h and width w , is associated with each atom. Summing over the contributions of all atoms, one obtains a continuous density field, which is usually sampled on a regular grid. The isosurface where the density is in the neighborhood of a given, target density (typically half between the maximum and the minimum) is then used as the definition of the interface itself. In the present implementation, we use both an exact method, by deriving a class from scipy.stat.gaussian_kde that implements periodic boundary conditions, and an approximated one that calculates the Gaussian contribution from points closer than 2.5. We make use of the topologically consistent version28 of the marching cubes algorithm29 as implemented in scikit‐image 30 to extract the isosurface from the kernel density estimate. The computation of the kernel density estimate scales like . Here, N is either the number of grid points in the simulation box (in the exact implementation) or the number of grid points within the sphere of radius 2.5 (in the approximated implementation). The marching cube algorithm scales like , so that is the global scaling. In the most common scenario, the grid spacing is kept constant, independently on the system size, meaning that the computational complexity of the exact implementation scales like , while the approximated one scales like .

The algorithm has in principle three free parameters, namely w , h , and the target density. However, by choosing the target density as half of the maximum one, h becomes irrelevant. The algorithm does not directly determine surface atoms, but one could define them as those within a certain distance from the isodensity surface. The Willard–Chandler algorithm is not restricted to planar surfaces.

Lee‐Richards SASA algorithm

The algorithm first creates, for each of the atoms, a list of neighbors within a cutoff radius, which is twice the sum of the largest atom radius and the probe sphere radius. The spherical volume associated to each atom that has as radius the sum of the atomic and of the probe sphere radii is then sliced in several discs along a chosen direction. Using simple geometrical considerations, one can calculate the total arc length of each disc, which happens to overlap with the corresponding discs associated to neighboring atoms. In this way, it is possible to compute an approximation of the area that is (or is not) accessible by a probe sphere. Atoms with nonzero accessible surface area are considered to belong to the interfacial layer. As the neighbor search algorithm scales like , and the calculation of the exposed area scales like for each atom, the global scaling is .

Algorithms to determine phases

To determine the interface between two phases containing the same chemical compound, it is, of course, necessary to assign which atoms belong to that phase. This task can be in principle trivial for liquid/vapor interfaces far from the critical temperature, as typical simulated system sizes are small, and often the vapor phase is just empty. For liquid/liquid interfaces, the equivalent condition is that of a perfectly demixing system. In these cases, it is in principle not necessary to perform any particular prefiltering of the particles before proceeding to the identification of surfaces or surface molecules using any of the above algorithms. One has to be careful, however, that the presence of a single molecule either in the vapor phase or solvated in the opposite phase, can jeopardize the whole interfacial determination. This problem is, of course, even more serious when getting closer to the critical temperature.

One of the most straightforward ways to separate molecules in the vapor phase from those in the liquid is to group all atoms in clusters and consider the biggest cluster as the liquid phase.

Simple clustering

A simple yet effective strategy is to consider molecules to be in the same cluster if any of their atoms are closer than a chosen cutoff distance.14 Usually, a reasonable choice for the cutoff value is the position of the minimum after the first peak of the pair distribution function in the bulk liquid. This way, one is sure to include in the cluster all atoms belonging to the first solvation shell. In this case, the liquid fraction will form a percolating network of connected molecules, whereas the molecules in the vapor phase will be part of smaller, separated clusters.

DBSCAN

The simple cutoff clustering scheme is limited to the case of low temperatures and low miscibilities. At high enough temperatures, for example, the density of the vapor phase can become high that the whole system percolates using any choice of cutoff. In this case, it is more effective to look at the local density rather than at the connectivity. The DBSCAN density‐based cluster algorithm31 is well‐suited for this task, and has been already used to study supercritical fluids32 and a two phase system of highly miscible fluids.15 The main idea of the algorithm is to include in a cluster all atoms with overlapping neighborhoods where the density is higher than a given threshold. In Pytim, we implemented an approach for the automatic determination of the density threshold, based on a k‐mean analysis33, 34 of the density distribution.15

Algorithms to compute interfacial properties

Once the surface atoms are determined, these are available through the python interface for further analysis. In the Pytim package, we provide the user with an extendable set of classes for the computation of observables, suited both for the calculation of bulk properties, and interfacial ones. In particular, one can combine each observable with the calculation of intrinsic and non‐intrinsic density profiles, as well as pair distribution and time correlation functions of atoms within the interfacial layers. The (non‐intrinsic) profile of observable O across a macroscopically planar interface is computed in the simulation box‐fixed reference frame aswhere 〈…〉 represents an ensemble average, is the Dirac delta function, and the sum is extended over the set of atoms of interest. This corresponds to the ordinary notion of a density profile, where the bin width is given by . Here the discretized version of the Dirac delta function is if |.

The intrinsic profile, on the contrary, is obtained by referring the position of the particles to the surface capillary waves, obtaining in this way the local, interfacial structure of the fluid without the smearing effect of the capillary waves themselves. If a particle is located at in the box reference frame, we denote the position of the local reference frame on the corrugated surface as (. In this way, the intrinsic profile can be written as7, 35

Notice that z can be either negative or positive, depending on the location with respect to the interface, and can be, therefore, thought of as a signed distance.

This concept can be generalized to nonplanar macroscopic interfaces by choosing the distance of an atom i located at from the surface as its distance to the closest surface atom, j. However, care has to be taken in determining the bin width and its sign. Although there is, in the general case, no analytical expression for the area of the iso‐distance surfaces, we estimate the bin volume numerically, by using a simple Monte Carlo scheme.13 The second point one has to take care of is the determination of the sign of the distance that determines whether a point is located below or above the interface. While in the planar case this is straightforward, in the general case, we resort to the following. We compute the center of mass of the local environment around the closest surface atom, , and determine the sign of the distance as the sign of the scalar product (. This choice is made because if both and are on the same side with respect to the atom representing the location of the local surface, one could consider the i‐th atom to be within the liquid phase.

Results

Computational performances

We have analyzed the computational complexity of the three main algorithms implemented in Pytim, and of the cluster search. As a test system, we have taken a macroscopically planar water/vapor interface. Starting from an initial configuration of 216 water molecules, we generated larger systems by replicating the simulation box in the two directions of the surface plane, obtaining systems up to 55.296 water molecules. In this way, we increased by the same proportion both the system size and the surface area determined by the different algorithms.

In Figures 4, 5, 6, we report the time needed to perform the surface analysis with itim, gitim, and the Willard–Chandler algorithm. In Figure 4, we report also the time needed to assign molecules to the liquid cluster. All algorithms are roughly following the linear scaling for large system sizes, as expected from the theoretical analysis of the average cases. For small system sizes, sublinear scaling is sometimes observed. The exact calculation of the kernel density estimate in the Willard–Chandler algorithm scales also, as expected, worse than its fast, approximated version, although slightly better than quadratic.

Figure 4

Computational complexity of the itim algorithm as implemented in Pytim. The time needed to identify the surface atoms in one configuration is separated into the time to perform the cluster analysis and the actual time for the identification. [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 5

Computational complexity of the gitim algorithm as implemented in Pytim. The calculation of the intrinsic surface atoms has been performed on the same systems as the in Figure 4: gitim is only about 2–3 times slower than itim. [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 6

Computational complexity of the Willard–Chandler algorithm, with the variant of the scipy.stats.gaussian_kde (squares) and with the approximated fast kernel (circles). In both cases, the grid spacing was set to 2 Å. [Color figure can be viewed at http://wileyonlinelibrary.com]

The probe sphere radius was the same (2 Å) in all tests with itim and gitim; for the instantaneous interface calculation using the Willard–Chandler algorithm, we used a gaussian width w = 3 Å and a grid spacing of 2 Å.

Regarding the prefactors, it is worth noticing that, while itim turned out to be the fastest one, gitim performed comparatively well, being only a factor 2–3 slower than itim. The performances of the Willard–Chandler algorithm are affected by the density field discretization (i.e., linear scaling in the number of grid points, or, equivalently, cubic scaling in the inverse grid spacing), but for the reasonable choice of 2 Å for the grid spacing, the performances were only slightly worse than those of gitim.

All tests were performed on a laptop with an Intel(R) Core(TM) i7‐4650 U CPU at 1.70 GHz.

Installation instructions

The latest development code can be retrieved from the address https://github.com/Marcello-Sega/pytim as a git repository, or as a zipped archive. From the source code, the command python setup.py install ‐–user will install the package in the user's directories. Requirements for the installation are the setuptools and numpy python package and Cython. The setup.py script will automatically download and install other required dependencies.

On the other hand, the user can install the latest stable version from the Python package index (https://pypi.python.org) by invoking pip install ‐–user pytim, or from Anaconda (https://anaconda.org/conda-forge) with conda install ‐c conda‐forge pytim.

Usage examples

The simplest script that can be used to compute the interfacial atoms is the following:

import MDAnalysis as mda

import pytim

from pytim.datafiles import WATER_GRO

u = mda.Universe(WATER_GRO)

inter = pytim.ITIM(u)

The script above starts with importing the MDAnalysis and Pytim modules, and one of the example files present in the distribution. Next, the configuration is loaded, and the interfacial atoms are computed using with default values. The interfacial atoms are accessible in a number of ways, one of which being the array inter.atoms that lists all of them, irrespective of the layer/side in which they are located.

In addition to the usual properties that are available for the Atom and AtomGroup classes of MDAnalysis, Pytim introduces new ones, namely, layers, sides, and clusters. These properties can be accessed from each atom in the system, for example, as u.atoms.layers, and give access to which layer the atoms have been assigned to (layers), to which side of a planar interface the atoms are (sides, only) and whether the atoms are in the main cluster or in one of the smaller ones (clusters). These properties can be used to quickly select the atoms of interest. For example, one can access the atoms in the first layer as u.atoms[u.atoms.layers==1].

Some usage examples involving the different methods have been already presented in Figures 1, 2, 3. Here, we describe two typical scripts that can be used to compute intrinsic and non‐intrinsic profiles.

Figure 3

Usage of the WillardChandler class to identify the isodensity continuous surface of a solvated dodecylphosphocholine (DPC) micelle. The image (water not shown) is a section cut of the micelle along with the continuous surface, as visualized using paraview. [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 7 shows an example of how one can compute an intrinsic density profile in Pytim. Note that the Pytim distribution provides several sample configurations and trajectories (in this case, the configuration of the Lennard‐Jones system is pointed at by the label LJ_GRO). Some larger trajectories, which would be too big for the distribution, can be retrieved from the online repository using the pytim_data.fetch() function, as, in this case, for the LJ_BIG_XTC binary trajectory. This script initializes first the interface calculation using , then instantiates a Profile object. The interface option signals that the profile has to be computed relative to the interface itself, that is, it will be an intrinsic profile. Because no observable is passed to the Profile constructor, this defaults to the number density profile calculation. In the successive lines, the code iterates over the trajectory frames, computing the surface atoms automatically every time a new frame is loaded, and sampling the intrinsic profile. The profile is stored internally with 0.01 Å resolution, and can be accessed through the get_values() function with the desired binning, resulting in the plot shown in Figure 7, with the liquid part located on the left (negative position values) of the interface (a delta function centered at zero, not shown) and the vapor part on its right (positive values). The distinction between liquid and vapor is performed by prefiltering the configurations used by using the simple cluster search, which is switched on by using the cluster_cut option. The cluster search is, instead, switched on by using, additionally, the cluster_threshold option.

Figure 7

Calculation of the intrinic density profile of a Lennard‐Jones liquid/vapor interface. [Color figure can be viewed at http://wileyonlinelibrary.com]

Another example, on how to compute the distribution of the surface layers along the interface normal is shown in Figure 8. In this case, one instantiates five different profiles (without specifying an interface option, the code will compute the non‐intrinsic profile), and passes to the sample() function a selection of atoms corresponding to those in the different layers.

Figure 8

Calculation of the distribution of the first four molecular layers in a water/vapor interface. [Color figure can be viewed at http://wileyonlinelibrary.com]

As a final example, in Figure 9 we show the analysis of the intrinsic curvature at the surface of the DPC micelle, together with the calculation of the local normal vector, using the method of Yesylevskyy and Ramseyer.11

Figure 9

Calculation of the local intrinsic curvature for surface atoms (red: positive; blue: negative) and local surface normal for the DPC micelle. [Color figure can be viewed at http://wileyonlinelibrary.com]

Other examples are provided on the github page of the project, in the online manual, and in the example files in the Pytim distribution, including details on how to compute time correlation functions or how to use the pytim tools on top of MDTraj or in an OpenMM simulation run.

List of functionalities

The code, thanks to the use of the MDAnalysis, allows to perform several surface‐related analyses of many popular trajectories and configuration file formats. The full list is available on the website of MDAnalysis (http://www.mdanalysis.org). Both Python 2.7 and Python 3 are supported. The functionalities of the code can be grouped in the following classes:

Interface/surface identification

itim algorithm for macroscopically planar interfaces

gitim algorithm for arbitrary interfaces

Willard–Chandler algorithm for arbitrary interfaces

Lee–Richards SASA algorithm for arbitrary interfaces

Clustering/filtering

simple clustering of nearest neighbors

DBSCAN density based clustering with automatic threshold tuning

Output/Export

In addition to all the file formats available through MDAnalysis, Pytim can save information about the interfaces in some additional ones, namely:

pdb format for particles with information on layering

vtk format for particles

vtk format for volumetric densities

vtk format for surface triangulation

cube format for particles and volumetric density

Wavefront obj format for surface triangulation

Observables

Pytim provides also an abstract class for coding new observables. A set of basic atomic/molecular properties (density, charge, molecular orientation, ...) are already provided in the code, and we plan to expand the set of available observables with time. All observables derived from the abstract class are suitable to be used in the calculation of the following correlations/distributions:

Intrinsic and non‐intrinsic profiles

Three‐ and two‐dimensional radial distribution functions

Time correlation functions, also for variable sets of atoms (both in the continuous and intermittent variants)

Test systems

Several configurations and small trajectory files are provided for testing purposes: the present set of systems includes so far as follows:

Water/vapor interfaces at different temperatures

Water/CCl4 interface

Methanol/vapor interface

Ionic liquid/benzene mixture

A glucose molecule in water

DPC micelle

DPPC bilayer

Fullerene

API documentation and tests

Pytim code has been written using docstrings for all public classes/functions, and for the most important private ones. The code reference manual, generated automatically from the code, is available at https://marcello-sega.github.io/pytim.

Along with the reference manual come several tutorials, which cover several topics ranging from the installation to troubleshooting, and include common usage of the main modules. In the source distribution, there is a directory that includes several python example files that show how to use the main modules and the analysis tools. In addition, jupyter notebooks are also present, and can be visualized (statically) directly from a web browser, without needing to run the code, from https://github.com/Marcello-Sega/pytim.

Code snippets within the documentation of the code are used to provide a testing framework using the doctest module in combination with Sphinx (http://www.sphinx-doc.org/). Any new update of the code, to be incorporated in the main branch, has to pass all tests (more than 200, with a code coverage larger than 85%) using different versions of Python, of the underlying libraries (MDAnalysis) on Linux and Mac OS‐X systems. Automated testing is performed using Travis‐CI (https://travis-ci.org/Marcello-Sega/pytim).

Conclusions

We presented Pytim, a new package for the interfacial analysis of molecular systems. Pytim, is based on the python scripting language and provides the user with a number of analysis tools dedicated to the properties of fluid interfaces. Being built on top of the MDAnalysis library, Pytim allows to analyze trajectory files in the native formats of the most common simulation packages, freeing the user from the burden of conversion. The package implements three main surface analysis algorithms, namely, itim, gitim, and the Willard–Chandler instantaneous liquid interface algorithm, as well as two different clustering algorithms to be used for identifying different phases. All algorithms have been carefully optimized using either fast, vectorized functions from the numpy and scipy libraries, or by writing, when necessary, Cython code, as well as parallel, multiprocess code. All algorithm have a linear or quasi‐linear scaling with respect to the system size. The use of python as the main coding language for this package also allows the user to easily modify and extend the package with minimum effort.