I'm a chemistry student working on a research team and part of my job is to produce simple computational scripts that can be understood and used by members unfamiliar with Python or other languages. I was recently tasked with producing a script for creating crystal lattices of any desired 3-dimensional permutation of AxBxC in reduced coordinates (aka fractional coordinates) using a dictionary of the primitive cell, which is essentially a minimum volume unit cell. The user sets their desired parameters, which are found in the beginning of the original script and in the userParam function of the second script, and runs the script.
Example of output:
*******************************************************
1 x 1 x 2 PbCO3 14 supercell
x-shift: 0.000
y-shift: 0.000
z-shift: 0.000
Sorting priority: z-y-x
*******************************************************
x y z
1: O | 0.010000000 | -0.800000000 | -0.195000000
2: Pb | -0.251000000 | -0.574000000 | -0.113500000
3: O | -0.440000000 | -0.310000000 | -0.055000000
4: O | -0.230000000 | -0.060000000 | -0.055000000
5: C | -0.210000000 | -0.220000000 | -0.044000000
6: C | 0.210000000 | 0.220000000 | 0.044000000
7: O | 0.010000000 | 1.300000000 | 0.045500000
8: O | 0.230000000 | 0.060000000 | 0.055000000
9: O | 0.440000000 | 0.310000000 | 0.055000000
10: Pb | 0.251000000 | 0.574000000 | 0.113500000
11: Pb | -0.251000000 | 1.074000000 | 0.136500000
12: O | -0.230000000 | 0.560000000 | 0.195000000
13: C | -0.210000000 | 0.720000000 | 0.195000000
14: O | -0.010000000 | 0.800000000 | 0.195000000
15: O | -0.440000000 | 0.810000000 | 0.195000000
16: O | 0.010000000 | -0.800000000 | 0.305000000
17: O | 0.440000000 | 0.190000000 | 0.305000000
18: C | 0.210000000 | 0.280000000 | 0.305000000
19: O | 0.230000000 | 0.440000000 | 0.305000000
20: Pb | 0.251000000 | -0.074000000 | 0.363500000
21: Pb | -0.251000000 | -0.574000000 | 0.386500000
22: O | -0.440000000 | -0.310000000 | 0.445000000
23: O | -0.010000000 | -0.300000000 | 0.445000000
24: O | -0.230000000 | -0.060000000 | 0.445000000
25: C | -0.210000000 | -0.220000000 | 0.456000000
26: C | 0.210000000 | 0.220000000 | 0.544000000
27: O | 0.010000000 | 1.300000000 | 0.545500000
28: O | 0.230000000 | 0.060000000 | 0.555000000
29: O | 0.440000000 | 0.310000000 | 0.555000000
30: Pb | 0.251000000 | 0.574000000 | 0.613500000
31: Pb | -0.251000000 | 1.074000000 | 0.636500000
32: O | -0.230000000 | 0.560000000 | 0.695000000
33: C | -0.210000000 | 0.720000000 | 0.695000000
34: O | -0.010000000 | 0.800000000 | 0.695000000
35: O | -0.440000000 | 0.810000000 | 0.695000000
36: O | 0.440000000 | 0.190000000 | 0.805000000
37: C | 0.210000000 | 0.280000000 | 0.805000000
38: O | 0.230000000 | 0.440000000 | 0.805000000
39: Pb | 0.251000000 | -0.074000000 | 0.863500000
40: O | -0.010000000 | -0.300000000 | 0.945000000
***Possible inversion center found at 5-6***
Enter two numbers e.g. "10-20" corresponding to the atomic positions you would like to cleave a surface between:
1-20
O 0.010000000 -0.800000000 -0.195000000
Pb -0.251000000 -0.574000000 -0.113500000
O -0.440000000 -0.310000000 -0.055000000
O -0.230000000 -0.060000000 -0.055000000
C -0.210000000 -0.220000000 -0.044000000
C 0.210000000 0.220000000 0.044000000
O 0.010000000 1.300000000 0.045500000
O 0.230000000 0.060000000 0.055000000
O 0.440000000 0.310000000 0.055000000
Pb 0.251000000 0.574000000 0.113500000
Pb -0.251000000 1.074000000 0.136500000
O -0.230000000 0.560000000 0.195000000
C -0.210000000 0.720000000 0.195000000
O -0.010000000 0.800000000 0.195000000
O -0.440000000 0.810000000 0.195000000
O 0.010000000 -0.800000000 0.305000000
O 0.440000000 0.190000000 0.305000000
C 0.210000000 0.280000000 0.305000000
O 0.230000000 0.440000000 0.305000000
Pb 0.251000000 -0.074000000 0.363500000
Net charge of cleaved surface: 0
Here is the original script:
#!/usr/bin/python3
x=1; y=2; z=3;
#-------------------------------------------------------------------------------------------
#enter input dictionaries here
PbCO3_14 = {
1 : ['Pb',(0.251, 0.574, 0.227)], 2 : [ 'Pb',(-0.251, 1.074, 0.273)],
3 : ['Pb',(-0.251, -0.574, -0.227)], 4 : [ 'Pb',(0.251, -0.074, 0.727)],
5 : ['C',(0.21, 0.22, 0.088)], 6 : [ 'C',(-0.21, 0.72, 0.39)],
7 : ['C',( -0.21, -0.22, -0.088)], 8 : [ 'C',(0.21, 0.28, 0.61)],
9 : ['O',(0.23, 0.06, 0.11)], 10 : [ 'O',(-0.23, 0.56, 0.39)],
11 : ['O',(-0.23, -0.06, -0.11)], 12 : [ 'O',(0.23, 0.44, 0.61)],
13 : ['O',(-0.01, 0.8, 0.39)], 14 : [ 'O',(0.01, 1.3, 0.091)],
15 : ['O',(0.01, -0.8, -0.39)], 16 : [ 'O',(-0.01, -0.30, 0.89)],
17 : ['O',( 0.44, 0.31, 0.11)], 18 : [ 'O',(-0.44, 0.81, 0.39)],
19 : ['O',(-0.44, -0.31, -0.11)], 20 : [ 'O',(0.44, 0.19, 0.61)]
}
#END OF DICTIONARIES
#*****************************************************USER-PARAMETERS*****************************************************
Name = ['PbCO3 #14'] #enter string corresponding to name of structure
output_key = ['Quantum'] #enter either 'Abinit' or 'Quantum' for desired output format of cleaved surface, default is Quantum Espresso format
#dimensions of supercell, 1x1x1 returns the primitive cell
X = 1
Y = 1
Z = 2
Sort_by = [z,y,x] #sorting priority for supercell coordinates - must be some permutation of x, y, and z
charges = { 'H' : +1, 'Pb' : +2, 'O' : -2, 'C' : +4 }
#shifts for x y and z coordinates of the cell - use to find inversion plane
x_shift = 0
y_shift = 0
z_shift = 0
STRUCTURE = PbCO3_14 #name of input dictionary for primitive cell atoms and coordinates
#****************************************************END-OF-PARAMETERS****************************************************
def Expand(prim_cell,Atoms,prim,X,Y,Z,x_shift,y_shift,z_shift):
Atom = [ prim[1][0]/X, prim[1][1]/Y, prim[1][2]/Z ]
L_X = [ [ Atom[0] + i/X, Atom[1], Atom[2] ] for i in range(X) ];
L_XY = [ [ Coord[0], Coord[1] + i/Y, Coord[2] ] for Coord in L_X for i in range(Y) ]
L_XYZ = [ [ prim[0], round(Coord[0] + x_shift,9), round(Coord[1] + y_shift,9),
round(Coord[2] + i/Z + z_shift, 9) ] for Coord in L_XY for i in range(Z) ]
for i in range(len(L_XYZ)): Atoms = { max( Atoms, key=int ) + 1 : L_XYZ[i], **Atoms }
return Atoms
#-------------------------------------------------------------------------------------------
def Layer(prim_cell,Atoms,Layers):
inversion = False; origin = False; k_inv = False; k_zero = False;
Layered_Atoms = sorted(Atoms.items(), key = lambda x:(x[1][Sort_by[0]],x[1][Sort_by[1]],x[1][Sort_by[2]]))
for i in range(len(Layered_Atoms)):[ Layers.update( { i + 1 : Layered_Atoms[i][1] } ) for i in range(len(Layered_Atoms)) ]
for key, value in Layers.items():
if (key != list(Layers)[-1]) and (Layers[key+1] == [ value[0],-value[1],-value[2],-value[3] ]):inversion = True; k_inv = key
if value[1] == value[2] == value[3] == 0.0: origin = True; k_zero = key
print('{:4} {:2} {:^12} {:^12} {:^12}\n'.format('','','x','y','z'))
for key, value in Layers.items():
if ((inversion == True ) and (key in [k_inv,k_inv+1])):print(fmt.inversion.format(key,*value));
if ((inversion == True) and (key not in [k_inv,k_inv+1])) or inversion == False: print(fmt.main.format(key,*value))
if inversion == True: print((color.BOLD + '\n***Possible inversion center found at {}-{}***\n' + color.END).format(k_inv,k_inv+1))
if origin == True: print('\n***Possible inversion center found at {}**\n'.format(k_zero))
return Layers
#-------------------------------------------------------------------------------------------
def Supercell(prim_cell,Atoms,Layers,X,Y,Z,x_shift,y_shift,z_shift):
for i in range(len(prim_cell)): Atoms = Expand(prim_cell,Atoms,prim_cell[i+1],*params()[1],*params()[2])
del Atoms[0]
Layers = Layer(prim_cell,Atoms,Layers)
print('\nEnter two numbers e.g. "10-20" corresponding to the atomic positions you would like to cleave a surface between')
Surface = input('\n: ').split('-');print('\n');
if Surface != [''] and (len(Surface) == 2):
[ Cleaved_Surface.update( { i + 1 : [ *Layers[i] ] } ) for i in range(int(Surface[0]),int(Surface[1])+1) ]
for i in range(len(Cleaved_Surface)):
if output_key[0] == 'Abinit':
Atom_Sort = sorted(Cleaved_Surface.items(), key = lambda x: (x[1][0],x[1][Sort_by[0]]));
print(fmt.abinit.format(Atom_Sort[i][1][1],Atom_Sort[i][1][2],Atom_Sort[i][1][3],Atom_Sort[i][1][0]))
if output_key[0] in ['','Quantum'] :
Atom_Sort = sorted(Cleaved_Surface.items(), key = lambda x: (x[1][Sort_by[0]],x[1][Sort_by[1]],x[1][Sort_by[2]]))
print(fmt.quantum.format(*Atom_Sort[i][1]))
Net_charge(charges,Cleaved_Surface)
#-------------------------------------------------------------------------------------------
def params():
Dim = [X,Y,Z]
if [(isinstance(dim,int)) for dim in Dim] != [True]*3: print('\nCheck supercell dimension parameters.\n');raise SystemExit;
Shift = [x_shift,y_shift,z_shift]
for shift in Shift:
if isinstance(shift,float) == False and shift != 0: print('\nInvalid shift parameters.\n'); raise SystemExit
key_map = { 1: 'x', 2: 'y', 3: 'z' }; Sorter = [key_map.get(key, 'Null') for key in Sort_by]
if [Sort_by.count(sort) for sort in Sort_by] != [1]*3: print('Check Sort_by parameters.'); raise SystemExit;
if output_key[0] not in ['','Quantum','Abinit']: print('Check output_key.'); raise SystemExit
return Sorter,Dim,Shift
#-------------------------------------------------------------------------------------------
def Net_charge(charge,cleaved_surface):
net_charge = int(sum([charge.get(v[0], 'Null') for k,v in cleaved_surface.items()]));
def colorize(col,net_charge):
print(fmt.charge.format(color.BOLD,col,net_charge,color.END))
if net_charge > 0: colorize(color.GREEN,'+'+str(net_charge))
elif net_charge == 0: colorize(color.CYAN,net_charge)
else: colorize(color.RED,net_charge)
#-------------------------------------------------------------------------------------------
class color:
CYAN ='\033[96m';BOLD ='\033[1m';END ='\033[0m';HIGHLIGHT ='\033[01;97;105m';
FLASH ='\033[5m';PINK ='\033[95m';GREEN ='\033[32m';RED ='\033[31m';
class fmt:
main = '{:4}: {:2} | {:12.9f} | {:12.9f} | {:12.9f}';abinit = ' {:12.9f} {:12.9f} {:12.9f} #{:3}';
quantum = '{:2} {:12.9f} {:12.9f} {:12.9f}';charge = '\nNet charge of cleaved surface: {}{}{}{}\n';
inversion = color.HIGHLIGHT + '{:4}: {:2} | {:12.9f} | {:12.9f} | {:12.9f}' + color.END;
#-------------------------------------------------------------------------------------------
if __name__ == "__main__":
Atoms = { 0 : 'Null' }; Layers = {}; Cleaved_Surface = {};
if X == Y == Z == 1: print('{}\n{} unaltered cell\n'.format('\n'+'*'*55+'\n',Name[0]))
else: print('{}\n{} x {} x {} {} supercell\n'.format('\n'+'*'*55+'\n',*params()[1],Name[0]))
print('x-shift: {} \ny-shift: {} \nz-shift: {} \n\nSorting priority: {}-{}-{}\n{}'.format(*params()[2],*params()[0],'\n'+'*'*55+'\n'))
Supercell(STRUCTURE,Atoms,Layers,*params()[1],*params()[2])
I'm inexperienced at programming and wanted to teach myself about classes and basic methods by incorporating them into this script in a meaningful way, and so this is how I updated it:
#!/usr/bin/python3
def userParam(x,y,z):
#USER DICTIONARIES
PbCO3_14 = {
1 : ['Pb',(0.251, 0.574, 0.227)], 2 : [ 'Pb',(-0.251, 1.074, 0.273)],
3 : ['Pb',(-0.251, -0.574, -0.227)], 4 : [ 'Pb',(0.251, -0.074, 0.727)],
5 : ['C',(0.21, 0.22, 0.088)], 6 : [ 'C',(-0.21, 0.72, 0.39)],
7 : ['C',( -0.21, -0.22, -0.088)], 8 : [ 'C',(0.21, 0.28, 0.61)],
9 : ['O',(0.23, 0.06, 0.11)], 10 : [ 'O',(-0.23, 0.56, 0.39)],
11 : ['O',(-0.23, -0.06, -0.11)], 12 : [ 'O',(0.23, 0.44, 0.61)],
13 : ['O',(-0.01, 0.8, 0.39)], 14 : [ 'O',(0.01, 1.3, 0.091)],
15 : ['O',(0.01, -0.8, -0.39)], 16 : [ 'O',(-0.01, -0.30, 0.89)],
17 : ['O',( 0.44, 0.31, 0.11)], 18 : [ 'O',(-0.44, 0.81, 0.39)],
19 : ['O',(-0.44, -0.31, -0.11)], 20 : [ 'O',(0.44, 0.19, 0.61)]
}
#END OF DICTIONARIES
#*****************************************************USER-PARAMETERS*****************************************************
STRUCTURE = PbCO3_14 # name of input dictionary for primitive cell atoms and coordinates
Name = ['PbCO3 #14'] # enter string corresponding to name of structure
# dimensions of supercell, 1x1x1 returns the primitive cell
X = 1
Y = 1
Z = 2
sortBy = [z,y,x] #sorting priority for supercell coordinates - must be some permutation of x, y, and z
charges = { 'H' : +1, 'Pb' : +2, 'O' : -2, 'C' : +4 }
# shifts for x y and z coordinates of the cell - use to find inversion plane
x_shift = 0
y_shift = 0
z_shift = 0
outputKey = ['Quantum'] # enter either 'Abinit' or 'Quantum' for desired output format of cleaved surface, default is Quantum Espresso format
#****************************************************END-OF-PARAMETERS****************************************************
return [STRUCTURE,Name,outputKey,[X,Y,Z],sortBy,[x_shift,y_shift,z_shift],charges]
#-------------------------------------------------------------------------------------------
class color:
CYAN ='\033[96m';BOLD ='\033[1m';END ='\033[0m';HIGHLIGHT ='\033[01;97;105m';
FLASH ='\033[5m';PINK ='\033[95m';GREEN ='\033[32m';RED ='\033[31m';
class fmt:
main = '{:4}: {:2} | {:12.9f} | {:12.9f} | {:12.9f}';
abinit = ' {:12.9f} {:12.9f} {:12.9f} #{:3}';
quantum = '{:2} {:12.9f} {:12.9f} {:12.9f}';
charge = '\nNet charge of cleaved surface: {}{}{}{}\n';
inversion = color.HIGHLIGHT + '{:4}: {:2} | {:12.9f} | {:12.9f} | {:12.9f}' + color.END;
header = '{:4} {:2} {:^12} {:^12} {:^12}\n'
showParam = 'x-shift: {:4.3f} \ny-shift: {:4.3f} \nz-shift: {:4.3f} \n\nSorting priority: {}-{}-{}\n{}'
class Supercell:
def __init__(self,structure,name,outputKey,dim,sortBy,shift,charges):
self.primitive = structure;
self.name = name[0]
self.format = outputKey[0]
self.sortBy = sortBy
self.dim = dim
self.X = dim[0]; self.Y = dim[1]; self.Z = dim[2];
self.shift = shift
self.x_shift = shift[0]; self.y_shift = shift[1]; self.z_shift = shift[2]
self.charges = charges
self.inversion = False; self.inversionKey = False
self.origin = False; self.zeroKey = False
self.layers = {}
self.cleaved = {}
#-------------------------------------------------------------------------------------------
def checkParam(self):
if [(isinstance(dim,int)) for dim in self.dim] != [True]*3: print('\nCheck supercell dimension parameters.\n');raise SystemExit;
for shift in self.shift:
if isinstance(shift,float) == False and shift != 0: print('\nInvalid shift parameters.\n'); raise SystemExit
if [self.sortBy.count(sort) for sort in self.sortBy] != [1]*3: print('Check sortBy parameters.'); raise SystemExit;
keyMap = { 1: 'x', 2: 'y', 3: 'z' }; sortKeys = [keyMap.get(key, 'Null') for key in self.sortBy]
if self.format not in ['','Quantum','Abinit']: print('Check outputKey.'); raise SystemExit
return sortKeys
#-------------------------------------------------------------------------------------------
def constructCell(self):
emptyCell = { 0 : 'Null' }
for i in range(len(self.primitive)):
atom = self.primitive[i+1][0]
basis = [ self.primitive[i+1][1][0]/self.X, self.primitive[i+1][1][1]/self.Y, self.primitive[i+1][1][2]/self.Z ]
X_RED = [ [ basis[0] + i/self.X, basis[1], basis[2] ] for i in range(self.X) ];
XY_RED = [ [ xred[0], xred[1] + i/self.Y, xred[2] ] for xred in X_RED for i in range(self.Y) ]
XYZ_RED = [ [ atom, round(xred[0] + self.x_shift,9), round(xred[1] + self.y_shift,9), round(xred[2] + i/self.Z + self.z_shift, 9) ] for xred in XY_RED for i in range(self.Z) ]
for i in range(len(XYZ_RED)): emptyCell = { max( emptyCell, key=int ) + 1 : XYZ_RED[i], **emptyCell }
del emptyCell[0]; return emptyCell
#-------------------------------------------------------------------------------------------
def layerCell(self):
cell = self.constructCell()
layeredCell = sorted(cell.items(), key = lambda x:(x[1][self.sortBy[0]],x[1][self.sortBy[1]],x[1][self.sortBy[2]]))
for i in range(len(layeredCell)):[ self.layers.update( { i + 1 : layeredCell[i][1] } ) for i in range(len(layeredCell)) ]
return self.layers
#-------------------------------------------------------------------------------------------
def displayParam(self):
sortKeys = self.checkParam()
if self.X == self.Y == self.Z == 1: print('{}\n{} unaltered cell\n'.format('\n'+'*'*55+'\n',self.name))
else: print('{}\n{} x {} x {} {} supercell\n'.format('\n'+'*'*55+'\n',*self.dim,self.name))
print(fmt.showParam.format(*self.shift,*sortKeys,'\n'+'*'*55+'\n'))
#-------------------------------------------------------------------------------------------
def displayCell(self):
self.displayParam()
layers = self.layerCell()
for key, value in layers.items():
if (key != list(layers)[-1]) and (layers[key+1] == [value[0],-value[1],-value[2],-value[3]]): self.inversion = True; self.inversionKey = key
if value[1] == value[2] == value[3] == 0.0: self.origin = True; self.zeroKey = key
print(fmt.header.format('','','x','y','z'))
for key, value in layers.items():
if ((self.inversion == True ) and (key in [self.inversionKey,self.inversionKey+1])):print(fmt.inversion.format(key,*value));
if ((self.inversion == True) and (key not in [self.inversionKey,self.inversionKey+1])) or self.inversion == False: print(fmt.main.format(key,*value))
if self.inversion == True: print((color.BOLD + '\n***Possible inversion center found at {}-{}***\n' + color.END).format(self.inversionKey,self.inversionKey+1))
if self.origin == True: print('\n***Possible inversion center found at {}**\n'.format(self.zeroKey))
while True:
Surface = input('\nEnter two numbers e.g. "10-20" corresponding to the atomic positions you would like to cleave a surface between:\n\n').split('-');print('\n')
if Surface != [''] and (len(Surface) == 2):
[ self.cleaved.update( { i + 1 : [ *layers[i] ] } ) for i in range(int(Surface[0]),int(Surface[1])+1) ]
else: print('Exiting...\n'); raise SystemExit
for i in range(len(self.cleaved)):
if self.format == 'Abinit':
cellSort = sorted(self.cleaved.items(), key = lambda x: (x[1][0],x[1][self.sortBy[0]]));
print(fmt.abinit.format(cellSort[i][1][1],cellSort[i][1][2],cellSort[i][1][3],cellSort[i][1][0]))
if self.format in ['','Quantum']:
cellSort = sorted(self.cleaved.items(), key = lambda x: (x[1][self.sortBy[0]],x[1][self.sortBy[1]],x[1][self.sortBy[2]]))
print(fmt.quantum.format(*cellSort[i][1]))
netCharge = int(sum([self.charges.get(data[0], 'Null') for data in self.cleaved.values()]));
if netCharge > 0: print(fmt.charge.format(color.BOLD,color.GREEN,'+'+str(netCharge),color.END))
elif netCharge == 0: print(fmt.charge.format(color.BOLD,color.CYAN,netCharge,color.END))
else: print(fmt.charge.format(color.BOLD,color.RED,netCharge,color.END))
print('Press enter to exit')
self.cleaved = {}
#-------------------------------------------------------------------------------------------
if __name__ == "__main__":
Supercell(*userParam(1,2,3)).displayCell()
Was this an improvement? Any criticisms are welcome.