PbcPy is no longer under active development. Please use its fork DFTpy instead, which is mantained by the Pavanello Research Group and has many additional features.
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 Nrank 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
andCoord
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
andReciprocalGrid
classes, which are derived fromDirectCell
andReciprocalCell
and provide space discretization;DirectField
andReciprocalField
, classes to represent a scalar (such as an electron density or a potential) and/or vector fields associated to either aDirectGrid
or aReciprocalGrid
;
Installation
Install pbcpy
through PyPI
pip install pbcpy
Install the dev version from gitlab
git clone [email protected]: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 cellorigin
: 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 twoCell
objects 
get_reciprocal
: returns a newReciprocalCell
object that is the "reciprocal" cell of self (if self is aDirectCell
) 
get_direct
: returns a newDirectCell
object that is the "direct" cell of self (if self is aReciprocalCell
)
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
: theDirectCell
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()
: convertself
to crystal or cartesian basis (returns a newCoord
object).d_mic(other)
: Calculate the vector connecting two coordinates (from self to other), using the minimum image convention (MIC). The result is itself aCoord
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 aCoord
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. Grid
s inherit all the attributes and methods of their respective Cell
s, 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 quantitiesnr
: array, number of grid point for each directionnnr
: total number of points in the gridr
: cartesian coordinates at each grid point. A rank 3 array of typeCoord
(DirectGrid
only)s
: crystal coordinates at each grid point. A rank 3 array of typeCoord
(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 aDirectGrid
orReciprocalGrid
object)span
: The number of dimensions of the grid for which the number of points is larger than 1rank
: The number of dimensions of the quantity at each grid point1
: scalar field (e.g. the rank of rho is1
)>1
: vector field (e.g. the rank of the gradient of rho is3
)
>>> 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 newDirectField
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 ofReciprocalField
, which contains the appropriateReciprocalGrid
# 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 ofDirectField
, which contains the appropriateDirectGrid
# 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 nameions
: collection of atoms and their coordinatescell
: the unit cell of the system (DirectCell
orDirectGrid
)field
: an optionalDirectField
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 postprocessing .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)