Repository: JustinYKC/FEPOPS (language: python)

Issue:
[Software paper] Lacking mention of how features are computed
Refers to the JOSS submission review: https://github.com/openjournals/joss-reviews/issues/5763

The paper is missing one or two sentences on how the features are computed:

- partial charges: Gasteiger charges calculated with RDKit
- Atomic logP contribution: Crippen method as implemented in RDKit
- Hbond donors and acceptors: SMARTS patterns taken from DAYLIGHT

Since the SMARTS pattern are taken from the [DAYLIGHT website](https://www.daylight.com/dayhtml_tutorials/languages/smarts/smarts_examples.html#H_BOND), they must be added to the references. (and ideally the [RDKit](https://www.rdkit.org/docs/Overview.html#citing-the-rdkit) as well)

[Optional] Small fixes and enhancements
These are just small enhancements/fixes that I noticed while going through the code or documentation, feel free to ignore them as I don't consider them blocking for the acceptance of the paper in JOSS:

- [x] Reorganize the benchmark files in a separate folder. Currently the notebooks and the `mol_filter.py` file are at the root of the repository, but to make it clear what they are for (especially the filtering file) they should be in a `benchmark` directory. Also the type hints for return types in `mol_file.py` should be `Optional[X]` instead of `X or None` in some places.
- [x] The `DudePreprocessor` class in the notebook is already defined in package, no need to redefine it in the notebook.
- [x] typo: `descriptor_sds` should be `descriptor_stds` [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L77)
- [x] Use the standard `warnings` or `logging` modules instead of `print` for warnings/errors/info. Also [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L555) the message will not be displayed, and it is a duplicate of [this print call](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L842-L844).
- [ ] The custom K-medoid implementation could be replaced by the one in the `sklearn-extra` package, which includes both PAM and CLARA algorithms described in the original FEPOPS implementation.
- [x] The output of CLI commands is not in a standard format (python list, numpy array…etc) which makes the `calc_sim` and `get_fepops` CLI commands hardly usable since they require manual editing to be parsed. They should be printed in a tab-separated format or JSON instead, or written to a file.
- [x] On the [welcome page of the docs](https://justinykc.github.io/FEPOPS/#:~:text=and%20is%20available%20here%3A%20readme), the link to the readme is broken (it works in the sidebar though)
[Optional] Small fixes and enhancements
These are just small enhancements/fixes that I noticed while going through the code or documentation, feel free to ignore them as I don't consider them blocking for the acceptance of the paper in JOSS:

- [x] Reorganize the benchmark files in a separate folder. Currently the notebooks and the `mol_filter.py` file are at the root of the repository, but to make it clear what they are for (especially the filtering file) they should be in a `benchmark` directory. Also the type hints for return types in `mol_file.py` should be `Optional[X]` instead of `X or None` in some places.
- [x] The `DudePreprocessor` class in the notebook is already defined in package, no need to redefine it in the notebook.
- [x] typo: `descriptor_sds` should be `descriptor_stds` [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L77)
- [x] Use the standard `warnings` or `logging` modules instead of `print` for warnings/errors/info. Also [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L555) the message will not be displayed, and it is a duplicate of [this print call](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L842-L844).
- [ ] The custom K-medoid implementation could be replaced by the one in the `sklearn-extra` package, which includes both PAM and CLARA algorithms described in the original FEPOPS implementation.
- [x] The output of CLI commands is not in a standard format (python list, numpy array…etc) which makes the `calc_sim` and `get_fepops` CLI commands hardly usable since they require manual editing to be parsed. They should be printed in a tab-separated format or JSON instead, or written to a file.
- [x] On the [welcome page of the docs](https://justinykc.github.io/FEPOPS/#:~:text=and%20is%20available%20here%3A%20readme), the link to the readme is broken (it works in the sidebar though)
[Documentation] API docs section not compiled correctly
Refers to the JOSS submission review: https://github.com/openjournals/joss-reviews/issues/5763

The [API section](https://justinykc.github.io/FEPOPS/fepops.html) seem to not be compiled correctly as it is not showing the content of the docstrings for the different classes/functions which is a shame since these are very well documented.

Without this it's hard for me to justify that the core functionality of the software is documented to a satisfactory level as users would have trouble accessing it unless they read the source code.
[Optional] Small fixes and enhancements
These are just small enhancements/fixes that I noticed while going through the code or documentation, feel free to ignore them as I don't consider them blocking for the acceptance of the paper in JOSS:

- [x] Reorganize the benchmark files in a separate folder. Currently the notebooks and the `mol_filter.py` file are at the root of the repository, but to make it clear what they are for (especially the filtering file) they should be in a `benchmark` directory. Also the type hints for return types in `mol_file.py` should be `Optional[X]` instead of `X or None` in some places.
- [x] The `DudePreprocessor` class in the notebook is already defined in package, no need to redefine it in the notebook.
- [x] typo: `descriptor_sds` should be `descriptor_stds` [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L77)
- [x] Use the standard `warnings` or `logging` modules instead of `print` for warnings/errors/info. Also [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L555) the message will not be displayed, and it is a duplicate of [this print call](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L842-L844).
- [ ] The custom K-medoid implementation could be replaced by the one in the `sklearn-extra` package, which includes both PAM and CLARA algorithms described in the original FEPOPS implementation.
- [x] The output of CLI commands is not in a standard format (python list, numpy array…etc) which makes the `calc_sim` and `get_fepops` CLI commands hardly usable since they require manual editing to be parsed. They should be printed in a tab-separated format or JSON instead, or written to a file.
- [x] On the [welcome page of the docs](https://justinykc.github.io/FEPOPS/#:~:text=and%20is%20available%20here%3A%20readme), the link to the readme is broken (it works in the sidebar though)
[Optional] Small fixes and enhancements
These are just small enhancements/fixes that I noticed while going through the code or documentation, feel free to ignore them as I don't consider them blocking for the acceptance of the paper in JOSS:

- [x] Reorganize the benchmark files in a separate folder. Currently the notebooks and the `mol_filter.py` file are at the root of the repository, but to make it clear what they are for (especially the filtering file) they should be in a `benchmark` directory. Also the type hints for return types in `mol_file.py` should be `Optional[X]` instead of `X or None` in some places.
- [x] The `DudePreprocessor` class in the notebook is already defined in package, no need to redefine it in the notebook.
- [x] typo: `descriptor_sds` should be `descriptor_stds` [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L77)
- [x] Use the standard `warnings` or `logging` modules instead of `print` for warnings/errors/info. Also [here](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L555) the message will not be displayed, and it is a duplicate of [this print call](https://github.com/JustinYKC/FEPOPS/blob/master/src/fepops/fepops.py#L842-L844).
- [ ] The custom K-medoid implementation could be replaced by the one in the `sklearn-extra` package, which includes both PAM and CLARA algorithms described in the original FEPOPS implementation.
- [x] The output of CLI commands is not in a standard format (python list, numpy array…etc) which makes the `calc_sim` and `get_fepops` CLI commands hardly usable since they require manual editing to be parsed. They should be printed in a tab-separated format or JSON instead, or written to a file.
- [x] On the [welcome page of the docs](https://justinykc.github.io/FEPOPS/#:~:text=and%20is%20available%20here%3A%20readme), the link to the readme is broken (it works in the sidebar though)


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
new file mode 100644
index 0000000..6889ce8
--- /dev/null
+++ b/.github/workflows/ci.yaml
@@ -0,0 +1,40 @@
+# Runs the black formatter and pytest
+name: OpenFEPOPS CI
+
+on:
+  push:
+    branches: [main, development]
+  pull_request:
+    branches: [main, development]
+
+jobs:
+  build:
+    strategy:
+      matrix:
+        python-version: [3.11]
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v2
+        with:
+          fetch-depth: 0
+          
+      - name: Checkout the appropriate branch
+        run: git checkout ${{ env.BRANCH }}
+              
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v1
+        with:
+          python-version: ${{ matrix.python-version }}
+          
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e .[dev]
+          
+      - name: Run black formatter
+        run: |
+          black .
+      - name: Run tests
+        run: |
+          pytest
diff --git a/Explore_DUDE_diversity_set.ipynb b/Explore_DUDE_diversity_set.ipynb
index 5579775..74830ed 100644
--- a/Explore_DUDE_diversity_set.ipynb
+++ b/Explore_DUDE_diversity_set.ipynb
@@ -14,149 +14,32 @@
     "Perform preprocessing on the DUDE database.  Assumes that the DUDE database has been downloaded and placed in the path \"data/dude/unprocessed\". Here, we expect directories for each target containing 'actives_final.ism' and 'decoys_final.ism'.\n",
     "It may be beneficial to split the below out of this Jupyter notebook and deploy to HPC or a non interactive environment. To pre-cache descriptors, run cache_mols_from_csv() on the DudePreprocessor object below.\n",
     "\n",
-    "The step may also be skipped by downloading the OpenFEPOPS-DUDE data directory created for analysis of the DUDE diversity set from FigShare at https://doi.org/10.6084/m9.figshare.23951445.v1 .  See the OpenFEPOPS README.md for further information."
+    "The step may also be skipped by downloading the OpenFEPOPS-DUDE data directory created for analysis of the DUDE diversity set from FigShare at https://doi.org/10.6084/m9.figshare.23951445.v2."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "data/dude/\n",
+      "Processing the following DUDE targets: []\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Preparing targets: 0it [00:00, ?it/s]\n"
+     ]
+    }
+   ],
    "source": [
-    "from typing import Union, Optional\n",
-    "from pathlib import Path\n",
-    "from tqdm import tqdm\n",
-    "from fepops import OpenFEPOPS\n",
-    "import pandas as pd\n",
-    "import multiprocessing as mp\n",
-    "from rdkit import Chem\n",
-    "\n",
-    "\n",
-    "class DudePreprocessor:\n",
-    "    def __init__(\n",
-    "        self,\n",
-    "        *,\n",
-    "        dude_directory: Union[Path, str] = \"data/dude/\",\n",
-    "    ) -> None:\n",
-    "        print(dude_directory)\n",
-    "        print(type(dude_directory))\n",
-    "        self.dude_path = Path(dude_directory)\n",
-    "        self.dude_unprocessed_path = self.dude_path / Path(\"unprocessed\")\n",
-    "        if not self.dude_path.exists():\n",
-    "            raise FileNotFoundError(f\"Dude dataset not found in path: {self.dude_path}\")\n",
-    "        self.dude_processed_path = self.dude_path / Path(\"processed\")\n",
-    "        self.dude_processed_path.parent.mkdir(parents=True, exist_ok=True)\n",
-    "        self.fepops_ob = OpenFEPOPS()\n",
-    "\n",
-    "    def __call__(\n",
-    "        self,\n",
-    "    ):\n",
-    "        self.process()\n",
-    "\n",
-    "    def process(\n",
-    "        self,\n",
-    "        targets: Optional[Union[str, list[str]]] = None,\n",
-    "        skip_existing: bool = True,\n",
-    "    ):\n",
-    "        if targets is None:\n",
-    "            dude_targets = [\n",
-    "                t.parent.name\n",
-    "                for t in self.dude_path.glob(\"unprocessed/*/actives_final.ism\")\n",
-    "            ]\n",
-    "        else:\n",
-    "            if isinstance(targets, str):\n",
-    "                targets = [targets]\n",
-    "        print(f\"Processing the following DUDE targets: {dude_targets}\")\n",
-    "        for target in tqdm(dude_targets, desc=f\"Preparing targets\"):\n",
-    "            self.create_dude_target_csv_data(target, skip_existing=skip_existing)\n",
-    "\n",
-    "    @staticmethod\n",
-    "    def _parallel_init_worker_desc_gen_shared_fepops_ob():\n",
-    "        global shared_fepops_ob\n",
-    "        shared_fepops_ob = OpenFEPOPS()\n",
-    "\n",
-    "    @staticmethod\n",
-    "    def _parallel_get_rdkit_cansmi(s):\n",
-    "        global shared_fepops_ob\n",
-    "        mol = shared_fepops_ob._mol_from_smiles(s)\n",
-    "        if mol is None:\n",
-    "            return \"\"\n",
-    "        return Chem.MolToSmiles(mol)\n",
-    "\n",
-    "    def create_dude_target_csv_data(\n",
-    "        self,\n",
-    "        dude_target: Path,\n",
-    "        actives_file: Path = Path(\"actives_final.ism\"),\n",
-    "        decoys_file: Path = Path(\"decoys_final.ism\"),\n",
-    "        seperator: str = \" \",\n",
-    "        skip_existing: bool = True,\n",
-    "    ):\n",
-    "        target_output_file = self.dude_processed_path / f\"dude_target_{dude_target}.csv\"\n",
-    "        if skip_existing and target_output_file.exists():\n",
-    "            print(\n",
-    "                f\"Found existing {target_output_file}, skipping due to skip_existing = True, rerun as False to regenerate\"\n",
-    "            )\n",
-    "        actives = pd.read_csv(\n",
-    "            self.dude_unprocessed_path / Path(dude_target) / actives_file,\n",
-    "            sep=seperator,\n",
-    "            header=None,\n",
-    "            names=[\"SMILES\", \"DUDEID\", \"CHEMBLID\"],\n",
-    "        )\n",
-    "        actives[\"Active\"] = 1\n",
-    "        decoys = pd.read_csv(\n",
-    "            self.dude_unprocessed_path / Path(dude_target) / decoys_file,\n",
-    "            sep=seperator,\n",
-    "            header=None,\n",
-    "            names=[\"SMILES\", \"DUDEID\"],\n",
-    "        )\n",
-    "        decoys[\"Active\"] = 0\n",
-    "        df = pd.concat([actives, decoys]).reset_index().drop(columns=\"index\")\n",
-    "        df[\"rdkit_canonical_smiles\"] = tqdm(\n",
-    "            mp.Pool(\n",
-    "                initializer=self._parallel_init_worker_desc_gen_shared_fepops_ob\n",
-    "            ).imap(self._parallel_get_rdkit_cansmi, df.SMILES, chunksize=100),\n",
-    "            desc=f\"Generating {dude_target} benchmark file\",\n",
-    "            total=len(df),\n",
-    "        )\n",
-    "        df.to_csv(target_output_file, index=False)\n",
-    "\n",
-    "    def cache_mols_from_csv(\n",
-    "        self,\n",
-    "        csv_path: Union[Path, str],\n",
-    "        rdkit_canonical_smiles_column_header: str = \"rdkit_canonical_smiles\",\n",
-    "    ):\n",
-    "        \"\"\"Cache mols from a CSV into a db for faster recall later\n",
-    "\n",
-    "        Parameters\n",
-    "        ----------\n",
-    "        csv_path : Union[Path, str]\n",
-    "            Path of CSV file. If None, then all CSV files in the DUDE datasets\n",
-    "            processed path are used.\n",
-    "        rdkit_canonical_smiles_column_header : str, optional\n",
-    "            _description_, by default \"rdkit_canonical_smiles\"\n",
-    "        \"\"\"\n",
-    "\n",
-    "        from fepops.fepops_persistent import get_persistent_fepops_storage_object\n",
-    "\n",
-    "        for csv_path in (\n",
-    "            [Path(csv_path)]\n",
-    "            if csv_path is not None\n",
-    "            else self.dude_processed_path.glob(\"dude_target_*.csv\")\n",
-    "        ):\n",
-    "            df = pd.read_csv(csv_path)\n",
-    "            if df[rdkit_canonical_smiles_column_header].isnull().values.any():\n",
-    "                print(\n",
-    "                    f\"Whilst working on caching {csv_path}, the following mol rows did not contain RDKit canonical SMILES:\"\n",
-    "                )\n",
-    "                print(df[df[rdkit_canonical_smiles_column_header].isnull()])\n",
-    "            smiles = [\n",
-    "                s\n",
-    "                for s in df[rdkit_canonical_smiles_column_header].tolist()\n",
-    "                if not pd.isnull(s)\n",
-    "            ]\n",
-    "            with get_persistent_fepops_storage_object(csv_path.with_suffix(\".db\")) as f:\n",
-    "                f.save_descriptors(smiles)\n",
-    "\n",
+    "from fepops.utils import DudePreprocessor\n",
     "dude_preprocessor=DudePreprocessor()\n",
     "dude_preprocessor()"
    ]
@@ -188,144 +71,134 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[PosixPath('data/dude/processed/diversity_set/dude_target_kif11.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_hivrt.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_cxcr4.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_ampc.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_gcr.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_cp3a4.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_akt1.csv'), PosixPath('data/dude/processed/diversity_set/dude_target_hivpr.csv')]\n",
-      "Working on dude_target_kif11\n"
+      "Working on data/dude/processed/diversity_set/dude_target_kif11.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "  0%|          | 0/6966 [00:00<?, ?it/s]"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "100%|██████████| 6966/6966 [00:05<00:00, 1255.33it/s]\n"
+      "100%|██████████| 6966/6966 [00:02<00:00, 2522.41it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_hivrt\n"
+      "Working on data/dude/processed/diversity_set/dude_target_hivrt.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 19229/19229 [00:13<00:00, 1378.07it/s]\n"
+      "100%|██████████| 19229/19229 [00:07<00:00, 2640.52it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_cxcr4\n"
+      "Working on data/dude/processed/diversity_set/dude_target_cxcr4.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 3446/3446 [00:02<00:00, 1555.36it/s]\n"
+      "100%|██████████| 3446/3446 [00:01<00:00, 2641.72it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_ampc\n"
+      "Working on data/dude/processed/diversity_set/dude_target_ampc.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 2898/2898 [00:01<00:00, 1903.67it/s]\n"
+      "100%|██████████| 2898/2898 [00:00<00:00, 3033.43it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_gcr\n"
+      "Working on data/dude/processed/diversity_set/dude_target_gcr.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      " 66%|██████▌   | 10066/15258 [00:06<00:03, 1631.98it/s][22:53:46] Explicit valence for atom # 12 N, 5, is greater than permitted\n",
-      "[22:53:46] Explicit valence for atom # 12 N, 5, is greater than permitted\n",
-      " 68%|██████▊   | 10394/15258 [00:06<00:03, 1605.42it/s]"
+      " 66%|██████▋   | 10142/15258 [00:04<00:02, 2546.86it/s][22:46:45] Explicit valence for atom # 12 N, 5, is greater than permitted\n",
+      " 70%|██████▉   | 10663/15258 [00:04<00:01, 2571.43it/s]"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Could not parse smiles to a valid molecule, smiles was: N#CC(NC(=O)c1ccccc1)=N1=C(c2ccccc2)C=C(c2ccccc2)C=C1c1ccccc1\n",
-      "Failed to make a molecule from N#CC(NC(=O)c1ccccc1)=N1=C(c2ccccc2)C=C(c2ccccc2)C=C1c1ccccc1\n"
+      "Could not parse smiles to a valid molecule, smiles was: N#CC(NC(=O)c1ccccc1)=N1=C(c2ccccc2)C=C(c2ccccc2)C=C1c1ccccc1\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 15258/15258 [00:10<00:00, 1503.96it/s]\n"
+      "100%|██████████| 15258/15258 [00:06<00:00, 2433.03it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_cp3a4\n"
+      "Working on data/dude/processed/diversity_set/dude_target_cp3a4.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 11970/11970 [00:08<00:00, 1350.43it/s]\n"
+      "100%|██████████| 11970/11970 [00:05<00:00, 2326.70it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_akt1\n"
+      "Working on data/dude/processed/diversity_set/dude_target_akt1.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 16743/16743 [00:11<00:00, 1400.16it/s]\n"
+      "100%|██████████| 16743/16743 [00:06<00:00, 2414.06it/s]\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Working on dude_target_hivpr\n"
+      "Working on data/dude/processed/diversity_set/dude_target_hivpr.csv\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 36286/36286 [00:30<00:00, 1176.21it/s]\n"
+      "100%|██████████| 36286/36286 [00:16<00:00, 2257.40it/s]\n"
      ]
     },
     {
@@ -340,11 +213,11 @@
     {
      "data": {
       "text/plain": [
-       "array([-0.28932319,  0.5166312 ,  0.37458883,  0.99913668, -0.04193182,\n",
-       "        1.03616917,  0.27327129,  0.99839024,  0.09701198,  1.12969387,\n",
-       "        0.23718642,  0.99865705,  0.35968991,  0.6649304 ,  0.4123743 ,\n",
-       "        0.99893657,  5.70852885,  6.3707943 ,  6.47354071,  6.26385429,\n",
-       "        6.19229367,  6.22946713])"
+       "array([-0.28971602,  0.5181022 ,  0.37487135,  0.99922747, -0.04187301,\n",
+       "        1.03382471,  0.27407036,  0.99853436,  0.09725517,  1.12824307,\n",
+       "        0.23735556,  0.99882914,  0.35977538,  0.66653514,  0.41238282,\n",
+       "        0.99902545,  5.71261449,  6.37716992,  6.47293777,  6.26134733,\n",
+       "        6.20354385,  6.23201498])"
       ]
      },
      "metadata": {},
@@ -362,11 +235,11 @@
     {
      "data": {
       "text/plain": [
-       "array([0.35067291, 1.00802116, 0.48380817, 0.02926675, 0.15400475,\n",
-       "       0.86220776, 0.44542581, 0.03999429, 0.16085455, 0.92042695,\n",
-       "       0.42515847, 0.03655217, 0.35778578, 1.36108994, 0.49210665,\n",
-       "       0.03252466, 1.96446927, 2.30792259, 2.5024708 , 2.4155645 ,\n",
-       "       2.29434487, 2.31437527])"
+       "array([0.35110473, 1.00839329, 0.4838859 , 0.02769204, 0.15418035,\n",
+       "       0.86446056, 0.44583626, 0.0381767 , 0.16095862, 0.92079483,\n",
+       "       0.42526185, 0.03413741, 0.35756229, 1.36093993, 0.4921059 ,\n",
+       "       0.0311619 , 1.9668792 , 2.31266486, 2.50699385, 2.41269982,\n",
+       "       2.30018205, 2.31527129])"
       ]
      },
      "metadata": {},
@@ -375,17 +248,15 @@
    ],
    "source": [
     "from pathlib import Path\n",
-    "from fepops import OpenFEPOPS\n",
     "from tqdm import tqdm\n",
     "from fepops.fepops_persistent import get_persistent_fepops_storage_object\n",
     "import numpy as np\n",
     "dude_diversity_set_path=Path(\"data/dude/processed/diversity_set/\")\n",
     "diversity_target_files=list(dude_diversity_set_path.glob(\"dude_target_*.csv\"))\n",
-    "print(diversity_target_files)\n",
     "descriptors=[]\n",
     "for diversity_target in diversity_target_files:\n",
     "    f=get_persistent_fepops_storage_object(diversity_target.with_suffix(\".db\"))\n",
-    "    print(f\"Working on {diversity_target.stem}\")\n",
+    "    print(f\"Working on {diversity_target}\")\n",
     "    for (orig_smi, dude_id, chemblid, active_flag, can_smi) in tqdm([l.strip().split(\",\") for l in open(diversity_target).readlines()[1:] if len(l)>3]):\n",
     "        status, retrieved_descriptors=f.get_fepops(can_smi)\n",
     "        if status.value ==1:\n",
@@ -415,26 +286,23 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Got 8 diversity_set_csv_files : ['dude_target_kif11', 'dude_target_hivrt', 'dude_target_cxcr4', 'dude_target_ampc', 'dude_target_gcr', 'dude_target_cp3a4', 'dude_target_akt1', 'dude_target_hivpr']\n"
+      "Got 1 diversity_set_csv_files : ['dude_target_ampc']\n"
      ]
     }
    ],
    "source": [
-    "import time\n",
     "import numpy as np\n",
     "from tqdm import tqdm\n",
-    "from rdkit.Chem import rdMolDescriptors\n",
     "\n",
-    "from fepops import OpenFEPOPS\n",
     "from dataclasses import dataclass\n",
-    "from typing import Callable, Union, Optional\n",
+    "from typing import Callable, Optional\n",
     "from pathlib import Path\n",
     "import pandas as pd\n",
     "from rdkit import Chem\n",
@@ -443,13 +311,10 @@
     "from rdkit.Chem import AllChem\n",
     "from rdkit.Chem import DataStructs\n",
     "from sklearn.metrics import roc_auc_score\n",
-    "from fepops.fepops_persistent import get_persistent_fepops_storage_object\n",
-    "from typing import Union, Optional\n",
     "from fepops import OpenFEPOPS\n",
     "from fepops.fepops_persistent import get_persistent_fepops_storage_object\n",
     "from fepops.fepops import GetFepopStatusCode\n",
     "import json\n",
-    "import sys\n",
     "\n",
     "open_fepops_object=OpenFEPOPS()\n",
     "\n",
@@ -484,7 +349,7 @@
     "    ),\n",
     "}\n",
     "\n",
-    "diversity_set_csv_files = list(Path(\"data/dude/processed/diversity_set/\").glob(\"dude_target*.csv\"))\n",
+    "diversity_set_csv_files = list(Path(\"data/dude/processed/diversity_set/\").glob(\"dude_target_hiv*.csv\"))\n",
     "\n",
     "# Write all CSVs to SMILES files\n",
     "print(f\"Got {len(diversity_set_csv_files)} diversity_set_csv_files : {[f.stem for f in diversity_set_csv_files]}\")\n",
@@ -505,9 +370,88 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 10,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "                                              SMILES DUDEID  CHEMBLID  Active  \\\n",
+      "0     S(Nc1c(O)cc(C(=O)O)cc1)(c2c(scc2)C(=O)O)(=O)=O    116       NaN       1   \n",
+      "1       S(Nc1ccc(C(=O)O)cc1)(c2ccsc2C(=O)[O-])(=O)=O    113       NaN       1   \n",
+      "2       S(Nc1cccc(c1)C(=O)O)(c2ccsc2C(=O)[O-])(=O)=O    114       NaN       1   \n",
+      "3   S(Nc1ccc(C(=O)O)c(c1)Cl)(c2ccsc2C(=O)[O-])(=O)=O    115       NaN       1   \n",
+      "4  OC(=O)[C@@H](CC1=CC=CC2=CC=CC=C12)N1C(=O)C2=C(...    211       NaN       1   \n",
+      "\n",
+      "                              rdkit_canonical_smiles  \n",
+      "0         O=C(O)c1ccc(NS(=O)(=O)c2ccsc2C(=O)O)c(O)c1  \n",
+      "1         O=C(O)c1ccc(NS(=O)(=O)c2ccsc2C(=O)[O-])cc1  \n",
+      "2         O=C(O)c1cccc(NS(=O)(=O)c2ccsc2C(=O)[O-])c1  \n",
+      "3       O=C(O)c1ccc(NS(=O)(=O)c2ccsc2C(=O)[O-])cc1Cl  \n",
+      "4  O=C(O)c1ccc2c(c1)C(=O)N([C@H](Cc1cccc3ccccc13)...  \n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Caching Morgan 2 descriptors for dude_target_ampc: 2898it [00:00, 6261.77it/s]\n",
+      "Caching MACCS descriptors for dude_target_ampc: 2898it [00:01, 1745.90it/s]\n",
+      "Caching RDKit descriptors for dude_target_ampc: 2898it [00:00, 3455.80it/s]\n",
+      "Caching OpenFEPOPS descriptors for dude_target_ampc: 585it [00:00, 2932.53it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Problem with Cc1ccc([C@H]2[C@@]([O-])([C@H](O)CO)[C@]3(O)[C@@H](c4ccc(C)cc4)[C@H]([O-])[C@@]23[O-])cc1, (<GetFepopStatusCode.FAILED_RETRIEVED_NONE: 4>, None)\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Caching OpenFEPOPS descriptors for dude_target_ampc: 2898it [00:00, 3002.50it/s]\n",
+      "Assessing active recall (AUROC) for Morgan 2: 100%|██████████| 48/48 [00:00<00:00, 135.08it/s]\n",
+      "Assessing active recall (AUROC) for MACCS: 100%|██████████| 48/48 [00:00<00:00, 203.02it/s]\n",
+      "Assessing active recall (AUROC) for RDKit: 100%|██████████| 48/48 [00:00<00:00, 200.67it/s]\n",
+      "Assessing active recall (AUROC) for OpenFEPOPS: 100%|██████████| 48/48 [00:06<00:00,  7.04it/s]"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Writing to  data/dude/processed/diversity_set\n",
+      "  target similarity_method smiles_count actives_count failed_smiles  \\\n",
+      "0   ampc          Morgan 2         2898            48             1   \n",
+      "1   ampc             MACCS         2898            48             1   \n",
+      "2   ampc             RDKit         2898            48             1   \n",
+      "3   ampc        OpenFEPOPS         2898            48             1   \n",
+      "\n",
+      "  failed_active_smiles average_auroc_score median_auroc_score q1_auroc_score  \\\n",
+      "0                    0            0.783629           0.831864       0.315866   \n",
+      "1                    0             0.67331           0.711149       0.365732   \n",
+      "2                    0            0.659837           0.662738       0.349989   \n",
+      "3                    0            0.639413           0.736706        0.22802   \n",
+      "\n",
+      "  q3_auroc_score  \n",
+      "0       0.321303  \n",
+      "1       0.377163  \n",
+      "2       0.361342  \n",
+      "3        0.23266  \n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\n"
+     ]
+    }
+   ],
    "source": [
     "auroc_scores_info_df=pd.DataFrame()\n",
     "\n",
@@ -602,7 +546,7 @@
     "    print(\"Writing to \", csv_file_path.parent)\n",
     "\n",
     "    json.dump(auroc_scores,open(csv_file_path.parent/Path(f\"res_scores_{csv_file_path.stem}.json\"),\"w\"))\n",
-    "auroc_scores_info_df.to_csv(csv_file_path.parent/Path(f\"res_df_{csv_file_path.stem}.csv\"))\n",
+    "    auroc_scores_info_df.to_csv(csv_file_path.parent/Path(f\"res_df_{csv_file_path.stem}.csv\"))\n",
     "print(auroc_scores_info_df)\n",
     "    "
    ]
