Repository: JustinYKC/FEPOPS (language: python)

Issue:
Feature/benchmark data preprocessing
Add data pre-processing. 


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 7e2fd2d..d9a4021 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,169 @@
+# Initial part of this .gitignore is taken from
+# https://github.com/github/gitignore/blob/main/Python.gitignore
+# ---------------------------------------------------------------
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
-.vscode
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/#use-with-ide
+.pdm.toml
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
+
+# Fepops specific folders to ignore
+data/
+.vscode/
+/*.prof
diff --git a/benchmark.py b/benchmark.py
deleted file mode 100644
index f28c406..0000000
--- a/benchmark.py
+++ /dev/null
@@ -1,190 +0,0 @@
-import time
-import fire
-import numpy as np
-from sklearn.datasets import make_classification
-from tqdm import tqdm
-from fepops import Fepops
-from dataclasses import dataclass
-from typing import Callable, Union
-from pathlib import Path
-import pandas as pd
-from rdkit.Chem import AllChem
-from rdkit.Chem import DataStructs
-from sklearn.metrics import roc_auc_score
-from fepops.fepops_persistent import get_persistent_fepops_storage_object
-
-
-@dataclass
-class SimilarityMethod:
-    name: str
-    score: Callable
-
-
-class ROCScorer:
-    def __init__(self, methods: list[SimilarityMethod]) -> None:
-        self.similarity_methods = methods
-
-    def get_auroc_scores(
-        self,
-        df: pd.DataFrame,
-        smiles_column_title="SMILES",
-        active_flag_column_title="Active",
-    ) -> dict[str, float]:
-        if active_flag_column_title not in df.columns:
-            raise ValueError(
-                f"Could not find a '{active_flag_column_title}' column in the DataFrame to indicate if the compound is active, or a decoy. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
-            )
-        if smiles_column_title not in df.columns:
-            raise ValueError(
-                f"Could not find a '{smiles_column_title}' column in the DataFrame. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
-            )
-        labels = np.array(df[active_flag_column_title].tolist(), dtype=int)
-        scores_dict = {sm.name: [] for sm in self.similarity_methods}
-        for active in df.query(f"{active_flag_column_title}==1")[
-            smiles_column_title
-        ].tolist():
-            scores = np.array(
-                [
-                    [
-                        sim_method.score(active, s)
-                        for sim_method in self.similarity_methods
-                    ]
-                    for s in df[smiles_column_title].tolist()
-                ]
-            )
-            roc_scores = [roc_auc_score(labels, s) for s in scores.T]
-            [
-                scores_dict[m.name].append(rs)
-                for m, rs in zip(self.similarity_methods, roc_scores)
-            ]
-        return scores_dict
-
-
-class FepopsBenchmarker:
-    def __init__(self, database_file: str = "benchmark.db"):
-        """FEPOPS benchmarker
-
-        Contains test data useful in assessing the peformance of the FEPOPS object
-        """
-        self.fepops = get_persistent_fepops_storage_object(database_file=database_file)
-
-    def get_1k_x_1024_fepops(self):
-        n_useful_features = 22
-        return np.array(
-            [
-                make_classification(
-                    n_samples=1024,
-                    n_features=n_useful_features,
-                    n_informative=n_useful_features,
-                    n_redundant=0,
-                    n_repeated=0,
-                    n_classes=21,
-                    n_clusters_per_class=1,
-                    weights=None,
-                    flip_y=0,
-                    class_sep=1.0,
-                    hypercube=True,
-                    shift=0.0,
-                    scale=1.0,
-                    shuffle=True,
-                )[0]
-                for _ in tqdm(range(1000), "Generating test data")
-            ]
-        )
-
-    def kmedoids(self):
-        """Benchmark kmedoid code using the standard classification dataset
-
-        All benchmarks run using github codepace running Python 3.10.9 on a
-        2cpu instance (Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz)
-
-        Github commit 0e9f282, using the first kmeans implementation gives
-        a median time of 0.04102 seconds to extract 7 centroid fepops from
-        1024 fepops.
-
-        """
-
-        cached_1k_1024_fepops = self.get_1k_x_1024_fepops()
-        timings = np.empty(cached_1k_1024_fepops.shape[0])
-        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
-            start = time.time()
-            self._get_k_medoids(f_1024)
-            timings[i] = time.time() - start
-        print(
-            f"Median time to compute medoids with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
-        )
-
-    def kmeans(self):
-        """Benchmark kmedoid code using the standard classification dataset
-
-        All benchmarks run using github codepace running Python 3.10.9 on a
-        2cpu instance (Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz)
-
-        Github commit 0e9f282, using the first kmeans implementation gives
-        a median time of 0.04102 seconds to extract 7 centroid fepops from
-        1024 fepops.
-
-        """
-        fepops_sklearn = Fepops(kmeans_method="sklearn")
-
-        cached_1k_1024_fepops = self.get_1k_x_1024_fepops()
-        timings = np.empty(cached_1k_1024_fepops.shape[0])
-        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
-            start = time.time()
-            fepops_sklearn._perform_kmeans(f_1024, 7, kmeans_method="sklearn")
-            timings[i] = time.time() - start
-        print(
-            f"Median time to compute kmeans with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
-        )
-
-        fepops_ptcpu = Fepops(kmeans_method="pytorch-cpu")
-
-        timings = np.empty(cached_1k_1024_fepops.shape[0])
-        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
-            start = time.time()
-            fepops_ptcpu._perform_kmeans(f_1024, 7, kmeans_method="pytorch-cpu")
-            timings[i] = time.time() - start
-        print(
-            f"Median time to compute kmeans with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
-        )
-
-    @staticmethod
-    def _score_morgan(m1: str, m2: str, smi_to_mol_func) -> float:
-        m1 = smi_to_mol_func(m1)
-        m2 = smi_to_mol_func(m2)
-        fp1 = AllChem.GetMorganFingerprint(m1, 2)
-        fp2 = AllChem.GetMorganFingerprint(m2, 2)
-        return DataStructs.DiceSimilarity(fp1, fp2)
-
-    def auroc_performance(
-        self,
-        data_tsv: Union[Path, str],
-        smiles_column_title: str = "Std_SMILES",
-        active_flag_column_title: str = "Activity",
-    ):
-        roc_scorer = ROCScorer(
-            [
-                SimilarityMethod(
-                    "Morgan",
-                    lambda x, y: self._score_morgan(
-                        x, y, self.fepops.fepops_object._mol_from_smiles
-                    ),
-                ),
-                SimilarityMethod("FEPOPS", self.fepops.calc_similarity),
-            ]
-        )
-        scores_dict = roc_scorer.get_auroc_scores(
-            pd.read_csv(
-                Path(data_tsv),
-                sep="\t",
-                index_col=[0],
-                header=0,
-            ).reset_index(),
-            smiles_column_title="Std_SMILES",
-            active_flag_column_title="Activity",
-        )
-        print(pd.DataFrame.from_dict(scores_dict).describe())
-
-
-if __name__ == "__main__":
-    fire.Fire(FepopsBenchmarker)
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..d805f12
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,172 @@
+name: fepops
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - _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
+  - 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
+  - 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
+  - rdkit=2023.03.1=py311hc21844d_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=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:
+      - asttokens==2.2.1
+      - attrs==21.4.0
+      - backcall==0.2.0
+      - black==23.3.0
+      - certifi==2023.5.7
+      - 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
+      - easydict==1.10
+      - executing==1.2.0
+      - fast-pytorch-kmeans==0.2.0.1
+      - filelock==3.12.0
+      - fire==0.5.0
+      - h5py==3.8.0
+      - idna==3.4
+      - ipykernel==6.23.1
+      - ipython==8.13.2
+      - itsdangerous==2.1.2
+      - jedi==0.18.2
+      - jinja2==3.1.2
+      - jupyter-client==8.2.0
+      - jupyter-core==5.3.0
+      - lit==16.0.5
+      - markdown-it-py==2.2.0
+      - markupsafe==2.1.2
+      - matplotlib-inline==0.1.6
+      - 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
+      - prompt-toolkit==3.0.38
+      - psutil==5.9.5
+      - ptyprocess==0.7.0
+      - pure-eval==0.2.2
+      - pygments==2.15.1
+      - pynvml==11.5.0
+      - pyyaml==6.0
+      - pyzmq==25.1.0
+      - requests==2.31.0
+      - requests-cache==0.7.5
+      - rich==13.4.0
+      - scikit-learn==1.2.2
+      - scipy==1.10.1
+      - stack-data==0.6.2
+      - sympy==1.12
+      - termcolor==2.3.0
+      - threadpoolctl==3.1.0
+      - 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
+prefix: /home/sshave/miniconda3/envs/fepops
diff --git a/fepops/__main__.py b/fepops/__main__.py
deleted file mode 100644
index f8410ab..0000000
--- a/fepops/__main__.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import argparse
-from .fepops import Fepops
-from .fepops_persistent import get_persistent_fepops_storage_object
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="Generate the FEPOPS descriptors")
-    subparsers = parser.add_subparsers(
-        dest="subcmd", help="subcommands", metavar="SUBCOMMAND"
-    )
-    subparsers.required = True
-
-    f1_parser = subparsers.add_parser(
-        "get_fepops", help="Generate the FEPOPS descriptors from an input SMILES string"
-    )
-    f1_parser.add_argument(
-        "-ismi",
-        "--input_smiles",
-        help="An input SMILES string of a molecule",
-        dest="insmile",
-        required=True,
-    )
-    f1_parser.add_argument(
-        "-db",
-        "--database_file",
-        help="sqlite3 database file to use as a descriptor cache",
-        dest="database_file",
-        required=False,
-    )
-
-    f2_parser = subparsers.add_parser(
-        "calc_sim", help="Calculate FEPOPS similarity between two molecules"
-    )
-    f2_parser.add_argument(
-        "-ismi1",
-        "--input_smiles_1",
-        help="An input SMILES string of the first molecule",
-        dest="insmiles1",
-        required=True,
-    )
-    f2_parser.add_argument(
-        "-ismi2",
-        "--input_smiles_2",
-        help="An input SMILES string of the second molecule",
-        dest="insmiles2",
-        required=True,
-    )
-    f2_parser.add_argument(
-        "-db",
-        "--database_file",
-        help="sqlite3 database file to use as a descriptor cache",
-        dest="database_file",
-        required=False,
-    )
-
-    f3_parser = subparsers.add_parser(
-        "save_descriptors",
-        help="Pregenerate FEPOPS descriptors for a set of SMILES strings",
-    )
-    f3_parser.add_argument(
-        "-ismi",
-        "--input_smiles",
-        help="Input file containing SMILES strings",
-        dest="smiles_file",
-        required=True,
-    )
-    f3_parser.add_argument(
-        "-db",
-        "--database_file",
-        help="sqlite3 database file to use as a descriptor cache",
-        dest="database_file",
-        required=True,
-    )
-
-    args = parser.parse_args()
-    if args.subcmd == "get_fepops":
-        if args.database_file is not None:
-            f_persistent = get_persistent_fepops_storage_object(args.database_file)
-            fepops_features = f_persistent.get_fepop(args.insmile)
-        else:
-            f = Fepops()
-            fepops_features = f.get_fepops(args.insmile)
-        print(fepops_features)
-
-    if args.subcmd == "calc_sim":
-        if args.database_file is not None:
-            with get_persistent_fepops_storage_object(
-                args.database_file
-            ) as f_persistent:
-                fepops_features1 = f_persistent.get_fepop(args.insmiles1)
-                fepops_features2 = f_persistent.get_fepop(args.insmiles2)
-                print(f_persistent.calc_similarity(fepops_features1, fepops_features2))
-        else:
-            f = Fepops()
-            fepops_features_1 = f.get_fepops(args.insmiles1)
-            fepops_features_2 = f.get_fepops(args.insmiles2)
-            print(f.calc_similarity(fepops_features_1, fepops_features_2))
-
-    if args.subcmd == "save_descriptors":
-        f_persistent = get_persistent_fepops_storage_object(args.database_file)
-        f_persistent.save_descriptors(args.smiles_file)
-        if args.database_file.endswith((".json")):
-            f_persistent.write()
diff --git a/fepops/fepops_persistent/fepops_persistent_abc.py b/fepops/fepops_persistent/fepops_persistent_abc.py
deleted file mode 100755
index 18c0dfc..0000000
--- a/fepops/fepops_persistent/fepops_persistent_abc.py
+++ /dev/null
@@ -1,236 +0,0 @@
-from abc import ABCMeta, abstractmethod
-from pathlib import Path
-from typing import Union
-
-import numpy as np
-import torch
-from fast_pytorch_kmeans import KMeans as _FastPTKMeans
-from joblib import Parallel, delayed
-from rdkit import Chem
-from rdkit.Chem import AllChem, Crippen, MolStandardize, rdMolTransforms
-from scipy.spatial.distance import cdist, pdist, squareform
-from sklearn.cluster import KMeans as _SKLearnKMeans
-from sklearn.preprocessing import StandardScaler
-from tqdm import tqdm
-
-from ..fepops import Fepops
-
-
-class FepopsPersistentAbstractBaseClass(metaclass=ABCMeta):
-    """Abstract base class for persistent fepops storage
-
-    New storage methods may be implemented as demonstrated in fepopsdb_json.py
-    by extending this abstract base class which provides some required
-    functionality like:
-
-    - save_descriptors(smiles: Union[str, Path, list[str]]), to save a smiles file/list of smiles to the persistent storage
-    - get_cansmi_to_mol_dict_not_in_database(smiles: Union[str, Path, list[str]]), to retrieve a unique dictionary with canonical smiles as keys not already stored in the database and rdkit mol objects as values.
-
-    When writing your own persistent storage methods, you must override the following methods:
-
-    add_fepop(rdkit_canonical_smiles: str, fepops: np.ndarray)
-    ----------------------------------------------------------
-    Add the fepop to persistent storage. super().add_fepop may be called by the overridden function to perform type checks on arguments.
-
-    fepop_exists(rdkit_canonical_smiles: str)
-    -----------------------------------------
-    Return True if the canonical smiles is already in the database, and False if not. super().fepop_exists may be called by the overridden
-    function to perform type checks on arguments.
-
-    get_fepop(rdkit_canonical_smiles: str)
-    --------------------------------------
-    Return a fepop from persistent storage. If it does not exist, then generate it by calling self.fepops_object.get_fepop which is supplied
-    by this base class. super().get_fepop may be called by the overridden function to perform type checks on arguments.
-    With this function in place it allows interface compatibility with a standard Fepops object.
-
-    Inheriting functions may also define __enter__ and __exit__ methods for use with context handlers. If none are defined, then empty ones
-    are provided. This can be useful in doing things like writing out large files after descriptor generation if incremental writes are not
-    possible, like in the case of the FepopsDBJSON child class.
-
-
-    Parameters
-    ----------
-    database_file : Union[str, Path]
-            File to use for persistent storage.
-    kmeans_method : str, optional
-            Method which should be used for kmeans calculation by
-            fepops objects, can be one of "sklearn", "pytorch-gpu",
-            or "pytorch-cpu".
-    parallel : bool, optional
-            Run in parallel (using joblib), by default True
-    n_jobs : int, optional
-            Number of jobs to be spawned with joblib. If -1, then use
-            all available cores. By default -1.
-    """
-
-    @abstractmethod
-    def __init__(
-        self,
-        database_file: Union[str, Path],
-        kmeans_method: str = "pytorch-cpu",
-        parallel: bool = True,
-        n_jobs: int = -1,
-    ):
-        self.database_file = Path(database_file)
-        self.fepops_object = Fepops(kmeans_method=kmeans_method)
-        self.parallel = parallel
-        self.n_jobs = n_jobs
-
-    def save_descriptors(
-        self,
-        smiles: Union[str, Path, list[str]],
-    ):
-        canonical_smiles_to_mol_dict = self.get_cansmi_to_mol_dict_not_in_database(
-            smiles
-        )
-        if not self.parallel:
-            for rdkit_canonical_smiles, mol in canonical_smiles_to_mol_dict.items():
-                self.add_fepop(
-                    rdkit_canonical_smiles, self.fepops_object.get_fepops(mol)
-                )
-            print(
-                f"Added {len(canonical_smiles_to_mol_dict)} new molecues to the database ({self.database_file})"
-            )
-        else:  # Do it in parallel
-            cansmi_fepops_tuples = Parallel(n_jobs=self.n_jobs, prefer="threads")(
-                delayed(lambda sm: (sm[0], self.fepops_object.get_fepops(sm[1])))(
-                    (cs, m)
-                )
-                for cs, m in canonical_smiles_to_mol_dict.items()
-            )
-            for rdkit_canonical_smiles, new_fepop in cansmi_fepops_tuples:
-                self.add_fepop(rdkit_canonical_smiles, new_fepop)
-            print(
-                f"Added {len(canonical_smiles_to_mol_dict)} new molecues to the database ({self.database_file})"
-            )
-
-    @abstractmethod
-    def add_fepop(self, rdkit_canonical_smiles: str, fepops: np.ndarray):
-        """Add canonical smiles and fepop to database.
-
-        Must be overridden
-        """
-        if not isinstance(rdkit_canonical_smiles, str):
-            raise ValueError(
-                f"Expected an rdkit canonical smiles string, but a {type(rdkit_canonical_smiles)} was passed"
-            )
-        if not isinstance(fepops, np.ndarray):
-            raise ValueError(f"Expected a fepop, but a {type(fepops)} was passed")
-
-    @abstractmethod
-    def get_fepop(self, smiles: str) -> Union[np.ndarray, None]:
-        if not isinstance(smiles, str):
-            raise ValueError(
-                f"Expected an rdkit canonical smiles string, but a {type(smiles)} was passed"
-            )
-
-    @abstractmethod
-    def fepop_exists(self, rdkit_canonical_smiles: str) -> bool:
-        """Return True if canonical smiles already exist in the database"""
-        if not isinstance(rdkit_canonical_smiles, str):
-            raise ValueError(
-                f"Expected an rdkit canonical smiles string, but a {type(rdkit_canonical_smiles)} was passed"
-            )
-
-    @staticmethod
-    def _get_can_smi_mol_tuple(s: str, is_canonical: bool = False):
-        try:
-            mol = Chem.MolFromSmiles(s)
-        except:
-            try:
-                mol = Chem.MolFromSmiles(s, sanitize=False)
-            finally:
-                mol = None
-        if mol is None:
-            print(f"Could not parse smiles to a valid molecule, smiles was: {s}")
-            return (s, mol)
-        if is_canonical:
-            return (s, mol)
-        else:
-            return (Chem.CanonSmiles(Chem.MolToSmiles(mol)), mol)
-
-    def get_cansmi_to_mol_dict_not_in_database(
-        self, smiles: Union[str, Path, list[str]]
-    ):
-        if isinstance(smiles, str):
-            smiles = Path(smiles)
-        if isinstance(smiles, Path):
-            if smiles.exists():
-                smiles = [
-                    s.strip() for s in open(smiles).readlines() if len(s.strip()) > 0
-                ]
-            else:
-                raise ValueError(
-                    f"smiles file ({smiles}) not found. If you are passing smiles, make it into a list"
-                )
-        if not isinstance(smiles, list):
-            raise ValueError(
-                "smiles should be a str or Path denoting the location of a smiles file, or a list of smiles"
-            )
-        smiles = list(set(smiles))
-        if not self.parallel:
-            # Ensure unique (canonical, also storing intermediate mol)
-            canonical_smiles_to_mol_dict = dict(
-                self._get_can_smi_mol_tuple(s)
-                for s in tqdm(smiles, desc="Ensuring unique")
-            )
-        else:
-            # Ensure unique (canonical, also storing intermediate mol)
-            canonical_smiles_to_mol_dict = dict(
-                Parallel(n_jobs=self.n_jobs)(
-                    delayed(FepopsPersistentAbstractBaseClass._get_can_smi_mol_tuple)(s)
-                    for s in smiles
-                )
-            )
-        # Make sure none are already in the database
-        canonical_smiles_to_mol_dict = {
-            cansmi: mol
-            for cansmi, mol in canonical_smiles_to_mol_dict.items()
-            if not self.fepop_exists(cansmi)
-        }
-        return canonical_smiles_to_mol_dict
-
-    def calc_similarity(
-        self,
-        fepops_features_1: Union[np.ndarray, str, None],
-        fepops_features_2: Union[np.ndarray, str, None],
-    ):
-        """Calculate FEPOPS similarity
-
-        A static method for calculating molecular similarity based on their FEPOPS descriptors.
-
-        Parameters
-        ----------
-        fepops_features_1 : Union[np.ndarray, str]
-                A Numpy array containing the FEPOPS descriptors of the query molecule
-                or a smiles string from which to generate FEPOPS descriptors for the
-                query molecule.
-        fepops_features_2 : Union[np.ndarray, str]
-                A Numpy array containing the FEPOPS descriptors of the candidate
-                molecule or a smiles string from which to generate FEPOPS descriptors
-                for the query molecule.
-
-        Returns
-        -------
-        float
-                Fepops similarity between two molecules
-        """
-
-        if isinstance(fepops_features_1, str):
-            fepops_features_1 = self.get_fepop(fepops_features_1)
-        if isinstance(fepops_features_2, str):
-            fepops_features_2 = self.get_fepop(fepops_features_2)
-        if any(x is None for x in (fepops_features_1, fepops_features_2)):
-            raise ValueError(
-                f"Unable to calculate similarity due to NoneType found in the fepops features:(fepops_features_1, fepops_features_2)=({type(fepops_features_1)}, {type(fepops_features_2)})"
-            )
-        return self.fepops_object.calc_similarity(fepops_features_1, fepops_features_2)
-
-    def write(self):
-        pass
-
-    def __enter__(self):
-        return self
-
-    def __exit__(self, exc_type, exc_value, traceback):
-        self.write()
diff --git a/fepops/fepops_persistent/utils.py b/fepops/fepops_persistent/utils.py
deleted file mode 100644
index 29061aa..0000000
--- a/fepops/fepops_persistent/utils.py
+++ /dev/null
@@ -1,8 +0,0 @@
-from . import FepopsDBSqlite, FepopsDBJSON
-
-
-def get_persistent_fepops_storage_object(database_file: str):
-    if database_file.endswith((".sqlite", ".sqlite3", ".db", ".db3", ".s3db", ".sl3")):
-        return FepopsDBSqlite(database_file)
-    if database_file.endswith((".json")):
-        return FepopsDBJSON(database_file)
diff --git a/pyproject.toml b/pyproject.toml
index 68b1f49..f2a7d5c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,53 @@
+[project]
+name="Fepops"
+version="1.4.0"
+description="Python implementation of the FEPOPS molecular descriptors"
+authors = [
+    {name = "Yan-Kai Chen"},
+    {name = "Steven Shave"},
+]
+readme = "README.rst"
+requires-python = ">=3.7"
+keywords = ["molecular similarity", "ligand based virtual screening"]
+
 [build-system]
-requires = [
-    "setuptools>=42",
-    "wheel",
+requires=[
+    "setuptools>=61.0.0",
     "torch",
+    "numpy>=1.19.0",
+    "pandas>=1.5.0",
+    "rdkit>=2019.09.0.0",
+    "scikit-learn>=0.20.4",
+    "scipy>=1.7.0",
+    "fast_pytorch_kmeans>=0.1.9",
+    "tqdm>=4.48.0",
+    "rdkit",
+    ]
+build-backend="setuptools.build_meta"
+
+[tool.black]
+skip-string-normalization = true
+include = '''
+(
+    ^/tests/
+  | ^/src/
+  | ^/setup[.]py
+)
+'''
+exclude = '''
+(
+    __pycache__
+    |.*\.egg-info
+)
+'''
+
+dependencies = [
+    "requests",
+    'importlib-metadata; python_version<"3.8"',
 ]
-build-backend = 'setuptools.build_meta'
+
+[project.scripts]
+fepops = "fepops.__main__:fepops_entrypoint"
+
+[tool.setuptools.packages.find]
+where = ["src"]
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index 15b33ee..0000000
--- a/setup.cfg
+++ /dev/null
@@ -1,32 +0,0 @@
-[metadata]
-name = fepops
-version = 1.3.0
-author = Yan-Kai Chen
-description = Python implementation of the FEPOPS molecular descriptors
-long_description = file: README.md
-long_description_content_type = text/markdown
-url = https://github.com/JustinYKC/FEPOPS
-classifiers =
-    Programming Language :: Python :: 3
-    License :: OSI Approved :: MIT License
-    Operating System :: OS Independent
-
-[options]
-package_dir =
-    = .
-packages = 
-    fepops
-
-python_requires = >=3.9
-install_requires =
-    numpy>=1.19.0
-    pandas>=1.5.0
-    rdkit>=2019.09.0.0
-    scikit-learn>=0.20.4
-    scipy>=1.7.0
-    Pytorch>=1.0.0
-    fast_pytorch_kmeans>=0.1.9
-    tqdm>=4.48.0
-
-[tool:pytest]
-testpaths = test
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..6068493
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,3 @@
+from setuptools import setup
+
+setup()
diff --git a/fepops/__init__.py b/src/fepops/__init__.py
similarity index 100%
rename from fepops/__init__.py
rename to src/fepops/__init__.py
diff --git a/src/fepops/__main__.py b/src/fepops/__main__.py
new file mode 100644
index 0000000..09fef2a
--- /dev/null
+++ b/src/fepops/__main__.py
@@ -0,0 +1,93 @@
+import numpy as np
+import fire
+from .fepops import Fepops, GetFepopStatusCode
+from .fepops_persistent import get_persistent_fepops_storage_object
+
+
+class FepopsCMDLineInterface:
+    def __init__(self, database_file: str = None):
+        self.database_file = database_file
+
+    def get_fepops(self, smi: str):
+        """Get Fepops descriptors from a SMILES string
+
+        Parameters
+        ----------
+        smi : str
+            Molecule as a SMILES string. Can also be None, in which case a failure error
+            status is returned along with None in place of the requested Fepops
+            descriptors.
+
+        Returns
+        -------
+        Tuple[GetFepopStatusCode, Union[np.ndarray, None]]
+            Returns a tuple, with the first value being a GetFepopStatusCode
+            (enum) denoting SUCCESS or FAILED_TO_GENERATE. The second tuple
+            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)
+        else:
+            f = Fepops()
+            return f.get_fepops(smi)
+
+    def calc_sim(self, smi1: str, smi2: str):
+        """Calculate FEPOPS similarity between two smiles strings
+
+        Parameters
+        ----------
+        smi1 : str
+            Smiles string representing molecule 1
+        smi2 : str
+            Smiles string representing molecule 2
+
+        Returns
+        -------
+        float
+            Fepops similarity between the two supplied molecules
+        """
+        if self.database_file is not None:
+            with get_persistent_fepops_storage_object(self.database_file) as f:
+                fepops_status1, fepops_features1 = f.get_fepops(smi1)
+                fepops_status2, fepops_features2 = f.get_fepops(smi2)
+                if (
+                    fepops_status1 is not GetFepopStatusCode.SUCCESS
+                    or fepops_status2 is not GetFepopStatusCode.SUCCESS
+                ):
+                    return np.nan
+                else:
+                    return f.calc_similarity(fepops_features1, fepops_features2)
+        else:
+            f = Fepops()
+            print(f.calc_similarity(smi1, smi2))
+
+    def save_descriptors(self, smi):
+        """Pregenerate FEPOPS descriptors for a set of SMILES strings"""
+        f_persistent = get_persistent_fepops_storage_object(self.database)
+        f_persistent.save_descriptors(smi)
+        if self.database_file.endswith((".json")):
+            f_persistent.write()
+
+    def dude_preprocessor(self, dude_directory: str = "data/dude/"):
+        """Access functions for preparation and benchmarking of the DUDE dataset
+
+        The DUDE dataset consists of 102 targets, 22,886 active compounds along with
+        decoys. More information is available here:
+        https://dude.docking.org/
+
+        """
+        from .dude_preprocessor import DudePreprocessor
+
+        print("Dude DIR========", dude_directory)
+        dude_preprop = DudePreprocessor(dude_directory=dude_directory)
+        return dude_preprop
+
+
+def fepops_entrypoint():
+    fire.Fire(FepopsCMDLineInterface)
+
+
+if __name__ == '__main__':
+    fire.Fire(FepopsCMDLineInterface)
diff --git a/src/fepops/benchmark.py b/src/fepops/benchmark.py
new file mode 100644
index 0000000..977d93a
--- /dev/null
+++ b/src/fepops/benchmark.py
@@ -0,0 +1,440 @@
+import time
+import numpy as np
+from sklearn.datasets import make_classification
+from tqdm import tqdm
+from fepops import Fepops
+from dataclasses import dataclass
+from typing import Callable, Union, Optional
+from pathlib import Path
+import pandas as pd
+from rdkit import Chem
+from rdkit.Chem import AllChem
+from rdkit.Chem import DataStructs
+from rdkit.Chem.MolStandardize import rdMolStandardize
+from sklearn.metrics import roc_auc_score
+from fepops.fepops_persistent import get_persistent_fepops_storage_object
+from typing import Union, Optional
+
+
+@dataclass
+class SimilarityMethod:
+    name: str
+    supports_multiple_candidates: bool
+    score: Callable
+    descriptor_calc_func: Optional[Callable] = None
+    descriptor_score_func: Optional[Callable] = None
+
+
+class Filter:
+    def __init__(self, min_atoms: int = 4) -> None:
+        self.min_atoms = min_atoms
+        self.lfc = rdMolStandardize.LargestFragmentChooser()
+
+    def _remove_salt(self, mol: Chem.rdchem.Mol) -> Union[Chem.rdchem.Mol, 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.
+        """
+        try:
+            cleaned_molecule = rdMolStandardize.Cleanup(mol)
+            smiles = Chem.MolToSmiles(cleaned_molecule)
+            if "." in smiles:
+                cleaned_molecule = rdMolStandardize.Cleanup(
+                    self.lfc.choose(cleaned_molecule)
+                )
+        except:
+            return None
+        return cleaned_molecule
+
+    def _filter_by_atom_num(self, mol: Chem.rdchem.Mol) -> Union[Chem.rdchem.Mol, 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 > self.min_atoms:
+            return mol
+        else:
+            return None
+
+    def __call__(self, mol: Chem.rdchem.Mol) -> Union[Chem.rdchem.Mol, 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) -> Union[Chem.rdchem.Mol, None]:
+        """Apply all filters
+
+        Perform two filters (minimum number of atoms, 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._remove_salt(mol)
+        if mol is None:
+            return None
+        else:
+            mol = self._filter_by_atom_num(mol)
+            if mol is None:
+                return None
+            else:
+                return mol
+
+
+class ROCScorer:
+    def __init__(self, methods: list[SimilarityMethod]) -> None:
+        self.similarity_methods = methods
+
+    def get_auroc_scores(
+        self,
+        df: pd.DataFrame,
+        smiles_column_title="SMILES",
+        active_flag_column_title="Active",
+    ) -> dict[str, float]:
+        if active_flag_column_title not in df.columns:
+            raise ValueError(
+                f"Could not find a '{active_flag_column_title}' column in the DataFrame to indicate if the compound is active, or a decoy. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
+            )
+        if smiles_column_title not in df.columns:
+            raise ValueError(
+                f"Could not find a '{smiles_column_title}' column in the DataFrame. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
+            )
+
+        smiles_list = df[smiles_column_title].tolist()
+        labels_list = df[active_flag_column_title].tolist()
+        descriptors = []
+        for sm in self.similarity_methods:
+            descriptors.append([])
+            for smiles in tqdm(smiles_list, desc=f"Caching descriptors for {sm.name}"):
+                res = sm.descriptor_calc_func(smiles)
+                if isinstance(res, tuple):
+                    descriptors[-1].append(res[-1])
+                else:
+                    descriptors[-1].append(res)
+
+        problematic_molecule_ilocs = []
+        assert (
+            len(set([len(d) for d in descriptors])) == 1
+        ), "Generation of descriptors yielded different lengths, something went very wrong"
+        for j in range(len(descriptors[0])):
+            vals = [descriptors[i][j] is not None for i in range(len(descriptors))]
+            if not all(vals):
+                problematic_molecule_ilocs.append(j)
+        problematic_molecule_smiles = df.iloc[problematic_molecule_ilocs][
+            smiles_column_title
+        ].tolist()
+        if len(problematic_molecule_smiles) > 0:
+            print('─' * 80)
+            print("Warning")
+            print(
+                "Problematic molecules for which no descriptors could be retrieved were found. "
+            )
+            print("The following molecules will not be included in AUROC calculations:")
+            for pm in problematic_molecule_smiles:
+                print(pm)
+            print('─' * 80)
+            problematic_molecule_ilocs.sort(reverse=True)
+            for pmi in problematic_molecule_ilocs:
+                del labels_list[pmi]
+                for d_i in range(len(descriptors)):
+                    del descriptors[d_i][pmi]
+
+        labels = np.array(labels_list, dtype=int)
+        scores_dict = {sm.name: [] for sm in self.similarity_methods}
+        for sm_i, sm in enumerate(self.similarity_methods):
+            for active_i in tqdm(
+                np.argwhere(np.array(labels_list) == 1).flatten(),
+                desc=f"Assessing active recall (AUROC) for {sm.name}",
+            ):
+                if sm.supports_multiple_candidates:
+                    scores = np.array(
+                        sm.descriptor_score_func(
+                            descriptors[sm_i][active_i], descriptors[sm_i]
+                        ),
+                        dtype=float,
+                    )
+                else:
+                    scores = np.array(
+                        [
+                            sm.descriptor_score_func(
+                                descriptors[sm_i][active_i], descriptors[sm_i][smiles_i]
+                            )
+                            for smiles_i in range(len(descriptors[0]))
+                        ],
+                        dtype=float,
+                    )
+                scores_dict[sm.name].append(
+                    roc_auc_score(
+                        labels[np.argwhere(~np.isnan(scores))],
+                        scores[np.argwhere(~np.isnan(scores))],
+                    )
+                )
+        return scores_dict
+
+    def get_auroc_scores_fast(
+        self,
+        df: pd.DataFrame,
+        smiles_column_title="SMILES",
+        active_flag_column_title="Active",
+    ) -> dict[str, float]:
+        if active_flag_column_title not in df.columns:
+            raise ValueError(
+                f"Could not find a '{active_flag_column_title}' column in the DataFrame to indicate if the compound is active, or a decoy. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
+            )
+        if smiles_column_title not in df.columns:
+            raise ValueError(
+                f"Could not find a '{smiles_column_title}' column in the DataFrame. Consider passing a column title to the 'active_flag_column_title' argument to get_auroc_scores"
+            )
+        labels = np.array(df[active_flag_column_title].tolist(), dtype=int)
+        descriptors = {
+            sm.name: {m: sm.descriptor_calc_func(m) for m in df[smiles_column_title]}
+            for sm in self.similarity_methods
+        }
+        scores_dict = {sm.name: [] for sm in self.similarity_methods}
+        for active in tqdm(
+            df.query(f"{active_flag_column_title}==1")[smiles_column_title].tolist(),
+            desc="Assessing active recall (AUROC)",
+        ):
+            scores = np.array(
+                [
+                    [
+                        sim_method.score(
+                            descriptors[sim_method.name][active],
+                            descriptors[sim_method.name][s],
+                        )
+                        for sim_method in self.similarity_methods
+                    ]
+                    for s in df[smiles_column_title].tolist()
+                ]
+            )
+            roc_scores = [
+                roc_auc_score(
+                    labels[np.argwhere(~np.isnan(s))], s[np.argwhere(~np.isnan(s))]
+                )
+                for s in scores
+            ]
+            [
+                scores_dict[m.name].append(rs)
+                for m, rs in zip(self.similarity_methods, roc_scores)
+            ]
+        return scores_dict
+
+
+class FepopsBenchmarker:
+    def __init__(self, database_file: str = "benchmark.db"):
+        """FEPOPS benchmarker
+
+        Contains test data useful in assessing the peformance of the FEPOPS object
+        """
+        self.fepops = get_persistent_fepops_storage_object(database_file=database_file)
+
+    def get_1k_x_1024_fepops(self):
+        n_useful_features = 22
+        return np.array(
+            [
+                make_classification(
+                    n_samples=1024,
+                    n_features=n_useful_features,
+                    n_informative=n_useful_features,
+                    n_redundant=0,
+                    n_repeated=0,
+                    n_classes=21,
+                    n_clusters_per_class=1,
+                    weights=None,
+                    flip_y=0,
+                    class_sep=1.0,
+                    hypercube=True,
+                    shift=0.0,
+                    scale=1.0,
+                    shuffle=True,
+                )[0]
+                for _ in tqdm(range(1000), "Generating test data")
+            ]
+        )
+
+    def kmedoids(self):
+        """Benchmark kmedoid code using the standard classification dataset
+
+        All benchmarks run using github codepace running Python 3.10.9 on a
+        2cpu instance (Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz)
+
+        Github commit 0e9f282, using the first kmeans implementation gives
+        a median time of 0.04102 seconds to extract 7 centroid fepops from
+        1024 fepops.
+
+        """
+
+        cached_1k_1024_fepops = self.get_1k_x_1024_fepops()
+        timings = np.empty(cached_1k_1024_fepops.shape[0])
+        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
+            start = time.time()
+            self._get_k_medoids(f_1024)
+            timings[i] = time.time() - start
+        print(
+            f"Median time to compute medoids with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
+        )
+
+    def kmeans(self):
+        """Benchmark kmedoid code using the standard classification dataset
+
+        All benchmarks run using github codepace running Python 3.10.9 on a
+        2cpu instance (Intel(R) Xeon(R) Platinum 8272CL CPU @ 2.60GHz)
+
+        Github commit 0e9f282, using the first kmeans implementation gives
+        a median time of 0.04102 seconds to extract 7 centroid fepops from
+        1024 fepops.
+
+        """
+        fepops_sklearn = Fepops(kmeans_method="sklearn")
+
+        cached_1k_1024_fepops = self.get_1k_x_1024_fepops()
+        timings = np.empty(cached_1k_1024_fepops.shape[0])
+        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
+            start = time.time()
+            fepops_sklearn._perform_kmeans(f_1024, 7, kmeans_method="sklearn")
+            timings[i] = time.time() - start
+        print(
+            f"Median time to compute kmeans with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
+        )
+
+        fepops_ptcpu = Fepops(kmeans_method="pytorch-cpu")
+
+        timings = np.empty(cached_1k_1024_fepops.shape[0])
+        for i, f_1024 in enumerate(tqdm(cached_1k_1024_fepops, "Benchmarking")):
+            start = time.time()
+            fepops_ptcpu._perform_kmeans(f_1024, 7, kmeans_method="pytorch-cpu")
+            timings[i] = time.time() - start
+        print(
+            f"Median time to compute kmeans with sklearn from 1024 fepops={np.median(timings)}, mean={np.mean(timings)}"
+        )
+
+    @staticmethod
+    def _score_morgan(m1: str, m2: Union[str, None], smi_to_mol_func) -> float:
+        m1 = smi_to_mol_func(m1)
+        if m2 is None:
+            return AllChem.GetMorganFingerprint(m1, 2)
+        m2 = smi_to_mol_func(m2)
+        fp1 = AllChem.GetMorganFingerprint(m1, 2)
+        fp2 = AllChem.GetMorganFingerprint(m2, 2)
+        return DataStructs.DiceSimilarity(fp1, fp2)
+
+    def auroc_performance(
+        self,
+        data_tsv: Union[Path, str],
+        data_tsv_contains_canonical_smiles: bool = True,
+        smiles_column_title: str = "SMILES",
+        active_flag_column_title: str = "Active",
+        results_output_dir: Optional[Union[str, Path]] = None,
+    ):
+        roc_scorer = ROCScorer(
+            [
+                SimilarityMethod(
+                    "Morgan",
+                    False,
+                    lambda x, y: self._score_morgan(
+                        x, y, self.fepops.fepops_object._mol_from_smiles
+                    ),
+                    lambda x: self._score_morgan(
+                        x, None, self.fepops.fepops_object._mol_from_smiles
+                    ),
+                    lambda x, y: DataStructs.DiceSimilarity(x, y),
+                ),
+                SimilarityMethod(
+                    "FEPOPS",
+                    True,
+                    self.fepops.calc_similarity,
+                    lambda x: self.fepops.get_fepops(
+                        x, is_canonical=data_tsv_contains_canonical_smiles
+                    ),
+                    self.fepops.calc_similarity,
+                ),
+            ]
+        )
+        scores_dict = roc_scorer.get_auroc_scores(
+            pd.read_csv(
+                Path(data_tsv),
+                sep="\t",
+                index_col=[0],
+                header=0,
+            ).reset_index(),
+            smiles_column_title=smiles_column_title,
+            active_flag_column_title=active_flag_column_title,
+        )
+        auroc_results_df = pd.DataFrame.from_dict(scores_dict)
+        if results_output_dir is None:
+            results_output_dir = Path(data_tsv).parent
+        results_output_dir.mkdir(parents=True, exist_ok=True)
+        print(data_tsv + "\n" + auroc_results_df.describe())
+        auroc_results_df.to_csv(
+            results_output_dir
+            / f"results_benchmarking_{Path(data_tsv).stem.replace('benchmarking_','').replace('_molecules','')}.csv"
+        )
+
+    def auroc_performance_on_dir(
+        self,
+        working_dir: Union[Path, str] = Path("."),
+        data_tsv_contains_canonical_smiles: bool = True,
+        smiles_column_title: str = "SMILES",
+        active_flag_column_title: str = "Active",
+        results_output_dir: Optional[Union[Path, str]] = None,
+    ):
+        working_dir = Path(working_dir)
+        csv_files = [Path(f) for f in working_dir.glob("benchmarking_*_molecules.csv")]
+        db_and_tsv = []
+        for csvf in csv_files:
+            target = str(csvf.stem).split("_")[1]
+            db_file_path = (
+                working_dir / f"benchmarking_dud_e_diverse_mols_by_target_{target}.db"
+            )
+            if db_file_path.exists():
+                db_and_tsv.append((db_file_path, csvf))
+        print(db_and_tsv)
+        for db_file_path, csv_file_path in tqdm(db_and_tsv):
+            self.auroc_performance(
+                csv_file_path,
+                data_tsv_contains_canonical_smiles=data_tsv_contains_canonical_smiles,
+                smiles_column_title=smiles_column_title,
+                active_flag_column_title=active_flag_column_title,
+                results_output_dir=results_output_dir,
+            )
diff --git a/src/fepops/dude_preprocessor.py b/src/fepops/dude_preprocessor.py
new file mode 100644
index 0000000..b826947
--- /dev/null
+++ b/src/fepops/dude_preprocessor.py
@@ -0,0 +1,133 @@
+from typing import Union, Optional
+from pathlib import Path
+from tqdm import tqdm
+from fepops import Fepops
+import pandas as pd
+import multiprocessing as mp
+from rdkit import Chem
+
+
+class DudePreprocessor:
+    def __init__(
+        self,
+        *,
+        dude_directory: Union[Path, str] = "data/dude/",
+    ) -> None:
+        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():
+            raise FileNotFoundError(f"Dude dataset not found in path: {self.dude_path}")
+        self.dude_processed_path = self.dude_path / Path("processed")
+        self.dude_processed_path.parent.mkdir(parents=True, exist_ok=True)
+        self.fepops_ob = Fepops()
+
+    def __call__(
+        self,
+    ):
+        self.process()
+
+    def process(
+        self,
+        targets: Optional[Union[str, list[str]]] = None,
+        skip_existing: bool = True,
+    ):
+        if targets is None:
+            dude_targets = [
+                t.parent.name
+                for t in self.dude_path.glob("unprocessed/*/actives_final.ism")
+            ]
+        else:
+            if isinstance(targets, str):
+                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)
+
+    @staticmethod
+    def _parallel_init_worker_desc_gen_shared_fepops_ob():
+        global shared_fepops_ob
+        shared_fepops_ob = Fepops()
+
+    @staticmethod
+    def _parallel_get_rdkit_cansmi(s):
+        global shared_fepops_ob
+        mol = shared_fepops_ob._mol_from_smiles(s)
+        if mol is None:
+            return ""
+        return Chem.MolToSmiles(mol)
+
+    def create_dude_target_csv_data(
+        self,
+        dude_target: Path,
+        actives_file: Path = Path("actives_final.ism"),
+        decoys_file: Path = Path("decoys_final.ism"),
+        seperator: str = " ",
+        skip_existing: bool = True,
+    ):
+        target_output_file = self.dude_processed_path / f"dude_target_{dude_target}.csv"
+        if skip_existing and target_output_file.exists():
+            print(
+                f"Found existing {target_output_file}, skipping due to skip_existing = True, rerun as False to regenerate"
+            )
+        actives = pd.read_csv(
+            self.dude_unprocessed_path / Path(dude_target) / actives_file,
+            sep=seperator,
+            header=None,
+            names=["SMILES", "DUDEID", "CHEMBLID"],
+        )
+        actives["Active"] = 1
+        decoys = pd.read_csv(
+            self.dude_unprocessed_path / Path(dude_target) / decoys_file,
+            sep=seperator,
+            header=None,
+            names=["SMILES", "DUDEID"],
+        )
+        decoys["Active"] = 0
+        df = pd.concat([actives, decoys]).reset_index().drop(columns="index")
+        df["rdkit_canonical_smiles"] = tqdm(
+            mp.Pool(
+                initializer=self._parallel_init_worker_desc_gen_shared_fepops_ob
+            ).imap(self._parallel_get_rdkit_cansmi, df.SMILES, chunksize=100),
+            desc=f"Generating {dude_target} benchmark file",
+            total=len(df),
+        )
+        df.to_csv(target_output_file, index=False)
+
+    def cache_mols_from_csv(
+        self,
+        csv_path: Union[Path, str],
+        rdkit_canonical_smiles_column_header: str = "rdkit_canonical_smiles",
+    ):
+        """Cache mols from a CSV into a db for faster recall later
+
+        Parameters
+        ----------
+        csv_path : Union[Path, str]
+            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"
+        """
+
+        from .fepops_persistent import get_persistent_fepops_storage_object
+
+        for csv_path in (
+            [Path(csv_path)]
+            if csv_path is not None
+            else self.dude_processed_path.glob("dude_target_*.csv")
+        ):
+            df = pd.read_csv(csv_path)
+            if df[rdkit_canonical_smiles_column_header].isnull().values.any():
+                print(
+                    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()])
+            smiles = [
+                s
+                for s in df[rdkit_canonical_smiles_column_header].tolist()
+                if not pd.isnull(s)
+            ]
+            with get_persistent_fepops_storage_object(csv_path.with_suffix(".db")) as f:
+                f.save_descriptors(smiles)
diff --git a/fepops/fepops.py b/src/fepops/fepops.py
similarity index 69%
rename from fepops/fepops.py
rename to src/fepops/fepops.py
index b8a8919..d6278df 100644
--- a/fepops/fepops.py
+++ b/src/fepops/fepops.py
@@ -5,29 +5,59 @@
 from rdkit.Chem import rdMolTransforms
 from sklearn.cluster import KMeans as _SKLearnKMeans
 from fast_pytorch_kmeans import KMeans as _FastPTKMeans
