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

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

## Project details

## Download files

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