PbcPy

pbcpy is a Python3 package providing some useful abstractions to deal with molecules and materials under periodic boundary conditions (PBC).

In addition, pbcpy exposes a fully periodic N-rank array, the pbcarray, which is derived from the numpy.ndarray.

Finally, pbcpy provides IO support to some common file formats:

• Quantum Espresso .pp format (read only)
• XCrySDen .xsf format (write only)

Index

• Authors
• Fundamentals
• Installation
• DirectCell and ReciprocalCell class
• Coord class
• DirectGrid and ReciprocalGrid class
• DirectField and ReciprocalField class
• System class
• pbcarray class
• File Formats

Authors

pbcpy has been developed @ Pavanello Research Group by:

• Alessandro Genova

with contributions from:

• Tommaso Pavanello
• Michele Pavanello

Fundamentals

• DirectCell and Coord classes which define a unit cell under PBC in real space, and a cartesian/crystal coordinate respectively;
• ReciprocalCell class which defines a cell in reciprocal space;
• DirectGrid and ReciprocalGrid classes, which are derived from DirectCell and ReciprocalCell and provide space discretization;
• DirectField and ReciprocalField, classes to represent a scalar (such as an electron density or a potential) and/or vector fields associated to either a DirectGrid or a ReciprocalGrid;

Installation

Install pbcpy through PyPI

pip install pbcpy

Install the dev version from gitlab

git clone git@gitlab.com:ales.genova/pbcpy.git

NOTE: pbcpy is in the early stages of development, classes and APIs are bound to be changed without prior notice.

DirectCell and ReciprocalCell class

A unit cell is defined by its lattice vectors. To create a DirectCell object we need to provide it a 3x3 matrix containing the lattice vectors (as columns). pbcpy expects atomic units, a flexible units system might be addedd in the future.

>>> from pbcpy.base import DirectCell, ReciprocalCell
>>> import numpy as np
>>> lattice = np.identity(3)*10 # Make sure that at1 is of type numpy array.
>>> cell1 = DirectCell(lattice=lattice, origin=[0,0,0]) # 10 Bohr cubic cell

DirectCell and ReciprocalCell properties

• lattice : the lattice vectors (as columns)
• volume : the volume of the cell
• origin : the origin of the Cartesian reference frame

# the lattice
>>> cell1.lattice
array([[ 10.,   0.,   0.],
       [  0.,  10.,   0.],
       [  0.,   0.,  10.]])

# the volume
>>> cell1.volume
1000.0

DirectCell and ReciprocalCell methods

• == operator : compare two Cell objects

• get_reciprocal: returns a new ReciprocalCell object that is the "reciprocal" cell of self (if self is a DirectCell)

• get_direct: returns a new DirectCell object that is the "direct" cell of self (if self is a ReciprocalCell)

Note, by default the physics convention is used when converting between direct and reciprocal lattice:

\big[\text{reciprocal.lattice}\big]^T = 2\pi \cdot \big[\text{direct.lattice}\big]^{-1}

>>> reciprocal_cell1 = cell1.get_reciprocal()
>>> print(reciprocal_cell1.lattice)
array([[ 0.62831853,  0. ,  0. ],
       [ 0. ,  0.62831853,  0. ],
       [ 0. ,  0. ,  0.62831853]])

>>> cell2 = reciprocal_cell1.get_direct()
>>> print(cell2.lattice)
array([[ 10.,  0.,  0.],
       [ 0.,  10.,  0.],
       [ 0.,  0.,  10.]])

>>> cell1 == cell2
True

Coord class

Coord is a numpy.array derived class, with some additional attributes and methods. Coordinates in a periodic system are meaningless without the reference unit cell, this is why a Coord object also has an embedded DirectCell attribute. Also, coordinates can be either expressed in either a "Cartesian" or "Crystal" basis.

>>> from pbcpy.base import Coord
>>> pos1 = Coord(pos=[0.5,0.6,0.3], cell=cell1, ctype="Cartesian")

Coord attributes

• basis : the coordinate type: 'Cartesian' or 'Crystal'.
• cell : the DirectCell object associated to the coordinates.

# the coordinate type (Cartesian or Crystal)
>>> pos1.basis
'Cartesian'

# the cell attribute is a Cell object
>>> type(pos1.cell)
pbcpy.base.DirectCell

Coord methods

• to_crys(), to_cart() : convert self to crystal or cartesian basis (returns a new Coord object).
• d_mic(other) : Calculate the vector connecting two coordinates (from self to other), using the minimum image convention (MIC). The result is itself a Coord object.
• dd_mic(other) : Calculate the distance between two coordinates, using the MIC.
• +/- operators : Calculate the difference/sum between two coordinates without using the MIC. basis conversions are automatically performed when needed. The result is itself a Coord object.

>>> pos1 = Coord(pos=[0.5,0.0,1.0], cell=cell1, ctype="Crystal")
>>> pos2 = Coord(pos=[0.6,-1.0,3.0], cell=cell1, ctype="Crystal")

