diff --git a/qcengine/programs/gamess/harvester.py b/qcengine/programs/gamess/harvester.py
index 1597dbd0..aee61361 100644
--- a/qcengine/programs/gamess/harvester.py
+++ b/qcengine/programs/gamess/harvester.py
@@ -72,7 +72,12 @@ def harvest(
             return_hess = mill.align_hessian(np.array(calc_hess))
 
     else:
-        raise ValueError("""No coordinate information extracted from gamess output.""")
+        if module == "fmo":
+            return_mol = in_mol
+            return_grad = None
+            return_hess = None
+        else:
+            raise ValueError("""No coordinate information extracted from gamess output.""")
 
     return qcvars, return_hess, return_grad, return_mol, module
 
@@ -114,6 +119,7 @@ def harvest_output(outtext):
         pass_grad.append(gamessgrad)
 
     retindx = -1 if pass_coord[-1] else -2
+    if module == "fmo": retindx = -1
     return pass_qcvar[retindx], pass_coord[retindx], pass_grad[retindx], module
 
 
@@ -626,6 +632,19 @@ def harvest_outfile_pass(outtext):
                     [Decimal(mobj_dipole.group(1)), Decimal(mobj_dipole.group(2)), Decimal(mobj_dipole.group(3))]
                 )
 
+        # Process EFMO
+        mobj_efmo = re.search(
+            # fmt: off
+            r'The best FMO energy is' + r'\s+' + NUMBER,
+            # fmt: on
+            outtext,
+            re.MULTILINE | re.DOTALL,
+        )
+        if mobj_efmo:
+            qcvar["CURRENT ENERGY"] = Decimal(mobj_efmo.group(1))
+            return qcvar, qcvar_coord, qcvar_grad, "fmo"
+
+
     # Process CURRENT Energies
     if "HF TOTAL ENERGY" in qcvar:
         qcvar["CURRENT REFERENCE ENERGY"] = qcvar["HF TOTAL ENERGY"]
diff --git a/qcengine/programs/gamess/keywords.py b/qcengine/programs/gamess/keywords.py
index eb127617..3bed9d8c 100644
--- a/qcengine/programs/gamess/keywords.py
+++ b/qcengine/programs/gamess/keywords.py
@@ -14,6 +14,9 @@ def format_keyword(keyword: str, val: Any, lop_off: bool = True) -> Tuple[str, s
     elif str(val) == "False":
         text += ".false."
 
+    elif isinstance(val, list):
+        text += ",".join(map(str, val))
+
     # No Transform
     else:
         text += str(val).lower()
@@ -34,11 +37,19 @@ def format_keywords(keywords: Dict[str, Any]) -> str:
 
     grouped_lines = {}
     for group, opts in sorted(grouped_options.items()):
+        is_block = False
         line = []
         line.append(f"${group.lower()}")
         for key, val in sorted(grouped_options[group].items()):
-            line.append("=".join(format_keyword(key, val, lop_off=False)))
+            if key == "":
+                is_block = True
+                line.extend(val.split("\n"))
+            else:
+                line.append("=".join(format_keyword(key, val, lop_off=False)))
         line.append("$end\n")
-        grouped_lines[group] = textwrap.fill(" ".join(line), initial_indent=" ", subsequent_indent="  ")
+        if is_block:
+            grouped_lines[group] = "\n".join(" " + ln for ln in line)
+        else:
+            grouped_lines[group] = textwrap.fill(" ".join(line), initial_indent=" ", subsequent_indent="  ")
 
     return "\n".join(grouped_lines.values()) + "\n"
diff --git a/qcengine/programs/gamess/runner.py b/qcengine/programs/gamess/runner.py
index b456dbc4..e757f6ef 100644
--- a/qcengine/programs/gamess/runner.py
+++ b/qcengine/programs/gamess/runner.py
@@ -7,6 +7,7 @@
 
 from qcelemental.models import AtomicInput, AtomicResult, BasisSet, Provenance
 from qcelemental.util import safe_version, which
+from qcelemental import periodictable
 
 from ...exceptions import InputError, UnknownError
 from ...util import execute
@@ -102,6 +103,15 @@ def build_input(
         molcmd, moldata = input_model.molecule.to_string(dtype="gamess", units="Bohr", return_data=True)
         opts.update(moldata["keywords"])
 
+        # Split molecule $data for fmo/efmo
+        if opts.get("fmo__iefmo", None) is not None:
+            molcmd = molcmd.split("\n")
+            datacmd, fmoxyzcmd = molcmd[:3], molcmd[3:]
+            datacmd.extend(f" {el} {periodictable.to_atomic_number(el)}" for el in set(input_model.molecule.symbols))
+            datacmd.append(" $end")
+            fmoxyzcmd.insert(0, " $fmoxyz")
+            molcmd = "\n".join(fmoxyzcmd) + "\n".join(datacmd) + "\n"
+
         # Handle calc type and quantum chemical method
         opts.update(muster_modelchem(input_model.model.method, input_model.driver.derivative_int()))