@@ -616,7 +560,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -659,56 +603,56 @@
        "      <td>0.835717</td>\n",
        "      <td>0.833538</td>\n",
        "      <td>0.741306</td>\n",
-       "      <td>0.828947</td>\n",
+       "      <td>0.831254</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ampc</th>\n",
        "      <td>0.783629</td>\n",
        "      <td>0.659837</td>\n",
        "      <td>0.673310</td>\n",
-       "      <td>0.639115</td>\n",
+       "      <td>0.639413</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>cp3a4</th>\n",
        "      <td>0.602774</td>\n",
        "      <td>0.613335</td>\n",
        "      <td>0.581545</td>\n",
-       "      <td>0.649807</td>\n",
+       "      <td>0.647462</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>cxcr4</th>\n",
        "      <td>0.697226</td>\n",
        "      <td>0.592485</td>\n",
        "      <td>0.854251</td>\n",
-       "      <td>0.898971</td>\n",
+       "      <td>0.899027</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>gcr</th>\n",
-       "      <td>0.670138</td>\n",
-       "      <td>0.708350</td>\n",
-       "      <td>0.665839</td>\n",
-       "      <td>0.616173</td>\n",
+       "      <td>0.670082</td>\n",
+       "      <td>0.708346</td>\n",
+       "      <td>0.665830</td>\n",
+       "      <td>0.616228</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>hivpr</th>\n",
-       "      <td>0.779689</td>\n",
-       "      <td>0.759308</td>\n",
-       "      <td>0.681237</td>\n",
-       "      <td>0.677882</td>\n",
+       "      <td>0.779684</td>\n",
+       "      <td>0.759303</td>\n",
+       "      <td>0.681231</td>\n",
+       "      <td>0.678182</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>hivrt</th>\n",
-       "      <td>0.651011</td>\n",
-       "      <td>0.660054</td>\n",
-       "      <td>0.669939</td>\n",
-       "      <td>0.583981</td>\n",
+       "      <td>0.651023</td>\n",
+       "      <td>0.660070</td>\n",
+       "      <td>0.669954</td>\n",
+       "      <td>0.582240</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>kif11</th>\n",
-       "      <td>0.763058</td>\n",
-       "      <td>0.672460</td>\n",
-       "      <td>0.667887</td>\n",
-       "      <td>0.713152</td>\n",
+       "      <td>0.763057</td>\n",
+       "      <td>0.672469</td>\n",
+       "      <td>0.667852</td>\n",
+       "      <td>0.713142</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -717,14 +661,14 @@
       "text/plain": [
        "similarity_method  Morgan 2     RDKit     MACCS  OpenFEPOPS\n",
        "target                                                     \n",
-       "akt1               0.835717  0.833538  0.741306    0.828947\n",
-       "ampc               0.783629  0.659837  0.673310    0.639115\n",
-       "cp3a4              0.602774  0.613335  0.581545    0.649807\n",
-       "cxcr4              0.697226  0.592485  0.854251    0.898971\n",
-       "gcr                0.670138  0.708350  0.665839    0.616173\n",
-       "hivpr              0.779689  0.759308  0.681237    0.677882\n",
-       "hivrt              0.651011  0.660054  0.669939    0.583981\n",
-       "kif11              0.763058  0.672460  0.667887    0.713152"
+       "akt1               0.835717  0.833538  0.741306    0.831254\n",
+       "ampc               0.783629  0.659837  0.673310    0.639413\n",
+       "cp3a4              0.602774  0.613335  0.581545    0.647462\n",
+       "cxcr4              0.697226  0.592485  0.854251    0.899027\n",
+       "gcr                0.670082  0.708346  0.665830    0.616228\n",
+       "hivpr              0.779684  0.759303  0.681231    0.678182\n",
+       "hivrt              0.651023  0.660070  0.669954    0.582240\n",
+       "kif11              0.763057  0.672469  0.667852    0.713142"
       ]
      },
      "metadata": {},
