Performance + bugfix with not_in_frag keyword

crystal_mec.py

@@ -84,21 +84,21 @@ if __name__=='__main__':
         xVec = find_center(X,coordinates)
         M = rotation_matrix(k,xVec)
 
-        rotate(M, coordinates)
+        coordinates = rotate(M, coordinates)
     if Y != []:
         k = [0,1,0]
 
         yVec = find_center(Y,coordinates)
         M = rotation_matrix(k,yVec)
 
-        rotate(M, coordinates)
+        coordinates = rotate(M, coordinates)
     if Z != []:
         k = [0,0,1]
 
         zVec = find_center(Z,coordinates)
         M = rotation_matrix(k,zVec)
 
-        rotate(M, coordinates)
+        coordinates = 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))
@@ -113,6 +113,8 @@ if __name__=='__main__':
 
     # Finding the fragment
 
+    coordinates = sorted(coordinates, key=lambda x:distance(x,[0,0,0]))
+
     nAt, coordinates = find_fragment(coordinates,pattern,npattern,notInFrag)
 
     if verbose > 2 or showFrag:

src/out.py

@@ -70,7 +70,7 @@ def out_input_param(rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, np
         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(" % 8.6f   % 8.6f    % 8.6f"%(i[0],i[1],i[2]))           
             print()
 
 # Write the coordinates to a file using the string

src/read_input.py

@@ -181,7 +181,7 @@ def read_input(inputFile):
             while line.strip().casefold() != 'end':
                 ls = line.split()
                 try:
-                    notInFrag.append([ls[0],float(ls[1]),float(ls[2]),float(ls[3])])
+                    notInFrag.append([float(ls[0]),float(ls[1]),float(ls[2])])
                 except ValueError:
                     print("Error while parsing the input file : bad value for the atom %s"%line)
                     sys.exit()

src/utils.py

@@ -204,6 +204,8 @@ def find_fragment(coordinates, patterns, npatterns,notInFrag):
             for at in coordinates:
                 if at[3] == pattern[1]:
                     d = distance([0,0,0],at)
+                    if d > 10:
+                        break
                     if d < dc :
                         accept = True
                         for exc in notInFrag:
@@ -222,12 +224,20 @@ def find_fragment(coordinates, patterns, npatterns,notInFrag):
                     atIn.append([100,100,100,d])
 
                 for at in coordinates:
+                    if distance(at,[0,0,0]) > 10:
+                        break
                     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
+                            accept = True
+                            for exc in notInFrag:
+                                d = distance(exc,trial)
+                                if d < 1e-5:
+                                    accept = False
+                            if accept:
+                                atIn[-1] = trial
             for at in atIn:
                 inPattern.append(at[4])