- Notifications
You must be signed in to change notification settings - Fork1
Python package providing some useful tools when dealing with molecules and materials under periodic boundary conditions and uniform grids. This is a mirror ofhttps://gitlab.com/ales.genova/pbcpy
License
alesgenova/pbcpy
Folders and files
| Name | Name | Last commit message | Last commit date | |
|---|---|---|---|---|
Repository files navigation
PbcPy is no longer under active development. Please use its forkDFTpy instead, which is mantained by the Pavanello Research Group and has many additional features.
pbcpy is a Python3 package providing some useful abstractions to deal withmolecules and materials under periodic boundary conditions (PBC).
In addition,pbcpy exposes a fully periodic N-rank array, thepbcarray, which is derived from thenumpy.ndarray.
Finally,pbcpy provides IO support to some common file formats:
- Quantum Espresso
.ppformat (read only) - XCrySDen
.xsfformat (write only)
- Authors
- Fundamentals
- Installation
- DirectCell andReciprocalCell class
- Coord class
- DirectGrid andReciprocalGrid class
- DirectField andReciprocalField class
- System class
- pbcarray class
- File Formats
pbcpy has been developed @Pavanello Research Group by:
- Alessandro Genova
with contributions from:
- Tommaso Pavanello
- Michele Pavanello
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;
Installpbcpy throughPyPI
pip install pbcpyInstall the dev version fromgitlab
git clone git@gitlab.com:ales.genova/pbcpy.gitNOTE:pbcpy is in the early stages of development, classes and APIs are bound to be changed without prior notice.
A unit cell is defined by its lattice vectors. To create aDirectCell object we need to provide it a3x3 matrix containing the lattice vectors (as columns).pbcpy expects atomic units, a flexible units system might be addedd in the future.
>>>frompbcpy.baseimportDirectCell,ReciprocalCell>>>importnumpyasnp>>>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
lattice: the lattice vectors (as columns)volume: the volume of the cellorigin: the origin of the Cartesian reference frame
# the lattice>>>cell1.latticearray([[10.,0.,0.], [0.,10.,0.], [0.,0.,10.]])# the volume>>>cell1.volume1000.0
==operator : compare twoCellobjectsget_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==cell2True
Coord is anumpy.array derived class, with some additional attributes and methods.Coordinates in a periodic system are meaningless without the reference unit cell, this is why aCoord object also has an embeddedDirectCell attribute.Also, coordinates can be either expressed in either a"Cartesian" or"Crystal" basis.
>>>frompbcpy.baseimportCoord>>>pos1=Coord(pos=[0.5,0.6,0.3],cell=cell1,ctype="Cartesian")
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.DirectCell
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-pos1Coord([0.1,-1. ,2. ])>>> (pos2-pos1).length()22.383029285599392
DirectGrid andReciprocalGrid are subclasses ofDirectGrid andReciprocalGrid respectively.Grids inherit all the attributes and methods of their respectiveCells, and have a few of their own to deal with quantities represented on a equally spaced grid.
>>>frompbcpy.gridimportDirectGrid# A 10x10x10 Bohr Grid, with 100x100x100 gridpoints>>>lattice=np.identity(3)*10>>>grid1=DirectGrid(lattice=lattice,nr=[100,100,100],origin=[0,0,0])
- 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.dV0.001# Grid points for each direction>>>grid1.nrarray([100,100,100])# Total number of grid points>>>grid1.nnr1000000# Cartesian coordinates at each grid point>>>grid1.rCoord([[[[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.sCoord([[[[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 methosreciprocal_grid1=grid1.get_reciprocal()# reciprocal_grid1 is an instance of ReciprocalGrid>>>reciprocal_grid1.garray([[[[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.ggarray([[[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)
TheDirectField andReciprocalField classes represent a scalar field on aDirectGrid andReciprocalGrid respectively. These classes are extensions of thenumpy.ndarray.
Operations such as interpolations, fft and invfft, and taking arbitrary 1D/2D/3D cuts are made very easy.
ADirectField can be generated directly from Quantum Espresso postprocessing.pp files (see below).
# A DirectField example>>>frompbcpy.fieldimportDirectField>>>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>>>frompbcpy.formats.qeppimportPP>>>water_dimer=PP(filepp="/path/to/density.pp").read()>>>rho=water_dimer.field>>>type(rho)pbcpy.field.DirectField
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.span3>>>rho.rank1# the density is a scalar field
- 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.span2>>>rho_cut.grid.latticearray([[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
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.gridTrue>>>np.isclose(rho1,rho).all()True# as expected ifft(fft(rho)) = rho
System is simply a class containing aDirectCell (orDirectGrid), a set of atomsions, and aDirectField
name: arbitrary nameions: collection of atoms and their coordinatescell: the unit cell of the system (DirectCellorDirectGrid)field: an optionalDirectFieldobject.
pbcarray is a sublass ofnumpy.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.>>>frompbcpy.baseimportpbcarray>>>importmatplotlib.pyplotasplt>>>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
pbcpy can read a Quantum Espresso post-processing.pp file into aSystem object.
>>>water_dimer=PP(filepp='/path/to/density.pp').read()# the output of PP.read() is a System object.
pbcpy can write aSystem 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)
About
Python package providing some useful tools when dealing with molecules and materials under periodic boundary conditions and uniform grids. This is a mirror ofhttps://gitlab.com/ales.genova/pbcpy