@@ -744,7 +688,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -794,7 +738,7 @@
        "      <th>akt1</th>\n",
        "      <td>0.741306</td>\n",
        "      <td>0.835717</td>\n",
-       "      <td>0.828947</td>\n",
+       "      <td>0.831254</td>\n",
        "      <td>0.833538</td>\n",
        "      <td>16743</td>\n",
        "      <td>293</td>\n",
@@ -805,7 +749,7 @@
        "      <th>ampc</th>\n",
        "      <td>0.673310</td>\n",
        "      <td>0.783629</td>\n",
-       "      <td>0.639115</td>\n",
+       "      <td>0.639413</td>\n",
        "      <td>0.659837</td>\n",
        "      <td>2898</td>\n",
        "      <td>48</td>\n",
@@ -816,7 +760,7 @@
        "      <th>cp3a4</th>\n",
        "      <td>0.581545</td>\n",
        "      <td>0.602774</td>\n",
-       "      <td>0.649807</td>\n",
+       "      <td>0.647462</td>\n",
        "      <td>0.613335</td>\n",
        "      <td>11970</td>\n",
        "      <td>170</td>\n",
@@ -827,7 +771,7 @@
        "      <th>cxcr4</th>\n",
        "      <td>0.854251</td>\n",
        "      <td>0.697226</td>\n",
-       "      <td>0.898971</td>\n",
+       "      <td>0.899027</td>\n",
        "      <td>0.592485</td>\n",
        "      <td>3446</td>\n",
        "      <td>40</td>\n",
@@ -836,46 +780,46 @@
        "    </tr>\n",
        "    <tr>\n",
        "      <th>gcr</th>\n",
-       "      <td>0.665839</td>\n",
-       "      <td>0.670138</td>\n",
-       "      <td>0.616173</td>\n",
-       "      <td>0.708350</td>\n",
+       "      <td>0.665830</td>\n",
+       "      <td>0.670082</td>\n",
+       "      <td>0.616228</td>\n",
+       "      <td>0.708346</td>\n",
        "      <td>15258</td>\n",
        "      <td>258</td>\n",
-       "      <td>4</td>\n",
+       "      <td>10</td>\n",
        "      <td>0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>hivpr</th>\n",
-       "      <td>0.681237</td>\n",
-       "      <td>0.779689</td>\n",
-       "      <td>0.677882</td>\n",
-       "      <td>0.759308</td>\n",
+       "      <td>0.681231</td>\n",
+       "      <td>0.779684</td>\n",
+       "      <td>0.678182</td>\n",
+       "      <td>0.759303</td>\n",
        "      <td>36286</td>\n",
        "      <td>535</td>\n",
-       "      <td>1</td>\n",
+       "      <td>2</td>\n",
        "      <td>1</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>hivrt</th>\n",
-       "      <td>0.669939</td>\n",
-       "      <td>0.651011</td>\n",
-       "      <td>0.583981</td>\n",
-       "      <td>0.660054</td>\n",
+       "      <td>0.669954</td>\n",
+       "      <td>0.651023</td>\n",
+       "      <td>0.582240</td>\n",
+       "      <td>0.660070</td>\n",
        "      <td>19229</td>\n",
        "      <td>338</td>\n",
-       "      <td>9</td>\n",
+       "      <td>8</td>\n",
        "      <td>0</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>kif11</th>\n",
-       "      <td>0.667887</td>\n",
-       "      <td>0.763058</td>\n",
-       "      <td>0.713152</td>\n",
-       "      <td>0.672460</td>\n",
+       "      <td>0.667852</td>\n",
+       "      <td>0.763057</td>\n",
+       "      <td>0.713142</td>\n",
+       "      <td>0.672469</td>\n",
        "      <td>6966</td>\n",
        "      <td>116</td>\n",
-       "      <td>0</td>\n",
+       "      <td>1</td>\n",
        "      <td>0</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
@@ -885,14 +829,14 @@
       "text/plain": [
        "           MACCS  Morgan 2  OpenFEPOPS     RDKit  smiles_count  actives_count  \\\n",
        "target                                                                          \n",
-       "akt1    0.741306  0.835717    0.828947  0.833538         16743            293   \n",
-       "ampc    0.673310  0.783629    0.639115  0.659837          2898             48   \n",
-       "cp3a4   0.581545  0.602774    0.649807  0.613335         11970            170   \n",
-       "cxcr4   0.854251  0.697226    0.898971  0.592485          3446             40   \n",
-       "gcr     0.665839  0.670138    0.616173  0.708350         15258            258   \n",
-       "hivpr   0.681237  0.779689    0.677882  0.759308         36286            535   \n",
-       "hivrt   0.669939  0.651011    0.583981  0.660054         19229            338   \n",
-       "kif11   0.667887  0.763058    0.713152  0.672460          6966            116   \n",
+       "akt1    0.741306  0.835717    0.831254  0.833538         16743            293   \n",
+       "ampc    0.673310  0.783629    0.639413  0.659837          2898             48   \n",
+       "cp3a4   0.581545  0.602774    0.647462  0.613335         11970            170   \n",
+       "cxcr4   0.854251  0.697226    0.899027  0.592485          3446             40   \n",
+       "gcr     0.665830  0.670082    0.616228  0.708346         15258            258   \n",
+       "hivpr   0.681231  0.779684    0.678182  0.759303         36286            535   \n",
+       "hivrt   0.669954  0.651023    0.582240  0.660070         19229            338   \n",
+       "kif11   0.667852  0.763057    0.713142  0.672469          6966            116   \n",
        "\n",
        "        failed_smiles  failed_active_smiles  \n",
        "target                                       \n",
@@ -900,10 +844,10 @@
        "ampc                1                     0  \n",
        "cp3a4               4                     0  \n",
        "cxcr4               0                     0  \n",
-       "gcr                 4                     0  \n",
-       "hivpr               1                     1  \n",
-       "hivrt               9                     0  \n",
-       "kif11               0                     0  "
+       "gcr                10                     0  \n",
+       "hivpr               2                     1  \n",
+       "hivrt               8                     0  \n",
+       "kif11               1                     0  "
       ]
      },
      "metadata": {},
diff --git a/README.md b/README.md
index 33c5047..c0762e9 100644
--- a/README.md
+++ b/README.md
@@ -15,27 +15,27 @@ OpenFEPOPS has been benchmarked alongside some commonly used fingerprints agains
 
 | Target | Morgan 2 | MACCS | RDKit |OpenFEPOPS|
 |--------:|----------:|-------:|-------:|------------:|
