A toolbox to make it easier to deal with materials under periodc boundary conditions.

## Project description

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

## Project details

## Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.