Repository: xchem/fragalysis-api (language: python)

Issue:
Convert Small PDB Chains into HETATM LIGS
Small Cheat to enable users to specify a certain Chain Length to be considered a ligand. This should be used with caution as not all edge cases are handled.

Just creating the PR now to see the conflicts (if the GH Actions are set-up) - still need thorough testing.

Do Not Merge Yet...
New Resample Method
Added Fast mode


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 92e82d43..fd91bf5c 100755
--- a/README.md
+++ b/README.md
@@ -6,14 +6,15 @@
 
 Documentation: https://xchem.github.io/fragalysis-api/
 
-This api aims to allow any user to upload pdb files from the pdb or that they have created themselves, 
+This api aims to allow any user to upload pdb files from the pdb or that they have created themselves,
 and analyse the ligand binding using the fragalysis webpage (https://fragalysis.diamond.ac.uk).
 
 ## Installation
 
-Starting out by initialising an environment and activating it. 
-Clone the repository and cd to the relevant directory. 
+Starting out by initialising an environment and activating it.
+Clone the repository and cd to the relevant directory.
 Install rdkit via conda, and the other dependencies via the setup.py file:
+
 ```bash
 conda create -n fragalysis_env anaconda -y
 conda activate fragalysis_env
@@ -38,114 +39,92 @@ pip install -e .
 cd ..
 ```
 
-You can check if it has installed using:    ```conda list```
+You can check if it has installed using: `conda list`
 
 ### How to use API
 
 1. Set up environment
 2. Download PDB files and query the PDB for structures of the same protein bound to the same or different ligands
-3. Submit PDB files - you will be given a query ID 
+3. Submit PDB files - you will be given a query ID
 4. Push your files into fragalysis and view them online :construction:
-5. Analyse the binding of ligands to your target protein! :construction:
+5. Analyse the binding of ligands to your target protein!
 
 Other functionalities that are available:
+
 - Import pdb files directly through the API
 - Query the pdb for similar structures that also have ligands bound and have the option to import these structures
 - Align files as individual monomers
 - Align volume density files
-
-
-### Enforced rules :scroll:
-
-* The pdb file shall not be greater than 5mb.
-* The pdb filename shall not contain non English language ascii characters 
-and shall be between 4 and 20 characters in length.
-* Each pdb file for alignment shall contain the same number of chains.
-* All pdb files to be aligned must be in the same directory.
-* If manually selecting a file for reference it must be in the same directory as all pdb files for alignment. 
-* PDB file must abide by best practices set out at https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#note5
-* In PDB file, ligands must be referred to by same code in 'HET' lines in header and in 'HETATM' lines in main part of file
-* If providing cpp4 map files: All maps should be of the '.map' or '.cpp4' variety and labelled as:
-    > ```
-    > Aprot.pdb
-    > Aprot_smiles.txt
-    > Aprot_fofc.map
-    > Aprot_2fofc.map
-    > Aprot_event.cpp4
-    > ```
-* If you have multiple events, you can add use a `_event_0.cpp4` type naming convention and increment the number accordingly. Sorry! Here is an example of how you can do it though!
-```bash
-module load ccp4
-cd /staging_dir/target/aligned
-# The directory with the fragalysis stuff. ls should return folders for each ligand
-folders=$(ls)
-for i in $folders
-  do
-    str=$(echo ./$i/)
-    maps=$(find $str -name *.map)
-    ccp4s=$(find $str -name *.ccp4)
-    for j in $maps
-      do
-        mapmask mapin $j mapout $j xyzin ./$i/$i.pdb << eof
-        border 6
-        end
-        eof
-      done
-
-    for j in $ccp4s
-      do
-        mapmask mapin $j mapout $j xyzin ./$i/$i.pdb << eof
-        border 6
-        end
-        eof
-      done
-  done
+- Automatically generate Site annotation and ligand-pair relationships between ligands in same crystal model.
+
+## Usage in Python
+
+To prepare input data-files using python you api can import the `xcimporer` or `import_single_file` functions and then provide the necessary values to the functions.
+
+e.g
+
+```python
+from fragalysis_api import xcimporter, import_single_file
+
+xcimporter(
+    in_dir,
+    out_dir,
+    target,
+    metadata=False,
+    validate=False,
+    reduce_reference_frame=False,
+    biomol=None,
+    covalent=False,
+    pdb_ref="",
+    max_lig_len=0
+)
+# Process a single file:
+import_single_file(
+    in_file,
+    out_dir,
+    target,
+    reduce_reference_frame,
+    reference_pdb,
+    biomol=None,
+    covalent=False,
+    self_ref=False,
+    max_lig_len=0
+)
 ```
 
-### 2.1 How to download PDB files
-
-Move to the fragalysis-api/xcimporter directory. 
-
-You will need two bits of information:
-1. Your 'user ID' - this is your name followed by your protein name. For example, Anna looking at protein ATAD2 wil have username 'Anna_ATAD2'.
-2. The PDB code you would like to download. For example, '6epu'.
-
-To download the PDB file 
+To analyse the context of the sites we can then do:
+
+```python
+from fragalysis_api import Sites, contextualize_crystal_ligands
+import os
+sites_obj = Sites.from_folder(folder=os.path.join(out_dir, target), recalculate=True)
+sites_obj.to_json()
+sites_obj.contextualize_crystal_ligands(folder=os.path.join(out_dir, target))
+
+# To run the recluster using existing sites after import_single_file you can use:
+sites_obj = Sites.from_folder(folder=os.path.join(out_dir, target), recalculate=True)
+sites_obj.cluster_missing_mols(folder=os.path.join(out_dir, target), com_tolerance=5.00, other_tolerance=1.00)
+sites_obj.to_json()
+# Then rerun individual crystal contextualiser
+sites_obj.contextualize_crystal_ligands(folder=os.path.join(out_dir, target))
 ```
-python pdbimporter.py -id [user_id] -pdb [pdb code]
-```
-In our example, this would be 
-```python pdbimporter.py -id Anna_ATAD2 -pdb 6epu```
 
-Alternatively, you can upload your own PDB files. :construction:
+See below for information about the arguments.
 
-### 2.2 How to query the PDB file
+## Command-line Usage
 
-Note: you don't need to manually download a pdb file before querying the PDB for structures of the same protein
+### How to submit PDB files for conversion to a fragalysis friendly format (fff)
 
-Similar to downloading PDB files, however initially all you need is the PDB code and the chain you would like to query. The command is 
 ```
-python pdbquery.py -pdb [pdb code] -chain [required chain]
-```
-The API will then query the PDB for structures of the same protein that also have ligands bound. You will be asked if you would like to see a list of these structures and ligands, and if you would like to download all of these pdb files in bulk. If you choose to download these pdb files, you will be asked for your user id. This is as before: your name followed by the protein name (e.g. Anna_ATAD2).
-
-### 3.How to submit PDB files for conversion to a fragalysis friendly format (fff)
-
-Once you have the files downloaded, they need to be processed before they can be visualised in fragalysis. This is done using
-```
-python xcimporter -id [user id] 
+python fragalysis-api/fragalysis_api/xcimporter/xcimporter.py  -i [input directory] -o [output directory] -t [target name] -m
 ```
-Default directories are used. These can however be changed by using ```-i [input directory] ``` or   ```-o [output directory] ``` if this is required.
 
-For example a more elaborate conversion would be:
-```
-python fragalysis-api/fragalysis_api/xcimporter/xcimporter.py  -i [input directory] -o [output directory] -t [target name] -m 
-```
 Which will align the input files with respects to individual chains inside the .pdb files (`-m`) and save the output
 to a folder specified by `-t [target name]`.
 
 A description of the command line arguments are as follows:
-- `-i`, `--in_dir` :  Input Directory
+
+- `-i`, `--in_dir` : Input Directory
 - `-o`, `--out_dir` : Output Directory
 - `-v`, `--validate`: (Optional) Validate the Inputs
 - `-m`, `--monomerize`: (Optional) Split the input PDBs into separate chains. E.G If a pdb has A and B chains it will create files pdb_name_A.pdb and pdb_name_B.pdb
@@ -154,6 +133,7 @@ A description of the command line arguments are as follows:
 - `-b`, `--biomol`: (Optional) File path to plain text file that contains an optional header that you would like to be added to PDB files.
 - `-r`, `--reference`: (Optional) The name/filepath of the pdb file you which to use as reference (can be PDB ID)
 - `-c`, `--covalent`: (Optional) Handle Covalent attachments by extending output .mol file to include covalent attachment atoms. Requires modified smiles strings.
+- `-mll`, `--max_lig_len`: (Optional) **EXPERIMENTAL** Integer, if >0 will convert all chains with residue length <mll to `HETATM LIG` - useful for converting short chain amino acids to ligands for example.
 
 The terminal will let you know when the conversion has been successful and if there are any files that have been found to be incompatible with the API. We are working to minimize any incompatibilities.
 
@@ -162,19 +142,23 @@ contain two main folders `aligned` and `crystallographic`. The `aligned` folder
 Inside each of these subfolders contains all the information that is required by the fragalysis platform to display and elaborate a given ligand.
 The `crystallographic` folder contains an exact copy of the data supplied to the above command.
 
-#### 3.1 How to submit single PDB files for conversion to a fragalysis friendly format
+### How to submit single PDB files for conversion to a fragalysis friendly format
+
 If you would like to add individual pdb files to the results of the fragalysis api we can use:
+
 ```
 python fragalysis-api/fragalysis_api/xcimporter/single_import.py --in_file=[pdbtobealigned.pdb] --out_dir=[output directory] --target [targetname] -m --reference[output directory/targetname/reference.pdb]
 ```
+
 Where `in_file` is the filepath of a pdb file we would like to align to the rest of our results (with `_smiles.txt` and other files to be aligned in the same directory)
 and `reference` is the pdb file that you would like to align files associated to `in_file`. If you had previous run xcimporter then a reference.pdb file should be located at `output direct/targetname/reference.pdb`
 If `--reference` or `-r` are not specified then `output direct/targetname/reference.pdb` will be used according to what arguments you have specified.
 
 You should be able to use `single_import` without having to delete any preexisting files. As the new outputs should overwrite what previously exists. Nice.
 
-A description of the command line arguments for `single_import.py` are as follows: 
-- `-i`, `--in_file` :  Input File
+A description of the command line arguments for `single_import.py` are as follows:
+
+- `-i`, `--in_file` : Input File
 - `-o`, `--out_dir` : Output Directory
 - `-v`, `--validate`: (Optional) Validate the Inputs
 - `-m`, `--monomerize`: (Optional) Split the input PDBs into separate chains. E.G If a pdb has A and B chains it will create files pdb_name_A.pdb and pdb_name_B.pdb
@@ -184,12 +168,15 @@ A description of the command line arguments for `single_import.py` are as follow
 - `-r`, `--reference`: (Optional) The name/filepath of the pdb file you which to use as reference (can be PDB ID)
 - `-c`, `--covalent`: (Optional) Handle Covalent attachments by extending output .mol file to include covalent attachment atoms. Requires modified smiles strings.
 - `-sr`, `--self_reference`: (Optional) Indicate whether you want pdb files to align to themselves (for testing purposes)
+- `-mll`, `--max_lig_len`: (Optional) **EXPERIMENTAL** Integer, if >0 will convert all chains with residue length <mll to `HETATM LIG` - useful for converting short chain amino acids to ligands for example.
+
+#### Running The Fragalysis API without alignment
 
-#### 3.2 Running The Fragalysis API without alignment 
 This is only available to the `single_import.py` method
-If for whatever reason you decide that you would like to simply split a pdb file without the need of an alignment step you can use the `-sr` or `--selfreference` flags when using the fragalysis API. Or specify the `-r` or `--reference` to be itself. 
+If for whatever reason you decide that you would like to simply split a pdb file without the need of an alignment step you can use the `-sr` or `--selfreference` flags when using the fragalysis API. Or specify the `-r` or `--reference` to be itself.
 
 An example in bash to process a folder of pbds without aligned in bash would be:
+
 ```bash
 $input=/path/to/input/folder
 pdblist=$(ls input)
@@ -201,44 +188,150 @@ done
 
 Notably the `-sr` flag will take precedence over the `-r` flag.
 
+### Site Annotation
+
+We have now introduced automatic site annotations. This can be performed by using the calling sites.py
+
+```
+python fragalysis-api/fragalysis_api/xcimporter/sites.py -i [outputfolder/target] --mode=new -om -cc
+```
+
+This will create a `_sites.json` file in each fragalysis ligand folder.
+
+A description of the arguments are as follows:
+
+- `-i`, `--input` : The input directory for the sites functionality, this will typically be the output fragalyis folder (e.g. `out_dir/target`)
+- `-m`, `--mode` : The mode to run the annotation, either `new` or `missing`. `new` will create entirely new annotations considering all ligands in the input directory, `missing` will preserve previous annotations and append new ligands (such as those brought into existance by single_import)
+- `-om`, `--overwritemeta`: (Optional) Will replace contents of \_meta.csv files with the center of mass site from the `_sites.json` file
+- `-cc`, `--cocrystal`: (Optional) Contextualise multiple ligands from the same crystal together, based on ligand similarity and distance from one another. Information is output to `_relationships.json` file
+- `-ct`, `--comtolerance`: The tolerance value for creating new clusters for centre of mass sites, default is 5.00.
+- `-ot`, `--othertolerance`: The tolerance values for creating new clusters for non-centre of mass sites, default is 1.00.
+
+### Enforced rules :scroll:
+
+- The pdb file shall not be greater than 5mb.
+- The pdb filename shall not contain non English language ascii characters
+  and shall be between 4 and 20 characters in length.
+- Each pdb file for alignment shall contain the same number of chains.
+- All pdb files to be aligned must be in the same directory.
+- If manually selecting a file for reference it must be in the same directory as all pdb files for alignment.
+- PDB file must abide by best practices set out at https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#note5
+- In PDB file, ligands must be referred to by same code in 'HET' lines in header and in 'HETATM' lines in main part of file
+- If providing cpp4 map files: All maps should be of the '.map' or '.cpp4' variety and labelled as:
+  > ```
+  > Aprot.pdb
+  > Aprot_smiles.txt
+  > Aprot_fofc.map
+  > Aprot_2fofc.map
+  > Aprot_event.cpp4
+  > ```
+
+### Cutting maps
+
+For optimal performance in fragalysis-api and localXCR and fragalysis we recommend cutting maps to be around the ligand as to avoid large file sizes.
+You can do whichever way you like but here is an example using ccp4 mapmask.
+
+```bash
+module load ccp4
+cd /staging_dir/target/aligned
+# The directory with the fragalysis stuff. ls should return folders for each ligand
+folders=$(ls)
+for i in $folders
+  do
+    str=$(echo ./$i/)
+    maps=$(find $str -name *.map)
+    ccp4s=$(find $str -name *.ccp4)
+    for j in $maps
+      do
+        mapmask mapin $j mapout $j xyzin ./$i/$i.pdb << eof
+        border 6
+        end
+        eof
+      done
+
+    for j in $ccp4s
+      do
+        mapmask mapin $j mapout $j xyzin ./$i/$i.pdb << eof
+        border 6
+        end
+        eof
+      done
+  done
+```
+
+### How to download PDB files
+
+Move to the fragalysis-api/xcimporter directory.
+
+You will need two bits of information:
+
+1. Your 'user ID' - this is your name followed by your protein name. For example, Anna looking at protein ATAD2 wil have username 'Anna_ATAD2'.
+2. The PDB code you would like to download. For example, '6epu'.
+
+To download the PDB file
+
+```
+python pdbimporter.py -id [user_id] -pdb [pdb code]
+```
+
+In our example, this would be
+`python pdbimporter.py -id Anna_ATAD2 -pdb 6epu`
+
+Alternatively, you can upload your own PDB files. :construction:
+
+### How to query the PDB file
+
+Note: you don't need to manually download a pdb file before querying the PDB for structures of the same protein
+
+Similar to downloading PDB files, however initially all you need is the PDB code and the chain you would like to query. The command is
+
+```
+python pdbquery.py -pdb [pdb code] -chain [required chain]
+```
+
+The API will then query the PDB for structures of the same protein that also have ligands bound. You will be asked if you would like to see a list of these structures and ligands, and if you would like to download all of these pdb files in bulk. If you choose to download these pdb files, you will be asked for your user id. This is as before: your name followed by the protein name (e.g. Anna_ATAD2).
 
 ### Preparing metadata for fragalysis upload
-Each fragalysis upload requires a .csv file (called `metadata.csv`) to be stored inside the output `aligned`. It should have the following structure:
 
+We recommend using our own tools [localXCR](https://github.com/xchem/localXCR) to annotate and visualise the results of the fragalysis-api.
+
+However you are able to do this manually if you want to.
+
+Each fragalysis upload requires a .csv file (called `metadata.csv`) to be stored inside the output `aligned`. It should have the following structure:
 
-|<empty>|crystal_name|RealCrystalName|smiles|new_smiles|alternate_name|site_name|pdb_entry|
-|---|---|---|---|---|---|---|---|
-|1|prot-x0001_0A|prot-x0001|CCCC|<empty>|<empty>|Binding Site|<empty>|
-|...|...|...|...|...|...|...|...|
+| <empty> | crystal_name  | RealCrystalName | smiles | new_smiles | alternate_name | site_name    | pdb_entry |
+| ------- | ------------- | --------------- | ------ | ---------- | -------------- | ------------ | --------- |
+| 1       | prot-x0001_0A | prot-x0001      | CCCC   | <empty>    | <empty>        | Binding Site | <empty>   |
+| ...     | ...           | ...             | ...    | ...        | ...            | ...          | ...       |
 
 This file will be generated by the api if using the `-md` flag or otherwise can be constructed manually.
-- [First Column]: This doesn't have a column name, but just indicates a rownumber, fill 1 to n  
-- __crystal_name__ (Required): The name of the crystal as it appear inside the aligned folder. e.g. prot-x0001_0A or prot-x0001_1A
-- __RealCrystalName__ (Required): The name of the crystal that the ligand is derived from. e.g. prot-x0001, prot-x0002
-- __smiles__ (Required): The SMILES string for the crystal as described in the `aligned/crystal_name/crystal_name_smiles.txt` file
-- __new_smiles__ (can be left empty): If the SMILES string is slightly wrong or you want to specify a specific enantiomer to render, declare the new SMILES string here
-- __alternate_name__ (can be left empty): If you have a particular title that you would like to render alongside the crystal_name
-- __site_name__ (Recommended): This is the most important field. All rows in this column need a value! e.g. Binding Site. So that Fragalysis will know where to cluster the fragments together. This will be rendered alphabetically so if you want to preserved a specific order of sites, prefix the site labels with alphanumeric combinations e.g. "A1 - Binding Site", "B2 - Allosteric Site". If entire column is left empty then fragalysis will cluster ligands together by center of mass and will use uninformative site labels.
-- __pdb_entry__ (can be left empty): The PDB code relating to the Crystal.
+
+- [First Column]: This doesn't have a column name, but just indicates a rownumber, fill 1 to n
+- **crystal_name** (Required): The name of the crystal as it appear inside the aligned folder. e.g. prot-x0001_0A or prot-x0001_1A
+- **RealCrystalName** (Required): The name of the crystal that the ligand is derived from. e.g. prot-x0001, prot-x0002
+- **smiles** (Required): The SMILES string for the crystal as described in the `aligned/crystal_name/crystal_name_smiles.txt` file
+- **new_smiles** (can be left empty): If the SMILES string is slightly wrong or you want to specify a specific enantiomer to render, declare the new SMILES string here
+- **alternate_name** (can be left empty): If you have a particular title that you would like to render alongside the crystal_name
+- **site_name** (Recommended): This is the most important field. All rows in this column need a value! e.g. Binding Site. So that Fragalysis will know where to cluster the fragments together. This will be rendered alphabetically so if you want to preserved a specific order of sites, prefix the site labels with alphanumeric combinations e.g. "A1 - Binding Site", "B2 - Allosteric Site". If entire column is left empty then fragalysis will cluster ligands together by center of mass and will use uninformative site labels.
+- **pdb_entry** (can be left empty): The PDB code relating to the Crystal.
 
 The rendered table line in plain text will look like this:
+
 ```
 ,crystal_name,RealCrystalName,smiles,new_smiles,alternate_name,site_name,pdb_entry
 1,prot-x0001_0A,prot-x0001,CCCC,,,Binding Site,
 ```
+
 ### Who we are
+
 We are the Fragment 5, a group of students at the University of Oxford.
 
-Anna :whale:   
+Anna :whale:  
 Maranga :fire:  
 George a.k.a Joff Boff :poop:  
 Tobias :cow:  
-Alister :panda_face:  
+Alister :panda_face:
 
-We are looked after by Rachael :crown:. 
+We are looked after by Rachael :crown:.
 
 Finally, the code quality was ruined by Tyler.
-
-
-
-
diff --git a/fragalysis_api/__init__.py b/fragalysis_api/__init__.py
index 76213984..935c9896 100755
--- a/fragalysis_api/__init__.py
+++ b/fragalysis_api/__init__.py
@@ -1,10 +1,12 @@
 from .xcglobalscripts.set_config import ConfigSetup
 
 from .xcimporter.validate import Validate, ValidatePDB
-from .xcimporter.conversion_pdb_mol import Ligand, pdb_apo, set_up
-from .xcimporter.align import Align, Monomerize
+from .xcimporter.conversion_pdb_mol import set_up, convert_small_AA_chains, copy_extra_files
+from .xcimporter.align import Align
 from .xcimporter.xc_utils import to_fragalysis_dir
+from .xcimporter.sites import Sites, contextualize_crystal_ligands
 from .xcimporter.xcimporter import xcimporter
+from .xcimporter.single_import import import_single_file
 from .xcextracter.getdata import GetTargetsData, GetMoleculesData, GetPdbData, GetMolgroupData
 from .xcextracter.frag_web_live import can_connect
 from .xcextracter.xcextracter import xcextracter
diff --git a/fragalysis_api/xcimporter/align.py b/fragalysis_api/xcimporter/align.py
index f121b097..6d579b7c 100755
--- a/fragalysis_api/xcimporter/align.py
+++ b/fragalysis_api/xcimporter/align.py
@@ -10,30 +10,34 @@
 import os
 import warnings
 import Bio.PDB.PDBExceptions as bpp
-import json
 import shutil
 import gemmi  # Oh boy...
-import numpy
+import numpy as np
+import json
 
 warnings.simplefilter('ignore', bpp.PDBConstructionWarning)
 
 
 class Align:
 
-    def __init__(self, directory, pdb_ref='', mono=False):
+    def __init__(self, directory, pdb_ref='', rrf=False, refset=True):
         '''
         :param directory: Directory path contain pdbs to be aligned.
         :param pdb_ref: The String Reference of a pdb that you want to
             optionally align the input pdbs to. This should be the name of the PDB without
             .pbd and should be located inside the specified directory.
-        :param mono: Bool, To indicate if chains have been split into seperate pdb files.
+        :param rrf: Bool, To indicate if chains have been split into seperate pdb files.
             If True, this will change the behaviour of the alignment to solely align the first
             chain of each structure together, and not consider other chains.
+        :param refset: Boolean, Indicate whether or not to automatically set a reference when constructing the alignment class
+            Can be set to False if only using align_to_reference function.
         '''
 
         self.directory = directory
-        self._get_ref = pdb_ref
-        self.mono = mono
+        if refset:
+            self._get_ref = pdb_ref
+        self.rrf = rrf
+        self.ref_map = ''
 
     @property
     def _get_files(self):
@@ -72,13 +76,16 @@ def _get_ref(self, pdb_ref):
         """
         if pdb_ref != '':
             try:
-                assert (os.path.isfile(os.path.join(self.directory, pdb_ref + ".pdb")) == 1)
+                assert (os.path.isfile(os.path.join(
+                    self.directory, pdb_ref + ".pdb")) == 1)
                 self.__pdb_ref = pdb_ref
             except AssertionError:
                 print('pdb desired as reference does not exists. Default pdb chosen.')
-                self.__pdb_ref = self.__best_length_and_resolution([i for i in self._get_files if 'pdb' in i])
+                self.__pdb_ref = self.__best_length_and_resolution(
+                    [i for i in self._get_files if 'pdb' in i])
         else:
-            self.__pdb_ref = self.__best_length_and_resolution([i for i in self._get_files if 'pdb' in i])
+            self.__pdb_ref = self.__best_length_and_resolution(
+                [i for i in self._get_files if 'pdb' in i])
 
     def __best_length_and_resolution(self, pdb_files):
         """
@@ -88,7 +95,8 @@ def __best_length_and_resolution(self, pdb_files):
         :return: str of filename with best .pdb file
         """
         a_df = pd.DataFrame(pdb_files, columns=['file'])
-        a_df[['resolution', 'p_len', 'pdb']] = a_df.file.apply(self.__get_length_and_resolution)
+        a_df[['resolution', 'p_len', 'pdb']] = a_df.file.apply(
+            self.__get_length_and_resolution)
 
         return a_df.sort_values(by=['p_len', 'resolution'], ascending=[False, True]).pdb[0]
 
@@ -100,65 +108,105 @@ def __get_length_and_resolution(self, file):
         """
         parser = bp.PDBParser()
         ppb = bp.PPBuilder()
-        structure = parser.get_structure(os.path.splitext(os.path.basename(file))[0], file)
+        structure = parser.get_structure(
+            os.path.splitext(os.path.basename(file))[0], file)
 
         seq_len = 0
 
-        for pp in ppb.build_peptides(structure):  # Retrieve length by looping through each chain in the protein
+        # Retrieve length by looping through each chain in the protein
+        for pp in ppb.build_peptides(structure):
             seq_len += len(pp.get_sequence())
 
         # using a functions from PDBParser parser class to get the resolution and protein id from the pdb file
         return pd.Series([structure.header['resolution'], seq_len, structure.id])
 
-    def align_to_reference(self, in_file, reference, out_dir):
+    def align_to_reference(self, in_file, reference_pdb, out_dir):
         '''
         Aligns a single pdb file to a reference and adds it to the specified output directory.
         :param in_file: filepath to corresponding pdb file to align. Accompanying map files should be located within
             the same directory as input file.
-        :param reference: the reference pdb file to align to.
+        :param reference_pdb: the reference pdb file to align to.
         :param out_dir: The desired output directory for the aligned pdb file.
         :return: an aligned pdb file in the output directory with the same name as input.
         '''
         input_files = in_file
         map_list = self._get_maplist
         base_names = os.path.splitext(os.path.basename(input_files))[0]
-        crystals = [y for y in [x for x in [base_names] if 'event' not in x] if 'fofc' not in y]
-        ref = reference
+        crystals = [y for y in [x for x in [base_names]
+                                if 'event' not in x] if 'fofc' not in y]
+        ref = reference_pdb
         dir = self.directory
-        mono = self.mono
+        rrf = self.rrf
+
+        reference_pdb = Structure.from_file(file=Path(ref))
+        s = time.time()
         for num, name in enumerate(crystals):
-            reference_pdb = Structure.from_file(file=Path(ref))
             all_maps = [j for j in map_list if name in j]
             current_pdb = Structure.from_file(file=Path(in_file))
-            try:
-                current_pdb, transform = current_pdb.align_to(other=reference_pdb, monomerized=mono)
-                print(transform)
-                print(transform.transform.mat.tolist())
-                print(transform.transform.vec.tolist())
-            except Exception as e:
-                # Poorly Documented. Use better stuff...
-                print(f'{e}')
-                continue
+            if rrf:
+                # current_pdb.structure, chains = split_chain_str(
+                chains = split_chain_str(os.path.join(dir, f'{name}.pdb'))
+            else:
+                chains = ['']
+            for chain in chains:
+                if rrf:
+                    print(
+                        f'Aligning Chain {chain} of {name} to first chain of {ref}')
+                current_pdb = Structure.from_file(file=Path(in_file))
+                try:
+                    current_pdb, transform = current_pdb.align_to(
+                        other=reference_pdb, rrf=rrf, chain_id=chain
+                    )
+                except Exception as e:
+                    print(f'{e}')
+                    continue
+                # Write new structure according to chain-name?
+                if rrf:
+                    current_pdb.structure.write_pdb(
+                        os.path.join(out_dir, f'{name}_{chain}_bound.pdb')
+                    )
+                    transform.to_json(filename=os.path.join(
+                        out_dir, f'{name}_{chain}_transform.json'))
+                    if os.path.exists(os.path.join(self.directory, f'{name}_smiles.txt')):
+                        shutil.copyfile(os.path.join(self.directory, f'{name}_smiles.txt'), os.path.join(
+                            out_dir, f'{name}_{chain}_smiles.txt'))
+                else:
+                    current_pdb.structure.write_pdb(
+                        os.path.join(out_dir, f'{name}_bound.pdb')
+                    )
+                    transform.to_json(filename=os.path.join(
+                        out_dir, f'{name}_transform.json'))
+                    if os.path.exists(os.path.join(self.directory, f'{name}_smiles.txt')):
+                        shutil.copyfile(os.path.join(self.directory, f'{name}_smiles.txt'), os.path.join(
+                            out_dir, f'{name}_smiles.txt'))
 
-            current_pdb.structure.write_pdb(os.path.join(out_dir, f'{name}_bound.pdb'))
-            # Align Xmaps + save!
-            s = time.time()
-            print(all_maps)
-            for i in all_maps:
-                base, ext = os.path.splitext(os.path.basename(i))
-                print(i)
-                s2 = time.time()
-                map = Xmap.from_file(file=Path(os.path.join(dir, f'{base}{ext}')))
-                array = np.array(map.xmap, copy=False)
-                array[~np.isfinite(array)] = 0                
-                #print(map.to_array())
-                map.resample(xmap=map, transform=transform)
-                #print(map.to_array())
-                map.save(path=Path(os.path.join(out_dir, f'{base}{ext}')))
-                e2 = time.time()
-                print(f'{int(e2 - s2)} seconds to transform map...')
-                e = time.time()
-                print(f'Total Running time: {int(e - s) / 60} minutes.')
+                # Align Xmaps + save!
+                for i in all_maps:
+                    base, ext = os.path.splitext(os.path.basename(i))
+                    s2 = time.time()
+                    map = Xmap.from_file(
+                        file=Path(os.path.join(dir, f'{base}{ext}')))
+                    array = np.array(map.xmap, copy=False)
+                    array[~np.isfinite(array)] = 0
+                    array_mean = np.mean(array)
+                    array_sd = np.std(array)
+                    array[:, :, :] = (array[:, :, :] - array_mean) / array_sd
+                    newmap = resample(
+                        moving_xmap=map, transform=transform, reference_structure=reference_pdb)
+                    template = Path(
+                        os.path.join(dir, f'{base}{ext}'))
+                    if rrf:
+                        base = base.replace(name, f'{name}_{chain}')
+                    fn = f'{base}{ext}'
+                    referenceSave(
+                        template_map_path=template,
+                        xmap=newmap,
+                        path_to_save=Path(os.path.join(out_dir, fn))
+                    )
+                    e2 = time.time()
+                    print(f'{int(e2 - s2)} seconds to transform map...')
+                    e = time.time()
+                    print(f'Total Running time: {int(e - s) / 60} minutes.')
 
     def write_align_ref(self, output):
         '''
@@ -167,7 +215,7 @@ def write_align_ref(self, output):
         :return: the pdbfile that was chosen as reference copied to the located specified by output
         '''
         fn = os.path.join(self.directory, f'{self._get_ref}.pdb')
-        shutil.copyfile(fn, output)
+        shutil.copyfile(fn, os.path.join(output, 'reference.pdb'))
 
     def align(self, out_dir):
         """
@@ -178,53 +226,89 @@ def align(self, out_dir):
         # load ref
         input_files = self._get_files
         map_list = self._get_maplist
-        crystals = [os.path.splitext(os.path.basename(f))[0] for f in input_files]
+        crystals = [os.path.splitext(os.path.basename(f))[0]
+                    for f in input_files]
         ref = self._get_ref
         dir = self.directory
-        mono = self.mono
+        rrf = self.rrf
+
+        # Reference stuff
+        reference_pdb = Structure.from_file(
+            file=Path(os.path.join(dir, f'{ref}.pdb')))
 
         s = time.time()
         for num, name in enumerate(crystals):
-            reference_pdb = Structure.from_file(file=Path(os.path.join(dir, f'{ref}.pdb')))
             all_maps = [j for j in map_list if name in j]
-            if not name == ref:
-                # Do an alignment + save
-                current_pdb = Structure.from_file(file=Path(os.path.join(dir, f'{name}.pdb')))
+            # Logic to do...
+            # Align Chain N to First Chain in Reference
+            current_pdb = Structure.from_file(
+                file=Path(os.path.join(dir, f'{name}.pdb'))
+            )
+            if rrf:
+                chains = split_chain_str(
+                    os.path.join(dir, f'{name}.pdb'))
+            else:
+                chains = ['']
+            for chain in chains:
+                current_pdb = Structure.from_file(
+                    file=Path(os.path.join(dir, f'{name}.pdb')))
                 try:
-                    current_pdb, transform = current_pdb.align_to(other=reference_pdb, monomerized=mono)
+                    current_pdb, transform = current_pdb.align_to(
+                        other=reference_pdb, rrf=rrf, chain_id=chain
+                    )
                 except Exception as e:
-                    # Poorly Documented. Use better stuff...
                     print(f'{e}')
                     continue
 
-                current_pdb.structure.write_pdb(os.path.join(out_dir, f'{name}_bound.pdb'))
+                # Write new structure according to chain-name?
+                if rrf:
+                    current_pdb.structure.write_pdb(
+                        os.path.join(out_dir, f'{name}_{chain}_bound.pdb')
+                    )
+                    transform.to_json(filename=os.path.join(
+                        out_dir, f'{name}_{chain}_transform.json'))
+                    if os.path.exists(os.path.join(self.directory, f'{name}_smiles.txt')):
+                        shutil.copyfile(os.path.join(self.directory, f'{name}_smiles.txt'), os.path.join(
+                            out_dir, f'{name}_{chain}_smiles.txt'))
+                else:
+                    current_pdb.structure.write_pdb(
+                        os.path.join(out_dir, f'{name}_bound.pdb')
+                    )
+                    transform.to_json(filename=os.path.join(
+                        out_dir, f'{name}_transform.json'))
+                    if os.path.exists(os.path.join(self.directory, f'{name}_smiles.txt')):
+                        shutil.copyfile(os.path.join(self.directory, f'{name}_smiles.txt'), os.path.join(
+                            out_dir, f'{name}_smiles.txt'))
+
                 # Align Xmaps + save!
                 for i in all_maps:
                     base, ext = os.path.splitext(os.path.basename(i))
-                    print(i)
-                    # self.read_reshape_resave(name=base, out_dir=out_dir, ext=ext, transform=transform)
                     s2 = time.time()
-                    map = Xmap.from_file(file=Path(os.path.join(dir, f'{base}{ext}')))
-                    map.resample(xmap=map, transform=transform)
-                    map.save(path=Path(os.path.join(out_dir, f'{base}{ext}')))
+                    map = Xmap.from_file(
+                        file=Path(os.path.join(dir, f'{base}{ext}')))
+                    array = np.array(map.xmap, copy=False)
+                    array[~np.isfinite(array)] = 0
+                    array_mean = np.mean(array)
+                    array_sd = np.std(array)
+                    array[:, :, :] = (array[:, :, :] - array_mean) / array_sd
+                    newmap = resample(
+                        moving_xmap=map, transform=transform, reference_structure=reference_pdb)
+                    template = Path(
+                        os.path.join(dir, f'{base}{ext}'))
+                    if rrf:
+                        base = base.replace(name, f'{name}_{chain}')
+                    fn = f'{base}{ext}'
+                    referenceSave(
+                        template_map_path=template,
+                        xmap=newmap,
+                        path_to_save=Path(os.path.join(out_dir, fn))
+                    )
                     e2 = time.time()
                     print(f'{int(e2 - s2)} seconds to transform map...')
-
                     e = time.time()
                     print(f'Total Running time: {int(e - s) / 60} minutes.')
 
-            else:
-                shutil.copyfile(os.path.join(self.directory, f'{name}.pdb'),
-                                os.path.join(out_dir, f'{name}_bound.pdb'))
-                if os.path.exists(os.path.join(self.directory, f'{name}_smiles.txt')):
-                    shutil.copyfile(os.path.join(self.directory, f'{name}_smiles.txt'),
-                                    os.path.join(out_dir, f'{name}_smiles.txt'))
-                for i in all_maps:
-                    base, ext = os.path.splitext(os.path.basename(i))
-                    shutil.copyfile(i, os.path.join(out_dir, f"{base}{ext}"))
 
-
-# Conor's stuff
 @dataclasses.dataclass()
 class ResidueID:
     model: str
@@ -245,7 +329,8 @@ class RFree:
 
     @staticmethod
     def from_structure(structure):
-        rfree = structure.structure.make_mmcif_document()[0].find_loop("_refine.ls_R_factor_R_free")[0]
+        rfree = structure.structure.make_mmcif_document(
+        )[0].find_loop("_refine.ls_R_factor_R_free")[0]
         return RFree(float(rfree))
 
     def to_float(self):
@@ -255,29 +340,52 @@ def to_float(self):
 @dataclasses.dataclass()
 class Transform:
     transform: gemmi.Transform
-    com_reference: numpy.array
-    com_moving: numpy.array
+    com_reference: np.array
+    com_moving: np.array
 
     def apply(self, position: gemmi.Position) -> gemmi.Position:
         rotation_frame_position = gemmi.Position(position[0] - self.com_reference[0],
-                                                 position[1] - self.com_reference[1],
+                                                 position[1] -
+                                                 self.com_reference[1],
                                                  position[2] - self.com_reference[2])
         transformed_vector = self.transform.apply(rotation_frame_position)
         transformed_position = gemmi.Position(transformed_vector[0] + self.com_moving[0],
-                                              transformed_vector[1] + self.com_reference[1],
-                                              transformed_vector[2] + self.com_reference[2])
+                                              transformed_vector[1] +
+                                              self.com_moving[1],
+                                              transformed_vector[2] + self.com_moving[2])
         return transformed_position
 
     def apply_inverse(self, position: gemmi.Position) -> gemmi.Position:
         rotation_frame_position = gemmi.Position(position[0] - self.com_moving[0],
-                                                 position[1] - self.com_moving[1],
+                                                 position[1] -
+                                                 self.com_moving[1],
                                                  position[2] - self.com_moving[2])
         transformed_vector = self.transform.inverse().apply(rotation_frame_position)
         transformed_position = gemmi.Position(transformed_vector[0] + self.com_reference[0],
-                                              transformed_vector[1] + self.com_reference[1],
+                                              transformed_vector[1] +
+                                              self.com_reference[1],
                                               transformed_vector[2] + self.com_reference[2])
         return transformed_position
 
+    @staticmethod
+    def from_json(json_file):
+        with open(json_file) as jfile:
+            data = json.load(jfile)
+        transform = gemmi.Transform()
+        transform.mat.fromlist(data['transform_mat'])
+        transform.vec.fromlist(data['transform_vec'])
+        return Transform(transform, data['com_reference'], data['com_moving'])
+
+    def to_json(self, filename):
+        data = {
+            'com_moving': self.com_moving.tolist(),
+            'com_reference': self.com_reference.tolist(),
+            'transform_mat': self.transform.mat.tolist(),
+            'transform_vec': self.transform.vec.tolist()
+        }
+        with open(filename, 'w') as f:
+            json.dump(data, f)
+
     @staticmethod
     def from_translation_rotation(translation, rotation, com_reference, com_moving):
         transform = gemmi.Transform()
@@ -330,238 +438,151 @@ def all_atoms(self):
                     for atom in residue:
                         yield atom
 
-    def align_to(self, other, monomerized=False):
+    def align_to(self, other, rrf=False, chain_id=''):
+        # TODO: CHANGE
         # Warning: inplace!
         # Aligns structures usings carbon alphas and transform self into the frame of the other
 
         ca_self = []
         ca_other = []
-
+        ca_res_self_names = []
+        ca_res_other_names = []
         # Get CAs
         for model in self.structure:
             for chain in model:
-                for res_self in chain.get_polymer():
-                    if 'LIG' in str(res_self):
-                        print('Skipping Ligand...')
-                        continue
-
-                    try:
-                        current_res_id = ResidueID.from_residue_chain(model, chain, res_self)
-                        if monomerized:
-                            # print(other.structure[current_res_id.model])
-                            # print(len(other.structure[current_res_id.model]))
-                            res_other = other.structure[current_res_id.model][0][current_res_id.insertion][0]
-                        else:
-                            res_other = \
-                                other.structure[current_res_id.model][current_res_id.chain][current_res_id.insertion][0]
-                        # print(f'{self.structure}|{res_self}')
-                        # print(f'{other.structure}|{res_other}')
-                        self_ca_pos = res_self["CA"][0].pos
-                        other_ca_pos = res_other["CA"][0].pos
-
-                    except:
-                        print('Skipping, Residue not found in chain')
-                        continue
-
-                    ca_list_self = Transform.pos_to_list(self_ca_pos)
-                    ca_list_other = Transform.pos_to_list(other_ca_pos)
-
-                    ca_self.append(ca_list_self)
-                    ca_other.append(ca_list_other)
+                if rrf and chain.name not in chain_id:
+                    continue
+                else:
+                    for res_self in chain.get_polymer():
+                        if 'LIG' in str(res_self):
+                            continue
+                        try:
+                            current_res_id = ResidueID.from_residue_chain(
+                                model, chain, res_self)
+                            if rrf:  # TODO CHANGE?
+                                res_other = other.structure[current_res_id.model][0][current_res_id.insertion][0]
+                            else:
+                                res_other = \
+                                    other.structure[current_res_id.model][current_res_id.chain][current_res_id.insertion][0]
+                            self_ca_pos = res_self["CA"][0].pos
+                            other_ca_pos = res_other["CA"][0].pos
+                            ca_res_self_names.append(res_self)
+                            ca_res_other_names.append(res_other)
+                        except:
+                            continue
+
+                        ca_list_self = Transform.pos_to_list(self_ca_pos)
+                        ca_list_other = Transform.pos_to_list(other_ca_pos)
+
+                        ca_self.append(ca_list_self)
+                        ca_other.append(ca_list_other)
 
         # Make coord matricies
-        matrix_self = numpy.array(ca_self)
-        matrix_other = numpy.array(ca_other)
+        matrix_self = np.array(ca_self)
+        matrix_other = np.array(ca_other)
 
         # Find means
-        mean_self = numpy.mean(matrix_self, axis=0)
-        mean_other = numpy.mean(matrix_other, axis=0)
+        mean_self = np.mean(matrix_self, axis=0)
+        mean_other = np.mean(matrix_other, axis=0)
 
         # demaen
         de_meaned_self = matrix_self - mean_self
         de_meaned_other = matrix_other - mean_other
 
         # Align
-        rotation, rmsd = spatial.transform.Rotation.align_vectors(de_meaned_self, de_meaned_other)
+        rotation, rmsd = spatial.transform.Rotation.align_vectors(
+            de_meaned_self, de_meaned_other)
 
         # Get transform
-        vec = numpy.array([0.0, 0.0, 0.0])
+        vec = np.array([0.0, 0.0, 0.0])
         # Transform is from other frame to self frame
         transform = Transform.from_translation_rotation(vec,
                                                         rotation,
                                                         mean_other,
                                                         mean_self,
                                                         )
-
         # Transform positions
         for atom in self.all_atoms():
-            # print(str(atom))
             atom.pos = transform.apply_inverse(atom.pos)
 
         return self, transform
 
 
-class Monomerize:
-
-    def __init__(self, directory, outdir):
-        '''
-        :param directory: Path to folder containing .pdb files to monomerize
-        :param outdir: Output folder where monomerized pdbs will be saved
-        :return: Initialises the Monomerize class that has the ability to split pdb files by chain.
-        '''
-        self.directory = directory
-        self.outdir = outdir
-        self.non_ligs = json.load(
-            open(os.path.join(os.path.dirname(__file__), "non_ligs.json"), "r")
-        )
-
-    def get_filelist(self):
-        '''
-        Get a list of pdb files within the input directory
-        :return: Returns a list of pdb files
-        '''
-        return glob.glob(os.path.join(self.directory, '*.pdb'))
-
-    def get_maplist(self):
-        '''
-        Get a list of .map or .ccp4 files from within the input directory
-        :return: Returns a list of filepaths corresponding to said volume density files.
-        '''
-        map_files = glob.glob(os.path.join(self.directory, '*.map'))
-        cpp4_files = glob.glob(os.path.join(self.directory, '*.ccp4'))
-        return map_files + cpp4_files
-
-    def split_chains(self, f):
-        '''
-        Split a pdb file according to chain names. While preserving chains that contain waters and assigned small, non-amino acid containing chains to the nearest chain by center of mass.
-        :param f: File name of the pdb file to split
-        :return: A list of pdb_files that have been created with the naming convention f_[chain_name]
-        '''
-        aa_codes = {'V': 'VAL', 'I': 'ILE', 'L': 'LEU', 'E': 'GLU', 'Q': 'GLN', 'D': 'ASP', 'N': 'ASN', 'H': 'HIS',
-                        'W': 'TRP', 'F': 'PHE', 'Y': 'TYR', 'R': 'ARG', 'K': 'LYS', 'S': 'SER', 'T': 'THR', 'M': 'MET',
-                        'A': 'ALA', 'G': 'GLY', 'P': 'PRO', 'C': 'CYS'}
-        filenames = []
-        base_structure = gemmi.read_structure(f)
-        base_models = base_structure[0]
-        all_chains = [x.name for x in base_models]
-        HOH_chains = find_water_chains(f)
-
-        nonHOH_chains = list(set(all_chains) - set(HOH_chains))
-
-        chain_centers = {}
-        chain_names = []
-        # From nonHOH chains, identify if they contain any aminoacids, if they do. Assume they are full chains.
-        # else, later on they will be attached to nearest chain by center.
-        for i in nonHOH_chains:
-            chain_centers[i] = get_chain_center(chain_name=i, file=f)
-            span = base_structure[0][i].whole()
-            residues = [x.name for x in span]
-            # Possible to convert this to a %age
-            if any([True for v in residues if v in aa_codes.values()]):
-                chain_names.append(i)
-
-        alt_chains = list(set(nonHOH_chains) - set(chain_names))
-        for i in chain_names:
-            # For each chain, convert all ligands,
-            temp_structure = gemmi.read_structure(f)
-            for j in alt_chains:
-                chain_dists = {}
-                for z in chain_names:
-                    chain_dists[z] = chain_centers[j].dist(chain_centers[z])
-                print(chain_dists)
-                temp_structure[0][j].name = min(chain_dists, key=chain_dists.get)
-
-            leftover_chains = list(set(nonHOH_chains) - set(i))
-            # Remove remaining chains
-            for j in leftover_chains:
-                if j == i:
-                    # Just in case.
-                    continue
-                else:
-                    temp_structure[0].remove_chain(j)
-                    print([x.name for x in temp_structure[0]])
-
-            # Rename Chain to corresponding chain then save!
-            name = os.path.splitext(os.path.basename(f))[0] + '_' + str(i)
-            filename = os.path.join(self.outdir, f'{name}_mono.pdb')
-            print(f'Writing to: {filename}')
-            temp_structure.write_pdb(filename)
-            filenames.append(filename)
-
-        return filenames
-
-    def process_pdb(self, filename, maplist):
-        '''
-        Processes a pdb file.
-        Firstly it splits the pdb file according to chains and then makes a copy all accompanying files (e.g. map files) to share the same name as the new split pdb file.
-        :param filename: The filepath of the pdb file to process
-        :param maplist: A list of mapfiles to be copied for the number of chains the filename may have.
-        :return: A list of filenames that were generated from the input pdb.
-        '''
-        out_names = self.split_chains(filename)
-        for o in out_names:
-            if os.path.isfile(filename.replace('.pdb', '_smiles.txt')):
-                shutil.copy(filename.replace('.pdb', '_smiles.txt'), o.replace('_mono.pdb', '_smiles.txt'))
-            base = os.path.splitext(os.path.basename(filename))[0]
-            new = os.path.splitext(os.path.basename(o))[0].replace('_mono', '')
-            allmaps = [j for j in maplist if base in j]
-            for map in allmaps:
-                if os.path.isfile(map):
-                    mapbase = os.path.basename(map)
-                    shutil.copy(map, os.path.join(self.outdir, mapbase.replace(base, new)))
-        return out_names
-
-    def write_bound(self, inname, outname):
-        '''
-        Appends header to pdb file and writes to new location
-        :param inname: input pdb file
-        :param outname: filepath of desired output
-        :return: Formally returns nothing, but a file will be written to outname
-        '''
-        with open(inname, 'r') as handle:
-            switch = 0
-            header_front, header_end = [], []
-
-            for line in handle:
-
-                if line.startswith('ATOM'): switch = 1
-
-                if line.startswith('HETATM'): switch = 2
-
-                if switch == 0 and not line.startswith('REMARK 350'): header_front.append(line)
-
-                if (switch == 2) and not line.startswith('HETATM'): header_end.append(line)
-
-        for o in outname:
-            newfile_contents = open(o, 'r').readlines()
+def split_chain_str(f):
+    aa_codes = {'V': 'VAL', 'I': 'ILE', 'L': 'LEU', 'E': 'GLU', 'Q': 'GLN', 'D': 'ASP', 'N': 'ASN', 'H': 'HIS',
+                'W': 'TRP', 'F': 'PHE', 'Y': 'TYR', 'R': 'ARG', 'K': 'LYS', 'S': 'SER', 'T': 'THR', 'M': 'MET',
+                'A': 'ALA', 'G': 'GLY', 'P': 'PRO', 'C': 'CYS'}
+    base_structure = gemmi.read_structure(f)
+    base_models = base_structure[0]
+    all_chains = [x.name for x in base_models]
+    HOH_chains = find_water_chains(f)
+    nonHOH_chains = list(set(all_chains) - set(HOH_chains))
+    chain_centers = {}
+    chain_names = []
+    for i in nonHOH_chains:
+        chain_centers[i] = get_chain_center(chain_name=i, file=f)
+        span = base_structure[0][i].whole()
+        residues = [x.name for x in span]
+        # Possible to convert this to a %age
+        if any([True for v in residues if v in aa_codes.values()]):
+            chain_names.append(i)
+    alt_chains = list(set(nonHOH_chains) - set(chain_names))
+    return list(set(nonHOH_chains) - set(alt_chains))
+
+
+def split_chains(f):
+    '''
+    Split a pdb file according to chain names. While preserving chains that contain waters and assigned small, non-amino acid containing chains to the nearest chain by center of mass.
+    :param f: File name of the pdb file to split
+    :return: A list of pdb_files that have been created with the naming convention f_[chain_name]
+    '''
+    aa_codes = {'V': 'VAL', 'I': 'ILE', 'L': 'LEU', 'E': 'GLU', 'Q': 'GLN', 'D': 'ASP', 'N': 'ASN', 'H': 'HIS',
+                'W': 'TRP', 'F': 'PHE', 'Y': 'TYR', 'R': 'ARG', 'K': 'LYS', 'S': 'SER', 'T': 'THR', 'M': 'MET',
+                'A': 'ALA', 'G': 'GLY', 'P': 'PRO', 'C': 'CYS'}
+    filenames = []
+    base_structure = gemmi.read_structure(f)
+    base_models = base_structure[0]
+    all_chains = [x.name for x in base_models]
+    HOH_chains = find_water_chains(f)
+    nonHOH_chains = list(set(all_chains) - set(HOH_chains))
+
+    chain_centers = {}
+    chain_names = []
+    # From nonHOH chains, identify if they contain any aminoacids, if they do. Assume they are full chains.
+    # else, later on they will be attached to nearest chain by center.
+    for i in nonHOH_chains:
+        chain_centers[i] = get_chain_center(chain_name=i, file=f)
+        span = base_structure[0][i].whole()
+        residues = [x.name for x in span]
+        # Possible to convert this to a %age
+        if any([True for v in residues if v in aa_codes.values()]):
+            chain_names.append(i)
+
+    alt_chains = list(set(nonHOH_chains) - set(chain_names))
+    for i in chain_names:
+        # For each chain, convert all ligands,
+        temp_structure = gemmi.read_structure(f)
+        for j in alt_chains:
+            chain_dists = {}
+            for z in chain_names:
+                chain_dists[z] = chain_centers[j].dist(chain_centers[z])
+            temp_structure[0][j].name = min(chain_dists, key=chain_dists.get)
+        leftover_chains = list(set(nonHOH_chains) - set(i))
+        # Remove remaining chains
+        for j in leftover_chains:
+            if j == i:
+                # Just in case.
+                continue
+            else:
+                temp_structure[0].remove_chain(j)
 
-            with open(o.replace('_mono.pdb', '.pdb'), 'w') as handle:
-                remark = ['REMARK warning: chains may be ommitted for alignment\n']
-                new_pdb = ''.join(remark + header_front + newfile_contents)
-                # print(new_pdb)
-                handle.write(new_pdb)
+        # Rename Chain to corresponding chain then save!
+        name = os.path.splitext(os.path.basename(f))[0] + '_' + str(i)
+        filename = os.path.join(self.outdir, f'{name}_mono.pdb')
+        temp_structure.write_pdb(filename)
+        filenames.append(filename)
 
-    def monomerize_single(self, file):
-        '''
-        Convert a pdb file into it's constituent chains
-        :param file: Input pdb filepath
-        :return: Formally nothing, but will create n number of monomerized pdbs per number of AA containing chains.
-        '''
-        maplist = self.get_maplist()
-        outnames = self.process_pdb(file, maplist)
-        print(outnames)
-        self.write_bound(file, outnames)
-        for o in outnames:
-            if os.path.isfile(o):
-                os.remove(o)
-
-    def monomerize_all(self):
-        '''
-        Converts a set of pdbs into constituent chains
-        :return: Formally nothing, but will separate pdb files according to chains.
-        '''
-        for f in self.get_filelist():
-            self.monomerize_single(file=f)
+    return filenames
 
 
 def get_chain_center(chain_name, file):
@@ -578,7 +599,6 @@ def get_chain_center(chain_name, file):
             continue
         else:
             structure.remove_chain(i)
-
     return structure.calculate_center_of_mass()
 
 
@@ -594,3 +614,147 @@ def find_water_chains(file):
     structure.remove_empty_chains()
     new = [x.name for x in structure[0]]
     return list(set(old) - set(new))
+
+
+def resample(
+        moving_xmap: Xmap,
+        transform: Transform,
+        reference_structure: Structure
+):
+    interpolated_grid = gridFromTemplate(moving_xmap)
+    mask = gemmi.FloatGrid(moving_xmap.xmap.nu,
+                           moving_xmap.xmap.nv,
+                           moving_xmap.xmap.nw, )
+    mask.set_unit_cell(moving_xmap.xmap.unit_cell)
+    mask.spacegroup = gemmi.find_spacegroup_by_name("P 1")
+
+    for model in reference_structure.structure:
+        for chain in model:
+            for residue in chain.get_polymer():
+                for atom in residue:
+                    mask.set_points_around(atom.pos, 5.0, 1.0)
+
+    mask_array = np.array(mask)
+    mask_indicies = np.hstack([x.reshape((len(x), 1))
+                              for x in np.nonzero(mask)])
+    fractional_coords = []
+    for model in reference_structure.structure:
+        for chain in model:
+            for residue in chain.get_polymer():
+                for atom in residue:
+                    fractional = moving_xmap.xmap.unit_cell.fractionalize(
+                        atom.pos)
+                    fractional_coords.append(
+                        [fractional.x, fractional.y, fractional.z])
+
+    fractional_coords_array = np.array(fractional_coords)
+    max_coord = np.max(fractional_coords_array, axis=0)
+    min_coord = np.min(fractional_coords_array, axis=0)
+
+    min_index = np.floor(
+        min_coord * np.array([interpolated_grid.nu, interpolated_grid.nv, interpolated_grid.nw]))
+    max_index = np.floor(
+        max_coord * np.array([interpolated_grid.nu, interpolated_grid.nv, interpolated_grid.nw]))
+
+    points = itertools.product(range(int(min_index[0]), int(max_index[0])),
+                               range(int(min_index[1]), int(max_index[1])),
+                               range(int(min_index[2]), int(max_index[2])),
+                               )
+
+    # Unpack the points, poitions and transforms
+    point_list: List[Tuple[int, int, int]] = []
+    position_list: List[Tuple[float, float, float]] = []
+    transform_list: List[gemmi.transform] = []
+    com_moving_list: List[np.array] = []
+    com_reference_list: List[np.array] = []
+
+    transform_rotate_reference_to_moving = transform.transform
+    transform_rotate_reference_to_moving.vec.fromlist([0.0, 0.0, 0.0])
+
+    transform_reference_to_centered = gemmi.Transform()
+    transform_reference_to_centered.vec.fromlist(
+        (-transform.com_reference).tolist())
+    transform_reference_to_centered.mat.fromlist(np.eye(3).tolist())
+
+    transform_centered_to_moving = gemmi.Transform()
+    transform_centered_to_moving.vec.fromlist(transform.com_moving.tolist())
+    transform_centered_to_moving.mat.fromlist(np.eye(3).tolist())
+
+    # indicies to positions
+    #   mask_array = np.array(mask)
+    # mask_indicies = np.hstack([x.reshape((len(x), 1))
+    #                          for x in np.nonzero(mask)])
+    #points2 = [x for x in points]
+    #ma, mb, mc = zip(*mask_indicies)
+    ##sma = set(ma)
+    #smb = set(mb)
+    #smc = set(mc)
+    #points_filtered = [(x, y, z) for x, y, z in points2 if x in sma and y in smb and z in smc]
+    for point in points:  # points_filtered:
+
+        position = interpolated_grid.point_to_position(
+            interpolated_grid.get_point(point[0], point[1], point[2]))
+
+        if mask.interpolate_value(position) <= 0:
+            continue
+        # Tranform to origin frame
+        trtc = transform_reference_to_centered.apply(position)
+        position_origin_reference = gemmi.Position(trtc[0], trtc[1], trtc[2])
+
+        # Rotate
+        trrtm = transform_rotate_reference_to_moving.apply(
+            position_origin_reference)
+        position_origin_moving = gemmi.Position(trrtm[0], trrtm[1], trrtm[2])
+
+        # Transform to moving frame
+        tctm = transform_centered_to_moving.apply(position_origin_moving)
+        position_moving = gemmi.Position(tctm[0], tctm[1], tctm[2])
+
+        # Interpolate moving map
+        interpolated_map_value = moving_xmap.xmap.interpolate_value(
+            position_moving)
+
+        # Set original point
+        interpolated_grid.set_value(
+            point[0], point[1], point[2], interpolated_map_value)
+
+    # interpolated_grid.symmetrize_max()
+    interpolated_array = np.array(interpolated_grid)
+
+    interpolated_grid_neg = gridFromTemplate(moving_xmap)
+    interpolated_array_neg = np.array(interpolated_grid_neg, copy=False)
+    interpolated_array_neg[:, :, :] = -interpolated_array[:, :, :]
+    interpolated_grid_neg.symmetrize_max()
+
+    interpolated_grid_pos = gridFromTemplate(moving_xmap)
+    interpolated_array_pos = np.array(interpolated_grid_pos, copy=False)
+    interpolated_array_pos[:, :, :] = interpolated_array[:, :, :]
+    interpolated_grid_pos.symmetrize_max()
+
+    interpolated_grid_sym = gridFromTemplate(moving_xmap)
+    interpolated_array_sym = np.array(interpolated_grid_sym, copy=False)
+    interpolated_array_sym[:, :, :] = interpolated_array_pos[:,
+                                                             :, :] - interpolated_array_neg[:, :, :]
+
+    return Xmap(interpolated_grid_sym)
+
+
+def referenceSave(template_map_path, xmap, path_to_save):
+    # Open Template map
+    ccp4 = gemmi.read_ccp4_map(str(template_map_path))
+    ccp4.setup()
+    # Replace Template map data
+    ccp4.grid = xmap.xmap
+    ccp4.setup()
+    # Save Template Map...
+    ccp4.write_ccp4_map(str(path_to_save))
+
+
+def gridFromTemplate(template):
+    interpolated_grid = gemmi.FloatGrid(
+        template.xmap.nu,
+        template.xmap.nv,
+        template.xmap.nw, )
+    interpolated_grid.set_unit_cell(template.xmap.unit_cell)
+    interpolated_grid.spacegroup = template.xmap.spacegroup
+    return interpolated_grid
diff --git a/fragalysis_api/xcimporter/conversion_pdb_mol.py b/fragalysis_api/xcimporter/conversion_pdb_mol.py
index da193e9c..6d4033fa 100755
--- a/fragalysis_api/xcimporter/conversion_pdb_mol.py
+++ b/fragalysis_api/xcimporter/conversion_pdb_mol.py
@@ -2,15 +2,17 @@
 import glob
 
 from rdkit import Chem
-from rdkit.Chem import AllChem
+from rdkit.Chem import AllChem, Draw
 from rdkit.Geometry import Point3D
 import json
 import os
 import shutil
 import warnings
 import csv
-import pypdb
+#import pypdb
 import numpy as np
+import gemmi
+import shutil
 
 
 class Ligand:
@@ -61,7 +63,7 @@ def remove_nonligands(self):
                 self.final_hets.append(line)
         return self.final_hets
 
-    def find_ligand_names_new(self):
+    def find_ligand_names_new(self, rrf=False):
         """
         Finds list of ligands contained in the structure, including solvents and ions
         :return: A listed of ligands that are not listed in the non_ligs.json file!
@@ -72,15 +74,15 @@ def find_ligand_names_new(self):
                 all_ligands.append(line)
 
         for lig in all_ligands:
-            if (
-                    lig.split()[3][-3:] not in self.non_ligs
-            ):  # this takes out the solvents and ions a.k.a non-ligands
+            if rrf:
+                # Take First Character?
+                # ? #lig.split()[4][0]: #??? This is breaking stuff??
+                if not os.path.basename(self.infile).rsplit('_', 2)[1] == lig[21]:
+                    continue
+            if lig[17:20].strip() not in self.non_ligs:
                 self.wanted_ligs.append(lig[16:20].strip() + lig[20:26])
-                # print(lig[16:20].strip() + lig[20:26])
 
         self.wanted_ligs = list(set(self.wanted_ligs))
-        # print(self.wanted_ligs)
-
         return self.wanted_ligs
 
     def get_3d_distance(self, coord_a, coord_b):
@@ -90,10 +92,11 @@ def get_3d_distance(self, coord_a, coord_b):
         :param coord_b: List of len 3 corresponding to xyz coords
         :return: The distance.
         '''
-        sum_ = (sum([(float(coord_a[i]) - float(coord_b[i])) ** 2 for i in range(3)]))
+        sum_ = (
+            sum([(float(coord_a[i]) - float(coord_b[i])) ** 2 for i in range(3)]))
         return np.sqrt(sum_)
 
-    def handle_covalent_mol(self, lig_res_name, non_cov_mol):
+    def handle_covalent_mol(self, lig_res_name, non_cov_mol, file_base):
         '''
         Do some magic if we think the molecule has a covalent attachment
         :param lig_res_name: Name of the covalent ligand
@@ -102,9 +105,12 @@ def handle_covalent_mol(self, lig_res_name, non_cov_mol):
         '''
         # original pdb = self.pdbfile (already aligned)
         # lig res name = name of ligand to find link for
+        fb = file_base.rsplit('_', 1)[1]
 
         covalent = False
 
+        chain = ''
+
         for line in self.pdbfile:
             if 'LINK' in line:
                 zero = line[13:27]
@@ -112,12 +118,19 @@ def handle_covalent_mol(self, lig_res_name, non_cov_mol):
 
                 if lig_res_name in zero:
                     res = one
+                    chain = one[8]
                     covalent = True
 
                 if lig_res_name in one:
                     res = zero
+                    chain = zero[8]
                     covalent = True
 
+        if(len(fb) > 1):
+            basechain = fb[-1]
+            if not chain == basechain:
+                return None
+
         if covalent:
             for line in self.pdbfile:
                 if 'ATOM' in line and line[13:27] == res:
@@ -125,26 +138,25 @@ def handle_covalent_mol(self, lig_res_name, non_cov_mol):
                     res_y = float(line[39:47])
                     res_z = float(line[47:55])
                     res_coords = [res_x, res_y, res_z]
-                    print(res_coords)
                     atm = Chem.MolFromPDBBlock(line)
                     atm_trans = atm.GetAtomWithIdx(0)
 
             orig_pdb_block = Chem.MolToPDBBlock(non_cov_mol)
 
-            lig_block = '\n'.join([l for l in orig_pdb_block.split('\n') if 'COMPND' not in l])
+            lig_block = '\n'.join(
+                [l for l in orig_pdb_block.split('\n') if 'COMPND' not in l])
             lig_lines = [l for l in lig_block.split('\n') if 'HETATM' in l]
             j = 0
             old_dist = 100
             for line in lig_lines:
                 j += 1
-                #                 print(line)
                 if 'HETATM' in line:
-                    coords = [line[31:39].strip(), line[39:47].strip(), line[47:55].strip()]
+                    coords = [line[31:39].strip(), line[39:47].strip(),
+                              line[47:55].strip()]
                     dist = self.get_3d_distance(coords, res_coords)
 
                     if dist < old_dist:
                         ind_to_add = j
-                        print(dist)
                         old_dist = dist
 
             i = non_cov_mol.GetNumAtoms()
@@ -153,7 +165,8 @@ def handle_covalent_mol(self, lig_res_name, non_cov_mol):
             edmol.AddBond(ind_to_add - 1, i, Chem.BondType.SINGLE)
             new_mol = edmol.GetMol()
             conf = new_mol.GetConformer()
-            conf.SetAtomPosition(i, Point3D(res_coords[0], res_coords[1], res_coords[2]))
+            conf.SetAtomPosition(i, Point3D(
+                res_coords[0], res_coords[1], res_coords[2]))
 
             return new_mol
 
@@ -167,84 +180,47 @@ def create_pdb_mol(self, file_base, lig_out_dir, smiles_file, handle_cov=False):
         :return: mol object that attempts to correct bond order if PDB entry
                 or mol object extracted from pdb file
         """
-        pdb_block = open(os.path.join(lig_out_dir, (file_base + ".pdb")), 'r').read()
-
-        lig_line = open(os.path.join(lig_out_dir, (file_base + ".pdb")), 'r').readline()
+        pdb_block = open(os.path.join(
+            lig_out_dir, (file_base + ".pdb")), 'r').read()
+        lig_line = open(os.path.join(
+            lig_out_dir, (file_base + ".pdb")), 'r').readline()
         res_name = lig_line[16:20].replace(' ', '')
 
-        # Look for PDB entries in PDB bank and use residue name to get bond order
-        if not smiles_file:
-            try:
-
-                mol = Chem.MolFromPDBBlock(pdb_block)
-                chem_desc = pypdb.describe_chemical(f"{res_name}")
-                new_smiles = chem_desc["describeHet"]["ligandInfo"]["ligand"]["smiles"]
-
-                template = Chem.MolFromSmiles(new_smiles)
-                new_mol = AllChem.AssignBondOrdersFromTemplate(template, mol)
-
-                if handle_cov:
-                    cov_mol = self.handle_covalent_mol(lig_res_name=res_name, non_cov_mol=new_mol)
-                    if cov_mol is not None:
-                        new_mol = cov_mol
-
-                return new_mol
-
-            except Exception as e:
-                new_pdb_block = ''
-
-                for lig in pdb_block.split('\n'):
-                    if 'ATM' in lig:
-                        pos = 16
-                        s = lig[:pos] + ' ' + lig[pos + 1:]
-                        new_pdb_block += s
-                    else:
-                        new_pdb_block += lig
-
-                    new_pdb_block += '\n'
-
-                mol = Chem.rdmolfiles.MolFromPDBBlock(new_pdb_block)
-
-                return mol
-
-        # Look for new XChem data - new XChem data must have associated smile.txt file
-        # Need to do this to catch corner case - x0685 from mArh residue
-        # name NHE was found ----> yielded wrong mol/smiles
-        if smiles_file:
-            new_pdb_block = ''
-
-            for lig in pdb_block.split('\n'):
-                if 'ATM' in lig:
-                    pos = 16
-                    s = lig[:pos] + ' ' + lig[pos + 1:]
-                    new_pdb_block += s
-                else:
-                    new_pdb_block += lig
-
-                new_pdb_block += '\n'
-
-            mol = Chem.rdmolfiles.MolFromPDBBlock(new_pdb_block)
-            if handle_cov:
-                cov_mol = self.handle_covalent_mol(lig_res_name=res_name, non_cov_mol=mol)
-                if cov_mol is not None:
-                    mol = cov_mol
-
-            return mol
-
-    def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent=False):
+        # Create new pdb_block and create a mol file regardles... (bond order gets assigned in create_mol_file...
+        new_pdb_block = ''
+        for lig in pdb_block.split('\n'):
+            if 'ATM' in lig:
+                pos = 16
+                s = lig[:pos] + ' ' + lig[pos + 1:]
+                new_pdb_block += s
+            else:
+                new_pdb_block += lig
+            new_pdb_block += '\n'
+        mol = Chem.rdmolfiles.MolFromPDBBlock(new_pdb_block)
+        # if not smiles_file:  # create a new template?
+        #    new_smiles = ''
+        #    template = Chem.MolFromSmiles(new_smiles)
+        #    mol = AllChem.AssignBondOrdersFromTemplate(template, mol)
+        if handle_cov:
+            cov_mol = self.handle_covalent_mol(
+                lig_res_name=res_name, non_cov_mol=mol, file_base=file_base)
+            if cov_mol is not None:
+                mol = cov_mol
+        return mol
+
+    def create_pdb_for_ligand(self, ligand, count, reduce, smiles_file, covalent=False):
         """
         A pdb file is produced for an individual ligand, containing atomic and connection information
         :param ligand: Name of the Ligand
         :param count: The index of the ligand
-        :param monomerize: Bool, if the file needs to be named using the chain name of the PDB
+        :param reduce: Bool, if the file needs to be named using the chain name of the PDB
         :param smiles_file: File path of smiles_file (if any)
         :param covalent: Bool, indicate whether or not covalent attach should be sought.
         :return: .pdb file for ligand.
         """
-
         # out directory and filename for lig pdb
         if not self.target_name in os.path.abspath(self.infile):
-            if not monomerize:
+            if not reduce:
                 file_base = str(
                     self.target_name
                     + "-"
@@ -255,7 +231,7 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
                     + "_"
                     + str(count)
                 )
-            if monomerize:
+            if reduce:
                 file_base = str(self.target_name
                                 + "-"
                                 + os.path.abspath(self.infile)
@@ -267,7 +243,7 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
                 file_base = file_base[:-2] + "_" + str(count) + chain
 
         else:
-            if not monomerize:
+            if not reduce:
                 file_base = str(
                     os.path.abspath(self.infile)
                     .split("/")[-1]
@@ -276,7 +252,7 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
                     + "_"
                     + str(count)
                 )
-            if monomerize:
+            if reduce:
                 file_base = str(
                     os.path.abspath(self.infile)
                     .split("/")[-1]
@@ -287,7 +263,6 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
                 file_base = file_base[:-2] + "_" + str(count) + chain
 
         lig_out_dir = os.path.join(self.RESULTS_DIRECTORY, file_base)
-
         individual_ligand = []
         individual_ligand_conect = []
         # adding atom information for each specific ligand to a list
@@ -299,10 +274,7 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
         for atom in individual_ligand:
             atom_number = atom.split()[1]
             for conection in self.conects:
-                if (
-                        atom_number in conection
-                        and conection not in individual_ligand_conect
-                ):
+                if (atom_number in conection and conection not in individual_ligand_conect):
                     individual_ligand_conect.append(conection)
                     con_num += 1
 
@@ -310,12 +282,14 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
         # (taking into account ligands that are covalently bound to the protein
 
         # assert 0 <= con_num - len(individual_ligand) <= 1
-
         # making into one list that is compatible with conversion to mol object
         ligand_het_con = individual_ligand + individual_ligand_conect
-
         # make a pdb file for the ligand molecule
 
+        if os.path.isdir(lig_out_dir):
+            # This is stupid but will correctly spec the files... is there a better solution??
+            shutil.rmtree(lig_out_dir)
+
         if not os.path.isdir(lig_out_dir):
             os.makedirs(lig_out_dir)
 
@@ -325,14 +299,14 @@ def create_pdb_for_ligand(self, ligand, count, monomerize, smiles_file, covalent
         for line in ligand_het_con:
             ligands_connections.write(str(line))
         ligands_connections.close()
-
         # making pdb file into mol object
-        mol = self.create_pdb_mol(file_base=file_base, lig_out_dir=lig_out_dir, smiles_file=smiles_file, handle_cov=covalent)
-
+        mol = self.create_pdb_mol(
+            file_base=file_base, lig_out_dir=lig_out_dir, smiles_file=smiles_file, handle_cov=covalent)
         # Move Map files into lig_out_dir
 
         if not mol:
-            print(f'WARNING: {file_base} did not produce a mol object from its pdb lig file!')
+            print(
+                f'WARNING: {file_base} did not produce a mol object from its pdb lig file!')
         else:
             try:
                 Chem.AddHs(mol)
@@ -365,20 +339,26 @@ def create_mol_file(self, directory, file_base, mol_obj, smiles_file=None):
             try:
                 smiles = open(smiles_file, 'r').readlines()[0].rstrip()
                 template = AllChem.MolFromSmiles(smiles)
-                new_mol = AllChem.AssignBondOrdersFromTemplate(template, mol_obj)
-
-                return Chem.rdmolfiles.MolToMolFile(new_mol, out_file)
+                mol_obj = AllChem.AssignBondOrdersFromTemplate(
+                    template, mol_obj)
+                #Draw.MolToFile(new_mol, out_png)
+                # return Chem.rdmolfiles.MolToMolFile(new_mol, out_file)
             except Exception as e:
                 print(e)
                 print('failed to fit template ' + smiles_file)
                 print(f'template smiles: {smiles}')
-                return Chem.rdmolfiles.MolToMolFile(mol_obj, out_file)
+                #Draw.MolToFile(mol_obj, out_png)
+                # return Chem.rdmolfiles.MolToMolFile(mol_obj, out_file)
 
         else:
             print(f'Warning: No smiles file: {file_base}')
 
-        # creating mol file
-        return Chem.rdmolfiles.MolToMolFile(mol_obj, out_file)
+        Chem.rdmolfiles.MolToMolFile(mol_obj, out_file)
+        # IS this enough... do we lost stuff?
+        out_png = os.path.join(directory, str(file_base + ".png"))
+        draw_mol = Chem.Mol(mol_obj)
+        AllChem.Compute2DCoords(draw_mol)
+        Draw.MolToFile(draw_mol, out_png, imageType='png')
 
     def create_sd_file(self, mol_obj, writer):
         """
@@ -398,7 +378,8 @@ def create_metadata_file(self, directory, file_base, mol_obj, smiles_file=None):
         """
 
         meta_out_file = os.path.join(directory, str(file_base + "_meta.csv"))
-        smiles_out_file = os.path.join(directory, str(file_base + "_smiles.txt"))
+        smiles_out_file = os.path.join(
+            directory, str(file_base + "_smiles.txt"))
 
         if smiles_file:
             try:
@@ -460,7 +441,8 @@ def make_apo_file(self, keep_headers=False):
         if keep_headers:
             include = ['CONECT', 'SEQRES', 'TITLE', 'ANISOU']
         else:
-            include = ['CONECT', 'REMARK', 'CRYST', 'SEQRES', 'HEADER', 'TITLE', 'ANISOU']
+            include = ['CONECT', 'REMARK', 'CRYST',
+                       'SEQRES', 'HEADER', 'TITLE', 'ANISOU']
 
         for line in self.pdbfile:
             if (
@@ -473,7 +455,8 @@ def make_apo_file(self, keep_headers=False):
                 lines += line
 
         apo_file = open(
-            os.path.join(self.RESULTS_DIRECTORY, str(self.filebase + "_apo.pdb")), "w+"
+            os.path.join(self.RESULTS_DIRECTORY, str(
+                self.filebase + "_apo.pdb")), "w+"
         )
         apo_file.write(str(lines))
         apo_file.close()
@@ -491,22 +474,24 @@ def add_biomol_remark(self):
         Add contents of biomol/additional text file to a .pdb file
         '''
         biomol_remark = open(self.biomol).readlines()
-        print(biomol_remark)
         f = self.apo_file
         with open(f) as handle:
             switch = 0
             header_front, header_end = [], []
             pdb = []
             for line in handle:
-                if line.startswith('ATOM'): switch = 1
-                if line.startswith('HETATM'): switch = 2
+                if line.startswith('ATOM'):
+                    switch = 1
+                if line.startswith('HETATM'):
+                    switch = 2
                 if switch == 0:
                     header_front.append(line)
                 elif (switch == 2) and not line.startswith('HETATM'):
                     header_end.append(line)
                 else:
                     pdb.append(line)
-            full_file = ''.join(header_front) + ''.join(biomol_remark) + ''.join(pdb) + ''.join(header_end)
+            full_file = ''.join(
+                header_front) + ''.join(biomol_remark) + ''.join(pdb) + ''.join(header_end)
             with open(f, 'w') as w:
                 w.write(full_file)
 
@@ -525,7 +510,8 @@ def make_apo_desol_files(self):
             "w+",
         )
         solv_file = open(
-            os.path.join(self.RESULTS_DIRECTORY, str(self.filebase + "_apo-solv.pdb")),
+            os.path.join(self.RESULTS_DIRECTORY, str(
+                self.filebase + "_apo-solv.pdb")),
             "w+",
         )
         if not self.apo_file:
@@ -542,40 +528,35 @@ def make_apo_desol_files(self):
         prot_file.close()
 
 
-def set_up(target_name, infile, out_dir, monomerize, smiles_file=None, biomol=None, covalent=False, keep_headers=False):
+def set_up(target_name, infile, out_dir, rrf, smiles_file=None, biomol=None, covalent=False, keep_headers=False):
     """
     For each ligand inside a pdb file, process each ligand seperately and create own outputs in individual folders.
     :param target_name: Name of the folder in out_dir
     :param infile: pdb file to be processed
     :param out_dir: Overall location of outputs
-    :param monomerize: Bool, indicate whether or not data was monomerized prior to set_up
+    :param rrf: Bool, indicate whether or not data is to be set to a single reference frame
     :param smiles_file: Filepath pointing to text file containing smiles string (if exists)
     :param biomol: Filepath pointing to text file containing biomol/header information for pdbs (if exists)
     :param covalent: Bool, indicate whether or not output mol files should find covalent attachment.
     :param keep_headers: Bool, indicate whether or not keep headers on apo files.
     :return: for each ligand: pdb, mol, sdf and _apo.pdb in seperate directorys inside out_dir/target_name
     """
-
     RESULTS_DIRECTORY = os.path.join(out_dir, target_name, 'aligned')
-
     if not os.path.isdir(RESULTS_DIRECTORY):
         os.makedirs(RESULTS_DIRECTORY)
 
-    print(RESULTS_DIRECTORY)
-
     new = Ligand(
         target_name, infile, RESULTS_DIRECTORY
     )  # takes in pdb file and returns specific ligand files
     new.hets_and_cons()  # takes only hetatm and conect file lines from pdb file
     new.remove_nonligands()  # removes ions and solvents from list of ligands
-    new.find_ligand_names_new()  # finds the specific name and locations of desired ligands
+    # finds the specific name and locations of desired ligands
+    new.find_ligand_names_new(rrf=rrf)
     for i in range(len(new.wanted_ligs)):
         new.create_pdb_for_ligand(
-            new.wanted_ligs[i], count=i, monomerize=monomerize, smiles_file=smiles_file, covalent=covalent
+            new.wanted_ligs[i], count=i, reduce=rrf, smiles_file=smiles_file, covalent=covalent
         )  # creates pdb file and mol object for specific ligand
-
     for i in range(len(new.mol_dict["directory"])):
-
         if not new.mol_dict["mol"][i]:
             warnings.warn(
                 str(
@@ -585,7 +566,6 @@ def set_up(target_name, infile, out_dir, monomerize, smiles_file=None, biomol=No
                 )
             )
             continue
-
         shutil.copy(infile,
                     os.path.join(new.mol_dict["directory"][i], str(new.mol_dict["file_base"][i] + "_bound.pdb")))
 
@@ -593,27 +573,26 @@ def set_up(target_name, infile, out_dir, monomerize, smiles_file=None, biomol=No
         basebase = os.path.basename(inpath)
         fofcmap_files = glob.glob(f'{inpath}_*.map')
         event_files = glob.glob(f'{inpath}_*.ccp4')
-        map_files = fofcmap_files + event_files
-        for map_file in map_files:
-            map_base = os.path.basename(map_file)
-            map_base = map_base.replace(basebase, new.mol_dict["file_base"][i])
-            shutil.copy(map_file,
-                        os.path.join(new.mol_dict["directory"][i], map_base))
-
+        json_files = glob.glob(f'{inpath}_*.json')
+        other_files = fofcmap_files + event_files + json_files
+        for other_file in other_files:
+            other_base = os.path.basename(other_file)
+            other_base = other_base.replace(
+                basebase, new.mol_dict["file_base"][i])
+            shutil.copy(other_file,
+                        os.path.join(new.mol_dict["directory"][i], other_base))
         new.create_mol_file(
             directory=new.mol_dict["directory"][i],
             file_base=new.mol_dict["file_base"][i],
             mol_obj=new.mol_dict["mol"][i],
             smiles_file=smiles_file,
         )  # creates mol file for each ligand
-
         writer = Chem.rdmolfiles.SDWriter(
             os.path.join(
                 new.mol_dict["directory"][i],
                 str(new.mol_dict["file_base"][i] + ".sdf"),
             )
         )
-
         new.create_sd_file(
             new.mol_lst[i], writer
         )  # creates sd file containing all mol files
@@ -625,7 +604,6 @@ def set_up(target_name, infile, out_dir, monomerize, smiles_file=None, biomol=No
             file_base=new.mol_dict["file_base"][i],
             smiles_file=smiles_file,
         )  # create metadata csv file for each ligand
-
         new_apo = pdb_apo(
             infile,
             target_name,
@@ -633,7 +611,49 @@ def set_up(target_name, infile, out_dir, monomerize, smiles_file=None, biomol=No
             new.mol_dict["file_base"][i],
             biomol=biomol
         )
-        new_apo.make_apo_file(keep_headers=keep_headers)  # creates pdb file that doesn't contain any ligand information
-        new_apo.make_apo_desol_files()  # makes apo file without solvent, ions and buffers, and file with just those
+        # creates pdb file that doesn't contain any ligand information
+        new_apo.make_apo_file(keep_headers=keep_headers)
+        # makes apo file without solvent, ions and buffers, and file with just those
+        new_apo.make_apo_desol_files()
 
     return new
+
+
+def convert_small_AA_chains(in_file, out_file, max_len=15):
+    pdb_file = gemmi.read_structure(in_file)
+    structure = pdb_file[0]
+    chain_lens = [j.get_polymer().length() for j in structure]
+    for i, j in enumerate(chain_lens):
+        if int(j) <= int(max_len):
+            for x in structure[i]:
+                x.name = 'LIG'
+                x.het_flag = 'H'
+                x.seqid.num = i + 1
+    w_chain = structure.find_chain('W')
+    nhoh = 1
+    if w_chain is not None:
+        newchain = gemmi.Chain('W')
+        for res in w_chain:
+            if res.name == 'LIG' or res.name == 'HOH':  # Should be a ligand in the water chain right??
+                for atom in res:
+                    emptyres = gemmi.Residue()
+                    emptyres.add_atom(atom)
+                    emptyres.name = 'HOH'
+                    emptyres.seqid.num = nhoh
+                    nhoh += 1
+                    newchain.add_residue(emptyres)
+            else:
+                newchain.add_residues(res)
+        structure.remove_chain('W')
+        structure.add_chain(newchain)
+    pdb_file.write_pdb(out_file)
+
+
+def copy_extra_files(in_file, out_dir):
+    bn = os.path.basename(in_file).replace('.pdb', '')
+    dirname = os.path.dirname(in_file)
+    other_files = set(glob.glob(os.path.join(
+        dirname, f'{bn}*'))) - set([in_file])
+    new_loc = [os.path.join(out_dir, os.path.basename(x)) for x in other_files]
+    for i, j in zip(other_files, new_loc):
+        shutil.copyfile(i, j)
diff --git a/fragalysis_api/xcimporter/non_ligs.json b/fragalysis_api/xcimporter/non_ligs.json
index 45d16f27..8121895f 100755
--- a/fragalysis_api/xcimporter/non_ligs.json
+++ b/fragalysis_api/xcimporter/non_ligs.json
@@ -1,1 +1,1 @@
-["PJE","DMS", "CL", "MG", "EDO", "PEG", "MLI", "NI", "SO4", "WAT", "LI", "NA", "K", "RB", "BE", "MG", "CA", "SR", "F", "CL", "BR", "I", "MN", "FE", "CO", "NI", "CU", "ZN", "3CO", "BUQ", "NAG", "NAD", "CR", "SF4", "EOH", "ZNO", "NAO", "EOM", "EHN", "NAP", "CCN", "NAW", "BR", "EGL", "NI2", "GSH", "NI1", "O2", "BA", "RU", "SAH", "GU7", "SAM", "TAS", "DCE", "2BM", "TP7", "OF3", "OF1", "RB", "IOH", "MW1", "IOD", "C2O", "BNZ", "TCN", "ARS", "NH4", "GD", "PER", "GA", "TPP", "CHX", "CME", "THB", "IPA", "CD1", "OH", "SO4", "DTT", "PQN", "CYN", "PQQ", "PYJ", "PEO", "NA6", "MBR", "NA5", "OS", "MAN", "CMO", "OCL", "DMF", "OCN", "MO3", "NGN", "ACT", "U1", "HDZ", "MO5", "MO4", "VO4", "DMS", "FUC", "PCL", "CB5", "EEE", "HG", "NO2", "CMP", "PR", "BMA", "IUM", "PT", "ZN2", "TTP", "NO3", "YT3", "TYS", "PB", "M2M", "ZO3", "PD", "AMP", "PI", "MH3", "AF3", "ZN", "MN3", "OXY", "NI", "CSD", "OX", "PS5", "MN5", "MN6", "S", "W", "SB", "FOL", "OXE", "PT4", "PBM", "O", "MW2", "MG", "543", "MSM", "C5P", "ANL", "MTO", "NO", "TBU", "OPY", "PC4", "GU3", "GU2", "GU1", "MOH", "ANP", "GU6", "GU5", "GU4", "AU", "OC3", "BTN", "I42", "OC4", "OC7", "OC6", "TMP", "RE", "GD3", "CTP", "ACE", "3OF", "ETZ", "MM4", "IN", "ACN", "DOD", "AST", "COA", "EU", "DOX", "COB", "B12", "REO", "ATP", "CD3", "U10", "ACY", "PEG", "YB", "NDP", "NBZ", "ETI", "SER", "C2C", "NA", "FMT", "ASC", "AU3", "FE2", "LNK", "SEK", "MO1", "EU3", "1BO", "AUC", "CLO", "FE", "DUM", "ADP", "OF2", "BEF", "FEL", "BF4", "HEX", "CUZ", "NDG", "XE", "FMN", "YAN", "CUA", "V", "CUO", "HEM", "GMP", "CU", "MGF", "GDP", "CFT", "SBT", "PLP", "SR", "FU1", "EDN", "EDO", "H2S", "ND4", "BRO", "KR", "CS", "NME", "CDP", "HGI", "SM", "ALY", "NMO", "TDP", "SE", "HO", "3CN", "AZI", "F42", "FLO", "6MO", "EMC", "Y1", "MO7", "SE4", "BF2", "CO", "NGD", "2MO", "202", "DIS", "MBN", "LA", "PGO", "CL", "HP6", "SO2", "LI", "PPS", "TPO", "POL", "GU0", "SGM", "DTU", "MOO", "TE", "TB", "CA", "FAD", "CNV", "GOL", "SCN", "AG", "PO4", "IR", "DIO", "NH2", "8CL", "3NI", "IRI", "UTP", "AR", "N4M", "CE", "NH3", "MN", "CNN", "HGC", "GU8", "GTP", "UDP", "OC2", "ART", "TFH", "MCH", "2NO", "6WO", "CD5", "KCX", "E1H", "ARF", "TL", "DXE", "GU9", "IDO", "KO4", "NRU", "4MO", "HOH", "010"]
\ No newline at end of file
+["F3F", "PJE","DMS", "CL", "MG", "EDO", "PEG", "MLI", "NI", "SO4", "WAT", "LI", "NA", "K", "RB", "BE", "MG", "CA", "SR", "F", "CL", "BR", "I", "MN", "FE", "CO", "NI", "CU", "ZN", "3CO", "BUQ", "NAG", "NAD", "CR", "SF4", "EOH", "ZNO", "NAO", "EOM", "EHN", "NAP", "CCN", "NAW", "BR", "EGL", "NI2", "GSH", "NI1", "O2", "BA", "RU", "SAH", "GU7", "SAM", "TAS", "DCE", "2BM", "TP7", "OF3", "OF1", "RB", "IOH", "MW1", "IOD", "C2O", "BNZ", "TCN", "ARS", "NH4", "GD", "PER", "GA", "TPP", "CHX", "CME", "THB", "IPA", "CD1", "OH", "SO4", "DTT", "PQN", "CYN", "PQQ", "PYJ", "PEO", "NA6", "MBR", "NA5", "OS", "MAN", "CMO", "OCL", "DMF", "OCN", "MO3", "NGN", "ACT", "U1", "HDZ", "MO5", "MO4", "VO4", "DMS", "FUC", "PCL", "CB5", "EEE", "HG", "NO2", "CMP", "PR", "BMA", "IUM", "PT", "ZN2", "TTP", "NO3", "YT3", "TYS", "PB", "M2M", "ZO3", "PD", "AMP", "PI", "MH3", "AF3", "ZN", "MN3", "OXY", "NI", "CSD", "OX", "PS5", "MN5", "MN6", "S", "W", "SB", "FOL", "OXE", "PT4", "PBM", "O", "MW2", "MG", "543", "MSM", "C5P", "ANL", "MTO", "NO", "TBU", "OPY", "PC4", "GU3", "GU2", "GU1", "MOH", "ANP", "GU6", "GU5", "GU4", "AU", "OC3", "BTN", "I42", "OC4", "OC7", "OC6", "TMP", "RE", "GD3", "CTP", "ACE", "3OF", "ETZ", "MM4", "IN", "ACN", "DOD", "AST", "COA", "EU", "DOX", "COB", "B12", "REO", "ATP", "CD3", "U10", "ACY", "PEG", "YB", "NDP", "NBZ", "ETI", "SER", "C2C", "NA", "FMT", "ASC", "AU3", "FE2", "LNK", "SEK", "MO1", "EU3", "1BO", "AUC", "CLO", "FE", "DUM", "ADP", "OF2", "BEF", "FEL", "BF4", "HEX", "CUZ", "NDG", "XE", "FMN", "YAN", "CUA", "V", "CUO", "HEM", "GMP", "CU", "MGF", "GDP", "CFT", "SBT", "PLP", "SR", "FU1", "EDN", "EDO", "H2S", "ND4", "BRO", "KR", "CS", "NME", "CDP", "HGI", "SM", "ALY", "NMO", "TDP", "SE", "HO", "3CN", "AZI", "F42", "FLO", "6MO", "EMC", "Y1", "MO7", "SE4", "BF2", "CO", "NGD", "2MO", "202", "DIS", "MBN", "LA", "PGO", "CL", "HP6", "SO2", "LI", "PPS", "TPO", "POL", "GU0", "SGM", "DTU", "MOO", "TE", "TB", "CA", "FAD", "CNV", "GOL", "SCN", "AG", "PO4", "IR", "DIO", "NH2", "8CL", "3NI", "IRI", "UTP", "AR", "N4M", "CE", "NH3", "MN", "CNN", "HGC", "GU8", "GTP", "UDP", "OC2", "ART", "TFH", "MCH", "2NO", "6WO", "CD5", "KCX", "E1H", "ARF", "TL", "DXE", "GU9", "IDO", "KO4", "NRU", "4MO", "HOH", "010"]
diff --git a/fragalysis_api/xcimporter/single_import.py b/fragalysis_api/xcimporter/single_import.py
index 383a044e..6548386c 100644
--- a/fragalysis_api/xcimporter/single_import.py
+++ b/fragalysis_api/xcimporter/single_import.py
@@ -3,10 +3,10 @@
 import os
 import shutil
 
-from fragalysis_api import Align, Monomerize, set_up
+from fragalysis_api import Align, set_up, convert_small_AA_chains, copy_extra_files, Sites, contextualize_crystal_ligands
 
 
-def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=None, covalent=False, self_ref=False):
+def import_single_file(in_file, out_dir, target, reduce_reference_frame, reference_pdb, biomol=None, covalent=False, self_ref=False, max_lig_len=0):
     '''Formats a PDB file into fragalysis friendly format.
     1. Validates the naming of the pdbs.
     2. It aligns the pdbs (_bound.pdb file).
@@ -16,12 +16,13 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
     :param in_file: Filepath of pdb you with to align (where additional files are stored in same directory)
     :param out_dir: Directory containing processed pdbs (will be created if it doesn't exists)
     :param target: Name of the folder to be created inside out_dir
-    :param monomerize: Bool, if True, will attempt to split pdb files into seperate chains
-    :param reference: Name of the Reference pdb to align in_file to. If called from command line this will default to out_dir/target/reference.pdb
+    :param reduce_reference_frame: Bool, if True, will attempt to split pdb files into seperate chains
+    :param reference_pdb: Name of the Reference pdb to align in_file to. If called from command line this will default to out_dir/target/reference.pdb
     :param biomol: plain-text file containing header information about the bio-molecular
         context of the pdb structures. If provided the contents will be appended to the top of the _apo.pdb files
     :param covalent: Bool, if True, will attempt to convert output .mol files to account for potential covalent attachments
     :param self_ref: Bool, if True, the import single file will align to itself.
+    :max_lig_len: Integer, If >0 will convert all chains with fewer than max_lig_len residues to HETATM with the name LIG. [Currently broken, yikes]
     :return: Hopefully, beautifully aligned files that be used with the fragalysis loader :)
     '''
 
@@ -35,49 +36,49 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
 
     pdb_smiles_dict = {'pdb': [], 'smiles': []}
 
-    if monomerize:
-        print("Monomerizing input structure")
-        out = os.path.join(out_dir, f'mono{target}/')
+    # Experimental - if option is used then chains are converted.
+    if int(max_lig_len) > int(0):
+        print(
+            f'EXPERIMENTAL: Converting all chains with less than {max_lig_len} residues to HETATM LIG')
+        out = os.path.join(out_dir, f'maxliglen{target}/')
         if not os.path.isdir(out):
             os.makedirs(out)
-
-        mono = Monomerize(directory=in_dir, outdir=out)
-        mono.monomerize_single(file=in_file)
+        convert_small_AA_chains(in_file=in_file, out_file=os.path.join(
+            out, in_file), max_len=max_lig_len)
+        copy_extra_files(in_file=in_file, out_dir=out)
         in_dir = out
-        for f in os.listdir(in_dir):
-            if '.pdb' in f:
-                pdb_smiles_dict['pdb'].append(os.path.join(in_dir, f))
-                if os.path.isfile(os.path.join(in_dir, f).replace('.pdb', '_smiles.txt')):
-                    pdb_smiles_dict['smiles'].append(os.path.join(in_dir, f).replace('.pdb', '_smiles.txt'))
-                else:
-                    pdb_smiles_dict['smiles'].append(None)
+
+    pdb_smiles_dict['pdb'].append(in_file)
+    if os.path.isfile(in_file.replace('.pdb', '_smiles.txt')):
+        pdb_smiles_dict['smiles'].append(
+            in_file.replace('.pdb', '_smiles.txt'))
     else:
-        pdb_smiles_dict['pdb'].append(in_file)
-        if os.path.isfile(in_file.replace('.pdb', '_smiles.txt')):
-            pdb_smiles_dict['smiles'].append(in_file.replace('.pdb', '_smiles.txt'))
-        else:
-            pdb_smiles_dict['smiles'].append(None)
+        pdb_smiles_dict['smiles'].append(None)
 
-    print("Aligning to Reference")
-    structure = Align(in_dir, "", monomerize)
+    print(f"Aligning to Reference: {reference_pdb}")
+    structure = Align(in_dir, "", rrf=reduce_reference_frame,
+                      refset=False)
 
     for i in pdb_smiles_dict['pdb']:
         if self_ref:
-            structure.align_to_reference(i, i, out_dir=os.path.join(out_dir, f"tmp{target}"))
+            structure.align_to_reference(
+                i, i, out_dir=os.path.join(out_dir, f"tmp{target}"))
         else:
-            structure.align_to_reference(i, reference, out_dir=os.path.join(out_dir, f"tmp{target}"))
+            structure.align_to_reference(
+                i, reference_pdb=reference_pdb, out_dir=os.path.join(out_dir, f"tmp{target}"))
 
     for smiles_file in pdb_smiles_dict['smiles']:
         if smiles_file:
-            print(smiles_file)
-            shutil.copyfile(smiles_file, os.path.join(os.path.join(out_dir, f"tmp{target}", smiles_file.split('/')[-1])))
-            print(os.path.join(out_dir, f"tmp{target}", smiles_file.split('/')[-1]))
+            shutil.copyfile(smiles_file, os.path.join(os.path.join(
+                out_dir, f"tmp{target}", smiles_file.split('/')[-1])))
+            print(os.path.join(
+                out_dir, f"tmp{target}", smiles_file.split('/')[-1]))
 
     aligned_dict = {'bound_pdb': [], 'smiles': []}
-
     for f in os.listdir(os.path.join(out_dir, f"tmp{target}")):
         if '.pdb' in f:
-            aligned_dict['bound_pdb'].append(os.path.join(out_dir, f"tmp{target}", f))
+            aligned_dict['bound_pdb'].append(
+                os.path.join(out_dir, f"tmp{target}", f))
             if os.path.isfile(os.path.join(out_dir, f"tmp{target}", f).replace('_bound.pdb', '_smiles.txt')):
                 aligned_dict['smiles'].append(os.path.join(out_dir, f"tmp{target}", f).replace('_bound.pdb',
                                                                                                '_smiles.txt'))
@@ -90,20 +91,20 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
                 _ = set_up(target_name=target,
                            infile=os.path.abspath(aligned),
                            out_dir=out_dir,
-                           monomerize=monomerize,
+                           rrf=reduce_reference_frame,
                            smiles_file=os.path.abspath(smiles),
                            biomol=biomol,
                            covalent=covalent,
-                           keep_headers=self_ref)
+                           keep_headers=True)
 
             else:
                 _ = set_up(target_name=target,
                            infile=os.path.abspath(aligned),
                            out_dir=out_dir,
-                           monomerize=monomerize,
+                           rrf=reduce_reference_frame,
                            biomol=biomol,
                            covalent=covalent,
-                           keep_headers=self_ref)
+                           keep_headers=True)
 
         except AssertionError:
             print(aligned, "is not suitable, please consider removal or editing")
@@ -115,15 +116,14 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
     dest_dir = os.path.join(out_dir, target, 'crystallographic')
     globstr = in_file.replace('.pdb', '*')
     for file in glob.glob(globstr):
-        print(file)
         shutil.copy(file, dest_dir)
 
-    if os.path.exists(os.path.join(out_dir, f'mono{target}')):
-        shutil.rmtree(os.path.join(out_dir, f'mono{target}'))
-
-    if os.path.exists(os.path.join(out_dir, f'tmp{target}')):
-        shutil.rmtree(os.path.join(out_dir, f'tmp{target}'))
+    # Time to use a for loop?
+    clean_up = [os.path.join(out_dir, f'maxliglen{target}'), os.path.join(
+        out_dir, f'mono{target}'), os.path.join(out_dir, f'tmp{target}')]
+    [shutil.rmtree(x) for x in clean_up if os.path.exists(x)]
 
+    # Cut Maps
     print("Files are now in a fragalysis friendly format!")
 
 
@@ -144,14 +144,17 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
         required=True,
     )
     parser.add_argument(
-        "-m", "--monomerize", action="store_true", default=False, help="Monomerize input"
+        "-rrf", "--reduce_reference_frame", action="store_true", default=False, help="Reduce Reference Frame"
     )
 
     parser.add_argument("-t", "--target", help="Target name", required=True)
 
-    parser.add_argument('-r', '--reference', help='Reference Structure', required=False, default=None)
-    parser.add_argument('-sr', '--self_reference', action="store_true", help='Include this flag if you want each pdb file to only align to itself', required=False, default=False)
-    parser.add_argument("-b", "--biomol_txt", help="Biomol Input txt file", required=False, default=None)
+    parser.add_argument('-r', '--reference_pdb',
+                        help='Reference PDB Structure', required=False, default=None)
+    parser.add_argument('-sr', '--self_reference', action="store_true",
+                        help='Include this flag if you want each pdb file to only align to itself', required=False, default=False)
+    parser.add_argument(
+        "-b", "--biomol_txt", help="Biomol Input txt file", required=False, default=None)
     parser.add_argument("-c",
                         "--covalent",
                         action="store_true",
@@ -159,39 +162,81 @@ def import_single_file(in_file, out_dir, target, monomerize, reference, biomol=N
                         required=False,
                         default=False
                         )
+    parser.add_argument("-mll",
+                        "--max_lig_len",
+                        help="Int, Convert all chains shorter than max_lig_len to HETATM LIG",
+                        required=False,
+                        default=0)
+
+    parser.add_argument(
+        "-cs",
+        "--cluster_sites",
+        action="store_true",
+        help='Include this flag if you would like to automatically assign center of mass sites as site labels',
+        required=False,
+        default=False
+    )
+
+    parser.add_argument(
+        "-cs_com",
+        "--cluster_sites_com",
+        help="Tolerance value for creating new clusters for centre of mass sites",
+        type=float,
+        default=5.00
+    )
+
+    parser.add_argument(
+        "-cs_other",
+        "--cluster_sites_other",
+        help="Tolerance value for creating new clusters for non centre of mass sites",
+        type=float,
+        default=1.00
+    )
 
     args = vars(parser.parse_args())
 
     in_file = args["in_file"]
     out_dir = args["out_dir"]
-    monomerize = args["monomerize"]
+    reduce_reference_frame = args["reduce_reference_frame"]
     target = args["target"]
     biomol = args["biomol_txt"]
     covalent = args["covalent"]
     self_ref = args['self_reference']
+    mll = args['max_lig_len']
+    cs = args['cluster_sites']
+    cs_com = args['cluster_sites_com']
+    cs_other = args['cluster_sites_other']
 
     # Will this work?
     if self_ref:
-        reference = in_file
-    elif args['reference'] is None:
-        reference = os.path.join(out_dir, target, 'reference.pdb')
+        reference_pdb = in_file
+    elif args['reference_pdb'] is None:
+        reference_pdb = os.path.join(out_dir, target, 'reference.pdb')
     else:
-        reference = args['reference']
+        reference_pdb = args['reference_pdb']
 
     if out_dir == os.path.join("..", "..", "data", "xcimporter", "output"):
         print("Using the default input directory ", out_dir)
 
-    if not os.path.isfile(reference):
-        print(f'Cannot find file called {reference}, please make sure the path is correct (or specify another reference using -r)!')
+    if not os.path.isfile(reference_pdb):
+        print(
+            f'Cannot find file called {reference_pdb}, please make sure the path is correct (or specify another reference using -r)!')
     else:
         import_single_file(in_file=in_file,
                            out_dir=out_dir,
                            target=target,
-                           monomerize=monomerize,
-                           reference=reference,
+                           reduce_reference_frame=reduce_reference_frame,
+                           reference_pdb=reference_pdb,
                            biomol=biomol,
                            covalent=covalent,
-                           self_ref=self_ref)
-        print(f'File has been aligned to {reference}')
-
-
+                           self_ref=self_ref,
+                           max_lig_len=mll)
+        print(f'File has been aligned to {reference_pdb}')
+        if cs:
+            folder = os.path.join(out_dir, target)
+            site_obj = Sites.from_folder(folder, recalculate=False)
+            site_obj.cluster_missing_mols(
+                com_tolerance=cs_com, other_tolerance=cs_other)
+            site_obj.to_json()
+            contextualize_crystal_ligands(folder=folder)
+            site_obj.apply_to_metadata()
diff --git a/fragalysis_api/xcimporter/sites.py b/fragalysis_api/xcimporter/sites.py
new file mode 100644
index 00000000..9c02f822
--- /dev/null
+++ b/fragalysis_api/xcimporter/sites.py
@@ -0,0 +1,353 @@
+from rdkit import Chem, DataStructs
+from rdkit.Chem import Descriptors
+import numpy as np
+import glob
+import os
+from frag.alysis.run_clustering import run_lig_cluster
+import json
+import argparse
+import dataclasses
+import csv
+
+# This needs refactoring.
+
+
+def get_distance(coord_a, coord_b):
+    sum_ = (sum([(float(coord_a[i])-float(coord_b[i]))**2 for i in range(3)]))
+    return np.sqrt(sum_)
+
+
+def json_file_generator(folder):
+    jsons = glob.glob(os.path.join(folder, 'aligned', '*', '*_sites.json'))
+    if len(jsons) == 0:
+        jsons = glob.glob(os.path.join(folder, '*', '*_sites.json'))
+    if len(jsons) == 0:
+        raise RuntimeError(f'Cannot find any .json files in {folder}')
+    return jsons
+
+
+def mol_file_generator(folder):
+    mols = glob.glob(os.path.join(folder, 'aligned', '*', '*.mol'))
+    if len(mols) == 0:
+        mols = glob.glob(os.path.join(folder, '*', '*.mol'))
+    if len(mols) == 0:
+        raise RuntimeError(f'Cannot find any .mol files in {folder}')
+    return mols
+
+
+def read_jfile(j_file):
+    with open(j_file) as jfile:
+        data = json.load(jfile)
+    data['fn'] = j_file.replace('_sites.json', '.mol')  # Needed?
+    return data
+
+
+def recluster_from_jsons(folder):
+    # Fetch list of all jsons
+    jsons = [read_jfile(x) for x in json_file_generator(folder)]
+    # Curate dictionary of sites from that
+    out_data = {}
+    for jf in jsons:
+        for keys in jf.keys():
+            if keys == 'fn':
+                continue
+            if out_data.get(keys) is None:
+                out_data[keys] = {}
+            for n, id in enumerate(jf.get(keys).get('site_id')):
+                if out_data[keys].get(id) is None:
+                    out_data[keys][id] = {
+                        'centre_of_mass': jf.get(keys).get('site_centre_of_mass')[n],
+                        'mol_ids': [jf.get('fn')]
+                    }
+                else:
+                    out_data[keys][id]['mol_ids'].append(jf.get('fn'))
+    return out_data
+
+
+def cluster_all_mols(folder):
+    mols = mol_file_generator(folder)
+    molecules = [Chem.MolFromMolFile(x) for x in mols]
+    out_data = run_lig_cluster(molecules, mols)
+    return out_data
+
+
+@dataclasses.dataclass()
+class Sites:
+    out_data: dict
+    mol_files: list
+    folder: str
+    missing_mols: list
+
+    @staticmethod
+    def from_folder(folder, recalculate=True):
+        missing_mols = []
+        if recalculate:
+            outdata = cluster_all_mols(folder)
+        else:
+            try:
+                outdata = recluster_from_jsons(folder)
+                missing_mols = [x for x in mol_file_generator(
+                    folder) if not os.path.exists(x.replace('.mol', '_sites.json'))]
+            except RuntimeError:
+                outdata = cluster_all_mols(folder)
+        return Sites(out_data=outdata, missing_mols=missing_mols, mol_files=mol_file_generator(folder), folder=folder)
+
+    def to_json(self):
+        molids = []
+        for key in self.out_data.keys():
+            for site in self.out_data.get(key).keys():
+                molids = molids + \
+                    self.out_data.get(key).get(site).get('mol_ids')
+        for x in set(molids):
+            out_dict = {}
+            for y in self.out_data.keys():
+                site_ids = [z for z in self.out_data.get(
+                    y).keys() if x in self.out_data.get(y).get(z).get('mol_ids')]
+                site_centre_of_mass = [self.out_data.get(y).get(z).get('centre_of_mass') for z in self.out_data.get(
+                    y).keys() if x in self.out_data.get(y).get(z).get('mol_ids')]
+                if(len(site_ids) > 0):
+                    out_dict[y] = {'site_id': site_ids,
+                                   'site_centre_of_mass': site_centre_of_mass}
+            fn = x.replace('.mol', '_sites.json')
+            print(f'Writing Sites to {fn}...')
+            with open(fn, 'w') as f:
+                json.dump(out_dict, f)
+
+    def apply_to_metadata(self):
+        for mf in self.mol_files:
+            mf_csv = mf.replace('.mol', '_meta.csv')
+            mf_json = mf.replace('.mol', '_sites.json')
+            jf = read_jfile(mf_json)
+            site_str = str(jf['c_of_m']['site_id'][0] + 1).zfill(2)
+            with open(mf_csv, 'r') as csvfile:
+                reader = csv.reader(csvfile)
+                row = [x for x in reader][0]
+            row[6] = f'Site {site_str}'
+            with open(mf_csv, 'w', newline='') as csvfile:
+                writer = csv.writer(csvfile, delimiter=',')
+                writer.writerow(row)
+
+    def curate_tags(self):
+        # Currently undefined, placeholder
+        pass
+
+    def cluster_get_or_create(self, key):
+        created = False
+        if self.out_data.get(key) is None:
+            created = True
+            self.out_data[key] = {}
+        return (self.out_data.get(key), created)
+
+    def next_site_number(self, key):
+        vals = list(self.out_data[key].keys())
+        site_num = 0
+        while site_num in vals:
+            site_num += 1
+        return site_num
+
+    def tolerate_sites(self, key, com, tolerance, fn):
+        cluster, created = self.cluster_get_or_create(key=key)
+        if not created:
+            if len(cluster.keys()) > 0:
+                coms = [cluster.get(x).get('centre_of_mass')
+                        for x in cluster.keys()]
+                dists = [get_distance(com, y) for y in coms]
+                trues = [x < tolerance for x in dists]
+                if any(trues):
+                    idx = np.argmin(dists)
+                    key_idx = list(cluster.keys())[idx]
+                    cluster.get(key_idx).get('mol_ids').append(fn)
+                else:
+                    site_num = self.next_site_number(key=key)
+                    cluster[site_num] = {
+                        'centre_of_mass': com, 'mol_ids': [fn]}
+            else:
+                cluster[0] = {'centre_of_mass': com, 'mol_ids': [fn]}
+        else:
+            cluster[0] = {'centre_of_mass': com, 'mol_ids': [fn]}
+
+    def append_data(self, dict, filename, com_tolerance=5, other_tolerance=1):
+        for key in dict.keys():
+            if key == 'c_of_m':
+                tolerance = com_tolerance
+            else:
+                tolerance = other_tolerance
+            for s in dict.get(key).keys():
+                for com in [dict.get(key).get(s).get('centre_of_mass')]:
+                    self.tolerate_sites(
+                        key=key, com=com, tolerance=tolerance, fn=filename)
+
+    def append_mol(self, mol, com_tolerance=5, other_tolerance=1):
+        x_out_data = run_lig_cluster([Chem.MolFromMolFile(mol)], [mol])
+        self.append_data(dict=x_out_data, com_tolerance=com_tolerance,
+                         other_tolerance=other_tolerance, filename=mol)
+
+    def append_json(self, jsonfile, com_tolerance=5, other_tolerance=1):
+        jf = read_jfile(jsonfile)
+        for key in jf.keys():
+            if key == 'fn':
+                continue
+            if key == 'c_of_m':
+                tolerance = com_tolerance
+            else:
+                tolerance = other_tolerance
+                for com in dict.get(key).get('site_centre_of_mass'):
+                    self.tolerate_sites(
+                        key=key, com=com, tolerance=tolerance, fn=jf.get('fn'))
+
+    def cluster_missing_mols(self, com_tolerance=5, other_tolerance=1):
+        for x in self.missing_mols:
+            self.append_mol(mol=x, com_tolerance=com_tolerance,
+                            other_tolerance=other_tolerance)
+
+
+def centre_of_mass(mol):
+    numatoms = mol.GetNumAtoms()
+    conf = mol.GetConformer()
+    if not conf.Is3D():
+        return 0
+    # get coordinate of each atoms
+    pts = np.array([list(conf.GetAtomPosition(atmidx))
+                   for atmidx in range(numatoms)])
+    atoms = [atom for atom in mol.GetAtoms()]
+    mass = Descriptors.MolWt(mol)
+    # get center of mass
+    center_of_mass = np.array(
+        sum(atoms[i].GetMass() * pts[i] for i in range(numatoms))) / mass
+    return center_of_mass
+
+
+def compare_2_mols(molfile, pair):
+    # Imagine being this obtuse.
+    # , isomericSmiles=False, sanitize=True, canonicalize=True)
+    mol1, mol2 = (Chem.MolFromMolFile(x) for x in [molfile, pair])
+    smi1, smi2 = (Chem.MolToSmiles(m) for m in [mol1, mol2])
+    fps = [Chem.RDKFingerprint(Chem.MolFromSmiles(x)) for x in [smi1, smi2]]
+    similarity = DataStructs.FingerprintSimilarity(fps[0], fps[1])
+    distance = get_distance(centre_of_mass(mol1), centre_of_mass(mol2))
+    num1, num2 = (x.GetNumAtoms() for x in [mol1, mol2])
+    conf1, conf2 = (x.GetConformer() for x in [mol1, mol2])
+    pts1 = [list(conf1.GetAtomPosition(atmidx)) for atmidx in range(num1)]
+    pts2 = [list(conf2.GetAtomPosition(atmidx)) for atmidx in range(num2)]
+    dist_list = []
+    for p in pts1:
+        for q in pts2:
+            dist_list.append(get_distance(p, q))
+            min_d = min(dist_list)
+    # Some rules:
+    # If similarity >= .95 + closeatom_Dictance = Really Small) : Alternative Conformation # If same chain Conf, else Pose
+    # If similarity >= .95 + closeatom_Distance = Big : Alternative Position
+    # If similarity < % + closeatomDistance = big, Primary
+    # If similarity < % + closeatomDistance = big, Secondary
+    if similarity >= 0.95:
+        if min_d <= 5:  # 5 Angstroms enough?
+            if molfile.replace('.mol', '')[-1] == pair.replace('.mol', '')[-1]:
+                relationship = 'Alternative Conformation'
+            else:
+                relationship = 'Alternative Pose'
+        else:
+            relationship = 'Alternative Position'
+    else:
+        if min_d <= 5:
+            relationship = 'Primary'
+        else:
+            relationship = 'Secondary'
+    return (similarity, distance, min_d, relationship)
+
+
+def contextualize_a_crystal(crystal_list):
+    for molfile in crystal_list:
+        d = {}
+        for pair in list(set(crystal_list) - set([molfile])):
+            similarity, distance, min_d, relationship = compare_2_mols(
+                molfile, pair)
+            d[os.path.basename(pair).replace('.mol', '')] = {
+                'similarity': similarity,
+                'com_Distance': distance,
+                'closeatom_Distance': min_d,
+                'relationship': relationship
+            }
+        fn = molfile.replace('.mol', '_relationships.json')
+        with open(fn, 'w') as f:
+            json.dump(d, f)
+
+
+def contextualize_crystal_ligands(folder):
+    mols = mol_file_generator(folder)
+    unique_mols = set([os.path.basename(x).rsplit('_', 1)[0] for x in mols])
+    dic = {}
+    for i in unique_mols:
+        dic[i] = [x for x in mols if i in x]
+    for crys in dic.keys():
+        if len(dic.get(crys)) > 1:
+            contextualize_a_crystal(dic.get(crys))
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "-i",
+        "--input",
+        help="Input folder",
+        required=True
+    )
+
+    parser.add_argument(
+        "-m",
+        "--mode",
+        choices=['new', 'missing'],
+        required=True
+    )
+
+    parser.add_argument(
+        "-om",
+        "--overwritemeta",
+        action="store_true",
+        help='Include this flag if you would like to automatically assign center of mass sites as site labels',
+        required=False,
+        default=False
+    )
+
+    parser.add_argument(
+        "-cc",
+        "--cocrystal",
+        action="store_true",
+        help='Include this flag if you want to define relationships between all ligands within each crystal',
+        required=False,
+        default=False
+    )
+
+    parser.add_argument(
+        "-ct",
+        "--comtolerance",
+        help="Tolerance value for creating new clusters for centre of mass sites",
+        type=float,
+        default=5.00
+    )
+
+    parser.add_argument(
+        "-ot",
+        "--othertolerance",
+        help="Tolerance value for creating new clusters for non centre of mass sites",
+        type=float,
+        default=1.00
+    )
+
+    args = vars(parser.parse_args())
+    folder = args["input"]
+    mode = args["mode"]
+    ot = args["othertolerance"]
+    ct = args["comtolerance"]
+    cc = args["cocrystal"]
+    om = args["overwritemeta"]
+
+    # If recalculate is True then stuff happens, if False it reads from jsons.
+    site_obj = Sites.from_folder(folder, recalculate=mode == 'new')
+    if mode == 'missing':
+        site_obj.cluster_missing_mols(com_tolerance=ct, other_tolerance=ot)
+    site_obj.to_json()
+    if cc:
+        contextualize_crystal_ligands(folder=folder)  # Add arguments?
+    if om:
+        site_obj.apply_to_metadata()
diff --git a/fragalysis_api/xcimporter/xcimporter.py b/fragalysis_api/xcimporter/xcimporter.py
index 072d09aa..2d4ae76e 100755
--- a/fragalysis_api/xcimporter/xcimporter.py
+++ b/fragalysis_api/xcimporter/xcimporter.py
@@ -1,17 +1,18 @@
 import argparse
-import shutil
 from sys import exit
 import os
+import glob
 
 from shutil import copyfile, rmtree
 
-from fragalysis_api import Validate, Align, Monomerize
-from fragalysis_api import set_up
+from fragalysis_api import Validate, Align, Sites, contextualize_crystal_ligands
+from fragalysis_api import set_up, convert_small_AA_chains, copy_extra_files
 
 from distutils.dir_util import copy_tree
 
 
-def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomerize=False, biomol=None, covalent=False, pdb_ref=""):
+def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, reduce_reference_frame=False, biomol=None, covalent=False,
+               pdb_ref="", max_lig_len=0):
     """Formats a lists of PDB files into fragalysis friendly format.
     1. Validates the naming of the pdbs.
     2. It aligns the pdbs (_bound.pdb file).
@@ -25,11 +26,12 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
     :param target: Name of the folder to be created inside out_dir
     :param metadata: If set to 1 will create a csv file called metadata.csv in target directory
     :param validate: Validates, and explicitly, warns if input PDB files are not suitable
-    :param monomerize: Bool, if True, will attempt to split pdb files into seperate chains
+    :param reduce_reference_frame: Boolean, if True, will attempt to the number of reference frames by aligning to the first chain in the reference file.
     :param biomol: plain-text file containing header information about the bio-molecular
-        context of the pdb structures. If provided the contents will be appended to the top of the _apo.pdb files
+        context of the pdb structures. If provided the contents will be appended to the top of the _apo.pdb files, Now defunct?
     :param covalent: Bool, if True, will attempt to convert output .mol files to account for potential covalent attachments
     :pdb_ref: String, if provided, all pdb files will be aligned to the name of the file (sans extnesion) that is specified.
+    :max_lig_len: Integer, If >0 will convert all chains with fewer than max_lig_len residues to HETATM with the name LIG. [Currently broken, yikes]
     :return: Hopefully, beautifully aligned files that be used with the fragalysis loader :)
     """
 
@@ -57,13 +59,18 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
 
     in_dir2 = in_dir
 
-    if monomerize:
-        print("Monomerizing input structures")
-        out = os.path.join(out_dir, f'mono{target}/')
+    # Experimental - if option is used then chains are converted. All files need to be moved????
+    if int(max_lig_len) > int(0):
+        print(
+            f'EXPERIMENTAL: Converting all chains with less than {max_lig_len} residues to HETATM LIG')
+        out = os.path.join(out_dir, f'maxliglen{target}/')
         if not os.path.isdir(out):
             os.makedirs(out)
-        mono = Monomerize(directory=in_dir, outdir=out)
-        mono.monomerize_all()
+        infiles = glob.glob(os.path.join(in_dir, '*.pdb'))
+        for i in infiles:
+            convert_small_AA_chains(in_file=i, out_file=os.path.join(
+                out, os.path.basename(i)), max_len=max_lig_len)
+            copy_extra_files(in_file=i, out_dir=out)
         in_dir = out
 
     pdb_smiles_dict = {'pdb': [], 'smiles': []}
@@ -73,33 +80,37 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
             pdb_smiles_dict['pdb'].append(os.path.join(in_dir, f))
             print(os.path.join(in_dir, f).replace('.pdb', '_smiles.txt'))
             if os.path.isfile(os.path.join(in_dir, f).replace('.pdb', '_smiles.txt')):
-                pdb_smiles_dict['smiles'].append(os.path.join(in_dir, f).replace('.pdb', '_smiles.txt'))
+                pdb_smiles_dict['smiles'].append(os.path.join(
+                    in_dir, f).replace('.pdb', '_smiles.txt'))
             else:
                 pdb_smiles_dict['smiles'].append(None)
 
     print(pdb_smiles_dict['smiles'])
     print("Aligning protein structures")
-    print('New Stuff')
-    structure = Align(directory=in_dir, pdb_ref=pdb_ref, mono=monomerize)
+    structure = Align(directory=in_dir, pdb_ref=pdb_ref,
+                      rrf=reduce_reference_frame)
     structure.align(out_dir=os.path.join(out_dir, f"tmp{target}"))
 
     for smiles_file in pdb_smiles_dict['smiles']:
         if smiles_file:
             print(smiles_file)
-            copyfile(smiles_file, os.path.join(os.path.join(out_dir, f"tmp{target}", smiles_file.split('/')[-1])))
-            print(os.path.join(out_dir, f"tmp{target}", smiles_file.split('/')[-1]))
+            copyfile(smiles_file, os.path.join(os.path.join(
+                out_dir, f"tmp{target}", smiles_file.split('/')[-1])))
+            print(os.path.join(
+                out_dir, f"tmp{target}", smiles_file.split('/')[-1]))
 
-    aligned_dict = {'bound_pdb':[], 'smiles':[]}
+    aligned_dict = {'bound_pdb': [], 'smiles': []}
 
     for f in os.listdir(os.path.join(out_dir, f"tmp{target}")):
         if '.pdb' in f:
-            aligned_dict['bound_pdb'].append(os.path.join(out_dir, f"tmp{target}",f))
-            if os.path.isfile(os.path.join(out_dir, f"tmp{target}",f).replace('_bound.pdb', '_smiles.txt')):
-                aligned_dict['smiles'].append(os.path.join(out_dir, f"tmp{target}",f).replace('_bound.pdb', '_smiles.txt'))
+            aligned_dict['bound_pdb'].append(
+                os.path.join(out_dir, f"tmp{target}", f))
+            if os.path.isfile(os.path.join(out_dir, f"tmp{target}", f).replace('_bound.pdb', '_smiles.txt')):
+                aligned_dict['smiles'].append(
+                    os.path.join(out_dir, f"tmp{target}", f).replace('_bound.pdb', '_smiles.txt'))
             else:
                 aligned_dict['smiles'].append(None)
 
-
     print(aligned_dict['smiles'])
     print("Identifying ligands")
 
@@ -109,21 +120,21 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
                 _ = set_up(target_name=target,
                            infile=os.path.abspath(aligned),
                            out_dir=out_dir,
-                           monomerize=monomerize,
+                           rrf=reduce_reference_frame,
                            smiles_file=os.path.abspath(smiles),
                            biomol=biomol,
                            covalent=covalent,
-                           keep_headers=False)
-                
+                           keep_headers=True)
+
             else:
                 _ = set_up(target_name=target,
                            infile=os.path.abspath(aligned),
                            out_dir=out_dir,
-                           monomerize=monomerize,
+                           rrf=reduce_reference_frame,
                            biomol=biomol,
                            covalent=covalent,
-                           keep_headers=False)
-                
+                           keep_headers=True)
+
         except AssertionError:
             print(aligned, "is not suitable, please consider removal or editing")
             for file in os.listdir(os.path.join(out_dir, f"tmp{target}")):
@@ -144,18 +155,18 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
                             f.write(line)
 
     # Copy reference pdb to aligned folder as: reference.pdb, so single_import can file off with ease.
-    structure.write_align_ref(os.path.join(out_dir, target, 'reference.pdb'))
+    structure.write_align_ref(output=os.path.join(out_dir, target))
 
     # Move input files into Target/crystallographic folder
     copy_tree(in_dir2, os.path.join(out_dir, target, 'crystallographic'))
 
-    if os.path.exists(os.path.join(out_dir, f'mono{target}')):
-        shutil.rmtree(os.path.join(out_dir, f'mono{target}'))
-    if os.path.exists(os.path.join(out_dir, f'tmp{target}')):
-        shutil.rmtree(os.path.join(out_dir, f'tmp{target}'))
-
+    # Time to use a for loop?
+    clean_up = [os.path.join(out_dir, f'maxliglen{target}'), os.path.join(
+        out_dir, f'mono{target}'), os.path.join(out_dir, f'tmp{target}')]
+    [rmtree(x) for x in clean_up if os.path.exists(x)]
     print("Files are now in a fragalysis friendly format!")
 
+
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
 
@@ -177,30 +188,66 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
         "-v", "--validate", action="store_true", default=False, help="Validate input"
     )
     parser.add_argument(
-        "-m", "--monomerize", action="store_true", default=False, help="Monomerize input"
+        "-rrf", "--reduce_reference_frame", action="store_true", default=False, help="Reduce Reference frame to single chain?"
     )
     parser.add_argument("-t", "--target", help="Target name", required=True)
-    parser.add_argument("-md", "--metadata",action="store_true", help="Metadata output", default=False)
-    parser.add_argument("-b", "--biomol_txt", help="Biomol Input txt file", required=False, default=None)
-    parser.add_argument('-r', '--reference', help='Reference Structure', required=False, default=None)
+    parser.add_argument("-md", "--metadata", action="store_true",
+                        help="Metadata output", default=False)
+    parser.add_argument(
+        "-b", "--biomol_txt", help="Biomol Input txt file", required=False, default=None)
+    parser.add_argument(
+        '-r', '--reference', help='Reference Structure', required=False, default=None)
     parser.add_argument("-c",
                         "--covalent",
                         action="store_true",
                         help="Handle covalent bonds between ligand and target",
                         required=False,
                         default=False)
+    parser.add_argument("-mll",
+                        "--max_lig_len",
+                        help="Int, Convert all chains shorter than max_lig_len to HETATM LIG",
+                        required=False,
+                        default=0)
 
+    parser.add_argument(
+        "-cs",
+        "--cluster_sites",
+        action="store_true",
+        help='Include this flag if you would like to automatically assign center of mass sites as site labels',
+        required=False,
+        default=False
+    )
+
+    parser.add_argument(
+        "-cs_com",
+        "--cluster_sites_com",
+        help="Tolerance value for creating new clusters for centre of mass sites",
+        type=float,
+        default=5.00
+    )
+
+    parser.add_argument(
+        "-cs_other",
+        "--cluster_sites_other",
+        help="Tolerance value for creating new clusters for non centre of mass sites",
+        type=float,
+        default=1.00
+    )
     args = vars(parser.parse_args())
 
     # user_id = args['user_id']
     in_dir = args["in_dir"]
     out_dir = args["out_dir"]
     validate = args["validate"]
-    monomerize = args["monomerize"]
+    reduce_reference_frame = args["reduce_reference_frame"]
     target = args["target"]
     metadata = args["metadata"]
     biomol = args["biomol_txt"]
     covalent = args["covalent"]
+    mll = args['max_lig_len']
+    cs = args['cluster_sites']
+    cs_com = args['cluster_sites_com']
+    cs_other = args['cluster_sites_other']
 
     if args['reference'] is None:
         print('Reference not set')
@@ -218,20 +265,18 @@ def xcimporter(in_dir, out_dir, target, metadata=False, validate=False, monomeri
                out_dir=out_dir,
                target=target,
                validate=validate,
-               monomerize=monomerize,
+               reduce_reference_frame=reduce_reference_frame,
                metadata=metadata,
                biomol=biomol,
                covalent=covalent,
-               pdb_ref=reference)
-
-    fix_pdb = open(os.path.join(out_dir, target, 'aligned', 'pdb_file_failures.txt'), 'w')
-
-    for target_file in os.listdir(os.path.join(out_dir, target)):
-        if target_file != 'pdb_file_failures.txt' and len(os.listdir(os.path.join(out_dir, target, target_file))) < 2:
-            rmtree(os.path.join(out_dir, target, target_file))
-            fix_pdb.write(target_file.split('-')[1]+'\n')
-
-    fix_pdb.close()
-    print('For files that we were unable to process, look at the pdb_file_failures.txt file in your results directory.'
-          ' These files were unable to produce RDKit molecules, so the error likely lies in the way the ligand atoms or'
-          'the conect files have been written in the pdb file')
+               pdb_ref=reference,
+               max_lig_len=mll
+               )
+    if cs:
+        folder = os.path.join(out_dir, target)
+        site_obj = Sites.from_folder(folder, recalculate=False)
+        site_obj.cluster_missing_mols(
+            com_tolerance=cs_com, other_tolerance=cs_other)
+        site_obj.to_json()
+        contextualize_crystal_ligands(folder=folder)
+        site_obj.apply_to_metadata()
diff --git a/fragalysis_api/xcimporter_utils.py b/fragalysis_api/xcimporter_utils.py
index 8fa6b3ff..1b5d249d 100755
--- a/fragalysis_api/xcimporter_utils.py
+++ b/fragalysis_api/xcimporter_utils.py
@@ -1,4 +1,5 @@
-from fragalysis_api.xcimporter.align import Align, Monomerize
-from fragalysis_api.xcimporter.conversion_pdb_mol import Ligand, pdb_apo, set_up
+from fragalysis_api.xcimporter.align import Align
+from fragalysis_api.xcimporter.sites import Sites, contextualize_crystal_ligands
+from fragalysis_api.xcimporter.conversion_pdb_mol import Ligand, pdb_apo, set_up, convert_small_AA_chains, copy_extra_files
 from fragalysis_api.xcimporter.validate import Validate, ValidatePDB
-from fragalysis_api.xcimporter.xc_utils import to_fragalysis_dir
\ No newline at end of file
+from fragalysis_api.xcimporter.xc_utils import to_fragalysis_dir
diff --git a/setup.py b/setup.py
index ceb4ec3d..932a787c 100755
--- a/setup.py
+++ b/setup.py
@@ -6,7 +6,7 @@
 
 setuptools.setup(
     name='fragalysis_api',
-    version='0.0.2.2',
+    version='0.0.5.1',
     author='Fragment 5 + Tyler',
     author_email="",
     description="A package to load PDBs into fragalysis format.",
@@ -19,7 +19,8 @@
                       'pandas',
                       'pypdb',
                       'matplotlib',
-                      'scipy'],  # Install requirements extracted from requirements.txt
+                      'scipy',
+                      'fragalysis'],  # Install requirements extracted from requirements.txt
     include_package_data=True,  # Allow to include other files than .py in package
     package_data={
         '': ['fragalysis_api/xcimporter/non_ligs.json',
@@ -32,4 +33,4 @@
         "Operating System :: OS Independent",
     ],
     python_requires='>=3.6',
-    )
+)


Test patch:
diff --git a/tests/xcimporter/test_align.py b/tests/xcimporter/test_align.py
index f53ad471..b505adf4 100755
--- a/tests/xcimporter/test_align.py
+++ b/tests/xcimporter/test_align.py
@@ -3,7 +3,7 @@
 import os
 from pathlib import Path
 
-from fragalysis_api import Align, Monomerize, set_up
+from fragalysis_api import Align, set_up
 from glob import glob
 from shutil import rmtree
 
@@ -23,16 +23,16 @@ class EasyAlign(AlignTest):
     @classmethod
     def setUpClass(cls):
         super(EasyAlign, cls).setUpClass()
-        cls.align_obj = Align(os.path.join(cls.dir_input, 'examples_to_test0'), pdb_ref='')
-        cls.align_obj_w_ref = Align(os.path.join(cls.dir_input, 'examples_to_test0'), pdb_ref='6hi3')
-        cls.align_obj_w_wrong_ref = Align(os.path.join(cls.dir_input, 'examples_to_test0'), pdb_ref='wrong_pdb')
-        cls.align_obj_w_maps = Align(os.path.join(cls.dir_input, 'examples_to_test5'), pdb_ref='')
-
-        # Do a monomerize
-        out = os.path.join(cls.dir_output, f'mono')
-        if not os.path.isdir(out):
-            os.makedirs(out)
-        cls.monomerize_obj = Monomerize(directory=os.path.join(cls.dir_input, 'examples_to_test5'), outdir=out)
+        cls.align_obj = Align(os.path.join(
+            cls.dir_input, 'examples_to_test0'), pdb_ref='')
+        cls.align_obj_w_ref = Align(os.path.join(
+            cls.dir_input, 'examples_to_test0'), pdb_ref='6hi3')
+        cls.align_obj_w_wrong_ref = Align(os.path.join(
+            cls.dir_input, 'examples_to_test0'), pdb_ref='wrong_pdb')
+        cls.align_obj_w_maps = Align(os.path.join(
+            cls.dir_input, 'examples_to_test5'), pdb_ref='', rrf=False)
+        cls.align_obj_w_maps_rrf = Align(os.path.join(
+            cls.dir_input, 'examples_to_test5'), pdb_ref='', rrf=True)
 
     @classmethod
     def tearDownClass(cls):
@@ -41,125 +41,65 @@ def tearDownClass(cls):
 
     def test_get_files(self):
         files = ['6epu.pdb', '6epv.pdb', '6epx.pdb', '6hi3.pdb']
-        whole_files = sorted([os.path.join(self.align_obj.directory, x) for x in files])
+        whole_files = sorted(
+            [os.path.join(self.align_obj.directory, x) for x in files])
         self.assertCountEqual(self.align_obj._get_files, whole_files)
 
     def test_get_mapfiles(self):
         self.assertCountEqual(self.align_obj._get_maplist, [])
 
-    def test_a_monomerize_single(self):
-        dir = os.path.join('tests', 'data_for_tests', 'mono')
-        if not os.path.exists(dir):
-            os.makedirs(dir)
-        self.monomerize_obj.monomerize_single(
-            file=os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x2097.pdb'))
-        single_test_cases = ['Mpro-x2097_A.pdb', 'Mpro-x2097_A_event.ccp4']
-        single_test_exists = [os.path.exists(os.path.join('tests', 'data_for_tests', 'mono', x)) for x in
-                              single_test_cases]
-        print(single_test_exists)
-        [self.assertTrue(x) for x in single_test_exists]
-
-    def test_b_monomerize_all(self):
-        dir = os.path.join('tests', 'data_for_tests', 'mono')
-        if not os.path.exists(dir):
-            os.makedirs(dir)
-        self.monomerize_obj.monomerize_all()
-        all_test_cases = ['Mpro-x0978_A.pdb', 'Mpro-x0981_A.pdb',
-                          'Mpro-x2097_A.pdb', 'Mpro-x2097_A_event.ccp4',
-                          'Mpro-x2119_A.pdb', 'Mpro-x2119_A_event.ccp4']
-        all_test_exists = [os.path.exists(os.path.join('tests', 'data_for_tests', 'mono', x)) for x in all_test_cases]
-        print(all_test_exists)
-        [self.assertTrue(x) for x in all_test_exists]
-
-    def test_c_align_from_mono(self):
-        print(os.path.exists(os.path.join('tests', 'data_for_tests', 'mono')))
-        a = Align(os.path.join('tests', 'data_for_tests', 'mono'), pdb_ref='', mono=True)
-        dir = os.path.join('tests', 'data_for_tests', 'tmp')
-        if not os.path.exists(dir):
-            os.makedirs(dir)
-        a.align(out_dir=dir)
-
-        a_test_cases = ['Mpro-x0978_A_bound.pdb', 'Mpro-x0981_A_bound.pdb', 'Mpro-x2097_A_bound.pdb',
-                        'Mpro-x2119_A_bound.pdb', 'Mpro-x2097_A_event.ccp4', 'Mpro-x2119_A_event.ccp4']
-        a_test_exists = [os.path.exists(os.path.join(dir, x)) for x in a_test_cases]
-        print(a_test_exists)
-        [self.assertTrue(x) for x in a_test_exists]
-
-        dir = os.path.join('tests', 'data_for_tests', 'tmpa2')
-        if not os.path.exists(dir):
-            os.makedirs(dir)
-
-        a2 = Align(os.path.join('tests', 'data_for_tests', 'mono'), pdb_ref='Mpro-x2119', mono=True)
-        a2.align(out_dir=dir)
-
-        a2_test_cases = ['Mpro-x0978_A_bound.pdb', 'Mpro-x0981_A_bound.pdb', 'Mpro-x2097_A_bound.pdb',
-                         'Mpro-x2119_A_bound.pdb', 'Mpro-x2097_A_event.ccp4', 'Mpro-x2119_A_event.ccp4']
-        a2_test_exists = [os.path.exists(os.path.join(dir, x)) for x in a2_test_cases]
-        print(a2_test_exists)
-        [self.assertTrue(x) for x in a2_test_exists]
-
-        # Also test if error state ruins things
-        # error = Align(os.path.join('tests', 'data_for_tests', 'mono'), pdb_ref='Mpro-x2119', mono=False)
-
-    def test_d_align_with_maps(self):
+    def test_a_align_with_maps(self):
         dir = os.path.join('tests', 'data_for_tests', 'tmp_map')
         if not os.path.exists(dir):
             os.makedirs(dir)
         self.align_obj_w_maps.align(out_dir=dir)
         map_test_cases = ['Mpro-x0978_bound.pdb', 'Mpro-x0981_bound.pdb', 'Mpro-x2097_bound.pdb',
                           'Mpro-x2119_bound.pdb', 'Mpro-x2097_event.ccp4', 'Mpro-x2119_event.ccp4']
-        map_test_exists = [os.path.exists(os.path.join(dir, x)) for x in map_test_cases]
-        print(map_test_exists)
+        map_test_exists = [os.path.exists(
+            os.path.join(dir, x)) for x in map_test_cases]
         [self.assertTrue(x) for x in map_test_exists]
 
-    def test_e_align_to_reference_with_map(self):
+    def test_b_align_to_reference_with_map(self):
         dir = os.path.join('tests', 'data_for_tests', 'tmp_map_single')
-        file = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x2119.pdb')
-        ref = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x0978.pdb')
+        file = os.path.join('tests', 'data_for_tests',
+                            'examples_to_test5', 'Mpro-x2119.pdb')
+        ref = os.path.join('tests', 'data_for_tests',
+                           'examples_to_test5', 'Mpro-x0978.pdb')
         if not os.path.exists(dir):
             os.makedirs(dir)
         self.align_obj_w_maps.align_to_reference(in_file=file,
-                                                 reference=ref,
+                                                 reference_pdb=ref,
                                                  out_dir=dir)
         test_cases = ['Mpro-x2119_bound.pdb', 'Mpro-x2119_event.ccp4']
-        test_exists = [os.path.exists(os.path.join(dir, x)) for x in test_cases]
+        test_exists = [os.path.exists(os.path.join(dir, x))
+                       for x in test_cases]
         [self.assertTrue(x) for x in test_exists]
 
+    def test_align_rrf(self):
+        dir = os.path.join('tests', 'data_for_tests', 'tmp_rrf')
+        if not os.path.exists(dir):
+            os.makedirs(dir)
+        self.align_obj_w_maps_rrf.align(out_dir=dir)
+        rrf_test_cases = ['Mpro-x0978_A_bound.pdb', 'Mpro-x0981_A_bound.pdb', 'Mpro-x2097_A_bound.pdb',
+                          'Mpro-x2119_A_bound.pdb', 'Mpro-x2097_A_event.ccp4', 'Mpro-x2119_A_event.ccp4']
+        rrf_test_exists = [os.path.exists(
+            os.path.join(dir, x)) for x in rrf_test_cases]
+        [self.assertTrue(x) for x in rrf_test_exists]
+
     def test_z_structure_class_tests(self):
-        path = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x2119.pdb')
+        path = os.path.join('tests', 'data_for_tests',
+                            'examples_to_test5', 'Mpro-x2119.pdb')
         struc = Structure.from_file(Path(path))
         self.assertEqual(struc.rfree().to_float(), 0.2173)
         self.assertEqual(len(struc.residue_ids()), 304)
         self.assertEqual(len([i for i in struc.protein_atoms()]), 2374)
         self.assertEqual(len([i for i in struc.all_atoms()]), 2732)
 
-    def test_f_monomerize_multichain(self):
-        dir = os.path.join('tests', 'data_for_tests', 'monomc')
-        if not os.path.exists(dir):
-            os.makedirs(dir)
-
-        in_dir = os.path.join('tests', 'data_for_tests', 'examples_to_test6')
-        mono = Monomerize(directory=in_dir, outdir=dir)
-        mono.monomerize_all()
-        test_cases = [
-            'mArh-x0091_A.pdb',
-            'mArh-x0104_A.pdb',
-            'mArh-x0128_A.pdb',
-            'mArh-x0142_A.pdb',
-            'mArh-x0158_A.pdb',
-            'mArh-x0091_B.pdb',
-            'mArh-x0104_B.pdb',
-            'mArh-x0128_B.pdb',
-            'mArh-x0142_B.pdb',
-            'mArh-x0158_B.pdb'
-        ]
-        test_exists = [os.path.exists(os.path.join(dir, x)) for x in test_cases]
-        [self.assertTrue(x) for x in test_exists]
-
     def test_g_conversion_pdb_mol(self):
         dir = os.path.join('tests', 'data_for_tests', 'conv')
         dir2 = os.path.join(dir, 'target')
-        insmiles = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x2119_smiles.txt')
+        insmiles = os.path.join('tests', 'data_for_tests',
+                                'examples_to_test5', 'Mpro-x2119_smiles.txt')
         # Requires aligned pdb file and text file...
         for mono in [False, True]:
             for covalent in [False, True]:
@@ -171,7 +111,8 @@ def test_g_conversion_pdb_mol(self):
                             os.makedirs(dir2)
 
                         if mono:
-                            inpdb = os.path.join('tests', 'data_for_tests', 'tmp', 'Mpro-x2119_A_bound.pdb')
+                            inpdb = os.path.join(
+                                'tests', 'data_for_tests', 'tmp_rrf', 'Mpro-x2119_A_bound.pdb')
                             final_path_dir = os.path.join('tests', 'data_for_tests', 'conv', 'target', 'aligned',
                                                           'target-Mpro-x2119_0A')
                             test_cases = ['target-Mpro-x2119_0A.mol',
@@ -184,7 +125,8 @@ def test_g_conversion_pdb_mol(self):
                                           'target-Mpro-x2119_0A_apo-desolv.pdb',
                                           'target-Mpro-x2119_0A_event.ccp4']
                         else:
-                            inpdb = os.path.join('tests', 'data_for_tests', 'tmp_map', 'Mpro-x2119_bound.pdb')
+                            inpdb = os.path.join(
+                                'tests', 'data_for_tests', 'tmp_map', 'Mpro-x2119_bound.pdb')
                             final_path_dir = os.path.join('tests', 'data_for_tests', 'conv', 'target', 'aligned',
                                                           'target-Mpro-x2119_0')
                             test_cases = ['target-Mpro-x2119_0.mol',
@@ -201,27 +143,20 @@ def test_g_conversion_pdb_mol(self):
                             target_name='target',
                             infile=os.path.abspath(inpdb),
                             out_dir=dir,
-                            monomerize=mono,
+                            rrf=mono,
                             smiles_file=smiles_file,
                             biomol=biomol,
                             covalent=covalent
                         )
                         # Assert the Files were made:
-                        test_exists = [os.path.exists(os.path.join(final_path_dir, x)) for x in test_cases]
-                        print(mono)
-                        print(smiles_file)
-                        print(biomol)
-                        print(covalent)
-                        print(test_cases)
-                        print(test_exists)
-                        print([i for i in os.walk(final_path_dir)])
+                        test_exists = [os.path.exists(os.path.join(
+                            final_path_dir, x)) for x in test_cases]
                         [self.assertTrue(x) for x in test_exists]
                         # Clean-up the files.
                         shutil.rmtree(dir2)
 
-
     def test_write_ref(self):
-        fp = os.path.join('tests', 'data_for_tests' 'reference.pdb')
+        fp = os.path.join('tests', 'data_for_tests')
         self.align_obj_w_ref.write_align_ref(fp)
         self.assertTrue(os.path.exists(fp))
 
diff --git a/tests/xcimporter/test_conversion.py b/tests/xcimporter/test_conversion.py
index c548a0bc..96fed27a 100755
--- a/tests/xcimporter/test_conversion.py
+++ b/tests/xcimporter/test_conversion.py
@@ -8,7 +8,8 @@ class ConversionTest(unittest.TestCase):
 
     @classmethod
     def setUpClass(cls):
-        cls.dir_input = os.path.join('tests', 'data_for_tests', 'examples_to_test3')
+        cls.dir_input = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test3')
         cls.dir_output = os.path.join('tests', 'data_for_tests', 'tmp')
 
 
@@ -17,9 +18,9 @@ class PDBexample1(ConversionTest):
     @classmethod
     def setUpClass(cls):
         super(PDBexample1, cls).setUpClass()
-        cls.obj_5q1j = set_up(target_name="5q1j", infile=os.path.join(cls.dir_input,'5q1j.pdb'), monomerize = False,
+        cls.obj_5q1j = set_up(target_name="5q1j", infile=os.path.join(cls.dir_input, '5q1j.pdb'), rrf=False,
                               out_dir=cls.dir_output)
-        cls.obj_5qj7 = set_up(target_name="5qj7", infile=os.path.join(cls.dir_input,'5qj7.pdb'), monomerize = False,
+        cls.obj_5qj7 = set_up(target_name="5qj7", infile=os.path.join(cls.dir_input, '5qj7.pdb'), rrf=False,
                               out_dir=cls.dir_output)
 
     @classmethod
@@ -70,22 +71,26 @@ def test_make_pdb_file(self):
         tests that a pdb file has been made for a particular ligand with the correct
         number of lines (HETATM & CONECT)
         """
-        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY, "5qj7_0/5qj7_0.pdb")).readlines()
+        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY,
+                    "5qj7_0/5qj7_0.pdb")).readlines()
         self.assertEqual(len(file), 36)
 
     def test_make_mol_file(self):
         """
         tests that a mol file has been made with the correct number of lines
         """
-        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY, "5qj7_0/5qj7_0.mol")).readlines()
+        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY,
+                    "5qj7_0/5qj7_0.mol")).readlines()
         self.assertEqual(len(file), 43)
 
     def test_make_sdf_file(self):
         """
         test that a sdf file has been made that incorporates the different mol objects
         """
-        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY, "5qj7_0/5qj7_0.sdf")).readlines()
+        file = open(os.path.join(self.obj_5qj7.RESULTS_DIRECTORY,
+                    "5qj7_0/5qj7_0.sdf")).readlines()
         self.assertEqual(len(file), 44)
 
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/xcimporter/test_xcimporter.py b/tests/xcimporter/test_xcimporter.py
index cfa84b6e..f9d629fe 100755
--- a/tests/xcimporter/test_xcimporter.py
+++ b/tests/xcimporter/test_xcimporter.py
@@ -1,7 +1,6 @@
 import unittest
 import os
 from fragalysis_api import xcimporter
-from shutil import rmtree
 
 from fragalysis_api.xcimporter.single_import import import_single_file
 
@@ -10,25 +9,27 @@ class XcImporterTest(unittest.TestCase):
 
     @classmethod
     def setUpClass(cls):
-        cls.in_dir = os.path.join('tests', 'data_for_tests', 'examples_to_test5')
-        cls.in_file = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x0978.pdb')
+        cls.in_dir = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5')
+        cls.in_file = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x0978.pdb')
         cls.in_file_no_ext = 'Mpro-x0978'
         cls.out_dir = os.path.join('tests', 'data_for_tests')
         cls.target = 'TestTarget'
         cls.validate = False
-        cls.monomerize = True
-        cls.biomol = os.path.join('tests', 'data_for_tests', 'examples_to_test5', 'biomol.txt')
+        cls.rrf = True
+        cls.biomol = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5', 'biomol.txt')
         cls.metadata = True
         cls.covalent = True
         cls.reference = 'Mpro-x0981_A'
-
+        cls.mll = -1
 
     @classmethod
     def tearDownClass(cls):
         #rmtree(os.path.join(cls.out_dir, cls.target))
         pass
 
-
     def test_xcimporter(self):
 
         xcimporter(
@@ -36,28 +37,33 @@ def test_xcimporter(self):
             out_dir=self.out_dir,
             target=self.target,
             validate=self.validate,
-            monomerize=self.monomerize,
+            reduce_reference_frame=self.rrf,
             metadata=self.metadata,
             biomol=self.biomol,
             covalent=self.covalent,
-            pdb_ref=self.reference
+            pdb_ref=self.reference,
+            max_lig_len=self.mll
         )
-
         # Write more Tests
-        self.assertTrue(expr=os.path.exists(os.path.join(self.out_dir, self.target)))
-        self.assertTrue(expr=os.path.exists(os.path.join(self.out_dir, self.target, 'reference.pdb')))
+        self.assertTrue(expr=os.path.exists(
+            os.path.join(self.out_dir, self.target)))
+        self.assertTrue(expr=os.path.exists(os.path.join(
+            self.out_dir, self.target, 'reference.pdb')))
 
         import_single_file(in_file=self.in_file,
                            out_dir=self.out_dir,
                            target=self.target,
-                           monomerize=self.monomerize,
-                           reference=os.path.join(self.out_dir, self.target, 'reference.pdb'),
+                           reduce_reference_frame=self.rrf,
+                           reference_pdb=os.path.join(
+                               self.out_dir, self.target, 'reference.pdb'),
                            biomol=self.biomol,
-                           covalent=self.covalent)
+                           covalent=self.covalent,
+                           max_lig_len=self.mll)
 
         # Tests should check if a thing is correctly removed and readded etc...
         # Write Many More... Single import should use a pdb in a seperate test!
-        self.assertTrue(expr=os.path.exists(os.path.join(self.out_dir, self.target, 'reference.pdb')))
+        self.assertTrue(expr=os.path.exists(os.path.join(
+            self.out_dir, self.target, 'reference.pdb')))
 
 
 if __name__ == '__main__':
diff --git a/tests/xcimporter/test_xcsites.py b/tests/xcimporter/test_xcsites.py
new file mode 100644
index 00000000..69b85915
--- /dev/null
+++ b/tests/xcimporter/test_xcsites.py
@@ -0,0 +1,59 @@
+import unittest
+import os
+from fragalysis_api import xcimporter, Sites, contextualize_crystal_ligands
+
+
+class XcSitesTest(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.in_dir = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5')
+        cls.in_file = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5', 'Mpro-x0978.pdb')
+        cls.in_file_no_ext = 'Mpro-x0978'
+        cls.out_dir = os.path.join('tests', 'data_for_tests')
+        cls.target = 'TestTarget'
+        cls.validate = False
+        cls.rrf = True
+        cls.biomol = os.path.join(
+            'tests', 'data_for_tests', 'examples_to_test5', 'biomol.txt')
+        cls.metadata = True
+        cls.covalent = True
+        cls.reference = 'Mpro-x0981_A'
+        cls.mll = -1
+
+    @classmethod
+    def tearDownClass(cls):
+        #rmtree(os.path.join(cls.out_dir, cls.target))
+        pass
+
+    def test_sites(self):
+
+        xcimporter(
+            in_dir=self.in_dir,
+            out_dir=self.out_dir,
+            target=self.target,
+            validate=self.validate,
+            reduce_reference_frame=self.rrf,
+            metadata=self.metadata,
+            biomol=self.biomol,
+            covalent=self.covalent,
+            pdb_ref=self.reference,
+            max_lig_len=self.mll
+        )
+
+        folder = os.path.join(self.out_dir, self.target)
+        site_obj = Sites.from_folder(folder, recalculate=False)
+        # site_obj.cluster_missing_mols(
+        #        folder=folder, com_tolerance=5.00, other_tolerance=1.00)
+        site_obj.to_json()
+        contextualize_crystal_ligands(folder=folder)
+        site_obj.apply_to_metadata()
+
+        self.assertTrue(expr=os.path.exists(os.path.join(
+            self.out_dir, self.target, 'aligned', 'Mpro-x0978_0A', 'Mpro-x0978_0A_sites.json')))
+
+
+if __name__ == '__main__':
+    unittest.main()