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
.ppformat (read only) - XCrySDen
.xsfformat (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
DirectCellandCoordclasses which define a unit cell under PBC in real space, and a cartesian/crystal coordinate respectively;ReciprocalCellclass which defines a cell in reciprocal space;DirectGridandReciprocalGridclasses, which are derived fromDirectCellandReciprocalCelland provide space discretization;DirectFieldandReciprocalField, classes to represent a scalar (such as an electron density or a potential) and/or vector fields associated to either aDirectGridor 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 cellDirectCell 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.0DirectCell and ReciprocalCell methods
-
==operator : compare twoCellobjects -
get_reciprocal: returns a newReciprocalCellobject that is the "reciprocal" cell of self (if self is aDirectCell) -
get_direct: returns a newDirectCellobject 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
TrueCoord 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: theDirectCellobject 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.DirectCellCoord methods
to_crys(),to_cart(): convertselfto crystal or cartesian basis (returns a newCoordobject).d_mic(other): Calculate the vector connecting two coordinates (from self to other), using the minimum image convention (MIC). The result is itself aCoordobject.dd_mic(other): Calculate the distance between two coordinates, using the MIC.+/-operators : Calculate the difference/sum between two coordinates without using the MIC.basisconversions are automatically performed when needed. The result is itself aCoordobject.
>>> 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.383029285599392DirectGrid 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 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(DirectGridonly)s: crystal coordinates at each grid point. A rank 3 array of typeCoord(DirectGridonly)g: G vector at each grid point (ReciprocalGridonly)gg: Square of G vector at each grid point (ReciprocalGridonly)
# 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.DirectFieldDirectField attributes
grid: Represent the grid associated to the field (it's aDirectGridorReciprocalGridobject)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 fieldDirectField 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 newDirectFieldobject). 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.ReciprocalFieldReciprocalField 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)) = rhoSystem 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 (DirectCellorDirectGrid)field: an optionalDirectFieldobject.
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 resultFile 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)