What is Pytim | Examples | More info | How to Install | References

_{If you try this software out and have some suggestions, remarks, or bugfixes, feel free to comment here on GitHub and/or make a pull request.}

Jupyter Notebooks with more examples are available at Marcello-Sega/pytim-notebooks

The paper about pytim has been published on J. Comput. Chem. It is open access and you can download the pdf from Wiley <img src="https://licensebuttons.net/l/by-nc/4.0/80x15.png"> (see also the references)

What is Pytim

Pytim is a cross-platform python implementation of several methods for the detection of fluid interfaces in molecular simulations. It is based on MDAnalysis, but it integrates also seamlessly with MDTraj, and can be even used for online analysis during an OpenMM simulation (see further down for examples with MDTraj and with OpenMM).

So far the following interface/phase identification methods have been implemented: <img src="https://github.com/Marcello-Sega/pytim/raw/IMAGES/_images/micelle_cut.png" width="380" align="right" style="z-index:999;">

• ITIM
• GITIM
• SASA
• Willard Chandler
• DBSCAN filtering

Supported formats

Pytim relies on the MDAnalysis package for reading/writing trajectories, and work therefore seamlessly for a number of popular trajectory formats, including:

• GROMACS
• CHARMM/NAMD
• LAMMPS
• AMBER
• DL_Poly

as well as common structure file formats such as XYZ or PDB (have a look at the complete list)

Install from PyPi or Anaconda -

PyPi: pip install --user --upgrade pytim

Anaconda: conda install -c conda-forge pytim

NOTE: on Mac OS you might want to use CFLAGS='-stdlib=libc++' pip install --user --upgrade pytim

<a name="example"></a> Show me an example usage, now!

Ok, ok ... have a look below: this example is about computing molecular layers of a flat interface:

Step 1: interface identification

import MDAnalysis as mda
import pytim
from pytim.datafiles import WATER_GRO

# load the configuration in MDAnalysis
u = mda.Universe(WATER_GRO)

# compute the interface using ITIM. Identify 4 layers.
inter = pytim.ITIM(u,max_layers=4)

That's it. There's no step 2!

Now interfacial atoms are accessible in different ways, pick the one you like:

1. Through the atoms group accessible as

inter.atoms.positions # this is a numpy array holding the position of atoms in the layers

array([[ 22.10000038,  16.05999947,  94.19633484],
       [ 22.43999863,  16.97999954,  93.96632385],
       ...,
       [ 33.70999908,  49.02999878,  62.52632904],
       [ 34.06999969,  48.18000031,  61.16632843]], dtype=float32)

2. Using the label that each atom in the MDAnalysis universe now has, which specifies in which layer it is found:

u.atoms.layers  # -1 if not in any layer, 1 if in the first layer, ...

3. Using the layers groups, stored as a list (of lists, in case of upper/lower layers in flat interfaces) of groups:

inter.layers

array([[<AtomGroup with 780 atoms>, <AtomGroup with 690 atoms>,
        <AtomGroup with 693 atoms>, <AtomGroup with 660 atoms>],
       [<AtomGroup with 777 atoms>, <AtomGroup with 687 atoms>,
        <AtomGroup with 663 atoms>, <AtomGroup with 654 atoms>]], dtype=object)

Visualisation

Pytim can export in different file formats: the PDB format with the beta field used to tag the layers, VTK, cube for both continuous surfaces and particles, and, of course, all formats supported by MDAnalysis.

In VMD, for example, using beta == 1 allows you to select all atoms in the first interfacial layer. Just save your PDB file with layers information using

inter.writepdb('myfile.pdb')

In a jupyter notebook, you can use nglview like this:

import nglview
v = nglview.show_mdanalysis(u)
v.camera='orthographic'
v.center()
system = v.component_0
colors = ['','red','orange','yellow','white']

for n in [1,2,3,4]:
    system.add_spacefill(selection = inter.atoms[inter.atoms.layers==n].indices, color=colors[n] )