# convert to Crystal or Cartesian (returns new object)
>>> pos1.to_cart() 
Coord([  5.,   0.,  10.]) # the coordinate was already Cartesian, the result is still correct.
>>> pos1.to_crys()
Coord([ 0.5,  0. ,  1. ]) # the coordinate was already Crystal, the result is still correct.

## vector connecting two coordinates (using the minimum image convention), and distance
>>> pos1.d_mic(pos2)
Coord([ 0.1,  0. ,  0. ])
>>> pos1.dd_mic(pos2)
0.99999999999999978

## vector connecting two coordinates (without using the minimum image convention) and distance
>>> pos2 - pos1
Coord([ 0.1, -1. ,  2. ])
>>> (pos2 - pos1).length()
22.383029285599392

DirectGrid and ReciprocalGrid classes

DirectGrid and ReciprocalGrid are subclasses of DirectGrid and ReciprocalGrid respectively. Grids inherit all the attributes and methods of their respective Cells, and have a few of their own to deal with quantities represented on a equally spaced grid.

>>> from pbcpy.grid import DirectGrid
# A 10x10x10 Bohr Grid, with 100x100x100 gridpoints
>>> lattice = np.identity(3)*10
>>> grid1 = DirectGrid(lattice=lattice, nr=[100,100,100], origin=[0,0,0])

Grid attributes

• All the attributes inherited from Cell
• dV : the volume of a single point, useful when calculating integral quantities
• nr : array, number of grid point for each direction
• nnr : total number of points in the grid
• r : cartesian coordinates at each grid point. A rank 3 array of type Coord (DirectGrid only)
• s : crystal coordinates at each grid point. A rank 3 array of type Coord (DirectGrid only)
• g : G vector at each grid point (ReciprocalGrid only)
• gg : Square of G vector at each grid point (ReciprocalGrid only)

# The volume of each point
>>> grid1.dV
0.001

# Grid points for each direction
>>> grid1.nr
array([100, 100, 100])

# Total number of grid points
>>> grid1.nnr
1000000

# Cartesian coordinates at each grid point
>>> grid1.r
Coord([[[[ 0. ,  0. ,  0. ],
       	 [ 0. ,  0. ,  0.1],
         [ 0. ,  0. ,  0.2],
         [ 0. ,  0. ,  0.3],
                        ...]]])

>>> grid1.r.shape
(100, 100, 100, 3)

>>> grid1.r[0,49,99]
Coord([ 0. ,  4.9,  9.9])

# Crystal coordinates at each grid point
>>> grid1.s
Coord([[[[ 0.  ,  0.  ,  0.  ],
  	 [ 0.  ,  0.  ,  0.01],
       	 [ 0.  ,  0.  ,  0.02],
         [ 0.  ,  0.  ,  0.03],
			  ...]]]])

# Since DirectGrid inherits from DirectCell, we can still use the get_reciprocal methos
reciprocal_grid1 = grid1.get_reciprocal()

# reciprocal_grid1 is an instance of ReciprocalGrid
>>> reciprocal_grid1.g
array([[[[ 0.  ,  0.  ,  0.  ],
         [ 0.  ,  0.  ,  0.01],
         [ 0.  ,  0.  ,  0.02],
         ..., 
         [ 0.  ,  0.  , -0.03],
         [ 0.  ,  0.  , -0.02],
         [ 0.  ,  0.  , -0.01]],
         		   ...]]])

>>> reciprocal_grid1.g.shape
(100, 100, 100, 3)

>>> reciprocal_grid1.gg
array([[[ 0.    ,  0.0001,  0.0004, ...,  0.0009,  0.0004,  0.0001],
        [ 0.0001,  0.0002,  0.0005, ...,  0.001 ,  0.0005,  0.0002],
        [ 0.0004,  0.0005,  0.0008, ...,  0.0013,  0.0008,  0.0005],
        ..., 
        [ 0.0009,  0.001 ,  0.0013, ...,  0.0018,  0.0013,  0.001 ],
        [ 0.0004,  0.0005,  0.0008, ...,  0.0013,  0.0008,  0.0005],
        [ 0.0001,  0.0002,  0.0005, ...,  0.001 ,  0.0005,  0.0002]],
        ...,
                                                                  ]])

>>> reciprocal_grid1.gg.shape
(100, 100, 100)                                          

DirectField and ReciprocalField class

The DirectField and ReciprocalField classes represent a scalar field on a DirectGrid and ReciprocalGrid respectively. These classes are extensions of the numpy.ndarray.

Operations such as interpolations, fft and invfft, and taking arbitrary 1D/2D/3D cuts are made very easy.

A DirectField can be generated directly from Quantum Espresso postprocessing .pp files (see below).

# A DirectField example
>>> from pbcpy.field import DirectField
>>> griddata = np.random.random(size=grid1.nr)
>>> field1 = DirectField(grid=grid1, griddata_3d=griddata)

