Changed main.py name to crystel_mec.py

crystal_mec.py

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