system.add_spacefill(selection = (u.atoms - inter.atoms).indices, color='gray' )

# this must go in a separate cell
v.display()

<p align="center"> <img src="https://github.com/Marcello-Sega/pytim/raw/IMAGES/_images/output_13_0.png" width="60%" align="center" style="z-index:999;"> </p>

Analysing trajectories (MDAnalysis and mdtraj)

Once one of the pytim classes (ITIM,GITIM,WillardChandler,...) has been initialised, whenever a new frame is loaded, the interfacial properties will be calculated automatically without need of doing anything else

MDAnalysis example

import MDAnalysis as mda
import pytim 
from pytim.datafiles import WATER_GRO, WATER_XTC 

u = mda.Universe(WATER_GRO,WATER_XTC)
inter = pytim.ITIM(u)
for step in u.trajectory[:]:
    print "surface atoms:", repr(inter.atoms)

mdtraj example

Under the hood, pytim uses MDAnalysis, but this is made (almost completely) transparent to the user, so that interoperability with other software is easy to implement. For example, to analyse a trajectory loaded with MDTraj, it is enough to do the following:

import mdtraj
import pytim                     
from pytim.datafiles import WATER_GRO, WATER_XTC 

t = mdtraj.load_xtc(WATER_XTC,top=WATER_GRO) 
inter = pytim.ITIM(t) 
for step in t[:]:
        print "surface atoms:" , repr(inter.atoms.indices)

openmm example

Another example is using pytim to perform online interfacial analysis during an OpenMM simulation:

# openmm imports
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
# pytim
import pytim
from pytim.datafiles import WATER_PDB

# usual openmm setup, we load one of pytim's example files
pdb = PDBFile(WATER_PDB)
forcefield = ForceField('amber99sb.xml', 'spce.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
        nonbondedCutoff=1*nanometer)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)

# just pass the openmm Simulation object to pytim
inter = pytim.ITIM(simulation)
print repr(inter.atoms)

# the new interfacial atoms will be computed at the end
# of the integration cycle
simulation.step(10)
print repr(inter.atoms)

<a name="non-flat-interfaces"></a> What if the interface is not flat?

You could use GITIM, or the Willard-Chandler method. Let's have a look first at GITIM:

import MDAnalysis as mda
import pytim
from   pytim.datafiles import *

u       = mda.Universe(MICELLE_PDB)
g       = u.select_atoms('resname DPC')

interface = pytim.GITIM(u,group=g,molecular=False,
                        symmetry='spherical',alpha=2.5)
layer = interface.layers[0]
interface.writepdb('gitim.pdb',centered=False)
print repr(layer)
<AtomGroup with 872 atoms>

nglview can be used to show a section cut of the micelle in a jupyter notebook:

import nglview
import numpy as np

v = nglview.show_mdanalysis(u)
v.camera='orthographic'
v.center()
v.add_unitcell()
system = v.component_0
system.clear()

selection = g.atoms.positions[:,2]>30.
system.add_spacefill(selection =g.atoms.indices[selection])

selection = np.logical_and(inter.atoms.layers==1,inter.atoms.positions[:,2]>30.)
system.add_spacefill(selection =inter.atoms.indices[selection], color='yellow' )

v.display()

<p align="center"> <img src="https://github.com/Marcello-Sega/pytim/raw/IMAGES/_images/micelle-gitim.png" width="40%" align="middle"> </p>

The Willard-Chandler method can be used, instead to find out isodensity surfaces:

import MDAnalysis as mda
import pytim
from pytim.datafiles import MICELLE_PDB
import nglview

u = mda.Universe(MICELLE_PDB)
g = u.select_atoms('resname DPC')

interface = pytim.WillardChandler(u, group=g, mesh=1.5, alpha=3.0)
interface.writecube('data.cube')

Done, the density file is written in .cube format, we can now open it with programs such as [Paraview](https://ww