-from sklearn.preprocessing import StandardScaler
 from scipy.spatial.distance import cdist, squareform, pdist
 from scipy.special import softmax
 import numpy as np
 import itertools, zlib
 import torch
-from typing import Union, Optional
+from typing import Union, Optional, Tuple, Literal
+from enum import Enum
+import multiprocessing as mp
+from multiprocessing import SimpleQueue
+
+GetFepopStatusCode = Enum(
+    "GetFepopStatusCode",
+    ["SUCCESS", "FAILED_TO_GENERATE", "FAILED_TO_RETRIEVE", "FAILED_RETRIEVED_NONE"],
+)
 
 
 class Fepops:
-    """Fepops molecular similarity object
+    """Fepops (Feature Points) molecular similarity object
 
-    Fepops allows the comparison of molecules using feature points, see
-    the original publication for more information https://pubs.acs.org/doi/10.1021/jm049654z
-    In short, featurepoints reduce the number of points used to represent a molecule and can
-    be used to compare molecules in the hope of discovering biosimilars based on queries.
+    Fepops allows the comparison of molecules using feature points, see the original
+    publication for more information: https://doi.org/10.1021/jm049654z. In short,
+    featurepoints reduce the number of points used to represent a molecule by combining
+    atoms and their properties. Typically used to compare libraries of small molecules
+    against known actives in the hope of discovering biosimilars based on queries.
 
     Parameters
     ----------
     kmeans_method : str, optional