-|akt1  |0.836|0.741 |0.833 |0.829|
+|akt1  |0.836|0.741 |0.833 |0.831|
 |ampc  |0.784|0.673|0.660 |0.639|
-|cp3a4 |0.603|0.582|0.613 |0.650|
+|cp3a4 |0.603|0.582|0.613 |0.647|
 |cxcr4 |0.697|0.854|0.592 |0.899|
 |gcr   |0.670|0.666|0.708  |0.616|
 |hivpr |0.780|0.681|0.759 |0.678|
-|hivrt |0.651|0.670 |0.660 |0.584|
+|hivrt |0.651|0.670 |0.660 |0.582|
 |kif11 |0.763|0.668 |0.672  |0.713|
 
 *Macro averaged AUROC scores by target and molecular similarity technique for the DUDE diversity set. Across all datasets, 19 small molecules out of 112,796 were excluded from analysis mainly due to issues in parsing to valid structures using RDKit.  See 'Explore_DUDE_diversity_set.ipynb' for further information*
 
 ## Method description
 Whilst OpenFEPOPS has included functionality for descriptor caching and profiling of libraries, the core functionality of the package is descriptor generation and scoring.
-### _Descriptor generation:_
+### _Descriptor generation_
 
 1. Tautomer enumeration
     - For a given small molecule, OpenFEPOPS uses RDKit to iterate over molecular tautomers. By default, there is no limit to the number of recoverable tautomers but a limit may be imposed which may be necessary if adapting the OpenFEPOPS code to large macromolecules and not just small molecules.
 2.  Conformer enumeration
     - For each tautomer, up to 1024 conformers are sampled by either complete enumeration of rotatable bond states (at the literature reported optimum increment of 90 degrees) if there are five or less rotatable bonds, or through random sampling of 1024 possible states if there are more than 5 rotatable bonds.
 3.  Defining feature points
-    - The KMeans algorithm is applied to each conformer of each tautomer to identify four (by default) representative or central points, into which the atomic information of neighbouring atoms is collapsed. As standard, the atomic properties of charge, logP, hydrogen bond donor, and hydrogen bond acceptor status are collapsed into four feature points per unique tautomer conformation. These feature points are encoded to 22 numeric values (a FEPOP) comprising four points, each with four properties, and six pairwise distances between these points. With many FEPOPS descriptors collected from a single molecule through tautomer and conformer enumeration, this set of representative FEPOPS should capture every possible state of the original molecule.
+    - The KMeans algorithm is applied to each conformer of each tautomer to identify four (by default) representative or central points, into which the atomic information of neighbouring atoms is collapsed. As standard, the atomic properties of charge, logP, hydrogen bond donor, and hydrogen bond acceptor status are collapsed into four feature points per unique tautomer conformation.The RDKit package is used to calculate these properties with the iterative Gasteiger charges algorithm applied to assign atomic charges, the Crippen method used to assign atomic logP contributions, and hydrogen bond acceptors and donors identified with appropriate SMARTS substructure queries.These feature points are encoded to 22 numeric values (a FEPOP) comprising four points, each with four properties, and six pairwise distances between these points. With many FEPOPS descriptors collected from a single molecule through tautomer and conformer enumeration, this set of representative FEPOPS should capture every possible state of the original molecule.
 4.  Selection of diverse FEPOPS
     - From the collection of FEPOPS derived from every tautomer conformation of a molecule, the K-Medoid algorithm is applied to identify seven (by default) diverse FEPOPS which are thought to best capture a fuzzy representation of the molecule. These seven FEPOPS each comprise 22 descriptors each, totalling 154 32-bit floating point numbers or 616 bytes.
 
@@ -72,7 +72,22 @@ Or via conda :
 conda env create -f environment.yml
 ```
 
-### Requirements:
+Requirements for developing and contributing to OpenFEPOPS can be installed from source:
+```console
+git clone https://github.com/JustinYKC/FEPOPS.git
+cd FEPOPS
+pip install .[dev]
+```
+
+### Testing
+Tests may be run using the pytest framework. With the development dependencies installed (see above), simply running:
+
+```console
+pytest
+```
+
+will run the built in tests.
+### Requirements
 OpenFEPOPS requires the following packages:
 - rdkit (>=2021.9.5.1)
 - numpy (>=1.22.4)
@@ -88,7 +103,7 @@ OpenFEPOPS requires the following packages:
 ## Usage
 OpenFEPOPS has been uploaded to the Python Packaging Index under the name 'fepops' and as such is installable using the pip package manager and the command 'pip install fepops'. With the package installed, entrypoints are used to expose commonly used OpenFEPOPS tasks such as descriptor generation and calculation on molecular similarity, enabling simple command line access without the need to explicitly invoke a Python interpreter. Whilst OpenFEPOPS may be used solely via the command line interface, a robust API is available and may be used within other programs or integrated into existing pipelines to enable more complex workflows.  API documentation is available at https://justinykc.github.io/FEPOPS.
 
-### Command line usage:
+### Command line usage
 With OpenFEPOPS installed to a Python environment, entrypoints in the code and a command line interface control object allow omission of the python interpreter program when calling common OpenFEPOPS tasks. This allows command line usage of the form:
 
 ```console
@@ -109,6 +124,14 @@ fepops --database_file=<DB_file_path> <subcommand> <arguments>
 ```
 Depending on the supplied file extension of <DB_file_path>, either a SQLite file (".sqlite", ".sqlite3", ".db", ".db3", ".s3db", or ".sl3" file extensions) is created or loaded for the appending of data, or a JSON file is used (when the extension is ".json").
 
+Additionally, a file to which JSON formatted ouput should be written to can be supplied with the --json_file flag when using the calc_sim and get_fepops commands as follows:
+
+```console
+fepops --database_file=<DB_file_path> --json_file=<JSON_file_path> get_fepops <arguments>
+fepops --database_file=<DB_file_path> --json_file=<JSON_file_path> calc_sim <arguments>
+```
+where <JSON_file_path> is the path to where JSON formatted output should be written.
+
 Running OpenFEPOPS with the -h switch prints help to the terminal and it may also be used after sub-command switches to inspect required arguments.
 
 #### Example: Calculating molecular similarity between two molecules (ibuprofen and diclofenac)
@@ -117,10 +140,68 @@ fepops calc_sim "CC(Cc1ccc(cc1)C(C(=O)O)C)C" "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl"
 ```
 Note the use of quotes around the smiles strings. This is required as BASH and other shells will try to parse the brackets often present in smiles to denote branching.
 
+Output may be directed to a JSON file as follows:
+
+
+```console
+fepops --json_file="example_fepops_score.json" calc_sim "CC(Cc1ccc(cc1)C(C(=O)O)C)C" "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl"
+```
+
+The file 'example_fepops_score.json' then contains output of the form:
+
+```javascript
+{
+    "SMI1": "CC(Cc1ccc(cc1)C(C(=O)O)C)C",
+    "SMI2": "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl",
+    "FEPOPS Similarity Score": 0.14135371590471596
+}
+```
+
 #### Example: Get the FEPOPS descriptors for diclofenac
 ```console
 fepops get_fepops "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl"
 ```
+Again, output may be directed to a JSON file as follows:
+
+
+```console
+fepops --json_file="example_diclofenac_fepops_descriptors.json" get_fepops "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl"
+```
+
+The file 'example_fepops_score.json' then contains output of the form below. For brevity, output is truncated at the first of 7 FEPOPS, truncation marked by TRUNCATION:
+
+```javascript
+{
+    "SMILES": "OC(=O)Cc1ccccc1Nc1c(Cl)cccc1Cl",
+    "FepopStatusCode": "SUCCESS",
+    "Fepops": [
+        [
+            -0.5684914743317748,
+            1.3061,
+            1.0,
+            1.0,
+            0.02878293161743492,
+            0.8432999999999999,
+            0.0,
+            1.0,
+            0.1949443985093149,
+            1.0734,
+            0.0,
+            1.0,
+            0.3447641442050251,
+            1.1413,
+            0.0,
+            1.0,
+            2.52506114069371,
+            4.753703619492762,
+            3.216333750680169,
+            4.563318239587458,
+            2.9657848574016255,
+            7.096648170644156
+        ],
+    TRUNCATION
+}
+```
 
 #### Example: pregenerating descriptors for an in-house compound archive
 With a SMILES file called 'inhouse_compounds.smi', we may pre-generate their FEPOPS descriptors for faster use and comparison against new molecules.
diff --git a/docsource/conf.py b/docsource/conf.py
index 645451a..90b42bf 100644
--- a/docsource/conf.py
+++ b/docsource/conf.py
@@ -16,6 +16,7 @@
 # https://www.lieret.net/2021/05/20/include-readme-sphinx/
 # to also rewrite the path of Figure1.png
 import pathlib
+
 readme_path = pathlib.Path(__file__).parent.resolve().parent / "README.md"
 readme_target = pathlib.Path(__file__).parent / "readme.md"
 with readme_target.open("w") as outf:
@@ -30,9 +31,9 @@
     lines = []
     for line in readme_path.read_text().split("\n")[1:]:
         if line.startswith("#"):
-            line=line[1:]
+            line = line[1:]
         if "Figure1.png" in line:
-            line=line.replace("Figure1.png", "_static/Figure1.png")
+            line = line.replace("Figure1.png", "_static/Figure1.png")
         lines.append(line)
     outf.write("\n".join(lines))
 # End of README.md importer/parser
@@ -40,29 +41,31 @@
 
 # Copy Figure1.png to docsource/_static
 Path("_static").mkdir(exist_ok=True)
-shutil.copyfile("../Figure1.png", Path("_static")/"Figure1.png")
+shutil.copyfile("../Figure1.png", Path("_static") / "Figure1.png")
 
 project = 'OpenFEPOPS'
 copyright = '2023, Yan-Kai Chen'
 author = 'Yan-Kai Chen'
-release = '1.4.0'
+release = '1.6.0'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
 
-extensions = ["sphinx.ext.todo", "sphinx.ext.viewcode", "sphinx.ext.autodoc", "sphinx.ext.napoleon", "myst_parser"]
+extensions = [
+    "sphinx.ext.todo",
+    "sphinx.ext.viewcode",
+    "sphinx.ext.autodoc",
+    "sphinx.ext.napoleon",
+    "myst_parser",
+]
 
 templates_path = ['_templates']
 exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
 
 
-
 # -- Options for HTML output -------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
 
 html_theme = 'sphinx_rtd_theme'
 html_static_path = ['_static']
 source_suffix = ['.rst', '.md']
-
-
-
diff --git a/docsource/fepops.rst b/docsource/fepops.rst
index d5c76de..bee13f2 100644
--- a/docsource/fepops.rst
+++ b/docsource/fepops.rst
@@ -8,26 +8,11 @@ Subpackages
    :maxdepth: 4
 
    fepops.fepops_persistent
+   fepops.utils
 
 Submodules
 ----------
 
-fepops.benchmark module
------------------------
-
-.. automodule:: fepops.benchmark
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
-fepops.dude\_preprocessor module
---------------------------------
-
-.. automodule:: fepops.dude_preprocessor
-   :members:
-   :undoc-members:
-   :show-inheritance:
-
 fepops.fepops module
 --------------------
 
diff --git a/docsource/index.rst b/docsource/index.rst
index 7d3000c..21e627c 100644
--- a/docsource/index.rst
+++ b/docsource/index.rst
@@ -17,20 +17,19 @@ valuable small molecules. This implementation was recreated following the origin
 
 Source code is available on the GitHub repo here: `OpenFEPOPS GitHub Repository <https://github.com/JustinYKC/FEPOPS>`_
 
-An online version of the GitHub README.md can be found here which documents how to use the package for
-calculation of descriptors and subsequent molecular similarity calculations and is available here: `readme`
+An online version of the GitHub README.md can be found broken down into sections below, detailing how to use the package for calculation of descriptors and subsequent molecular similarity calculations.
 
 .. toctree::
    :maxdepth: 4
-   :caption: API documentation:
+   :caption: Further information:
 
-   modules
+   readme
 
 .. toctree::
    :maxdepth: 4
-   :caption: Further information:
+   :caption: API documentation:
 
-   readme
+   modules
 
 Indices and tables
 ==================
diff --git a/environment.yml b/environment.yml
index 1c742e7..f6a22d5 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,49 +1,185 @@
 name: openfepops
 channels:
+  - conda-forge
   - defaults
 dependencies:
-  - _libgcc_mutex=0.1=main
-  - ca-certificates=2023.05.30=h06a4308_0
+  - _libgcc_mutex=0.1=conda_forge
+  - _openmp_mutex=4.5=2_kmp_llvm
+  - blas=1.0=mkl
+  - boost=1.78.0=py311h59ea3da_4
+  - boost-cpp=1.78.0=h6582d0a_3
+  - bottleneck=1.3.5=py311hbed6279_0
+  - brotli=1.0.9=h166bdaf_8
+  - brotli-bin=1.0.9=h166bdaf_8
+  - bzip2=1.0.8=h7b6447c_0
+  - ca-certificates=2023.5.7=hbcca054_0
+  - cairo=1.16.0=hb05425b_3
+  - certifi=2023.5.7=pyhd8ed1ab_0
+  - contourpy=1.0.7=py311ha3edf6b_0
+  - cycler=0.11.0=pyhd8ed1ab_0
+  - expat=2.5.0=hcb278e6_1
+  - fontconfig=2.14.1=h4c34cd2_2
+  - fonttools=4.39.4=py311h459d7ec_0
+  - freetype=2.12.1=hca18f0e_1
+  - gettext=0.21.1=h27087fc_0
+  - glib=2.76.3=hfc55251_0
+  - glib-tools=2.76.3=hfc55251_0
+  - greenlet=2.0.2=py311hb755f60_1
+  - icu=72.1=hcb278e6_0
+  - intel-openmp=2023.1.0=hdb19cb5_46305
+  - kiwisolver=1.4.4=py311h4dd048b_1
+  - lcms2=2.15=haa2dc70_1
   - ld_impl_linux-64=2.38=h1181459_1
