PbcPy

Index

Authors

Fundamentals

Installation

pip install pbcpy

git clone git@gitlab.com:ales.genova/pbcpy.git

DirectCell andReciprocalCell class

>>>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

DirectCell andReciprocalCell properties

# the lattice>>>cell1.latticearray([[10.,0.,0.],       [0.,10.,0.],       [0.,0.,10.]])# the volume>>>cell1.volume1000.0

DirectCell andReciprocalCell methods

\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 class

>>>frompbcpy.baseimportCoord>>>pos1=Coord(pos=[0.5,0.6,0.3],cell=cell1,ctype="Cartesian")

Coord attributes

# the coordinate type (Cartesian or Crystal)>>>pos1.basis'Cartesian'# the cell attribute is a Cell object>>>type(pos1.cell)pbcpy.base.DirectCell

Coord methods

>>>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 classes

>>>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])

Grid attributes

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

DirectField andReciprocalField class

# 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

DirectField attributes

>>>type(rho.grid)pbcpy.grid.DirectGrid>>>rho.span3>>>rho.rank1# the density is a scalar field

DirectField methods

# 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

ReciprocalField methods

# 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 class

System attributes

pbcarray class

# 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

File Formats

PP class

>>>water_dimer=PP(filepp='/path/to/density.pp').read()# the output of PP.read() is a System object.

XSF class

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