Started final pres, code refactoring

Major improvements made to code, put the elements into classes (finally!)

.python-version

@@ -0,0 +1 @@
+3.10.12

config-gen.py

@@ -0,0 +1,216 @@
+import sympy
+import math
+import random
+import argparse
+import sys
+import yaml
+
+
+# Needs to be a list of 3d-points of type sympy.Matrix (3d colun vector)
+def list_crystal_mod(list):
+    for point in list:
+        for coord in point:
+            const = sum([term for term in coord.as_ordered_terms() if
+                        term.is_constant()])
+            if const > 0:
+                const = -sympy.floor(const)
+            else:
+                const = -sympy.ceiling(const)
+            coord = coord + const
+        point.simplify()
+
+
+def point_crystal_mod(point):
+    for coord in point:
+        const = sum([term for term in coord.as_ordered_terms() if
+                    term.is_constant()])
+        if const > 0:
+            const = -sympy.floor(const)
+        else:
+            const = -sympy.ceiling(const)
+        coord = coord + const
+    point.simplify()
+
+
+# Checks if the point is within a certain tolerance of any point contained
+# within the list.
+def is_point_in_list(my_point, my_list, tolerance):
+    for list_point in my_list:
+        distance = math.sqrt(((my_point[0] - list_point[0]) ** 2) +
+                             ((my_point[1] - list_point[1]) ** 2) +
+                             ((my_point[2] - list_point[2]) ** 2))
+        if distance <= tolerance:
+            my_point = list_point
+            return True
+    return False
+
+
+# Checks if the points in list1 are within a certain tolerance of any points
+# contained within the list2.
+def is_list_in_list(list1, list2, tolerance):
+    for my_point in list1:
+        for list_point in list2:
+            distance = math.sqrt(((my_point[0] - list_point[0]) ** 2) +
+                                 ((my_point[1] - list_point[1]) ** 2) +
+                                 ((my_point[2] - list_point[2]) ** 2))
+            if distance <= tolerance:
+                my_point = list_point
+                break
+        return False
+    return True
+
+
+class Element:
+    def __init__(self, symbol, spacegroup_num, sites):
+        self.symbol = symbol
+        self.spacegroup_num = spacegroup_num
+        self.sites = sites
+        self.unused_sites = self.sites
+        self.used_sites = []
+        self.symmetry_ops = [sympy.Matrix([x, y, z])]
+        self.inv_symmetry_ops = self.symmetry_ops
+        self.spacegroup = spacegroups[self.spacegroup_num]
+        i = 1
+        while i in self.spacegroup:
+            M_i = sympy.Matrix(self.spacegroup[i]["M"])
+            T_i = sympy.Matrix([0, 0, 0])
+            if "T" in self.spacegroup[i]:
+                T_i = sympy.Matrix(self.spacegroup[i]["T"])
+            self.symmetry_ops.append((M_i * self.symmetry_ops[0]) + T_i)
+            self.inv_symmetry_ops.append((M_i.T * self.symmetry_ops[0]) - T_i)
+            i += 1
+        if "+" in self.spacegroup:
+            for vec in self.spacegroup["+"]:
+                vector = sympy.Matrix(vec)
+                for i in range(len(self.symmetry_ops)):
+                    self.symmetry_ops.append(self.symmetry_ops[i] + vector)
+        # Now we need to remove all whole number constants from each
+        # component
+        list_crystal_mod(self.symmetry_ops)
+
+    def gen_first_point_set(self, attempt_num):
+
+        # Check for attempt timeout first
+        if attempt_num >= ATTEMPT_CREATION_TIMEOUT:
+            raise Exception("Attempt timeout reached.")
+
+        origin = random.choice(self.unused_sites)
+        temp_point_list = []
+        for op in self.symmetry_ops:
+            new_point = op.subs([(x, origin[0]), (y, origin[1]),
+                                 (z, origin[2])]).simplify()
+            point_crystal_mod(new_point)
+            temp_point_list.append(new_point)
+        if not is_list_in_list(temp_point_list, self.unused_sites):
+            self.gen_first_point_set(attempt_num + 1)
+        else:
+            for point in temp_point_list:
+                point
+
+
+DIST_TOL = 0.00000001
+CU_OCC = 18
+CU_AV = 108
+AG_OCC = 9
+BI_OCC = 9
+AG_BI_AV = 27
+I_AV_OCC = 54
+x, y, z = sympy.symbols('x, y, z')
+POSSIBLE_AG_BI_SG = [i for i in range(143, 162)]
+POSSIBLE_CU_SG = [143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,
+                  155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166,
+                  167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178,
+                  179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190,
+                  195, 198]
+ATTEMPT_CREATION_TIMEOUT = 100
+
+parser = argparse.ArgumentParser(
+        prog="Config Generator",
+        description="Given user-inputted parameters, generates VASP POSCAR\
+                files for Cu2AgBiI6 in a 3x3 crystal of a given space group\
+                from 143-161."
+        )
+
+parser.add_argument("-i", "--index", help="Filepath to the index file\
+        containing all possible positions for each of the elemental\
+        sublattices", required=True)
+
+# Check if positive later
+parser.add_argument("-n", "--number", type=int, help="Number of configurations\
+        that will be generated by the script of the specified spacegroup.",
+                    required=True)
+
+parser.add_argument("-s", "--spacegroup", type=int, choices=range(143, 162),
+                    help="The spacegroup that we will be generating\
+                            configurations to fall into", required=True)
+
+parser.add_argument("-y", "--symmetries", help="Filepath to the YAML file\
+        containing the symmetry operations for each of the spacegroups that\
+        in the specified range.", required=True)
+
+args = parser.parse_args()
+
+if args.number <= 0:
+    print("Please input a positiv number of configurations to generate.")
+    sys.exit(1)
+
+# Open index file
+index_file = open(args.index, "r")
+
+# Read in INDEX.vasp into linetable and parse components
+index = index_file.readlines()
+
+index_file.close()
+
+lattice_constant = index[1]
+# 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()])]
+# Just making sure that the file is in the correct format :)
+assert len(Elements) == 3
+assert index[7].strip() == "Direct"
+
+line_offset = 8
+AgBiIndex = []
+CuIndex = []
+IIndex = []
+
+# Load up the lists that index the possible positions for each element and
+# check to make sure that each element has the correct number of spots
+for element in Elements:
+    match element[0]:
+        case "Ag" | "Bi":
+            assert element[1] == AG_BI_AV
+            for i in range(element[1]):
+                position = sympy.Matrix([float(j) for j in index[line_offset +
+                                                                 i].split()])
+                AgBiIndex.append(position)
+            print(f"Found {element[1]} positions for {element[0]}.")
+            line_offset += element[1]
+        case "Cu":
+            assert element[1] == CU_AV
+            for i in range(element[1]):
+                position = sympy.Matrix([float(j) for j in index[line_offset +
+                                                                 i].split()])
+                CuIndex.append(position)
+            print(f"Found {element[1]} positions for {element[0]}.")
+            line_offset += element[1]
+        case "I":
+            assert element[1] == I_AV_OCC
+            for i in range(element[1]):
+                position = sympy.Matrix([float(j) for j in index[line_offset +
+                                                                 i].split()])
+                IIndex.append(position)
+            print(f"Found {element[1]} positions for {element[0]}.")
+            line_offset += element[1]
+        case _:
+            print("Syntax error on lines 6 and 7.")
+
+# Open symmetries file
+spacegroups = yaml.safe_load(open(args.symmetries, "r"))
+
+# Randomize order of filling as well as what spacegroup is assigned to each
+# element.
+