# When importing a Quantum Espresso .pp files a DirectField object is created
>>> from pbcpy.formats.qepp import PP
>>> water_dimer = PP(filepp="/path/to/density.pp").read()
>>> rho = water_dimer.field
>>> type(rho)
pbcpy.field.DirectField

DirectField attributes

• grid : Represent the grid associated to the field (it's a DirectGrid or ReciprocalGrid object)
• span : The number of dimensions of the grid for which the number of points is larger than 1
• rank : The number of dimensions of the quantity at each grid point
  • 1 : scalar field (e.g. the rank of rho is 1)
  • >1 : vector field (e.g. the rank of the gradient of rho is 3)

>>> type(rho.grid)
pbcpy.grid.DirectGrid

>>> rho.span
3

>>> rho.rank
1
# the density is a scalar field

DirectField methods

• Any method inherited from numpy.array.
• integral : returns the integral of the field.
• get_3dinterpolation : Interpolates the data to a different grid (returns a new DirectField object). 3rd order spline interpolation.
• get_cut(r0, [r1], [r2], [origin], [center], [nr]) : Get 1D/2D/3D cuts of the scalar field, by providing arbitraty vectors and an origin/center.
• fft : Calculates the Fourier transform of self, and returns an instance of ReciprocalField, which contains the appropriate ReciprocalGrid

# Integrate the field over the whole grid
>>> rho.integral()
16.000000002898673 # the electron density of a water dimer has 16 valence electrons as expected

# Interpolate the scalar field from one grid to another
>>> rho.shape
(125, 125, 125)

>>> rho_interp = rho.get_3dinterpolation([90,90,90])
>>> rho_interp.shape
(90, 90, 90)

>> rho_interp.integral()
15.999915251442873


# Get arbitrary cuts of the scalar field.
# In this example get the cut of the electron density in the plane of the water molecule
>>> ppfile = "/path/to/density.pp"
>>> water_dimer = PP(ppfile).read()

>>> o_pos = water_dimer.ions[0].pos
>>> h1_pos = water_dimer.ions[1].pos
>>> h2_pos = water_dimer.ions[2].pos

>>> rho_cut = rho.get_cut(r0=o_h1_vec*4, r1=o_h2_vec*4, center=o_pos, nr=[100,100])

# plot_cut is itself a DirectField instance, and it can be either exported to an xsf file (see next session)
# or its values can be analized/manipulated in place.
>>> rho_cut.shape
(100,100)
>>> rho_cut.span
2
>>> rho_cut.grid.lattice
array([[ 1.57225214, -6.68207161, -0.43149218],
       [-1.75366585, -3.04623853,  0.8479004 ],
       [-7.02978121,  0.97509868, -0.30802502]])

# plot_cut is itself a Grid_Function_Base instance, and it can be either exported to an xsf file (see next session)
# or its values can be analized/manipulated in place.
>>> plot_cut.values.shape
(200, 200)

# Fourier transform of the DirectField
>>> rho_g = rho.fft()
>>> type(rho_g)
pbcpy.field.ReciprocalField

ReciprocalField methods

• ifft : Calculates the inverse Fourier transform of self, and returns an instance of DirectField, which contains the appropriate DirectGrid

# inv fft:
# recall that rho_g = fft(rho)
>>> rho1 = rho_g.ifft()
>>> type(rho1)
pbcpy.field.DirectField

>>> rho1.grid == rho.grid
True

>>> np.isclose(rho1, rho).all()
True
# as expected ifft(fft(rho)) = rho

System class

System is simply a class containing a DirectCell (or DirectGrid), a set of atoms ions, and a DirectField

System attributes

• name : arbitrary name
• ions : collection of atoms and their coordinates
• cell : the unit cell of the system (DirectCell or DirectGrid)
• field : an optional DirectField object.

pbcarray class

pbcarray is a sublass of numpy.ndarray, and is suitable to represent periodic quantities, by including robust wrapping capabilities. pbcarray can be of any rank, and it can be freely sliced.

# 1D example, but it is valid for any rank.
>>> from pbcpy.base import pbcarray
>>> import  matplotlib.pyplot as plt
>>> x = np.linspace(0,2*np.pi, endpoint=False, num=100)
>>> y = np.sin(x)
>>> y_pbc = pbcarray(y)
>>> y_pbc.shape
(100,) 							# y_pbc only has 100 elements, but we can freely do operations such as:
>>> plt.plot(y_pbc[-100:200])	# and get the expected result

File Formats

PP class

pbcpy can read a Quantum Espresso post-processing .pp file into a System object.

>>> water_dimer = PP(filepp='/path/to/density.pp').read() 
# the output of PP.read() is a System object.

XSF class

pbcpy can write a System object into a XCrySDen .xsf file.

>>> XSF(filexsf='/path/to/output.xsf').write(system=water_dimer)

# an optional field parameter can be passed to XSF.write() in order to override the DirectField in system.
# This is especially useful if one wants to output one system and an arbitrary cut of the grid,
# such as the one we generated earlier
>>> XSF(filexsf='/path/to/output.xsf').write(system=water_dimer, field=rho_cut)