mirror of
https://github.com/NehZio/Crystal-MEC
synced 2024-09-27 12:01:01 +02:00
479 lines
19 KiB
Python
479 lines
19 KiB
Python
import numpy as np
|
|
import os
|
|
import operator
|
|
|
|
|
|
###### CONSTANTS #######
|
|
|
|
a = 5.48463
|
|
b = 5.48463
|
|
c = 25.7977
|
|
|
|
alpha = 90
|
|
beta = 90
|
|
gamma = 90
|
|
|
|
dr = 0.001
|
|
da = 0.01
|
|
DA = 1.5
|
|
|
|
|
|
chO = -2.00
|
|
chIr = 4.00
|
|
chSr = 2.00
|
|
dIrO = 2.50
|
|
dSrO = 3.10
|
|
#######################
|
|
|
|
def parse(fileName):
|
|
f = open(fileName,'r')
|
|
line = 0
|
|
dataTable = []
|
|
|
|
f.readline()
|
|
f.readline()
|
|
|
|
while True:
|
|
line = f.readline().strip().split()
|
|
if line == [] :
|
|
break
|
|
line[1] = float(line[1])
|
|
line[2] = float(line[2])
|
|
line[3] = float(line[3])
|
|
|
|
dataTable.append(line)
|
|
|
|
return dataTable
|
|
|
|
def distance(u,v):
|
|
return np.sqrt((u[0]-v[0])**2+(u[1]-v[1])**2+(u[2]-v[2])**2)
|
|
|
|
def vect_product(u,v):
|
|
return [u[1]*v[2]-u[2]*v[1],u[2]*v[0]-u[0]*v[2],u[0]*v[1]-u[1]*v[0]]
|
|
|
|
def dot_product(u,v):
|
|
return (u[0]*v[0]+u[1]*v[1]+u[2]*v[2])
|
|
|
|
def normalize(u):
|
|
norm = np.sqrt(u[0]**2+u[1]**2+u[2]**2)
|
|
return [u[0]/norm, u[1]/norm, u[2]/norm]
|
|
|
|
def rot_matrix(oldAxis,newAxis):
|
|
newAxis = normalize(newAxis)
|
|
vp = normalize(vect_product(oldAxis,newAxis))
|
|
angle = np.arccos(dot_product(oldAxis,newAxis))
|
|
|
|
rMat = [[np.cos(angle)+vp[0]*vp[0]*(1-np.cos(angle)),vp[0]*vp[1]*(1-np.cos(angle))-vp[2]*np.sin(angle),vp[0]*vp[2]*(1-np.cos(angle))+vp[1]*np.sin(angle)],[vp[0]*vp[1]*(1-np.cos(angle))+vp[2]*np.sin(angle),vp[1]*vp[1]*(1-np.cos(angle))+np.cos(angle),vp[1]*vp[2]*(1-np.cos(angle))-vp[0]*np.sin(angle)],[vp[0]*vp[2]*(1-np.cos(angle))-vp[1]*np.sin(angle),vp[1]*vp[2]*(1-np.cos(angle))+vp[0]*np.sin(angle),np.cos(angle)+vp[2]*vp[2]*(1-np.cos(angle))]]
|
|
|
|
##This is the formula for the rotation matrix to change axis
|
|
|
|
return rMat
|
|
|
|
def rotation(coord, rMat):
|
|
newCoord = []
|
|
|
|
for i in range(len(coord)):
|
|
newX = coord[i][0]*rMat[0][0] + coord[i][1]*rMat[1][0] + coord[i][2]*rMat[2][0]
|
|
newY = coord[i][0]*rMat[0][1] + coord[i][1]*rMat[1][1] + coord[i][2]*rMat[2][1]
|
|
newZ = coord[i][0]*rMat[0][2] + coord[i][1]*rMat[1][2] + coord[i][2]*rMat[2][2]
|
|
newCoord.append([newX,newY,newZ])
|
|
|
|
return newCoord
|
|
|
|
def cut_bath(rBath, coords):
|
|
bath = []
|
|
for i in range(len(coords)):
|
|
if distance([0,0,0],[coords[i][0],coords[i][1],coords[i][2]]) <= rBath:
|
|
bath.append([coords[i][0], coords[i][1], coords[i][2], coords[i][3], 'C'])
|
|
|
|
return bath
|
|
|
|
def set_pp(rPP,coords, notIn):
|
|
|
|
pp = []
|
|
|
|
for i in range(len(coords)):
|
|
for j in range(len(coords)):
|
|
if coords[i][4] == 'O':
|
|
if coords[j][4] == 'C' and coords[j][3] not in notIn:
|
|
if distance([coords[i][0],coords[i][1],coords[i][2]],[coords[j][0],coords[j][1],coords[j][2]]) <= rPP:
|
|
coords[j][4] == 'Cl'
|
|
pp.append(j)
|
|
for i in pp:
|
|
coords[i][4] = 'Cl'
|
|
return coords
|
|
|
|
def find_frag(pattern, n, coords): #We mark the atoms in the bath corresponding to
|
|
#the fragment according to user input
|
|
inFrag = []
|
|
for k in range(n):
|
|
closest = [100,100,100]
|
|
for j in coords:
|
|
if j[3] == pattern[1]:
|
|
if distance(j,[0,0,0]) < distance([0,0,0],closest) and [j[0],j[1],j[2],distance(j,j), coords.index(j)] not in inFrag:
|
|
closest = [j[0],j[1],j[2],distance(j,j), coords.index(j)]
|
|
for i in range(1,len(pattern)/2):
|
|
inPattern = [closest]
|
|
for j in range(int(pattern[2*i])):
|
|
inPattern.append([100,100,100,distance([100,100,100],closest)])
|
|
for j in coords:
|
|
inPattern = sorted(inPattern,key=operator.itemgetter(3))
|
|
if j[3] == pattern[2*i+1]:
|
|
if distance(j,closest) <= inPattern[-1][3]:
|
|
inPattern[-1] = [j[0],j[1],j[2], distance(j,closest), coords.index(j)]
|
|
for j in inPattern:
|
|
inFrag.append(j)
|
|
for j in inFrag:
|
|
coords[j[4]][4] = 'O'
|
|
|
|
return coords
|
|
|
|
def eivgen(coords):
|
|
charges = []
|
|
for i in coords:
|
|
nIr = 0
|
|
nSr = 0
|
|
nO = 0
|
|
for j in coords:
|
|
if i[3] == 'Ir' and j[3] == 'O' and distance(i,j) < dIrO:
|
|
nO += 1
|
|
elif i[3] == 'Sr' and j[3] == 'O' and distance(i,j) < dSrO:
|
|
nO += 1
|
|
if i[3] == 'Ir':
|
|
charges.append(chIr*(nO/6.00))
|
|
if i[3] == 'Sr':
|
|
charges.append(chSr*(nO/9.00))
|
|
|
|
typeO = 0
|
|
|
|
if i[3] == 'O':
|
|
for j in coords:
|
|
if j[3] == 'Ir' and distance(i,j) < dIrO:
|
|
nIr += 1
|
|
if distance(i,j) < 2.00:
|
|
typeO = 1
|
|
elif j[3] == 'Sr' and distance(i,j) < dSrO:
|
|
nSr += 1
|
|
if typeO == 0:
|
|
charges.append(chO*((chSr*nSr+chIr*nIr)/(5*chSr+chIr)))
|
|
else:
|
|
charges.append(chO*((chSr*nSr+chIr*nIr)/(4*chSr+2*chIr)))
|
|
|
|
|
|
|
|
|
|
coords = [[coords[i][0],coords[i][1],coords[i][2],coords[i][3],coords[i][4],charges[i]] for i in range(len(coords))]
|
|
|
|
return coords
|
|
|
|
def symmetry(coord,atoms,charges, operations):
|
|
newCoord = []
|
|
|
|
while coord != []:
|
|
toDel = []
|
|
newCoord.append(coord[0]) #Add the atom to a new list
|
|
name = atoms[0]
|
|
a = newCoord[-1][:] #label a = E
|
|
b = [-a[0],a[1],a[2]] #label b = yOz mirror plan
|
|
c = [-a[0],-a[1],a[2]] #label c = C2 rotation around z
|
|
d = [a[0],-a[1],a[2]] #label d = xOz mirror plan
|
|
e = [a[0],a[1],-a[2]] #label e = xOy mirror plan
|
|
f = [-a[0],a[1],-a[2]] #label f = C2 rotation around y
|
|
g = [a[0],-a[1],-a[2]] #label g = C2 rotation around x
|
|
h = [-a[0],-a[1],-a[2]] #label h = i
|
|
newCoord[-1].append(name+"a")
|
|
newCoord[-1].append(charges[0])
|
|
del atoms[0] #Delete from old list
|
|
del coord[0]
|
|
del charges[0]
|
|
|
|
for t in coord:
|
|
index = coord.index(t)
|
|
if name == atoms[index]: #Check if it is the same atom
|
|
if distance(t,a) == 0:
|
|
print("ERROR : Twice the same atom")
|
|
if 'xOz' in operations:
|
|
if distance(t,d) <= da:
|
|
newCoord.append(d)
|
|
newCoord[-1].append(name+'d')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,d) and distance(t,d) < DA:
|
|
print "Error : This atom should not be there",distance(t,d),t,d,a,charges[index]
|
|
print "Are you sure about the xOz symmetry operation ?"
|
|
break
|
|
if 'yOz' in operations:
|
|
if distance(t,b) <= da:
|
|
newCoord.append(b)
|
|
newCoord[-1].append(name+'b')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,b) and distance(t,b) < DA:
|
|
print "Error : This atom should not be there",distance(t,b),t,b,a,charges[index]
|
|
print "Are you sure about the yOz symmetry operation ?"
|
|
break
|
|
if 'C2z' in operations:
|
|
if distance(t,c) <= da:
|
|
newCoord.append(c)
|
|
newCoord[-1].append(name+'c')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,c) and distance(t,c) < DA:
|
|
print "Error : This atom should not be there",distance(t,c),t,c,a,charges[index]
|
|
print "Are you sure about the C2z axis ?"
|
|
break
|
|
if 'xOy' in operations:
|
|
if distance(t,e) <= da:
|
|
newCoord.append(e)
|
|
newCoord[-1].append(name+'e')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,e) and distance(t,e) < DA:
|
|
print "Error : This atom should not be there",distance(t,e),t,e,a,charges[index]
|
|
print "Are you sure about the xOy operation ?"
|
|
break
|
|
if 'C2y' in operations:
|
|
if distance(t,f) <= da:
|
|
newCoord.append(f)
|
|
newCoord[-1].append(name+'f')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,f) and distance(t,f) < DA:
|
|
print "Error : This atom should not be there",distance(t,f),t,f,a,charges[index]
|
|
print "Are you sure about the C2y axis ?"
|
|
break
|
|
if 'C2x' in operations:
|
|
if distance(t,g) <= da:
|
|
newCoord.append(g)
|
|
newCoord[-1].append(name+'g')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,g) and distance(t,g) < DA:
|
|
print "Error : This atom should not be there",distance(t,g),t,g,a,charges[index]
|
|
print "Are you sure about the C2x axis ?"
|
|
break
|
|
if 'i' in operations:
|
|
if distance(t,h) <= da:
|
|
newCoord.append(h)
|
|
newCoord[-1].append(name+'h')
|
|
newCoord[-1].append(charges[index])
|
|
toDel.append(index)
|
|
elif da < distance(t,h) and distance(t,h) < DA:
|
|
print "Error : This atom should not be there",distance(t,h),t,h,a,charges[index]
|
|
print "Are you sure about the i operation ?"
|
|
break
|
|
|
|
|
|
for m in range(len(toDel)): #We delete the atoms seen in the simmetry
|
|
del coord[toDel[m]-m]
|
|
del atoms[toDel[m]-m]
|
|
del charges[toDel[m]-m]
|
|
|
|
return newCoord
|
|
|
|
def write_input(fragCoord,ppCoord,bathCoord,fileName, sym):
|
|
g = open(fileName,'w')
|
|
if sym == 'y':
|
|
for i in range(len(fragCoord)):
|
|
if fragCoord[i][3][-1] == 'a':
|
|
g.write('%4s % 5.3f % 5.3f % 5.3f angstrom\n' %(fragCoord[i][3].replace('a','')+str(i+1),fragCoord[i][0],fragCoord[i][1],fragCoord[i][2]))
|
|
g.write("\n\n")
|
|
for i in range(len(ppCoord)):
|
|
if ppCoord[i][3][-1] == 'a':
|
|
g.write('%4s % 5.3f % 5.3f % 5.3f angstrom\n' %(ppCoord[i][3].replace('a','')+str(i+1),ppCoord[i][0],ppCoord[i][1],ppCoord[i][2]))
|
|
g.write("\n\n")
|
|
for i in range(len(bathCoord)):
|
|
if bathCoord[i][3][-1] == 'a':
|
|
g.write('% 5.3f % 5.3f % 5.3f % 5.3f 0. 0. 0. \n'%(bathCoord[i][0],bathCoord[i][1],bathCoord[i][2],bathCoord[i][4]))
|
|
if sym == 'n':
|
|
for i in range(len(fragCoord)):
|
|
g.write('%4s % 5.3f % 5.3f % 5.3f angstrom\n' %(fragCoord[i][3]+str(i+1),fragCoord[i][0],fragCoord[i][1],fragCoord[i][2]))
|
|
g.write('\n\n')
|
|
for i in range(len(ppCoord)):
|
|
g.write('%4s % 5.3f % 5.3f % 5.3f angstrom\n' %(ppCoord[i][3]+str(i+1),ppCoord[i][0],ppCoord[i][1],ppCoord[i][2]))
|
|
g.write('\n\n')
|
|
for i in range(len(bathCoord)):
|
|
g.write('% 5.3f % 5.3f % 5.3f % 5.3f 0. 0. 0. \n' %(bathCoord[i][0],bathCoord[i][1],bathCoord[i][2], bathCoord[i][4]))
|
|
g.close()
|
|
|
|
def translation(vec, coords):
|
|
for i in range(len(coords)):
|
|
coords[i][0] -= vec[0]
|
|
coords[i][1] -= vec[1]
|
|
coords[i][2] -= vec[2]
|
|
return coords
|
|
|
|
def main():
|
|
|
|
fileName = raw_input("Name of the cif file : ")
|
|
rBath = input("Chose the bath radius (in Angstrom) : ")
|
|
rPP = input("Chose the pseudopotential radius (in Angstrom) : ")
|
|
notIn = raw_input("What atoms should not be in the first shell ? ").split()
|
|
|
|
nA = int(np.floor(2*rBath/a)+2) #We chose the number of time we need to replicate
|
|
nB = int(np.floor(2*rBath/(b*np.sin(np.radians(gamma))))+2) #to be able to cut the bath
|
|
nC = int(np.floor(2*rBath/(c*np.sin(np.radians(beta))))+2)
|
|
|
|
cmd = 'atomsk '+ fileName + ' ' + '-duplicate ' + str(nA) + ' ' + str(nB) + ' ' + str(nC) + ' ' + fileName.replace('cif','xyz') + ' -v 0' #This is the command that calls the program that generates the big cell
|
|
os.system(cmd)
|
|
|
|
|
|
data = parse(fileName.replace('cif','xyz')) #Read the data from the xyz file
|
|
|
|
coords = [[data[i][1],data[i][2],data[i][3]] for i in range(len(data))]
|
|
labels = [data[i][0] for i in range(len(data))]
|
|
|
|
coords = [[coords[i][0],coords[i][1],coords[i][2],labels[i]] for i in range(len(coords))]
|
|
|
|
coords = translation([nA*a/2,nB*b/2,nC*c/2],coords) #Putting the origin at the center of the cell
|
|
|
|
|
|
labels = [i[3] for i in coords]
|
|
|
|
center = raw_input("Specify the center of the fragment (if between atoms, specify them all) ").split() #Searching for the center according to user
|
|
centers = [] #input, and translating the coordinates
|
|
for i in range(len(center)):
|
|
centers.append([100,100,100])
|
|
|
|
for i in centers:
|
|
i.append(distance([0,0,0],i))
|
|
|
|
for i in coords:
|
|
centers = sorted(centers,key=operator.itemgetter(3))
|
|
if i[3] in center:
|
|
if distance(i,[0,0,0]) <= centers[-1][3]:
|
|
centers[-1] = [i[0],i[1],i[2],distance(i,[0,0,0])]
|
|
|
|
newOgn = np.mean(np.array(centers),axis=0)
|
|
newOgn = [newOgn[0], newOgn[1], newOgn[2]]
|
|
|
|
coords = translation(newOgn,coords)
|
|
|
|
axis = ['x','y','z']
|
|
|
|
for k in axis:
|
|
nAxis = raw_input("Where should the %s axis be headed ? (if between atoms, specify them all), x if none "%k).split() #Searching for the new orientation
|
|
nAxiss = [] #according to user input, and
|
|
#rotating the coordinates
|
|
if nAxis[0] == 'x':
|
|
continue
|
|
for i in range(len(nAxis)):
|
|
nAxiss.append([100,100,100])
|
|
|
|
for i in nAxiss:
|
|
i.append(distance([0,0,0],i))
|
|
|
|
for i in coords:
|
|
nAxiss = sorted(nAxiss,key=operator.itemgetter(3))
|
|
if i[3] in nAxis:
|
|
if distance(i,[0,0,0]) <= nAxiss[-1][3]:
|
|
nAxiss[-1] = [i[0],i[1],i[2],distance(i,[0,0,0])]
|
|
newN = np.mean(np.array(nAxiss),axis=0)
|
|
newN = [newN[0],newN[1],newN[2]]
|
|
|
|
if k == 'x':
|
|
oldN = [1,0,0]
|
|
elif k == 'y':
|
|
oldN = [0,1,0]
|
|
elif k == 'z':
|
|
oldN = [0,0,1]
|
|
|
|
rMat = rot_matrix(oldN,newN)
|
|
coords = rotation(coords,rMat)
|
|
coords = [[coords[i][0],coords[i][1],coords[i][2],labels[i]] for i in range(len(coords))]
|
|
|
|
#We now have one big cell oriented and centered as we want
|
|
#The rest of the code will cut what we want in this big cell
|
|
|
|
pattern = raw_input("What pattern is in the fragment ? ").split()
|
|
|
|
n = input("How many time does the pattern repeat in the fragment ? ")
|
|
|
|
|
|
|
|
coords = cut_bath(rBath,coords)
|
|
coords = find_frag(pattern, n ,coords)
|
|
coords = set_pp(rPP,coords,notIn)
|
|
|
|
|
|
coords = eivgen(coords)
|
|
ch = 0
|
|
for i in range(len(coords)):
|
|
ch += coords[i][5]
|
|
print("Total charge : % 8.5f"%ch)
|
|
sym = 'x'
|
|
|
|
os.system('rm '+fileName.replace('cif','xyz'))
|
|
|
|
frag = sorted([i for i in coords if i[4] == 'O'],key=operator.itemgetter(3))
|
|
pp = sorted([i for i in coords if i[4] == 'Cl'],key=operator.itemgetter(3))
|
|
bath = sorted([i for i in coords if i[4] == 'C'],key=operator.itemgetter(3))
|
|
|
|
if raw_input("Do you want to see the fragment ? (y/n) ") == 'y':
|
|
g = open('tmp.xyz','w')
|
|
g.write('%i \n \n'%len(frag))
|
|
for j in frag:
|
|
g.write('%s % 6.2f % 6.2f % 6.2f \n'%(j[3],j[0],j[1],j[2]))
|
|
g.close()
|
|
os.system('avogadro tmp.xyz')
|
|
#os.system('rm tmp.xyz')
|
|
|
|
while sym != 'y' and sym != 'n':
|
|
sym = raw_input("Do you want to treat the symmetry ? (y/n) ")
|
|
|
|
if sym == 'y':
|
|
operation = raw_input("What are the symmetry operations you'd like to treat ? (C2(x,y,z), xOy xOz yOz, i) ").split()
|
|
rep = 0
|
|
for i in range(len(coords)-1):
|
|
for j in range(i+1,len(coords)):
|
|
rep += (coords[i][5]*coords[j][5])/distance(coords[i],coords[j])
|
|
print("Nuclear repulsion before symmetry : %f"%rep)
|
|
|
|
frag = symmetry([[i[0],i[1],i[2]] for i in frag],[i[3] for i in frag], [i[5] for i in frag], operation)
|
|
pp = symmetry([[i[0],i[1],i[2]] for i in pp],[i[3] for i in pp], [i[5] for i in pp], operation)
|
|
bath = symmetry([[i[0],i[1],i[2]] for i in bath],[i[3] for i in bath], [i[5] for i in bath], operation)
|
|
|
|
coords = frag+pp+bath
|
|
|
|
rep = 0
|
|
for i in range(len(coords)-1):
|
|
for j in range(i+1,len(coords)):
|
|
rep += (coords[i][4]*coords[j][4])/distance(coords[i],coords[j])
|
|
|
|
print("Nuclear repulsion after symmetry : %f"%rep)
|
|
|
|
if sym == 'n':
|
|
frag = [[i[0],i[1],i[2],i[3],i[5]] for i in frag]
|
|
pp = [[i[0],i[1],i[2],i[3],i[5]] for i in pp]
|
|
bath = [[i[0],i[1],i[2],i[3],i[5]] for i in bath]
|
|
fileName = raw_input("What name do you wish yo give to your file ? ")
|
|
|
|
write_input(frag,pp,bath,fileName,sym)
|
|
|
|
if raw_input("Do you want to visualize the bath ? (y/n) ") == 'y':
|
|
g = open('tmp.xyz','w')
|
|
if raw_input("Do you want it to be sorted with colors ? (y/n) ") == 'y':
|
|
g.write('%i \n\n'%(len(frag)+len(pp)+len(bath)))
|
|
for i in frag:
|
|
g.write('O % 5.3f % 5.3f % 5.3f \n'%(i[0],i[1],i[2]))
|
|
for i in pp:
|
|
g.write('Cl % 5.3f % 5.3f % 5.3f \n'%(i[0],i[1],i[2]))
|
|
for i in bath:
|
|
g.write('C % 5.3f % 5.3f % 5.3f \n'%(i[0],i[1],i[2]))
|
|
else:
|
|
coords = frag+pp+bath
|
|
if sym == 'y':
|
|
for i in coords:
|
|
l = ['a','b','c','d','e','f','g','h']
|
|
if i[3][-1] in l:
|
|
i[3] = i[3][:-1]
|
|
g.write('%i \n \n'%len(coords))
|
|
for i in coords:
|
|
g.write('%s % 5.3f % 5.3f % 5.3f \n'%(i[3],i[0],i[1],i[2]))
|
|
g.close()
|
|
os.system('avogadro tmp.xyz')
|
|
os.system('rm tmp.xyz')
|
|
|
|
|
|
main()
|