diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py
index fca61fd..b6f55a3 100644
--- a/src/qibochem/driver/molecule.py
+++ b/src/qibochem/driver/molecule.py
@@ -42,7 +42,7 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file
         else:
             # Check if xyz_file exists, then fill in the Molecule attributes
             assert Path(f"{xyz_file}").exists(), f"{xyz_file} not found!"
-            self._process_xyz_file(xyz_file)
+            self._process_xyz_file(xyz_file, charge, multiplicity)
         if basis is None:
             # Default bais is STO-3G
             self.basis = "sto-3g"
@@ -80,7 +80,7 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file
         self.n_active_e = None
         self.n_active_orbs = None  # Number of spin-orbitals in the active space
 
-    def _process_xyz_file(self, xyz_file):
+    def _process_xyz_file(self, xyz_file, charge, multiplicity):
         """
         Reads a .xyz file to obtain and set the molecular coordinates (in OpenFermion format),
             charge, and multiplicity
@@ -91,7 +91,15 @@ def _process_xyz_file(self, xyz_file):
         with open(xyz_file, encoding="utf-8") as file_handler:
             # First two lines: # atoms and comment line (charge, multiplicity)
             _n_atoms = int(file_handler.readline())  # Not needed/used
-            _charge, _multiplicity = (int(_num) for _num in file_handler.readline().split())
+
+            # Try to read charge and multiplicity from comment line
+            split_line = [int(_num) for _num in file_handler.readline().split()]
+            if len(split_line) == 2:
+                # Format of comment line matches (charge, multiplicity):
+                _charge, _multiplicity = split_line
+            else:
+                # Otherwise, use the default (from __init__) values of 0 and 1
+                _charge, _multiplicity = charge, multiplicity
 
             # Start reading xyz coordinates from the 3rd line onwards
             _geometry = []