diff --git a/docs/make.jl b/docs/make.jl
index 40228b6ea..1308ad6bf 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -39,7 +39,7 @@ using PorousMaterials
         repo = "github.com/SimonEnsemble/PorousMaterials.jl.git",
         # This is a link to the main repo and the master branch
         # target = "build",
-        julia = "1.0",
+        julia = "1.2",
         osname = "linux",
         deps = Deps.pip("mkdocs", "mkdocs-windmill", "pymdown-extensions") # These are dependencies for the site, not the package
     )
diff --git a/docs/src/guides/faq.md b/docs/src/guides/faq.md
index e828abb3e..01cdf6fae 100755
--- a/docs/src/guides/faq.md
+++ b/docs/src/guides/faq.md
@@ -15,8 +15,8 @@ Yes! See [here](https://github.com/JuliaLang/IJulia.jl).
 
 **How can I convert my `.cif` into P1 symmetry for `PorousMaterials.jl`?**
 
-We hope someone will contribute this feature to `PorousMaterials.jl` eventually. For now, you can use [OpenBabel](http://openbabel.org/wiki/Main_Page):
+`PorousMaterials.jl` will automatically do this for you! It looks for the
+`_symmetry_equiv_pos_as_xyz` tag in the `.cif` file and uses those symmetry operations to replicate the structure in a lower symmetry into P1 symmetry.
 
-```
-obabel -icif non-P1.cif -ocif -O P1.cif --fillUC strict
-```
+It is important to note that `PorousMaterials.jl` will read in the space group
+name, but it does **NOT** use this for converting your structure to P1.
diff --git a/docs/src/guides/help_wanted.md b/docs/src/guides/help_wanted.md
index bc984b09a..e908a56d5 100755
--- a/docs/src/guides/help_wanted.md
+++ b/docs/src/guides/help_wanted.md
@@ -8,7 +8,7 @@
 * consolidate `eikar`, `eikbr`, `eikcr` somehow without slowing down the Ewald sum.
 * more tests added to `tests/runtests.jl`, `tests/henry_tests.jl`, `tests/gcmc_tests.jl`
 * geometric-based pore size calculations (largest free and included spheres), surface area, and porosity calculations that take `Framework`'s as input
-* handle .cif's without P1 symmetry. i.e. convert any .cif to P1 symmetry
+* ~~handle .cif's without P1 symmetry. i.e. convert any .cif to P1 symmetry~~
 * extend `gcmc_simulation` to handle mixtures
 * better default rules for choosing Ewald sum parameters? alpha, kvectors required...
 * Henry coefficient code prints off Ewald sum params 5 times if run with one core...
diff --git a/docs/src/manual/boxes_crystals_grids.md b/docs/src/manual/boxes_crystals_grids.md
index 289e38dd6..509dafae0 100644
--- a/docs/src/manual/boxes_crystals_grids.md
+++ b/docs/src/manual/boxes_crystals_grids.md
@@ -1,6 +1,6 @@
 ## Loading in Crystal Structure Files
 
-Place `.cif` and `.cssr` crystal structure files in `PorousMaterials.PATH_TO_CRYSTALS`. `PorousMaterials.jl` currently takes crystals in P1 symmetry only. From here you can start julia and do the following to load a framework and start working with it.
+Place `.cif` and `.cssr` crystal structure files in `PorousMaterials.PATH_To_CRYSTALS`. `PorousMaterials.jl` can load in `.cif` files of any symmetry as long as the symmetry operations are included. From here you can start julia and do the following to load a framework and start working with it.
 
 ```julia
 using PorousMaterials
@@ -22,6 +22,91 @@ Number of charges = 0
 Chemical formula: Dict(:H=>8,:S=>1,:Ca=>1,:O=>6,:C=>14)
 ```
 
+If the file is not in P1 symmetry, it will be converted within the framework reader and this message will be displayed.
+
+```julia
+┌ Warning: Name_of_file.cif is not in P1 symmetry. It is being converted to P1 for use in PorousMaterials.jl.
+└ @ PorousMaterials ~/osu_undergrad/simon_ensemble/PorousMaterials.jl/src/Crystal.jl:284
+```
+
+PorousMaterials also gives the option to read in structures in the lower level
+symmetries and convert them to P1 before simulation.
+
+```julia
+framework = Framework("ORIVOC_clean.cif"; remove_overlap=true, convert_to_p1=false)
+# the remove_overlap argument is specific to this structure, not all frameworks need it
+
+###
+    Perform any operation on the structure while it is not in P1
+###
+
+framework = apply_symmetry_rules(framework)
+# the framework is now in P1 and it can be used in simulations
+```
+
+## Generating Bond Information for Frameworks
+
+The bonds are stored in a `SimpleGraph` from the `LightGraphs.jl` package, and
+can be accessed through the `bonds` attribute.  
+
+### Reading from a file
+
+`PorousMaterials` can read in bonds from `.cif` files if they have the tags
+`_geom_bond_atom_site_label_1` and `_geom_bond_atom_site_label_2`. To choose to
+read bonds from a file, pass `read_bonds_from_file=true` to the `Framework`
+constructor.
+
+```julia
+using PorousMaterials
+
+f = Framework("KAXQIL_clean.cif"; read_bonds_from_file=true, convert_to_p1=false)
+
+f.bonds
+```
+
+This example uses a structure that is not in P1 symmetry. `PorousMaterials`
+cannot replicate a structure or apply symmetry rules if it currently has bonds.
+However, this structure can be converted to P1 without bonds, and then bonds can
+be inferred for the full P1 structure.
+
+### Inferring bonds using `BondingRule`s
+
+`PorousMaterials` can infer bonds for a structure and populate the bond graph by
+using `BondingRule`s. Each `BondingRule` has two species of atoms that it works
+for. It also has a minimum and maximum distance that a bond can be defined for
+the two atoms.
+
+```julia
+using PorousMaterials
+
+f = Framework("SBMOF-1.cif")
+
+# define an array of BondingRule's that will be used to define bonds in the
+#   framework. These need to be in the order that they are applied
+bonding_rules = [BondingRule(:H, :*, 0.4, 1.2),
+                 BondingRule(:*, :*, 0.4, 1.9)]
+
+# Alternatively, you could get the above bonding rules with the following command
+bonding_rules = default_bondingrules()
+
+# infer the bonds for the framework f with bonds across periodic boundaries
+infer_bonds!(f, true, bonding_rules)
+
+# redefine bonding_rules to account for edge cases between Ca and O atoms. `pushfirst!` adds the newly
+#   defined Bondingrule to the front of `bonding_rules`
+pushfirst!(BondingRule(:Ca, :O, 0.4, 2.5), bonding_rules)
+
+# remove old bonds from framework before inferring bonds with new rules
+remove_bonds!(f)
+
+# re-infer bonds
+infer_bonds!(f, true, bonding_rules)
+
+# output the bond information to visualize it and double check
+write_bond_information(f, "SBMOF-1_bonds.vtk")
+```
+
+
 ## Building Blocks of PorousMaterials: Bravais lattice
 
 We later apply periodic boundary conditions to mimic a crystal of infinite extent. A `Box` describes a [Bravais lattice](https://en.wikipedia.org/wiki/Bravais_lattice).
@@ -111,13 +196,22 @@ write_cube(grid, "CH4_in_SBMOF1.cube")
     Framework
     remove_overlapping_atoms_and_charges
     strip_numbers_from_atom_labels!
+    wrap_atoms_to_unit_cell
     chemical_formula
     molecular_weight
     crystal_density
     replicate(::Framework, ::Tuple{Int, Int, Int})
     charged(::Framework; ::Bool)
-    write_cif
     assign_charges
+    apply_symmetry_rules
+    is_symmetry_equal
+    write_cif
+    BondingRule
+    write_bond_information
+    infer_bonds!
+    remove_bonds!
+    compare_bonds_in_framework
+    default_bondingrules
 ```
 
 ## Grids
diff --git a/src/Box.jl b/src/Box.jl
index 352a4ebac..88f23863d 100644
--- a/src/Box.jl
+++ b/src/Box.jl
@@ -1,6 +1,7 @@
 """
     box = Box(a, b, c, α, β, γ, volume, f_to_c, c_to_f, reciprocal_lattice)
     box = Box(a, b, c, α, β, γ)
+    box = Box(a, b, c) # α=β=γ=π/2 assumed.
     box = Box(f_to_c)
 
 Data structure to describe a unit cell box (Bravais lattice) and convert between
@@ -114,6 +115,7 @@ This function generates a unit cube, each side is 1.0 Angstrom long, and all the
 corners are right angles.
 """
 UnitCube() = Box(1.0, 1.0, 1.0, π/2, π/2, π/2)
+Box(a::Float64, b::Float64, c::Float64) = Box(a, b, c, π/2, π/2, π/2) # right angle box
 
 """
     new_box = replicate(original_box, repfactors)
diff --git a/src/Crystal.jl b/src/Crystal.jl
index 918516948..2d961421e 100644
--- a/src/Crystal.jl
+++ b/src/Crystal.jl
@@ -6,23 +6,82 @@ struct Framework
     box::Box
     atoms::Atoms
     charges::Charges
+    bonds::SimpleGraph
+    symmetry::Array{AbstractString, 2}
+    space_group::AbstractString
+    is_p1::Bool
 end
 
+function Framework(name::AbstractString, box::Box, atoms::Atoms, charges::Charges;
+                   bonds::SimpleGraph=SimpleGraph(atoms.n_atoms),
+                   symmetry::Array{AbstractString, 2}=[Array{AbstractString, 2}(undef, 3, 0) ["x", "y", "z"]],
+                   space_group::AbstractString="P1", is_p1::Bool=true)
+    return Framework(name, box, atoms, charges, bonds, symmetry, space_group, is_p1)
+end
+
+# struct for holding bonding information
+"""
+    bonding_rule = BondingRule(:Ca, :O, 0.4, 2.0)
+    bonding_rules = [BondingRule(:H, :*, 0.4, 1.2),
+                     BondingRule(:*, :*, 0.4, 1.9)]
+
+A rule for determining if two atoms within a framework are bonded. 
+
+# Attributes
+-`species_i::Symbol`: One of the atoms types for this bond rule
+-`species_j::Symbol`: The other atom type for this bond rule
+-`min_dist`: The minimum distance between the atoms for bonding to occur
+-`max_dist`: The maximum distance between the atoms for bonding to occur
+"""
+struct BondingRule
+    species_i::Symbol
+    species_j::Symbol
+    min_dist::Float64
+    max_dist::Float64
+end
+
+"""
+    default_bondingrules = default_bondingrules()
+
+Returns the default bonding rules. Using `append!` and/or `prepend!` to add to the default bonding rules:
+
+# Example
+```
+bond_rules = default_bondingrules()
+prepend!(bond_rules, BondingRule(:Cu, :*, 0.1, 2.6))
+```
+
+# Returns
+-`default_bondingrules::Array{BondingRule, 1}`: The default bonding rules: `[BondingRule(:*, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]`
+"""
+default_bondingrules() = [BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]
+
 """
     framework = Framework(filename, check_charge_neutrality=true,
                           net_charge_tol=0.001, check_atom_and_charge_overlap=true,
-                          remove_overlap=false)
-    framework = Framework(name, box, atoms, charges)
+                          remove_overlap=false, convert_to_p1=true,
+                          read_bonds_from_file=false, wrap_to_unit_cell=true)
+    framework = Framework(name, box, atoms, charges; bonds=SimpleGraph(atoms.n_atoms),
+                          symmetry=["x", "y", "z"], space_group="P1", is_p1=true)
 
 Read a crystal structure file (.cif or .cssr) and populate a `Framework` data structure,
 or construct a `Framework` data structure directly.
 
+If the framework is constructed using the `Framework(name, box, atoms, charges)`
+function it is assumed it is in P1 symmetry.
+
 # Arguments
 - `filename::AbstractString`: the name of the crystal structure file (include ".cif" or ".cssr") read from `PATH_TO_CRYSTALS`.
 - `check_charge_neutrality::Bool`: check for charge neutrality
 - `net_charge_tol::Float64`: when checking for charge neutrality, throw an error if the absolute value of the net charge is larger than this value.
 - `check_atom_and_charge_overlap::Bool`: throw an error if overlapping atoms are detected.
 - `remove_overlap::Bool`: remove identical atoms automatically. Identical atoms are the same element atoms which overlap.
+- `convert_to_p1::Bool`: If the structure is not in P1 it will be converted to
+    P1 symmetry using the symmetry rules from the `_symmetry_equiv_pos_as_xyz` list in the .cif file. 
+    (We do not use the space groups name to look up symmetry rules).
+- `read_bonds_from_file::Bool`: Whether or not to read bonding information from
+    cif file. If false, the bonds can be inferred later. note that, if the crystal is not in P1 symmetry, we cannot *both* read bonds and convert to P1 symmetry.
+- `wrap_to_unit_cell::Bool`: if true, enforce that fractional coords of atoms/charges are in [0,1]³ by mod(x, 1)
 
 # Returns
 - `framework::Framework`: A framework containing the crystal structure information
@@ -32,10 +91,23 @@ or construct a `Framework` data structure directly.
 - `box::Box`: unit cell (Bravais Lattice)
 - `atoms::Atoms`: list of Atoms in crystal unit cell
 - `charges::Charges`: list of point charges in crystal unit cell
+- `bonds::SimpleGraph`: Unweighted, undirected graph showing all of the atoms
+    that are bonded within the framework
+- `symmetry::Array{Function, 2}`: 2D array of anonymous functions that represent
+    the symmetry operations. If the structure is in P1 there will be one
+    symmetry operation.
+- `space_group::AbstractString`: The name of the space group. This is stored
+    so that it can be written out again in the write_cif function. The space
+    group is not used to verify the symmetry rules.
+- `is_p1::Bool`: Stores whether the framework is currently in P1 symmetry. Prior
+    to GCMC and Widom insertion simulations this is checked.
+- `wrap_to_unit_cell::Bool`: Whether the atom and charge positions will be
+    wrapped to the unit cell so their coordinates are in [0, 1]
 """
 function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
                    net_charge_tol::Float64=0.001, check_atom_and_charge_overlap::Bool=true,
-                   remove_overlap::Bool=false)
+                   remove_overlap::Bool=false, convert_to_p1::Bool=true,
+                   read_bonds_from_file::Bool=false, wrap_to_unit_cell::Bool=true)
     # Read file extension. Ensure we can read the file type
     extension = split(filename, ".")[end]
     if ! (extension in ["cif", "cssr"])
@@ -54,34 +126,211 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
     yf = Array{Float64, 1}()
     zf = Array{Float64, 1}()
     coords = Array{Float64, 2}(undef, 3, 0)
+    # default for symmetry rules is P1.
+    # These will be overwritten if the user chooses to read in non-P1
+    symmetry_rules = Array{AbstractString, 2}(undef, 3, 0)
+    # creating empty SimpleGraph, might not have any information read in
+    bonds = SimpleGraph()
+    # used for remembering whether fractional/cartesian coordinates are read in
+    # placed here so it will be defined for the if-stmt after the box is defined
+    fractional = false
+    cartesian = false
+    # used for determining if the framework is in P1 symmetry for simulations
+    p1_symmetry = false
+    space_group = ""
 
 
     # Start of .cif reader
+    ###################################
+    # CIF READER
+    ###################################
     if extension == "cif"
         data = Dict{AbstractString, Float64}()
         loop_starts = -1
-        for (i, line) in enumerate(lines)
-            line = split(line)
+        i = 1
+        # used for reading in symmetry options and replications
+        symmetry_info = false
+        atom_info = false
+        label_num_to_idx = Dict{AbstractString, Int}()
+        fractional = false
+        cartesian = false
+        while i <= length(lines)
+            line = split(lines[i])
             # Skip empty lines
             if length(line) == 0
+                i += 1
                 continue
             end
 
             # Make sure the space group is P1
             if line[1] == "_symmetry_space_group_name_H-M"
-                if length(line) == 3
-                    @assert (occursin("P1", line[2] * line[3]) ||
-                            occursin("P 1", line[2] * line[3]) ||
-                            occursin("-P1", line[2] * line[3]))
-                            "cif must have P1 symmetry.\n"
-                elseif length(line) == 2
-                    @assert (occursin("P1", line[2]) ||
-                            occursin("P 1", line[2]) ||
-                            occursin("-P1", line[2]))
-                            "cif must have P1 symmetry.\n"
-                else
-                    println(line)
-                    error("Does this .cif have P1 symmetry? Use `convert_cif_to_P1_symmetry` to convert to P1 symmetry")
+                # use anonymous function to combine all terms past the first
+                #   to extract space group name
+                space_group = reduce((x, y) -> x * " " * y, line[2:end])
+                space_group = split(space_group, [''', '"'], keepempty=false)[1]
+                if space_group == "P1" || space_group == "P 1" ||
+                        space_group == "-P1"
+                    # simplify by only having one P1 space_group name
+                    space_group = "P1"
+                    p1_symmetry = true
+                end
+            end
+
+            # checking for information about atom sites and symmetry
+            if line[1] == "loop_"
+                # creating dictionary of column names to determine what should be done
+                atom_column_name = ""
+                # name_to_column is a dictionary that e.g. returns which column contains x fractional coord
+                #   use example: name_to_column["_atom_site_fract_x"] gives 3
+                name_to_column = Dict{AbstractString, Int}()
+
+                i += 1
+                loop_starts = i
+                while length(split(lines[i])) == 1 && split(lines[i])[1][1] == '_'
+                    if i == loop_starts
+                        atom_column_name = split(lines[i])[1]
+                    end
+                    name_to_column[split(lines[i])[1]] = i + 1 - loop_starts
+                    # iterate to next line in file
+                    i += 1
+                end
+
+                fractional = fractional || haskey(name_to_column, "_atom_site_fract_x") &&
+                                haskey(name_to_column, "_atom_site_fract_y") &&
+                                haskey(name_to_column, "_atom_site_fract_z")
+                # if the file provides cartesian coordinates
+                cartesian = cartesian || ! fractional && haskey(name_to_column, "_atom_site_Cartn_x") &&
+                                haskey(name_to_column, "_atom_site_Cartn_y") &&
+                                haskey(name_to_column, "_atom_site_Cartn_z")
+                                             # if both are provided, will default
+                                             #  to using fractional, so keep cartesian
+                                             #  false
+
+                # =====================
+                # SYMMETRY READER
+                # =====================
+                if haskey(name_to_column, "_symmetry_equiv_pos_as_xyz")
+                    symmetry_info = true
+
+                    symmetry_count = 0
+                    # CSD stores symmetry as one column in a string that ends
+                    #   up getting split on the spaces between commas (i.e. its
+                    #   not really one column) the length(name_to_column) + 2
+                    #   should catch this hopefully there aren't other weird
+                    #   ways of writing cifs...
+                    while i <= length(lines) && length(lines[i]) > 0 && lines[i][1] != '_' && !occursin("loop_", lines[i])
+                        symmetry_count += 1
+                        line = lines[i]
+                        sym_funcs = split(line, [' ', ',', ''', '"'], keepempty=false)
+
+                        # store as strings so it can be written out later
+                        new_sym_rule = Array{AbstractString, 1}(undef, 3)
+
+                        sym_start = name_to_column["_symmetry_equiv_pos_as_xyz"] - 1
+                        for j = 1:3
+                            new_sym_rule[j] = sym_funcs[j + sym_start]
+                        end
+
+                        symmetry_rules = [symmetry_rules new_sym_rule]
+
+                        i += 1
+                    end
+
+                    @assert symmetry_count == size(symmetry_rules, 2) "number of symmetry rules must match the count"
+
+                    # finish reading in symmetry information, skip to next
+                    #   iteration of outer while-loop
+                    continue
+                # =====================
+                # FRACTIONAL READER
+                # =====================
+                elseif fractional && ! atom_info
+                    atom_info = true
+                    atom_column_name = [name for (name, column) in name_to_column if column == 1][end]
+                    
+                    while i <= length(lines) && length(split(lines[i])) == length(name_to_column)
+                        line = split(lines[i])
+
+                        push!(species, Symbol(line[name_to_column[atom_column_name]]))
+                        coords = [coords [mod(parse(Float64, split(line[name_to_column["_atom_site_fract_x"]], '(')[1]), 1.0),
+                                mod(parse(Float64, split(line[name_to_column["_atom_site_fract_y"]], '(')[1]), 1.0),
+                                mod(parse(Float64, split(line[name_to_column["_atom_site_fract_z"]], '(')[1]), 1.0)]]
+                        # If charges present, import them
+                        if haskey(name_to_column, "_atom_site_charge")
+                            push!(charge_values, parse(Float64, line[name_to_column["_atom_site_charge"]]))
+                        else
+                            push!(charge_values, 0.0)
+                        end
+                        # add to label_num_to_idx so that bonds can be converted later
+                        if read_bonds_from_file
+                            label_num_to_idx[line[name_to_column["_atom_site_label"]]] = length(species)
+                        end
+                        # iterate to next line in file
+                        i += 1
+                    end
+                    # set up graph of correct size
+                    bonds = SimpleGraph(length(species))
+                    # finish reading in atom_site information, skip to next
+                    #   iteration of outer while-loop
+                    # prevents skipping a line after finishing reading atoms
+                    continue
+                # =====================
+                # CARTESIAN READER
+                # =====================
+                elseif cartesian && ! atom_info
+                    atom_info = true
+                    atom_column_name = ""
+                    for (name, column) in name_to_column
+                        if column == 1
+                            atom_column_name = name
+                        end
+                    end
+
+                    while i <= length(lines) && length(split(lines[i])) == length(name_to_column)
+                        line = split(lines[i])
+
+                        push!(species, Symbol(line[name_to_column[atom_column_name]]))
+                        coords = [coords [parse(Float64, split(line[name_to_column["_atom_site_Cartn_x"]], '(')[1]),
+                                parse(Float64, split(line[name_to_column["_atom_site_Cartn_y"]], '(')[1]),
+                                parse(Float64, split(line[name_to_column["_atom_site_Cartn_z"]], '(')[1])]]
+                        # If charges present, import them
+                        if haskey(name_to_column, "_atom_site_charge")
+                            push!(charge_values, parse(Float64, line[name_to_column["_atom_site_charge"]]))
+                        else
+                            push!(charge_values, 0.0)
+                        end
+                        # add to label_num_to_idx so that bonds can be converted later
+                        if read_bonds_from_file
+                            label_num_to_idx[line[name_to_column["_atom_site_label"]]] = length(species)
+                        end
+                        # iterate to next line in file
+                        i += 1
+                    end
+                    # set up graph of correct size
+                    bonds = SimpleGraph(length(species))
+                    # finish reading in atom_site information, skip to next
+                    #   iteration of outer while-loop
+                    # prevents skipping a line after finishing reading atoms
+                    continue
+                # =====================
+                # BOND READER
+                # =====================
+                elseif read_bonds_from_file &&
+                       haskey(name_to_column, "_geom_bond_atom_site_label_1") &&
+                       haskey(name_to_column, "_geom_bond_atom_site_label_2")
+                    while i <= length(lines) && length(split(lines[i])) == length(name_to_column)
+                        line = split(lines[i])
+
+                        atom_one_idx = label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_1"]]]
+                        atom_two_idx = label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_2"]]]
+                        add_edge!(bonds, atom_one_idx, atom_two_idx) 
+
+                        # iterate to next line in file
+                        i += 1
+                    end
+
+                    # skip to next iteration in outer while loop
+                    continue
                 end
             end
 
@@ -99,51 +348,17 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
                 end
             end
 
-            # As soon as we reach the coordinate loop, we break this for-loop
-            # and replace it with a while-loop further down
-            if line[1] == "loop_"
-                next_line = split(lines[i+1])
-                if occursin("_atom_site", next_line[1])
-                    loop_starts = i + 1
-                    break
-                end
-            end
+            i += 1
         end # End loop over lines
 
-        if loop_starts == -1
+        if !atom_info
             error("Could not find _atom_site* after loop_ in .cif file\n")
         end
 
-        atom_column_name = ""
-        # name_to_column is a dictionary that e.g. returns which column contains x fractional coord
-        #   use example: name_to_column["_atom_site_fract_x"] gives 3
-        name_to_column = Dict{AbstractString, Int}()
-
-        i = loop_starts
-        while length(split(lines[i])) == 1
-            if i == loop_starts
-                atom_column_name = split(lines[i])[1]
-            end
-            name_to_column[split(lines[i])[1]] = i + 1 - loop_starts
-            i += 1
-        end
-
-        for i = loop_starts+length(name_to_column):length(lines)
-            line = split(lines[i])
-            if length(line) != length(name_to_column)
-                break
-            end
-
-            push!(species, Symbol(line[name_to_column[atom_column_name]]))
-            coords = [coords [mod(parse(Float64, line[name_to_column["_atom_site_fract_x"]]), 1.0),
-                    mod(parse(Float64, line[name_to_column["_atom_site_fract_y"]]), 1.0),
-                    mod(parse(Float64, line[name_to_column["_atom_site_fract_z"]]), 1.0)]]
-            # If charges present, import them
-            if haskey(name_to_column, "_atom_site_charge")
-                push!(charge_values, parse(Float64, line[name_to_column["_atom_site_charge"]]))
-            else
-                push!(charge_values, 0.0)
-            end
+        # Structure must either be in P1 symmetry or have replication information
+        if !p1_symmetry && !symmetry_info
+            error(@sprintf("%s is not in P1 symmetry and the .cif does not have a _symmetry_equiv_pos_as_xyz column 
+            for us to apply symmetry operations to convert into P1 symmetry.", filename))
         end
 
         a = data["a"]
@@ -153,7 +368,15 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
         β = data["beta"]
         γ = data["gamma"]
 
+        # redo coordinates if they were read in cartesian
+        if cartesian && ! fractional
+            coords = Box(a, b, c, α, β, γ).c_to_f * coords
+        end
+
     # Start of cssr reader #TODO make sure this works for different .cssr files!
+    ###################################
+    # CSSR READER
+    ###################################
     elseif extension == "cssr"
         # First line contains unit cell lenghts
         line = split(lines[1])
@@ -184,6 +407,11 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
         for i = 1:n_atoms
             coords = [ coords [xf[i], yf[i], zf[i]] ]
         end
+
+        # add P1 symmetry rules for consistency
+        symmetry_rules = [symmetry_rules ["x", "y", "z"]]
+        p1_symmetry = true
+        space_group = "P1"
     end
 
     # Construct the unit cell box
@@ -194,7 +422,8 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
     idx_nz = charge_values .!= 0.0
     charges = Charges(charge_values[idx_nz], coords[:, idx_nz])
 
-    framework = Framework(filename, box, atoms, charges)
+    framework = Framework(filename, box, atoms, charges; bonds=bonds, symmetry=symmetry_rules, space_group=space_group, is_p1=p1_symmetry)
+
 
     if check_charge_neutrality
         if ! charge_neutral(framework, net_charge_tol)
@@ -205,6 +434,19 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
         end
     end
 
+    if convert_to_p1 && ! p1_symmetry && ! read_bonds_from_file
+        @warn @sprintf("Framework %s has %s space group. We are converting it to P1 symmetry for use in molecular simulations.
+        To afrain from this, pass `convert_to_p1=false` to the `Framework` constructor.\n",
+            framework.name, framework.space_group)
+        return apply_symmetry_rules(framework; remove_overlap=remove_overlap,
+                                         check_charge_neutrality=check_charge_neutrality,
+                                         check_atom_and_charge_overlap=check_atom_and_charge_overlap)
+    end
+    
+    if wrap_to_unit_cell
+        wrap_atoms_to_unit_cell!(framework)
+    end
+
     if remove_overlap
         return remove_overlapping_atoms_and_charges(framework)
     end
@@ -219,6 +461,17 @@ function Framework(filename::AbstractString; check_charge_neutrality::Bool=true,
     return framework
 end
 
+"""
+    wrap_atoms_to_unit_cell!(framework)
+
+Wraps the coordinates of all atom and charge positions to be within the unit
+cell defined for the framework.
+"""
+function wrap_atoms_to_unit_cell!(framework::Framework)
+    framework.atoms.xf .= mod.(framework.atoms.xf, 1.0)
+    framework.charges.xf .= mod.(framework.charges.xf, 1.0)
+end
+
 """
     replicated_frame = replicate(framework, repfactors)
 
@@ -233,6 +486,8 @@ construct a new `Framework`. Note `replicate(framework, (1, 1, 1))` returns the
 - `replicated_frame::Framework`: Replicated framework
 """
 function replicate(framework::Framework, repfactors::Tuple{Int, Int, Int})
+    @assert ne(framework.bonds) == 0 @sprintf("The framework %s has bonds within it. Remove the bonds with `remove_bonds!` to replicate, and then use `infer_bonds(framework)` to recalculate bond information", framework.name)
+    assert_P1_symmetry(framework)
     # determine number of atoms in replicated framework
     n_atoms = size(framework.atoms.xf, 2) * repfactors[1] * repfactors[2] * repfactors[3]
 
@@ -266,9 +521,12 @@ function replicate(framework::Framework, repfactors::Tuple{Int, Int, Int})
 
     @assert (new_charges.n_charges == framework.charges.n_charges * prod(repfactors))
     @assert (new_atoms.n_atoms == framework.atoms.n_atoms * prod(repfactors))
-    return Framework(framework.name, new_box, new_atoms, new_charges)
+    return Framework(framework.name, new_box, new_atoms, new_charges,
+                     symmetry=deepcopy(framework.symmetry),
+                     space_group=framework.space_group, is_p1=framework.is_p1)
 end
 
+
 # doc string in Misc.jl
 function write_xyz(framework::Framework, filename::AbstractString;
                       comment::AbstractString="", center::Bool=false)
@@ -287,14 +545,14 @@ function write_xyz(framework::Framework, filename::AbstractString;
     write_xyz(atoms, x, filename, comment=comment)
 end
 write_xyz(framework::Framework; comment::AbstractString="", center::Bool=false) = write_xyz(
-    framework, 
+    framework,
     replace(replace(framework.name, ".cif" => ""), ".cssr" => "") * ".xyz",
     comment=comment, center=center)
 
 """
-    is_overlap = atom_overlap(framework; overlap_tol=0.1, verbose=false)
+    overlap = atom_overlap(framework; overlap_tol=0.1, verbose=false)
 
-Return true iff any two `Atoms` in the crystal overlap by calculating the distance
+Return true if any two `Atoms` in the crystal overlap by calculating the distance
 between every pair of atoms and ensuring distance is greater than
 `overlap_tol`. If verbose, print the pair of atoms which are culprits.
 
@@ -317,7 +575,7 @@ function atom_overlap(framework::Framework; overlap_tol::Float64=0.1, verbose::B
                         framework.box, overlap_tol)
                 overlap = true
                 if verbose
-                    @warn @sprintf("Atoms %d and %d in %s are less than %d Å apart.", i, j,
+                    @warn @sprintf("Atoms %d and %d in %s are less than %f Å apart.", i, j,
                         framework.name, overlap_tol)
                 end
             end
@@ -337,7 +595,7 @@ function charge_overlap(framework::Framework; overlap_tol::Float64=0.1, verbose:
                         framework.box, overlap_tol)
                 overlap = true
                 if verbose
-                    @warn @sprintf("Charges %d and %d in %s are less than %d Å apart.", i, j,
+                    @warn @sprintf("Charges %d and %d in %s are less than %f Å apart.", i, j,
                         framework.name, overlap_tol)
                 end
             end
@@ -350,7 +608,7 @@ end
 #   do overlap, and can then use that number to determine if they overlap or are repeats
 function _overlap(xf_1::Array{Float64, 1}, xf_2::Array{Float64, 1},
                   box::Box, overlap_tol::Float64)
-    dxf = xf_1 .- xf_2
+    dxf = mod.(xf_1, 1.0) .- mod.(xf_2, 1.0)
     nearest_image!(dxf)
     dxc = box.f_to_c * dxf
     return norm(dxc) < overlap_tol
@@ -362,7 +620,8 @@ end
 
 #TODO write tests for this! one with diff elements
 """
-    new_framework = remove_overlapping_atoms_and_charges(framework, overlap_tol=0.1, verbose=true)
+    new_framework = remove_overlapping_atoms_and_charges(framework, atom_overlap_tol=0.1, 
+                                                        charge_overlap_tol=0.1, verbose=true)
 
 Takes in a framework and returns a new framework with where overlapping atoms and overlapping
 charges were removed. i.e. if there is an overlapping pair, one in the pair is removed.
@@ -373,6 +632,7 @@ must be identical.
 - `framework::Framework`: The framework containing the crystal structure information
 - `atom_overlap_tol::Float64`: The minimum distance between two atoms that is tolerated
 - `charge_overlap_tol::Float64`: The minimum distance between two charges that is tolerated
+- `verbose::Bool`: Will print out information regarding the function call
 
 # Returns
 - `new_framework::Framework`: A new framework where identical atoms have been removed.
@@ -381,7 +641,7 @@ function remove_overlapping_atoms_and_charges(framework::Framework;
     atom_overlap_tol::Float64=0.1, charge_overlap_tol::Float64=0.1, verbose::Bool=true)
 
     atoms_to_keep = trues(framework.atoms.n_atoms)
-    charges_to_keep = trues(framework.atoms.n_atoms)
+    charges_to_keep = trues(framework.charges.n_charges)
 
     for i = 1:framework.atoms.n_atoms
         for j =  1:framework.atoms.n_atoms
@@ -441,7 +701,9 @@ function remove_overlapping_atoms_and_charges(framework::Framework;
     atoms = Atoms(framework.atoms.species[atoms_to_keep], atom_coords_to_keep)
     charges = Charges(framework.charges.q[charges_to_keep], charge_coords_to_keep)
 
-    new_framework = Framework(framework.name, framework.box, atoms, charges)
+    new_framework = Framework(framework.name, framework.box, atoms, charges,
+                              symmetry=deepcopy(framework.symmetry),
+                              space_group=framework.space_group, is_p1=framework.is_p1)
 
     @assert (! atom_overlap(new_framework, overlap_tol=atom_overlap_tol))
     @assert (! charge_overlap(new_framework, overlap_tol=charge_overlap_tol))
@@ -582,15 +844,271 @@ function crystal_density(framework::Framework)
 end
 
 """
-    write_cif(framework, filename)
+    simulation_ready_framework = apply_symmetry_rules(non_p1_framework;
+                                                check_charge_neutrality=true,
+                                                net_charge_tol=0.001,
+                                                check_atom_and_charge_overlap=true,
+                                                remove_overlap=false,
+                                                wrap_to_unit_cell=true)
+
+Convert a framework to P1 symmetry based on internal symmetry rules. This will
+return the new framework.
+
+# Arguments
+- `f::Framework`: The framework to be converted to P1 symmetry
+- `check_charge_neutrality::Bool`: check for charge neutrality
+- `net_charge_tol::Float64`: when checking for charge neutrality, throw an error if the absolute value of the net charge is larger than this value.
+- `check_atom_and_charge_overlap::Bool`: throw an error if overlapping atoms are detected.
+- `remove_overlap::Bool`: remove identical atoms automatically. Identical atoms are the same element atoms which overlap.
+- `wrap_to_unit_cell::Bool`: if true, enforce that fractional coords of atoms/charges are in [0,1]³ by mod(x, 1)
+
+# Returns
+- `P1_framework::Framework`: The framework after it has been converted to P1
+    symmetry. The new symmetry rules will be the P1 symmetry rules
+"""
+function apply_symmetry_rules(framework::Framework; check_charge_neutrality::Bool=true,
+                              net_charge_tol::Float64=0.001, check_atom_and_charge_overlap::Bool=true,
+                              remove_overlap::Bool=false, wrap_to_unit_cell::Bool=true)
+    if framework.is_p1
+        return framework
+    end
+    new_atom_xfs = Array{Float64, 2}(undef, 3, framework.atoms.n_atoms * size(framework.symmetry, 2))
+    new_charge_xfs = Array{Float64, 2}(undef, 3, framework.charges.n_charges * size(framework.symmetry, 2))
+    new_atom_species = Array{Symbol, 1}(undef, 0)
+    new_charge_qs = Array{Float64, 1}(undef, 0)
+
+    # for each symmetry rule
+    for i in 1:size(framework.symmetry, 2)
+        # loop over all atoms in lower level symmetry
+        sym_rule = eval.(Meta.parse.("(x, y, z) -> " .* framework.symmetry[:, i]))
+        for j in 1:size(framework.atoms.xf, 2)
+            # apply current symmetry rule to current atom for x, y, and z coordinates
+            current_atom_idx = (i - 1) * framework.atoms.n_atoms + j
+            new_atom_xfs[:, current_atom_idx] .= [Base.invokelatest.(
+                        sym_rule[k], framework.atoms.xf[:, j]...) for k in 1:3]
+        end
+        # loop over all charges in lower level symmetry
+        for j in 1:size(framework.charges.xf, 2)
+            # apply current symmetry rule to current atom for x, y, and z coordinates
+            current_charge_idx = (i - 1) * framework.charges.n_charges + j
+            new_charge_xfs[:, current_charge_idx] .= [Base.invokelatest.(
+                        sym_rule[k], framework.charges.xf[:, j]...) for k in 1:3]
+        end
+        # repeat charge_qs and atom_species for every symmetry applied
+        new_atom_species = [new_atom_species; framework.atoms.species]
+        new_charge_qs = [new_charge_qs; framework.charges.q]
+    end
+
+    new_framework = Framework(framework.name, framework.box,
+        Atoms(new_atom_species, new_atom_xfs),
+        Charges(new_charge_qs, new_charge_xfs))
+
+    if check_charge_neutrality
+        if ! charge_neutral(new_framework, net_charge_tol)
+            error(@sprintf("Framework %s is not charge neutral; net charge is %f e. Ignore
+            this error message by passing check_charge_neutrality=false or increasing the
+            net charge tolerance `net_charge_tol`\n",
+                            new_framework.name, total_charge(new_framework)))
+        end
+    end
+        
+    if wrap_to_unit_cell
+        wrap_atoms_to_unit_cell!(framework)
+    end
+
+    if remove_overlap
+        return remove_overlapping_atoms_and_charges(new_framework)
+    end
+
+    if check_atom_and_charge_overlap
+        if atom_overlap(new_framework) | charge_overlap(new_framework)
+            error(@sprintf("At least one pair of atoms/charges overlap in %s.
+            Consider passing `remove_overlap=true`\n", new_framework.name))
+        end
+    end
+
+    return new_framework
+end
+
+"""
+    symmetry_equal = is_symmetry_equal(framework1.symmetry, framework2.symmetry)
+
+Returns true if both symmetry rules can create the same set from the same set
+of coordinates. Returns false if they don't contain the same number of rules or
+if they create different sets of points.
+
+# Arguments
+- `sym1::Array{AbstractString, 2}`: Array of strings that represent
+    symmetry operations
+- `sym2::Array{AbstractString, 2}`: Array of strings that represent
+    symmetry operations
+
+# Returns
+- `is_equal::Bool`: True if they are the same set of symmetry rules
+    False if they are different
+"""
+function is_symmetry_equal(sym1::Array{AbstractString, 2}, sym2::Array{AbstractString, 2})
+    # need same number of symmetry operations
+    if size(sym1, 2) != size(sym2, 2)
+        return false
+    end
+    # define a test array that operations will be performed on
+    test_array = [0.0 0.25 0.0  0.0  0.0  0.25 0.25 0.25;
+                  0.0 0.0  0.25 0.0  0.25 0.0  0.25 0.25;
+                  0.0 0.0  0.0  0.25 0.25 0.25 0.25 0.25]
+    # set up both arrays for storing replicated coords
+    sym1_applied_to_test = Array{Float64, 2}(undef, 3, 0)
+    sym2_applied_to_test = Array{Float64, 2}(undef, 3, 0)
+
+    # loop over all positions in the test_array
+    for i in 1:size(test_array, 2)
+        # loop over f1 symmetry rules
+        for j in 1:size(sym1, 2)
+            sym1_applied_to_test = [sym1_applied_to_test [Base.invokelatest.(
+                eval(Meta.parse("(x, y, z) -> " * sym1[k, j])), test_array[:, i]...) for k in 1:3]]
+        end
+        # loop over f2 symmetry rules
+        for j in 1:size(sym2, 2)
+            sym2_applied_to_test = [sym2_applied_to_test [Base.invokelatest.(
+                eval(Meta.parse("(x, y, z) -> " * sym2[k, j])), test_array[:, i]...) for k in 1:3]]
+        end
+    end
+
+    # convert to sets for using issetequal, symmetry rules might be in a a different order
+    sym1_set = Set([sym1_applied_to_test[:, i] for i in 1:size(sym1_applied_to_test, 2)])
+    sym2_set = Set([sym2_applied_to_test[:, i] for i in 1:size(sym2_applied_to_test, 2)])
+
+    # return if the sets of coords are equal
+    return issetequal(sym1_set, sym2_set)
+end
+
+"""
+    assert_P1_symmetry(framework_to_be_tested)
+
+Used for making sure that a framework is in P1 symmetry before running
+simulations on it. If the structure is not in P1, this will throw an
+AssertionError.
+-`framework::Framework`: The framework to be tested for P1 symmetry
+"""
+function assert_P1_symmetry(framework::Framework)
+    @assert framework.is_p1 @sprintf("The framework %s is not in P1 symmetry.\n
+                                     try running:\n
+                                     \tframework_p1 = apply_symmetry_rules(framework)\n
+                                     and pass `framework_p1` into this simulation",
+                                    framework.name)
+end
+
+"""
+    remove_bonds!(framework)
+
+Remove all bonds from a framework structure.
+
+# Arguments
+-`framework::Framework`: the framework that bonds wil be removed from
+"""
+function remove_bonds!(framework::Framework)
+    while ne(framework.bonds) > 0
+        rem_edge!(framework.bonds, collect(edges(framework.bonds))[1].src, collect(edges(framework.bonds))[1].dst)
+    end
+end
+
+"""
+    infer_bonds!(framework, include_bonds_across_periodic_boundaries, 
+                    bonding_rules=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)])
+
+Populate the bonds in the framework object based on the bonding rules. If a
+pair doesn't have a suitable rule then they will not be considered bonded. 
+
+`:*` is considered a wildcard and can be substituted for any species. It is a
+good idea to include a bonding rule between two `:*` to allow any atoms to bond
+as long as they are close enough.
+
+The bonding rules are hierarchical, i.e. the first bonding rule takes precedence over the latter ones.
+
+# Arguments
+-`framework::Framework`: The framework that bonds will be added to
+-`include_bonds_across_periodic_boundaries::Bool`: Whether to check across the
+    periodic boundary when calculating bonds
+-`bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will
+    be used to fill the bonding information. They are applied in the order that
+    they appear.
+"""
+function infer_bonds!(framework::Framework, include_bonds_across_periodic_boundaries::Bool,
+                      bonding_rules::Array{BondingRule, 1}=
+                      [BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)])
+    @assert ne(framework.bonds) == 0 @sprintf("The framework %s already has bonds. Remove them with the `remove_bonds!` function before inferring new ones.", framework.name)
+
+    # loop over every atom
+    for i in 1:framework.atoms.n_atoms
+        # loop over every unique pair of atoms
+        for j in i+1:framework.atoms.n_atoms
+            if is_bonded(framework, i, j, bonding_rules; include_bonds_across_periodic_boundaries=include_bonds_across_periodic_boundaries)
+                add_edge!(framework.bonds, i, j)
+            end
+        end
+    end
+end
+
+"""
+    bonds_equal = compare_bonds_in_framework(framework1, framework2, atol=0.0)
+
+Returns whether the bonds defined in framework1 are the same as the bonds
+defined in framework2. It checks whether the atoms in the same positions
+have the same bonds.
+
+# Arguments
+-`framework1::Framework`: The first framework
+-`framework2::Framework`: The second framework
+-`atol::Float64`: absolute tolerance for the comparison of coordinates in the framework
+
+# Returns
+-`bonds_equal::Bool`: Wether the bonds in framework1 and framework2 are equal
+"""
+function compare_bonds_in_framework(fi::Framework, fj::Framework; atol::Float64=0.0)
+    if ne(fi.bonds) != ne(fj.bonds)
+        return false
+    end
+
+    num_in_common = 0
+    for edge_i in collect(edges(fi.bonds))
+        for edge_j in collect(edges(fj.bonds))
+            # either the bond matches going src-src dst-dst
+            if  (fi.atoms.species[edge_i.src] == fj.atoms.species[edge_j.src] &&
+                 fi.atoms.species[edge_i.dst] == fj.atoms.species[edge_j.dst] &&
+                 isapprox(fi.atoms.xf[:, edge_i.src], fj.atoms.xf[:, edge_j.src]; atol=atol) &&
+                 isapprox(fi.atoms.xf[:, edge_i.dst], fj.atoms.xf[:, edge_j.dst]; atol=atol)) ||
+                # or the bond matches going src-dst dst-src
+                (fi.atoms.species[edge_i.src] == fj.atoms.species[edge_j.dst] &&
+                 fi.atoms.species[edge_i.dst] == fj.atoms.species[edge_j.src] &&
+                 isapprox(fi.atoms.xf[:, edge_i.src], fj.atoms.xf[:, edge_j.dst]; atol=atol) &&
+                 isapprox(fi.atoms.xf[:, edge_i.dst], fj.atoms.xf[:, edge_j.src]; atol=atol))
+                num_in_common += 1
+                break
+            end
+        end
+    end
+    return num_in_common == ne(fi.bonds) && num_in_common == ne(fj.bonds)
+end
+
+"""
+    write_cif(framework, filename; fractional=true)
 
 Write a `framework::Framework` to a .cif file with `filename::AbstractString`. If `filename` does
 not include the .cif extension, it will automatically be added.
 """
-function write_cif(framework::Framework, filename::AbstractString)
+function write_cif(framework::Framework, filename::AbstractString; fractional::Bool=true)
     if charged(framework) && (framework.atoms.n_atoms != framework.charges.n_charges)
         error("write_cif assumes equal numbers of Charges and Atoms (or zero charges)")
     end
+
+    # create dictionary for tracking label numbers
+    label_numbers = Dict{Symbol, Int}()
+    for atom in framework.atoms.species
+        if !haskey(label_numbers, atom)
+            label_numbers[atom] = 1
+        end
+    end
+
     # append ".cif" to filename if it doesn't already have the extension
     if ! occursin(".cif", filename)
         filename *= ".cif"
@@ -604,23 +1122,32 @@ function write_cif(framework::Framework, filename::AbstractString)
         @printf(cif_file, "data_%s_PM\n", split(framework.name, ".")[1])
     end
 
-    @printf(cif_file, "_symmetry_space_group_name_H-M   'P 1'\n")
+    @printf(cif_file, "_symmetry_space_group_name_H-M\t'%s'\n", framework.space_group)
 
-    @printf(cif_file, "_cell_length_a %f\n", framework.box.a)
-    @printf(cif_file, "_cell_length_b %f\n", framework.box.b)
-    @printf(cif_file, "_cell_length_c %f\n", framework.box.c)
+    @printf(cif_file, "_cell_length_a\t%f\n", framework.box.a)
+    @printf(cif_file, "_cell_length_b\t%f\n", framework.box.b)
+    @printf(cif_file, "_cell_length_c\t%f\n", framework.box.c)
 
-    @printf(cif_file, "_cell_angle_alpha %f\n", framework.box.α * 180.0 / pi)
-    @printf(cif_file, "_cell_angle_beta %f\n", framework.box.β * 180.0 / pi)
-    @printf(cif_file, "_cell_angle_gamma %f\n", framework.box.γ * 180.0 / pi)
+    @printf(cif_file, "_cell_angle_alpha\t%f\n", framework.box.α * 180.0 / pi)
+    @printf(cif_file, "_cell_angle_beta\t%f\n", framework.box.β * 180.0 / pi)
+    @printf(cif_file, "_cell_angle_gamma\t%f\n", framework.box.γ * 180.0 / pi)
 
     @printf(cif_file, "_symmetry_Int_Tables_number 1\n\n")
-    @printf(cif_file, "loop_\n_symmetry_equiv_pos_as_xyz\n 'x, y, z'\n\n")
+    @printf(cif_file, "loop_\n_symmetry_equiv_pos_as_xyz\n")
+    for i in 1:size(framework.symmetry, 2)
+        @printf(cif_file, "'%s,%s,%s'\n", framework.symmetry[:, i]...)
+    end
+    @printf(cif_file, "\n")
 
-    @printf(cif_file, "loop_\n_atom_site_label\n")
-    @printf(cif_file, "_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n")
+    @printf(cif_file, "loop_\n_atom_site_label\n_atom_site_type_symbol\n")
+    if fractional
+        @printf(cif_file, "_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n")
+    else
+        @printf(cif_file, "_atom_site_Cartn_x\n_atom_site_Cartn_y\n_atom_site_Cartn_z\n")
+    end
     @printf(cif_file, "_atom_site_charge\n")
 
+    idx_to_label = Array{AbstractString, 1}(undef, framework.atoms.n_atoms)
     for i = 1:framework.atoms.n_atoms
         q = 0.0
         if charged(framework)
@@ -629,13 +1156,73 @@ function write_cif(framework::Framework, filename::AbstractString)
                 error("write_cif assumes charges correspond to LJspheres")
             end
         end
-        @printf(cif_file, "%s %f %f %f %f\n", framework.atoms.species[i], 
-                framework.atoms.xf[1, i], framework.atoms.xf[2, i], 
-                framework.atoms.xf[3, i], q)
-     end
-     close(cif_file)
+        # print label and type symbol
+        @printf(cif_file, "%s\t%s\t", string(framework.atoms.species[i]) *
+                string(label_numbers[framework.atoms.species[i]]),
+                framework.atoms.species[i])
+        # store label for this atom idx
+        idx_to_label[i] = string(framework.atoms.species[i]) *
+                    string(label_numbers[framework.atoms.species[i]])
+        # increment label
+        label_numbers[framework.atoms.species[i]] += 1
+        if fractional
+            @printf(cif_file, "%f\t%f\t%f\t%f\n", framework.atoms.xf[:, i]..., q)
+        else
+            
+            @printf(cif_file, "%f\t%f\t%f\t%f\n", (framework.box.f_to_c * framework.atoms.xf[:, i])..., q)
+        end
+    end
+
+    # only print bond information if it is in the framework
+    if ne(framework.bonds) > 0
+        # print column names for bond information
+        @printf(cif_file, "\nloop_\n_geom_bond_atom_site_label_1\n_geom_bond_atom_site_label_2\n_geom_bond_distance\n")
+
+        for edge in collect(edges(framework.bonds))
+            dxf = framework.atoms.xf[:, edge.src] - framework.atoms.xf[:, edge.dst]
+            nearest_image!(dxf)
+            @printf(cif_file, "%s\t%s\t%0.5f\n", idx_to_label[edge.src], idx_to_label[edge.dst],
+                    norm(dxf))
+        end
+    end
+    close(cif_file)
 end
 
+"""
+    write_bond_information(framework, filename)
+    write_bond_information(framework)
+
+Writes the bond information from a framework to the selected filename.
+
+# Arguments
+-`framework::Framework`: The framework to have its bonds written to a vtk file
+-`filename::AbstractString`: The filename the bond information will be saved to. If left out, will default to framework name.
+"""
+function write_bond_information(framework::Framework, filename::AbstractString)
+    if ne(framework.bonds) == 0
+        @warn("Framework %s has no bonds present. To get bonding information for this framework run `infer_bonds!` with an array of bonding rules\n", framework.name)
+    end
+    if ! occursin(".vtk", filename)
+        filename *= ".vtk"
+    end
+
+    vtk_file = open(filename, "w")
+
+    @printf(vtk_file, "# vtk DataFile Version 2.0\n%s bond information\nASCII\nDATASET POLYDATA\nPOINTS %d double\n", framework.name, nv(framework.bonds))
+
+    for i = 1:framework.atoms.n_atoms
+        @printf(vtk_file, "%0.5f\t%0.5f\t%0.5f\n", (framework.box.f_to_c * framework.atoms.xf[:, i])...)
+    end
+    @printf(vtk_file, "\nLINES %d %d\n", ne(framework.bonds), 3 * ne(framework.bonds))
+    for edge in collect(edges(framework.bonds))
+        @printf(vtk_file, "2\t%d\t%d\n", edge.src - 1, edge.dst - 1)
+    end
+    close(vtk_file)
+    @printf("Saving bond information for framework %s to %s.\n", framework.name, joinpath(pwd(), filename))
+end
+
+write_bond_information(framework::Framework) = write_bond_information(framework, framework.name * "_bonds.vtk")
+
 """
     new_framework = assign_charges(framework, charges, net_charge_tol=1e-5)
 
@@ -704,7 +1291,9 @@ function assign_charges(framework::Framework, charges::Union{Dict{Symbol, Float6
     charges = Charges(charge_vals, charge_coords)
 
     # construct new framework
-    new_framework = Framework(framework.name, framework.box, framework.atoms, charges)
+    new_framework = Framework(framework.name, framework.box, framework.atoms, charges,
+                              symmetry=deepcopy(framework.symmetry),
+                              space_group=framework.space_group, is_p1=framework.is_p1)
 
     # check for charge neutrality
     if abs(total_charge(new_framework)) > net_charge_tol
@@ -716,19 +1305,96 @@ function assign_charges(framework::Framework, charges::Union{Dict{Symbol, Float6
     return new_framework
 end
 
+"""
+    are_atoms_bonded = is_bonded(framework, i, j, bonding_rules=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)],
+                                 include_bonds_across_periodic_boundaries=true)
+
+Checks to see if atoms `i` and `j` in `framework` are bonded according to the `bonding_rules`.
+
+# Arguments
+-`framework::Framework`: The framework that bonds will be added to
+-`i::Int`: Index of the first atom
+-`j::Int`: Index of the second atom
+-`bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will
+    be used to fill the bonding information. They are applied in the order that
+    they appear.
+-`include_bonds_across_periodic_boundaries::Bool`: Whether to check across the
+    periodic boundary when calculating bonds
+
+# Returns
+-`are_atoms_bonded::Bool`: Whether atoms `i` and `j` are bonded according to `bonding_rules`
+
+"""
+function is_bonded(framework::Framework, i::Int64, j::Int64,
+                   bonding_rules::Array{BondingRule, 1}=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)];
+                   include_bonds_across_periodic_boundaries::Bool=true)
+    species_i = framework.atoms.species[i]
+    species_j = framework.atoms.species[j]
+    x_i = framework.atoms.xf[:, i]
+    x_j = framework.atoms.xf[:, j]
+    # loop over possible bonding rules
+    for br in bonding_rules
+        # determine if the atom species correspond to the species in `bonding_rules`
+        species_match = false
+        if br.species_i == :* && br.species_j == :*
+            species_match = true
+        elseif br.species_i == :* && (species_i == br.species_j || species_j == br.species_j)
+            species_match = true
+        elseif br.species_j == :* && (species_i == br.species_i || species_j == br.species_j)
+            species_match = true
+        elseif (species_i == br.species_i && species_j == br.species_j) || (species_j == br.species_i && species_i == br.species_j)
+            species_match = true
+        end
+
+        if species_match
+            # determine if the atoms are close enough to bond
+            dxf = x_i - x_j
+            if include_bonds_across_periodic_boundaries
+                nearest_image!(dxf)
+            end
+            cartesian_dist_between_atoms = norm(framework.box.f_to_c * dxf)
+            if br.min_dist < cartesian_dist_between_atoms && br.max_dist > cartesian_dist_between_atoms
+                return true
+            end
+        end
+    end
+    return false
+end
+
+
 function Base.show(io::IO, framework::Framework)
     println(io, "Name: ", framework.name)
     println(io, framework.box)
 	@printf(io, "Number of atoms = %d\n", framework.atoms.n_atoms)
 	@printf(io, "Number of charges = %d\n", framework.charges.n_charges)
     println(io, "Chemical formula: ", chemical_formula(framework))
+    @printf(io, "Space Group: %s\n", framework.space_group)
+    @printf(io, "Symmetry Operations:\n")
+    for i in 1:size(framework.symmetry, 2)
+        @printf(io, "\t'%s, %s, %s'\n", framework.symmetry[:, i]...)
+    end
 end
 
-function Base.isapprox(f1::Framework, f2::Framework; checknames::Bool=false)
+function has_same_sets_of_atoms_and_charges(f1::Framework, f2::Framework; atol::Float64=1e-6, checknames::Bool=false)
     names_flag = f1.name == f2.name
     if checknames && (! names_flag)
         return false
     end
+    box_flag = isapprox(f1.box, f2.box)
+    if f1.charges.n_charges != f2.charges.n_charges
+        return false
+    end
+    if f1.atoms.n_atoms != f2.atoms.n_atoms
+        return false
+    end
+    charges_flag = has_same_set_of_charges(f1.charges, f2.charges; atol=atol)
+    atoms_flag = has_same_set_of_atoms(f1.atoms, f2.atoms; atol=atol)
+    symmetry_flag = is_symmetry_equal(f1.symmetry, f2.symmetry)
+    return box_flag && charges_flag && atoms_flag && symmetry_flag
+end
+
+
+function Base.isapprox(f1::Framework, f2::Framework)
     box_flag = isapprox(f1.box, f2.box)
     if f1.charges.n_charges != f2.charges.n_charges
         return false
@@ -738,5 +1404,43 @@ function Base.isapprox(f1::Framework, f2::Framework; checknames::Bool=false)
     end
     charges_flag = isapprox(f1.charges, f2.charges)
     atoms_flag = isapprox(f1.atoms, f2.atoms)
-    return box_flag && charges_flag && atoms_flag
+    symmetry_flag = is_symmetry_equal(f1.symmetry, f2.symmetry)
+    return box_flag && charges_flag && atoms_flag && symmetry_flag
+end
+
+function Base.:+(frameworks::Framework...; check_overlap::Bool=true)
+    new_framework = deepcopy(frameworks[1])
+    for (i, f) in enumerate(frameworks)
+        if i == 1
+            continue
+        end
+        @assert isapprox(new_framework.box, f.box) @sprintf("Framework %s has a different box\n", f.name)
+        @assert is_symmetry_equal(new_framework.symmetry, f.symmetry) @sprintf("Framework %s has different symmetry rules\n", f.name)
+        @assert new_framework.space_group == f.space_group
+
+        new_atoms = new_framework.atoms + f.atoms
+        new_charges = new_framework.charges + f.charges
+
+        nf_n_atoms = new_framework.atoms.n_atoms
+        add_vertices!(new_framework.bonds, nf_n_atoms)
+        for edge in collect(edges(f.bonds))
+            add_edge!(new_framework.bonds, nf_n_atoms + edge.src, nf_n_atoms + edge.dst)
+        end
+
+        new_framework = Framework(split(new_framework.name, ".")[1] * "_" * split(f.name, ".")[1],
+                                 new_framework.box, new_atoms, new_charges,
+                                 symmetry=new_framework.symmetry,space_group=new_framework.space_group,
+                                 is_p1=new_framework.is_p1, bonds=new_framework.bonds)
+    end
+    if check_overlap
+        if atom_overlap(new_framework)
+            @warn "This new framework has overlapping atoms, use:\n`remove_overlapping_atoms_and_charges(framework)`\nto remove them"
+        end
+
+        if charge_overlap(new_framework)
+            @warn "This new framework has overlapping charges, use:\n`remove_overlapping_atoms_and_charges(framework)`\nto remove them"
+        end
+    end
+
+    return new_framework
 end
diff --git a/src/GCMC.jl b/src/GCMC.jl
index 71ab230ce..e56c94b87 100644
--- a/src/GCMC.jl
+++ b/src/GCMC.jl
@@ -129,6 +129,9 @@ end
                                   verbose=true, ewald_precision=1e-6, eos=:ideal,
                                   load_checkpoint_file=false, checkpoint=Dict(),
                                   write_checkpoints=false, checkpoint_frequency=50,
+                                  write_adsorbate_snapshots=false,
+                                  snapshot_frequency=1, calculate_density_grid=false,
+                                  density_grid_dx=1.0, density_grid_species=nothing,
                                   filename_comment="", show_progress_bar=false)
 
 Run a set of grand-canonical (μVT) Monte Carlo simulations in series. Arguments are the
@@ -157,6 +160,10 @@ function stepwise_adsorption_isotherm(framework::Framework,
                                       density_grid_dx::Float64=1.0, 
                                       density_grid_species::Union{Nothing, Symbol}=nothing,
                                       filename_comment::AbstractString="")
+
+    # simulation only works if framework is in P1
+    assert_P1_symmetry(framework)
+
     results = Dict{String, Any}[] # push results to this array
     molecules = Molecule[] # initiate with empty framework
     for (i, pressure) in enumerate(pressures)
@@ -189,6 +196,9 @@ end
                                   verbose=true, ewald_precision=1e-6, eos=:ideal, 
                                   load_checkpoint_file=false, checkpoint=Dict(), 
                                   write_checkpoints=false, checkpoint_frequency=50,
+                                  write_adsorbate_snapshots=false,
+                                  snapshot_frequency=1, calculate_density_grid=false,
+                                  density_grid_dx=1.0, density_grid_species=nothing,
                                   filename_comment="", show_progress_bar=false)
 
 Run a set of grand-canonical (μVT) Monte Carlo simulations in parallel. Arguments are the
@@ -213,6 +223,10 @@ function adsorption_isotherm(framework::Framework,
                              density_grid_dx::Float64=1.0, 
                              density_grid_species::Union{Nothing, Symbol}=nothing,
                              filename_comment::AbstractString="")
+
+    # simulation only works if framework is in P1
+    assert_P1_symmetry(framework)
+
     # make a function of pressure only to facilitate uses of `pmap`
     run_pressure(pressure::Float64) = gcmc_simulation(framework, molecule, temperature,
                                                       pressure, ljforcefield,
@@ -252,6 +266,9 @@ end
                                          load_checkpoint_file=false,
                                          show_progress_bar=false, checkpoint=Dict(),
                                          write_checkpoints=false, checkpoint_frequency=100,
+                                         write_adsorbate_snapshots=false,
+                                         snapshot_frequency=1, calculate_density_grid=false,
+                                         density_grid_dx=1.0, density_grid_species=nothing,
                                          filename_comment="")
 
 Runs a grand-canonical (μVT) Monte Carlo simulation of the adsorption of a molecule in a
@@ -307,6 +324,9 @@ function gcmc_simulation(framework::Framework, molecule_::Molecule, temperature:
     calculate_density_grid::Bool=false, density_grid_dx::Float64=1.0, 
     density_grid_species::Union{Nothing, Symbol}=nothing, filename_comment::AbstractString="")
 
+    # simulation only works if framework is in P1
+    assert_P1_symmetry(framework)
+
     start_time = time()
     # to avoid changing the outside object `molecule_` inside this function, we make
     #  a deep copy of it here. this serves as a template to copy when we insert a new molecule.
diff --git a/src/Grid.jl b/src/Grid.jl
index b0b0427b4..bd5f7286e 100644
--- a/src/Grid.jl
+++ b/src/Grid.jl
@@ -246,6 +246,7 @@ function energy_grid(framework::Framework, molecule::Molecule, ljforcefield::LJF
                      n_pts::Tuple{Int, Int, Int}=(50,50,50), n_rotations::Int=1000, 
                      temperature::Float64=NaN, units::Symbol=:kJ_mol, center::Bool=false,
                      verbose::Bool=true)
+    assert_P1_symmetry(framework)
     if ! (units in [:kJ_mol, :K])
         error("Pass :kJ_mol or :K for units of kJ/mol or K, respectively.")
     end
@@ -626,6 +627,7 @@ function compute_accessibility_grid(framework::Framework, probe::Molecule, force
                                     energy_tol::Float64=10.0,  energy_units::Symbol=:kJ_mol,
                                     verbose::Bool=true, write_b4_after_grids::Bool=true, 
                                     block_inaccessible_pockets::Bool=true)
+    assert_P1_symmetry(framework)
     if verbose
         printstyled(@sprintf("Computing accessibility grid of %s using %f %s potential energy tol and %s probe...\n",
             framework.name, energy_tol, energy_units, probe.species), color=:green)
diff --git a/src/Henry.jl b/src/Henry.jl
index a08216944..bbbd245be 100644
--- a/src/Henry.jl
+++ b/src/Henry.jl
@@ -46,6 +46,10 @@ function henry_coefficient(framework::Framework, molecule_::Molecule, temperatur
                            write_checkpoint::Bool=false, load_checkpoint::Bool=false,
                            checkpoint_frequency::Int=10000,
                            accessibility_grid::Union{Nothing, Grid{Bool}}=nothing)
+
+    # simulation only works if framework is in P1
+    assert_P1_symmetry(framework)
+
     time_start = time()
     if verbose
         print("Simulating Henry coefficient of ")
diff --git a/src/Matter.jl b/src/Matter.jl
index e54e06513..f8282e489 100644
--- a/src/Matter.jl
+++ b/src/Matter.jl
@@ -22,6 +22,14 @@ Atoms(species::Array{Symbol, 1}, xf::Array{Float64, 2}) = Atoms(size(xf, 2), spe
 
 Base.isapprox(a1::Atoms, a2::Atoms) = (a1.species == a2.species) && isapprox(a1.xf, a2.xf)
 
+function has_same_set_of_atoms(a1::Atoms, a2::Atoms; atol::Float64=1e-6)
+    return issetequal(
+    Set([(round.(a1.xf[:, i], digits=Int(abs(log10(atol)))), a1.species[i]) for i in 1:a1.n_atoms]),
+    Set([(round.(a2.xf[:, i], digits=Int(abs(log10(atol)))), a2.species[i]) for i in 1:a2.n_atoms])) 
+end
+
+Base.:+(a1::Atoms, a2::Atoms) = Atoms(a1.n_atoms + a2.n_atoms, [a1.species; a2.species], [a1.xf a2.xf])
+
 """
 Data structure holds a set of point charges and their positions in fractional coordinates.
 
@@ -44,4 +52,12 @@ end
 # compute n_charges automatically from array sizes
 Charges(q::Array{Float64, 1}, xf::Array{Float64, 2}) = Charges(size(xf, 2), q, xf)
 
-Base.isapprox(c1::Charges, c2::Charges) = (isapprox(c1.q, c2.q) && isapprox(c1.xf, c2.xf))
+Base.isapprox(c1::Charges, c2::Charges) = isapprox(c1.q, c2.q) && isapprox(c1.xf, c2.xf)
+
+function has_same_set_of_charges(c1::Charges, c2::Charges; atol::Float64=1e-6)
+    return issetequal(
+    Set([(round.(c1.xf[:, i], digits=Int(abs(log10(atol)))), c1.q[i]) for i in 1:c1.n_charges]),
+    Set([(round.(c2.xf[:, i], digits=Int(abs(log10(atol)))), c2.q[i]) for i in 1:c2.n_charges]))
+end
+
+Base.:+(c1::Charges, c2::Charges) = Charges(c1.n_charges + c2.n_charges, [c1.q; c2.q], [c1.xf c2.xf])
diff --git a/src/PorousMaterials.jl b/src/PorousMaterials.jl
index 4b7542cd7..994c1228a 100644
--- a/src/PorousMaterials.jl
+++ b/src/PorousMaterials.jl
@@ -104,7 +104,7 @@ export
     Box, replicate, UnitCube, write_vtk, inside,
 
     # Matter.jl
-    Atoms, Charges,
+    Atoms, Charges, has_same_set_of_atoms, has_same_set_of_charges,
 
     # NearestImage.jl
     nearest_image!, nearest_r², nearest_r,
@@ -113,9 +113,13 @@ export
     read_xyz, read_cpk_colors, read_atomic_radii, write_xyz, fit_adsorption_isotherm,
 
     # Crystal.jl
-    Framework, read_crystal_structure_file, remove_overlapping_atoms_and_charges,
-    strip_numbers_from_atom_labels!, chemical_formula, molecular_weight, crystal_density,
-    construct_box, replicate, read_atomic_masses, charged, write_cif, assign_charges,
+    Framework, BondingRule, read_crystal_structure_file,
+    remove_overlapping_atoms_and_charges, strip_numbers_from_atom_labels!,
+    chemical_formula, molecular_weight, crystal_density, construct_box,
+    replicate, read_atomic_masses, charged, write_cif, assign_charges,
+    is_symmetry_equal, apply_symmetry_rules, assert_P1_symmetry, infer_bonds!,
+    remove_bonds!, compare_bonds_in_framework, wrap_atoms_to_unit_cell!,
+    write_bond_information, is_bonded, default_bondingrules, has_same_sets_of_atoms_and_charges,
 
     # Molecules.jl
     Molecule, n_atoms, set_fractional_coords!, translate_by!, outside_box,
diff --git a/test/crystal_test.jl b/test/crystal_test.jl
index 88a07619c..a10774dc9 100644
--- a/test/crystal_test.jl
+++ b/test/crystal_test.jl
@@ -3,6 +3,7 @@ module Crystal_Test
 using PorousMaterials
 using OffsetArrays
 using LinearAlgebra
+using LightGraphs
 using Test
 using JLD2
 using Statistics
@@ -14,34 +15,167 @@ using Random
     @test framework.name == "test_structure2.cif"
     @test isapprox(framework.box, Box(10.0, 20.0, 30.0, 90*π/180, 45*π/180, 120*π/180))
     @test framework.atoms.n_atoms == 2
-    @test isapprox(framework.atoms, Atoms([:Ca, :O], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
-    @test isapprox(framework.charges, Charges([1.0, -1.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
+    @test has_same_set_of_atoms(framework.atoms, Atoms([:Ca, :O], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
+    @test has_same_set_of_charges(framework.charges, Charges([1.0, -1.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
     new_frame = assign_charges(framework, Dict(:Ca => -2.0, :O => 2.0))
-    @test isapprox(new_frame.charges, Charges([-2.0, 2.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
+    @test has_same_set_of_charges(new_frame.charges, Charges([-2.0, 2.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
     new_frame = assign_charges(framework, [4.0, -4.0])
-    @test isapprox(new_frame.charges, Charges([4.0, -4.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
+    @test has_same_set_of_charges(new_frame.charges, Charges([4.0, -4.0], [0.2 0.6; 0.5 0.3; 0.7 0.1]))
     @test charged(framework)
     @test chemical_formula(framework) == Dict(:Ca => 1, :O => 1)
     @test molecular_weight(framework) ≈ 15.9994 + 40.078
     # same as test_structure.cif but with overlapping atoms.
     framework2 = Framework("test_structure2B.cif", remove_overlap=true, check_charge_neutrality=false)
     strip_numbers_from_atom_labels!(framework2)
-    @test isapprox(framework.atoms, framework2.atoms) && isapprox(framework.charges, framework2.charges)
+    @test has_same_set_of_atoms(framework.atoms, framework2.atoms) && has_same_set_of_charges(framework.charges, framework2.charges)
 
     # test .cif writer; write, read in, assert equal
+    #   more write_cif tests below with symmetry tests
     write_cif(framework, joinpath("data", "crystals", "rewritten_test_structure2.cif"))
     framework_rewritten = Framework("rewritten_test_structure2.cif")
-    @test isapprox(framework, framework_rewritten)
+    strip_numbers_from_atom_labels!(framework_rewritten)
+    write_cif(framework, joinpath("data", "crystals", "rewritten_test_structure2_cartn.cif"); fractional=false)
+    framework_rewritten_cartn = Framework("rewritten_test_structure2_cartn.cif")
+    strip_numbers_from_atom_labels!(framework_rewritten_cartn)
+    @test has_same_sets_of_atoms_and_charges(framework, framework_rewritten)
+    @test has_same_sets_of_atoms_and_charges(framework, framework_rewritten_cartn)
+
+    # test .cif reader for non-P1 symmetry
+    #   no atoms should overlap
+    #   should place atoms in the same positions as the P1 conversion using
+    #       openBabel
+    non_P1_framework = Framework("symmetry_test_structure.cif")
+    strip_numbers_from_atom_labels!(non_P1_framework)
+    non_P1_cartesian = Framework("symmetry_test_structure_cartn.cif")
+    strip_numbers_from_atom_labels!(non_P1_cartesian)
+    P1_framework = Framework("symmetry_test_structure_P1.cif")
+    strip_numbers_from_atom_labels!(P1_framework)
+
+    # wrap all atoms and charges to be within the unit cell
+    non_P1_framework.atoms.xf .= mod.(non_P1_framework.atoms.xf, 1.0)
+    non_P1_framework.charges.xf .= mod.(non_P1_framework.charges.xf, 1.0)
+
+    non_P1_cartesian.atoms.xf .= mod.(non_P1_cartesian.atoms.xf, 1.0)
+    non_P1_cartesian.charges.xf .= mod.(non_P1_cartesian.charges.xf, 1.0)
+
+    P1_framework.atoms.xf .= mod.(P1_framework.atoms.xf, 1.0)
+    P1_framework.charges.xf .= mod.(P1_framework.charges.xf, 1.0)
+
+    @test has_same_sets_of_atoms_and_charges(non_P1_framework, P1_framework; atol=1e-2) 
+    # test that fractional and cartesian produce same results
+    @test has_same_sets_of_atoms_and_charges(non_P1_framework, non_P1_cartesian; atol=1e-2)
+    # test that cartesian and P1 produce same results
+    @test has_same_sets_of_atoms_and_charges(non_P1_cartesian, P1_framework; atol=1e-2)
+    # test that incorrect file formats throw proper errors
+    @test_throws ErrorException Framework("non_P1_no_symmetry.cif")
+    # test that a file with no atoms throws error
+    @test_throws ErrorException Framework("no_atoms.cif")
+    # test that these carry the is_p1 flag
+    @test non_P1_framework.is_p1
+    @test non_P1_cartesian.is_p1
+    @test P1_framework.is_p1
+
+    # test reading in non-P1 then applying symmetry later
+    # read in the same files as above, then convert to P1, then compare
+    non_P1_framework_symmetry = Framework("symmetry_test_structure.cif", convert_to_p1=false)
+    strip_numbers_from_atom_labels!(non_P1_framework_symmetry)
+    non_P1_cartesian_symmetry = Framework("symmetry_test_structure_cartn.cif", convert_to_p1=false)
+    strip_numbers_from_atom_labels!(non_P1_cartesian_symmetry)
+
+    # make sure these frameworks are not in P1 symmetry when convert_to_p1 is
+    #   set to false
+    @test ! non_P1_framework_symmetry.is_p1
+    @test ! non_P1_cartesian_symmetry.is_p1
+
+    # test write_cif in non_p1 symmetry
+    write_cif(non_P1_framework_symmetry, joinpath("data", "crystals", "rewritten_symmetry_test_structure.cif"))
+    # keep this in cartesian to test both
+    write_cif(non_P1_cartesian_symmetry, joinpath("data", "crystals", "rewritten_symmetry_test_structure_cartn.cif"), fractional=false)
+    rewritten_non_p1_fractional = Framework("rewritten_symmetry_test_structure.cif"; convert_to_p1=false)
+    strip_numbers_from_atom_labels!(rewritten_non_p1_fractional)
+    rewritten_non_p1_cartesian = Framework("rewritten_symmetry_test_structure_cartn.cif"; convert_to_p1=false)
+    strip_numbers_from_atom_labels!(rewritten_non_p1_cartesian)
+
+    @test isapprox(rewritten_non_p1_fractional, non_P1_framework_symmetry)
+    @test isapprox(rewritten_non_p1_cartesian, non_P1_cartesian_symmetry)
+
+    non_P1_framework_symmetry = apply_symmetry_rules(non_P1_framework_symmetry)
+    non_P1_cartesian_symmetry = apply_symmetry_rules(non_P1_cartesian_symmetry)
+
+    # wrap all atoms and charges to be within the unit cell
+    non_P1_framework_symmetry.atoms.xf .= mod.(non_P1_framework_symmetry.atoms.xf, 1.0)
+    non_P1_framework_symmetry.charges.xf .= mod.(non_P1_framework_symmetry.charges.xf, 1.0)
+
+    non_P1_cartesian_symmetry.atoms.xf .= mod.(non_P1_cartesian_symmetry.atoms.xf, 1.0)
+    non_P1_cartesian_symmetry.charges.xf .= mod.(non_P1_cartesian_symmetry.charges.xf, 1.0)
+
+    # make sure frameworks are now recorded as being in P1 with their is_p1 flag
+    @test non_P1_framework_symmetry.is_p1
+    @test non_P1_cartesian_symmetry.is_p1
+
+    # test that same structure is created when reading and converting to P1 and
+    #   when reading then converting to P1
+    @test isapprox(non_P1_framework, non_P1_framework_symmetry)
+    @test isapprox(non_P1_cartesian, non_P1_cartesian_symmetry)
 
     # test .cssr reader too; test_structure2.{cif,cssr} designed to be the same.
     framework_from_cssr = Framework("test_structure2.cssr")
     strip_numbers_from_atom_labels!(framework_from_cssr)
-    @test isapprox(framework_from_cssr, framework, checknames=false)
+    @test has_same_sets_of_atoms_and_charges(framework_from_cssr, framework, checknames=false)
 
     # test replicate framework
     sbmof = Framework("SBMOF-1.cif")
+    strip_numbers_from_atom_labels!(sbmof)
     replicated_sbmof = replicate(sbmof, (1, 1, 1))
     @test isapprox(sbmof, replicated_sbmof)
+    # test replication no bonds assertion
+    sbmof_bonds = Framework("SBMOF-1.cif")
+    strip_numbers_from_atom_labels!(sbmof_bonds)
+    bonding_rules = [BondingRule(:H, :*, 0.4, 1.2),
+                     BondingRule(:Ca, :O, 0.4, 2.5),
+                     BondingRule(:*, :*, 0.4, 1.9)]
+    infer_bonds!(sbmof_bonds, true, bonding_rules)
+    @test_throws AssertionError replicate(sbmof_bonds, (2, 2, 2))
+    # write out and compare to inferred bonds
+    write_cif(sbmof_bonds, joinpath(pwd(), "data", "crystals", "SBMOF-1_inferred_bonds.cif"))
+    sbmof_inferred_bonds = Framework("SBMOF-1_inferred_bonds.cif"; read_bonds_from_file=true)
+    strip_numbers_from_atom_labels!(sbmof_inferred_bonds)
+    @test compare_bonds_in_framework(sbmof_bonds, sbmof_inferred_bonds)
+    # other bond info tests
+    # TODO find more robust test/confirm these are the correct numbers
+    # replacing this test with the one below comparing pdb bond info to inferred
+    #   bond info
+    sbmof_bonds_copy = Framework("SBMOF-1.cif")
+    strip_numbers_from_atom_labels!(sbmof_bonds_copy)
+    # reverse the order of the atoms and bond info should still be the same
+    sbmof_bonds_copy.atoms.xf .= reverse(sbmof_bonds_copy.atoms.xf; dims=2)
+    sbmof_bonds_copy.atoms.species .= reverse(sbmof_bonds_copy.atoms.species)
+    infer_bonds!(sbmof_bonds_copy, true, bonding_rules)
+    @test compare_bonds_in_framework(sbmof_bonds, sbmof_bonds_copy)
+    remove_bonds!(sbmof_bonds)
+    @test ne(sbmof_bonds.bonds) == 0
+    @test !compare_bonds_in_framework(sbmof_bonds, sbmof_bonds_copy)
+
+    # test reading in bonds as part of `Framework()`
+    sbmof_read_bonds = Framework("test_bond_viz.cif"; read_bonds_from_file=true, check_atom_and_charge_overlap=false)
+    strip_numbers_from_atom_labels!(sbmof_read_bonds)
+    @test ne(sbmof_read_bonds.bonds) == 5
+    write_cif(sbmof_read_bonds, joinpath(pwd(), "data", "crystals", "rewritten_sbmof_read_bonds.cif"))
+    reloaded_sbmof_read_bonds = Framework("rewritten_sbmof_read_bonds.cif"; read_bonds_from_file=true, check_atom_and_charge_overlap=false)
+    strip_numbers_from_atom_labels!(reloaded_sbmof_read_bonds)
+    @test compare_bonds_in_framework(sbmof_read_bonds, reloaded_sbmof_read_bonds)
+
+    # Test that reading in bonding information is the same as inferring the
+    #   bonding info
+    # Bonding information is from a pdb file saved by avogadro
+    # using non-p1 because it meant copying over fewer bonds
+    read_bonds = Framework("KAXQIL_clean_cartn.cif"; convert_to_p1=false, read_bonds_from_file=true)
+    strip_numbers_from_atom_labels!(read_bonds)
+    inferred_bonds = Framework("KAXQIL_clean.cif"; convert_to_p1=false)
+    strip_numbers_from_atom_labels!(inferred_bonds)
+    # Using same bonding rules as above
+    infer_bonds!(inferred_bonds, true, bonding_rules)
+    @test compare_bonds_in_framework(read_bonds, inferred_bonds; atol=1e-6)
 
     repfactors = replication_factors(sbmof.box, 14.0)
     replicated_sbmof = replicate(sbmof, repfactors)
@@ -51,6 +185,58 @@ using Random
     @test chemical_formula(sbmof) == chemical_formula(replicated_sbmof)
     @test isapprox(crystal_density(sbmof), crystal_density(replicated_sbmof), atol=1e-7)
 
+    # test symmetry_rules
+    # define default symmetry_rules
+    symmetry_rules = [Array{AbstractString, 2}(undef, 3, 0) ["x", "y", "z"]]
+    other_symmetry_rules = [Array{AbstractString, 2}(undef, 3, 0) ["y + z", "x + z", "x + y"]]
+    symmetry_rules_two = [Array{AbstractString, 2}(undef, 3, 0) ["x" "y + z";
+                                                                 "y" "x + z";
+                                                                 "z" "x + y"]]
+    symmetry_rules_two_cpy = deepcopy(symmetry_rules_two)
+    @test ! is_symmetry_equal(symmetry_rules, symmetry_rules_two)
+    @test ! is_symmetry_equal(symmetry_rules, other_symmetry_rules)
+    @test is_symmetry_equal(symmetry_rules, symmetry_rules)
+    @test is_symmetry_equal(symmetry_rules_two, symmetry_rules_two_cpy)
+
+    # test framework addition
+    f1 = Framework("framework 1", UnitCube(), Atoms([:a, :b],
+                                                    [1.0 4.0;
+                                                     2.0 5.0;
+                                                     3.0 6.0]),
+                                              Charges([0.1, 0.2],
+                                                      [1.0 4.0;
+                                                       2.0 5.0;
+                                                       3.0 6.0]))
+    f2 = Framework("framework 2", UnitCube(), Atoms([:c, :d],
+                                                    [7.0 10.0;
+                                                     8.0 11.0;
+                                                     9.0 12.0]),
+                                              Charges([0.3, 0.4],
+                                                      [7.0 10.0;
+                                                       8.0 11.0;
+                                                       9.0 12.0]))
+    f3 = f1 + f2
+    addition_bonding_rules = [BondingRule(:a, :b, 4.5, 5.3),
+                              BondingRule(:c, :d, 4.5, 5.3)]
+    @test is_bonded(f1, 1, 2, [BondingRule(:a, :b, 1.0, 5.5)]; include_bonds_across_periodic_boundaries=false)
+    @test ! is_bonded(f2, 1, 2, [BondingRule(:c, :d, 1.0, 4.5)]; include_bonds_across_periodic_boundaries=false)
+    infer_bonds!(f1, false, addition_bonding_rules)
+    infer_bonds!(f2, false, addition_bonding_rules)
+    @test ! compare_bonds_in_framework(f1 + f2, f3)
+    infer_bonds!(f3, false, addition_bonding_rules)
+    @test compare_bonds_in_framework(f1 + f2, f3)
+    @test_throws AssertionError f1 + sbmof # only allow frameworks with same box
+    @test isapprox(f1.box, f3.box)
+    @test isapprox(f2.box, f3.box)
+    @test isapprox(f1.atoms + f2.atoms, f3.atoms)
+    @test isapprox(f1.charges + f2.charges, f3.charges)
+
+    # test infinite framework_addition
+    f4 = +(f1, f2, f3; check_overlap=false)
+    @test isapprox(f3.box, f4.box)
+    @test isapprox(f4.atoms, f1.atoms + f2.atoms + f3.atoms)
+    @test isapprox(f4.charges, f1.charges + f2.charges + f3.charges)
+
     # more xtal tests
     sbmof1 = Framework("SBMOF-1.cif")
     @test !charged(sbmof1)
@@ -64,6 +250,16 @@ using Random
     @test rbox.Ω ≈ sbmof1.box.Ω * 2 * 3 * 4
     @test all(rbox.c_to_f * sbmof1.box.f_to_c * [1.0, 1.0, 1.0] .≈ [1/2, 1/3, 1/4])
 
+    # write bond information to manually inspect if bonds are in order
+    hkust1 = Framework("FIQCEN_clean.cif")
+    strip_numbers_from_atom_labels!(hkust1)
+    br = default_bondingrules()
+    pushfirst!(br, BondingRule(:Cu, :O, 0.3, 2.4))
+    infer_bonds!(hkust1, true, br)
+    reference_bonds = loadgraph("hkust1_reference_bonds.lgz")
+    @test reference_bonds == hkust1.bonds
+
+
 #=    ## .cssr reader test
     # test replicate framework
     sbmof = Framework("SBMOF-1.cif")
diff --git a/test/data/crystals/FIQCEN_clean.cif b/test/data/crystals/FIQCEN_clean.cif
new file mode 100755
index 000000000..d7d2adad3
--- /dev/null
+++ b/test/data/crystals/FIQCEN_clean.cif
@@ -0,0 +1,180 @@
+data_FIQCEN_clean
+_audit_creation_date              2014-07-02
+_audit_creation_method            'Materials Studio'
+_symmetry_space_group_name_H-M    'P1'
+_symmetry_Int_Tables_number       1
+_symmetry_cell_setting            triclinic
+loop_
+_symmetry_equiv_pos_as_xyz
+  x,y,z
+_cell_length_a                    18.6273
+_cell_length_b                    18.6273
+_cell_length_c                    18.6273
+_cell_angle_alpha                 60.0000
+_cell_angle_beta                  60.0000
+_cell_angle_gamma                 60.0000
+loop_
+_atom_site_label
+_atom_site_type_symbol
+_atom_site_fract_x
+_atom_site_fract_y
+_atom_site_fract_z
+_atom_site_U_iso_or_equiv
+_atom_site_adp_type
+_atom_site_occupancy
+Cu1    Cu    0.57057  -0.00000   0.42943   0.01267  Uiso   1.00
+Cu2    Cu    0.00000   0.57057   0.00000   0.01267  Uiso   1.00
+Cu3    Cu    0.42943   0.00000   0.57057   0.01267  Uiso   1.00
+Cu4    Cu    0.00000   0.42943   0.00000   0.01267  Uiso   1.00
+Cu5    Cu    0.57057   0.42943  -0.00000   0.01267  Uiso   1.00
+Cu6    Cu    0.00000   1.00000   0.42943   0.01267  Uiso   1.00
+Cu7    Cu    0.42943   0.57057  -0.00000   0.01267  Uiso   1.00
+Cu8    Cu    0.00000   0.00000   0.57057   0.01267  Uiso   1.00
+Cu9    Cu    0.57057   1.00000  -0.00000   0.01267  Uiso   1.00
+Cu10   Cu    0.00000   0.42943   0.57057   0.01267  Uiso   1.00
+Cu11   Cu    0.42943   1.00000   0.00000   0.01267  Uiso   1.00
+Cu12   Cu    0.00000   0.57057   0.42943   0.01267  Uiso   1.00
+H1     H     0.72800   0.72800   0.51160   0.01267  Uiso   1.00
+H2     H     0.72800   0.72800   0.03240   0.01267  Uiso   1.00
+H3     H     0.51160   0.03240   0.72800   0.01267  Uiso   1.00
+H4     H     0.03240   0.51160   0.72800   0.01267  Uiso   1.00
+H5     H     0.72800   0.51160   0.03240   0.01267  Uiso   1.00
+H6     H     0.72800   0.03240   0.51160   0.01267  Uiso   1.00
+H7     H     0.51160   0.72800   0.72800   0.01267  Uiso   1.00
+H8     H     0.03240   0.72800   0.72800   0.01267  Uiso   1.00
+H9     H     0.72800   0.03240   0.72800   0.01267  Uiso   1.00
+H10    H     0.72800   0.51160   0.72800   0.01267  Uiso   1.00
+H11    H     0.51160   0.72800   0.03240   0.01267  Uiso   1.00
+H12    H     0.03240   0.72800   0.51160   0.01267  Uiso   1.00
+H13    H     0.27200   0.27200   0.48840   0.01267  Uiso   1.00
+H14    H     0.27200   0.27200   0.96760   0.01267  Uiso   1.00
+H15    H     0.96760   0.48840   0.27200   0.01267  Uiso   1.00
+H16    H     0.48840   0.96760   0.27200   0.01267  Uiso   1.00
+H17    H     0.48840   0.27200   0.96760   0.01267  Uiso   1.00
+H18    H     0.96760   0.27200   0.48840   0.01267  Uiso   1.00
+H19    H     0.27200   0.48840   0.27200   0.01267  Uiso   1.00
+H20    H     0.27200   0.96760   0.27200   0.01267  Uiso   1.00
+H21    H     0.96760   0.27200   0.27200   0.01267  Uiso   1.00
+H22    H     0.48840   0.27200   0.27200   0.01267  Uiso   1.00
+H23    H     0.27200   0.48840   0.96760   0.01267  Uiso   1.00
+H24    H     0.27200   0.96760   0.48840   0.01267  Uiso   1.00
+C1     C     0.25700   0.96900   0.38700   0.01267  Uiso   1.00
+C2     C     0.96900   0.25700   0.38700   0.01267  Uiso   1.00
+C3     C     0.38700   0.38700   0.25700   0.01267  Uiso   1.00
+C4     C     0.38700   0.38700   0.96900   0.01267  Uiso   1.00
+C5     C     0.25700   0.38700   0.38700   0.01267  Uiso   1.00
+C6     C     0.96900   0.38700   0.38700   0.01267  Uiso   1.00
+C7     C     0.38700   0.25700   0.96900   0.01267  Uiso   1.00
+C8     C     0.38700   0.96900   0.25700   0.01267  Uiso   1.00
+C9     C     0.25700   0.38700   0.96900   0.01267  Uiso   1.00
+C10    C     0.96900   0.38700   0.25700   0.01267  Uiso   1.00
+C11    C     0.38700   0.96900   0.38700   0.01267  Uiso   1.00
+C12    C     0.38700   0.25700   0.38700   0.01267  Uiso   1.00
+C13    C     0.03100   0.74300   0.61300   0.01267  Uiso   1.00
+C14    C     0.74300   0.03100   0.61300   0.01267  Uiso   1.00
+C15    C     0.61300   0.61300   0.74300   0.01267  Uiso   1.00
+C16    C     0.61300   0.61300   0.03100   0.01267  Uiso   1.00
+C17    C     0.61300   0.74300   0.61300   0.01267  Uiso   1.00
+C18    C     0.61300   0.03100   0.61300   0.01267  Uiso   1.00
+C19    C     0.74300   0.61300   0.03100   0.01267  Uiso   1.00
+C20    C     0.03100   0.61300   0.74300   0.01267  Uiso   1.00
+C21    C     0.61300   0.74300   0.03100   0.01267  Uiso   1.00
+C22    C     0.61300   0.03100   0.74300   0.01267  Uiso   1.00
+C23    C     0.03100   0.61300   0.61300   0.01267  Uiso   1.00
+C24    C     0.74300   0.61300   0.61300   0.01267  Uiso   1.00
+C25    C     0.56870   0.56870   0.83770   0.01267  Uiso   1.00
+C26    C     0.56870   0.56870   0.02490   0.01267  Uiso   1.00
+C27    C     0.83770   0.02490   0.56870   0.01267  Uiso   1.00
+C28    C     0.02490   0.83770   0.56870   0.01267  Uiso   1.00
+C29    C     0.56870   0.83770   0.02490   0.01267  Uiso   1.00
+C30    C     0.56870   0.02490   0.83770   0.01267  Uiso   1.00
+C31    C     0.83770   0.56870   0.56870   0.01267  Uiso   1.00
+C32    C     0.02490   0.56870   0.56870   0.01267  Uiso   1.00
+C33    C     0.56870   0.02490   0.56870   0.01267  Uiso   1.00
+C34    C     0.56870   0.83770   0.56870   0.01267  Uiso   1.00
+C35    C     0.83770   0.56870   0.02490   0.01267  Uiso   1.00
+C36    C     0.02490   0.56870   0.83770   0.01267  Uiso   1.00
+C37    C     0.43130   0.43130   0.16230   0.01267  Uiso   1.00
+C38    C     0.43130   0.43130   0.97510   0.01267  Uiso   1.00
+C39    C     0.97510   0.16230   0.43130   0.01267  Uiso   1.00
+C40    C     0.16230   0.97510   0.43130   0.01267  Uiso   1.00
+C41    C     0.16230   0.43130   0.97510   0.01267  Uiso   1.00
+C42    C     0.97510   0.43130   0.16230   0.01267  Uiso   1.00
+C43    C     0.43130   0.16230   0.43130   0.01267  Uiso   1.00
+C44    C     0.43130   0.97510   0.43130   0.01267  Uiso   1.00
+C45    C     0.97510   0.43130   0.43130   0.01267  Uiso   1.00
+C46    C     0.16230   0.43130   0.43130   0.01267  Uiso   1.00
+C47    C     0.43130   0.16230   0.97510   0.01267  Uiso   1.00
+C48    C     0.43130   0.97510   0.16230   0.01267  Uiso   1.00
+C49    C     0.30060   0.96840   0.43040   0.01267  Uiso   1.00
+C50    C     0.96840   0.30060   0.30060   0.01267  Uiso   1.00
+C51    C     0.43040   0.30060   0.30060   0.01267  Uiso   1.00
+C52    C     0.30060   0.43040   0.96840   0.01267  Uiso   1.00
+C53    C     0.30060   0.43040   0.30060   0.01267  Uiso   1.00
+C54    C     0.96840   0.30060   0.43040   0.01267  Uiso   1.00
+C55    C     0.43040   0.30060   0.96840   0.01267  Uiso   1.00
+C56    C     0.30060   0.96840   0.30060   0.01267  Uiso   1.00
+C57    C     0.30060   0.30060   0.96840   0.01267  Uiso   1.00
+C58    C     0.96840   0.43040   0.30060   0.01267  Uiso   1.00
+C59    C     0.43040   0.96840   0.30060   0.01267  Uiso   1.00
+C60    C     0.30060   0.30060   0.43040   0.01267  Uiso   1.00
+C61    C     0.03160   0.69940   0.56960   0.01267  Uiso   1.00
+C62    C     0.69940   0.03160   0.69940   0.01267  Uiso   1.00
+C63    C     0.69940   0.56960   0.69940   0.01267  Uiso   1.00
+C64    C     0.56960   0.69940   0.03160   0.01267  Uiso   1.00
+C65    C     0.56960   0.69940   0.69940   0.01267  Uiso   1.00
+C66    C     0.69940   0.03160   0.56960   0.01267  Uiso   1.00
+C67    C     0.69940   0.56960   0.03160   0.01267  Uiso   1.00
+C68    C     0.03160   0.69940   0.69940   0.01267  Uiso   1.00
+C69    C     0.69940   0.69940   0.03160   0.01267  Uiso   1.00
+C70    C     0.56960   0.03160   0.69940   0.01267  Uiso   1.00
+C71    C     0.03160   0.56960   0.69940   0.01267  Uiso   1.00
+C72    C     0.69940   0.69940   0.56960   0.01267  Uiso   1.00
+O1     O     0.61201   0.49247   0.87425   0.01267  Uiso   1.00
+O2     O     0.49247   0.61201   0.02127   0.01267  Uiso   1.00
+O3     O     0.87425   0.02127   0.61201   0.01267  Uiso   1.00
+O4     O     0.02127   0.87425   0.49247   0.01267  Uiso   1.00
+O5     O     0.61201   0.87425   0.02127   0.01267  Uiso   1.00
+O6     O     0.49247   0.02127   0.87425   0.01267  Uiso   1.00
+O7     O     0.87425   0.61201   0.49247   0.01267  Uiso   1.00
+O8     O     0.02127   0.49247   0.61201   0.01267  Uiso   1.00
+O9     O     0.61201   0.02127   0.49247   0.01267  Uiso   1.00
+O10    O     0.49247   0.87425   0.61201   0.01267  Uiso   1.00
+O11    O     0.87425   0.49247   0.02127   0.01267  Uiso   1.00
+O12    O     0.02127   0.61201   0.87425   0.01267  Uiso   1.00
+O13    O     0.50753   0.38799   0.12575   0.01267  Uiso   1.00
+O14    O     0.38799   0.50753   0.97873   0.01267  Uiso   1.00
+O15    O     0.97873   0.12575   0.38799   0.01267  Uiso   1.00
+O16    O     0.12575   0.97873   0.50753   0.01267  Uiso   1.00
+O17    O     0.12575   0.38799   0.97873   0.01267  Uiso   1.00
+O18    O     0.97873   0.50753   0.12575   0.01267  Uiso   1.00
+O19    O     0.38799   0.12575   0.50753   0.01267  Uiso   1.00
+O20    O     0.50753   0.97873   0.38799   0.01267  Uiso   1.00
+O21    O     0.97873   0.38799   0.50753   0.01267  Uiso   1.00
+O22    O     0.12575   0.50753   0.38799   0.01267  Uiso   1.00
+O23    O     0.50753   0.12575   0.97873   0.01267  Uiso   1.00
+O24    O     0.38799   0.97873   0.12575   0.01267  Uiso   1.00
+O25    O     0.38799   0.50753   0.12575   0.01267  Uiso   1.00
+O26    O     0.50753   0.38799   0.97873   0.01267  Uiso   1.00
+O27    O     0.12575   0.97873   0.38799   0.01267  Uiso   1.00
+O28    O     0.97873   0.12575   0.50753   0.01267  Uiso   1.00
+O29    O     0.38799   0.12575   0.97873   0.01267  Uiso   1.00
+O30    O     0.50753   0.97873   0.12575   0.01267  Uiso   1.00
+O31    O     0.12575   0.38799   0.50753   0.01267  Uiso   1.00
+O32    O     0.97873   0.50753   0.38799   0.01267  Uiso   1.00
+O33    O     0.38799   0.97873   0.50753   0.01267  Uiso   1.00
+O34    O     0.50753   0.12575   0.38799   0.01267  Uiso   1.00
+O35    O     0.12575   0.50753   0.97873   0.01267  Uiso   1.00
+O36    O     0.97873   0.38799   0.12575   0.01267  Uiso   1.00
+O37    O     0.49247   0.61201   0.87425   0.01267  Uiso   1.00
+O38    O     0.61201   0.49247   0.02127   0.01267  Uiso   1.00
+O39    O     0.02127   0.87425   0.61201   0.01267  Uiso   1.00
+O40    O     0.87425   0.02127   0.49247   0.01267  Uiso   1.00
+O41    O     0.87425   0.61201   0.02127   0.01267  Uiso   1.00
+O42    O     0.02127   0.49247   0.87425   0.01267  Uiso   1.00
+O43    O     0.61201   0.87425   0.49247   0.01267  Uiso   1.00
+O44    O     0.49247   0.02127   0.61201   0.01267  Uiso   1.00
+O45    O     0.02127   0.61201   0.49247   0.01267  Uiso   1.00
+O46    O     0.87425   0.49247   0.61201   0.01267  Uiso   1.00
+O47    O     0.49247   0.87425   0.02127   0.01267  Uiso   1.00
+O48    O     0.61201   0.02127   0.87425   0.01267  Uiso   1.00
diff --git a/test/data/crystals/FIQCEN_clean_min_charges.cif b/test/data/crystals/FIQCEN_clean_min_charges.cif
deleted file mode 100644
index 0767e85b6..000000000
--- a/test/data/crystals/FIQCEN_clean_min_charges.cif
+++ /dev/null
@@ -1,185 +0,0 @@
-#======================================================================
-
-# CRYSTAL DATA
-
-#----------------------------------------------------------------------
-
-data_VESTA_phase_1
-
-
-_pd_phase_name                         MOF
-_cell_length_a                         18.74300
-_cell_length_b                         18.76500
-_cell_length_c                         18.76200
-_cell_angle_alpha                      59.98000
-_cell_angle_beta                       60.08300
-_cell_angle_gamma                      60.10900
-_symmetry_space_group_name_H-M         "P 1"
-_symmetry_Int_Tables_number            1
-
-loop_
-_symmetry_equiv_pos_as_xyz
-   "x, y, z"
-
-loop_
-   _atom_site_label
-   _atom_site_fract_x
-   _atom_site_fract_y
-   _atom_site_fract_z
- _atom_site_charge
- Cu1 0.565372 0.000174 0.43460  0.936049294872
- Cu1 0.000235 0.565345 0.99996  0.936049294872
- Cu1 0.434599 0.999819 0.56541  0.936049294872
- Cu1 0.999705 0.434699 8.1e-0  0.936049294872
- Cu1 0.565375 0.434645 0.00013  0.936049294872
- Cu1 0.999771 6.2e-05 0.4346  0.936049294872
- Cu1 0.434607 0.565305 0.9998  0.936049294872
- Cu1 0.000199 0.999914 0.56538  0.936049294872
- Cu1 0.565388 9.1e-05 7.8e-0  0.936049294872
- Cu1 0.999913 0.434586 0.565  0.936049294872
- Cu1 0.434519 0.99988 0.99996  0.936049294872
- Cu1 4.4e-05 0.565377 0.43451  0.936049294872
- H1 0.734134 0.734163 0.50424  0.114709294872
- H1 0.734064 0.734487 0.0273  0.114709294872
- H1 0.504265 0.027895 0.73400  0.114709294872
- H1 0.027289 0.504629 0.73401  0.114709294872
- H1 0.734049 0.505152 0.02700  0.114709294872
- H1 0.734177 0.026584 0.50518  0.114709294872
- H1 0.50476 0.733824 0.73365  0.114709294872
- H1 0.028097 0.733942 0.73364  0.114709294872
- H1 0.734049 0.02771 0.73423  0.114709294872
- H1 0.734257 0.504742 0.73405  0.114709294872
- H1 0.504338 0.733979 0.02771  0.114709294872
- H1 0.028018 0.733945 0.50412  0.114709294872
- H1 0.265827 0.26583 0.49577  0.114709294872
- H1 0.265807 0.265538 0.97267  0.114709294872
- H1 0.972646 0.49535 0.26605  0.114709294872
- H1 0.495719 0.972071 0.26597  0.114709294872
- H1 0.495633 0.265874 0.97229  0.114709294872
- H1 0.97194 0.26599 0.49590  0.114709294872
- H1 0.265742 0.495204 0.26591  0.114709294872
- H1 0.265919 0.972219 0.26582  0.114709294872
- H1 0.971869 0.266065 0.2663  0.114709294872
- H1 0.49523 0.266113 0.26639  0.114709294872
- H1 0.266005 0.494784 0.97303  0.114709294872
- H1 0.265833 0.973392 0.49482  0.114709294872
- C1 0.256691 0.971502 0.38557  -0.170020705128
- C1 0.97067 0.256844 0.38632  -0.170020705128
- C1 0.385733 0.385918 0.25691  -0.170020705128
- C1 0.386069 0.385604 0.97132  -0.170020705128
- C1 0.256608 0.385792 0.38610  -0.170020705128
- C1 0.970887 0.385957 0.38624  -0.170020705128
- C1 0.38592 0.256519 0.97127  -0.170020705128
- C1 0.386037 0.970908 0.25670  -0.170020705128
- C1 0.256702 0.385432 0.97152  -0.170020705128
- C1 0.970944 0.385983 0.25697  -0.170020705128
- C1 0.386037 0.971371 0.38567  -0.170020705128
- C1 0.385749 0.256766 0.38629  -0.170020705128
- C1 0.029296 0.743126 0.61368  -0.170020705128
- C1 0.743301 0.028449 0.61445  -0.170020705128
- C1 0.614256 0.614029 0.74308  -0.170020705128
- C1 0.61394 0.614298 0.02870  -0.170020705128
- C1 0.614222 0.743203 0.61372  -0.170020705128
- C1 0.613962 0.028602 0.61431  -0.170020705128
- C1 0.743261 0.614553 0.02851  -0.170020705128
- C1 0.029006 0.614014 0.74305  -0.170020705128
- C1 0.614 0.743424 0.02876  -0.170020705128
- C1 0.613941 0.029042 0.74331  -0.170020705128
- C1 0.029057 0.614 0.61380  -0.170020705128
- C1 0.743368 0.614181 0.61388  -0.170020705128
- C2 0.570015 0.569768 0.83672  0.692409294872
- C2 0.569676 0.569984 0.02361  0.692409294872
- C2 0.836951 0.023115 0.570  0.692409294872
- C2 0.024044 0.836744 0.56946  0.692409294872
- C2 0.569789 0.836929 0.02358  0.692409294872
- C2 0.569723 0.023923 0.83683  0.692409294872
- C2 0.836882 0.56996 0.5696  0.692409294872
- C2 0.023885 0.569762 0.56960  0.692409294872
- C2 0.569727 0.023395 0.57008  0.692409294872
- C2 0.569984 0.83677 0.56948  0.692409294872
- C2 0.836871 0.57037 0.02334  0.692409294872
- C2 0.023722 0.569814 0.83672  0.692409294872
- C2 0.429962 0.430184 0.16327  0.692409294872
- C2 0.430347 0.429912 0.97640  0.692409294872
- C2 0.975946 0.163219 0.43052  0.692409294872
- C2 0.163042 0.976837 0.42973  0.692409294872
- C2 0.163094 0.429662 0.97669  0.692409294872
- C2 0.976236 0.430203 0.163  0.692409294872
- C2 0.429978 0.163203 0.43056  0.692409294872
- C2 0.430275 0.976585 0.42990  0.692409294872
- C2 0.976078 0.430178 0.43045  0.692409294872
- C2 0.163096 0.430024 0.43030  0.692409294872
- C2 0.430112 0.163016 0.97643  0.692409294872
- C2 0.430247 0.976024 0.16318  0.692409294872
- C3 0.29982 0.971615 0.42857  0.0341588782051
- C3 0.970626 0.299953 0.30018  0.0341588782051
- C3 0.428786 0.299926 0.30017  0.0341588782051
- C3 0.299903 0.428483 0.97147  0.0341588782051
- C3 0.299666 0.428844 0.29993  0.0341588782051
- C3 0.970644 0.299906 0.42939  0.0341588782051
- C3 0.429112 0.299661 0.97110  0.0341588782051
- C3 0.299851 0.971008 0.29968  0.0341588782051
- C3 0.299767 0.299475 0.97131  0.0341588782051
- C3 0.97101 0.428978 0.30005  0.0341588782051
- C3 0.429169 0.970912 0.29978  0.0341588782051
- C3 0.299698 0.299764 0.42928  0.0341588782051
- C3 0.029306 0.700049 0.57063  0.0341588782051
- C3 0.700131 0.028932 0.70035  0.0341588782051
- C3 0.700321 0.57111 0.70005  0.0341588782051
- C3 0.570856 0.700243 0.02892  0.0341588782051
- C3 0.5712 0.700028 0.69984  0.0341588782051
- C3 0.700179 0.028352 0.57144  0.0341588782051
- C3 0.700107 0.571457 0.02856  0.0341588782051
- C3 0.029328 0.70004 0.69982  0.0341588782051
- C3 0.70015 0.700508 0.02872  0.0341588782051
- C3 0.570817 0.02905 0.7002  0.0341588782051
- C3 0.028924 0.571002 0.70000  0.0341588782051
- C3 0.700272 0.700216 0.5707  0.0341588782051
- O1 0.613937 0.492317 0.87335  -0.569640705128
- O1 0.492152 0.613822 0.0206  -0.569640705128
- O1 0.873625 0.020109 0.61422  -0.569640705128
- O1 0.020827 0.873315 0.4920  -0.569640705128
- O1 0.61378 0.873494 0.02066  -0.569640705128
- O1 0.492202 0.020866 0.87325  -0.569640705128
- O1 0.873398 0.613903 0.49210  -0.569640705128
- O1 0.020853 0.492246 0.61366  -0.569640705128
- O1 0.613694 0.020561 0.49257  -0.569640705128
- O1 0.492503 0.873086 0.61345  -0.569640705128
- O1 0.873137 0.492907 0.02046  -0.569640705128
- O1 0.020922 0.613731 0.87318  -0.569640705128
- O1 0.507458 0.386317 0.12688  -0.569640705128
- O1 0.386412 0.50739 0.97921  -0.569640705128
- O1 0.978811 0.126777 0.38655  -0.569640705128
- O1 0.126736 0.97978 0.50720  -0.569640705128
- O1 0.126374 0.385808 0.97965  -0.569640705128
- O1 0.979462 0.507647 0.12672  -0.569640705128
- O1 0.386036 0.12656 0.50805  -0.569640705128
- O1 0.507793 0.979634 0.38598  -0.569640705128
- O1 0.978996 0.386216 0.50802  -0.569640705128
- O1 0.126687 0.507544 0.38625  -0.569640705128
- O1 0.50764 0.126662 0.97942  -0.569640705128
- O1 0.386241 0.978913 0.1266  -0.569640705128
- O1 0.386029 0.507632 0.12664  -0.569640705128
- O1 0.507866 0.386077 0.97939  -0.569640705128
- O1 0.126366 0.979858 0.38580  -0.569640705128
- O1 0.979179 0.12664 0.5079  -0.569640705128
- O1 0.386101 0.126452 0.97939  -0.569640705128
- O1 0.507756 0.979114 0.12675  -0.569640705128
- O1 0.126582 0.386094 0.5078  -0.569640705128
- O1 0.979096 0.5077 0.38639  -0.569640705128
- O1 0.386304 0.979437 0.50740  -0.569640705128
- O1 0.507457 0.126882 0.38659  -0.569640705128
- O1 0.126852 0.50713 0.97955  -0.569640705128
- O1 0.979018 0.386304 0.12683  -0.569640705128
- O1 0.492524 0.613647 0.87311  -0.569640705128
- O1 0.613621 0.492521 0.02076  -0.569640705128
- O1 0.02117 0.873189 0.61343  -0.569640705128
- O1 0.873262 0.02016 0.49283  -0.569640705128
- O1 0.873573 0.614239 0.02038  -0.569640705128
- O1 0.020513 0.492364 0.87330  -0.569640705128
- O1 0.613928 0.873424 0.49199  -0.569640705128
- O1 0.492215 0.02033 0.61401  -0.569640705128
- O1 0.020989 0.613713 0.49202  -0.569640705128
- O1 0.873291 0.492436 0.61372  -0.569640705128
- O1 0.492277 0.873288 0.02055  -0.569640705128
- O1 0.613735 0.020997 0.87337  -0.569640705128
diff --git a/test/data/crystals/KAXQIL_clean.cif b/test/data/crystals/KAXQIL_clean.cif
new file mode 100644
index 000000000..e539e59fa
--- /dev/null
+++ b/test/data/crystals/KAXQIL_clean.cif
@@ -0,0 +1,72 @@
+
+#######################################################################
+#
+#                 Cambridge Crystallographic Data Centre
+#                                CCDC 
+#
+#######################################################################
+#
+# If this CIF has been generated from an entry in the Cambridge 
+# Structural Database, then it will include bibliographic, chemical, 
+# crystal, experimental, refinement or atomic coordinate data resulting 
+# from the CCDC's data processing and validation procedures.
+#
+#######################################################################
+
+data_casdb
+_symmetry_cell_setting           monoclinic
+_symmetry_space_group_name_H-M   'P 21/n'
+_symmetry_Int_Tables_number      14
+_space_group_name_Hall           '-P 2yn'
+loop_
+_symmetry_equiv_pos_site_id
+_symmetry_equiv_pos_as_xyz
+1 x,y,z
+2 1/2-x,1/2+y,1/2-z
+3 -x,-y,-z
+4 1/2+x,1/2-y,1/2+z
+_cell_length_a                   11.8214(3)
+_cell_length_b                   5.56730(13)
+_cell_length_c                   22.7603(5)
+_cell_angle_alpha                90.00
+_cell_angle_beta                 100.356(2)
+_cell_angle_gamma                90.00
+_cell_volume                     1473.53
+loop_
+_atom_site_label
+_atom_site_type_symbol
+_atom_site_fract_x
+_atom_site_fract_y
+_atom_site_fract_z
+Ca1 Ca 0.36318(3) 0.24415(5) 0.277126(14)
+S1 S 0.60165(3) 0.03637(7) 0.392320(16)
+O6 O 0.50690(10) 0.1950(2) 0.36904(5)
+O1 O 0.73774(11) 0.0848(2) 0.69064(5)
+O5 O 0.58466(12) -0.2189(2) 0.38824(6)
+O3 O 1.01577(12) 0.4743(2) 0.25622(6)
+O4 O 1.04305(12) 0.0809(2) 0.24715(7)
+C4 C 0.99546(14) 0.2586(3) 0.26511(8)
+C4A C 0.72326(14) 0.2607(3) 0.65507(7)
+C6 C 0.64494(14) 0.1092(3) 0.46864(7)
+C7 C 0.71894(14) 0.1104(3) 0.35781(7)
+C8 C 0.80704(16) -0.0567(3) 0.36004(8)
+H8 H 0.8055 -0.1996 0.3811
+C9 C 0.72975(17) -0.0114(3) 0.56772(8)
+H9 H 0.7682 -0.1249 0.5939
+C10 C 0.69872(14) 0.2073(3) 0.58931(7)
+C11 C 0.61656(16) 0.3313(3) 0.48920(8)
+H11 H 0.5794 0.4459 0.4628
+C12 C 0.72051(15) 0.3272(3) 0.32810(8)
+H12 H 0.6612 0.4377 0.3271
+C13 C 0.81276(15) 0.3756(3) 0.29975(8)
+H13 H 0.8161 0.5217 0.2803
+C14 C 0.89738(16) -0.0071(3) 0.33040(9)
+H14 H 0.9561 -0.1186 0.3309
+C15 C 0.70345(18) -0.0610(3) 0.50695(8)
+H15 H 0.7249 -0.2066 0.4922
+C16 C 0.90018(14) 0.2083(3) 0.30006(7)
+C17 C 0.64454(16) 0.3796(3) 0.54979(8)
+H17 H 0.6269 0.5287 0.5642
+O2 O 0.73242(11) 0.4749(2) 0.67168(6)
+
+#END
diff --git a/test/data/crystals/KAXQIL_clean_cartn.cif b/test/data/crystals/KAXQIL_clean_cartn.cif
new file mode 100644
index 000000000..b05f05bc4
--- /dev/null
+++ b/test/data/crystals/KAXQIL_clean_cartn.cif
@@ -0,0 +1,89 @@
+data_KAXQIL_clean_PM
+_symmetry_space_group_name_H-M	'P 21/n'
+_cell_length_a	11.821400
+_cell_length_b	5.567300
+_cell_length_c	22.760300
+_cell_angle_alpha	90.000000
+_cell_angle_beta	100.356000
+_cell_angle_gamma	90.000000
+_symmetry_Int_Tables_number 1
+
+loop_
+_symmetry_equiv_pos_as_xyz
+'x,y,z'
+'1/2-x,1/2+y,1/2-z'
+'-x,-y,-z'
+'1/2+x,1/2-y,1/2+z'
+
+loop_
+_atom_site_label
+_atom_site_type_symbol
+_atom_site_Cartn_x
+_atom_site_Cartn_y
+_atom_site_Cartn_z
+_atom_site_charge
+Ca1	Ca	3.159441	1.359256	6.204721	0.000000
+S1	S	5.507177	0.202483	8.783861	0.000000
+O1	O	4.482349	1.085624	8.262632	0.000000
+O2	O	5.895382	0.472107	15.463105	0.000000
+O3	O	5.323025	4.348618	8.692511	0.000000
+O4	O	-0.861895	2.640570	5.736645	0.000000
+O5	O	-0.502297	0.450395	5.533572	0.000000
+C1	C	10.683039	1.439704	5.935688	0.000000
+C2	C	5.869742	1.451395	14.666710	0.000000
+C3	C	5.706664	0.607949	10.492630	0.000000
+C4	C	7.034906	0.614630	8.011198	0.000000
+C5	C	8.067247	5.251634	8.061127	0.000000
+H1	H	7.962876	4.456067	8.532650	0.000000
+C6	C	6.303853	5.503833	12.710984	0.000000
+H2	H	6.651271	4.871944	13.297142	0.000000
+C7	C	5.848700	1.154101	13.194374	0.000000
+C8	C	5.287052	1.844446	10.952958	0.000000
+H3	H	4.955784	2.482459	10.361875	0.000000
+C9	C	7.175023	1.821621	7.346005	0.000000
+H4	H	6.477987	2.436807	7.323615	0.000000
+C10	C	8.381541	2.091078	6.711262	0.000000
+H5	H	8.500603	2.904460	6.275785	0.000000
+C11	C	9.256464	5.527772	7.397501	0.000000
+H6	H	9.948571	4.907018	7.408696	0.000000
+C12	C	6.241589	5.227695	11.350372	0.000000
+H7	H	6.555508	4.417096	11.020127	0.000000
+C13	C	9.413699	1.159669	6.718202	0.000000
+C14	C	5.369912	2.113347	12.309540	0.000000
+H8	H	5.102424	2.943432	12.632173	0.000000
+O6	O	5.910066	2.643911	15.038600	0.000000
+
+loop_
+_geom_bond_atom_site_label_1
+_geom_bond_atom_site_label_2
+Ca1     O1
+S1      C4
+S1      O1
+S1      O3
+S1      C3
+O2      C2
+O4      C1
+O5      C1
+C1      C13
+C2      C7
+C2      O6
+C3      C8
+C3      C12
+C4      C9
+C4      C5
+C5      C11
+C5      H1
+C6      C12
+C6      C7
+C6      H2
+C7      C14
+C8      H3
+C8      C14
+C9      C10
+C9      H4
+C10     H5
+C10     C13
+C11     C13
+C11     H6
+C12     H7
+C14     H8
diff --git a/test/data/crystals/no_atoms.cif b/test/data/crystals/no_atoms.cif
new file mode 100644
index 000000000..e98b208ee
--- /dev/null
+++ b/test/data/crystals/no_atoms.cif
@@ -0,0 +1,36 @@
+
+#######################################################################
+#
+#                 Cambridge Crystallographic Data Centre
+#                                CCDC 
+#
+#######################################################################
+#
+# If this CIF has been generated from an entry in the Cambridge 
+# Structural Database, then it will include bibliographic, chemical, 
+# crystal, experimental, refinement or atomic coordinate data resulting 
+# from the CCDC's data processing and validation procedures.
+#
+#######################################################################
+
+data_casdb
+_symmetry_cell_setting           monoclinic
+_symmetry_space_group_name_H-M   'P 21/n'
+_symmetry_Int_Tables_number      14
+_space_group_name_Hall           '-P 2yn'
+loop_
+_symmetry_equiv_pos_site_id
+_symmetry_equiv_pos_as_xyz
+1 x,y,z
+2 1/2-x,1/2+y,1/2-z
+3 -x,-y,-z
+4 1/2+x,1/2-y,1/2+z
+_cell_length_a                   11.8214(3)
+_cell_length_b                   5.56730(13)
+_cell_length_c                   22.7603(5)
+_cell_angle_alpha                90.00
+_cell_angle_beta                 100.356(2)
+_cell_angle_gamma                90.00
+_cell_volume                     1473.53
+
+#END
diff --git a/test/data/crystals/non_P1_no_symmetry.cif b/test/data/crystals/non_P1_no_symmetry.cif
new file mode 100644
index 000000000..17f5aa09f
--- /dev/null
+++ b/test/data/crystals/non_P1_no_symmetry.cif
@@ -0,0 +1,66 @@
+
+#######################################################################
+#
+#                 Cambridge Crystallographic Data Centre
+#                                CCDC 
+#
+#######################################################################
+#
+# If this CIF has been generated from an entry in the Cambridge 
+# Structural Database, then it will include bibliographic, chemical, 
+# crystal, experimental, refinement or atomic coordinate data resulting 
+# from the CCDC's data processing and validation procedures.
+#
+#######################################################################
+
+data_casdb
+_symmetry_cell_setting           monoclinic
+_symmetry_space_group_name_H-M   'P 21/n'
+_symmetry_Int_Tables_number      14
+_space_group_name_Hall           '-P 2yn'
+
+_cell_length_a                   11.8214(3)
+_cell_length_b                   5.56730(13)
+_cell_length_c                   22.7603(5)
+_cell_angle_alpha                90.00
+_cell_angle_beta                 100.356(2)
+_cell_angle_gamma                90.00
+_cell_volume                     1473.53
+loop_
+_atom_site_label
+_atom_site_type_symbol
+_atom_site_fract_x
+_atom_site_fract_y
+_atom_site_fract_z
+Ca1 Ca 0.36318(3) 0.24415(5) 0.277126(14)
+S1 S 0.60165(3) 0.03637(7) 0.392320(16)
+O6 O 0.50690(10) 0.1950(2) 0.36904(5)
+O1 O 0.73774(11) 0.0848(2) 0.69064(5)
+O5 O 0.58466(12) -0.2189(2) 0.38824(6)
+O3 O 1.01577(12) 0.4743(2) 0.25622(6)
+O4 O 1.04305(12) 0.0809(2) 0.24715(7)
+C4 C 0.99546(14) 0.2586(3) 0.26511(8)
+C4A C 0.72326(14) 0.2607(3) 0.65507(7)
+C6 C 0.64494(14) 0.1092(3) 0.46864(7)
+C7 C 0.71894(14) 0.1104(3) 0.35781(7)
+C8 C 0.80704(16) -0.0567(3) 0.36004(8)
+H8 H 0.8055 -0.1996 0.3811
+C9 C 0.72975(17) -0.0114(3) 0.56772(8)
+H9 H 0.7682 -0.1249 0.5939
+C10 C 0.69872(14) 0.2073(3) 0.58931(7)
+C11 C 0.61656(16) 0.3313(3) 0.48920(8)
+H11 H 0.5794 0.4459 0.4628
+C12 C 0.72051(15) 0.3272(3) 0.32810(8)
+H12 H 0.6612 0.4377 0.3271
+C13 C 0.81276(15) 0.3756(3) 0.29975(8)
+H13 H 0.8161 0.5217 0.2803
+C14 C 0.89738(16) -0.0071(3) 0.33040(9)
+H14 H 0.9561 -0.1186 0.3309
+C15 C 0.70345(18) -0.0610(3) 0.50695(8)
+H15 H 0.7249 -0.2066 0.4922
+C16 C 0.90018(14) 0.2083(3) 0.30006(7)
+C17 C 0.64454(16) 0.3796(3) 0.54979(8)
+H17 H 0.6269 0.5287 0.5642
+O2 O 0.73242(11) 0.4749(2) 0.67168(6)
+
+#END
diff --git a/test/data/crystals/symmetry_test_structure.cif b/test/data/crystals/symmetry_test_structure.cif
new file mode 100644
index 000000000..abd2caf8a
--- /dev/null
+++ b/test/data/crystals/symmetry_test_structure.cif
@@ -0,0 +1,393 @@
+data_2P0_publ
+
+_pd_block_id
+  2011-01-13T16:00|FE-B||Overall
+
+_audit_creation_method  "from EXP file using GSAS2CIF"
+_audit_creation_date                   2011-01-13T16:00
+_audit_author_name                     ""
+_audit_update_record
+; 2011-01-13T16:00  Initial CIF as created by GSAS2CIF
+;
+
+#=============================================================================
+# this information describes the project, paper etc. for the CIF             #
+# Acta Cryst. Section C papers and editorial correspondence is generated     #
+# from the information in this section                                       #
+#                                                                            #
+#   (from)   CIF submission form for Rietveld refinements (Acta Cryst. C)    #
+#                                                 Version 14 December 1998   #
+#=============================================================================
+# 1. SUBMISSION DETAILS
+
+_publ_contact_author_name            ?   # Name of author for correspondence
+_publ_contact_author_address             # Address of author for correspondence
+; ?
+;
+_publ_contact_author_email           ?
+_publ_contact_author_fax             ?
+_publ_contact_author_phone           ?
+
+_publ_contact_letter
+; ?
+;
+
+_publ_requested_journal              ?
+_publ_requested_coeditor_name        ?
+_publ_requested_category             ?   # Acta C: one of CI/CM/CO/FI/FM/FO
+
+#==============================================================================
+
+# 2. PROCESSING SUMMARY (IUCr Office Use Only)
+
+_journal_data_validation_number      ?
+
+_journal_date_recd_electronic        ?
+_journal_date_to_coeditor            ?
+_journal_date_from_coeditor          ?
+_journal_date_accepted               ?
+_journal_date_printers_first         ?
+_journal_date_printers_final         ?
+_journal_date_proofs_out             ?
+_journal_date_proofs_in              ?
+_journal_coeditor_name               ?
+_journal_coeditor_code               ?
+_journal_coeditor_notes
+; ?
+;
+_journal_techeditor_code             ?
+_journal_techeditor_notes
+; ?
+;
+_journal_coden_ASTM                  ?
+_journal_name_full                   ?
+_journal_year                        ?
+_journal_volume                      ?
+_journal_issue                       ?
+_journal_page_first                  ?
+_journal_page_last                   ?
+_journal_paper_category              ?
+_journal_suppl_publ_number           ?
+_journal_suppl_publ_pages            ?
+
+#==============================================================================
+
+# 3. TITLE AND AUTHOR LIST
+
+_publ_section_title
+; ?
+;
+_publ_section_title_footnote
+; ?
+;
+
+# The loop structure below should contain the names and addresses of all
+# authors, in the required order of publication. Repeat as necessary.
+
+loop_
+	_publ_author_name
+        _publ_author_footnote
+	_publ_author_address
+ ?                                   #<--'Last name, first name'
+; ?
+;
+; ?
+;
+
+#==============================================================================
+
+# 4. TEXT
+
+_publ_section_synopsis
+;  ?
+;
+_publ_section_abstract
+; ?
+;
+_publ_section_comment
+; ?
+;
+_publ_section_exptl_prep      # Details of the preparation of the sample(s)
+                              # should be given here.
+; ?
+;
+_publ_section_exptl_refinement
+; ?
+;
+_publ_section_references
+; ?
+;
+_publ_section_figure_captions
+; ?
+;
+_publ_section_acknowledgements
+; ?
+;
+
+#=============================================================================
+# 5. OVERALL REFINEMENT & COMPUTING DETAILS
+
+_refine_special_details
+; ?
+;
+_pd_proc_ls_special_details
+; ?
+;
+
+# The following items are used to identify the programs used.
+_computing_molecular_graphics     ?
+_computing_publication_material   ?
+
+_refine_ls_weighting_scheme       ?
+_refine_ls_weighting_details      ?
+_refine_ls_hydrogen_treatment     ?
+_refine_ls_extinction_method      ?
+_refine_ls_extinction_coef        ?
+_refine_ls_number_constraints     ?
+
+_refine_ls_restrained_S_all       ?
+_refine_ls_restrained_S_obs       ?
+
+#==============================================================================
+# 6. SAMPLE PREPARATION DATA
+
+# (In the unusual case where multiple samples are used in a single
+#  Rietveld study, this information should be moved into the phase
+#  blocks)
+
+# The following three fields describe the preparation of the material.
+# The cooling rate is in K/min.  The pressure at which the sample was
+# prepared is in kPa.  The temperature of preparation is in K.
+
+_pd_prep_cool_rate                ?
+_pd_prep_pressure                 ?
+_pd_prep_temperature              ?
+
+_pd_char_colour                   ?       # use ICDD colour descriptions
+data_FE-B_overall
+
+_refine_ls_shift/su_max                0.27
+_refine_ls_shift/su_mean               0.11
+_computing_structure_refinement        GSAS
+_refine_ls_number_parameters           45
+_refine_ls_goodness_of_fit_all         0.98
+_refine_ls_number_restraints           2
+_refine_ls_matrix_type                 full
+
+# pointers to the phase blocks
+loop_   _pd_phase_block_id
+      2011-01-13T16:00|FE-B_phase1|||
+# pointers to the diffraction patterns
+loop_   _pd_block_diffractogram_id
+      ?
+
+# Information for phase 1
+data_FE-B_phase_1
+
+_pd_block_id
+  2011-01-13T16:00|FE-B_phase1|||
+
+#==============================================================================
+# 7. CHEMICAL, STRUCTURAL AND CRYSTAL DATA
+
+_pd_char_particle_morphology      ?
+
+_chemical_name_systematic
+; ?
+;
+_chemical_name_common             ?
+_chemical_formula_moiety          ?
+_chemical_formula_structural      ?
+_chemical_formula_analytical      ?
+_chemical_melting_point           ?
+_chemical_compound_source         ?       # for minerals and
+                                          # natural products
+_symmetry_space_group_name_Hall   ?
+
+_exptl_crystal_F_000               ?
+_exptl_crystal_density_diffrn      ?
+_exptl_crystal_density_meas        ?
+_exptl_crystal_density_method      ?
+
+_cell_measurement_temperature     ?
+
+_cell_special_details
+; ?
+;
+
+_geom_special_details             ?
+
+# The following item identifies the program(s) used (if appropriate).
+_computing_structure_solution     ?
+
+#==============================================================================
+
+# 8. Phase information from GSAS
+
+_pd_phase_name
+                   "from /data/people/craigy/MgMOF/BT1/1108/MOF74.xtl"
+_cell_length_a                         25.5177(9)
+_cell_length_b                         25.5177
+_cell_length_c                         6.9661(4)
+_cell_angle_alpha                      90.0
+_cell_angle_beta                       90.0
+_cell_angle_gamma                      120.0
+_cell_volume                           3928.31(24)
+_symmetry_cell_setting                 trigonal
+_symmetry_space_group_name_H-M         "R -3"
+loop_
+_symmetry_equiv_pos_site_id
+_symmetry_equiv_pos_as_xyz
+       1 +x,+y,+z 
+       2 -y,x-y,+z 
+       3 y-x,-x,+z 
+      -1 -x,-y,-z 
+      -2 +y,y-x,-z 
+      -3 x-y,+x,-z 
+     101 +x+1/3,+y+2/3,+z+2/3 
+     102 -y+1/3,x-y+2/3,+z+2/3 
+     103 y-x+1/3,-x+2/3,+z+2/3 
+    -101 -x+2/3,-y+1/3,-z+1/3 
+    -102 +y+2/3,y-x+1/3,-z+1/3 
+    -103 x-y+2/3,+x+1/3,-z+1/3 
+     201 +x+2/3,+y+1/3,+z+1/3 
+     202 -y+2/3,x-y+1/3,+z+1/3 
+     203 y-x+2/3,-x+1/3,+z+1/3 
+    -201 -x+1/3,-y+2/3,-z+2/3 
+    -202 +y+1/3,y-x+2/3,-z+2/3 
+    -203 x-y+1/3,+x+2/3,-z+2/3 
+
+# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS
+
+
+loop_
+      _atom_site_type_symbol
+      _atom_site_label
+      _atom_site_fract_x
+      _atom_site_fract_y
+      _atom_site_fract_z
+      _atom_site_occupancy
+      _atom_site_thermal_displace_type
+      _atom_site_U_iso_or_equiv
+      _atom_site_symmetry_multiplicity
+Fe Fe      0.38873(29)  0.35262(31)  0.1477(10)   1.0        Uiso   0.0093       18
+O O1      0.3182(6)    0.2951(5)    0.3500(16)   1.0        Uiso   0.01896      18
+O O2      0.3068(5)    0.2318(6)    0.5894(17)   1.0        Uiso   0.0267       18
+O O3      0.3545(5)    0.2758(5)    0.0009(16)   1.0        Uiso   0.00934      18
+C C1      0.3189(5)    0.2489(5)    0.4144(13)   1.0        Uiso   0.00914      18
+C C2      0.3276(5)    0.2071(5)    0.2780(16)   1.0        Uiso   0.00989      18
+C C3      0.3444(5)    0.2224(5)    0.0839(15)   1.0        Uiso   0.0201       18
+C C4      0.3512(4)    0.1806(5)   -0.0277(12)   1.0        Uiso   0.01662      18
+H H       0.3596(9)    0.1861(8)   -0.1755(22)   1.0        Uiso   0.01289      18
+O O1a     0.5591(6)    0.6699(6)    0.6775(20)   0.901(8)   Uiso   0.01063      18
+O O1b     0.5236(6)    0.6317(6)    0.7955(17)   0.901(8)   Uiso   0.01161      18
+O O2a     0.5134(6)    0.6747(7)    0.2470(20)   0.866(9)   Uiso   0.03843      18
+O O2b     0.5295(7)    0.7244(6)    0.3087(26)   0.866(9)   Uiso   0.06115      18
+O O3a     0.059(4)     0.1362(25)   0.967(11)    0.221(9)   Uiso   0.05492      18
+O O3b     0.081(13)    0.183(5)     0.879(28)    0.221(9)   Uiso   0.44935      18
+
+loop_ _atom_type_symbol
+      _atom_type_number_in_cell
+                   Fe 18.0        
+                   O  125.32      
+                   C  72.0        
+                   H  18.0        
+ 
+# If you change Z, be sure to change all 3 of the following
+_chemical_formula_sum                  "C4 H Fe O6.96"
+_chemical_formula_weight               216.29
+_cell_formula_units_Z                  18
+
+# MOLECULAR GEOMETRY
+
+loop_
+      _geom_bond_atom_site_label_1
+      _geom_bond_atom_site_label_2
+      _geom_bond_distance
+      _geom_bond_site_symmetry_1
+      _geom_bond_site_symmetry_2
+      _geom_bond_publ_flag
+  Fe        O1               2.177(13)  .      1_555  N
+  Fe        O1               2.155(12)  .    103_554  N
+  Fe        O2               2.047(14)  .    202_554  N
+  Fe        O3               1.985(12)  .      1_555  N
+  Fe        O3               1.978(12)  .    202_555  N
+  Fe        O1a              2.087(15)  .     -1_666  N
+  Fe        O1b              2.103(14)  .     -1_666  N
+  O1        Fe               2.177(13)  .      1_555  N
+  O1        Fe               2.155(12)  .    202_555  N
+  O1        C1               1.269(13)  .      1_555  N
+  O2        Fe               2.047(14)  .    103_555  N
+  O2        C1               1.280(12)  .      1_555  N
+  O3        Fe               1.985(12)  .      1_555  N
+  O3        Fe               1.978(12)  .    103_554  N
+  O3        C3               1.380(17)  .      1_555  N
+  C1        O1               1.269(13)  .      1_555  N
+  C1        O2               1.280(12)  .      1_555  N
+  C1        C2               1.525(13)  .      1_555  N
+  C2        C1               1.525(13)  .      1_555  N
+  C2        C3               1.413(12)  .      1_555  N
+  C2        C4               1.389(14)  .   -201_444  N
+  C2        H                2.098(19)  .   -201_444  N
+  C3        O3               1.380(17)  .      1_555  N
+  C3        C2               1.413(12)  .      1_555  N
+  C3        C4               1.400(14)  .      1_555  N
+  C4        C2               1.389(14)  .   -201_444  N
+  C4        C3               1.400(14)  .      1_555  N
+  C4        H                1.047(17)  .      1_555  N
+  H         C2               2.098(19)  .   -201_444  N
+  H         C4               1.047(17)  .      1_555  N
+  O1a       Fe               2.087(15)  .     -1_666  N
+  O1a       O1b              1.250(16)  .      1_555  N
+  O1b       Fe               2.103(14)  .     -1_666  N
+  O1b       O1a              1.250(16)  .      1_555  N
+  O2a       O2b             1.1999(10)  .      1_555  N
+  O2b       O2a             1.1999(10)  .      1_555  N
+  O3a       O3b               1.201(4)  .      1_555  N
+  O3b       O3a               1.201(4)  .      1_555  N
+
+loop_
+      _geom_angle_atom_site_label_1
+      _geom_angle_atom_site_label_2
+      _geom_angle_atom_site_label_3
+      _geom_angle
+      _geom_angle_site_symmetry_1
+      _geom_angle_site_symmetry_2
+      _geom_angle_site_symmetry_3
+      _geom_angle_publ_flag
+  O1        Fe        O2                 84.7(5)    103_554  .    202_554  N
+  O1        Fe        O3                 78.2(4)    103_554  .      1_555  N
+  O1        Fe        O3                 89.1(6)    103_554  .    202_555  N
+  O1        Fe        O1a               160.0(6)    103_554  .     -1_666  N
+  O1        Fe        O1b               164.3(5)    103_554  .     -1_666  N
+  O2        Fe        O3                 96.2(5)    202_554  .      1_555  N
+  O2        Fe        O3                100.6(6)    202_554  .    202_555  N
+  O2        Fe        O1a               114.4(6)    202_554  .     -1_666  N
+  O2        Fe        O1b                79.9(5)    202_554  .     -1_666  N
+  O3        Fe        O3                157.9(6)      1_555  .    202_555  N
+  O3        Fe        O1a                93.1(6)      1_555  .     -1_666  N
+  O3        Fe        O1b               100.7(6)      1_555  .     -1_666  N
+  O3        Fe        O1a                92.9(5)    202_555  .     -1_666  N
+  O3        Fe        O1b                96.4(6)    202_555  .     -1_666  N
+  O1a       Fe        O1b                34.7(4)     -1_666  .     -1_666  N
+  Fe        O1        C1               134.2(10)    202_555  .      1_555  N
+  Fe        O2        C1                127.6(9)    103_555  .      1_555  N
+  Fe        O3        Fe                106.1(7)      1_555  .    103_554  N
+  Fe        O3        C3                122.1(9)      1_555  .      1_555  N
+  Fe        O3        C3                121.8(9)    103_554  .      1_555  N
+  O1        C1        O2               122.8(12)      1_555  .      1_555  N
+  O1        C1        C2               120.2(11)      1_555  .      1_555  N
+  O2        C1        C2               116.7(10)      1_555  .      1_555  N
+  C1        C2        C3               122.3(11)      1_555  .      1_555  N
+  C1        C2        C4                113.7(9)      1_555  .   -201_555  N
+  C3        C2        C4               124.0(10)      1_555  .   -201_555  N
+  O3        C3        C2               123.7(11)      1_555  .      1_555  N
+  O3        C3        C4               119.1(10)      1_555  .      1_555  N
+  C2        C3        C4               117.2(11)      1_555  .      1_555  N
+  C2        C4        C3                118.7(8)   -201_555  .      1_555  N
+  C2        C4        H                118.3(13)   -201_555  .      1_555  N
+  C3        C4        H                122.6(14)      1_555  .      1_555  N
+  Fe        O1a       O1b                73.4(9)     -1_666  .      1_555  N
+  Fe        O1b       O1a                71.9(8)     -1_666  .      1_555  N
+#--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--eof--#
+

diff --git a/test/data/crystals/symmetry_test_structure_P1.cif b/test/data/crystals/symmetry_test_structure_P1.cif
new file mode 100644
index 000000000..5d4d9ffe4
--- /dev/null
+++ b/test/data/crystals/symmetry_test_structure_P1.cif
@@ -0,0 +1,291 @@
+# CIF file generated by openbabel 2.4.1, see http://openbabel.sf.net
+data_I
+_chemical_name_common '?'
+_cell_length_a 25.5177
+_cell_length_b 25.5177
+_cell_length_c 6.9661
+_cell_angle_alpha 90
+_cell_angle_beta 90
+_cell_angle_gamma 120
+_space_group_name_H-M_alt 'P 1'
+_space_group_name_Hall 'P 1'
+loop_
+    _symmetry_equiv_pos_as_xyz
+    x,y,z
+loop_
+    _atom_site_label
+    _atom_site_type_symbol
+    _atom_site_fract_x
+    _atom_site_fract_y
+    _atom_site_fract_z
+    _atom_site_occupancy
+    Fe      Fe   0.38873   0.35262   0.14770   1.000
+    O1      O    0.31820   0.29510   0.35000   1.000
+    O2      O    0.30680   0.23180   0.58940   1.000
+    O3      O    0.35450   0.27580   0.00090   1.000
+    C1      C    0.31890   0.24890   0.41440   1.000
+    C2      C    0.32760   0.20710   0.27800   1.000
+    C3      C    0.34440   0.22240   0.08390   1.000
+    C4      C    0.35120   0.18060   0.97230   1.000
+    H       H    0.35960   0.18610   0.82450   1.000
+    O1a     O    0.55910   0.66990   0.67750   0.901
+    O1b     O    0.52360   0.63170   0.79550   0.901
+    O2a     O    0.51340   0.67470   0.24700   0.866
+    O2b     O    0.52950   0.72440   0.30870   0.866
+    O3a     O    0.05900   0.13620   0.96700   0.221
+    O3b     O    0.08100   0.18300   0.87900   0.221
+    Fe      Fe   0.64738   0.03611   0.14770   1.000
+    Fe      Fe   0.96389   0.61127   0.14770   1.000
+    Fe      Fe   0.61127   0.64738   0.85230   1.000
+    Fe      Fe   0.35262   0.96389   0.85230   1.000
+    Fe      Fe   0.03611   0.38873   0.85230   1.000
+    Fe      Fe   0.05540   0.68595   0.48103   1.000
+    Fe      Fe   0.31405   0.36944   0.48103   1.000
+    Fe      Fe   0.63056   0.94460   0.48103   1.000
+    Fe      Fe   0.27794   0.98071   0.18563   1.000
+    Fe      Fe   0.01929   0.29722   0.18563   1.000
+    Fe      Fe   0.70278   0.72206   0.18563   1.000
+    Fe      Fe   0.72206   0.01929   0.81437   1.000
+    Fe      Fe   0.98071   0.70278   0.81437   1.000
+    Fe      Fe   0.29722   0.27794   0.81437   1.000
+    Fe      Fe   0.94460   0.31405   0.51897   1.000
+    Fe      Fe   0.68595   0.63056   0.51897   1.000
+    Fe      Fe   0.36944   0.05540   0.51897   1.000
+    O1      O    0.70490   0.02310   0.35000   1.000
+    O1      O    0.97690   0.68180   0.35000   1.000
+    O1      O    0.68180   0.70490   0.65000   1.000
+    O1      O    0.29510   0.97690   0.65000   1.000
+    O1      O    0.02310   0.31820   0.65000   1.000
+    O1      O    0.98487   0.62843   0.68333   1.000
+    O1      O    0.37157   0.35643   0.68333   1.000
+    O1      O    0.64357   0.01513   0.68333   1.000
+    O1      O    0.34847   0.03823   0.98333   1.000
+    O1      O    0.96177   0.31023   0.98333   1.000
+    O1      O    0.68977   0.65153   0.98333   1.000
+    O1      O    0.65153   0.96177   0.01667   1.000
+    O1      O    0.03823   0.68977   0.01667   1.000
+    O1      O    0.31023   0.34847   0.01667   1.000
+    O1      O    0.01513   0.37157   0.31667   1.000
+    O1      O    0.62843   0.64357   0.31667   1.000
+    O1      O    0.35643   0.98487   0.31667   1.000
+    O2      O    0.76820   0.07500   0.58940   1.000
+    O2      O    0.92500   0.69320   0.58940   1.000
+    O2      O    0.69320   0.76820   0.41060   1.000
+    O2      O    0.23180   0.92500   0.41060   1.000
+    O2      O    0.07500   0.30680   0.41060   1.000
+    O2      O    0.97347   0.56513   0.92273   1.000
+    O2      O    0.43487   0.40833   0.92273   1.000
+    O2      O    0.59167   0.02653   0.92273   1.000
+    O2      O    0.35987   0.10153   0.74393   1.000
+    O2      O    0.89847   0.25833   0.74393   1.000
+    O2      O    0.74167   0.64013   0.74393   1.000
+    O2      O    0.64013   0.89847   0.25607   1.000
+    O2      O    0.10153   0.74167   0.25607   1.000
+    O2      O    0.25833   0.35987   0.25607   1.000
+    O2      O    0.02653   0.43487   0.07727   1.000
+    O2      O    0.56513   0.59167   0.07727   1.000
+    O2      O    0.40833   0.97347   0.07727   1.000
+    O3      O    0.72420   0.07870   0.00090   1.000
+    O3      O    0.92130   0.64550   0.00090   1.000
+    O3      O    0.64550   0.72420   0.99910   1.000
+    O3      O    0.27580   0.92130   0.99910   1.000
+    O3      O    0.07870   0.35450   0.99910   1.000
+    O3      O    0.02117   0.60913   0.33423   1.000
+    O3      O    0.39087   0.41203   0.33423   1.000
+    O3      O    0.58797   0.97883   0.33423   1.000
+    O3      O    0.31217   0.05753   0.33243   1.000
+    O3      O    0.94247   0.25463   0.33243   1.000
+    O3      O    0.74537   0.68783   0.33243   1.000
+    O3      O    0.68783   0.94247   0.66757   1.000
+    O3      O    0.05753   0.74537   0.66757   1.000
+    O3      O    0.25463   0.31217   0.66757   1.000
+    O3      O    0.97883   0.39087   0.66577   1.000
+    O3      O    0.60913   0.58797   0.66577   1.000
+    O3      O    0.41203   0.02117   0.66577   1.000
+    C1      C    0.75110   0.07000   0.41440   1.000
+    C1      C    0.93000   0.68110   0.41440   1.000
+    C1      C    0.68110   0.75110   0.58560   1.000
+    C1      C    0.24890   0.93000   0.58560   1.000
+    C1      C    0.07000   0.31890   0.58560   1.000
+    C1      C    0.98557   0.58223   0.74773   1.000
+    C1      C    0.41777   0.40333   0.74773   1.000
+    C1      C    0.59667   0.01443   0.74773   1.000
+    C1      C    0.34777   0.08443   0.91893   1.000
+    C1      C    0.91557   0.26333   0.91893   1.000
+    C1      C    0.73667   0.65223   0.91893   1.000
+    C1      C    0.65223   0.91557   0.08107   1.000
+    C1      C    0.08443   0.73667   0.08107   1.000
+    C1      C    0.26333   0.34777   0.08107   1.000
+    C1      C    0.01443   0.41777   0.25227   1.000
+    C1      C    0.58223   0.59667   0.25227   1.000
+    C1      C    0.40333   0.98557   0.25227   1.000
+    C2      C    0.79290   0.12050   0.27800   1.000
+    C2      C    0.87950   0.67240   0.27800   1.000
+    C2      C    0.67240   0.79290   0.72200   1.000
+    C2      C    0.20710   0.87950   0.72200   1.000
+    C2      C    0.12050   0.32760   0.72200   1.000
+    C2      C    0.99427   0.54043   0.61133   1.000
+    C2      C    0.45957   0.45383   0.61133   1.000
+    C2      C    0.54617   0.00573   0.61133   1.000
+    C2      C    0.33907   0.12623   0.05533   1.000
+    C2      C    0.87377   0.21283   0.05533   1.000
+    C2      C    0.78717   0.66093   0.05533   1.000
+    C2      C    0.66093   0.87377   0.94467   1.000
+    C2      C    0.12623   0.78717   0.94467   1.000
+    C2      C    0.21283   0.33907   0.94467   1.000
+    C2      C    0.00573   0.45957   0.38867   1.000
+    C2      C    0.54043   0.54617   0.38867   1.000
+    C2      C    0.45383   0.99427   0.38867   1.000
+    C3      C    0.77760   0.12200   0.08390   1.000
+    C3      C    0.87800   0.65560   0.08390   1.000
+    C3      C    0.65560   0.77760   0.91610   1.000
+    C3      C    0.22240   0.87800   0.91610   1.000
+    C3      C    0.12200   0.34440   0.91610   1.000
+    C3      C    0.01107   0.55573   0.41723   1.000
+    C3      C    0.44427   0.45533   0.41723   1.000
+    C3      C    0.54467   0.98893   0.41723   1.000
+    C3      C    0.32227   0.11093   0.24943   1.000
+    C3      C    0.88907   0.21133   0.24943   1.000
+    C3      C    0.78867   0.67773   0.24943   1.000
+    C3      C    0.67773   0.88907   0.75057   1.000
+    C3      C    0.11093   0.78867   0.75057   1.000
+    C3      C    0.21133   0.32227   0.75057   1.000
+    C3      C    0.98893   0.44427   0.58277   1.000
+    C3      C    0.55573   0.54467   0.58277   1.000
+    C3      C    0.45533   0.01107   0.58277   1.000
+    C4      C    0.81940   0.17060   0.97230   1.000
+    C4      C    0.82940   0.64880   0.97230   1.000
+    C4      C    0.64880   0.81940   0.02770   1.000
+    C4      C    0.18060   0.82940   0.02770   1.000
+    C4      C    0.17060   0.35120   0.02770   1.000
+    C4      C    0.01787   0.51393   0.30563   1.000
+    C4      C    0.48607   0.50393   0.30563   1.000
+    C4      C    0.49607   0.98213   0.30563   1.000
+    C4      C    0.31547   0.15273   0.36103   1.000
+    C4      C    0.84727   0.16273   0.36103   1.000
+    C4      C    0.83727   0.68453   0.36103   1.000
+    C4      C    0.68453   0.84727   0.63897   1.000
+    C4      C    0.15273   0.83727   0.63897   1.000
+    C4      C    0.16273   0.31547   0.63897   1.000
+    C4      C    0.98213   0.48607   0.69437   1.000
+    C4      C    0.51393   0.49607   0.69437   1.000
+    C4      C    0.50393   0.01787   0.69437   1.000
+    H       H    0.81390   0.17350   0.82450   1.000
+    H       H    0.82650   0.64040   0.82450   1.000
+    H       H    0.64040   0.81390   0.17550   1.000
+    H       H    0.18610   0.82650   0.17550   1.000
+    H       H    0.17350   0.35960   0.17550   1.000
+    H       H    0.02627   0.51943   0.15783   1.000
+    H       H    0.48057   0.50683   0.15783   1.000
+    H       H    0.49317   0.97373   0.15783   1.000
+    H       H    0.30707   0.14723   0.50883   1.000
+    H       H    0.85277   0.15983   0.50883   1.000
+    H       H    0.84017   0.69293   0.50883   1.000
+    H       H    0.69293   0.85277   0.49117   1.000
+    H       H    0.14723   0.84017   0.49117   1.000
+    H       H    0.15983   0.30707   0.49117   1.000
+    H       H    0.97373   0.48057   0.84217   1.000
+    H       H    0.51943   0.49317   0.84217   1.000
+    H       H    0.50683   0.02627   0.84217   1.000
+    O1a     O    0.33010   0.88920   0.67750   0.901
+    O1a     O    0.11080   0.44090   0.67750   0.901
+    O1a     O    0.44090   0.33010   0.32250   0.901
+    O1a     O    0.66990   0.11080   0.32250   0.901
+    O1a     O    0.88920   0.55910   0.32250   0.901
+    O1a     O    0.22577   0.00323   0.01083   0.901
+    O1a     O    0.99677   0.22253   0.01083   0.901
+    O1a     O    0.77747   0.77423   0.01083   0.901
+    O1a     O    0.10757   0.66343   0.65583   0.901
+    O1a     O    0.33657   0.44413   0.65583   0.901
+    O1a     O    0.55587   0.89243   0.65583   0.901
+    O1a     O    0.89243   0.33657   0.34417   0.901
+    O1a     O    0.66343   0.55587   0.34417   0.901
+    O1a     O    0.44413   0.10757   0.34417   0.901
+    O1a     O    0.77423   0.99677   0.98917   0.901
+    O1a     O    0.00323   0.77747   0.98917   0.901
+    O1a     O    0.22253   0.22577   0.98917   0.901
+    O1b     O    0.36830   0.89190   0.79550   0.901
+    O1b     O    0.10810   0.47640   0.79550   0.901
+    O1b     O    0.47640   0.36830   0.20450   0.901
+    O1b     O    0.63170   0.10810   0.20450   0.901
+    O1b     O    0.89190   0.52360   0.20450   0.901
+    O1b     O    0.19027   0.96503   0.12883   0.901
+    O1b     O    0.03497   0.22523   0.12883   0.901
+    O1b     O    0.77477   0.80973   0.12883   0.901
+    O1b     O    0.14307   0.70163   0.53783   0.901
+    O1b     O    0.29837   0.44143   0.53783   0.901
+    O1b     O    0.55857   0.85693   0.53783   0.901
+    O1b     O    0.85693   0.29837   0.46217   0.901
+    O1b     O    0.70163   0.55857   0.46217   0.901
+    O1b     O    0.44143   0.14307   0.46217   0.901
+    O1b     O    0.80973   0.03497   0.87117   0.901
+    O1b     O    0.96503   0.77477   0.87117   0.901
+    O1b     O    0.22523   0.19027   0.87117   0.901
+    O2a     O    0.32530   0.83870   0.24700   0.866
+    O2a     O    0.16130   0.48660   0.24700   0.866
+    O2a     O    0.48660   0.32530   0.75300   0.866
+    O2a     O    0.67470   0.16130   0.75300   0.866
+    O2a     O    0.83870   0.51340   0.75300   0.866
+    O2a     O    0.18007   0.00803   0.58033   0.866
+    O2a     O    0.99197   0.17203   0.58033   0.866
+    O2a     O    0.82797   0.81993   0.58033   0.866
+    O2a     O    0.15327   0.65863   0.08633   0.866
+    O2a     O    0.34137   0.49463   0.08633   0.866
+    O2a     O    0.50537   0.84673   0.08633   0.866
+    O2a     O    0.84673   0.34137   0.91367   0.866
+    O2a     O    0.65863   0.50537   0.91367   0.866
+    O2a     O    0.49463   0.15327   0.91367   0.866
+    O2a     O    0.81993   0.99197   0.41967   0.866
+    O2a     O    0.00803   0.82797   0.41967   0.866
+    O2a     O    0.17203   0.18007   0.41967   0.866
+    O2b     O    0.27560   0.80510   0.30870   0.866
+    O2b     O    0.19490   0.47050   0.30870   0.866
+    O2b     O    0.47050   0.27560   0.69130   0.866
+    O2b     O    0.72440   0.19490   0.69130   0.866
+    O2b     O    0.80510   0.52950   0.69130   0.866
+    O2b     O    0.19617   0.05773   0.64203   0.866
+    O2b     O    0.94227   0.13843   0.64203   0.866
+    O2b     O    0.86157   0.80383   0.64203   0.866
+    O2b     O    0.13717   0.60893   0.02463   0.866
+    O2b     O    0.39107   0.52823   0.02463   0.866
+    O2b     O    0.47177   0.86283   0.02463   0.866
+    O2b     O    0.86283   0.39107   0.97537   0.866
+    O2b     O    0.60893   0.47177   0.97537   0.866
+    O2b     O    0.52823   0.13717   0.97537   0.866
+    O2b     O    0.80383   0.94227   0.35797   0.866
+    O2b     O    0.05773   0.86157   0.35797   0.866
+    O2b     O    0.13843   0.19617   0.35797   0.866
+    O3a     O    0.86380   0.92280   0.96700   0.221
+    O3a     O    0.07720   0.94100   0.96700   0.221
+    O3a     O    0.94100   0.86380   0.03300   0.221
+    O3a     O    0.13620   0.07720   0.03300   0.221
+    O3a     O    0.92280   0.05900   0.03300   0.221
+    O3a     O    0.72567   0.46953   0.30033   0.221
+    O3a     O    0.53047   0.25613   0.30033   0.221
+    O3a     O    0.74387   0.27433   0.30033   0.221
+    O3a     O    0.60767   0.19713   0.36633   0.221
+    O3a     O    0.80287   0.41053   0.36633   0.221
+    O3a     O    0.58947   0.39233   0.36633   0.221
+    O3a     O    0.39233   0.80287   0.63367   0.221
+    O3a     O    0.19713   0.58947   0.63367   0.221
+    O3a     O    0.41053   0.60767   0.63367   0.221
+    O3a     O    0.27433   0.53047   0.69967   0.221
+    O3a     O    0.46953   0.74387   0.69967   0.221
+    O3a     O    0.25613   0.72567   0.69967   0.221
+    O3b     O    0.81700   0.89800   0.87900   0.221
+    O3b     O    0.10200   0.91900   0.87900   0.221
+    O3b     O    0.91900   0.81700   0.12100   0.221
+    O3b     O    0.18300   0.10200   0.12100   0.221
+    O3b     O    0.89800   0.08100   0.12100   0.221
+    O3b     O    0.74767   0.51633   0.21233   0.221
+    O3b     O    0.48367   0.23133   0.21233   0.221
+    O3b     O    0.76867   0.25233   0.21233   0.221
+    O3b     O    0.58567   0.15033   0.45433   0.221
+    O3b     O    0.84967   0.43533   0.45433   0.221
+    O3b     O    0.56467   0.41433   0.45433   0.221
+    O3b     O    0.41433   0.84967   0.54567   0.221
+    O3b     O    0.15033   0.56467   0.54567   0.221
+    O3b     O    0.43533   0.58567   0.54567   0.221
+    O3b     O    0.25233   0.48367   0.78767   0.221
+    O3b     O    0.51633   0.76867   0.78767   0.221
+    O3b     O    0.23133   0.74767   0.78767   0.221
diff --git a/test/data/crystals/symmetry_test_structure_cartn.cif b/test/data/crystals/symmetry_test_structure_cartn.cif
new file mode 100644
index 000000000..555eabbee
--- /dev/null
+++ b/test/data/crystals/symmetry_test_structure_cartn.cif
@@ -0,0 +1,53 @@
+# CIF file generated by openbabel 2.3.2, see http://openbabel.sf.net
+data_I
+_chemical_name_common ''
+_cell_length_a 25.5177
+_cell_length_b 25.5177
+_cell_length_c 6.9661
+_cell_angle_alpha 90
+_cell_angle_beta 90
+_cell_angle_gamma 120
+_space_group_name_H-M_alt 'R -3'
+_space_group_name_Hall '-R 3'
+_symmetry_space_group_name_H-M 'R -3'
+loop_
+    _symmetry_equiv_pos_as_xyz
+    'x,y,z'
+    '-y,x-y,z'
+    '-x+y,-x,z'
+    '-x,-y,-z'
+    'y,-x+y,-z'
+    'x-y,x,-z'
+    '2/3+x,1/3+y,1/3+z'
+    '2/3-y,1/3+x-y,1/3+z'
+    '2/3-x+y,1/3-x,1/3+z'
+    '2/3-x,1/3-y,1/3-z'
+    '2/3+y,1/3-x+y,1/3-z'
+    '2/3+x-y,1/3+x,1/3-z'
+    '1/3+x,2/3+y,2/3+z'
+    '1/3-y,2/3+x-y,2/3+z'
+    '1/3-x+y,2/3-x,2/3+z'
+    '1/3-x,2/3-y,2/3-z'
+    '1/3+y,2/3-x+y,2/3-z'
+    '1/3+x-y,2/3+x,2/3-z'
+loop_
+    _atom_site_type_symbol
+    _atom_site_label
+    _atom_site_Cartn_x
+    _atom_site_Cartn_y
+    _atom_site_Cartn_z
+     Fe   Fe     5.42047    7.79254    1.02889
+      O   O1     4.35459    6.52141    2.43814
+      O   O2     4.87133    5.12254    4.10582
+      O   O3     5.52713    6.09490    0.00627
+      C   C1     4.96192    5.50043    2.88675
+      C   C2     5.71724    4.57670    1.93658
+      C   C3     5.95073    4.91481    0.58446
+      C   C4     6.65757    3.99108   -0.19296
+      H    H     6.80174    4.11262   -1.22255
+      O  O1a     5.71979   14.80410    4.71953
+      O  O1b     5.30130   13.95992    5.54153
+      O  O2a     4.49239   14.91018    1.72063
+      O  O2b     4.26911   16.00850    2.15044
+      O  O3a    -0.23221    3.00988    6.73622
+      O  O3b    -0.26794    4.04411    6.12320
diff --git a/test/data/crystals/test_bond_viz.cif b/test/data/crystals/test_bond_viz.cif
new file mode 100644
index 000000000..e4ce798aa
--- /dev/null
+++ b/test/data/crystals/test_bond_viz.cif
@@ -0,0 +1,36 @@
+data_MOF_publ
+_cell_length_a  10.0
+_cell_length_b  10.0
+_cell_length_c  10.0
+_cell_angle_alpha   90.0
+_cell_angle_beta    90.0
+_cell_angle_gamma   90.0
+
+_symmetry_space_group_name_H-M  'P1'
+
+loop_
+_symmetry_equiv_pos_as_xyz
+    'x,y,z'
+
+loop_
+_atom_site_type_symbol
+_atom_site_label
+_atom_site_fract_x
+_atom_site_fract_y
+_atom_site_fract_z
+    C   C1  0.1     0.1     0.1
+    H   H1  0.05    0.05    0.05
+    O   O1  0.5     0.5     0.5
+    O   O2  1.0     1.0     1.0
+    O   O3  1.0     1.0     0.0
+    O   O4  1.0     0.0     1.0
+    O   O5  0.0     1.0     1.0
+
+loop_
+_geom_bond_atom_site_label_1
+_geom_bond_atom_site_label_2
+    C1  O1
+    C1  O2
+    C1  O3
+    C1  O4
+    C1  O5
diff --git a/test/henry_checkpoint_test.jl b/test/henry_checkpoint_test.jl
index 6aa35e8c4..f5f7ff215 100644
--- a/test/henry_checkpoint_test.jl
+++ b/test/henry_checkpoint_test.jl
@@ -4,7 +4,7 @@ using FileIO
 using Test
 
 @testset "Henry Checkpoint Tests" begin
-    framework = Framework("FIQCEN_clean_min_charges.cif")
+    framework = Framework("FIQCEN_clean.cif")
     strip_numbers_from_atom_labels!(framework)
     co2 = Molecule("CO2")
     ljff = LJForceField("UFF.csv")
diff --git a/test/hkust1_reference_bonds.lgz b/test/hkust1_reference_bonds.lgz
new file mode 100644
index 000000000..586e8e3a2
--- /dev/null
+++ b/test/hkust1_reference_bonds.lgz
@@ -0,0 +1,193 @@
+156,192,u,graph,2,Int64,simplegraph
+1,117
+1,128
+1,142
+1,151
+2,120
+2,126
+2,143
+2,149
+3,118
+3,127
+3,141
+3,152
+4,119
+4,125
+4,144
+4,150
+5,109
+5,121
+5,134
+5,146
+6,112
+6,123
+6,135
+6,148
+7,110
+7,122
+7,133
+7,145
+8,111
+8,124
+8,136
+8,147
+9,113
+9,131
+9,138
+9,156
+10,116
+10,129
+10,139
+10,154
+11,114
+11,132
+11,137
+11,155
+12,115
+12,130
+12,140
+12,153
+13,108
+14,105
+15,106
+16,107
+17,103
+18,102
+19,101
+20,104
+21,98
+22,99
+23,100
+24,97
+25,96
+26,93
+27,94
+28,95
+29,91
+30,90
+31,89
+32,92
+33,86
+34,87
+35,88
+36,85
+37,76
+37,85
+37,92
+38,75
+38,86
+38,90
+39,73
+39,87
+39,89
+40,74
+40,88
+40,91
+41,82
+41,89
+41,96
+42,81
+42,90
+42,94
+43,83
+43,91
+43,93
+44,84
+44,92
+44,95
+45,77
+45,88
+45,93
+46,78
+46,86
+46,94
+47,80
+47,85
+47,95
+48,79
+48,87
+48,96
+49,64
+49,97
+49,104
+50,63
+50,98
+50,102
+51,61
+51,99
+51,101
+52,62
+52,100
+52,103
+53,70
+53,101
+53,108
+54,69
+54,102
+54,106
+55,71
+55,103
+55,105
+56,72
+56,104
+56,107
+57,65
+57,100
+57,105
+58,66
+58,98
+58,106
+59,68
+59,97
+59,107
+60,67
+60,99
+60,108
+61,109
+61,145
+62,110
+62,146
+63,111
+63,148
+64,112
+64,147
+65,113
+65,155
+66,114
+66,156
+67,115
+67,154
+68,116
+68,153
+69,117
+69,152
+70,118
+70,151
+71,119
+71,149
+72,120
+72,150
+73,121
+73,133
+74,122
+74,134
+75,123
+75,136
+76,124
+76,135
+77,125
+77,143
+78,126
+78,144
+79,127
+79,142
+80,128
+80,141
+81,129
+81,140
+82,130
+82,139
+83,131
+83,137
+84,132
+84,138
diff --git a/test/matter_test.jl b/test/matter_test.jl
new file mode 100644
index 000000000..4f11544bc
--- /dev/null
+++ b/test/matter_test.jl
@@ -0,0 +1,55 @@
+module Matter_Test
+
+using PorousMaterials
+using Random
+using Test
+
+@testset "Matter Tests" begin
+    # testing addition for atoms type
+    a1 = Atoms([:a, :b], [1.0 4.0;
+                          2.0 5.0;
+                          3.0 6.0])
+    a2 = Atoms([:c, :d], [7.0 10.0;
+                          8.0 11.0;
+                          9.0 12.0])
+    a3 = a1 + a2
+    # total atoms in addition is number of atoms combined
+    @test a3.n_atoms == a1.n_atoms + a2.n_atoms
+    # atoms from a1 & a2 should be within the atoms for a3
+    @test issubset(Set(a1.xf), Set(a3.xf))
+    @test issubset(Set(a2.xf), Set(a3.xf))
+    @test issubset(Set(a1.species), Set(a3.species))
+    @test issubset(Set(a2.species), Set(a3.species))
+    # no atoms that weren't in a1 or a2 are in a3
+    @test issetequal(setdiff(Set(a3.xf), union(Set(a1.xf), Set(a2.xf))), Set())
+    @test issetequal(setdiff(Set(a3.species), union(Set(a1.species), Set(a2.species))), Set())
+    a1_mismatch = Atoms([:b, :a], [1.0 4.0;
+                                   2.0 5.0;
+                                   3.0 6.0])
+    @test ! has_same_set_of_atoms(a1, a1_mismatch)
+
+    # testing addition for charges type
+    c1 = Charges([0.1, 0.2], [1.0 4.0;
+                              2.0 5.0;
+                              3.0 6.0])
+    c2 = Charges([0.3, 0.4], [7.0 10.0;
+                              8.0 11.0;
+                              9.0 12.0])
+    c3 = c1 + c2
+    # total charges in c3 must be the number of charges in c1 and c2
+    @test c3.n_charges == c1.n_charges + c2.n_charges
+    # charges from c1 & c2 should be within c3
+    @test issubset(Set(c1.xf), Set(c3.xf))
+    @test issubset(Set(c2.xf), Set(c3.xf))
+    @test issubset(Set(c1.q), Set(c3.q))
+    @test issubset(Set(c2.q), Set(c3.q))
+    # no charges are added to c3 that weren't within c1 or c2
+    @test issetequal(setdiff(Set(c3.xf), union(Set(c1.xf), Set(c2.xf))), Set())
+    @test issetequal(setdiff(Set(c3.q), union(Set(c1.q), Set(c2.q))), Set())
+    c1_mismatch = Charges([0.2, 0.1], [1.0 4.0;
+                                       2.0 5.0;
+                                       3.0 6.0])
+    @test ! has_same_set_of_charges(c1, c1_mismatch)
+
+end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index be7cff947..2233e2176 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@
 # Start Test Script
 
 include("box_test.jl")
+include("matter_test.jl")
 include("crystal_test.jl")
 include("molecule_test.jl")
 include("nearest_image_test.jl")
@@ -12,6 +13,7 @@ include("electrostatics_energetics_test.jl")
 include("mchelpers_test.jl")
 include("guest_guest_energetics_test.jl")
 include("eos_test.jl")
+include("simulation_rules_test.jl")
 include("gcmc_checkpoints_test.jl")
 include("henry_checkpoint_test.jl")
 include("grid_test.jl")
diff --git a/test/simulation_rules_test.jl b/test/simulation_rules_test.jl
new file mode 100644
index 000000000..a2d7ee10c
--- /dev/null
+++ b/test/simulation_rules_test.jl
@@ -0,0 +1,44 @@
+module Simulation_Rules_Test
+
+using PorousMaterials
+using Test
+
+# Test set for making sure simulations meet certain rules
+#   i.e. frameworks musts be in P1 symmetry to be used in GCMC or Henry
+#   coefficient test
+@testset "Simulation Rules Tests" begin
+    non_P1_framework = Framework("symmetry_test_structure.cif", convert_to_p1=false)
+    molecule = Molecule("CO2")
+    ljff = LJForceField("UFF.csv")
+    temp = 298.0
+    pressure = 5.0
+
+    # make sure that the framework is not in P1 before running tests
+
+    # assertion errors should be thrown when trying to run gcmc simulations
+    #   with non-P1 frameworks
+    @test_throws AssertionError gcmc_simulation(non_P1_framework, molecule,
+                                               temp, pressure, ljff)
+    pressures = [1.0, 5.0, 10.0, 15.0, 20.0]
+    @test_throws AssertionError adsorption_isotherm(non_P1_framework, molecule,
+                                                   temp, pressures, ljff)
+    @test_throws AssertionError stepwise_adsorption_isotherm(non_P1_framework,
+                                            molecule, temp, pressures, ljff)
+
+    # assertion error should be thrown when attempting to run a henry
+    #   coefficient with a non-P1 framework
+    @test_throws AssertionError henry_coefficient(non_P1_framework, molecule,
+                                                 temp, ljff)
+
+    # Test that an assertion is thrown when a non-P1 structure is passed into
+    #   energy_grid(), and compute_accessibility_grid()
+    @test_throws AssertionError energy_grid(non_P1_framework, molecule, ljff)
+    @test_throws AssertionError compute_accessibility_grid(non_P1_framework, molecule, ljff)
+
+    # Test that an assertion is thrown when trying to replicate a non-P1
+    #   structure
+    @test_throws AssertionError replicate(non_P1_framework, (2, 2, 2))
+    @test_throws AssertionError replicate(non_P1_framework, (1, 1, 1))
+
+end
+end