manual-config-generator.py

@@ -240,50 +240,61 @@ def growSeed(origin, symops):
                     # remaining points if they are unoccupied. If they are
                     # occupied, try the remaining possible origins before
                     # declaring it failed and starting over.
-                    if multiplier > 1 and overlapping_sites > 0:
-                        subsets = list(itertools.permutations(
-                                       used_ag_sites), overlapping_sites)
-                        num_subsets = len(subsets)
-                        possible_origins = []
-                        for i in range(num_subsets):
-                            # systematically assign each point in the subset to
-                            # a point in the set of inverse operations and try
-                            # to derive a matching origin that is within the
-                            # set of possible atomic sites.
-                            inv_op_subs = list(itertools.permuatations(
-                                               inv_symmetry_ops),
-                                               overlapping_sites)
-                            # Remove so that we don't get the original origin
-                            # in the set of new possible origins.
-                            inv_op_subs.remove(i)
-                            for inv_op_sub in inv_op_subs:
-                                # Need to now iterate through each point in
-                                # each subset and use it to find its possible
-                                # origin
-                                poss_origin = []
-                                for j in range(len(subsets[i])):
-                                    poss_origin.append(inv_op_sub[j].subs([
-                                        (x, subsets[i][j][0]),
-                                        (y, subsets[i][j][1]),
-                                        (z, subsets[i][j][2])]).simplify())
-                                # Check to see if all points in poss_origin are
-                                # the same
-                                remove_whole_number_constants(poss_origin)
-                                if listIsAllIdentical(poss_origin):
-                                    possible_origins.append(poss_origin[0])
-                        # If number of possible origins is 0, then mark as
-                        # failed and break
-                        if len(possible_origins) == 0:
-                            failed = True
-                            break
-                        elif len(possible_origins >= multiplier - num_groups):
-                            new_origins = random.sample(possible_origins,
-                                                        multiplier -
-                                                        num_groups)
-                            for org in new_origins:
-
-
-
+                    if multiplier <= 1 and overlapping_sites <= 0:
+                        break
+
+                    subsets = list(itertools.permutations(
+                                   used_ag_sites), overlapping_sites)
+                    num_subsets = len(subsets)
+                    possible_origins = []
+                    for i in range(num_subsets):
+                        # systematically assign each point in the subset to
+                        # a point in the set of inverse operations and try
+                        # to derive a matching origin that is within the
+                        # set of possible atomic sites.
+                        inv_op_subs = list(itertools.permuatations(
+                                           inv_symmetry_ops),
+                                           overlapping_sites)
+                        # Remove so that we don't get the original origin
+                        # in the set of new possible origins.
+                        inv_op_subs.remove(i)
+                        for inv_op_sub in inv_op_subs:
+                            # Need to now iterate through each point in
+                            # each subset and use it to find its possible
+                            # origin
+                            poss_origin = []
+                            for j in range(len(subsets[i])):
+                                poss_origin.append(inv_op_sub[j].subs([
+                                    (x, subsets[i][j][0]),
+                                    (y, subsets[i][j][1]),
+                                    (z, subsets[i][j][2])]).simplify())
+                            # Check to see if all points in poss_origin are
+                            # the same
+                            remove_whole_number_constants(poss_origin)
+                            if listIsAllIdentical(poss_origin):
+                                possible_origins.append(poss_origin[0])
+                    # If number of possible origins is 0, then mark as
+                    # failed and break
+                    if len(possible_origins) == 0:
+                        failed = True
+                        break
+                    elif len(possible_origins >= multiplier - num_groups):
+                        new_origins = random.sample(possible_origins,
+                                                    multiplier -
+                                                    num_groups)
+                        for org in new_origins:
+                            new_point_set = [op.subs([(x, origin[0]),
+                                                      (y, origin[1]),
+                                                      (z, origin[2])])
+                                             .simplify()
+                                             for op in symmetry_ops]
+                            remove_whole_number_constants(new_point_set)
+                            for loc in new_point_set:
+                                if loc not in unused_ag_sites and\
+                                        loc not in used_ag_sites:
+                                    failed = True
+                                    break
+                                elif loc in 
 
 
                 if failed:

