I got the script working I think?!

Now we know that only space groups 1 and trigonal groups will satisfy
the requirements and are possible candidates.

first-gen-filler.py

@@ -28,10 +28,12 @@
 import pymatgen.io.vasp
 import pymatgen.symmetry.analyzer
 
+
 def chi(x):
-    out = integrate.quad(lambda t: math.exp(-((t ** 2)/2)), 0, x) 
+    out = integrate.quad(lambda t: math.exp(-((t ** 2)/2)), 0, x)
     return out[0]
 
+
 CU_OCC = 18
 CU_AV = 108
 AG_OCC = 9
@@ -56,11 +58,11 @@ def chi(x):
 arg_parser.add_argument("-n", "--number", type=int, default=100,
                         help="Specifies how many configurations there will be\
                         in the generation. Default 100.")
-arg_parser.add_argument("-a", "--aggressiveness", type=float, default=3, help=\
-                        "Specifies how aggressively biased the binning will be\
-                        towards higher order space groups. Mathematically, this\
-                        is specifying a z-score as a cutoff on the curve\
-                        that is being sampled. Default 3.")
+arg_parser.add_argument("-a", "--aggressiveness", type=float, default=3,
+                        help="Specifies how aggressively biased the binning\
+                        will be towards higher order space groups.\
+                        Mathematically, this is specifying a z-score as a\
+                        cutoff on the curve that is being sampled. Default 3.")
 
 args = arg_parser.parse_args()
 
@@ -81,7 +83,8 @@ def chi(x):
 # A: 0, B: 1, C: 2
 LatticeMatrix = [index[i] for i in range(2, 5)]
 
-Elements = [i for i in zip(index[5].split(), [int(x) for x in index[6].split()])]
+Elements = [i for i in zip(index[5].split(), [int(x) for x in index[6]
+                                              .split()])]
 # Just making sure that the file is in the correct format :)
 assert len(Elements) == 3
 assert index[7].strip() == "Direct"
@@ -119,13 +122,13 @@ def chi(x):
 # Calculate the distribution of space groups into bins on a normal
 # distribution
 if args.aggressiveness == 0:
-    Bins = [round((SPACE_GROUPS * i) / args.bins) for i in\
+    Bins = [round((SPACE_GROUPS * i) / args.bins) for i in
             range(1, args.bins + 1)]
 else:
-    Bins = [round(SPACE_GROUPS *\
-                       ((chi((args.aggressiveness * i) / args.bins)) /\
-                        (chi(args.aggressiveness)))) for i in\
-                            range(1, args.bins + 1)]
+    Bins = [round(SPACE_GROUPS *
+                  ((chi((args.aggressiveness * i) / args.bins)) /
+                   (chi(args.aggressiveness)))) for i in
+            range(1, args.bins + 1)]
 BinSize = [0] * args.bins
 MAX_BIN_SIZE = args.number // args.bins
 
@@ -167,7 +170,7 @@ def chi(x):
     # returned 1 ever. Am I just not scanning enough input or is there actually
     # a problem with the system?
     spacegroup = pymatgen.symmetry.analyzer.\
-        SpacegroupAnalyzer(poscar.structure, symprec=R_TOL,\
+        SpacegroupAnalyzer(poscar.structure, symprec=R_TOL,
                            angle_tolerance=A_TOL).get_space_group_number()
 
     vasp_file = open(file_path, "r")
@@ -192,5 +195,5 @@ def chi(x):
                 break
             os.remove(file_path)
             break
-    
+
     attempt += 1

spacegroup-sym.py

@@ -1,7 +1,7 @@
 import yaml
 import argparse
 import sympy
-import pdb
+import math
 
 x, y, z = sympy.symbols('x, y, z')
 
@@ -25,6 +25,7 @@
 input = yaml.safe_load(input)
 
 # Big ass for loop - determine symmetry operations
+passing_list = [1]
 for spacegroup in input.items():
     if spacegroup[0] == 1:
         print(f"Group {spacegroup[0]} is closed")
@@ -59,11 +60,11 @@
     # Now run each site through the same set of ops and see if we get the
     # same site back
     for op in symmetry_ops:
-        # TODO: make sure that this will result in N vectors (okay bc will run N times)
+        if op == sympy.Matrix([[x], [y], [z]]):
+            continue
         sym_op_alt = []
         compose = [(x, op[0]), (y, op[1]), (z, op[2])]
         for op2 in symmetry_ops:
-            breakpoint()
             new_vec = sympy.Matrix([op2[0].subs(compose), op2[1].subs(compose),
                                     op2[2].subs(compose)])
             sym_op_alt.append(new_vec)
@@ -80,14 +81,31 @@
                 op.simplify()
         # Now check to see if all operations are contained within the original
         # list
-        while isClosed:
-            for op2 in sym_op_alt:
-                if op2 not in symmetry_ops:
-                    isClosed = False
-                    break
+        for op2 in sym_op_alt:
+            if op2 not in symmetry_ops:
+                isClosed = False
+                break
+
+        if not isClosed:
+            break
 
-        # Print whether or not the space group is closed
-        if isClosed:
-            print(f"Group {spacegroup[0]} is closed")
+    # Print whether or not the space group is closed
+    if isClosed:
+        print(f"Group {spacegroup[0]} is closed")
+        if len(symmetry_ops) > 9:
+            print(f"Group {spacegroup[0]} fails")
         else:
-            print(f"Group {spacegroup[0]} is open")
+            if 9 % len(symmetry_ops) == 0:
+                print(f"Group {spacegroup[0]} passes")
+                passing_list.append(spacegroup[0])
+            else:
+                print(f"Group {spacegroup[0]} fails")
+    else:
+        print(f"Group {spacegroup[0]} is open")
+        if math.gcd(len(symmetry_ops), 9) != 1:
+            print(f"Group {spacegroup[0]} passes")
+            passing_list.append(spacegroup[0])
+        else:
+            print(f"Group {spacegroup[0]} fails")
+
+print(passing_list)