diff --git a/tangelo/toolboxes/molecular_computation/frozen_orbitals.py b/tangelo/toolboxes/molecular_computation/frozen_orbitals.py
index dc81038e9..6bf5cecb9 100644
--- a/tangelo/toolboxes/molecular_computation/frozen_orbitals.py
+++ b/tangelo/toolboxes/molecular_computation/frozen_orbitals.py
@@ -42,7 +42,7 @@ def get_frozen_core(molecule):
 
     # Counting how many of each element is in the molecule.
     elements = {i: molecule.elements.count(i) for i in molecule.elements}
-    frozen_core = sum([v * core_orbitals[k] for k, v in elements.items()])
+    frozen_core = sum([v * core_orbitals.get(k, 0) for k, v in elements.items()])
 
     return frozen_core
 
diff --git a/tangelo/toolboxes/molecular_computation/molecule.py b/tangelo/toolboxes/molecular_computation/molecule.py
index 1265f7329..7306e1de1 100644
--- a/tangelo/toolboxes/molecular_computation/molecule.py
+++ b/tangelo/toolboxes/molecular_computation/molecule.py
@@ -94,7 +94,7 @@ class Molecule:
     n_min_orbitals: int = field(init=False)
 
     def __post_init__(self):
-        mol = self.to_pyscf(basis="sto-3g", symmetry=False)
+        mol = self.to_pyscf()
         self.n_atoms = mol.natm
         self.n_electrons = mol.nelectron
         self.n_min_orbitals = mol.nao_nr()
@@ -107,11 +107,13 @@ def elements(self):
     def coords(self):
         return np.array([self.xyz[i][1] for i in range(self.n_atoms)])
 
-    def to_pyscf(self, basis="sto-3g", symmetry=False):
+    def to_pyscf(self, basis="CRENBL", symmetry=False, ecp=None):
         """Method to return a pyscf.gto.Mole object.
 
         Args:
             basis (string): Basis set.
+            symmetry (bool): Flag to turn symmetry on
+            ecp (dict): Dictionary with ecp definition for each atom e.g. {"Cu": "crenbl"}
 
         Returns:
             pyscf.gto.Mole: PySCF compatible object.
@@ -122,6 +124,7 @@ def to_pyscf(self, basis="sto-3g", symmetry=False):
         mol.charge = self.q
         mol.spin = self.spin
         mol.symmetry = symmetry
+        mol.ecp = ecp if ecp else dict()
         mol.build()
 
         self.xyz = list()
@@ -137,7 +140,7 @@ def to_file(self, filename, format=None):
             filename (str): The name of the file to output the geometry.
             format (str): The output type of "raw", "xyz", or "zmat". If None, will be inferred by the filename
         """
-        mol = self.to_pyscf(basis="sto-3g")
+        mol = self.to_pyscf()
         mol.tofile(filename, format)
 
     def to_openfermion(self, basis="sto-3g"):
@@ -161,6 +164,8 @@ class SecondQuantizedMolecule(Molecule):
 
     Attributes:
         basis (string): Basis set.
+        ecp (dict): The effective core potential (ecp) for any atoms in the molecule.
+            e.g. {"C": "crenbl"} use CRENBL ecp for Carbon atoms.
         symmetry (bool or str): Whether to use symmetry in RHF or ROHF calculation.
             Can also specify point group using pyscf allowed string.
             e.g. "Dooh", "D2h", "C2v", ...
@@ -195,6 +200,7 @@ class SecondQuantizedMolecule(Molecule):
         actives_mos (list): Active MOs indexes.
     """
     basis: str = "sto-3g"
+    ecp: dict = field(default_factory=dict)
     symmetry: bool = False
     frozen_orbitals: list or int = field(default="frozen_core", repr=False)
 
@@ -273,7 +279,7 @@ def _compute_mean_field(self):
         (mf_energy, mo_energies, mo_occ, n_mos and n_sos).
         """
 
-        molecule = self.to_pyscf(self.basis, self.symmetry)
+        molecule = self.to_pyscf(self.basis, self.symmetry, self.ecp)
 
         self.mean_field = scf.RHF(molecule)
         self.mean_field.verbose = 0
@@ -332,7 +338,7 @@ def _convert_frozen_orbitals(self, frozen_orbitals):
         """
 
         if frozen_orbitals == "frozen_core":
-            frozen_orbitals = get_frozen_core(self.to_pyscf(self.basis))
+            frozen_orbitals = get_frozen_core(self.to_pyscf(self.basis)) if not self.ecp else 0
         elif frozen_orbitals is None:
             frozen_orbitals = 0
 
@@ -456,7 +462,7 @@ def get_integrals(self, mo_coeff=None, consider_frozen=True):
         """
 
         # Pyscf molecule to get integrals.
-        pyscf_mol = self.to_pyscf(self.basis, self.symmetry)
+        pyscf_mol = self.to_pyscf(self.basis, self.symmetry, self.ecp)
         if mo_coeff is None:
             mo_coeff = self.mean_field.mo_coeff
 
diff --git a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py
index 347b55368..cf0a8af12 100644
--- a/tangelo/toolboxes/molecular_computation/tests/test_molecule.py
+++ b/tangelo/toolboxes/molecular_computation/tests/test_molecule.py
@@ -184,6 +184,20 @@ def test_symmetry_label(self):
         assert(mo_symm_labels == molecule.mo_symm_labels)
         assert(mo_symm_ids == molecule.mo_symm_ids)
 
+    def test_ecp(self):
+        """Verify that the number of electrons is reduced when ecp is called."""
+
+        molecule = SecondQuantizedMolecule(xyz="Yb", q=0, spin=0, basis="crenbl", ecp="crenbl")
+        # "Yb" has 70 electrons but the ecp reduces this to 16
+        assert(molecule.n_active_electrons == 16)
+        assert(molecule.n_active_mos == 96)
+
+        molecule = SecondQuantizedMolecule(xyz="Cu", q=0, spin=1, basis="cc-pvdz", ecp="crenbl",
+                                           frozen_orbitals=list(range(8)))
+        # "Cu" has 29 electrons but the ecp reduces this to 19. The active electrons are 19 - 8 * 2 = 3
+        assert(molecule.n_active_electrons == 3)
+        assert(molecule.n_active_mos == 35)
+
 
 if __name__ == "__main__":
     unittest.main()