diff --git a/conformers.py b/conformers.py
index 464bd30..c493372 100644
--- a/conformers.py
+++ b/conformers.py
@@ -133,7 +133,7 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # If provided crest path different to that stored, use it and save it
         if Path(avo_input["crest_bin"]) != easyxtb.CREST_BIN:
@@ -157,8 +157,6 @@
         # Run calculation; returns set of conformers as well as Calculation object
         conformers, calc = easyxtb.calc.conformers(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=solvation,
             method=2,
             ewin=ewin_kcal,
diff --git a/deprotonate.py b/deprotonate.py
index e1f0023..a0e8bcb 100644
--- a/deprotonate.py
+++ b/deprotonate.py
@@ -39,13 +39,11 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns set of tautomers as well as Calculation object
         tautomers, calc = easyxtb.calc.deprotonate(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=2,
             return_calc=True,
diff --git a/easyxtb/src/easyxtb/calc.py b/easyxtb/src/easyxtb/calc.py
index 0968623..0455ce8 100644
--- a/easyxtb/src/easyxtb/calc.py
+++ b/easyxtb/src/easyxtb/calc.py
@@ -172,6 +172,10 @@ def run(self):
                     *self.runtype_args,
                     "--",
                 ]
+            if "chrg" not in self.options and self.input_geometry.charge != 0:
+                command.extend(["--chrg", str(self.input_geometry.charge)])
+            if "uhf" not in self.options and self.input_geometry.multiplicity != 1:
+                command.extend(["--uhf", str(self.input_geometry.multiplicity - 1)])
             for flag, value in self.options.items():
                 # Add appropriate number of minuses to flags
                 if len(flag) == 1:
@@ -216,7 +220,11 @@ def process_xtb(self):
         # If there's an output geometry, read it
         if self.output_file.with_name("xtbopt.xyz").exists():
             logger.debug(f"Output geometry found at {self.output_file.with_name('xtbopt.xyz')}")
-            self.output_geometry = Geometry.from_file(self.output_file.with_name("xtbopt.xyz"))
+            self.output_geometry = Geometry.from_file(
+                self.output_file.with_name("xtbopt.xyz"),
+                charge=self.input_geometry.charge,
+                multiplicity=self.input_geometry.multiplicity,
+            )
             logger.debug("Read output geometry")
         else:
             # Assume geom was the same at end of calc as at beginning
@@ -248,7 +256,11 @@ def process_crest(self):
                 # Conformer search
                 logger.debug("A conformer search was requested, so checking for files")
                 # Get energy and geom of lowest conformer
-                best = Geometry.from_file(self.output_file.with_name("crest_best.xyz"))
+                best = Geometry.from_file(
+                    self.output_file.with_name("crest_best.xyz"),
+                    charge=self.input_geometry.charge,
+                    multiplicity=self.input_geometry.multiplicity,
+                )
                 logger.debug(f"Geometry of lowest energy conformer read from {self.output_file.with_name('crest_best.xyz')}")
                 self.output_geometry = best
                 self.energy = float(best._comment)
@@ -257,6 +269,8 @@ def process_crest(self):
                 conformer_geoms = Geometry.from_file(
                     self.output_file.with_name("crest_conformers.xyz"),
                     multi=True,
+                    charge=self.input_geometry.charge,
+                    multiplicity=self.input_geometry.multiplicity,
                 )
                 logger.debug(f"Geometries of conformers read from {self.output_file.with_name('crest_conformers.xyz')}")
                 self.conformers = [
@@ -273,6 +287,8 @@ def process_crest(self):
                 tautomer_geoms = Geometry.from_file(
                     self.output_file.with_name("protonated.xyz"),
                     multi=True,
+                    charge=self.input_geometry.charge + 1,
+                    multiplicity=self.input_geometry.multiplicity,
                 )
                 logger.debug(f"Geometries of tautomers read from {self.output_file.with_name('protonated.xyz')}")
                 self.tautomers = [
@@ -294,6 +310,8 @@ def process_crest(self):
                 tautomer_geoms = Geometry.from_file(
                     self.output_file.with_name("deprotonated.xyz"),
                     multi=True,
+                    charge=self.input_geometry.charge - 1,
+                    multiplicity=self.input_geometry.multiplicity,
                 )
                 logger.debug(f"Geometries of tautomers read from {self.output_file.with_name('deprotonated.xyz')}")
                 self.tautomers = [
@@ -314,8 +332,6 @@ def process_crest(self):
 
 def energy(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     n_proc: int | None = None,
@@ -326,8 +342,6 @@ def energy(
     calc = Calculation(
         input_geometry=input_geometry,
         options={
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "gfn": method,
             "alpb": solvation,
             "p": n_proc if n_proc else config["n_proc"],
@@ -342,8 +356,6 @@ def energy(
 
 def optimize(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     level: str = "normal",
@@ -358,8 +370,6 @@ def optimize(
         runtype="opt",
         runtype_args=[level],
         options={
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "gfn": method,
             "alpb": solvation,
             "p": n_proc if n_proc else config["n_proc"],
@@ -377,8 +387,6 @@ def optimize(
 
 def frequencies(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     n_proc: int | None = None,
@@ -390,8 +398,6 @@ def frequencies(
         input_geometry=input_geometry,
         runtype="hess",
         options={
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "gfn": method,
             "alpb": solvation,
             "p": n_proc if n_proc else config["n_proc"],
@@ -406,8 +412,6 @@ def frequencies(
 
 def opt_freq(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     level: str = "normal",
@@ -426,8 +430,6 @@ def opt_freq(
         runtype="ohess",
         runtype_args=[level],
         options={
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "gfn": method,
             "alpb": solvation,
             "p": n_proc if n_proc else config["n_proc"],
@@ -441,12 +443,10 @@ def opt_freq(
     distorted_geom_file = calc.output_file.with_name("xtbhess.xyz")
     while distorted_geom_file.exists():
         calc = Calculation(
-            input_geometry=Geometry.from_file(distorted_geom_file),
+            input_geometry=Geometry.from_file(distorted_geom_file, charge=input_geometry.charge, multiplicity=input_geometry.multiplicity),
             runtype="ohess",
             runtype_args=[level],
             options={
-                "chrg": charge,
-                "uhf": multiplicity - 1,
                 "gfn": method,
                 "alpb": solvation,
                 "p": n_proc if n_proc else config["n_proc"],
@@ -462,8 +462,6 @@ def opt_freq(
 
 def orbitals(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     n_proc: int | None = None,
@@ -479,8 +477,6 @@ def orbitals(
         input_geometry=input_geometry,
         options={
             "molden": True,
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "gfn": method,
             "alpb": solvation,
             "p": n_proc if n_proc else config["n_proc"],
@@ -495,8 +491,6 @@ def orbitals(
 
 def conformers(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     ewin: int | float = 6,
@@ -519,8 +513,6 @@ def conformers(
         options={
             "xnam": XTB_BIN,
             method_flag: True,
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "alpb": solvation,
             "ewin": ewin,
             "T": n_proc if n_proc else config["n_proc"],
@@ -537,8 +529,6 @@ def conformers(
 
 def protonate(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     n_proc: int | None = None,
@@ -547,8 +537,6 @@ def protonate(
     """Screen possible protonation sites and return set of tautomer Geometries and energies.
 
     The returned tautomers are ordered from lowest to highest energy.
-    
-    Note that `charge` should be the charge on the molecule before protonation.
     """
     method_flag = f"gfn{method}"
     calc = Calculation(
@@ -558,8 +546,6 @@ def protonate(
         options={
             "xnam": XTB_BIN,
             method_flag: True,
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "alpb": solvation,
             "T": n_proc if n_proc else config["n_proc"],
         },
@@ -573,8 +559,6 @@ def protonate(
 
 def deprotonate(
     input_geometry: Geometry,
-    charge: int = 0,
-    multiplicity: int = 1,
     solvation: str | None = None,
     method: int = 2,
     n_proc: int | None = None,
@@ -583,8 +567,6 @@ def deprotonate(
     """Screen possible deprotonation sites and return set of tautomer Geometries and energies.
 
     The returned tautomers are ordered from lowest to highest energy.
-    
-    Note that `charge` should be the charge on the molecule before deprotonation.
     """
     method_flag = f"gfn{method}"
     calc = Calculation(
@@ -594,8 +576,6 @@ def deprotonate(
         options={
             "xnam": XTB_BIN,
             method_flag: True,
-            "chrg": charge,
-            "uhf": multiplicity - 1,
             "alpb": solvation,
             "T": n_proc if n_proc else config["n_proc"],
         },
diff --git a/easyxtb/src/easyxtb/geometry.py b/easyxtb/src/easyxtb/geometry.py
index 120d6cd..9ff3000 100644
--- a/easyxtb/src/easyxtb/geometry.py
+++ b/easyxtb/src/easyxtb/geometry.py
@@ -26,9 +26,12 @@ class Geometry:
     def __init__(
         self,
         atoms: list[Atom],
+        charge: int = 0,
+        multiplicity: int = 1,
         _comment: str | None = None,
     ):
-        """A set of atoms within a 3D space for use in calculations.
+        """A set of atoms within a 3D space for use in calculations with an associated
+        charge and multiplicity.
 
         Provides class methods for creation from, and instance methods for writing to,
         XYZ and CJSON formats.
@@ -36,6 +39,11 @@ def __init__(
         The coordinates should be provided as a list of `Atom` objects.
         """
         self.atoms = atoms
+        self.charge = charge
+        if multiplicity >= 1:
+            self.multiplicity = multiplicity
+        else:
+            raise ValueError("Multiplicity must be at least 1 (a singlet).")
         self._comment = _comment
     
     def __iter__(self):
@@ -78,6 +86,10 @@ def to_cjson(self) -> dict:
                     "number": all_element_numbers,
                 },
             },
+            "properties": {
+                "totalCharge": self.charge,
+                "totalSpinMultiplicity": self.multiplicity,
+            }
         }
         return cjson
 
@@ -117,7 +129,7 @@ def to_file(self, dest: os.PathLike, format: str = None) -> os.PathLike:
             self.write_cjson(dest)
 
     @classmethod
-    def from_xyz(cls, xyz_lines: list[str]):
+    def from_xyz(cls, xyz_lines: list[str], charge: int = 0, multiplicity: int = 1):
         """Create a `Geometry` object from an XYZ in the form of a list of lines."""
         logger.debug("Instantiating a geometry from an xyz")
         atoms = []
@@ -130,10 +142,10 @@ def from_xyz(cls, xyz_lines: list[str]):
                 atoms.append(
                     Atom(atom_parts[0], *[float(n) for n in atom_parts[1:4]])
                 )
-        return Geometry(atoms, _comment=xyz_lines[1])
+        return Geometry(atoms, charge, multiplicity, _comment=xyz_lines[1])
     
     @classmethod
-    def from_multi_xyz(cls, xyz_lines: list[str]):
+    def from_multi_xyz(cls, xyz_lines: list[str], charge: int = 0, multiplicity: int = 1):
         """Create a set of `Geometry` objects from an XYZ (in the form of a list of
         lines) that contains multiple different structures.
         
@@ -149,16 +161,20 @@ def from_multi_xyz(cls, xyz_lines: list[str]):
         for i, l in enumerate(xyz_lines):
             current_xyz.append(l)
             if i == len(xyz_lines) - 1 or xyz_lines[i+1] == atom_count_line:
-                geometries.append(cls.from_xyz(current_xyz))
+                geometries.append(cls.from_xyz(current_xyz, charge, multiplicity))
                 current_xyz = []
         return geometries
     
     @classmethod
-    def from_cjson(cls, cjson_dict: dict):
-        """Create a `Geometry` object from an CJSON in the form of a Python dict."""
+    def from_cjson(cls, cjson_dict: dict, charge: int = None, multiplicity: int = None):
+        """Create a `Geometry` object from an CJSON in the form of a Python dict.
+        
+        If the CJSON does not specify the overall charge and multiplicity, a neutral
+        singlet is assumed, regardless of the chemical feasibility of that.
+        """
         logger.debug("Instantiating a geometry from a cjson")
         atoms = []
-        for i, atomic_no in enumerate(cjson_dict["atoms"]["elements"]):
+        for i, atomic_no in enumerate(cjson_dict["atoms"]["elements"]["number"]):
             atoms.append(
                 Atom(
                     get_element_symbol(atomic_no),
@@ -167,16 +183,28 @@ def from_cjson(cls, cjson_dict: dict):
                     cjson_dict["atoms"]["coords"]["3d"][3*i + 2],
                 )
             )
-        return Geometry(atoms)
+        charge = charge if charge else cjson_dict.get("properties", {}).get("totalCharge", 0)
+        multiplicity = multiplicity if multiplicity else cjson_dict.get("properties", {}).get("totalSpinMultiplicity", 1)
+        return Geometry(atoms, charge, multiplicity)
 
     @classmethod
-    def from_file(cls, file: os.PathLike, format: str = None, multi: bool = False):
+    def from_file(
+        cls,
+        file: os.PathLike,
+        format: str = None,
+        multi: bool = False,
+        charge: int = None,
+        multiplicity: int = None,
+    ):
         """Create a `Geometry` object from an XYZ or CJSON file.
         
         The format can be specified by passing either ".xyz" or ".cjson" as the `format`
         argument, or it can be left to automatically be detected based on the filename
         ending.
 
+        If the file is a CJSON, charge and multiplicity will be read from the file if
+        present. Passing them as arguments will override this, but is not required.
+
         The method attempts to automatically detect an XYZ file containing multiple
         structures, parse as appropriate, and return a list of `Geometry` objects.
         If you know that a file contains multiple structures, or you wish to make sure
@@ -192,15 +220,17 @@ def from_file(cls, file: os.PathLike, format: str = None, multi: bool = False):
         if format == ".xyz":
             with open(filepath, encoding="utf-8") as f:
                 xyz_lines = f.read().split("\n")
+            charge = charge if charge else 0
+            multiplicity = multiplicity if multiplicity else 1
             # Detect multiple structures in single xyz file
             # Make the reasonable assumption that all structures have the same number of
             # atoms and that therefore the first line of the file repeats itself
             atom_count_line = xyz_lines[0]
             n_structures = xyz_lines.count(atom_count_line)
             if n_structures > 1 or multi:
-                return cls.from_multi_xyz(xyz_lines)
+                return cls.from_multi_xyz(xyz_lines, charge=charge, multiplicity=multiplicity)
             else:
-                return cls.from_xyz(xyz_lines)
+                return cls.from_xyz(xyz_lines, charge=charge, multiplicity=multiplicity)
 
         if format == ".cjson":
             with open(filepath, encoding="utf-8") as f:
diff --git a/energy.py b/energy.py
index a9f8af4..b5583cf 100644
--- a/energy.py
+++ b/energy.py
@@ -38,14 +38,12 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns energy as float in hartree
         logger.debug("avo_xtb is requesting a single point energy calculation")
         energy_hartree = easyxtb.calc.energy(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=easyxtb.config["method"],
         )
diff --git a/freq.py b/freq.py
index a622918..eaf771b 100644
--- a/freq.py
+++ b/freq.py
@@ -38,14 +38,12 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns set of frequency data
         logger.debug("avo_xtb is requesting a frequency calculation")
         freqs = easyxtb.calc.frequencies(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=easyxtb.config["method"],
         )
diff --git a/ohess.py b/ohess.py
index 8efdf7b..33d5dec 100644
--- a/ohess.py
+++ b/ohess.py
@@ -39,14 +39,12 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns path to Gaussian file containing frequencies
         logger.debug("avo_xtb is requesting an opt+freq calculation")
         opt_geom, freqs, calc = easyxtb.calc.opt_freq(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=easyxtb.config["method"],
             level=easyxtb.config["opt_lvl"],
diff --git a/opt.py b/opt.py
index fc0a766..f0ae138 100644
--- a/opt.py
+++ b/opt.py
@@ -57,14 +57,12 @@ def cleanup_after_opt(cjson: dict) -> dict:
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns optimized geometry as well as Calculation object
         logger.debug("avo_xtb is requesting a geometry optimization")
         opt_geom, calc = easyxtb.calc.optimize(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=easyxtb.config["method"],
             level=easyxtb.config["opt_lvl"],
diff --git a/orbitals.py b/orbitals.py
index 9e948c3..fd40cb0 100644
--- a/orbitals.py
+++ b/orbitals.py
@@ -38,14 +38,12 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns Molden output file as string
         logger.debug("avo_xtb is requesting a molecular orbitals calculation")
         molden_string = easyxtb.calc.orbitals(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=easyxtb.config["method"],
         )
diff --git a/protonate.py b/protonate.py
index 80f8346..9d2f7b5 100644
--- a/protonate.py
+++ b/protonate.py
@@ -52,13 +52,11 @@ def cleanup_after_taut(cjson: dict) -> dict:
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         # Run calculation; returns set of tautomers as well as Calculation object
         tautomers, calc = easyxtb.calc.protonate(
             geom,
-            charge=avo_input["charge"],
-            multiplicity=avo_input["spin"],
             solvation=easyxtb.config["solvent"],
             method=2,
             return_calc=True,
diff --git a/run.py b/run.py
index 37b51a6..ddc39b2 100644
--- a/run.py
+++ b/run.py
@@ -94,7 +94,7 @@
         # Read input from Avogadro
         avo_input = json.loads(sys.stdin.read())
         # Extract the coords
-        geom = easyxtb.Geometry.from_xyz(avo_input["xyz"].split("\n"))
+        geom = easyxtb.Geometry.from_cjson(avo_input["cjson"])
 
         command = avo_input["command"].split()
         # Remove xtb and <input_geometry> placeholders