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 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)

alesgenova/pbcpy

PbcPy

Index

Authors

Fundamentals

Installation

DirectCell and ReciprocalCell class

DirectCell and ReciprocalCell properties

DirectCell and ReciprocalCell methods

Coord class

Coord attributes

Coord methods

DirectGrid and ReciprocalGrid classes

Grid attributes

DirectField and ReciprocalField class

DirectField attributes

DirectField methods

ReciprocalField methods

System class

System attributes

pbcarray class