Using pip install to get numpy & scipy linked against the Intel MKL

The pip package management system is a very convenient way of managing custom python installations that are not managed by the operating system (such as virtual python environments).

pip install numpy scipy essentially works out of the box. Yet, when it comes to numpy and scipy, significant speed can be gained by linking these modules against the Intel Math Kernel Library (MKL) for linear algebra operations.

You can tell pip to do that in the following way:

(1) Tell numpy where to find the MKL

Create a file .numpy-site.cfg in your home directory with a content similar to 

[mkl]
library_dirs = /opt/intel/composer_xe_2013.1.117/mkl/lib/intel64
include_dirs = /opt/intel/composer_xe_2013.1.117/mkl/include
mkl_libs = mkl_rt

lapack_libs =

(2) Tell pip to use the Intel compiler

pip install numpy \

--global-option="config" \

--global-option="--compiler=intelem" \

--global-option="build_clib" \

--global-option="--compiler=intelem" \

--global-option="build_ext" \

--global-option="--compiler=intelem"

pip install scipy \

--global-option="config" \

--global-option="--compiler=intelem" \

--global-option="--fcompiler=intelem" \

--global-option="build_clib" \

--global-option="--compiler=intelem" \

--global-option="--fcompiler=intelem" \

--global-option="build_ext" \

--global-option="--compiler=intelem" \

--global-option="--fcompiler=intelem"

Alternatively, you may also create a file .pydistutils.cfg in your home directory with content

[config]

compiler=intelem
fcompiler=intelem

[build_clib]

compiler=intelem

fcompiler=intelem
[build]

compiler=intelem
[build_ext]
compiler=intelem

fcompiler=intelem

and then simply run

pip install numpy
pip install scipy

Note: If you don't provide this info to pip, it will look for and use the GNU compiler collection by default.

Note: While numpy installed without problems I did receive a compilation error with icc 14.0.0. Possibly related to this 

Parts (1) and (2) relied on the following StackExchange posts

• http://stackoverflow.com/questions/13769936/supplying-numpy-site-cfg-arguments-to-pip
• http://stackoverflow.com/questions/15725869/pass-build-ext-options-to-pip-install

Importing cube files into Blender

In this post, we are going to import the grid data from a Gaussian cube file into Blender, an open-source 3d graphics software.

Why this effort? In the atomistic simulations community, we are used to displaying wave functions and electron densities as isosurfaces and lots of software exists that makes this process very easy. But an isosurface doesn't tell the whole story. Blender gives us the opportunity to visualize the electron density as a kind of "smoke cloud", where the density determines how light scatters going through the cloud (we may then add an intrinsic glow to the cloud as well). While this is not needed in the daily routine of understanding the basic features of an electron distribution, it allows to convey the inherent fuzzyness in a single image.

Of course, there is many more reasons, why one may want to use Blender. Among its many features, some of the ones relevant to atomistic simulations are

1. a great graphical user interface [1]
2. ray-tracing capabilities
3. advanced animation capabilities (including physics!)
4. a Python interface that allows for automating workflows

However, Blender was originally developed for animation studios, and so the import of file formats customary in atomistic simulations is not always straightforward.

For atomic structure data, several routes are already available

• The XYZ addon for directly importing XYZ (.xyz) files into Blender
• Via VMD, exporting to the Wavefront Objects (.obj/.mtl), which can be imported by Blender
  This allows not only to import atomic structure models, but also isosurfaces etc.
• Via ASE, exporting to the X3D (.x3d) format, which can be imported by Blender

But how would we go about importing a function f(x,y,z) defined on a regular 3d grid, such as an orbital wave function or electron density? In this post, we are going to import the grid data from a Gaussian cube file.

First, we have to note that the capability to render this kind data (called "voxels" in the Blender community) is currently available only for the internal "Blender Render" engine. While the "Cycles Render" engine is being used and developed more actively, we will have to abandon it for now. "Blender Internal" offers two alternative file formats for importing voxel data

• "8bit raw", in which f(x,y,z) is stored as integers ranging from 0 to 255
• "Blender voxel", in which f(x,y,z) is stored as 32bit floats. 

The first format is basically unusable for wave functions/densities of molecules - a relative accuracy of 1/256 doesn't get you very far when sampling exponential decays. We thus choose the latter format.

So, we need to convert the .cube file into a .bvox file. Adapting a post by Matthias Meschede, this is easily achieved using a combination of ASE and numpy: 

#!/usr/bin/env python

import ase.io.cube

import numpy as np

import os

fname = 'PROJ-WFN_00183_1-1_0.cube'

print("Reading cube file {}".format(fname))

data, atoms = ase.io.cube.read_cube_data(fname)

# Here, I want the electron density, not the wave function

data = data**2

# If data is too large, just reduce it by striding with steps >1

sx, sy, sz = 1, 1, 1

data = data[::sx,::sy,::sz]

# Note the reversed order!!

nz, ny, nx = data.shape

nframes = 1

header = np.array([nx,ny,nz,nframes])

#open and write to file

vdata = data.flatten() / np.max(data)

vfname = os.path.splitext(fname)[0] + '.bvox'

vfile = open(vfname,'wb')

print("Writing Blender voxel file {}".format(vfname))

header.astype('<i4').tofile(vfile)

vdata.astype('<f4').tofile(vfile)

Now we need to load the file into Blender. The steps are essentially the following

1. Open Blender and select "Blender Render"
2. Select the default cube, set its material to "Volume" and its "Density" to 0
3. Add a new texture of type "Voxel Data", file format "Blender Voxel" and load the file
   Note: You don't get any confirmation of the file being loaded but if you don't get an error message, assume that it is.
4. Activate "Still frame" and set still frame number to 1 (not 0!)
   
5. Activate influence on density and emission and ramp up the intensity until you see something in the rendered output. Note: You may never see anything in the preview window - don't let this confuse you. should be able to continue by yourself and unleash your creativity on the voxels!
   

To give you one possible idea of what can be done, here is what I came up with for my first Blender project

End state of N=7 armchair graphene nanoribbon

Cheers!

[1] This may sound trivial, but alternatives, such as POV-ray may have no gui whatsoever.