-  - libffi=3.3=he6710b0_2
-  - libgcc-ng=9.1.0=hdf63c60_0
-  - libstdcxx-ng=9.1.0=hdf63c60_0
-  - ncurses=6.3=h7f8727e_2
-  - openssl=1.1.1v=h7f8727e_0
-  - pip=23.2.1=py39h06a4308_0
-  - python=3.9.12=h12debd9_1
-  - readline=8.1.2=h7f8727e_1
-  - setuptools=68.0.0=py39h06a4308_0
-  - sqlite=3.38.5=hc218d9a_0
+  - lerc=4.0.0=h27087fc_0
+  - libbrotlicommon=1.0.9=h166bdaf_8
+  - libbrotlidec=1.0.9=h166bdaf_8
+  - libbrotlienc=1.0.9=h166bdaf_8
+  - libdeflate=1.18=h0b41bf4_0
+  - libexpat=2.5.0=hcb278e6_1
+  - libffi=3.4.4=h6a678d5_0
+  - libgcc-ng=12.2.0=h65d4601_19
+  - libglib=2.76.3=hebfc3b9_0
+  - libiconv=1.17=h166bdaf_0
+  - libjpeg-turbo=2.1.5.1=h0b41bf4_0
+  - libpng=1.6.39=h753d276_0
+  - libstdcxx-ng=12.2.0=h46fd767_19
+  - libtiff=4.5.0=ha587672_6
+  - libuuid=1.41.5=h5eee18b_0
+  - libwebp-base=1.3.0=h0b41bf4_0
+  - libxcb=1.15=h0b41bf4_0
+  - libxml2=2.10.4=hfdac1af_0
+  - libzlib=1.2.13=h166bdaf_4
+  - llvm-openmp=16.0.4=h4dfa4b3_0
+  - matplotlib-base=3.7.1=py311h8597a09_0
+  - mkl=2023.1.0=h6d00ec8_46342
+  - mkl-service=2.4.0=py311h5eee18b_1
+  - mkl_fft=1.3.6=py311ha02d727_1
+  - mkl_random=1.2.2=py311ha02d727_1
+  - munkres=1.1.4=pyh9f0ad1d_0
+  - ncurses=6.4=h6a678d5_0
+  - numexpr=2.8.4=py311h65dcdc2_1
+  - numpy=1.24.3=py311h08b1b3b_1
+  - numpy-base=1.24.3=py311hf175353_1
+  - openjpeg=2.5.0=hfec8fc6_2
+  - openssl=1.1.1t=h0b41bf4_0
+  - packaging=23.1=pyhd8ed1ab_0
+  - pandas=1.5.3=py311hba01205_0
+  - pcre2=10.40=hc3806b6_0
+  - pillow=9.5.0=py311h0b84326_1
+  - pip=23.0.1=py311h06a4308_0
+  - pixman=0.40.0=h36c2ea0_0
+  - pthread-stubs=0.4=h36c2ea0_1001
+  - pycairo=1.23.0=py311hdece40b_0
+  - pyparsing=3.0.9=pyhd8ed1ab_0
+  - python=3.11.3=h7a1cb2a_0
+  - python-dateutil=2.8.2=pyhd3eb1b0_0
+  - python_abi=3.11=2_cp311
+  - pytz=2022.7=py311h06a4308_0
+  - readline=8.2=h5eee18b_0
+  - reportlab=3.6.13=py311hdc8528d_0
+  - setuptools=66.0.0=py311h06a4308_0
+  - six=1.16.0=pyhd3eb1b0_1
+  - sqlalchemy=2.0.15=py311h459d7ec_0
+  - sqlite=3.41.2=h5eee18b_0
+  - tbb=2021.8.0=hdb19cb5_0
   - tk=8.6.12=h1ccaba5_0
+  - typing-extensions=4.6.2=hd8ed1ab_0
+  - typing_extensions=4.6.2=pyha770c72_0
   - tzdata=2023c=h04d1e81_0
-  - wheel=0.38.4=py39h06a4308_0
-  - xz=5.2.5=h7f8727e_1
-  - zlib=1.2.12=h7f8727e_2
+  - wheel=0.38.4=py311h06a4308_0
+  - xorg-libxau=1.0.11=hd590300_0
+  - xorg-libxdmcp=1.1.3=h7f98852_0
+  - xz=5.4.2=h5eee18b_0
+  - zlib=1.2.13=h166bdaf_4
+  - zstd=1.5.2=h3eb15da_6
   - pip:
-    - build==0.10.0
-    - fast-pytorch-kmeans==0.1.9
-    - fepops
-    - fire==0.5.0
-    - joblib==1.3.2
-    - numpy==1.22.4
-    - packaging==23.1
-    - pandas==1.5.0
-    - pillow==10.0.0
-    - pynvml==11.5.0
-    - pyproject_hooks==1.0.0
-    - python-dateutil==2.8.2
-    - pytz==2023.3
-    - rdkit-pypi==2021.9.5.1
-    - scikit-learn==1.3.0
-    - scipy==1.7.0
-    - six==1.16.0
-    - termcolor==2.3.0
-    - threadpoolctl==3.2.0
-    - tomli==2.0.1
-    - torch==1.7.1+cu101
-    - torchaudio==0.7.2
-    - torchvision==0.8.2+cu101
-    - tqdm==4.48.0
-    - typing_extensions==4.7.1
+      - astroid==2.15.5
+      - asttokens==2.2.1
+      - attrs==21.4.0
+      - backcall==0.2.0
+      - black==23.3.0
+      - charset-normalizer==3.1.0
+      - chembl-webresource-client==0.10.8
+      - click==8.1.3
+      - cmake==3.26.3
+      - comm==0.1.3
+      - debugpy==1.6.7
+      - decorator==5.1.1
+      - dill==0.3.6
+      - easydict==1.10
+      - executing==1.2.0
+      - fast-pytorch-kmeans==0.2.0.1
+      - fepops
+      - filelock==3.12.0
+      - fire==0.5.0
+      - h5py==3.8.0
+      - idna==3.4
+      - iniconfig==2.0.0
+      - ipykernel==6.23.1
+      - ipython==8.13.2
+      - isort==5.12.0
+      - itsdangerous==2.1.2
+      - jedi==0.18.2
+      - jinja2==3.1.2
+      - joblib==1.2.0
+      - jupyter-client==8.2.0
+      - jupyter-core==5.3.0
+      - lazy-object-proxy==1.9.0
+      - lit==16.0.5
+      - markdown-it-py==2.2.0
+      - markupsafe==2.1.2
+      - matplotlib-inline==0.1.6
+      - mccabe==0.7.0
+      - mdurl==0.1.2
+      - mpmath==1.3.0
+      - mypy-extensions==1.0.0
+      - nest-asyncio==1.5.6
+      - networkx==3.1
+      - nvidia-cublas-cu11==11.10.3.66
+      - nvidia-cuda-cupti-cu11==11.7.101
+      - nvidia-cuda-nvrtc-cu11==11.7.99
+      - nvidia-cuda-runtime-cu11==11.7.99
+      - nvidia-cudnn-cu11==8.5.0.96
+      - nvidia-cufft-cu11==10.9.0.58
+      - nvidia-curand-cu11==10.2.10.91
+      - nvidia-cusolver-cu11==11.4.0.1
+      - nvidia-cusparse-cu11==11.7.4.91
+      - nvidia-nccl-cu11==2.14.3
+      - nvidia-nvtx-cu11==11.7.91
+      - parso==0.8.3
+      - pathspec==0.11.1
+      - pexpect==4.8.0
+      - pickleshare==0.7.5
+      - platformdirs==3.5.1
+      - pluggy==1.0.0
+      - prompt-toolkit==3.0.38
+      - psutil==5.9.5
+      - ptyprocess==0.7.0
+      - pure-eval==0.2.2
+      - pygments==2.15.1
+      - pylint==2.17.4
+      - pynvml==11.5.0
+      - pytest==7.3.1
+      - pyyaml==6.0
+      - pyzmq==25.1.0
+      - rdkit==2023.3.3
+      - requests==2.31.0
+      - requests-cache==0.7.5
+      - rich==13.4.0
+      - scikit-learn==1.3.0
+      - scipy==1.10.1
+      - stack-data==0.6.2
+      - sympy==1.12
+      - termcolor==2.3.0
+      - threadpoolctl==3.1.0
+      - tomlkit==0.11.8
+      - torch==2.0.1
+      - tornado==6.3.2
+      - tqdm==4.65.0
+      - traitlets==5.9.0
+      - triton==2.0.0
+      - url-normalize==1.4.3
+      - urllib3==2.0.2
+      - wcwidth==0.2.6
+      - wrapt==1.15.0
 
