mirror of https://github.com/NehZio/Crystal-MEC
All the code at once
This commit is contained in:
Léo Gaspard
parent
76d3adcb0a
commit
c5a5281810
|
|
@ -0,0 +1,163 @@
|
|||
import sys
|
||||
import numpy as np
|
||||
|
||||
from src.read_input import *
|
||||
from src.out import *
|
||||
from src.utils import *
|
||||
|
||||
if __name__=='__main__':
|
||||
|
||||
verbose = 2
|
||||
|
||||
if len(sys.argv) == 1:
|
||||
print("Please provide the input file")
|
||||
sys.exit()
|
||||
elif len(sys.argv) == 3 :
|
||||
if sys.argv[2] == 's':
|
||||
verbose = 0
|
||||
elif sys.argv[2] == 'v':
|
||||
verbose = 1
|
||||
elif sys.argv[2] == 'vv':
|
||||
verbose = 2
|
||||
elif sys.argv[2] == 'vvv':
|
||||
verbose = 3
|
||||
else:
|
||||
print("Unknown argument -%s-"%sys.argv[3])
|
||||
elif len(sys.argv) > 3:
|
||||
print("Too much arguments, please provide an input file and a verbose level (v, vv, vvv)")
|
||||
sys.exit()
|
||||
|
||||
inputFile = sys.argv[1]
|
||||
|
||||
# Reads all the parameters from the input file
|
||||
rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation = read_input(inputFile)
|
||||
|
||||
if verbose > 0:
|
||||
out_input_param(rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation)
|
||||
|
||||
# Converting the angles to radian
|
||||
alpha = alpha * np.pi / 180.0
|
||||
beta = beta * np.pi / 180.0
|
||||
gamma = gamma * np.pi / 180.0
|
||||
|
||||
|
||||
# Computing the number of replications needed in each directions using the interreticular distance
|
||||
# for the planes 100, 010 and 001.
|
||||
#
|
||||
# The condition is : d_{hkl} >= 2*bath_radius + |translation_vector|
|
||||
|
||||
fac = np.sqrt(1-np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma))
|
||||
|
||||
nA = int(np.ceil(np.sin(alpha)*(2*rB+np.linalg.norm(translation))/(a*fac)))+1
|
||||
nB = int(np.ceil(np.sin(beta)*(2*rB+np.linalg.norm(translation))/(b*fac)))+1
|
||||
nC = int(np.ceil(np.sin(gamma)*(2*rB+np.linalg.norm(translation))/(c*fac)))+1
|
||||
|
||||
if verbose > 1:
|
||||
print("The big cell will be of dimensions %2ix%2ix%2i\n"%(nA,nB,nC))
|
||||
|
||||
coordinates = big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC)
|
||||
|
||||
# Computing the translation vector by addition of the user translation vector and a translation
|
||||
# vector that puts the origin at the center of the big cell
|
||||
t = [-0.5,-0.5,-0.5]
|
||||
M = get_cell_matrix(nA*a,nB*b,nC*c,alpha,beta,gamma)
|
||||
t = np.matmul(M,t)
|
||||
t = [t[0]+translation[0],t[1]+translation[1],t[2]+translation[2]]
|
||||
|
||||
# Translating the coordinates
|
||||
coordinates = translate(t, coordinates)
|
||||
|
||||
# Finding the center and translating the coordinates
|
||||
# If this vector creates a displacment bigger than a, b or c
|
||||
# in any of the abc directions, this might result in an incomplete
|
||||
# sphere later, the user should provide a translation vector
|
||||
# to correct this
|
||||
if center != []:
|
||||
c = find_center(center,coordinates)
|
||||
coordinates = translate(-c,coordinates)
|
||||
|
||||
|
||||
# Orienting the big cell
|
||||
if X != []:
|
||||
k = [1,0,0]
|
||||
|
||||
xVec = find_center(X,coordinates)
|
||||
M = rotation_matrix(k,xVec)
|
||||
|
||||
rotate(M, coordinates)
|
||||
if Y != []:
|
||||
k = [0,1,0]
|
||||
|
||||
yVec = find_center(Y,coordinates)
|
||||
M = rotation_matrix(k,yVec)
|
||||
|
||||
rotate(M, coordinates)
|
||||
if Z != []:
|
||||
k = [0,0,1]
|
||||
|
||||
zVec = find_center(Z,coordinates)
|
||||
M = rotation_matrix(k,zVec)
|
||||
|
||||
rotate(M, coordinates)
|
||||
|
||||
if verbose > 2:
|
||||
print("The big cell contains %5i atoms and will be printed in the file big_cell.xyz\n"%len(coordinates))
|
||||
write_coordinates(coordinates,'big_cell.xyz',3)
|
||||
|
||||
# Cutting the sphere in the big cell
|
||||
coordinates = cut_sphere(coordinates,rB)
|
||||
|
||||
if verbose > 2:
|
||||
print("The sphere contains %5i atoms and will be printed in the file sphere.xyz\n"%len(coordinates))
|
||||
write_coordinates(coordinates,'sphere.xyz',3)
|
||||
|
||||
# Finding the fragment
|
||||
|
||||
nAt, coordinates = find_fragment(coordinates,pattern,npattern,notInFrag)
|
||||
|
||||
if verbose > 2 or showFrag:
|
||||
print("The fragment contains %3i atoms and will be printed in the file fragment.xyz\n"%nAt)
|
||||
write_coordinates(coordinates,'fragment.xyz',4,'O')
|
||||
|
||||
coordinates = find_pseudo(coordinates,rPP,notInPseudo)
|
||||
|
||||
if verbose > 2 or showBath:
|
||||
print("The bath will be printed in the file bath.xyz\n")
|
||||
write_coordinates(coordinates,'bath.xyz',3)
|
||||
print("The bath sorted with the fragment/pseudo/charge will be printed in the file bath_coloured.xyz\n")
|
||||
write_coordinates(coordinates,'bath_coloured.xyz',3,color='yes')
|
||||
|
||||
if evjen:
|
||||
charges = evjen_charges(coordinates,atoms)
|
||||
else:
|
||||
charges = []
|
||||
for i in range(len(coordinates)):
|
||||
li = coordinates[i][3]
|
||||
ii = np.where(atoms=li)[0]
|
||||
charges.append(atoms[ii+1])
|
||||
|
||||
if verbose > 1:
|
||||
print("The total charge fragment+pseudopotential+bath is : % 8.6f\n"%np.sum(charges))
|
||||
|
||||
if symmetry != []:
|
||||
nuc1 = nuclear_repulsion(coordinates,charges)
|
||||
if verbose > 1:
|
||||
print("Nuclear repulsion before the symmetry : % 8.6f\n"%nuc1)
|
||||
|
||||
coordinates,charges,indexList = compute_symmetry(coordinates,charges,symmetry)
|
||||
|
||||
nuc2 = nuclear_repulsion(coordinates,charges)
|
||||
if verbose > 1:
|
||||
print("Nuclear repulsion after the symmetry : % 8.6f\n"%nuc2)
|
||||
print("The total charge fragment+pseudopotential+bath after symmetry is : % 8.6f\n"%np.sum(charges))
|
||||
if verbose > 2:
|
||||
print("The symmetrized coordinates contain %5i atoms \n"%len(indexList))
|
||||
|
||||
else:
|
||||
indexList = [i for i in range(len(coordinates))]
|
||||
|
||||
write_output(outputFile,coordinates,charges,indexList)
|
||||
if verbose > 2:
|
||||
print("The output has been written to %s \n"%outputFile)
|
||||
out_interatomic_distances(coordinates)
|
||||
|
||||
|
|
@ -0,0 +1,157 @@
|
|||
import numpy as np
|
||||
import operator
|
||||
|
||||
from src.utils import distance
|
||||
|
||||
def out_input_param(rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation):
|
||||
print("Output file is : %s\n"%outputFile)
|
||||
print("Bath radius : % 16.10f\n"%rB)
|
||||
print("Pseudo-potential width : % 16.10f\n"%rPP)
|
||||
print("Lattice parameters :")
|
||||
print(" a = % 8.5f"%a)
|
||||
print(" b = % 8.5f"%b)
|
||||
print(" c = % 8.5f"%c)
|
||||
print(" alpha = % 6.5f"%alpha)
|
||||
print(" beta = % 6.5f"%beta)
|
||||
print(" gamma = % 6.5f\n"%gamma)
|
||||
print("Group symmetry operations :")
|
||||
for i in symGenerator:
|
||||
print(" %10s %10s %10s"%(i[0],i[1],i[2]))
|
||||
print("\nAtomic positions :")
|
||||
for i in generator:
|
||||
print(" %2s % 8.5f % 8.5f % 8.5f"%(i[0],i[1],i[2],i[3]))
|
||||
print("\nThere %3s %2i %8s :"%('is'*(len(pattern)==1)+'are'*(len(pattern)>1),len(pattern),'pattern'+'s'*(len(pattern)>1)))
|
||||
for i in range(len(pattern)):
|
||||
print("The following pattern will appear %2i %5s"%(npattern[i],'time'+'s'*(npattern[i]>1)))
|
||||
for j in range(0,len(pattern[i]),2):
|
||||
print(" %2i %2s\n"%(pattern[i][j],pattern[i][j+1]),end='')
|
||||
print("\nThere %3s %2i %5s"%('is'*(len(atoms)==1)+'are'*(len(atoms)>1), len(atoms), 'atom'+'s'*(len(atoms)>1)))
|
||||
print("Label Charge Number of neighbours Maximum bond length")
|
||||
for i in atoms:
|
||||
print(" %2s % 5.2f %2i % 6.4f\n"%(i[0],i[1],i[2],i[3]))
|
||||
if(len(center)>0):
|
||||
print("Centering between the following atom(s) :")
|
||||
for i in center:
|
||||
print(" %2s \n"%i,end='')
|
||||
print()
|
||||
if(len(X) > 0):
|
||||
print("Aligning the X axis with the following atom(s) :")
|
||||
for i in X:
|
||||
print(" %2s "%i,end='')
|
||||
print()
|
||||
if(len(Y) > 0):
|
||||
print("Aligning the Y axis with the following atom(s) :")
|
||||
for i in Y:
|
||||
print(" %2s "%i,end='')
|
||||
print()
|
||||
if(len(Z) > 0):
|
||||
print("Aligning the Z axis with the following atom(s) :")
|
||||
for i in Z:
|
||||
print(" %2s "%i,end='')
|
||||
print()
|
||||
if translation != [0.0,0.0,0.0]:
|
||||
print("The following translation will be applied : % 5.3f % 5.3f % 5.3f"%(translation[0],translation[1],translation[2]))
|
||||
if(len(symmetry) > 0):
|
||||
print("Treating the following symmetry operations :")
|
||||
for i in symmetry:
|
||||
print(" %3s "%i,end='')
|
||||
print()
|
||||
if(evjen):
|
||||
print("The program will reequilibrate the charges at the limits of the spheres using Evjen method\n")
|
||||
if(showFrag):
|
||||
print("The program will print the fragment coordinates in the file fragment.xyz\n")
|
||||
if(showFrag):
|
||||
print("The program will print the bath coordinates in the files bath.xyz and bath_coloured.xyz\n")
|
||||
if len(notInPseudo) > 0:
|
||||
print("The following %5s will not be considered in the pseudopotential shell :"%('atom'+'s'*(len(notInPseudo)>1)))
|
||||
for i in notInPseudo:
|
||||
print(" %2s "%i)
|
||||
print()
|
||||
if len(notInFrag) > 0:
|
||||
print("The following %5s will be excluded from the fragment :"%('atom'+'s'*(len(notInFrag)>1)))
|
||||
for i in notInFrag:
|
||||
print(" %2s % 8.6f % 8.6f % 8.6f"%(i[0],i[1],i[2],i[3]))
|
||||
print()
|
||||
|
||||
# Write the coordinates to a file using the string
|
||||
# at the position 'label' as atom name
|
||||
# if which is specified, it must be an array of labels
|
||||
# color = 'no' means that the label will be printed
|
||||
# color = 'yes' means that the color in the bath will be printed
|
||||
def write_coordinates(coordinates, fileName, label, which='all',color='no'):
|
||||
|
||||
f = open(fileName,'w')
|
||||
l = [i[label] for i in coordinates]
|
||||
|
||||
if color=='no':
|
||||
lab = 3
|
||||
else:
|
||||
lab = 4
|
||||
|
||||
unique, count = np.unique(l,return_counts=True)
|
||||
labelDic = dict(zip(unique,count))
|
||||
|
||||
if which=='all':
|
||||
f.write("%i \n\n"%len(coordinates))
|
||||
for i in coordinates:
|
||||
f.write("%2s % 10.6f % 10.6f % 10.6f\n"%(i[lab],i[0],i[1],i[2]))
|
||||
else:
|
||||
count = 0
|
||||
for la in which:
|
||||
count += labelDic[la]
|
||||
f.write("%i \n\n"%count)
|
||||
for i in coordinates:
|
||||
if i[label] in which:
|
||||
f.write("%2s % 10.6f % 10.6f % 10.6f\n"%(i[lab],i[0],i[1],i[2]))
|
||||
|
||||
f.close()
|
||||
|
||||
# Writes the output file
|
||||
def write_output(outputFile, coordinates, charges, indexList):
|
||||
|
||||
f = open(outputFile,'w')
|
||||
|
||||
d = {'O':'FRAGMENT','Cl':'PSEUDOPOTENTIAL','C':'BATH'}
|
||||
|
||||
count = 0
|
||||
|
||||
for t in d:
|
||||
f.write("%s\n"%d[t])
|
||||
|
||||
f.write(" Label x y z Charge\n")
|
||||
for i in indexList:
|
||||
if coordinates[i][4] == t:
|
||||
count += 1
|
||||
l = coordinates[i][3]+str(count)
|
||||
f.write("%6s % 16.10f % 16.10f % 16.10f % 8.5f\n"%(l,coordinates[i][0],coordinates[i][1],coordinates[i][2],charges[i]))
|
||||
f.write("\n")
|
||||
|
||||
|
||||
# Writes the interatomic distances inferior to 4 A
|
||||
def out_interatomic_distances(coordinates):
|
||||
distanceList = []
|
||||
|
||||
for i in range(len(coordinates)-1):
|
||||
for j in range(i+1,len(coordinates)):
|
||||
a = coordinates[i]
|
||||
b = coordinates[j]
|
||||
d = distance(a,b)
|
||||
|
||||
accept = True
|
||||
if d > 4:
|
||||
accept = False
|
||||
else:
|
||||
for k in distanceList:
|
||||
if ((a[3] == k[0] and b[3] == k[1]) or (a[3] == k[1] and b[3] == k[0])) and np.abs(d-k[2]) < 1e-5:
|
||||
accept = False
|
||||
|
||||
if accept:
|
||||
distanceList.append([a[3],b[3],d])
|
||||
|
||||
distanceList = sorted(distanceList,key=operator.itemgetter(2))
|
||||
|
||||
print("All the interatomic distances < 5 A :")
|
||||
print("Atom1 Atom2 Distance")
|
||||
for i in distanceList:
|
||||
print(" %2s %2s %6.4f"%(i[0],i[1],i[2]))
|
||||
|
||||
|
|
@ -0,0 +1,195 @@
|
|||
import sys
|
||||
|
||||
# Parses the input file
|
||||
def read_input(inputFile):
|
||||
|
||||
rB = 0.0
|
||||
rPP = 0.0
|
||||
center = []
|
||||
X = []
|
||||
Y = []
|
||||
Z = []
|
||||
symmetry = []
|
||||
outputFile = ""
|
||||
pattern = []
|
||||
npattern = []
|
||||
atoms = []
|
||||
dist = []
|
||||
a = 0.0
|
||||
b = 0.0
|
||||
c = 0.0
|
||||
alpha = 90.0
|
||||
beta = 90.0
|
||||
gamma = 90.0
|
||||
showBath = False
|
||||
evjen = False
|
||||
showFrag = False
|
||||
notInPseudo = []
|
||||
notInFrag = []
|
||||
symGenerator = []
|
||||
generator = []
|
||||
translate = [0.0,0.0,0.0]
|
||||
|
||||
checkInput = {'bath':False,'pseudo':False,'output':False,'pattern':False,'npattern':False,'a':False,'b':False,'c':False,'atoms':False,'symmetry_generator':False,'generator':False}
|
||||
|
||||
f = open(inputFile,'r')
|
||||
|
||||
line = 'x'
|
||||
|
||||
while line.casefold() != 'end_of_input':
|
||||
line = f.readline().strip()
|
||||
ls = line.split()
|
||||
if ls == []:
|
||||
continue
|
||||
if ls[0].casefold() in checkInput:
|
||||
checkInput[ls[0].casefold()] = True
|
||||
if ls[0].casefold() == 'bath':
|
||||
try:
|
||||
rB = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : %s is not a valid bath radius value"%(ls[1]))
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'pseudo':
|
||||
try:
|
||||
rPP = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : %s is not a valid pseudopotential radius value"%(ls[1]))
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'center':
|
||||
ls.pop(0)
|
||||
center = [i for i in ls]
|
||||
elif ls[0].casefold() == 'x_axis':
|
||||
ls.pop(0)
|
||||
X = [i for i in ls]
|
||||
elif ls[0].casefold() == 'y_axis':
|
||||
ls.pop(0)
|
||||
Y = [i for i in ls]
|
||||
elif ls[0].casefold() == 'z_axis':
|
||||
ls.pop(0)
|
||||
Z = [i for i in ls]
|
||||
elif ls[0].casefold() == 'symmetry':
|
||||
ls.pop(0)
|
||||
symmetry = [i for i in ls]
|
||||
elif ls[0].casefold() == 'output':
|
||||
outputFile = ls[1]
|
||||
elif ls[0].casefold() == 'pattern':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
pattern.append([])
|
||||
ls = line.split()
|
||||
for i in range(len(ls)):
|
||||
if i%2 == 0:
|
||||
pattern[-1].append(int(ls[i]))
|
||||
else:
|
||||
pattern[-1].append(ls[i])
|
||||
line = f.readline()
|
||||
elif ls[0].casefold() == 'npattern':
|
||||
ls.pop(0)
|
||||
try:
|
||||
npattern = [int(i) for i in ls]
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : the number of patterns is not valid %s"%(line))
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'lattice':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
ls = line.split()
|
||||
if ls[0].casefold() in checkInput:
|
||||
checkInput[ls[0].casefold()] = True
|
||||
if ls[0].casefold() == 'a':
|
||||
try:
|
||||
a = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'b':
|
||||
try:
|
||||
b = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'c':
|
||||
try:
|
||||
c = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'alpha':
|
||||
try:
|
||||
alpha = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'beta':
|
||||
try:
|
||||
beta = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'gamma':
|
||||
try:
|
||||
gamma = float(ls[1])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
|
||||
sys.exit()
|
||||
line = f.readline()
|
||||
elif ls[0].casefold() == 'atoms':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
ls = line.split()
|
||||
if(len(ls)) != 4:
|
||||
print("Error while parsing the input file : not enough values given for the atom in line %s"%line)
|
||||
sys.exit()
|
||||
try:
|
||||
atoms.append([ls[0], float(ls[1]), int(ls[2]), float(ls[3])])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the atom %s"%line)
|
||||
line = f.readline()
|
||||
elif ls[0].casefold() == 'show_bath':
|
||||
showBath = True
|
||||
elif ls[0].casefold() == 'translate':
|
||||
ls.pop(0)
|
||||
try:
|
||||
translate = [float(ls[0]),float(ls[1]),float(ls[2])]
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : the translation vector is not valid %s"%line)
|
||||
sys.exit()
|
||||
elif ls[0].casefold() == 'not_in_pseudo':
|
||||
ls.pop(0)
|
||||
notInPseudo = [i for i in ls]
|
||||
elif ls[0].casefold() == 'show_frag':
|
||||
showFrag = True
|
||||
elif ls[0].casefold() == 'evjen':
|
||||
evjen = True
|
||||
elif ls[0].casefold() == 'symmetry_generator':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
symGenerator.append(line.split(','))
|
||||
line = f.readline()
|
||||
elif ls[0].casefold() == 'generator':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
ls = line.split()
|
||||
try:
|
||||
generator.append([ls[0],float(ls[1]),float(ls[2]),float(ls[3])])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the generator atom %s"%line)
|
||||
sys.exit()
|
||||
line = f.readline()
|
||||
elif ls[0].casefold() == 'not_in_frag':
|
||||
line = f.readline()
|
||||
while line.strip().casefold() != 'end':
|
||||
ls = line.split()
|
||||
try:
|
||||
notInFrag.append([ls[0],float(ls[1]),float(ls[2]),float(ls[3])])
|
||||
except ValueError:
|
||||
print("Error while parsing the input file : bad value for the atom %s"%line)
|
||||
sys.exit()
|
||||
line = f.readline()
|
||||
f.close()
|
||||
for t in checkInput:
|
||||
if checkInput[t] == False:
|
||||
print("Bad input : missing the keyword -- %s --"%t)
|
||||
sys.exit()
|
||||
|
||||
return rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate
|
||||
|
|
@ -0,0 +1,373 @@
|
|||
import numpy as np
|
||||
import operator
|
||||
import sys
|
||||
|
||||
def distance(a,b):
|
||||
# Returns the 3D distance between a and b where
|
||||
# a and b are array where x, y and z are at the
|
||||
# position 0, 1 and 2
|
||||
|
||||
x = a[0]-b[0]
|
||||
y = a[1]-b[1]
|
||||
z = a[2]-b[2]
|
||||
|
||||
return np.sqrt(x**2+y**2+z**2)
|
||||
|
||||
def get_cell_matrix(a,b,c,alpha,beta,gamma):
|
||||
# Computing the volume of the primitive cell
|
||||
omega = a*b*c*np.sqrt(1-np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma))
|
||||
|
||||
# Computing the matrix
|
||||
M = [
|
||||
[a,b*np.cos(gamma),c*np.cos(beta)],
|
||||
[0,b*np.sin(gamma),c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/(np.sin(gamma))],
|
||||
[0,0,omega/(a*b*np.sin(gamma))]
|
||||
]
|
||||
return M
|
||||
|
||||
def big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC):
|
||||
coords = []
|
||||
|
||||
# Computing the matrix converting fractional to cartesian
|
||||
fracToCart = get_cell_matrix(a,b,c,alpha,beta,gamma)
|
||||
|
||||
for gen in generator:
|
||||
x = gen[1]
|
||||
y = gen[2]
|
||||
z = gen[3]
|
||||
|
||||
for sym in symGenerator:
|
||||
u = eval(sym[0])
|
||||
v = eval(sym[1])
|
||||
w = eval(sym[2])
|
||||
|
||||
# Making sure the value is within the range [0,1]
|
||||
u = u + 1*(u<0) - 1*(u>1)
|
||||
v = v + 1*(v<0) - 1*(v>1)
|
||||
w = w + 1*(w<0) - 1*(w>1)
|
||||
|
||||
coords.append([u,v,w,gen[0]])
|
||||
|
||||
# Deleting the redundant atoms
|
||||
toDel = []
|
||||
for i in range(len(coords)-1):
|
||||
for j in range(i+1,len(coords)):
|
||||
# Computing the distance using the minimum image convention
|
||||
# as described in Appendix B equation 9 of
|
||||
# "Statistical Mechanics : Theory and Molecular Simulations
|
||||
# Mark E. Tuckerman"
|
||||
|
||||
r1 = np.array(coords[i][:3])
|
||||
r2 = np.array(coords[j][:3])
|
||||
r12 = r1-r2
|
||||
r12 = r12 - np.round(r12)
|
||||
r12 = np.matmul(fracToCart,r12)
|
||||
d = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
|
||||
|
||||
if(d<1e-5):
|
||||
# We check if we don't already want to delete this atom
|
||||
if j not in toDel:
|
||||
toDel.append(j)
|
||||
|
||||
# We delete the atoms in the list
|
||||
for i in range(len(toDel)):
|
||||
coords.pop(toDel[i]-i)
|
||||
|
||||
newCoords = []
|
||||
|
||||
# We replicate the cell nA, nB, nC times
|
||||
for at in coords:
|
||||
newCoords.append([at[0],at[1],at[2],at[3]])
|
||||
for a in range(1,nA):
|
||||
newCoords.append([at[0]+a,at[1],at[2],at[3]])
|
||||
for b in range(1,nB):
|
||||
newCoords.append([at[0]+a,at[1]+b,at[2],at[3]])
|
||||
for c in range(1,nC):
|
||||
newCoords.append([at[0]+a,at[1]+b,at[2]+c,at[3]])
|
||||
for c in range(1,nC):
|
||||
newCoords.append([at[0]+a,at[1],at[2]+c,at[3]])
|
||||
for b in range(1,nB):
|
||||
newCoords.append([at[0],at[1]+b,at[2],at[3]])
|
||||
for c in range(1,nC):
|
||||
newCoords.append([at[0],at[1]+b,at[2]+c,at[3]])
|
||||
for c in range(1,nC):
|
||||
newCoords.append([at[0],at[1],at[2]+c,at[3]])
|
||||
|
||||
# Now we convert the fractionnal coordinates to cartesian coordinates
|
||||
coords = []
|
||||
|
||||
for at in newCoords:
|
||||
r = [at[0],at[1],at[2]]
|
||||
rxyz = np.matmul(fracToCart,r)
|
||||
coords.append([rxyz[0],rxyz[1],rxyz[2],at[3],'C'])
|
||||
|
||||
# Returns the list of the atoms [x,y,z,label,second_label]
|
||||
return coords
|
||||
|
||||
# Translates all the coordinates with the vector v
|
||||
def translate(v,coordinates):
|
||||
for c in coordinates:
|
||||
c[0] += v[0]
|
||||
c[1] += v[1]
|
||||
c[2] += v[2]
|
||||
return coordinates
|
||||
|
||||
# Finds the point at the center of the given atoms that are the
|
||||
# closest to the origin
|
||||
def find_center(centerList, coordinates):
|
||||
|
||||
centers = []
|
||||
for i in range(len(centerList)):
|
||||
centers.append([100,100,100]) # Setting a large value for each center
|
||||
|
||||
for c in centers:
|
||||
c.append(distance(c,[0,0,0])) # Computing the distance to the origin
|
||||
|
||||
for at in coordinates:
|
||||
if at[3] in centerList:
|
||||
centers = sorted(centers, key=operator.itemgetter(3)) # Sorting the list with respect to the distance to the origin
|
||||
d = distance(at,[0,0,0])
|
||||
if d <= centers[-1][-1]:
|
||||
centers[-1] = [at[0],at[1],at[2],d]
|
||||
|
||||
center = np.mean(centers,axis=0)[:3] # Computing the barycenter
|
||||
|
||||
return center
|
||||
|
||||
# Defines a rotation matrix that will put r1 at the position r2
|
||||
def rotation_matrix(r1,r2):
|
||||
|
||||
r1 = np.array(r1)/np.linalg.norm(r1)
|
||||
r2 = np.array(r2)/np.linalg.norm(r2)
|
||||
|
||||
# Computing the cross product which is the vector around which
|
||||
# the rotation is done
|
||||
crossProduct = np.cross(r1,r2)
|
||||
crossProduct = crossProduct/np.linalg.norm(crossProduct)
|
||||
|
||||
# Computing the angle of the rotation
|
||||
a = np.arccos(np.dot(r1,r2))
|
||||
|
||||
c = np.cos(a)
|
||||
s = np.sin(a)
|
||||
x = crossProduct[0]
|
||||
y = crossProduct[1]
|
||||
z = crossProduct[2]
|
||||
|
||||
M = [
|
||||
[x**2*(1-c)+c,x*y*(1-c)-z*s,x*z*(1-c)+y*s],
|
||||
[x*y*(1-c)+z*s,y**2*(1-c)+c,y*z*(1-c)-x*s],
|
||||
[x*z*(1-c)-y*s,y*z*(1-c)+x*s,z**2*(1-c)+c]
|
||||
]
|
||||
|
||||
return M
|
||||
|
||||
# Rotates all the coordinates using the rotation matric M
|
||||
def rotate(M,coordinates):
|
||||
for i in range(len(coordinates)):
|
||||
r = [coordinates[i][0],coordinates[i][1],coordinates[i][2]]
|
||||
rV = np.matmul(M,r)
|
||||
coordinates[i][0] = rV[0]
|
||||
coordinates[i][1] = rV[1]
|
||||
coordinates[i][2] = rV[2]
|
||||
|
||||
return coordinates
|
||||
|
||||
# Cuts a sphere centered on the origin in the coordinates
|
||||
def cut_sphere(coordinates,r):
|
||||
sphere = []
|
||||
for i in range(len(coordinates)):
|
||||
if distance(coordinates[i],[0,0,0]) <= r:
|
||||
sphere.append(coordinates[i])
|
||||
|
||||
return sphere
|
||||
|
||||
# Finds the fragment in the coordinates
|
||||
def find_fragment(coordinates, patterns, npatterns,notInFrag):
|
||||
|
||||
inFrag = []
|
||||
|
||||
for n in range(len(patterns)):
|
||||
pattern = patterns[n]
|
||||
npattern = npatterns[n]
|
||||
for i in range(npattern):
|
||||
c = [100,100,100]
|
||||
dc = distance([0,0,0],c)
|
||||
|
||||
inPattern = []
|
||||
# Finding the closest atom of the first type in the pattern
|
||||
for at in coordinates:
|
||||
if at[3] == pattern[1]:
|
||||
d = distance([0,0,0],at)
|
||||
if d < dc :
|
||||
accept = True
|
||||
for exc in notInFrag:
|
||||
d = distance(exc,at)
|
||||
if d < 1e-5:
|
||||
accept = False
|
||||
if accept and coordinates.index(at) not in inFrag:
|
||||
c = [at[0],at[1],at[2],0.0, coordinates.index(at)]
|
||||
dc = distance([0,0,0],c)
|
||||
# Finding the rest of the pattern around the atom previously found
|
||||
atIn = []
|
||||
for j in range(0,len(pattern),2):
|
||||
d = distance(c,[100,100,100])
|
||||
# Initializing the atoms
|
||||
for k in range(pattern[j]):
|
||||
atIn.append([100,100,100,d])
|
||||
|
||||
for at in coordinates:
|
||||
if at[3] == pattern[j+1]:
|
||||
atIn = sorted(atIn,key=operator.itemgetter(3))
|
||||
d = distance(at,c)
|
||||
trial = [at[0],at[1],at[2],d,coordinates.index(at)]
|
||||
if d < atIn[-1][3] and trial not in atIn:
|
||||
atIn[-1] = trial
|
||||
for at in atIn:
|
||||
inPattern.append(at[4])
|
||||
|
||||
for at in inPattern:
|
||||
if at not in inFrag:
|
||||
inFrag.append(at)
|
||||
|
||||
for at in inFrag:
|
||||
coordinates[at][4] = 'O'
|
||||
|
||||
return len(inFrag), coordinates
|
||||
|
||||
|
||||
# Finds the pseudopotential layer around
|
||||
# the fragment
|
||||
def find_pseudo(coordinates, rPP, notInPseudo):
|
||||
|
||||
for at in coordinates:
|
||||
if at[4] != 'O':
|
||||
continue
|
||||
for i in range(len(coordinates)):
|
||||
if coordinates[i][4] != 'C':
|
||||
continue
|
||||
d = distance(at,coordinates[i])
|
||||
if d < rPP:
|
||||
coordinates[i][4] = 'Cl'
|
||||
|
||||
return coordinates
|
||||
|
||||
# Creates lists containing the neighbours of each
|
||||
# atom
|
||||
def find_neighbours(coordinates, atoms):
|
||||
neighbourList = [[] for i in range(len(coordinates))]
|
||||
|
||||
atoms = np.array(atoms).flatten()
|
||||
|
||||
for i in range(len(coordinates)-1):
|
||||
for j in range(i+1,len(coordinates)):
|
||||
li = coordinates[i][3] # Label of the atom i
|
||||
lj = coordinates[j][3] # Label of the atom j
|
||||
|
||||
ii = np.where(atoms==li)[0]
|
||||
jj = np.where(atoms==lj)[0]
|
||||
|
||||
ci = float(atoms[ii+1]) # Charge of the atom i
|
||||
cj = float(atoms[jj+1]) # Charge of the atom j
|
||||
|
||||
if ci*cj < 0: # Checking if the charges have opposite signs
|
||||
d = distance(coordinates[i],coordinates[j])
|
||||
|
||||
if d < float(atoms[ii+3]) and d < float(atoms[jj+3]):
|
||||
neighbourList[i].append(j)
|
||||
neighbourList[j].append(i)
|
||||
return neighbourList
|
||||
|
||||
# For each atom, finds if it has the correct number of neighbours,
|
||||
# if not, modify its charge
|
||||
def evjen_charges(coordinates,atoms):
|
||||
neighbourList = find_neighbours(coordinates,atoms)
|
||||
|
||||
atoms = np.array(atoms).flatten()
|
||||
|
||||
charges = []
|
||||
|
||||
for i in range(len(coordinates)):
|
||||
li = coordinates[i][3]
|
||||
ii = np.where(atoms==li)[0]
|
||||
|
||||
nr = len(neighbourList[i])
|
||||
nt = int(atoms[ii+2])
|
||||
ci = float(atoms[ii+1])
|
||||
|
||||
if nr > nt:
|
||||
print("Error : too much neighbours for atom n°%i, count %i neighbours where it should have a maximum of %i"%(i,nr,nt))
|
||||
sys.exit()
|
||||
charges.append(ci*nr/nt)
|
||||
|
||||
return charges
|
||||
|
||||
# Computes the nuclear repulsion
|
||||
def nuclear_repulsion(coordinates,charges):
|
||||
|
||||
rep = 0.0
|
||||
|
||||
for i in range(len(coordinates)-1):
|
||||
for j in range(i+1,len(coordinates)):
|
||||
rij = distance(coordinates[i],coordinates[j])
|
||||
ci = charges[i]
|
||||
cj = charges[j]
|
||||
|
||||
if(rij < 1):
|
||||
print(i,j,"\n",coordinates[i],"\n",coordinates[j],"\n",rij)
|
||||
|
||||
rep += (ci*cj)/rij
|
||||
return rep
|
||||
|
||||
# Computes the symmetry in the whole system
|
||||
def compute_symmetry(coordinates,charges,symmetry):
|
||||
symmetrizedCoordinates = []
|
||||
symmetrizedCharges = []
|
||||
uniqueIndexList = [] # The list containing the indexes of the unique atoms
|
||||
|
||||
treated = [] # Will store the index of the atoms already treated
|
||||
|
||||
symOp = []
|
||||
|
||||
# Storing all the symmetry operations
|
||||
for s in symmetry:
|
||||
if s == 'C2x':
|
||||
symOp.append(np.array([1,-1,-1]))
|
||||
elif s == 'C2y':
|
||||
symOp.append(np.array([-1,1,-1]))
|
||||
elif s == 'C2z':
|
||||
symOp.append(np.array([-1,-1,1]))
|
||||
elif s == 'xOy':
|
||||
symOp.append(np.array([1,1,-1]))
|
||||
elif s == 'xOz':
|
||||
symOp.append(np.array([1,-1,1]))
|
||||
elif s == 'yOz':
|
||||
symOp.append(np.array([-1,1,1]))
|
||||
elif s == 'i':
|
||||
symOp.append(np.array([-1,-1,-1]))
|
||||
|
||||
for i in range(len(coordinates)):
|
||||
if i in treated:
|
||||
continue
|
||||
|
||||
treated.append(i)
|
||||
at1 = np.array(coordinates[i][:3])
|
||||
symmetrizedCoordinates.append(coordinates[i])
|
||||
symmetrizedCharges.append(charges[i])
|
||||
uniqueIndexList.append(len(symmetrizedCoordinates)-1)
|
||||
|
||||
for j in range(len(coordinates)):
|
||||
if j in treated or coordinates[i][3] != coordinates[j][3]:
|
||||
continue
|
||||
|
||||
at2 = np.array(coordinates[j][:3])
|
||||
|
||||
for s in symOp:
|
||||
if distance(at1,at1*s) > 1e-4 and distance(at2,at1*s) < 1e-4: # Checking if op.at1 != at1 and that op.at2 = at1
|
||||
p = at1*s
|
||||
treated.append(j)
|
||||
symmetrizedCoordinates.append([p[0],p[1],p[2],coordinates[i][3],coordinates[i][4]])
|
||||
symmetrizedCharges.append(charges[i])
|
||||
break
|
||||
|
||||
return symmetrizedCoordinates,symmetrizedCharges,uniqueIndexList
|
||||
Loading…
Reference in New Issue