notes/spacegroup-sym.md

@@ -78,3 +78,40 @@ units.
 If it is, then multiply and make sure they do not overlap.
 
 Copper space group is 221, iodine is 225. Ag/Bi is 166.
+
+## Algorithm for inverse
+
+Read in index file, spacegroup yml
+
+Read in all of the elements into sympy format
+
+Decide filling order, randomize the spacegroup
+
+Perform the filling operation for each spacegroup:
+1. get the spacegroup from the YML file into sympy form
+1. Figure out how many sites are needed left to fill, how many need to overlap. 
+
+```
+overlapping_sites = multiplier * len(symmetry_ops) - AG_OCC
+
+```
+1. Fill the first symmetry set
+    - We need to be careful about checking the tolerances on each of the
+        coordinates when doing this. Too strict and nothing will ever register,
+        too loose and we'll get error accumulating.
+1. Get all of the possible subsets of that set and subsets of the symmetry operations
+with `itertools.permuations`
+1. Permute the subsets together such that all possible pairs are generated
+    except for the ones that coincide indices (so we don't get the same set of
+    points back out.) For each of these pairs:
+    1. Call a function that calculates a potential origin for these points
+    On success, return point. On failure, return code.
+    1. If no successes, fail the attempt and skip to the next. If there is
+        a success, generate all of the points corresponding to that origin
+        and ensure that they are all valid.
+    1. Continue until we have the `multiplier` number of sites filled and
+        number of sites filled is what it's supposed to be.
+    1. Copy and save list as the final for element going into poscar, then if
+        Ag/Bi, copy use sites over to other element.
+1. Write the POSCAR file. Repeat for as many configurations as requested.
+