diff --git a/mol_filter.py b/mol_filter.py
deleted file mode 100644
index e83b004..0000000
--- a/mol_filter.py
+++ /dev/null
@@ -1,226 +0,0 @@
-from rdkit import Chem
-from rdkit.Chem import rdMolDescriptors
-from rdkit.Chem.MolStandardize import rdMolStandardize
-import pandas as pd
-
-
-class Filter:
-    """Filter object allows filtering of molecules by property
-
-    Filter object allowing filtering by min atoms, max rings, and max rotatable bonds in a molecule.
-    """
-
-    def __init__(
-        self, min_atoms: int = 4, max_rings: int = 9, max_n_rot: int = 40
-    ) -> None:
-        self.min_atoms = min_atoms
-        self.max_rings = max_rings
-        self.max_n_rot = max_n_rot
-
-    def __call__(self, mol: Chem.rdchem.Mol) -> Chem.rdchem.Mol or None:
-        """Apply all filters
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-
-        Returns
-        -------
-        filter_result
-          Return a Rdkit mol object if the input molecule meets all criteria, otherwise return None.
-        """
-        filter_result = self.filter_mol(mol)
-        return filter_result
-
-    def filter_mol(self, mol: Chem.rdchem.Mol) -> Chem.rdchem.Mol or None:
-        """Apply all filters
-
-        Perform four filters (minimum number of atoms, maximum number of rings, maximum number of rotatable bonds,
-        and salt/ion) in sequence to obtain the final molecules as required.
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-
-        Returns
-        -------
-        Chem.rdchem.Mol or None
-          Return a Rdkit mol object if the input molecule meets all criteria, otherwise return None.
-        """
-        mol = self.atom_num(mol, self.min_atoms)
-        if mol is None:
-            return None
-        else:
-            mol = self.ring_num(mol, self.max_rings)
-            if mol is None:
-                return None
-            else:
-                mol = self.rotation_bond_num(mol, self.max_n_rot)
-                if mol is None:
-                    return None
-                else:
-                    mol = self.remove_salt(mol)
-                    if mol is None:
-                        return None
-                    else:
-                        return mol
-
-    def read_as_mol(self, query_smiles: str) -> Chem.rdchem.Mol:
-        """Read SMILES strings from a file
-
-        Read molecular SMILES strings as Rdkit mol objects from a file.
-
-        Parameters
-        ----------
-        query_smiles : str
-          A molecular SMILES string.
-
-        Returns
-        -------
-        Chem.rdchem.Mol
-          Return a Rdkit mol object.
-        """
-        return Chem.MolFromSmiles(query_smiles)
-
-    def atom_num(
-        self, mol: Chem.rdchem.Mol, cutoff: int = 4
-    ) -> Chem.rdchem.Mol or None:
-        """Filter molecules by number of atoms
-
-        Filter out molecules by the number of atoms.
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-        cutoff : int
-          The criterion of atom numbers for filtering. By default 4.
-
-        Returns
-        -------
-        Chem.rdchem.Mol or None
-          Return a Rdkit mol object if the input molecule meets the criterion, otherwise return None.
-        """
-        atom_num = mol.GetNumAtoms()
-        if atom_num > cutoff:
-            return mol
-        else:
-            return None
-
-    def ring_num(
-        self, mol: Chem.rdchem.Mol, cutoff: int = 9
-    ) -> Chem.rdchem.Mol or None:
-        """Fliter molecules by number of rings
-
-        Fliter out molecules by the number of rings.
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-        cutoff : int
-          The criterion of ring numbers for filtering. By default, 9.
-
-        Returns
-        -------
-        Chem.rdchem.Mol or None
-          Return a Rdkit mol object if the input molecule meets the criterion, otherwise return None.
-        """
-        ring_num = mol.GetRingInfo().NumRings()
-        if ring_num < cutoff:
-            return mol
-        else:
-            return None
-
-    def rotation_bond_num(
-        self, mol: Chem.rdchem.Mol, cutoff: int = 40
-    ) -> Chem.rdchem.Mol or None:
-        """Fliter molecules by number of rotatable bonds
-
-        Fliter out molecules by the number of rotatable bonds.
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-        cutoff : int
-          The criterion of rotatable bond numbers for filtering. By default 40.
-
-        Returns
-        -------
-        Chem.rdchem.Mol or None
-          Return a Rdkit mol object if the input molecule meets the criterion, otherwise return None.
-        """
-        rotation_bond_num = rdMolDescriptors.CalcNumRotatableBonds(mol)
-        if rotation_bond_num < cutoff:
-            return mol
-        else:
-            return None
-
-    def remove_salt(self, mol: Chem.rdchem.Mol) -> Chem.rdchem.Mol or None:
-        """Romve salts or ions from molecules
-
-        Remove salts or ions from molecules
-
-        Parameters
-        ----------
-        mol : Chem.rdchem.Mol
-          The Rdkit mol object of the input molecule.
-
-        Returns
-        -------
-        Chem.rdchem.Mol or None
-          Return a Rdkit mol object if the input molecule meets the criterion, otherwise return None.
-        """
-        lfc = rdMolStandardize.LargestFragmentChooser()
-        try:
-            cleaned_molecule = rdMolStandardize.Cleanup(mol)
-            smiles = Chem.MolToSmiles(cleaned_molecule)
-            if "." in smiles:
-                cleaned_molecule = rdMolStandardize.Cleanup(
-                    lfc.choose(cleaned_molecule)
-                )
-        except:
-            return None
-        return cleaned_molecule
-
-    def print_valid_cpd(self, ori_df: pd.DataFrame) -> None:
-        """Print out the result of filtering
-
-        Parameters
-        ----------
-        ori_df : pd.DataFrame
-          The resulting dataframe after applying a filter.
-        """
-        print(f"Number of the original molecules: {ori_df.shape[0]}")
-        # print (f"{ori_df.shape[0] - filtered_df.shape[0]} invalid smiles found, dataset contains {filtered_df.shape[0]} valid molecules")
-        print(f"Number of molecules should be removed: {ori_df.Mol.isnull().sum()}")
-
-
-if __name__ == "__main__":
-    dataset_raw = pd.read_csv(
-        "./COCONUT_DB.smi", sep=" ", header=0, names=["SMILES", "ID"]
-    )
-    print(dataset_raw)
-
-    filter = Filter()
-    dataset_raw["Mol"] = dataset_raw.SMILES.apply(filter.read_as_mol)
-    filter.print_valid_cpd(dataset_raw)
-    dataset_vali = dataset_raw.loc[~dataset_raw["Mol"].isna()]
-
-    dataset_vali["Mol"] = dataset_vali.Mol.apply(filter)
-    filter.print_valid_cpd(dataset_vali)
-
-    dataset_vali = dataset_vali.loc[~dataset_vali["Mol"].isna()]
-    print(dataset_vali)
-
-    dataset_vali["Std_SMILES"] = dataset_vali.Mol.apply(Chem.MolToSmiles)
-    dataset_vali.to_csv(
-        "./COCONUT_after_filter.smi",
-        columns=["Std_SMILES", "ID"],
-        sep="\t",
-        index=False,
-        header=False,
-    )
diff --git a/paper.bib b/paper.bib
index f52acd7..1f72d43 100644
--- a/paper.bib
+++ b/paper.bib
@@ -93,13 +93,11 @@ @article{mysinger2012directory
 }
 
 @article{landrum2013rdkit,
-  title={RDKit: A software suite for cheminformatics, computational chemistry, and predictive modeling},
-  author={Landrum, Greg and others},
-  journal={Greg Landrum},
-  volume={8},
-  pages={31},
-  year={2013}
-}
+  title={RDKit: open-source cheminformatics},
+  author={Landrum, Greg},
+  year={2013},
+  note = {Available at \url{http://www.rdkit.org}, version 2023.3.3}
+  }
 
 @article{schreyer2012usrcat,
   title={USRCAT: real-time ultrafast shape recognition with pharmacophoric constraints},
@@ -132,3 +130,27 @@ @article{park2009simple
   publisher={Elsevier},
   doi = {10.1016/j.eswa.2008.01.039}
 }
+
+@article{gillet1998identification,
+  title={Identification of biological activity profiles using substructural analysis and genetic algorithms},
+  author={Gillet, Valerie J and Willett, Peter and Bradshaw, John},
+  journal={Journal of chemical information and computer sciences},
+  volume={38},
+  number={2},
+  pages={165--179},
+  year={1998},
+  publisher={ACS Publications},
+  doi={10.1021/ci970431+}
+}
+
+@article{wildman1999prediction,
+  title={Prediction of physicochemical parameters by atomic contributions},
+  author={Wildman, Scott A and Crippen, Gordon M},
+  journal={Journal of chemical information and computer sciences},
+  volume={39},
+  number={5},
+  pages={868--873},
+  year={1999},
+  publisher={ACS Publications},
+  doi={10.1021/ci990307l}
+}
\ No newline at end of file
diff --git a/paper.md b/paper.md
index b39a612..3958060 100644
--- a/paper.md
+++ b/paper.md
@@ -49,11 +49,11 @@ Whilst OpenFEPOPS has included functionality for descriptor caching and profilin
 The OpenFEPOPS descriptor generation process as outlined in \autoref{fig:descriptor_generation} follows;
 
 1. Tautomer enumeration
-    - For a given small molecule, OpenFEPOPS uses RDKit to iterate over molecular tautomers. By default, there is no limit to the number of recoverable tautomers but a limit may be imposed which may be necessary if adapting the OpenFEPOPS code to large macromolecules and not just small molecules.
+    - For a given small molecule, OpenFEPOPS uses RDKit [@landrum2013rdkit] to iterate over molecular tautomers. By default, there is no limit to the number of recoverable tautomers but a limit may be imposed which may be necessary if adapting the OpenFEPOPS code to large macromolecules and not just small molecules.
 2.  Conformer enumeration
     - For each tautomer, up to 1024 conformers are sampled by either complete enumeration of rotatable bond states (at the literature reported optimum increment of 90 degrees) if there are five or less rotatable bonds, or through random sampling of 1024 possible states if there are more than 5 rotatable bonds.
 3.  Defining feature points
-    - The KMeans algorithm [@arthur2007k] is applied to each conformer of each tautomer to identify four (by default) representative or central points, into which the atomic information of neighbouring atoms is collapsed. As standard, the atomic properties of charge, logP, hydrogen bond donor, and hydrogen bond acceptor status are collapsed into four feature points per unique tautomer conformation. These feature points are encoded to 22 numeric values (a FEPOP) comprising four points, each with four properties, and six pairwise distances between these points. With many FEPOPS descriptors collected from a single molecule through tautomer and conformer enumeration, this set of representative FEPOPS should capture every possible state of the original molecule.
+    - The KMeans algorithm [@arthur2007k] is applied to each conformer of each tautomer to identify four (by default) representative or central points, into which the atomic information of neighbouring atoms is collapsed. As standard, the atomic properties of charge, logP, hydrogen bond donor, and hydrogen bond acceptor status are collapsed into four feature points per unique tautomer conformation. The RDKit package is used to calculate these properties with the iterative Gasteiger charges algorithm [@wildman1999prediction] applied to assign atomic charges, the Crippen method [@wildman1999prediction] used to assign atomic logP contributions, and hydrogen bond acceptors and donors identified with appropriate [@gillet1998identification] SMARTS substructure queries. These feature points are encoded to 22 numeric values (a FEPOP) comprising four points, each with four properties, and six pairwise distances between these points. With many FEPOPS descriptors collected from a single molecule through tautomer and conformer enumeration, this set of representative FEPOPS should capture every possible state of the original molecule.
 4.  Selection of diverse FEPOPS
     - From the collection of FEPOPS derived from every tautomer conformation of a molecule, the K-Medoid algorithm [@park2009simple] is applied to identify seven (by default) diverse FEPOPS which are thought to best capture a fuzzy representation of the molecule. These seven FEPOPS each comprise 22 descriptors each, totalling 154 32-bit floating point numbers or 616 bytes.
 
@@ -74,15 +74,15 @@ Literature highlights that the choice of the Pearson correlation coefficient lea
 
 The predictive performance of OpenFEPOPS has been evaluated using the DUDE [@mysinger2012directory] diversity set. This dataset comprises eight protein targets accompanied by decoy ligands and known active ligands. Macro-averaged AUROC scores for each target were generated using  known actives to retrieve the entire set of actives for the target. Table 1 shows the average AUROC scores for DUDE diversity set targets along with scores obtained using the popular Morgan 2, MACCS, and RDKit fingerprints as implemented in RDKit and scored using the Tanimoto distance metric. See the Jupyter notebook 'Explore_DUDE_diversity_set.ipynb' in the source repository for further methods and data availability using the FigShare service. 
 
-|  | Morgan 2 | MACCS | RDKit |OpenFEPOPS|
+| Target | Morgan 2 | MACCS | RDKit |OpenFEPOPS|
 |--------:|----------:|-------:|-------:|------------:|
-|akt1  |0.836|0.741 |0.833 |0.829|
+|akt1  |0.836|0.741 |0.833 |0.831|
 |ampc  |0.784|0.673|0.660 |0.639|
-|cp3a4 |0.603|0.582|0.613 |0.650|
+|cp3a4 |0.603|0.582|0.613 |0.647|
 |cxcr4 |0.697|0.854|0.592 |0.899|
 |gcr   |0.670|0.666|0.708  |0.616|
 |hivpr |0.780|0.681|0.759 |0.678|
-|hivrt |0.651|0.670 |0.660 |0.584|
+|hivrt |0.651|0.670 |0.660 |0.582|
 |kif11 |0.763|0.668 |0.672  |0.713|
 
 **Table 1:** Macro averaged AUROC scores by target and molecular similarity technique for the DUDE diversity set. Across all datasets, 19 small molecules out of 112,796 were excluded from analysis mainly due to issues in parsing to valid structures using RDKit.
diff --git a/pyproject.toml b/pyproject.toml
index 160690d..53c4d31 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -9,10 +9,7 @@ authors = [
 readme = "README.md"
 requires-python = ">=3.9"
 keywords = ["molecular similarity", "ligand based virtual screening"]
-
-[build-system]
-requires=[
-    "setuptools>=61.0.0",
+dependencies = [
     "rdkit>=2021.9.5.1",
     "numpy>=1.22.4",
     "pandas>=1.5.0",
@@ -23,8 +20,15 @@ requires=[
     "fire>=0.5.0",
     "tqdm>=4.48.0",
     ]
+[build-system]
+requires=[
+    "setuptools>=61.0.0",
+    ]
 build-backend="setuptools.build_meta"
 
+
+
+
 [tool.black]
 skip-string-normalization = true
 include = '''
@@ -41,11 +45,6 @@ exclude = '''
 )
 '''
 
-dependencies = [
-    "requests",
-    'importlib-metadata; python_version<"3.10"',
-]
-
 [project.scripts]
 fepops = "fepops.__main__:fepops_entrypoint"
 
@@ -53,4 +52,14 @@ fepops = "fepops.__main__:fepops_entrypoint"
 where = ["src"]
 
 [project.urls]
-"Homepage" = "https://github.com/JustinYKC/FEPOPS"
\ No newline at end of file
+"Homepage" = "https://github.com/JustinYKC/FEPOPS"
+
+[project.optional-dependencies]
+dev = [
+    'black==23.3.0',
+    'pytest==7.4.0',
+    'build==0.10.0',
+]
+
+[tool.pytest.ini_options]
+"testpaths" = "tests/test*"
diff --git a/src/fepops/__init__.py b/src/fepops/__init__.py
index 4f23c13..fa78cb5 100644
--- a/src/fepops/__init__.py
+++ b/src/fepops/__init__.py
@@ -1,1 +1,3 @@
+"""FEPOPS module containing the OpenFEPOPS class"""
+
 from .fepops import OpenFEPOPS
diff --git a/src/fepops/__main__.py b/src/fepops/__main__.py
index a013ea2..f87960f 100644
--- a/src/fepops/__main__.py
+++ b/src/fepops/__main__.py
@@ -1,7 +1,8 @@
 import numpy as np
 import fire
+import json
 from pathlib import Path
-from typing import Union
+from typing import Optional
 from .fepops import OpenFEPOPS, GetFepopStatusCode
 from .fepops_persistent import get_persistent_fepops_storage_object
 
@@ -18,11 +19,11 @@ class FepopsCMDLineInterface:
 
     """
 
-    def __init__(self, database_file: str = None):
+    def __init__(self, database_file: str = None, json_file: Optional[str] = None):
         """Constructor for the command line inteface object
 
-        Allows singular definisition of a database file and subsequent use in
-        all subcommands.
+        Allows singular definition of a database file and subsequent use in
+        all subcommands and setting of JSON output file if required.
 
         Parameters
         ----------
@@ -30,8 +31,15 @@ def __init__(self, database_file: str = None):
             If a string, then this is used as the database/cache file used by
             all subprocesses. If None, then no database of cache file is used,
             by default None.
+        json_file : str, optional
+            If a string, then this is used as the JSON output file for results
+            from fepops generation and similarity calculations. If None, then
+            results are printed to the terminal, By default None
         """
         self.database_file = database_file
+        self.json_file = None if json_file is None else Path(json_file)
+        if self.json_file is not None and not self.json_file.parent.resolve().exists():
+            self.json_file.parent.resolve().mkdir(parents=True)
 
     def get_fepops(self, smi: str):
         """Get Fepops descriptors from a SMILES string
@@ -51,12 +59,27 @@ def get_fepops(self, smi: str):
             element is either None (if unsuccessful), or a np.ndarray containing
             the calculated Fepops descriptors of the requested input molecule.
         """
+
         if self.database_file is not None:
             with get_persistent_fepops_storage_object(self.database_file) as f:
-                return f.get_fepops(smi)
+                fepops_result = f.get_fepops(smi)
         else:
             f = OpenFEPOPS()
-            return f.get_fepops(smi)
+            fepops_result = f.get_fepops(smi)
+        if self.json_file is None:
+            return fepops_result
+        else:
+            fepops_status, fepops_features = fepops_result
+            with open(Path(self.json_file), 'w') as json_output:
+                json.dump(
+                    {
+                        "SMILES": smi,
+                        "FepopStatusCode": str(fepops_status.name),
+                        "Fepops": fepops_features.tolist(),
+                    },
+                    json_output,
+                    indent=4,
+                )
 
     def calc_sim(self, smi1: str, smi2: str):
         """Calculate FEPOPS similarity between two smiles strings
@@ -81,12 +104,32 @@ def calc_sim(self, smi1: str, smi2: str):
                     fepops_status1 is not GetFepopStatusCode.SUCCESS
                     or fepops_status2 is not GetFepopStatusCode.SUCCESS
                 ):
-                    return np.nan
+                    similarity_score = np.nan
                 else:
-                    return f.calc_similarity(fepops_features1, fepops_features2)
+                    similarity_score = f.calc_similarity(
+                        fepops_features1, fepops_features2
+                    )
         else:
             f = OpenFEPOPS()
-            print(f.calc_similarity(smi1, smi2))
+            similarity_score = f.calc_similarity(smi1, smi2)
+
+        if self.json_file is not None:
+            with open(Path(self.json_file), 'w') as json_output:
+                json.dump(
+                    {
+                        "SMI1": smi1,
+                        "SMI2": smi2,
+                        "FEPOPS Similarity Score": similarity_score
+                        if len(similarity_score) > 1
+                        else similarity_score[0],
+                    },
+                    json_output,
+                    indent=4,
+                )
+        else:
+            return (
+                similarity_score if len(similarity_score) > 1 else similarity_score[0]
+            )
 
     def save_descriptors(self, smi: str):
         """Pregenerate FEPOPS descriptors for a set of SMILES strings
@@ -110,7 +153,7 @@ def dude_preprocessor(self, dude_directory: str = "data/dude/"):
         decoys, bringing the total number of molecules to 1,434,022.
         More information is available here: https://dude.docking.org/
         """
-        from .dude_preprocessor import DudePreprocessor
+        from .utils import DudePreprocessor
 
         dude_preprop = DudePreprocessor(dude_directory=dude_directory)
         return dude_preprop
diff --git a/src/fepops/fepops.py b/src/fepops/fepops.py
index 767c729..12b5485 100644
--- a/src/fepops/fepops.py
+++ b/src/fepops/fepops.py
@@ -14,6 +14,7 @@
 from enum import Enum
 import multiprocessing as mp
 from multiprocessing import SimpleQueue
+import logging
 
 GetFepopStatusCode = Enum(
     "GetFepopStatusCode",
@@ -102,52 +103,52 @@ def __init__(
         num_fepops_per_mol: int = 7,
         num_centroids_per_fepop: int = 4,
         descriptor_means: Tuple[float, ...] = (
-            -0.28932319,
-            0.5166312,
-            0.37458883,
-            0.99913668,
-            -0.04193182,
-            1.03616917,
-            0.27327129,
-            0.99839024,
-            0.09701198,
-            1.12969387,
-            0.23718642,
-            0.99865705,
-            0.35968991,
-            0.6649304,
-            0.4123743,
-            0.99893657,
-            5.70852885,
-            6.3707943,
-            6.47354071,
-            6.26385429,
-            6.19229367,
-            6.22946713,
+            -0.28971602,
+            0.5181022,
+            0.37487135,
+            0.99922747,
+            -0.04187301,
+            1.03382471,
+            0.27407036,
+            0.99853436,
+            0.09725517,
+            1.12824307,
+            0.23735556,
+            0.99882914,
+            0.35977538,
+            0.66653514,
+            0.41238282,
+            0.99902545,
+            5.71261449,
+            6.37716992,
+            6.47293777,
+            6.26134733,
+            6.20354385,
+            6.23201498,
         ),
         descriptor_stds: Tuple[float, ...] = (
-            0.35067291,
-            1.00802116,
-            0.48380817,
-            0.02926675,
-            0.15400475,
-            0.86220776,
-            0.44542581,
-            0.03999429,
-            0.16085455,
-            0.92042695,
-            0.42515847,
-            0.03655217,
-            0.35778578,
-            1.36108994,
-            0.49210665,
-            0.03252466,
-            1.96446927,
-            2.30792259,
-            2.5024708,
-            2.4155645,
-            2.29434487,
-            2.31437527,
+            0.35110473,
+            1.00839329,
+            0.4838859,
+            0.02769204,
+            0.15418035,
+            0.86446056,
+            0.44583626,
+            0.0381767,
+            0.16095862,
+            0.92079483,
+            0.42526185,
+            0.03413741,
+            0.35756229,
+            1.36093993,
+            0.4921059,
+            0.0311619,
+            1.9668792,
+            2.31266486,
+            2.50699385,
+            2.41269982,
+            2.30018205,
+            2.31527129,
         ),
     ):
         """OpenFEPOPS (Feature Points) molecular similarity object
@@ -202,7 +203,7 @@ def __init__(
             0.09701198,1.12969387,0.23718642,0.99865705,0.35968991,0.6649304,
             0.4123743,0.99893657,5.70852885,6.3707943,6.47354071,6.26385429,
             6.19229367,6.22946713)
-        descriptor_sds : Tuple[float, ...], optional
+        descriptor_stds : Tuple[float, ...], optional
             Due to the need to apply scaling to FEPOPS, the DUDE diversity set has
             been profiled and the means collected for all contained FEPOPS. This
             this allows centering and scaling of FEPOPS before scoring. This field
@@ -219,9 +220,12 @@ def __init__(
         ValueError
             Invalid kmeans method
         """
-
-        self.descriptor_means = descriptor_means
-        self.descriptor_stds = descriptor_stds
+        # Descriptor stds may contain zeros. If they do, then we mimic Scikit-Learn's
+        # StandardScaler, whereby if unit variance is not achievable, no scaling is
+        # applied (value of 1.0)
+        self.descriptor_stds_no_zeros = np.array(descriptor_stds)
+        self.descriptor_stds_no_zeros[self.descriptor_stds_no_zeros == 0.0] = 1.0
+        self.descriptor_means = np.array(descriptor_means)
 
         try:
             self.kmeans_func = getattr(self, f"_perform_kmeans_{kmeans_method}")
@@ -282,31 +286,33 @@ def _get_k_medoids(
         if input_x.shape[0] <= k:
             return input_x
 
-        point_to_centroid_map = np.ones(input_x.shape[0])
+        # Apply standard scaling to FEPOP features. Behaviour when std dev is 0 mimics
+        # Scikit-Learn's StandardScaler, whereby if unit variance is not achievable, no
+        # scaling is applied (value of 1.0)
+        input_x_std = np.std(input_x, axis=0)
+        input_x_std[input_x_std == 0.0] = 1.0
+        X = (input_x - np.mean(input_x, axis=0)) / input_x_std
+        point_to_centroid_map = np.ones(X.shape[0])
         point_to_centroid_map_prev = np.zeros_like(point_to_centroid_map)
 
         np_rng = np.random.default_rng(seed=random_state)
-        medoids = input_x[
-            np_rng.choice(np.arange(input_x.shape[0]), size=k, replace=False), :
-        ]
+        medoids = X[np_rng.choice(np.arange(X.shape[0]), size=k, replace=False), :]
 
         while (point_to_centroid_map != point_to_centroid_map_prev).any():
             point_to_centroid_map_prev = point_to_centroid_map
-            point_to_centroid_map = np.argmin(
-                np.square(cdist(input_x, medoids)), axis=1
-            )
+            point_to_centroid_map = np.argmin(np.square(cdist(X, medoids)), axis=1)
             for i in range(k):
-                medoid_members = input_x[point_to_centroid_map == i]
+                medoid_members = X[point_to_centroid_map == i]
                 if len(medoid_members) == 0:
-                    chosen_x_point = np_rng.choice(np.arange(input_x.shape[0]))
-                    medoids[i] = input_x[chosen_x_point, :]
+                    chosen_x_point = np_rng.choice(np.arange(X.shape[0]))
+                    medoids[i] = X[chosen_x_point, :]
                     point_to_centroid_map[chosen_x_point] = i
-                medoids[i] = np.median(input_x[point_to_centroid_map == i], axis=0)
+                medoids[i] = np.median(X[point_to_centroid_map == i], axis=0)
         # Sorting at this stage for reproducibility with existing pregenerated
         # descriptor sets and convention with early FEPOPS versions which relied
         # upon FEPOPS within a molecule being sorted by charge (before moving to
         # the newer CombiAlign scoring algorithm)
-        return medoids[np.lexsort(medoids.T[::-1])]
+        return input_x[np.lexsort(medoids.T[::-1])]
 
     def _calculate_atomic_logPs(self, mol: Chem.rdchem.Mol) -> dict:
         """Calculate logP contribution for each of atom in a molecule
@@ -552,7 +558,6 @@ def _mol_from_smiles(self, smiles_string: str) -> Chem.rdchem.Mol:
             try:
                 mol = Chem.MolFromSmiles(smiles_string, sanitize=False)
             except:
-                f"Could not parse smiles to a valid molecule, smiles was: {smiles_string}"
                 return None
         return mol
 
@@ -673,7 +678,7 @@ def generate_conformers(self, mol: Chem.rdchem.Mol, random_state: int = 42) -> l
             params = AllChem.ETKDGv2()
             id = AllChem.EmbedMolecule(mol, params)
             if id == -1:
-                print(
+                logging.warning(
                     "Coords could not be generated without using random coords. using random coords now"
                 )
                 params.useRandomCoords = True
@@ -682,7 +687,7 @@ def generate_conformers(self, mol: Chem.rdchem.Mol, random_state: int = 42) -> l
                     AllChem.EmbedMolecule(mol, params)
                 )
             except ValueError:
-                print("Conformer embedding failed")
+                logging.warning("Conformer embedding failed")
                 return []
         dihedrals = self._get_dihedrals(mol)
         starting_angles = [
@@ -839,12 +844,12 @@ def get_fepops(
             original_smiles = mol
             mol = self._mol_from_smiles(mol)
         if mol is None:
-            print(
+            logging.error(
                 f"Failed to make a molecule{' from '+original_smiles if original_smiles is not None else ''}"
             )
             return GetFepopStatusCode.FAILED_TO_GENERATE, None
         if Lipinski.HeavyAtomCount(mol) < self.num_centroids_per_fepop:
-            print(
+            logging.error(
                 f"Number of heavy atoms ({Lipinski.HeavyAtomCount(mol)}) below requested feature points ({self.num_centroids_per_fepop}) for molecule {original_smiles if original_smiles is not None else ''}"
             )
             return GetFepopStatusCode.FAILED_TO_GENERATE, None
@@ -856,7 +861,7 @@ def get_fepops(
             conf_list = self.generate_conformers(t_mol)
             each_mol_with_all_confs_list.extend(conf_list)
         if each_mol_with_all_confs_list == []:
-            print(
+            logging.error(
                 f"Failed to generate conformers/tautomers {' for '+original_smiles if original_smiles is not None else ''}"
             )
             return GetFepopStatusCode.FAILED_TO_GENERATE, None
@@ -870,8 +875,8 @@ def get_fepops(
             )
         except ValueError as e:
             if original_smiles is not None:
-                print(f"Failed molecule had SMILES: {original_smiles}")
-            print(e)
+                logging.error(f"Failed molecule had SMILES: {original_smiles}")
+            logging.error(e)
             return GetFepopStatusCode.FAILED_TO_GENERATE, None
 
         medoids = self._get_k_medoids(
@@ -894,10 +899,10 @@ def pairwise_correlation(self, A: np.ndarray, B: np.ndarray):
         np.ndarray
             2D matrix containing A vs B feature correlations
         """
-        if len(A)<len(B):
-            A=np.pad(A, (0,len(B)-len(A)), mode='constant', constant_values=0)
-        if len(B)<len(A):
-            B=np.pad(B, (0,len(A)-len(B)), mode='constant', constant_values=0)
+        if len(A) < len(B):
+            A = np.pad(A, (0, len(B) - len(A)), mode='constant', constant_values=0)
+        if len(B) < len(A):
+            B = np.pad(B, (0, len(A) - len(B)), mode='constant', constant_values=0)
         am = A - np.mean(A, axis=0, keepdims=True)
         bm = B - np.mean(B, axis=0, keepdims=True)
         return (
@@ -961,12 +966,9 @@ def calc_similarity(
             raise ValueError(
                 "candidate was not, or could not be coerced into a np.ndarray"
             )
-        q = np.nan_to_num(
-            ((query - self.descriptor_means) / self.descriptor_stds), nan=1e-9
-        )
-        c = np.nan_to_num(
-            (candidate - self.descriptor_means) / self.descriptor_stds, nan=1e-9
-        )
+
+        q = (query - self.descriptor_means) / self.descriptor_stds_no_zeros
+        c = (candidate - self.descriptor_means) / self.descriptor_stds_no_zeros
         return self.pairwise_correlation(q.flatten(), c.flatten())
 
     def __call__(
diff --git a/src/fepops/fepops_persistent/__init__.py b/src/fepops/fepops_persistent/__init__.py
index cf51bd4..12be501 100755
--- a/src/fepops/fepops_persistent/__init__.py
+++ b/src/fepops/fepops_persistent/__init__.py
@@ -1,3 +1,4 @@
+"""fepops_persistent module contains functionality to cache/save FEPOPS descriptors to a file"""
 from .fepopsdb_json import FepopsDBJSON
 from .fepopsdb_sqlite import FepopsDBSqlite
 from .utils import get_persistent_fepops_storage_object
diff --git a/src/fepops/fepops_persistent/fepops_persistent_abc.py b/src/fepops/fepops_persistent/fepops_persistent_abc.py
index 5836025..988cbe8 100755
--- a/src/fepops/fepops_persistent/fepops_persistent_abc.py
+++ b/src/fepops/fepops_persistent/fepops_persistent_abc.py
@@ -9,6 +9,7 @@
 from fepops.fepops import GetFepopStatusCode
 from ..fepops import OpenFEPOPS
 import multiprocessing as mp
+import logging
 
 
 class FepopsPersistentAbstractBaseClass(metaclass=ABCMeta):
@@ -289,7 +290,9 @@ def _get_can_smi_mol_tuple(s: str, smiles_guaranteed_rdkit_canonical: bool = Fal
             except:
                 mol = None
         if mol is None:
-            print(f"Could not parse smiles to a valid molecule, smiles was: {s}")
+            logging.warning(
+                f"Could not parse smiles to a valid molecule, smiles was: {s}"
+            )
             return (s, mol)
         if smiles_guaranteed_rdkit_canonical:
             return (s, mol)
@@ -439,7 +442,9 @@ def calc_similarity(
                 return np.nan
         if any(x is None for x in (fepops_features_1, fepops_features_2)):
             return np.nan
-        score = self.openfepops_object.calc_similarity(fepops_features_1, fepops_features_2)
+        score = self.openfepops_object.calc_similarity(
+            fepops_features_1, fepops_features_2
+        )
         return score if score is not None else np.nan
 
     def write(self):
diff --git a/src/fepops/utils/__init__.py b/src/fepops/utils/__init__.py
new file mode 100644
index 0000000..eac32f0
--- /dev/null
+++ b/src/fepops/utils/__init__.py
@@ -0,0 +1,3 @@
+"""Utils to help with pre-filtering molecules and preparing benchmark files"""
+from .dude_preprocessor import DudePreprocessor
+from .mol_filter import Filter
diff --git a/src/fepops/dude_preprocessor.py b/src/fepops/utils/dude_preprocessor.py
similarity index 69%
rename from src/fepops/dude_preprocessor.py
rename to src/fepops/utils/dude_preprocessor.py
index 2d1acb9..742aa73 100644
--- a/src/fepops/dude_preprocessor.py
+++ b/src/fepops/utils/dude_preprocessor.py
@@ -5,6 +5,7 @@
 import pandas as pd
 import multiprocessing as mp
 from rdkit import Chem
+import logging
 
 
 class DudePreprocessor:
@@ -12,9 +13,21 @@ def __init__(
         self,
         *,
         dude_directory: Union[Path, str] = "data/dude/",
-    ) -> None:
+    ):
+        """Perform preprocessing on the DUDE database.
+
+        Assumes that the DUDE database has been downloaded and placed in the path
+        "data/dude/unprocessed". Here, we expect directories for each target containing
+        actives_final.ism and decoys_final.ism. These files are then processed into CSVs
+        representing the molecules (actives and decoys) for each DUDE target.
+
+        Parameters
+        ----------
+        dude_directory : Union[Path, str], optional
+            Directory containg , by default "data/dude/"
+
+        """
         print(dude_directory)
-        print(type(dude_directory))
         self.dude_path = Path(dude_directory)
         self.dude_unprocessed_path = self.dude_path / Path("unprocessed")
         if not self.dude_path.exists():
@@ -25,14 +38,38 @@ def __init__(
 
     def __call__(
         self,
+        targets: Optional[Union[str, list[str]]] = None,
+        skip_existing: bool = True,
     ):
-        self.process()
+        """Perform processing
+
+        Parameters
+        ----------
+        targets : Optional[Union[str, list[str]]], optional
+            Optionally provide a list of target names to process. If None, then
+            all targets are processed, by default None
+        skip_existing : bool, optional
+            If True, then existing processed target files are not regenerated, by
+            default True
+        """
+        self.process(targets=targets, skip_existing=skip_existing)
 
     def process(
         self,
-        targets: Optional[Union[str, list[str]]] = None,
+        targets: Optional[list[str]] = None,
         skip_existing: bool = True,
     ):
+        """Perform processing
+
+        Parameters
+        ----------
+        targets : Optional[Union[str, list[str]]], optional
+            Optionally provide a list of target names to process. If None, then
+            all targets are processed, by default None
+        skip_existing : bool, optional
+            If True, then existing processed target files are not regenerated, by
+            default True
+        """
         if targets is None:
             dude_targets = [
                 t.parent.name
@@ -43,7 +80,7 @@ def process(
                 targets = [targets]
         print(f"Processing the following DUDE targets: {dude_targets}")
         for target in tqdm(dude_targets, desc=f"Preparing targets"):
-            self.create_dude_target_csv_data(target, skip_existing=skip_existing)
+            self._create_dude_target_csv_data(target, skip_existing=skip_existing)
 
     @staticmethod
     def _parallel_init_worker_desc_gen_shared_fepops_ob():
@@ -58,7 +95,7 @@ def _parallel_get_rdkit_cansmi(s):
             return ""
         return Chem.MolToSmiles(mol)
 
-    def create_dude_target_csv_data(
+    def _create_dude_target_csv_data(
         self,
         dude_target: Path,
         actives_file: Path = Path("actives_final.ism"),
@@ -68,7 +105,7 @@ def create_dude_target_csv_data(
     ):
         target_output_file = self.dude_processed_path / f"dude_target_{dude_target}.csv"
         if skip_existing and target_output_file.exists():
-            print(
+            logging.warning(
                 f"Found existing {target_output_file}, skipping due to skip_existing = True, rerun as False to regenerate"
             )
         actives = pd.read_csv(
@@ -108,10 +145,10 @@ def cache_mols_from_csv(
             Path of CSV file. If None, then all CSV files in the DUDE datasets
             processed path are used.
         rdkit_canonical_smiles_column_header : str, optional
-            _description_, by default "rdkit_canonical_smiles"
+            Column header containing RDKit canonical SMILES, by default
+            "rdkit_canonical_smiles"
         """
-
-        from .fepops_persistent import get_persistent_fepops_storage_object
+        from ...fepops.fepops_persistent import get_persistent_fepops_storage_object
 
         for csv_path in (
             [Path(csv_path)]
@@ -120,7 +157,7 @@ def cache_mols_from_csv(
         ):
             df = pd.read_csv(csv_path)
             if df[rdkit_canonical_smiles_column_header].isnull().values.any():
-                print(
+                logging.warning(
                     f"Whilst working on caching {csv_path}, the following mol rows did not contain RDKit canonical SMILES:"
                 )
                 print(df[df[rdkit_canonical_smiles_column_header].isnull()])
diff --git a/src/fepops/utils/mol_filter.py b/src/fepops/utils/mol_filter.py
new file mode 100644
index 0000000..3e35306
--- /dev/null
+++ b/src/fepops/utils/mol_filter.py
@@ -0,0 +1,184 @@
+from rdkit import Chem
+from rdkit.Chem import rdMolDescriptors
+from rdkit.Chem.MolStandardize import rdMolStandardize
+from typing import Optional
+
+
+class Filter:
+    """Filter molecules by property
+
+    Filter object allows filtering by min atoms, max rings, and max rotatable
+    bonds present within a molecule.
+    """
+
+    def __init__(
+        self, min_atoms: int = 4, max_rings: int = 9, max_n_rot: int = 40
+    ) -> None:
+        self.min_atoms = min_atoms
+        self.max_rings = max_rings
+        self.max_n_rot = max_n_rot
+
+    def __call__(self, mol: Chem.rdchem.Mol) -> Optional[Chem.rdchem.Mol]:
+        """Apply all filters
+
+        Apply four filters, removing molecules which violate the defined minimum
+        cutoff for the number of atoms present and maximum cutoffs for rings and
+        rotatable bonds. Additionally remove salts/ions.
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          The Rdkit Mol object of the input molecule.
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Chem.rdchem.Mol if the molecule passes all filters, otherwise None.
+        """
+        return self.filter_mol(mol)
+
+    def filter_mol(self, mol: Chem.rdchem.Mol) -> Optional[Chem.rdchem.Mol]:
+        """Apply all filters
+
+        Apply four filters, removing molecules which violate the defined minimum
+        cutoff for the number of atoms present and maximum cutoffs for rings and
+        rotatable bonds. Additionally remove salts/ions.
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          The Rdkit Mol object of the input molecule.
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Chem.rdchem.Mol if the molecule passes all filters, otherwise None.
+        """
+        mol = self.filter_num_atoms(mol, self.min_atoms)
+        if mol is None:
+            return None
+        mol = self.filter_num_rings(mol, self.max_rings)
+        if mol is None:
+            return None
+        mol = self.rotation_bond_num(mol, self.max_n_rot)
+        if mol is None:
+            return None
+        mol = self.remove_salts(mol)
+        return mol
+
+    def smiles_to_mol(self, query_smiles: str) -> Chem.rdchem.Mol:
+        """Read SMILES strings from a file
+
+        Read molecular SMILES strings as Rdkit mol objects from a file.
+
+        Parameters
+        ----------
+        query_smiles : str
+          A molecular SMILES string.
+
+        Returns
+        -------
+        Chem.rdchem.Mol
+          Return a Rdkit mol object.
+        """
+        return Chem.MolFromSmiles(query_smiles)
+
+    def filter_num_atoms(
+        self, mol: Chem.rdchem.Mol, cutoff: int = 4
+    ) -> Optional[Chem.rdchem.Mol]:
+        """Filter out molecules with atom counts less than cutoff
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          Input molecule
+        cutoff : int
+          The number of atoms which the molecule must have to pass this filter
+          step. By default 4.
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Chem.rdchem.Mol if the input molecule meets the cutoff criteria,
+          otherwise None.
+        """
+        atom_num = mol.GetNumAtoms()
+        if atom_num > cutoff:
+            return mol
+        else:
+            return None
+
+    def filter_num_rings(
+        self, mol: Chem.rdchem.Mol, cutoff: int = 9
+    ) -> Optional[Chem.rdchem.Mol]:
+        """Filter out molecules with ring counts greater than cutoff
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          Input molecule
+        cutoff : int
+          The number of rings which the molecule must have less than to pass
+          this filter. By default 9
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Chem.rdchem.Mol if the input molecule meets the cutoff criteria,
+          otherwise None.
+        """
+        ring_num = mol.GetRingInfo().NumRings()
+        if ring_num < cutoff:
+            return mol
+        else:
+            return None
+
+    def rotation_bond_num(
+        self, mol: Chem.rdchem.Mol, cutoff: int = 40
+    ) -> Optional[Chem.rdchem.Mol]:
+        """Filter out molecules with rotatable bond counts greater than cutoff
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          Input molecule
+        cutoff : int
+          The number of rotatable bonds which the molecule must have less than
+          to pass this filter. By default 40
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Chem.rdchem.Mol if the input molecule meets the cutoff criteria,
+          otherwise None.
+        """
+        rotation_bond_num = rdMolDescriptors.CalcNumRotatableBonds(mol)
+        if rotation_bond_num < cutoff:
+            return mol
+        else:
+            return None
+
+    def remove_salts(self, mol: Chem.rdchem.Mol) -> Optional[Chem.rdchem.Mol]:
+        """Remove salts and ions from molecules
+
+        Parameters
+        ----------
+        mol : Chem.rdchem.Mol
+          Input molecule/salt
+
+        Returns
+        -------
+        Optional[Chem.rdchem.Mol]
+          Standardised molecule or None if an error in standardisation occurs
+        """
+        lfc = rdMolStandardize.LargestFragmentChooser()
+        try:
+            cleaned_molecule = rdMolStandardize.Cleanup(mol)
+            smiles = Chem.MolToSmiles(cleaned_molecule)
+            if "." in smiles:
+                cleaned_molecule = rdMolStandardize.Cleanup(
+                    lfc.choose(cleaned_molecule)
+                )
+        except:
+            return None
+        return cleaned_molecule


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/test/conftest.py b/tests/conftest.py
similarity index 100%
rename from test/conftest.py
rename to tests/conftest.py
diff --git a/test/tests.py b/tests/tests.py
similarity index 94%
rename from test/tests.py
rename to tests/tests.py
index 0188313..d689c9a 100644
--- a/test/tests.py
+++ b/tests/tests.py
@@ -1,6 +1,7 @@
 """Initial framework of FEPOPS test suite"""
 
 from pathlib import Path
+import numpy as np
 from fepops import OpenFEPOPS
 from fepops.fepops import GetFepopStatusCode
 from fepops.fepops_persistent import get_persistent_fepops_storage_object
@@ -60,3 +61,5 @@ def test_fepops_persistent_sqlite(test_mol_smiles, tmp_path):
     }
     file.unlink(missing_ok=True)
     assert fepops_dict.keys() == fepops_dict2.keys()
+    for k in fepops_dict.keys():
+        np.testing.assert_array_equal(fepops_dict[k], fepops_dict2[k])