/pbcpy

Python package providing some useful tools when dealing with molecules and materials under periodic boundary conditions and uniform grids. This is a mirror of https://gitlab.com/ales.genova/pbcpy

Primary LanguagePythonMIT LicenseMIT

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

PyPI version PyPI status pipeline status coverage report License: MIT

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

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)