-        Method which should be used for kmeans calculation, can be
-        one of "sklearn", "pytorch-gpu", or "pytorch-cpu". By
-        default "pytorch-cpu".
+        String literal denoting the method which should be used for kmeans calculations.
+        May be one of "sklearn", "pytorch-gpu", or "pytorch-cpu". If "sklearn" is passed
+        then Scikit-learn's kmeans implementation is used. However a faster
+        implementation from the fast_pytorch_kmeans package can also be used if Pytorch
+        is available and may be run in cpu-only mode, or GPU accelerated mode. Note:
+        GPU accelerated mode should only be used if you are stretching the capabilities
+        in terms of feature points for large molecules.  Small molecules will not
+        benefit from GPU acceleration due to overheads.  By default "pytorch-cpu"
+    max_tautomers : Optional[int], optional
+        Maximum number of tautomers which should be generated. Internally, this
+        implementation of FEPOPS relies upon RDKit's TautomerEnumerator to generate
+        tautomers and may optionally pass in a limit to the number of Tautomers to
+        generate. Unless the molecules (or macromolecules) you are working with generate
+        massive numbers of tautomers, this should be None implying that no limit should
+        be placed on tautomer generation. By default None
+    num_fepops_per_mol : int, optional
+            Number of feature points to use in the representation of a molecule.
+            Literature notes that 7 has been empirically found to be a good number of
+            feature points for performant representations of small molecules. This might
+            be increased if you are dealing with large and very flexible molecules. By
+            default 7
+    num_centroids_per_fepop : int, optional
+            Each fepop is represented by a number of centres, into which atom properties
+            are compressed. Literature notes that this has been empirically determined
+            to be 4 for a performant representation of small molecules. By default 4
 
     Raises
     ------
