diff --git a/aiida_gaussian_datatypes/basisset/cli.py b/aiida_gaussian_datatypes/basisset/cli.py
index 00ceb20..2ccc864 100644
--- a/aiida_gaussian_datatypes/basisset/cli.py
+++ b/aiida_gaussian_datatypes/basisset/cli.py
@@ -73,7 +73,8 @@ def cli():
     multiple=True,
     help="filter by a tag (all tags must be present if specified multiple times)")
 @click.option(
-    'fformat', '-f', '--format', type=click.Choice(['cp2k']), default='cp2k',
+    'fformat', '-f', '--format', type=click.Choice(['cp2k',
+                                                    'nwchem']), default='cp2k',
     help="the format of the basis set file")
 @click.option(
     '--duplicates',
@@ -92,6 +93,8 @@ def import_basisset(basisset_file, fformat, sym, tags, duplicates, group):
 
     loaders = {
         "cp2k": BasisSet.from_cp2k,
+        "nwchem": BasisSet.from_nwchem,
+        "gaussian": BasisSet.from_gaussian,
     }
 
     filters = {
@@ -176,7 +179,7 @@ def list_basisset(sym, name, tags):
               help="filter by name")
 @click.option('tags', '--tag', '-t', multiple=True,
               help="filter by a tag (all tags must be present if specified multiple times)")
-@click.option('output_format', '-f', '--format', type=click.Choice(['cp2k', ]), default='cp2k',
+@click.option('output_format', '-f', '--format', type=click.Choice(['cp2k', 'nwchem', 'gamess',]), default='cp2k',
               help="Chose the output format for the basiset: " + ', '.join(['cp2k', ]))
 # fmt: on
 @decorators.with_dbenv()
@@ -191,6 +194,8 @@ def dump_basisset(sym, name, tags, output_format, data):
 
     writers = {
         "cp2k": BasisSet.to_cp2k,
+        "nwchem" : BasisSet.to_nwchem,
+        "gamess" : BasisSet.to_gamess,
     }
 
     if data:
diff --git a/aiida_gaussian_datatypes/basisset/data.py b/aiida_gaussian_datatypes/basisset/data.py
index 5135337..77f6a52 100644
--- a/aiida_gaussian_datatypes/basisset/data.py
+++ b/aiida_gaussian_datatypes/basisset/data.py
@@ -14,12 +14,15 @@
     NotExistent,
     UniquenessError,
     ValidationError,
+    ParsingError
 )
+import re
 from aiida.orm import Data, Group
+from pathlib import Path
 from cp2k_input_tools.basissets import BasisSetData
 
 
-class BasisSet(Data):
+class BasisSetCommon(Data):
     """
     Provide a general way to store GTO basis sets from different codes within the AiiDA framework.
     """
@@ -46,7 +49,7 @@ def __init__(self, element=None, name=None, aliases=None, tags=None, n_el=None,
         if "label" not in kwargs:
             kwargs["label"] = name
 
-        super(BasisSet, self).__init__(**kwargs)
+        super(BasisSetCommon, self).__init__(**kwargs)
 
         self.set_attribute("name", name)
         self.set_attribute("element", element)
@@ -57,38 +60,24 @@ def __init__(self, element=None, name=None, aliases=None, tags=None, n_el=None,
         self.set_attribute("version", version)
 
     def store(self, *args, **kwargs):
-        """
-        Store the node, ensuring that the combination (element,name,version) is unique.
-        """
-        # TODO: this uniqueness check is not race-condition free.
-
-        try:
-            existing = self.get(self.element, self.name, self.version, match_aliases=False)
-        except NotExistent:
-            pass
-        else:
-            raise UniquenessError(
-                f"Gaussian Basis Set already exists for"
-                f" element={self.element}, name={self.name}, version={self.version}: {existing.uuid}"
-            )
-
-        return super(BasisSet, self).store(*args, **kwargs)
+        return super(BasisSetCommon, self).store(*args, **kwargs)
 
     def _validate(self):
-        super(BasisSet, self)._validate()
+        super(BasisSetCommon, self)._validate()
 
         try:
             # directly raises an exception for the data if something's amiss, extra fields are ignored
+            # BasisSetData.from_dict({"identifiers": self.aliases, **self.attributes})
             _dict2basissetdata(self.attributes)
 
-            assert isinstance(self.name, str) and self.name
+            #assert isinstance(self.name, str) and self.name
             assert (
                 isinstance(self.aliases, list)
                 and all(isinstance(alias, str) for alias in self.aliases)
                 and self.aliases
             )
-            assert isinstance(self.tags, list) and all(isinstance(tag, str) for tag in self.tags)
-            assert isinstance(self.version, int) and self.version > 0
+            #assert isinstance(self.tags, list) and all(isinstance(tag, str) for tag in self.tags)
+            #assert isinstance(self.version, int) and self.version > 0
         except Exception as exc:
             raise ValidationError("One or more invalid fields found") from exc
 
@@ -291,6 +280,240 @@ def exists(bset):
 
         return [cls(**bs) for bs in bsets]
 
+    @classmethod
+    def from_gaussian(cls, fhandle, filters=None, duplicate_handling="ignore", attrs = None):
+        """
+        Constructs a list with basis set objects from a Basis Set in Gaussian format
+
+        :param fhandle: open file handle
+        :param filters: a dict with attribute filter functions
+        :param duplicate_handling: how to handle duplicates ("ignore", "error", "new" (version))
+        :rtype: list
+        """
+
+        def exists(bset):
+            try:
+                cls.get(bset["element"], bset["name"], match_aliases=False)
+            except NotExistent:
+                return False
+
+            return True
+
+        """
+        Gaussian parser
+
+        TODO Maybe parser should move to "parsers"
+        """
+
+        element = None
+        data = []
+        blocks = []
+
+        if not attrs:
+            attrs = {}
+
+        def block_creator(b, orb, blocks = blocks):
+            orb_dict = {"s" : 0,
+                        "p" : 1,
+                        "d" : 2,
+                        "f" : 3,
+                        "g" : 4,
+                        "h" : 5,
+                        "i" : 6 }
+            block = { "n": 0, # I dont know how to setup main quantum number
+                      "l": [(orb_dict[orb], 1)],
+                      "coefficients" : [ [ d["exp"], d["cont"] ] for d in b ] }
+            blocks.append(block)
+
+        orb = "x"
+        for ii, line in enumerate(fhandle):
+            if ii == 1:
+                element = line.lower().split()[0]
+                continue
+            if re.match(r"^[A-z ]+[0-9\. ]*$", line):
+                if len(data) != 0:
+                    block_creator(data, orb)
+                data = []
+                orb = line.lower().split()[0]
+            if re.match(r"^[+-.0-9 ]+$", line):
+                exp, cont, = [ float(x) for x in line.split() ]
+                data.append({"exp" : exp,
+                             "cont" : cont })
+        if len(data) != 0:
+            block_creator(data, orb)
+            data = []
+
+        try:
+            basis = {"element" : element.capitalize(),
+                     "version" : 1,
+                     "name" : "unknown",
+                     "tags" : [],
+                     "aliases" : [],
+                     "blocks" : blocks }
+        except:
+            return []
+
+        basis["name"] = "NA"
+
+        if hasattr(fhandle, "name"):
+            basis["name"] = Path(fhandle.name).name.replace(".nwchem", "")
+            basis["aliases"].append(basis["name"].split(".")[-1])
+
+        if "name" in attrs:
+            basis["aliases"].append(basis["name"])
+            basis["name"] = attrs["name"]
+
+        for attr in ("n_el", "tags",):
+            if attr in attrs:
+                basis[attr] = attrs[attr]
+
+        if len(basis["aliases"]) == 0:
+            del basis["aliases"]
+
+        if duplicate_handling == "force-ignore":  # It will check at the store stage
+            pass
+
+        elif duplicate_handling == "ignore":  # simply filter duplicates
+            if exists(basis):
+                return []
+
+        elif duplicate_handling == "error":
+            if exists(basis):
+                raise UniquenessError( f"Gaussian Basis Set already exists for"
+                                       f" element={basis['element']}, name={basis['name']}: {latest.uuid}")
+
+        elif duplicate_handling == "new":
+                try:
+                    latest = cls.get(basis["element"], basis["name"], match_aliases=False)
+                except NotExistent:
+                    pass
+                else:
+                    basis["version"] = latest.version + 1
+
+        else:
+            raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
+
+        return [cls(**basis)]
+
+    @classmethod
+    def from_nwchem(cls, fhandle, filters=None, duplicate_handling="ignore", attrs = None):
+        """
+        Constructs a list with basis set objects from a Basis Set in NWCHEM format
+
+        :param fhandle: open file handle
+        :param filters: a dict with attribute filter functions
+        :param duplicate_handling: how to handle duplicates ("ignore", "error", "new" (version))
+        :rtype: list
+        """
+
+        def exists(bset):
+            try:
+                cls.get(bset["element"], bset["name"], match_aliases=False)
+            except NotExistent:
+                return False
+
+            return True
+
+        """
+        NWCHEM parser
+
+        TODO Maybe parser should move to "parsers"
+        """
+
+        element = None
+        data = []
+        blocks = []
+
+        if not attrs:
+            attrs = {}
+
+        def block_creator(b, orb, blocks = blocks):
+            orb_dict = {"s" : 0,
+                        "p" : 1,
+                        "d" : 2,
+                        "f" : 3,
+                        "g" : 4,
+                        "h" : 5,
+                        "i" : 6 }
+            block = { "n": 0, # I dont know how to setup main quantum number
+                      "l": [(orb_dict[orb], 1)],
+                      "coefficients" : [ [ d["exp"], d["cont"] ] for d in b ] }
+            blocks.append(block)
+
+        for line in fhandle:
+            """
+            Element symbol has to be every block
+            """
+            if re.match("^[A-z ]+$", line):
+                if len(data) != 0:
+                    block_creator(data, orb)
+                    data = []
+                el, orb = line.lower().split()
+                if element is None:
+                    """
+                    TODO check validity of element
+                    """
+                    element = el
+                elif element != el:
+                    raise ParsingError(f"Element previous {element}, and now {el}.") # Element cannot be changed
+            if re.match("^[+-.0-9 ]+$", line):
+                exp, cont, = [ float(x) for x in line.split() ]
+                data.append({"exp" : exp,
+                             "cont" : cont })
+        if len(data) != 0:
+            block_creator(data, orb)
+            data = []
+
+        try:
+            basis = {"element" : element.capitalize(),
+                     "version" : 1,
+                     "name" : "unknown",
+                     "tags" : [],
+                     "aliases" : [],
+                     "blocks" : blocks }
+        except:
+            return []
+
+        if hasattr(fhandle, "name"):
+            basis["name"] = Path(fhandle.name).name.replace(".nwchem", "")
+            basis["aliases"].append(basis["name"].split(".")[-1])
+
+        if "name" in attrs:
+            basis["aliases"].append(basis["name"])
+            basis["name"] = attrs["name"]
+
+        for attr in ("n_el", "tags",):
+            if attr in attrs:
+                basis[attr] = attrs[attr]
+
+        if len(basis["aliases"]) == 0:
+            del basis["aliases"]
+
+        if duplicate_handling == "force-ignore":  # It will check at the store stage
+            pass
+
+        elif duplicate_handling == "ignore":  # simply filter duplicates
+            if exists(basis):
+                return []
+
+        elif duplicate_handling == "error":
+            if exists(basis):
+                raise UniquenessError( f"Gaussian Basis Set already exists for"
+                                       f" element={basis['element']}, name={basis['name']}: {latest.uuid}")
+
+        elif duplicate_handling == "new":
+                try:
+                    latest = cls.get(basis["element"], basis["name"], match_aliases=False)
+                except NotExistent:
+                    pass
+                else:
+                    basis["version"] = latest.version + 1
+
+        else:
+            raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
+
+        return [cls(**basis)]
+
     def to_cp2k(self, fhandle):
         """
         Write the Basis Set to the passed file handle in the format expected by CP2K.
@@ -302,6 +525,85 @@ def to_cp2k(self, fhandle):
             fhandle.write(line)
             fhandle.write("\n")
 
+    def to_nwchem(self, fhandle):
+        """
+        Write the Basis Set to the passed file handle in the format expected by NWCHEM.
+
+        :param fhandle: A valid output file handle
+        """
+        orb_dict = {0 : "s",
+                    1 : "p",
+                    2 : "d",
+                    3 : "f",
+                    4 : "g",
+                    5 : "h",
+                    6 : "i" }
+
+        fhandle.write(f"# from AiiDA BasisSet<uuid: {self.uuid}>\n")
+        for block in self.blocks:
+            offset = 0
+            for orb, num, in block["l"]:
+                fhandle.write(f"{self.element} {orb_dict[orb]}\n")
+                for lnum in range(num):
+                    for entry in block["coefficients"]:
+                        exponent = entry[0]
+                        coefficient = entry[1 + lnum + offset]
+                        fhandle.write(f"  {exponent:15.7f} {coefficient:15.7f}\n")
+                offset = num
+
+    def to_gamess(self, fhandle):
+        """
+        Write the Basis Set to the passed file handle in the format expected by GAMESS.
+
+        :param fhandle: A valid output file handle
+        """
+        orb_dict = {0 : "s",
+                    1 : "p",
+                    2 : "d",
+                    3 : "f",
+                    4 : "g",
+                    5 : "h",
+                    6 : "i" }
+
+        fhandle.write(f"# from AiiDA BasisSet<uuid: {self.uuid}>\n")
+        for block in self.blocks:
+            offset = 0
+            for orb, num, in block["l"]:
+                fhandle.write(f" {orb_dict[orb].upper()}  {len(block['coefficients'])}\n")
+                for lnum in range(num):
+                    for ii, entry in enumerate(block["coefficients"]):
+                        exponent = entry[0]
+                        coefficient = entry[1 + lnum + offset]
+                        fhandle.write(f"  {ii + 1:3d} {exponent:15.7f} {coefficient:15.7f}\n")
+                offset = num
+
+    def to_gaussian(self, fhandle):
+        """
+        Write the Basis Set to the passed file handle in the format expected by Gaussian.
+
+        :param fhandle: A valid output file handle
+        """
+        orb_dict = {0 : "s",
+                    1 : "p",
+                    2 : "d",
+                    3 : "f",
+                    4 : "g",
+                    5 : "h",
+                    6 : "i" }
+
+        fhandle.write(f"# from AiiDA BasisSet<uuid: {self.uuid}>\n")
+        for block in self.blocks:
+            offset = 0
+            for orb, num, in block["l"]:
+                fhandle.write(f" {orb_dict[orb].upper()}  {len(block['coefficients'])}\n")
+                for lnum in range(num):
+                    for ii, entry in enumerate(block["coefficients"]):
+                        exponent = entry[0]
+                        coefficient = entry[1 + lnum + offset]
+                        fhandle.write(f"  {ii + 1:3d} {exponent:15.7f} {coefficient:15.7f}\n")
+                offset = num
+
+
     def get_matching_pseudopotential(self, *args, **kwargs):
         """
         Get a pseudopotential matching this basis set by at least element and number of valence electrons.
@@ -314,6 +616,36 @@ def get_matching_pseudopotential(self, *args, **kwargs):
         else:
             return Pseudopotential.get(element=self.element, *args, **kwargs)
 
+class BasisSet(BasisSetCommon):
+
+    def __init__(self, *args, **kwargs):
+        super(BasisSet, self).__init__(*args, **kwargs)
+
+    def store(self, *args, **kwargs):
+        """
+        Store the node, ensuring that the combination (element,name,version) is unique.
+        """
+        # TODO: this uniqueness check is not race-condition free.
+
+        try:
+            existing = self.get(self.element, self.name, self.version, match_aliases=False)
+        except NotExistent:
+            pass
+        else:
+            raise UniquenessError(
+                f"Gaussian Basis Set already exists for"
+                f" element={self.element}, name={self.name}, version={self.version}: {existing.uuid}"
+            )
+
+        return super(BasisSet, self).store(*args, **kwargs)
+
+class BasisSetFree(BasisSetCommon):
+
+    def __init__(self, *args, **kwargs):
+        super(BasisSetFree, self).__init__(*args, **kwargs)
+
+    def store(self, *args, **kwargs):
+        return super(BasisSetFree, self).store(*args, **kwargs)
 
 def _basissetdata2dict(data: BasisSetData) -> Dict[str, Any]:
     """
diff --git a/aiida_gaussian_datatypes/calc/__init__.py b/aiida_gaussian_datatypes/calc/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/aiida_gaussian_datatypes/calc/__init__py b/aiida_gaussian_datatypes/calc/__init__py
new file mode 100644
index 0000000..e69de29
diff --git a/aiida_gaussian_datatypes/calc/uncontract.py b/aiida_gaussian_datatypes/calc/uncontract.py
new file mode 100644
index 0000000..6bc4654
--- /dev/null
+++ b/aiida_gaussian_datatypes/calc/uncontract.py
@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+
+from aiida.plugins import DataFactory
+from aiida.engine import calcfunction
+"""
+
+"""
+
+BasisSet = DataFactory("gaussian.basisset")
+BasisSetFree = DataFactory("gaussian.basissetfree")
+
+@calcfunction
+def uncontract(basisset):
+    """
+
+    """
+    def disassemble(block):
+        n = block["n"]
+        l = block["l"]
+        for exp, cont in block["coefficients"]:
+            yield {"n" : n,
+                   "l" : l,
+                   "coefficients": [[exp, 1.0]]}
+    attr = basisset.attributes
+    blocks = []
+    for block in attr["blocks"]:
+        blocks.extend([ b for b in disassemble(block) ])
+    attr["blocks"] = blocks
+    attr["name"] += "-uncont"
+    ret = BasisSetFree(**attr)
+    return ret
+
diff --git a/aiida_gaussian_datatypes/fetcher/__init__.py b/aiida_gaussian_datatypes/fetcher/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/aiida_gaussian_datatypes/fetcher/cli.py b/aiida_gaussian_datatypes/fetcher/cli.py
new file mode 100644
index 0000000..37fc395
--- /dev/null
+++ b/aiida_gaussian_datatypes/fetcher/cli.py
@@ -0,0 +1,179 @@
+# -*- coding: utf-8 -*-
+
+import click
+import tabulate
+import pydriller
+from pathlib import Path
+from aiida.cmdline.utils import decorators, echo
+from aiida.cmdline.commands.cmd_data import verdi_data
+from aiida.orm import load_group
+from ..libraries import *
+from ..basisset.data import BasisSet
+from ..pseudopotential.data import Pseudopotential
+#from ..groups import (
+#    BasisSetGroup,
+#    PseudopotenialGroup,
+#)
+from ..groups import BasisSetGroup
+from ..groups import PseudopotentialGroup
+
+#from ..utils import (
+#    click_parse_range,  # pylint: disable=relative-beyond-top-level
+#    SYM2NUM,
+#)
+
+from ..utils import click_parse_range
+from ..utils import SYM2NUM
+from aiida.common.exceptions import UniquenessError
+
+def _formatted_table_import(elements):
+    """generates a formatted table (using tabulate) for importable basis and PPs"""
+
+    def _boldformater(f):
+
+        def fout(*args, **kwargs):
+            if args[1] % 2 == 1:
+                return ( f"\033[1m{x}\033[0m" for x in f(*args, **kwargs))
+            else:
+                return ( x for x in f(*args, **kwargs))
+        return fout
+
+    class row():
+
+        num = []
+        element = []
+        t = []
+
+        @_boldformater
+        def __new__(cls, num, element, t, p, tags, b):
+
+            if element in cls.element:
+                element = ""
+            else:
+                cls.element.append(element)
+                element = str(element)
+                cls.t = []
+
+            if num in cls.num:
+                num = ""
+            else:
+                cls.num.append(num)
+                num = str(num)
+
+            if t in cls.t:
+                t = ""
+            else:
+                cls.t.append(t)
+
+            if t == "":
+                p = ""
+                tags = []
+
+            name = ""
+            m = re.match("http:\/\/burkatzki\.com\|([A-z]+)", b)
+            if m:
+                name = m.group(1)
+            m = re.match("[A-z]{1,2}\.(.+).nwchem", b)
+            if m:
+                name = m.group(1)
+
+            return (
+                num,
+                element,
+                t,
+                p,
+                " ".join(sorted(tags)),
+                name,
+                b
+            )
+
+    table_content = []
+    for ii, (e, d) in enumerate(elements):
+        for t in d["types"]:
+            if len(d["types"][t]["pseudos"]) == 0:
+                continue
+            p = d["types"][t]["pseudos"][0]
+            for b in d["types"][t]["basis"]:
+                name = ""
+                if isinstance(b["path"], str):
+                    name = b["path"]
+                if hasattr(b["path"], "name"):
+                    name = b["path"].name
+                table_content.append(row(ii, e, t, name,
+                                         d["types"][t]["tags"],
+                                         name))
+
+    #table_content = [row(n, p, v) for n, (p, v) in enumerate(elements.items())]
+    return tabulate.tabulate(table_content, headers=["Nr.", "Element", "Type", "PseudoFile", "Tags", "Basis", "BasisFile"])
+
+@verdi_data.group("gaussian")
+def cli():
+    """Manage Pseudopotentials for GTO-based codes"""
+
+# fmt: off
+@cli.command('fetch')
+@click.argument('library',
+                type=click.Choice(LibraryBookKeeper.get_library_names()))
+@decorators.with_dbenv()
+# fmt: on
+def install_family(library):
+    """
+    Installs a family of pseudo potentials from a remote repository
+    """
+
+    basissetgname = f"{library}-basis"
+    try:
+        basisgroup = load_group(basissetgname)
+    except:
+        echo.echo_info("Creating library basis group ... ", nl = False)
+        basisgroup = BasisSetGroup(basissetgname)
+        basisgroup.store()
+        echo.echo("DONE")
+
+    pseudogname = f"{library}-pseudo"
+    try:
+        pseudogroup = load_group(pseudogname)
+    except:
+        echo.echo_info("Creating library pseudo group ... ", nl = False)
+        pseudogroup = PseudopotentialGroup(pseudogname)
+        pseudogroup.store()
+        echo.echo("DONE")
+
+    elements = LibraryBookKeeper.get_library_by_name(library).fetch()
+
+    elements = [ [el, p] for el, p in sorted(elements.items(), key = lambda x: SYM2NUM[x[0]]) ]
+    echo.echo_info(f"Found {len(elements)} elements")
+    echo.echo(_formatted_table_import(elements))
+    echo.echo("")
+    indexes = click.prompt(
+        "Which Elements do you want to add?"
+        " ('n' for none, 'a' for all, comma-seperated list or range of numbers)",
+        value_proc=lambda v: click_parse_range(v, len(elements)))
+    for idx in indexes:
+        e, v = elements[idx]
+        for t, o in v["types"].items():
+            for b in o["basis"]:
+                basis = b["obj"]
+                echo.echo_info(f"Adding Basis for: ", nl=False)
+                echo.echo(f"{basis.element} ({basis.name})...  ", nl=False)
+                try:
+                    basis.store()
+                    basisgroup.add_nodes([basis])
+                    echo.echo("Imported")
+                except UniquenessError:
+                    echo.echo("Skipping (already in)")
+                except Exception as e:
+                    echo.echo("Skipping (something went wrong)")
+            for p in o["pseudos"]:
+                pseudo = p["obj"]
+                echo.echo_info(f"Adding Pseudopotential for: ", nl=False)
+                echo.echo(f"{pseudo.element} ({pseudo.name})...  ", nl=False)
+                try:
+                    pseudo.store()
+                    pseudogroup.add_nodes([pseudo])
+                    echo.echo("Imported")
+                except UniquenessError:
+                    echo.echo("Skipping (already in)")
+                except Exception as e:
+                    echo.echo("Skipping (something went wrong)")
+
diff --git a/aiida_gaussian_datatypes/libraries.py b/aiida_gaussian_datatypes/libraries.py
new file mode 100644
index 0000000..4e58933
--- /dev/null
+++ b/aiida_gaussian_datatypes/libraries.py
@@ -0,0 +1,220 @@
+# -*- coding: utf-8 -*-
+# SPDX-License-Identifier: MIT
+#
+# Was there really a fish
+# That grants you that kind of wish
+#
+
+import os
+import re
+import git
+import tempfile
+import pathlib
+import pydriller
+from aiida_gaussian_datatypes import utils
+from typing import Dict, Generic, List, Optional, Sequence, Type, TypeVar
+from .basisset.data import BasisSet
+from .pseudopotential.data import Pseudopotential, ECPPseudopotential
+
+from aiida.common.exceptions import (
+    NotExistent,
+)
+
+class LibraryBookKeeper:
+
+    classes = []
+
+    @classmethod
+    def register_library(cls, cls_):
+        cls.classes.append(cls_)
+
+    @classmethod
+    def get_libraries(cls):
+        return cls.classes
+
+    @classmethod
+    def get_library_names(cls):
+        return [ re.match("<class '[0-9A-z_\.]*\.([A-z]+)'>", str(x)).group(1) for x in cls.classes ]
+
+    @classmethod
+    def get_library_by_name(cls, name):
+        for cls_ in cls.get_libraries():
+            if re.match(f"<class '[0-9A-z_\.]*\.({name})'>", str(cls_)) is not None:
+                return cls_
+        return None
+
+class _ExternalLibrary:
+
+    @classmethod
+    def fetch(cls):
+        pass
+
+@LibraryBookKeeper.register_library
+class EmptyLibrary(_ExternalLibrary):
+    pass
+
+@LibraryBookKeeper.register_library
+class QmcpackLibrary(_ExternalLibrary):
+
+    _URL = "https://github.com/QMCPACK/pseudopotentiallibrary.git"
+
+    @classmethod
+    def fetch(cls):
+
+        elements = {}
+        def add_data(p, tempdir, elements = elements):
+            element = str(p.parent.parent.name)
+            if element not in utils.SYM2NUM: # Check if element is valid
+                return
+            element_path = p.parent.parent
+
+            typ = str(p.parent.name)
+            typ_path = str(p.parent.name)
+
+            tags = ["ECP", typ, ]
+
+            """ Load Pseudopotential first """
+            with open(p, "r") as fhandle:
+                pseudo, = Pseudopotential.from_gamess(fhandle,
+                                                      duplicate_handling = "force-ignore",
+                                                      attrs = {"name" : typ })
+            commithash = ""
+            for version, commit in enumerate(pydriller.Repository(str(tempdir), filepath=str(p)).traverse_commits()):
+                commithash = commit.hash
+            if commithash == "": return
+            pseudo.extras["commithash"] = commithash
+            pseudo.attributes["version"] = version + 1
+
+            tags.append(f"q{pseudo.n_el_tot}")
+            tags.append(f"c{pseudo.core_electrons}")
+            pseudo.tags.extend(tags)
+
+            pseudos = [{"path": p,
+                        "obj": pseudo}]
+
+            """ Load Basis sets """
+            basis = []
+            for r in (p.parent).glob("**/*"):
+                if re.match("[A-z]{1,2}\.[A-z\-]*cc-.*\.nwchem", r.name):
+                    name = re.match("[A-z]{1,2}\.([A-z\-]*cc-.*)\.nwchem", r.name).group(1)
+                    name = f"{typ}-{name}"
+                    with open(r, "r") as fhandle:
+                        b = BasisSet.from_nwchem(fhandle,
+                                                 duplicate_handling = "force-ignore",
+                                                 attrs = {"n_el": pseudo.n_el_tot,
+                                                          "name": name,
+                                                          "tags": tags})
+                        if len(b) == 0: continue
+                        b, = b
+
+                        commithash = ""
+                        for version, commit in enumerate(pydriller.Repository(str(tempdir), filepath=str(r)).traverse_commits()):
+                            commithash = commit.hash
+                        if commithash == "": return
+                        b.extras["commithash"] = commithash
+                        b.attributes["version"] = version + 1
+
+                    basis.append({"path": r,
+                                  "obj": b})
+
+
+            if element not in elements:
+                elements[element] = {"path": element_path,
+                                     "types": {}}
+
+            if typ not in elements[element]["types"]:
+                elements[element]["types"][typ] = {"path": typ_path,
+                                                   "basis": basis,
+                                                   "pseudos": pseudos,
+                                                   "tags": tags}
+
+        tempdir = pathlib.Path(tempfile.mkdtemp())
+        git.Repo.clone_from(cls._URL, tempdir)
+
+        for p in (tempdir/"recipes").glob("**/*"):
+            if re.match("[A-z]{1,2}\.ccECP\.gamess", p.name):
+                add_data(p, tempdir)
+
+        return elements
+
+@LibraryBookKeeper.register_library
+class BFDLibrary(_ExternalLibrary):
+
+    _URL = "http://burkatzki.com/pseudos/step4.2.php?format=gaussian&element={e}&basis={b}"
+
+    @classmethod
+    def fetch(cls):
+
+        from ase.data import chemical_symbols
+        from ase.data import atomic_numbers
+        from urllib.request import urlopen
+        from time import sleep
+        import io
+
+        elements = {}
+        def add_data(source, e, b):
+
+            source = str(source)
+            pat=re.compile("^.*?(" + e + "\s0.*$)",re.M|re.DOTALL)
+            x = pat.sub("\g<1>", source)
+            x = re.sub("\<br\s*/\>", "\n", x)
+            x = re.sub("\&nbsp", "", x)
+            x = re.sub("\&nbsp", "", x)
+            x = re.sub(".*html.*$", "", x)
+            pat=re.compile("^.*?(" + e + "\s0.*)("+e+" 0.*)$",re.M|re.DOTALL)
+            m = pat.match(x)
+            if(m):
+                bas = m.group(1)
+                ecp = m.group(2)
+                if len(bas) < 15: return
+                typ = "BFD"
+                pseudo, = Pseudopotential.from_gaussian(io.StringIO(ecp),
+                                                        duplicate_handling = "force-ignore",
+                                                        attrs = {"name" : f"{typ}" })
+
+                tags = [typ]
+                if "ano" in b:
+                    tags.append("ANO")
+                basisset, = BasisSet.from_gaussian(io.StringIO(f"\n{bas}"),
+                                                   duplicate_handling = "force-ignore",
+                                                   attrs = {"name" : f"{typ}-{b}",
+                                                            "n_el" : pseudo.n_el_tot,
+                                                            "tags" : tags})
+                pseudos = [{"path": "",
+                            "obj": pseudo}]
+                version = 1
+                pseudo.attributes["version"] = version
+
+                tags = []
+                tags.append(f"q{pseudo.n_el_tot}")
+                tags.append(f"c{pseudo.core_electrons}")
+                pseudo.tags.extend(tags)
+
+                if e not in elements:
+                    elements[e] = {"path": "",
+                                   "types": {}}
+
+                if typ not in elements[e]["types"]:
+                    elements[e]["types"][typ] = {"path": "",
+                                                 "basis": [],
+                                                 "pseudos": pseudos,
+                                                 "tags": []}
+                elements[e]["types"][typ]["tags"].extend(tags)
+                elements[e]["types"][typ]["basis"].append({"path" : f"http://burkatzki.com|{b}",
+                                                           "obj"  : basisset})
+                elements[e]["types"][typ]["tags"].append("BFD")
+                tt = set(elements[e]["types"][typ]["tags"])
+                elements[e]["types"][typ]["tags"] = list(tt)
+
+
+        list_of_basis  =[ f"v{s}z" for s in "dtq56" ]
+        list_of_basis += [ f"{x}_ano" for x in list_of_basis ]
+
+        for b in list_of_basis:
+            for ie in range(1, 86):
+                l = cls._URL.format(b = b, e = chemical_symbols[ie])
+                add_data(urlopen(l).read(), chemical_symbols[ie], b)
+                """ Cool down """
+                sleep(0.5)
+
+        return elements
diff --git a/aiida_gaussian_datatypes/pseudopotential/cli.py b/aiida_gaussian_datatypes/pseudopotential/cli.py
index f2b9ba1..5da6678 100644
--- a/aiida_gaussian_datatypes/pseudopotential/cli.py
+++ b/aiida_gaussian_datatypes/pseudopotential/cli.py
@@ -31,15 +31,17 @@ def _formatted_table_import(pseudos):
     def row(num, pseudo):
         return (
             num + 1,
+            pseudo.__name__.replace("Pseudopotential", "") if hasattr(pseudo, "__name__") else "",
             pseudo.element,
             _names_column(pseudo.name, pseudo.aliases),
             ", ".join(pseudo.tags),
             ", ".join(f"{n:2d}" for n in pseudo.n_el + (3 - len(pseudo.n_el)) * [0]),
+            pseudo.n_el_tot,
             pseudo.version,
         )
 
     table_content = [row(n, p) for n, p in enumerate(pseudos)]
-    return tabulate.tabulate(table_content, headers=["Nr.", "Sym", "Names", "Tags", "Val. e⁻ (s, p, ..)", "Version"])
+    return tabulate.tabulate(table_content, headers=["Nr.", "Type", "Sym", "Names", "Tags", "Val. e⁻ (s, p, ..)", "Tot. val. e⁻", "Version"])
 
 
 def _formatted_table_list(pseudos):
@@ -48,15 +50,17 @@ def _formatted_table_list(pseudos):
     def row(pseudo):
         return (
             pseudo.uuid,
+            pseudo.__name__.replace("Pseudopotential", "") if hasattr(pseudo, "__name__") else "",
             pseudo.element,
             _names_column(pseudo.name, pseudo.aliases),
             ", ".join(pseudo.tags),
             ", ".join(f"{n:2d}" for n in pseudo.n_el + (3 - len(pseudo.n_el)) * [0]),
+            pseudo.n_el_tot,
             pseudo.version,
         )
 
     table_content = [row(p) for p in pseudos]
-    return tabulate.tabulate(table_content, headers=["ID", "Sym", "Names", "Tags", "Val. e⁻ (s, p, ..)", "Version"])
+    return tabulate.tabulate(table_content, headers=["Nr.", "Type", "Sym", "Names", "Tags", "Val. e⁻ (s, p, ..)", "Tot. val. e⁻", "Version"])
 
 
 @verdi_data.group("gaussian.pseudo")
@@ -74,7 +78,7 @@ def cli():
     help="filter by a tag (all tags must be present if specified multiple times)")
 @click.option(
     'fformat', '-f', '--format',
-    type=click.Choice(['cp2k', ]), default='cp2k',
+    type=click.Choice(['cp2k', 'gamess', 'turborvb' ]), default='cp2k',
     help="the format of the pseudopotential file")
 @click.option(
     '--duplicates',
@@ -97,6 +101,8 @@ def import_pseudo(pseudopotential_file, fformat, sym, tags, duplicates, ignore_i
 
     loaders = {
         "cp2k": Pseudopotential.from_cp2k,
+        "gamess": Pseudopotential.from_gamess,
+        "turborvb": Pseudopotential.from_turborvb,
     }
 
     filters = {
@@ -181,11 +187,15 @@ def list_pseudo(sym, name, tags):
               help="filter by name")
 @click.option('tags', '--tag', '-t', multiple=True,
               help="filter by a tag (all tags must be present if specified multiple times)")
-@click.option('output_format', '-f', '--format', type=click.Choice(['cp2k', ]), default='cp2k',
+@click.option('output_format', '-f', '--format', type=click.Choice(['cp2k',
+                                                                    'gamess',
+                                                                    'turborvb']), default='cp2k',
               help="Chose the output format for the pseudopotentials: " + ', '.join(['cp2k', ]))
+@click.option('-r', '--tolerance', type=float, default=1.0e-5,
+              help="set tolerance value for pseudo cutoff (default 1.0e-5, only for turborvb format)")
 @decorators.with_dbenv()
 # fmt: on
-def dump_pseudo(sym, name, tags, output_format, data):
+def dump_pseudo(sym, name, tags, output_format, data, tolerance):
     """
     Print specified Pseudopotentials
     """
@@ -196,6 +206,8 @@ def dump_pseudo(sym, name, tags, output_format, data):
 
     writers = {
         "cp2k": Pseudopotential.to_cp2k,
+        "gamess": Pseudopotential.to_gamess,
+        "turborvb": Pseudopotential.to_turborvb,
     }
 
     if data:
@@ -225,4 +237,4 @@ def dump_pseudo(sym, name, tags, output_format, data):
         if echo.is_stdout_redirected():
             echo.echo_report("Dumping {}/{} ({})...".format(pseudo.name, pseudo.element, pseudo.uuid), err=True)
 
-        writers[output_format](pseudo, sys.stdout)
+        writers[output_format](pseudo, sys.stdout, tolerance = tolerance)
diff --git a/aiida_gaussian_datatypes/pseudopotential/data.py b/aiida_gaussian_datatypes/pseudopotential/data.py
index c1fb524..11de26f 100644
--- a/aiida_gaussian_datatypes/pseudopotential/data.py
+++ b/aiida_gaussian_datatypes/pseudopotential/data.py
@@ -6,6 +6,11 @@
 Gaussian Pseudopotential Data class
 """
 
+import dataclasses
+from ..utils import SYM2NUM
+from decimal import Decimal
+import re
+import numpy as np
 from decimal import Decimal
 from typing import Any, Dict
 
@@ -25,6 +30,8 @@ class Pseudopotential(Data):
     fixme: extend to NLCC pseudos.
     """
 
+    __name__ = "Pseudopotential"
+
     def __init__(
         self,
         element=None,
@@ -32,9 +39,7 @@ def __init__(
         aliases=None,
         tags=None,
         n_el=None,
-        local=None,
-        non_local=None,
-        nlcc=None,
+        n_el_tot=None,
         version=1,
         **kwargs,
     ):
@@ -44,8 +49,6 @@ def __init__(
         :param aliases: alternative names
         :param tags: additional tags
         :param n_el: number of valence electrons covered by this basis set
-        :param local: see :py:attr:`~local`
-        :param local: see :py:attr:`~non_local`
         """
 
         if not aliases:
@@ -56,19 +59,18 @@ def __init__(
 
         if not n_el:
             n_el = []
-
-        if not non_local:
-            non_local = []
-
-        if not nlcc:
-            nlcc = []
+        else:
+            if not n_el_tot:
+                n_el_tot = sum(n_el)
+            else:
+                raise #TODO a propiate error here
 
         if "label" not in kwargs:
             kwargs["label"] = name
 
         super().__init__(**kwargs)
 
-        for attr in ("name", "element", "tags", "aliases", "n_el", "local", "non_local", "nlcc", "version"):
+        for attr in ("name", "element", "tags", "aliases", "n_el", "n_el_tot", "version"):
             self.set_attribute(attr, locals()[attr])
 
     def store(self, *args, **kwargs):
@@ -94,7 +96,6 @@ def _validate(self):
 
         try:
             # directly raises a ValidationError for the pseudo data if something's amiss
-            _dict2pseudodata(self.attributes)
 
             assert isinstance(self.name, str) and self.name
             assert (
@@ -104,6 +105,8 @@ def _validate(self):
             )
             assert isinstance(self.tags, list) and all(isinstance(tag, str) for tag in self.tags)
             assert isinstance(self.version, int) and self.version > 0
+            if len(self.n_el) != 0:
+                assert(sum(self.n_el) == self.n_el_tot)
         except Exception as exc:
             raise ValidationError("One or more invalid fields found") from exc
 
@@ -162,46 +165,13 @@ def n_el(self):
         return self.get_attribute("n_el", [])
 
     @property
-    def local(self):
-        """
-        Return the local part
-
-        The format of the returned dictionary::
-
-            {
-                'r': float,
-                'coeffs': [float, float, ...],
-            }
-
-        :rtype:dict
-        """
-        return self.get_attribute("local", None)
-
-    @property
-    def non_local(self):
+    def n_el_tot(self):
         """
-        Return a list of non-local projectors (for l=0,1...).
-
-        Each list element will have the following format::
-
-            {
-                'r': float,
-                'nproj': int,
-                'coeffs': [float, float, ...],  # only the upper-triangular elements
-            }
-
-        :rtype:list
-        """
-        return self.get_attribute("non_local", [])
-
-    @property
-    def nlcc(self):
+        Return the number of electrons per angular momentum
+        :rtype:int
         """
-        Return a list of the non-local core-corrections data
 
-        :rtype:list
-        """
-        return self.get_attribute("nlcc", [])
+        return self.get_attribute("n_el_tot", [])
 
     @classmethod
     def get(cls, element, name=None, version="latest", match_aliases=True, group_label=None, n_el=None):
@@ -223,7 +193,7 @@ def get(cls, element, name=None, version="latest", match_aliases=True, group_lab
             query.append(Group, filters={"label": group_label}, tag="group")
             params["with_group"] = "group"
 
-        query.append(Pseudopotential, **params)
+        query.append(cls, **params)
 
         filters = {"attributes.element": {"==": element}}
 
@@ -236,7 +206,7 @@ def get(cls, element, name=None, version="latest", match_aliases=True, group_lab
             else:
                 filters["attributes.name"] = {"==": name}
 
-        query.add_filter(Pseudopotential, filters)
+        query.add_filter(cls, filters)
 
         # SQLA ORM only solution:
         # query.order_by({Pseudopotential: [{"attributes.version": {"cast": "i", "order": "desc"}}]})
@@ -324,19 +294,420 @@ def exists(pseudo):
         else:
             raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
 
-        return [cls(**p) for p in pseudos]
+        return [GTHPseudopotential(**p) for p in pseudos]
+
+    @classmethod
+    def from_gaussian(cls, fhandle, filters=None, duplicate_handling="ignore", ignore_invalid=False, attrs = None):
+        """
+        Constructs a list with pseudopotential objects from a Pseudopotential in Gaussian format
+
+        :param fhandle: open file handle
+        :param filters: a dict with attribute filter functions
+        :param duplicate_handling: how to handle duplicates ("ignore", "error", "new" (version))
+        :param ignore_invalid: whether to ignore invalid entries silently
+        :rtype: list
+        """
+
+        def exists(pseudo):
+            try:
+                cls.get(pseudo["element"], pseudo["name"], match_aliases=False)
+            except NotExistent:
+                return False
+
+            return True
+
+        if not attrs:
+            attrs = {}
+
+        """
+        Parser for Gaussian format
+        """
+
+        block_counter = 0
+        functions = []
+        for ii, line in enumerate(fhandle):
+            ic(line.strip())
+            if ii == 0:
+                element, n, = line.split()
+                continue
+            if ii == 1:
+                qmc, n, core_electrons, = line.split()
+                continue
+            if ii == 2:
+                continue
+            else:
+                ic(block_counter)
+                if block_counter == 0:
+                    if line.strip() == "":
+                        continue
+                    m = re.match("[ ]*([0-9])+[ ]*$", line)
+                    if m:
+                        block_counter = int(m.group(1))
+                        functions.append({"prefactors" : [],
+                                          "polynoms"   : [],
+                                          "exponents"  : []})
+                else:
+                    functions[-1]["polynoms"].append(int(line.strip().split()[0]))
+                    functions[-1]["exponents"].append(float(line.strip().split()[1]))
+                    functions[-1]["prefactors"].append(float(line.strip().split()[2]))
+                    block_counter -= 1
+
+        """
+        Change the order of functions so they match orbital momentum
+
+        In Gaussian format first block represents upper most lmax
+        and then the rest s, p, d, ...
+        """
+        functions = functions[1:] + [functions[0]]
+
+        """
+        TODO properly extract name
+        """
+
+        lmax = len(functions) - 1
+        core_electrons = int(core_electrons)
+
+        data = {"functions"      : functions,
+                "element"        : element,
+                "aliases"        : [qmc],
+                "name"           : qmc,
+                "core_electrons" : core_electrons,
+                "lmax"           : lmax,
+                "version"        : 1,
+                "n_el"           : None,
+                "n_el_tot"       : SYM2NUM[element] - core_electrons}
+
+        if "name" in attrs:
+            data["aliases"].append(data["name"])
+            data["name"] = attrs["name"]
+
+        if duplicate_handling == "force-ignore":  # This will be checked at the store stage
+            pass
+
+        elif duplicate_handling == "ignore":  # simply filter duplicates
+            if exists(data):
+                return []
+
+        elif duplicate_handling == "error":
+            if exists(data):
+                raise UniquenessError(
+                    f"Gaussian Pseudopotential already exists for"
+                    f" element={data['element']}, name={data['name']}: {latest.uuid}"
+                )
+
+        elif duplicate_handling == "new":
+            if exists(data):
+                latest = cls.get(data["element"], data["name"], match_aliases=False)
+                data["version"] = latest.version + 1
+
+        else:
+            raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
+
+        return [ECPPseudopotential(**data)]
+
+    @classmethod
+    def from_gamess(cls, fhandle, filters=None, duplicate_handling="ignore", ignore_invalid=False, attrs = None):
+        """
+        Constructs a list with pseudopotential objects from a Pseudopotential in GAMESS format
+
+        :param fhandle: open file handle
+        :param filters: a dict with attribute filter functions
+        :param duplicate_handling: how to handle duplicates ("ignore", "error", "new" (version))
+        :param ignore_invalid: whether to ignore invalid entries silently
+        :rtype: list
+        """
+
+        def exists(pseudo):
+            try:
+                cls.get(pseudo["element"], pseudo["name"], match_aliases=False)
+            except NotExistent:
+                return False
+
+            return True
+
+        if not attrs:
+            attrs = {}
 
-    def to_cp2k(self, fhandle):
+        """
+        Parser for Gamess format
+        """
+
+        functions = []
+        ns = 0
+        for ii, line in enumerate(fhandle):
+            if len(line.strip()) == 0: continue
+            if ii == 0:
+                name, gen, core_electrons, lmax = line.split()
+                continue
+            if ns == 0:
+                ns = int(line)
+                functions.append({"prefactors" : [],
+                                  "polynoms"   : [],
+                                  "exponents"  : []})
+            else:
+                for key, value in zip(("prefactors", "polynoms", "exponents"), map(float, line.split())):
+                    functions[-1][key].append(value)
+                ns -= 1
+
+                """
+                Cast polynoms to Integers
+                """
+                functions[-1]["polynoms"] = [ int(x) for x in functions[-1]["polynoms"] ]
+
+        """
+        Change the order of functions so they match orbital momentum
+
+        In GAMESS format first block represents upper most lmax
+        and then the rest s, p, d, ...
+        """
+        functions = functions[1:] + [functions[0]]
+
+        """
+        TODO properly extract name
+        """
+
+        element = name.split("-")[0]
+        lmax = int(lmax)
+        core_electrons = int(core_electrons)
+
+
+        data = {"functions"      : functions,
+                "element"        : element,
+                "aliases"        : [name],
+                "name"           : name,
+                "core_electrons" : core_electrons,
+                "lmax"           : lmax,
+                "version"        : 1,
+                "n_el"           : None,
+                "n_el_tot"       : SYM2NUM[element] - core_electrons}
+
+        if "name" in attrs:
+            data["aliases"].append(data["name"])
+            data["name"] = attrs["name"]
+
+        if duplicate_handling == "force-ignore":  # This will be checked at the store stage
+            pass
+
+        elif duplicate_handling == "ignore":  # simply filter duplicates
+            if exists(data):
+                return []
+
+        elif duplicate_handling == "error":
+            if exists(data):
+                raise UniquenessError(
+                    f"Gaussian Pseudopotential already exists for"
+                    f" element={data['element']}, name={data['name']}: {latest.uuid}"
+                )
+
+        elif duplicate_handling == "new":
+            if exists(data):
+                latest = cls.get(data["element"], data["name"], match_aliases=False)
+                data["version"] = latest.version + 1
+
+        else:
+            raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
+
+        return [ECPPseudopotential(**data)]
+
+    @classmethod
+    def from_turborvb(cls, fhandle, filters=None, duplicate_handling="ignore", ignore_invalid=False, attrs = None, name = None):
+        """
+        Constructs a list with pseudopotential objects from a Pseudopotential in TurboRVB format
+
+        :param fhandle: open file handle
+        :param filters: a dict with attribute filter functions
+        :param duplicate_handling: how to handle duplicates ("ignore", "error", "new" (version))
+        :param ignore_invalid: whether to ignore invalid entries silently
+        :rtype: list
+        """
+
+        if hasattr(fhandle, "name"):
+            import re
+            if re.match(r"Z[0-9]{1,2}\_atomnumber[0-9]{1,2}\.[A-z]+",
+                        fhandle.name):
+                ret = re.match(r"Z[0-9]{1,2}\_atomnumber([0-9]{1,2})\.[A-z]+",
+                               fhandle.name)
+                atnum = int(ret.group(1))
+                element = list(SYM2NUM.keys())[list(SYM2NUM.values()).index(atnum)]
+                name = fhandle.name
+
+
+        def exists(pseudo):
+            try:
+                cls.get(pseudo["element"], pseudo["name"], match_aliases=False)
+            except NotExistent:
+                return False
+
+            return True
+
+        if not attrs:
+            attrs = {}
+
+        """
+        Parser for TurboRVB format
+        """
+
+        functions = []
+        ns = 0
+        for ii, line in enumerate(fhandle):
+            if ii == 0: continue
+            if ii == 1:
+                num, r0, lmax = [float(x) for x in line.split()]
+                continue
+            if ii == 2:
+                numf = [float(x) for x in line.split()]
+                for jj in range(len(numf)):
+                    functions.append({"prefactors" : [],
+                                      "polynoms"   : [],
+                                      "exponents"  : []})
+                continue
+            for jj in range(len(numf)):
+                if numf[jj] < 1: continue
+                numf[jj] -= 1
+                for key, value in zip(("prefactors", "polynoms", "exponents"), map(float, line.split())):
+                    functions[jj][key].append(value)
+
+                functions[jj]["polynoms"] = [ int(x) for x in functions[jj]["polynoms"] ]
+                break
+
+        """
+        TODO properly extract name
+        """
+
+        lmax = int(lmax)
+
+        data = {"functions"      : functions,
+                "element"        : element,
+                "aliases"        : [name],
+                "name"           : name,
+                "core_electrons" : 0,
+                "lmax"           : lmax,
+                "version"        : 1,
+                "n_el"           : None,
+                "n_el_tot"       : 0}
+
+        if "name" in attrs:
+            data["aliases"].append(data["name"])
+            data["name"] = attrs["name"]
+
+        if duplicate_handling == "force-ignore":  # This will be checked at the store stage
+            pass
+
+        elif duplicate_handling == "ignore":  # simply filter duplicates
+            if exists(data):
+                return []
+
+        elif duplicate_handling == "error":
+            if exists(data):
+                raise UniquenessError(
+                    f"Gaussian Pseudopotential already exists for"
+                    f" element={data['element']}, name={data['name']}: {latest.uuid}"
+                )
+
+        elif duplicate_handling == "new":
+            if exists(data):
+                latest = cls.get(data["element"], data["name"], match_aliases=False)
+                data["version"] = latest.version + 1
+
+        else:
+            raise ValueError(f"Specified duplicate handling strategy not recognized: '{duplicate_handling}'")
+
+        pp = ECPPseudopotential(**data)
+        pp.set_extra("r0", r0)
+        return [pp]
+
+    def to_cp2k(self, fhandle, **kwargs):
         """
         Write this Pseudopotential instance to a file in CP2K format.
 
         :param fhandle: open file handle
         """
 
-        fhandle.write(f"# from AiiDA Pseudopotential<uuid: {self.uuid}>\n")
-        for line in _dict2pseudodata(self.attributes).cp2k_format_line_iter():
-            fhandle.write(line)
+        if isinstance(self, GTHPseudopotential):
+
+            fhandle.write(f"# from AiiDA Pseudopotential<uuid: {self.uuid}>\n")
+            for line in _dict2pseudodata(self.attributes).cp2k_format_line_iter():
+                fhandle.write(line)
+                fhandle.write("\n")
+
+        else:
+            """
+            make an error
+            """
+            pass
+
+    def to_gamess(self, fhandle, **kwargs):
+        """
+        Write this Pseudopotential instance to a file in Gamess format.
+
+        :param fhandle: open file handle
+        """
+
+        if isinstance(self, ECPPseudopotential):
+            fhandle.write(f"{self.name} GEN {self.core_electrons} {self.lmax}\n")
+            functions = [self.functions[-1]] + self.functions[:-1]
+            for fun in functions:
+                fhandle.write(f"{len(fun['polynoms'])}\n")
+                for prefactor, polynom, exponent in zip(*[ fun[k] for k in ("prefactors", "polynoms", "exponents")]):
+                    fhandle.write(f"{prefactor:12.7f} {polynom:4d} {exponent:12.7f}\n")
+
+
+        else:
+            """
+            make an error
+            """
+            pass
+
+    def to_turborvb(self, fhandle, tolerance = 1.0e-5, index = 1, **kwargs):
+        """
+        Write this Pseudopotential instance to a file in TurboRVB format.
+
+        :param fhandle: open file handle
+        :param tolerance: tolerance for pseudopotential
+        """
+        def f(r, block):
+            nmax = len(block)
+            psip = np.zeros(nmax)
+            fun = 0.0
+            if r < 1.0e-9: r = 1.0e-9
+
+            for i in range(nmax):
+                psip[i] = np.exp(-block[i][2]*r*r + np.log(r)*block[i][1])
+
+            for i in range(nmax):
+                fun += psip[i] * block[i][0]
+
+            return fun/r/r
+
+        if isinstance(self, ECPPseudopotential):
+            fhandle.write(f"ECP\n")
+            r0 = 0.0
+            if "r0" in self.extras:
+                r0 = self.extras["r0"]
+            r0s = []
+            for fun in self.functions:
+                X = [ ii for ii in np.arange(0,10,0.01) ]
+                block = [ [prefactor, polynom, exponent] for prefactor, polynom, exponent in zip(*[ fun[k] for k in ("prefactors", "polynoms", "exponents")])]
+                Y = [ f(x, block) for x in X ]
+                for ii in reversed(range(len(X))):
+                    if Y[ii] > tolerance:
+                        r0s.append(X[ii])
+                        break
+            r0 = max(r0s)
+
+            fhandle.write(f"{index} {r0:4.2f} {len(self.functions)}\n")
+            fhandle.write(" ".join([ f"{len(x['polynoms'])}" for x in self.functions ]))
             fhandle.write("\n")
+            for fun in self.functions:
+                for prefactor, polynom, exponent in zip(*[ fun[k] for k in ("prefactors", "polynoms", "exponents")]):
+                    fhandle.write(f"{prefactor:12.7f} {polynom:4d} {exponent:12.7f}\n")
+
+
+        else:
+            """
+            make an error
+            """
+            pass
 
     def get_matching_basisset(self, *args, **kwargs):
         """
@@ -351,6 +722,141 @@ def get_matching_basisset(self, *args, **kwargs):
             return BasisSet.get(element=self.element, *args, **kwargs)
 
 
+class GTHPseudopotential(Pseudopotential):
+
+    __name__ = "GTHPseudopotential"
+
+    def __init__(
+        self,
+        local=None,
+        non_local=None,
+        nlcc=None,
+        **kwargs):
+        """
+        :param local: see :py:attr:`~local`
+        :param local: see :py:attr:`~non_local`
+        """
+
+        if not non_local:
+            non_local = []
+
+        if not nlcc:
+            nlcc = []
+
+        super().__init__(**kwargs)
+
+        for attr in ("local", "non_local", "nlcc"):
+            self.set_attribute(attr, locals()[attr])
+
+    @property
+    def local(self):
+        """
+        Return the local part
+
+        The format of the returned dictionary::
+
+            {
+                'r': float,
+                'coeffs': [float, float, ...],
+            }
+
+        :rtype:dict
+        """
+        return self.get_attribute("local", None)
+
+    @property
+    def non_local(self):
+        """
+        Return a list of non-local projectors (for l=0,1...).
+
+        Each list element will have the following format::
+
+            {
+                'r': float,
+                'nproj': int,
+                'coeffs': [float, float, ...],  # only the upper-triangular elements
+            }
+
+        :rtype:list
+        """
+        return self.get_attribute("non_local", [])
+
+    @property
+    def nlcc(self):
+        """
+        Return a list of the non-local core-corrections data
+
+        :rtype:list
+        """
+        return self.get_attribute("nlcc", [])
+
+    def _validate(self):
+        super()._validate()
+
+        try:
+            _dict2pseudodata(self.attributes)
+        except Exception as exc:
+            raise ValidationError("One or more invalid fields found") from exc
+
+
+class ECPPseudopotential(Pseudopotential):
+
+    __name__ = "ECPPseudopotential"
+
+    def __init__(
+        self,
+        functions=None,
+        lmax=1,
+        core_electrons=0,
+        **kwargs):
+        """
+        :param functions:
+        :param lmax: maximum angular momentum
+        """
+
+        if not functions:
+            functions = []
+
+        super().__init__(**kwargs)
+
+        for attr in ("functions", "lmax", "core_electrons"):
+            self.set_attribute(attr, locals()[attr])
+
+    @property
+    def lmax(self):
+        """
+        Return maximum angular momentum
+
+        :rtype:int
+        """
+        return self.get_attribute("lmax", [])
+
+    @property
+    def core_electrons(self):
+        """
+        Returns number of core electrons
+
+        :rtype:int
+        """
+        return self.get_attribute("core_electrons", [])
+
+    @property
+    def functions(self):
+        """
+        Returns list of basis functions
+
+        :rtype:list
+        """
+        return self.get_attribute("functions", [])
+
+
+def _dict2pseudodata(data):
+    from cp2k_input_tools.pseudopotentials import (
+        PseudopotentialData,
+        PseudopotentialDataLocal,
+        PseudopotentialDataNonLocal,
+    )
+
 def _pseudodata2dict(data: PseudopotentialData) -> Dict[str, Any]:
     """
     Convert a PseudopotentialData to a compatible dict with:
@@ -382,6 +888,6 @@ def _pseudodata2dict(data: PseudopotentialData) -> Dict[str, Any]:
 
 
 def _dict2pseudodata(data: Dict[str, Any]) -> PseudopotentialData:
-    obj = {k: v for k, v in data.items() if k not in ("name", "tags", "version")}
+    obj = {k: v for k, v in data.items() if k not in ("name", "tags", "version", "n_el_tot")}
     obj["identifiers"] = obj.pop("aliases")
     return PseudopotentialData.parse_obj(obj)
diff --git a/setup.json b/setup.json
index a7260d2..f5b1bd9 100644
--- a/setup.json
+++ b/setup.json
@@ -22,19 +22,30 @@
     "entry_points": {
         "aiida.data": [
             "gaussian.basisset = aiida_gaussian_datatypes.basisset.data:BasisSet",
-            "gaussian.pseudo = aiida_gaussian_datatypes.pseudopotential.data:Pseudopotential"
+            "gaussian.basissetfree = aiida_gaussian_datatypes.basisset.data:BasisSetFree",
+            "gaussian.pseudo = aiida_gaussian_datatypes.pseudopotential.data:Pseudopotential",
+            "gaussian.pseudo.gthpseudopotential = aiida_gaussian_datatypes.pseudopotential.data:GTHPseudopotential",
+            "gaussian.pseudo.ecppseudopotential = aiida_gaussian_datatypes.pseudopotential.data:ECPPseudopotential"
         ],
         "aiida.cmdline.data": [
             "gaussian.basisset = aiida_gaussian_datatypes.basisset.cli:cli",
-            "gaussian.pseudo = aiida_gaussian_datatypes.pseudopotential.cli:cli"
+            "gaussian.pseudo = aiida_gaussian_datatypes.pseudopotential.cli:cli",
+            "gaussian = aiida_gaussian_datatypes.fetcher.cli:cli"
         ],
         "aiida.groups": [
             "gaussian.basisset = aiida_gaussian_datatypes.groups:BasisSetGroup",
             "gaussian.pseudo = aiida_gaussian_datatypes.groups:PseudopotentialGroup"
+        ],
+        "aiida.calculations": [
+            "gaussian.uncontract = aiida_gaussian_datatypes.calc.uncontract:uncontract"
         ]
     },
     "scripts": [],
     "install_requires": [
+	"gitpython >= 3.1.24",    
+	"icecream >= 2.1.1",
+	"pydriller >= 2.0",
+        "pydantic >= 1.8.1",
         "aiida-core >= 1.6.2",
         "cp2k-input-tools >= 0.8.0"
     ],
diff --git a/tests/GAMESS_ECP.B b/tests/GAMESS_ECP.B
new file mode 100644
index 0000000..9ab1f03
--- /dev/null
+++ b/tests/GAMESS_ECP.B
@@ -0,0 +1,7 @@
+B-ccECP GEN 2 1
+3
+ 1.00000   1 30.0000
+ 100.000   3 22.0000
+-1.00000   2  5.0000
+1
+ 20.0000   2  4.0000
diff --git a/tests/test_group.py b/tests/test_group.py
index 5aceea9..afb7ab2 100644
--- a/tests/test_group.py
+++ b/tests/test_group.py
@@ -33,10 +33,11 @@ def test_pseudopotential_group_get():
     pseudogroup.add_nodes([pseudo.store() for pseudo in pseudos])
 
     retrieved_pseudos = pseudogroup.get_pseudos(elements=["Li", "H"])
+    retrieved_pseudos = {key: {x for x in val} for key, val in retrieved_pseudos.items()}
 
     assert retrieved_pseudos == {
-        "Li": [p for p in pseudos if p.element == "Li"],
-        "H": [p for p in pseudos if p.element == "H"],
+        "Li": {p for p in pseudos if p.element == "Li"},
+        "H": {p for p in pseudos if p.element == "H"},
     }
 
 
@@ -66,8 +67,9 @@ def test_pseudopotential_group_get_structure():
     structure.append_atom(position=(0.500, 0.500, 0.500), symbols="H")
 
     retrieved_pseudos = pseudogroup.get_pseudos(structure=structure)
+    retrieved_pseudos = {key: {x for x in val} for key, val in retrieved_pseudos.items()}
 
     assert retrieved_pseudos == {
-        "Li": [p for p in pseudos if p.element == "Li"],
-        "H": [p for p in pseudos if p.element == "H"],
+        "Li": {p for p in pseudos if p.element == "Li"},
+        "H": {p for p in pseudos if p.element == "H"},
     }
diff --git a/tests/test_pseudo_data.py b/tests/test_pseudo_data.py
index 03c6f9a..7539386 100644
--- a/tests/test_pseudo_data.py
+++ b/tests/test_pseudo_data.py
@@ -81,7 +81,7 @@ def test_validation_empty():
 
 
 def test_validation_invalid_local():
-    Pseudo = DataFactory("gaussian.pseudo")
+    Pseudo = DataFactory("gaussian.pseudo.gthpseudopotential")
     pseudo = Pseudo(name="test", element="H", local={"r": 1.23, "coeffs": [], "something": "else"})
 
     with pytest.raises(ValidationError):
@@ -98,6 +98,20 @@ def test_get_matching_empty():
         pseudos[0].get_matching_basisset()
 
 
+def test_import_from_gamess():
+    Pseudopotential = DataFactory("gaussian.pseudo")
+
+    with open(TEST_DIR.joinpath("GAMESS_ECP.B"), "r") as fhandle:
+        # get only the He PADE pseudo
+        pseudos = Pseudopotential.from_gamess( fhandle )
+
+    assert len(pseudos) == 1
+
+    pseudos[0].store()
+
+    # check that the name is used for the node label
+    assert pseudos[0].label == pseudos[0].name
+
 def test_to_cp2k():
     """Check whether writing a CP2K datafile works"""
     Pseudo = DataFactory("gaussian.pseudo")