@@ -37,13 +67,12 @@ class Fepops:
 
     def __init__(
         self,
-        kmeans_method: str = "pytorch-cpu",
-        max_tautomers: Optional[int] = None,
+        kmeans_method: Literal['sklearn', 'pytorch-cpu', 'pytorch-gpu'] = 'sklearn',
         *,
+        max_tautomers: Optional[int] = None,
         num_fepops_per_mol: int = 7,
         num_centroids_per_fepop: int = 4,
     ):
-
         self.implemented_kmeans_methods = ["sklearn", "pytorch-cpu", "pytorch-gpu"]
         self.sort_by_features_col_index_dict = {
             name: sort_order_index
@@ -70,15 +99,15 @@ def __init__(
         self.tautomer_enumerator = MolStandardize.tautomer.TautomerEnumerator(
             **{"max_tautomers": max_tautomers} if max_tautomers is not None else {}
         )
-        self.scaler = StandardScaler()
 
     def _get_k_medoids(
-        self, input_x: np.ndarray, k: int = 7, seed: int = 42
+        self, input_x: np.ndarray, k: int = 7, random_state: int = 42
     ) -> np.ndarray:
         """Select k Fepops from conformers
 
-        Gets k mediods from conformers (and conformers of tautomers) which
-        are representative of the molcule.
+        Gets k mediods from conformers (and tautomers) which are representative
+        of the molecule as a function of conformer and tautomer states by virtue
+        of chosen fepops being diverse.
 
         Parameters
         ----------
@@ -86,11 +115,15 @@ def _get_k_medoids(
             The pharmacophore features of all conformers.
         k : int
             The number of medoids for clustering. By default 7.
+        random_state : int
+            Integer to use as a random state when seeding the random number
+            generator.  By default 42.
 
         Returns
         -------
         np.ndarray
-            The final Fepops descriptors of the k representative conformers.
+            The final Fepops descriptors comprised of k representative
+            conformers/tautomers.
         """
         input_x = np.unique(input_x, axis=0)
 
@@ -100,7 +133,7 @@ def _get_k_medoids(
         point_to_centroid_map = np.ones(input_x.shape[0])
         point_to_centroid_map_prev = np.zeros_like(point_to_centroid_map)
 
-        np_rng = np.random.default_rng(seed=seed)
+        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), :
         ]
@@ -118,14 +151,11 @@ def _get_k_medoids(
                     point_to_centroid_map[chosen_x_point] = i
                 medoids[i] = np.median(input_x[point_to_centroid_map == i], axis=0)
         # Sorting at this stage for reproducibility with existing pregenerated
-        # descriptor sets.
+        # descriptor sets and convention with early fepops versions which relied
+        # upon fepops being sorted by charge (before moving to the newer CombiAlign
+        # scoring algorithm)
         return medoids[np.lexsort(medoids.T[::-1])]
 
-    def annotate_atom_idx(self, mol: Chem.rdchem.Mol):
-        for i, atom in enumerate(mol.GetAtoms()):
-            # For each atom, set the property "molAtomMapNumber" to a custom number, let's say, the index of the atom in the molecule
-            atom.SetProp("molAtomMapNumber", str(atom.GetIdx()))
-
     def _calculate_atomic_logPs(self, mol: Chem.rdchem.Mol) -> dict:
         """Calculate logP contribution for each of atom in a molecule
 
@@ -159,7 +189,7 @@ def _calculate_partial_charges(self, mol: Chem.rdchem.Mol) -> dict:
         dict
             A dictionary containing all atom symobls with their charges.
         """
-        Chem.rdPartialCharges.ComputeGasteigerCharges(mol)
+        Chem.rdPartialCharges.ComputeGasteigerCharges(mol, throwOnParamFailure=True)
         return {
             atom.GetIdx(): float(atom.GetProp("_GasteigerCharge"))
             for atom in mol.GetAtoms()
@@ -238,7 +268,9 @@ def _perform_kmeans(
         """
         if kmeans_method == "sklearn":
             kmeans = _SKLearnKMeans(
-                n_clusters=num_centroids, random_state=seed, n_init="auto"
+                n_clusters=num_centroids,
+                random_state=seed,
+                n_init="auto",
             ).fit(atom_coords)
             centroid_coors = kmeans.cluster_centers_
             instance_cluster_labels = kmeans.labels_
@@ -284,21 +316,13 @@ def _get_centroid_distances(
         np.ndarray
             Ordered centroid distances
         """
-
-        if is_distance_matrix:
-            distance_matrix = centroid_coords_or_distmat
+        if not is_distance_matrix:
+            dmat = squareform(pdist(centroid_coords_or_distmat))
         else:
-            distance_matrix = squareform(pdist(centroid_coords_or_distmat))
-        distances = np.array(
-            [distance_matrix[0, distance_matrix.shape[0] - 1]]
-            + [
-                ele
-                for arr in [
-                    distance_matrix.diagonal(i)
-                    for i in range(1, distance_matrix.shape[0] - 1)
-                ]
-                for ele in arr
-            ]
+            dmat = centroid_coords_or_distmat.copy()
+        distances = np.hstack(
+            [dmat[0, -1]]
+            + [np.diagonal(dmat, offset=k) for k in range(1, dmat.shape[0] - 1)]
         )
         return distances
 
@@ -326,11 +350,8 @@ def _mol_from_smiles(self, smiles_string: str) -> Chem.rdchem.Mol:
             try:
                 mol = Chem.MolFromSmiles(smiles_string, sanitize=False)
             except:
-                pass
-        if mol is None:
-            print(
                 f"Could not parse smiles to a valid molecule, smiles was: {smiles_string}"
-            )
+                return None
         return mol
 
     def _get_flanking_atoms(
@@ -542,7 +563,6 @@ def get_centroid_pharmacophoric_features(
         pharmacophore_features_arr = np.empty(shape=[self.num_centroids_per_fepop, 4])
         for centroid in range(self.num_centroids_per_fepop):
             centroid_atomic_id = np.where(instance_cluster_labels == centroid)[0]
-
             sum_of_logP = self._sum_of_atomic_features_by_centroids(
                 atomic_logP_dict, centroid_atomic_id
             )
@@ -579,8 +599,10 @@ def get_centroid_pharmacophoric_features(
         return pharmacophore_features_arr
 
     def get_fepops(
-        self, mol: Union[str, None, Chem.rdchem.Mol]
-    ) -> Union[np.ndarray, None]:
+        self,
+        mol: Union[str, None, Chem.rdchem.Mol],
+        is_canonical: bool = False,
+    ) -> Tuple[GetFepopStatusCode, Union[np.ndarray, None]]:
         """Get Fepops descriptors
 
         This method returns Fepops descriptors from a smiles string.
@@ -588,48 +610,101 @@ def get_fepops(
 
         Parameters
         ----------
-        smiles_string : str
-            SMILES string of an input molecule.
+        mol : Union[str, None, Chem.rdchem.Mol]
+            Molecule as a SMILES string or RDKit molecule. Can also be None,
+            in which case a failure error status is returned along with None
+            in place of the requested Fepops descriptors.
 
         Returns
         -------
-        np.ndarray
-            A Numpy array containing the calculated Fepops descriptors of an input molecule.
+        Tuple[GetFepopStatusCode, Union[np.ndarray, None]]
+            Returns a tuple, with the first value being a GetFepopStatusCode
+            (enum) denoting SUCCESS or FAILED_TO_GENERATE. The second tuple
+            element is either None (if unsuccessful), or a np.ndarray containing
+            the calculated Fepops descriptors of the requested input molecule.
         """
+        original_smiles = None
+        if isinstance(mol, np.ndarray):
+            return GetFepopStatusCode.SUCCESS, mol
         if isinstance(mol, str):
+            original_smiles = mol
             mol = self._mol_from_smiles(mol)
         if mol is None:
-            return None
-
-        mol = Chem.AddHs(mol)
+            print(
+                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(
-                f"Number of heavy atoms (:{Lipinski.HeavyAtomCount(mol)}) below requested feature points (:{self.num_centroids_per_fepop})"
+                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 None
+            return GetFepopStatusCode.FAILED_TO_GENERATE, None
+        mol = Chem.AddHs(mol)
 
         tautomers_list = self.tautomer_enumerator.enumerate(mol)
         each_mol_with_all_confs_list = []
         for index, t_mol in enumerate(tautomers_list):
             conf_list = self.generate_conformers(t_mol)
             each_mol_with_all_confs_list.extend(conf_list)
-
         if each_mol_with_all_confs_list == []:
-            return None
+            print(
+                f"Failed to generate conformers/tautomers {' for '+original_smiles if original_smiles is not None else ''}"
+            )
+            return GetFepopStatusCode.FAILED_TO_GENERATE, None
 
-        pharmacophore_feature_all_confs = np.array(
-            [
-                self.get_centroid_pharmacophoric_features(
-                    each_mol,
-                    kmeans_method_str=self.kmeans_method_str,
-                )
-                for each_mol in each_mol_with_all_confs_list
-            ]
-        )
+        try:
+            pharmacophore_feature_all_confs = np.array(
+                [
+                    self.get_centroid_pharmacophoric_features(
+                        each_mol,
+                        kmeans_method_str=self.kmeans_method_str,
+                    )
+                    for each_mol in each_mol_with_all_confs_list
+                ]
+            )
+        except ValueError as e:
+            if original_smiles is not None:
+                print(f"Failed molecule had SMILES: {original_smiles}")
+            print(e)
+            return GetFepopStatusCode.FAILED_TO_GENERATE, None
 
-        return self._get_k_medoids(
+        medoids = self._get_k_medoids(
             pharmacophore_feature_all_confs, self.num_fepops_per_mol
         )
+        return GetFepopStatusCode.SUCCESS, medoids
+
+    def _score_scaler(self, x1: np.ndarray, x2: np.ndarray) -> float:
+        """Score function for the similarity calculation
+
+        The score function for the similarity calculation on the FEPOPS descriptors.
+
+        Parameters
+        ----------
+        x1 : np.ndarray
+            A Numpy array containing the FEPOPS descriptors 1.
+        x2 : np.ndarray
+            A Numpy array containing the FEPOPS descriptors 2.
+
+        Returns
+        -------
+        float
+            The FEPOPS similarity score (Pearson correlation).
+        """
+        x1 = self.scaler.fit_transform(x1.reshape(-1, 1))
+        x2 = self.scaler.fit_transform(x2.reshape(-1, 1))
+        return np.corrcoef(x1.flatten(), x2.flatten())[0, 1]
+
+    def pairwise_correlation(self, A, B):
+        am = A - np.mean(A, axis=0, keepdims=True)
+        bm = B - np.mean(B, axis=0, keepdims=True)
+        return (
+            am.T
+            @ bm
+            / (
+                np.sqrt(np.sum(am**2, axis=0, keepdims=True)).T
+                * np.sqrt(np.sum(bm**2, axis=0, keepdims=True))
+            )
+        )
 
     def _score_combialign(self, x1: np.ndarray, x2: np.ndarray):
         """Score fepops using CombiAlign
@@ -698,43 +773,60 @@ def _score_combialign(self, x1: np.ndarray, x2: np.ndarray):
 
         # Rebuild x2 distances to squareform matrix, then reorder as per the best permutation and extract in required FEPOPS order.
         dmat = np.zeros((self.num_centroids_per_fepop, self.num_centroids_per_fepop))
+
         dmat[0, -1] = x2_dists[0]
-        dmat[-1, 0] = x2_dists[0]
-        rows, cols = np.diag_indices_from(dmat)
-        for (r, c), v in zip(
-            [
-                (x, y)
-                for d in [
-                    np.stack((rows[:-i], cols[i:]), axis=1)
-                    for i in range(1, dmat.shape[0] - 1)
-                ]
-                for x, y in d
-            ],
-            x2_dists,
-        ):
-            dmat[r, c] = v
-            dmat[c, r] = v
+        dist_position = 1
+        for offset in range(1, self.num_centroids_per_fepop - 1):
+            num_in_diagonal = self.num_centroids_per_fepop - offset
+            dmat += np.diag(
+                x2_dists[dist_position : dist_position + num_in_diagonal], k=offset
+            )
+            dist_position += num_in_diagonal
+        # Not required, but for completeness, copy upper triangle to lower
+        # triangle of dmat matrix
+        dmat = dmat + dmat.T - np.diag(np.diag(dmat))
 
         # Reorder the distance matrix using best permutation
+
         new_dmat = np.zeros_like(dmat)
         for i, p in enumerate(best_permutaion):
             for j in range(dmat.shape[0]):
                 if i == j:
                     continue
                 new_dmat[i, j] = dmat[p, best_permutaion[j]]
-
         distances = self._get_centroid_distances(new_dmat, is_distance_matrix=True)
 
         # Reform x2 with reordered medoids and medoid distances
         x2 = np.hstack([x2_desc[[best_permutaion]].flatten(), distances])
-
         # Apply softmax and return pearson correlation between the two
-        return np.corrcoef(softmax(x1), softmax(x2))[0, 1]
+
+        # return np.corrcoef(softmax(x1), softmax(x2))[0, 1]
+        # return np.corrcoef(x1, x2)[0, 1]
+        # x1 = self.scaler.fit_transform(x1.reshape(-1, 1))
+        # x2 = self.scaler.fit_transform(x2.reshape(-1, 1))
+        return self.pairwise_correlation(
+            (x1 - x1.mean()) / x1.std(), (x2 - x2.mean()) / x2.std()
+        )
+        # return np.corrcoef(x1.flatten(), x2.flatten())[0, 1]
+        # A= [0.4, 0.6]
+        # SoftA = []
+
+    def _init_worker_calc_similarity(self, query_descriptors_, candidate_descriptors_):
+        global shared_query_descriptors, shared_candidate_descriptors
+        shared_query_descriptors = query_descriptors_
+        shared_candidate_descriptors = candidate_descriptors_
+
+    def _work_calc_similarity(self, candidate_descriptor_i):
+        global shared_query_descriptors, shared_candidate_descriptors
+        return self.calc_similarity(
+            shared_query_descriptors,
+            shared_candidate_descriptors[candidate_descriptor_i],
+        )
 
     def calc_similarity(
         self,
         query: Union[np.ndarray, str, None],
-        candidate: Union[np.ndarray, str, None],
+        candidate: Union[np.ndarray, str, None, list[np.ndarray, str, None]],
     ) -> float:
         """Calculate FEPOPS similarity
 
@@ -745,25 +837,40 @@ def calc_similarity(
         query : Union[np.ndarray, str]
             A Numpy array containing the FEPOPS descriptors of the query molecule
             or a smiles string from which to generate FEPOPS descriptors for the
-            query molecule.
-        candidate : Union[np.ndarray, str]
+            query molecule. Can also be None, in which case, np.nan is returned
+            as a score.
+        candidate : Union[np.ndarray, str, None, list[np.ndarray, str, None]],
             A Numpy array containing the FEPOPS descriptors of the candidate
             molecule or a smiles string from which to generate FEPOPS descriptors
-            for the query molecule.
+            for the candidate molecule.  Can also be None, in which case, np.nan is
+            returned as a score, or a list of any of these. If it is a list,
+            then a list of scores against the single candidate is returned.
 
         Returns
         -------
         float
             Fepops similarity between two molecules
         """
-        if isinstance(query, str):
-            query = self.get_fepops(query)
-        if isinstance(candidate, str):
-            candidate = self.get_fepops(candidate)
-        if any(x is None for x in (query, candidate)):
-            raise ValueError(
-                f"Unable to calculate similarity due to NoneType found in the fepops features:(query, candidate)=({type(query)}, {type(candidate)})"
-            )
+        if not isinstance(query, np.ndarray):
+            query_status, query = self.get_fepops(query)
+            if query_status != GetFepopStatusCode.SUCCESS:
+                return np.nan
+
+        if isinstance(candidate, list):
+            shared_queue = SimpleQueue()
+            scores = []
+            with mp.Pool(
+                initializer=self._init_worker_calc_similarity,
+                initargs=(query, candidate),
+            ) as pool:
+                parallel_results = list(
+                    pool.imap(self._work_calc_similarity, (range(len(candidate))))
+                )
+            return parallel_results
+        if not isinstance(candidate, np.ndarray):
+            candidate_status, candidate = self.get_fepops(candidate)
+            if candidate_status != GetFepopStatusCode.SUCCESS:
+                return np.nan
         return np.max(cdist(query, candidate, metric=self._score_combialign))
 
     def __call__(
diff --git a/fepops/fepops_persistent/__init__.py b/src/fepops/fepops_persistent/__init__.py
similarity index 100%
rename from fepops/fepops_persistent/__init__.py
rename to src/fepops/fepops_persistent/__init__.py
diff --git a/src/fepops/fepops_persistent/fepops_persistent_abc.py b/src/fepops/fepops_persistent/fepops_persistent_abc.py
new file mode 100755
index 0000000..85cfabd
--- /dev/null
+++ b/src/fepops/fepops_persistent/fepops_persistent_abc.py
@@ -0,0 +1,324 @@
+from abc import ABCMeta, abstractmethod
+from pathlib import Path
+from typing import Union
+
+import numpy as np
+from rdkit import Chem
+from scipy.spatial.distance import cdist, pdist, squareform
+from tqdm import tqdm
+from fepops.fepops import GetFepopStatusCode
+from ..fepops import Fepops
+import multiprocessing as mp
+
+
+class FepopsPersistentAbstractBaseClass(metaclass=ABCMeta):
+    """Abstract base class for persistent fepops storage
+
+    New storage methods may be implemented as demonstrated in fepopsdb_json.py
+    by extending this abstract base class which provides some required
+    functionality like:
+
+    - save_descriptors(smiles: Union[str, Path, list[str]]), to save a smiles file/list of smiles to the persistent storage
+    - get_cansmi_to_mol_dict_not_in_database(smiles: Union[str, Path, list[str]]), to retrieve a unique dictionary with canonical smiles as keys not already stored in the database and rdkit mol objects as values.
+
+    When writing your own persistent storage methods, you must override the following methods:
+
+    add_fepop(rdkit_canonical_smiles: str, fepops: np.ndarray)
+    ----------------------------------------------------------
+    Add the fepop to persistent storage. super().add_fepop may be called by the overridden function to perform type checks on arguments.
+
+    fepop_exists(rdkit_canonical_smiles: str)
+    -----------------------------------------
+    Return True if the canonical smiles is already in the database, and False if not. super().fepop_exists may be called by the overridden
+    function to perform type checks on arguments.
+
+    get_fepops(rdkit_canonical_smiles: str)
+    --------------------------------------
+    Return a fepop from persistent storage. If it does not exist, then generate it by calling self.fepops_object.get_fepops which is supplied
+    by this base class. super().get_fepops may be called by the overridden function to perform type checks on arguments.
+    With this function in place it allows interface compatibility with a standard Fepops object.
+
+    Inheriting functions may also define __enter__ and __exit__ methods for use with context handlers. If none are defined, then empty ones
+    are provided. This can be useful in doing things like writing out large files after descriptor generation if incremental writes are not
+    possible, like in the case of the FepopsDBJSON child class.
+
+
+    Parameters
+    ----------
+    database_file : Union[str, Path]
+            File to use for persistent storage.
+    kmeans_method : str, optional
+            Method which should be used for kmeans calculation by
+            fepops objects, can be one of "sklearn", "pytorch-gpu",
+            or "pytorch-cpu".
+    parallel : bool, optional
+            Run in parallel (using joblib), by default True
+    n_jobs : int, optional
+            Number of jobs to be spawned with joblib. If -1, then use
+            all available cores. By default -1.
+    """
+
+    @staticmethod
+    def _parallel_init_worker_get_cansmi_mol_tuple(smiles_is_rdkit_canonical):
+        global sirdkc
+        sirdkc = smiles_is_rdkit_canonical
+
+    @staticmethod
+    def _parallel_get_cansmi_tuple(m):
+        global sirdkc
+        return FepopsPersistentAbstractBaseClass._get_can_smi_mol_tuple(
+            m, smiles_guaranteed_rdkit_canonical=sirdkc
+        )
+
+    @staticmethod
+    def _parallel_init_worker_desc_gen_shared_fepops_ob(fepops_object):
+        global shared_fepops_ob
+        shared_fepops_ob = fepops_object
+
+    @staticmethod
+    def _parallel_get_gen_fepops_descriptors(m):
+        global shared_fepops_ob
+        return m[0], shared_fepops_ob.get_fepops(m[1])
+
+    @abstractmethod
+    def __init__(
+        self,
+        database_file: Union[str, Path],
+        kmeans_method: str = "pytorch-cpu",
+        parallel: bool = True,
+        n_jobs: int = -1,
+    ):
+        self.database_file = Path(database_file)
+        self.fepops_object = Fepops(kmeans_method=kmeans_method)
+        self.parallel = parallel
+        self.n_jobs = n_jobs
+
+    def save_descriptors(
+        self,
+        smiles: Union[str, Path, list[str]],
+        add_failures_to_database: bool = True,
+        smiles_guaranteed_rdkit_canonical: bool = False,
+        fepops_object_constructor_kwargs: dict = {},
+    ):
+        """Pregenerate FEPOPS descriptors for a set of SMILES strings"""
+
+        canonical_smiles_to_mol_dict = self.get_cansmi_to_mol_dict_not_in_database(
+            smiles, smiles_guaranteed_rdkit_canonical=smiles_guaranteed_rdkit_canonical
+        )
+        if len(canonical_smiles_to_mol_dict) == 0:
+            print("Nothing to add to database")
+            return
+        if not self.parallel:
+            for rdkit_canonical_smiles, mol in tqdm(
+                canonical_smiles_to_mol_dict.items(), desc="Generating fepops"
+            ):
+                status, fepops_array = self.fepops_object.get_fepops(mol)
+                if status == GetFepopStatusCode.SUCCESS or add_failures_to_database:
+                    self.add_fepop(rdkit_canonical_smiles, fepops_array)
+            print(
+                f"Added {len(canonical_smiles_to_mol_dict)} new molecues to the database ({self.database_file})"
+            )
+        else:  # Do it in parallel
+            n_successes = 0
+            n_failures = 0
+            for rdkit_canonical_smiles, (status, new_fepop) in tqdm(
+                mp.Pool(
+                    # processes=min(len(canonical_smiles_to_mol_dict), mp.cpu_count()),
+                    processes=2,
+                    initializer=self._parallel_init_worker_desc_gen_shared_fepops_ob,
+                    initargs=(Fepops(**fepops_object_constructor_kwargs),),
+                ).imap(
+                    self._parallel_get_gen_fepops_descriptors,
+                    canonical_smiles_to_mol_dict.items(),
+                ),
+                desc="Generating descriptors (parallel)",
+                total=len(canonical_smiles_to_mol_dict),
+            ):
+                if status == GetFepopStatusCode.SUCCESS:
+                    n_successes += 1
+                else:
+                    n_failures += 1
+                if status == GetFepopStatusCode.SUCCESS or add_failures_to_database:
+                    self.add_fepop(rdkit_canonical_smiles, new_fepop)
+            print(
+                f"Successfully added {n_successes} new molecues to the database ({self.database_file}), {n_failures} failed"
+            )
+
+    @abstractmethod
+    def add_fepop(self, rdkit_canonical_smiles: str, fepops: np.ndarray):
+        """Add canonical smiles and fepop to database.
+
+        Must be overridden
+        """
+        if not isinstance(rdkit_canonical_smiles, str):
+            raise ValueError(
+                f"Expected an rdkit canonical smiles string, but a {type(rdkit_canonical_smiles)} was passed"
+            )
+        if not isinstance(fepops, np.ndarray):
+            raise ValueError(f"Expected a fepop, but a {type(fepops)} was passed")
+
+    @abstractmethod
+    def get_fepops(
+        self, smiles: Union[str, Chem.rdchem.Mol, np.ndarray], is_canonical: bool = True
+    ) -> None:
+        if not isinstance(smiles, (str, Chem.rdchem.Mol, np.ndarray)):
+            raise ValueError(
+                f"Expected an rdkit canonical smiles string, rdkit mol, or a numpy array of descriptors but a {type(smiles)} was passed"
+            )
+
+    @abstractmethod
+    def fepop_exists(self, rdkit_canonical_smiles: str) -> bool:
+        """Return True if canonical smiles already exist in the database"""
+        if not isinstance(rdkit_canonical_smiles, str):
+            raise ValueError(
+                f"Expected an rdkit canonical smiles string, but a {type(rdkit_canonical_smiles)} was passed"
+            )
+
+    @staticmethod
+    def _get_can_smi_mol_tuple(s: str, smiles_guaranteed_rdkit_canonical: bool = False):
+        try:
+            mol = Chem.MolFromSmiles(s)
+        except:
+            try:
+                mol = Chem.MolFromSmiles(s, sanitize=False)
+            finally:
+                mol = None
+        if mol is None:
+            print(f"Could not parse smiles to a valid molecule, smiles was: {s}")
+            return (s, mol)
+        if smiles_guaranteed_rdkit_canonical:
+            return (s, mol)
+        else:
+            return (Chem.MolToSmiles(mol), mol)
+
+    def get_cansmi_to_mol_dict_not_in_database(
+        self,
+        smiles: Union[str, Path, list[str]],
+        smiles_guaranteed_rdkit_canonical: bool = False,
+    ):
+        if isinstance(smiles, str):
+            smiles = Path(smiles)
+        if isinstance(smiles, Path):
+            if smiles.exists():
+                smiles = [
+                    s.strip() for s in open(smiles).readlines() if len(s.strip()) > 0
+                ]
+            else:
+                raise ValueError(
+                    f"smiles file ({smiles}) not found. If you are passing smiles, place it into a list first"
+                )
+        if not isinstance(smiles, list):
+            raise ValueError(
+                "smiles should be a str or Path denoting the location of a smiles file, or a list of smiles"
+            )
+        smiles = list(set(smiles))
+        print(f"Got {len(smiles)} unique SMILES strings")
+
+        if not self.parallel:
+            # Ensure unique (canonical, also storing intermediate mol)
+            canonical_smiles_to_mol_dict = dict(
+                self._get_can_smi_mol_tuple(
+                    s,
+                    smiles_guaranteed_rdkit_canonical=smiles_guaranteed_rdkit_canonical,
+                )
+                for s in tqdm(smiles, desc="Uniquifying input smiles (non-parallel)")
+            )
+        else:
+            tmp_res_list = []
+            # Ensure unique (canonical, also storing intermediate mol)
+            for res in tqdm(
+                mp.Pool(
+                    initializer=self._parallel_init_worker_get_cansmi_mol_tuple,
+                    initargs=(smiles_guaranteed_rdkit_canonical,),
+                    processes=max(1, min(len(smiles) / 50, mp.cpu_count())),
+                ).map(self._parallel_get_cansmi_tuple, smiles, chunksize=1),
+                desc="Uniquifying input smiles (parallel)",
+                total=len(smiles),
+            ):
+                tmp_res_list.append(res)
+
+            canonical_smiles_to_mol_dict = dict(tmp_res_list)
+
+            del tmp_res_list
+        # Make sure none are already in the database
+        canonical_smiles_to_mol_dict = {
+            cansmi: mol
+            for cansmi, mol in tqdm(
+                canonical_smiles_to_mol_dict.items(),
+                desc="Checking if mols already exist in the database",
+            )
+            if not self.fepop_exists(cansmi)
+        }
+        print(
+            f"Got {len(canonical_smiles_to_mol_dict)} unique molecules not already in the database"
+        )
+
+        return canonical_smiles_to_mol_dict
+
+    def calc_similarity(
+        self,
+        fepops_features_1: Union[np.ndarray, str, None],
+        fepops_features_2: Union[np.ndarray, str, None],
+        is_canonical=True,
+    ):
+        """Calculate FEPOPS similarity
+
+        A static method for calculating molecular similarity based on their FEPOPS descriptors.
+
+        Parameters
+        ----------
+        fepops_features_1 : Union[np.ndarray, str, None]
+            A Numpy array containing the FEPOPS descriptors of the query molecule
+            or a smiles string from which to generate FEPOPS descriptors for the
+            query molecule.
+        fepops_features_2 : Union[np.ndarray, str, None, list[np.ndarray, str, None]]
+            A Numpy array containing the FEPOPS descriptors of the candidate
+            molecule or a smiles string from which to generate FEPOPS descriptors
+            for the candidate molecule.  Can also be None, in which case, np.nan is
+            returned as a score, or a list of any of these. If it is a list,
+            then a list of scores against the single candidate is returned.
+
+        Returns
+        -------
+        float
+                Fepops similarity between two molecules
+        """
+        if fepops_features_1 is None:
+            return np.nan
+        if isinstance(fepops_features_1, (str, Chem.rdchem.Mol)):
+            status, fepops_features_1 = self.get_fepops(
+                fepops_features_1, is_canonical=is_canonical
+            )
+            if status != GetFepopStatusCode.SUCCESS:
+                return np.nan
+
+        if isinstance(fepops_features_2, list):
+            new_fepops_features_2 = []
+            for item in fepops_features_2:
+                status, fpop = self.get_fepops(item, is_canonical=is_canonical)
+                new_fepops_features_2.append(
+                    fpop if status == GetFepopStatusCode.SUCCESS else None
+                )
+            return self.fepops_object.calc_similarity(
+                fepops_features_1, new_fepops_features_2
+            )
+
+        if isinstance(fepops_features_2, (str, Chem.rdchem.Mol)):
+            status, fepops_features_2 = self.get_fepops(
+                fepops_features_2, is_canonical=is_canonical
+            )
+            if status != GetFepopStatusCode.SUCCESS:
+                return np.nan
+        if any(x is None for x in (fepops_features_1, fepops_features_2)):
+            return np.nan
+        score = self.fepops_object.calc_similarity(fepops_features_1, fepops_features_2)
+        return score if score is not None else np.nan
+
+    def write(self):
+        pass
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        self.write()
diff --git a/fepops/fepops_persistent/fepopsdb_json.py b/src/fepops/fepops_persistent/fepopsdb_json.py
similarity index 95%
rename from fepops/fepops_persistent/fepopsdb_json.py
rename to src/fepops/fepops_persistent/fepopsdb_json.py
index a830c67..c5db044 100755
--- a/fepops/fepops_persistent/fepopsdb_json.py
+++ b/src/fepops/fepops_persistent/fepopsdb_json.py
@@ -43,10 +43,10 @@ def add_fepop(self, rdkit_canonical_smiles: str, fepops: Union[np.ndarray, None]
             ).decode("ascii")
             self._db_changed = True
 
-    def get_fepop(
+    def get_fepops(
         self, smiles: str, is_canonical: bool = False
     ) -> Union[np.ndarray, None]:
-        super().get_fepop(smiles=smiles)
+        super().get_fepops(smiles=smiles)
         if not is_canonical:
             smiles, mol = self._get_can_smi_mol_tuple(smiles)
         if self.fepop_exists(rdkit_canonical_smiles=smiles):
diff --git a/fepops/fepops_persistent/fepopsdb_sqlite.py b/src/fepops/fepops_persistent/fepopsdb_sqlite.py
similarity index 69%
rename from fepops/fepops_persistent/fepopsdb_sqlite.py
rename to src/fepops/fepops_persistent/fepopsdb_sqlite.py
index b9cefd6..6f48348 100755
--- a/fepops/fepops_persistent/fepopsdb_sqlite.py
+++ b/src/fepops/fepops_persistent/fepopsdb_sqlite.py
@@ -2,9 +2,9 @@
 import sqlite3
 from pathlib import Path
 from typing import Union
-
+from rdkit import Chem
 import numpy as np
-
+from fepops.fepops import GetFepopStatusCode
 from .fepops_persistent_abc import FepopsPersistentAbstractBaseClass
 
 
@@ -89,20 +89,35 @@ def fepop_exists(self, rdkit_canonical_smiles: str) -> bool:
             return False
         return True
 
-    def get_fepop(self, smiles, is_canonical=False) -> Union[np.ndarray, None]:
-        super().get_fepop(smiles=smiles)
-        rdkit_canonical_smiles, mol = self._get_can_smi_mol_tuple(
-            smiles, is_canonical=is_canonical
-        )
+    def get_fepops(self, smiles, is_canonical=False) -> Union[np.ndarray, None]:
+        super().get_fepops(smiles=smiles)
+
+        if isinstance(smiles, str):
+            if is_canonical:
+                rdkit_canonical_smiles = smiles
+                mol = None
+            else:
+                rdkit_canonical_smiles, mol = self._get_can_smi_mol_tuple(
+                    smiles, smiles_guaranteed_rdkit_canonical=is_canonical
+                )
+        elif isinstance(smiles, Chem.rdchem.Mol):
+            mol = smiles
+            rdkit_canonical_smiles = Chem.MolToSmiles(mol)
+        else:
+            # At this point is is guaranteed to be np.ndarray (type checking
+            # performed by super), so just return the array (smiles) and
+            # success.
+            return GetFepopStatusCode.SUCCESS, smiles
+
         if self.fepop_exists(rdkit_canonical_smiles):
             res = self.cur.execute(
                 f"""SELECT fepops FROM fepops_lookup_table where cansmi="{rdkit_canonical_smiles}" """
             )
             fepop = res.fetchone()[0]
             if np.isnan(fepop).any():
-                return None
+                return GetFepopStatusCode.FAILED_RETRIEVED_NONE, None
             else:
-                return fepop.reshape(
+                return GetFepopStatusCode.SUCCESS, fepop.reshape(
                     -1,
                     (
                         self.fepops_object.num_centroids_per_fepop
@@ -111,8 +126,14 @@ def get_fepop(self, smiles, is_canonical=False) -> Union[np.ndarray, None]:
                     + self.fepops_object.num_distances_per_fepop,
                 )
         else:
-            fepops_descriptors = self.fepops_object.get_fepops(mol)
-            self.add_fepop(
-                rdkit_canonical_smiles=rdkit_canonical_smiles, fepops=fepops_descriptors
-            )
-            return fepops_descriptors
+            if mol is None:
+                mol = rdkit_canonical_smiles
+            status, fepops_descriptors = self.fepops_object.get_fepops(mol)
+            if status == GetFepopStatusCode.SUCCESS:
+                self.add_fepop(
+                    rdkit_canonical_smiles=rdkit_canonical_smiles,
+                    fepops=fepops_descriptors,
+                )
+                return status, fepops_descriptors
+            else:
+                return status, None
diff --git a/src/fepops/fepops_persistent/utils.py b/src/fepops/fepops_persistent/utils.py
new file mode 100644
index 0000000..6eabca5
--- /dev/null
+++ b/src/fepops/fepops_persistent/utils.py
@@ -0,0 +1,14 @@
+from . import FepopsDBSqlite, FepopsDBJSON
+from typing import Union
+from pathlib import Path
+
+
+def get_persistent_fepops_storage_object(
+    database_file: Union[str, Path], parallel: bool = True, n_jobs: int = -1
+):
+    if str(database_file).endswith(
+        (".sqlite", ".sqlite3", ".db", ".db3", ".s3db", ".sl3")
+    ):
+        return FepopsDBSqlite(database_file, parallel=parallel, n_jobs=n_jobs)
+    if str(database_file).endswith((".json")):
+        return FepopsDBJSON(database_file, parallel=parallel, n_jobs=n_jobs)


Test patch:
diff --git a/test/tests.py b/test/tests.py
index e6088b0..6afb4b5 100644
--- a/test/tests.py
+++ b/test/tests.py
@@ -1,6 +1,9 @@
 """Initial framework of FEPOPS test suite"""
 
-from fepops.fepops import Fepops
+from pathlib import Path
+from fepops import Fepops
+from fepops.fepops import GetFepopStatusCode
+from fepops.fepops_persistent import get_persistent_fepops_storage_object
 
 
 def test_k_medoids_simple(fepops_1024_unclustered):
@@ -9,21 +12,51 @@ def test_k_medoids_simple(fepops_1024_unclustered):
     assert f._get_k_medoids(fepops_1024_unclustered).shape[0] == 7
 
 
-def test_paracetamol_vs_all(test_mol_smiles):
+def test_all_test_mols_rank_themselves_top(test_mol_smiles):
     """Test similarity of paracetamol to itself, diclofenac, ibuprofen, nicotine and vanilin"""
     f = Fepops()
     fepops_dict = {
         mol_name: f.get_fepops(mol_smiles)
         for mol_name, mol_smiles in test_mol_smiles.items()
     }
-    scores = {
-        f"Paracetamol_{mol_name}_similarity": f.calc_similarity(
-            fepops_dict["Paracetamol"], fepop
-        )
-        for mol_name, fepop in fepops_dict.items()
+    assert all(
+        [t[0] == GetFepopStatusCode.SUCCESS for t in fepops_dict.values()]
+    ), "Could not generate valid fepops for all molecules within the test set"
+    scores = {}
+    for query in fepops_dict.keys():
+        scores[query] = {
+            f"{query}_{candidate_name}_similarity": f.calc_similarity(
+                fepops_dict[query][1], candidate_fepops
+            )
+            for candidate_name, (_, candidate_fepops) in fepops_dict.items()
+        }
+    sorted_results = {
+        query: {
+            k: v
+            for k, v in sorted(
+                scores[query].items(), key=lambda item: item[1], reverse=True
+            )
+        }
+        for query in fepops_dict
     }
-    sorted_scores = {
-        k: v for k, v in sorted(scores.items(), key=lambda item: item[1], reverse=True)
+    for query, results in sorted_results.items():
+        assert query in next(iter(results.keys())).split("_")[1]
+
+
+def test_fepops_persistent_sqlite(test_mol_smiles, tmp_path):
+    """Test creating and saving fepops descriptors to a sqlite3 database"""
+
+    d = tmp_path / "test_dir"
+    d.mkdir()
+    file = d / "test.db"
+    f = get_persistent_fepops_storage_object(database_file=str(file.name))
+    fepops_dict = {
+        mol_name: f.get_fepops(mol_smiles)[1]
+        for mol_name, mol_smiles in test_mol_smiles.items()
+    }
+    fepops_dict2 = {
+        mol_name: f.get_fepops(mol_smiles)[1]
+        for mol_name, mol_smiles in test_mol_smiles.items()
     }
-    it = iter(sorted_scores.keys())
-    assert "Paracetamol" in next(it).split("_")[1]
+    file.unlink(missing_ok=True)
+    assert fepops_dict.keys() == fepops_dict2.keys()