diff --git a/.gitignore b/.gitignore
index 0cb92af..8ae03ca 100644
--- a/.gitignore
+++ b/.gitignore
@@ -120,4 +120,5 @@ docs/images/logo.ai
 docs/images/logo2.ai
 
 # PyCharm files
-.idea/*
\ No newline at end of file
+.idea/*
+.idea/*
diff --git a/src/matchmaps/_compute_mr_diff.py b/src/matchmaps/_compute_mr_diff.py
index 9d8f453..6d13724 100644
--- a/src/matchmaps/_compute_mr_diff.py
+++ b/src/matchmaps/_compute_mr_diff.py
@@ -3,26 +3,21 @@
 import argparse
 import os
 import sys
-import subprocess
 import time
-from functools import partial
 from pathlib import Path
 
 import gemmi
-import numpy as np
 import reciprocalspaceship as rs
 
+from matchmaps._phenix_utils import rigid_body_refinement_wrapper, phaser_wrapper, _remove_waters
 from matchmaps._utils import (
     _handle_special_positions,
     make_floatgrid_from_mtz,
-    rigid_body_refinement_wrapper,
     _realspace_align_and_subtract,
     _rbr_selection_parser,
-    _remove_waters,
     _restore_ligand_occupancy,
     _validate_environment,
     _validate_inputs,
-    phaser_wrapper,
     _clean_up_files,
     _cif_or_pdb_to_pdb,
     _cif_or_mtz_to_mtz,
@@ -51,7 +46,8 @@ def compute_mr_difference_map(
     radius : float = 5,
     alpha : float = 0,
     no_bss = False,
-):
+    phenix_version: str = None,
+    ):
     """
     Compute a real-space difference map from mtzs in different spacegroups.
 
@@ -105,7 +101,12 @@ def compute_mr_difference_map(
         If True, skip bulk solvent scaling feature of phenix.refine
     """
     
-    _validate_environment(ccp4=False)
+    auto_phenix_version = _validate_environment(ccp4=False)
+
+    if phenix_version:
+        pass
+    else:
+        phenix_version = auto_phenix_version
 
     output_dir_contents = list(output_dir.glob("*"))
     
@@ -130,15 +131,8 @@ def compute_mr_difference_map(
         f"{time.strftime('%H:%M:%S')}: Running phenix.phaser to place 'off' model into 'on' data..."
     )
 
-    phaser_nickname = phaser_wrapper(
-        mtzfile=mtzon,
-        pdb=pdboff,
-        input_dir=input_dir,
-        output_dir=output_dir,
-        off_labels=f"{Fon},{SigFon}",
-        eff=None,
-        verbose=verbose,
-    )
+    phaser_nickname = phaser_wrapper(mtzfile=mtzon, pdb=pdboff, output_dir=output_dir, off_labels=f"{Fon},{SigFon}",
+                                     phenix_style=phenix_version, eff=None, verbose=verbose)
 
     # TO-DO: fix ligand occupancies in pdb_mr_to_on
     edited_mr_pdb = _restore_ligand_occupancy(
@@ -161,6 +155,7 @@ def compute_mr_difference_map(
         off_labels=f"{Fon},{SigFon}",  # workaround for compatibility
         mr_on=True,
         no_bss=no_bss,
+        phenix_style=phenix_version,
     )
 
     print(f"{time.strftime('%H:%M:%S')}: Running phenix.refine for the 'off' data...")
@@ -175,8 +170,8 @@ def compute_mr_difference_map(
         verbose=verbose,
         rbr_selections=rbr_phenix,
         off_labels=f"{Foff},{SigFoff}",
-        mr_off=True,
         no_bss=no_bss,
+        phenix_style=phenix_version,
     )
 
     # from here down I just copied over the stuff from the normal version
@@ -443,6 +438,16 @@ def parse_arguments():
         )
     )
 
+    parser.add_argument(
+        "--phenix-version",
+        required=False,
+        help=(
+            "Specify phenix version as a string, e.g. '1.20'. "
+            "If omitted, matchmaps will attempt to automatically detect the version in use "
+            "by analyzing the output of phenix.version"
+        )
+    )
+
     return parser
 
 
@@ -485,7 +490,8 @@ def main():
         alpha=args.alpha,
         on_as_stationary=args.on_as_stationary,
         keep_temp_files=args.keep_temp_files,
-        no_bss = args.no_bss
+        no_bss = args.no_bss,
+        phenix_version = args.phenix_version,
     )
     
     if args.script:
@@ -493,6 +499,7 @@ def main():
             utility = 'matchmaps.mr', 
             arguments = sys.argv[1:],
             script_name = args.script,
+            phenix_version=args.phenix_version,
             )
 
     return
diff --git a/src/matchmaps/_compute_ncs_diff.py b/src/matchmaps/_compute_ncs_diff.py
index fcd1b5d..15c6ad1 100644
--- a/src/matchmaps/_compute_ncs_diff.py
+++ b/src/matchmaps/_compute_ncs_diff.py
@@ -3,23 +3,17 @@
 import argparse
 import os
 import sys
-import subprocess
 import time
-from functools import partial
 from pathlib import Path
 
 import gemmi
-import numpy as np
 import reciprocalspaceship as rs
 
-
+from matchmaps._phenix_utils import rigid_body_refinement_wrapper, _renumber_waters
 from matchmaps._utils import (
     _handle_special_positions,
     make_floatgrid_from_mtz,
-    rigid_body_refinement_wrapper,
-    _realspace_align_and_subtract,
     _rbr_selection_parser,
-    _renumber_waters,
     _ncs_align_and_subtract,
     _validate_environment,
     _validate_inputs,
@@ -48,6 +42,7 @@ def compute_ncs_difference_map(
     eff : str = None,
     keep_temp_files : str = None,
     no_bss = False,
+    phenix_version: str = None,
 ):
     """
     Compute an internal difference map across non-crystallographic symmetry
@@ -91,8 +86,13 @@ def compute_ncs_difference_map(
     no_bss : bool, optional
         If True, skip bulk solvent scaling feature of phenix.refine
     """
-    _validate_environment(ccp4=False)
-    
+    auto_phenix_version = _validate_environment(ccp4=False)
+
+    if phenix_version:
+        pass
+    else:
+        phenix_version = auto_phenix_version
+
     output_dir_contents = list(output_dir.glob("*"))
     
     pdb = _cif_or_pdb_to_pdb(pdb, output_dir)
@@ -121,7 +121,8 @@ def compute_ncs_difference_map(
             verbose=verbose,
             rbr_selections=rbr_phenix,
             off_labels=f"{F},{SigF}",
-            no_bss=no_bss
+            no_bss=no_bss,
+            phenix_style=phenix_version,
         )
 
         # use phenix names for columns when computing FloatGrid
@@ -327,6 +328,16 @@ def parse_arguments():
             "Note that this file is written out in the current working directory, NOT the input or output directories"
         )
     )
+
+    parser.add_argument(
+        "--phenix-version",
+        required=False,
+        help=(
+            "Specify phenix version as a string, e.g. '1.20'. "
+            "If omitted, matchmaps will attempt to automatically detect the version in use "
+            "by analyzing the output of phenix.version"
+        )
+    )
     
     return parser
 
@@ -366,6 +377,7 @@ def main():
         spacing=args.spacing,
         keep_temp_files=args.keep_temp_files,
         no_bss=args.no_bss,
+        phenix_version=args.phenix_version,
     )
     
     if args.script:
@@ -373,6 +385,7 @@ def main():
             utility = 'matchmaps.ncs', 
             arguments = sys.argv[1:],
             script_name = args.script,
+            phenix_version=args.phenix_version,
             )
 
     return
diff --git a/src/matchmaps/_compute_realspace_diff.py b/src/matchmaps/_compute_realspace_diff.py
index 9b55bc6..580bd55 100755
--- a/src/matchmaps/_compute_realspace_diff.py
+++ b/src/matchmaps/_compute_realspace_diff.py
@@ -1,27 +1,20 @@
 """Compute unbiased real space difference map."""
 
 import argparse
-import os
-import sys
-import glob
 import subprocess
+import sys
 import time
-from functools import partial
 from pathlib import Path
 
 import gemmi
-import numpy as np
 import reciprocalspaceship as rs
 
-from IPython import embed
-
+from matchmaps._phenix_utils import rigid_body_refinement_wrapper, _renumber_waters
 from matchmaps._utils import (
     _handle_special_positions,
     make_floatgrid_from_mtz,
-    rigid_body_refinement_wrapper,
     _realspace_align_and_subtract,
     _rbr_selection_parser,
-    _renumber_waters,
     _clean_up_files,
     _validate_environment,
     _validate_inputs,
@@ -51,7 +44,8 @@ def compute_realspace_difference_map(
     keep_temp_files : str = None,
     radius : float = 5,
     alpha : float = 0,
-    no_bss = False
+    no_bss = False,
+    phenix_version: str = None,
 ):
     """
     Compute a real-space difference map from mtzs.
@@ -103,9 +97,16 @@ def compute_realspace_difference_map(
         Alpha to use in error weighting of F-obs prior to Fourier Transform. Defaults to 0, e.g. no weighting.
     no_bss : bool, optional
         If True, skip bulk solvent scaling feature of phenix.refine
+    phenix_version: str, optional
+        Phenix version string to override the automatically detected version. I don't know why this would be necessary.
     """
 
-    _validate_environment(ccp4=True)
+    auto_phenix_version = _validate_environment(ccp4=True)
+
+    if phenix_version:
+        pass
+    else:
+        phenix_version = auto_phenix_version
 
     output_dir_contents = list(output_dir.glob("*"))
 
@@ -173,6 +174,7 @@ def compute_realspace_difference_map(
         verbose=verbose,
         rbr_selections=rbr_phenix,
         no_bss=no_bss,
+        phenix_style=phenix_version,
     )
 
     print(f"{time.strftime('%H:%M:%S')}: Running phenix.refine for the 'off' data...")
@@ -188,6 +190,7 @@ def compute_realspace_difference_map(
         rbr_selections=rbr_phenix,
         off_labels=f"{Foff},{SigFoff}",
         no_bss=no_bss,
+        phenix_style=phenix_version,
     )
 
     # read back in the files created by phenix
@@ -446,6 +449,16 @@ def parse_arguments():
             "Note that this file is written out in the current working directory, NOT the input or output directories"
         )
     )
+    
+    parser.add_argument(
+        "--phenix-version",
+        required=False,
+        help=(
+            "Specify phenix version as a string, e.g. '1.20'. "
+            "If omitted, matchmaps will attempt to automatically detect the version in use "
+            "by analyzing the output of phenix.version"
+        )
+    )
 
     return parser
 
@@ -484,6 +497,7 @@ def main():
         on_as_stationary=args.on_as_stationary,
         keep_temp_files=args.keep_temp_files,
         no_bss = args.no_bss,
+        phenix_version = args.phenix_version,
     )
     
     if args.script:
@@ -491,6 +505,7 @@ def main():
             utility = 'matchmaps', 
             arguments = sys.argv[1:],
             script_name = args.script,
+            phenix_version=args.phenix_version,
             )
 
     return
diff --git a/src/matchmaps/_phenix_utils.py b/src/matchmaps/_phenix_utils.py
new file mode 100644
index 0000000..968ba15
--- /dev/null
+++ b/src/matchmaps/_phenix_utils.py
@@ -0,0 +1,381 @@
+import shutil
+import subprocess
+import time
+from pathlib import Path
+
+import reciprocalspaceship as rs
+
+
+def _auto_eff_refinement_template(phenix_style: str):
+    if phenix_style == '1.20':
+        eff_contents = """
+    refinement {
+      crystal_symmetry {
+        unit_cell = cell_parameters
+        space_group = sg
+      }
+      input {
+        pdb {
+          file_name = pdb_input
+        }
+        xray_data {
+          file_name = "mtz_input"
+          labels = columns
+          r_free_flags {
+            generate=True
+          }
+          force_anomalous_flag_to_be_equal_to = False
+        }
+        monomers {
+          ligands
+        }
+      }
+      output {
+        prefix = '''nickname'''
+        serial = 1
+        serial_format = "%d"
+        job_title = '''nickname'''
+        write_def_file = False
+        write_eff_file = False
+        write_geo_file = False
+      }
+      electron_density_maps {
+        map_coefficients {
+          map_type = "2mFo-DFc"
+          mtz_label_amplitudes = "2FOFCWT"
+          mtz_label_phases = "PH2FOFCWT"
+        }
+        map_coefficients {
+          map_type = "mFo-DFc"
+          mtz_label_amplitudes = "FOFCWT"
+          mtz_label_phases = "PHFOFCWT"
+        }
+      }
+      refine {
+        strategy = *rigid_body
+        sites {
+          rigid_body_sites
+        }
+      }
+      main {
+        number_of_macro_cycles = 1
+        nproc = 8
+        bulk_solvent_and_scale=bss
+        nqh_flips=False
+      }
+    }
+        """
+    elif phenix_style == '1.21':
+        eff_contents = """
+data_manager {
+  model {
+    file = pdb_input
+  }
+  miller_array {
+    file = mtz_input
+    labels {
+      name = columns
+    }
+  }
+  fmodel {
+    xray_data {
+        r_free_flags {
+            generate = True
+        }
+    }
+  }
+}
+refinement {
+  crystal_symmetry {
+    unit_cell = cell_parameters
+    space_group = sg
+  }
+  output {
+    write_def_file = False
+    write_eff_file = False
+    write_geo_file = False
+  }
+  electron_density_maps {
+    apply_default_maps = False
+    map_coefficients {
+      map_type = "2mFo-DFc"
+      mtz_label_amplitudes = "2FOFCWT"
+      mtz_label_phases = "PH2FOFCWT"
+    }
+    map_coefficients {
+      map_type = "mFo-DFc"
+      mtz_label_amplitudes = "FOFCWT"
+      mtz_label_phases = "PHFOFCWT"
+    }
+  }
+  refine {
+    strategy = *rigid_body
+    sites {
+      rigid_body_sites
+    }
+  }
+  main {
+    number_of_macro_cycles = 1
+    nproc = 8
+    bulk_solvent_and_scale=bss
+    nqh_flips=False
+  }
+  rigid_body {
+    bulk_solvent_and_scale=bss
+    high_resolution=None
+  }
+}
+output {
+  prefix = '''nickname'''
+  serial = 1
+  serial_format = "%d"
+}
+    """
+    else:
+        raise NotImplementedError('Unsupported phenix version')
+
+    return eff_contents
+
+
+def rigid_body_refinement_wrapper(
+        mtzon,
+        pdboff,
+        input_dir,
+        output_dir,
+        phenix_style,
+        off_labels=None,
+        ligands=None,
+        eff=None,
+        verbose=False,
+        rbr_selections=None,
+        mr_on=False,
+        no_bss=False,
+):
+    if eff is None:
+        eff_contents = _auto_eff_refinement_template(phenix_style=phenix_style)
+    else:
+        with open(input_dir + eff) as file:
+            eff_contents = file.read()
+
+    if (off_labels is None) or mr_on:
+        nickname = f"{mtzon.name.removesuffix('.mtz')}_rbr_to_{pdboff.name.removesuffix('.pdb')}"
+    else:
+        nickname = f"{mtzon.name.removesuffix('.mtz')}_rbr_to_self"
+    ####
+    # update this logic in the future if matchmaps.mr changes
+    # mtz_location = input_dir if (mr_on or mr_off) else output_dir
+    ####
+    similar_files = list(output_dir.glob(f"{nickname}_[0-9]_1.*"))
+    if len(similar_files) == 0:
+        nickname += "_0"
+    else:
+        nums = []
+        for s in similar_files:
+            try:
+                nums.append(int(str(s).split("_")[-2]))
+            except ValueError:
+                pass
+        nickname += f"_{max(nums) + 1}"
+    # read in mtz to access cell parameters and spacegroup
+    cell_string, sg = _parse_mtz(mtzfile=str(mtzon))
+
+    # name for modified refinement file
+    eff = output_dir / f"params_{nickname}.eff"
+    params = {
+        "sg": sg,
+        "cell_parameters": cell_string,
+        "bss": str(not no_bss),
+        "pdb_input": str(pdboff),
+        "mtz_input": str(mtzon),
+        "nickname": str(output_dir / nickname),
+    }
+    if off_labels is None:
+        params["columns"] = "FPH1,SIGFPH1"  # names from scaleit output
+    else:
+        params["columns"] = off_labels  # user-provided column nanes
+    # if selection is not None:
+    #     params["all"] = selection  # overwrite atom selection
+    for key, value in params.items():
+        eff_contents = eff_contents.replace(key, value)
+    # either add ligands to .eff file or delete "ligands" placeholder
+    if ligands is not None:
+        ligand_string = "\n".join([f"file_name = '{l}'" for l in ligands])
+        eff_contents = eff_contents.replace("ligands", ligand_string)
+    else:
+        eff_contents = eff_contents.replace("ligands", "")
+    if rbr_selections is not None:
+        selection_string = "\n".join(
+            [f"rigid_body = '{sel}'" for sel in rbr_selections]
+        )
+        eff_contents = eff_contents.replace("rigid_body_sites", selection_string)
+    else:
+        eff_contents = eff_contents.replace("rigid_body_sites", "rigid_body = all")
+    # write out customized .eff file for use by phenix
+    with open(eff, "w") as file:
+        file.write(eff_contents)
+    # run refinement!
+    # print refinement output to terminal if user supplied the --verbose flag
+    subprocess.run(
+        f"phenix.refine {eff}",
+        shell=True,
+        capture_output=(not verbose),
+    )
+
+    return output_dir / nickname
+
+
+def _auto_eff_phaser_template(phenix_style):
+    if (phenix_style == '1.20') or (phenix_style == '1.21'):
+        eff_contents = """
+phaser {
+  mode = ANO CCA EP_AUTO *MR_AUTO MR_FRF MR_FTF MR_PAK MR_RNP NMAXYZ SCEDS
+  hklin = mtz_input
+  labin = labels
+  model = pdb_input
+  model_identity = 100
+  component_copies = 1
+  search_copies = 1
+  chain_type = *protein dna rna
+  crystal_symmetry {
+    unit_cell = cell_parameters
+    space_group = sg
+  }
+  keywords {
+    general {
+      root = '''nickname'''
+      title = '''matchmaps_MR'''
+      mute = None
+      xyzout = True
+      xyzout_ensemble = True
+      hklout = True
+      jobs = 6
+    }
+  }
+}
+        """
+    elif phenix_style == '1.21':
+        eff_contents = """"""
+
+    else:
+        raise NotImplementedError(f"Phenix version {phenix_style} not supported")
+
+    return eff_contents
+
+def phaser_wrapper(
+        mtzfile,
+        pdb,
+        output_dir,
+        off_labels,
+        phenix_style,
+        eff=None,
+        verbose=False,
+):
+    """
+    Handle simple phaser run from the command line
+
+    Args:
+        phenix_style:
+    """
+
+    # this should never be needed; environment is already checked
+    if shutil.which("phenix.phaser") is None:
+        raise OSError(
+            "Cannot find executable, phenix.phaser. Please set up your phenix environment."
+        )
+
+    if eff is None:
+        eff_contents = _auto_eff_phaser_template(phenix_style=phenix_style)
+
+    else:
+        raise NotImplementedError("Custom phaser specifications are not yet supported")
+
+    nickname = f"{mtzfile.name.removesuffix('.mtz')}_phased_with_{pdb.name.removesuffix('.pdb')}"
+
+    similar_files = list(output_dir.glob(f"{nickname}_*"))
+    if len(similar_files) == 0:
+        nickname += "_0"
+    else:
+        nums = []
+        for s in similar_files:
+            try:
+                nums.append(int(str(s).split("_")[-1].split(".")[0]))
+            except ValueError:
+                pass
+        nickname += f"_{max(nums) + 1}"
+
+    cell_string, sg = _parse_mtz(mtzfile)
+
+    eff = output_dir / f"params_{nickname}.eff"
+
+    params = {
+        "sg": sg,
+        "cell_parameters": cell_string,
+        "pdb_input": str(pdb),
+        "mtz_input": str(mtzfile),
+        "nickname": str(output_dir / nickname),
+        "labels": off_labels,  # should be prepackaged as a string
+    }
+
+    for key, value in params.items():
+        eff_contents = eff_contents.replace(key, value)
+
+    with open(eff, "w") as file:
+        file.write(eff_contents)
+
+    subprocess.run(
+        f"phenix.phaser {eff}",
+        shell=True,
+        capture_output=(not verbose),
+    )
+
+    return output_dir / nickname
+
+
+def _parse_mtz(mtzfile):
+    mtz = rs.read_mtz(str(mtzfile))
+    cell_string = f"{mtz.cell.a} {mtz.cell.b} {mtz.cell.c} {mtz.cell.alpha} {mtz.cell.beta} {mtz.cell.gamma}"
+    sg = mtz.spacegroup.short_name()
+    return cell_string, sg
+
+
+def _renumber_waters(pdb):
+    """
+    Call phenix.sort_hetatms to place waters onto the nearest protein chain.
+    This ensures that rbr selections handle waters properly
+
+    Parameters
+    ----------
+    pdb : str
+        name of pdb file
+    dir : str
+        directory in which pdb file lives
+    """
+
+    pdb_renumbered = Path(str(pdb).removesuffix(".pdb") + "_renumbered.pdb")
+
+    subprocess.run(
+        f"phenix.sort_hetatms file_name={pdb} output_file={pdb_renumbered}",
+        shell=True,
+        capture_output=True,
+    )
+
+    print(f"{time.strftime('%H:%M:%S')}: Moved waters to nearest protein chains...")
+
+    return pdb_renumbered
+
+
+def _remove_waters(
+        pdb,
+        output_dir,
+):
+    pdb_dry = pdb.name.removesuffix(".pdb") + "_dry"
+
+    subprocess.run(
+        f"phenix.pdbtools {pdb} remove='water' \
+            output.prefix='{output_dir}/' \
+            output.suffix='{pdb_dry}'",
+        shell=True,
+        capture_output=True,
+    )
+
+    return output_dir / (pdb_dry + ".pdb")
diff --git a/src/matchmaps/_utils.py b/src/matchmaps/_utils.py
index 31afaab..2bfe8f1 100644
--- a/src/matchmaps/_utils.py
+++ b/src/matchmaps/_utils.py
@@ -6,14 +6,12 @@
 """
 
 import os
-import glob
 import shutil
 import subprocess
 import time
 import re
 from functools import partial
 from pathlib import Path
-from IPython import embed
 import gemmi
 import numpy as np
 import reciprocalspaceship as rs
@@ -22,6 +20,8 @@
 def _validate_environment(ccp4):
     """
     Check if the environment contains phenix (and if necessary, ccp4) and throw a helpful error if not
+    
+    If the function runs successfully, it returns a string of form "X.XX" denoting the major and minor version of the phenix detected, e.g. "1.20" or "1.21".
     """
 
     if shutil.which("phenix.refine") is None:
@@ -31,16 +31,7 @@ def _validate_environment(ccp4):
             "For more information, see https://rs-station.github.io/matchmaps/quickstart.html#additional-dependencies"
         )
     else:
-        version_printout = subprocess.run(
-            "phenix.version | grep Version", shell=True, capture_output=True
-        )
-        
-        version_string = str(version_printout.stdout)
-        
-        if version_string.find('21') > 0:
-            raise NotImplementedError("It seems that you are using phenix 1.21, which is not yet supported by matchmaps"
-                                      "\n"
-                                      "Please use phenix 1.20 or earlier.")
+        phenix_version = _detect_phenix_version()
 
     if ccp4:
         if shutil.which("scaleit") is None:
@@ -50,6 +41,25 @@ def _validate_environment(ccp4):
                 "For more information, see https://rs-station.github.io/matchmaps/quickstart.html#additional-dependencies"
             )
 
+    print(f'Detected phenix {phenix_version} in your environment.',
+            '\n',
+            'If this is not the version you are using, please specify the version directly via the --phenix-version flag')
+    
+    return phenix_version
+
+
+def _detect_phenix_version():
+    version_printout = subprocess.run(
+        "phenix.version | grep Version", shell=True, capture_output=True
+    )
+    version_string = str(version_printout.stdout)
+    # if version_string.find('21') > 0:
+    #     raise NotImplementedError("It seems that you are using phenix 1.21, which is not yet supported by matchmaps"
+    #                               "\n"
+    #                               "Please use phenix 1.20 or earlier.")
+    phenix_version = '.'.join(version_string.split(': ')[1].split('.')[:-1])
+    return phenix_version
+
 
 def _rbr_selection_parser(rbr_selections):
     # end early and return nones if this feature isn't being used
@@ -209,162 +219,6 @@ def make_floatgrid_from_mtz(
     return float_grid
 
 
-def rigid_body_refinement_wrapper(
-    mtzon,
-    pdboff,
-    input_dir,
-    output_dir,
-    off_labels=None,
-    ligands=None,
-    eff=None,
-    verbose=False,
-    rbr_selections=None,
-    mr_on=False,
-    mr_off=False,
-    no_bss=False,
-):
-    if eff is None:
-        eff_contents = """
-refinement {
-  crystal_symmetry {
-    unit_cell = cell_parameters
-    space_group = sg
-  }
-  input {
-    pdb {
-      file_name = pdb_input
-    }
-    xray_data {
-      file_name = "mtz_input"
-      labels = columns
-      r_free_flags {
-        generate=True
-      }
-      force_anomalous_flag_to_be_equal_to = False
-    }
-    monomers {
-      ligands
-    }
-  }
-  output {
-    prefix = '''nickname'''
-    serial = 1
-    serial_format = "%d"
-    job_title = '''nickname'''
-    write_def_file = False
-    write_eff_file = False
-    write_geo_file = False
-  }
-  electron_density_maps {
-    map_coefficients {
-      map_type = "2mFo-DFc"
-      mtz_label_amplitudes = "2FOFCWT"
-      mtz_label_phases = "PH2FOFCWT"
-    }
-    map_coefficients {
-      map_type = "mFo-DFc"
-      mtz_label_amplitudes = "FOFCWT"
-      mtz_label_phases = "PHFOFCWT"
-    }
-  }
-  refine {
-    strategy = *rigid_body
-    sites {
-      rigid_body_sites
-    }
-  }
-  main {
-    number_of_macro_cycles = 1
-    nproc = 8
-    bulk_solvent_and_scale=bss
-    nqh_flips=False
-  }
-}
-    """
-    else:
-        with open(input_dir + eff) as file:
-            eff_contents = file.read()
-
-    if (off_labels is None) or (mr_on):
-        nickname = f"{mtzon.name.removesuffix('.mtz')}_rbr_to_{pdboff.name.removesuffix('.pdb')}"
-    else:
-        nickname = f"{mtzon.name.removesuffix('.mtz')}_rbr_to_self"
-
-    ####
-    # update this logic in the future if matchmaps.mr changes
-    # mtz_location = input_dir if (mr_on or mr_off) else output_dir
-    ####
-
-    similar_files = list(output_dir.glob(f"{nickname}_[0-9]_1.*"))
-    if len(similar_files) == 0:
-        nickname += "_0"
-    else:
-        nums = []
-        for s in similar_files:
-            try:
-                nums.append(int(str(s).split("_")[-2]))
-            except ValueError:
-                pass
-        nickname += f"_{max(nums)+1}"
-
-    # read in mtz to access cell parameters and spacegroup
-    mtz = rs.read_mtz(str(mtzon))
-    cell_string = f"{mtz.cell.a} {mtz.cell.b} {mtz.cell.c} {mtz.cell.alpha} {mtz.cell.beta} {mtz.cell.gamma}"
-    sg = mtz.spacegroup.short_name()
-
-    # name for modified refinement file
-    eff = output_dir / f"params_{nickname}.eff"
-
-    params = {
-        "sg": sg,
-        "cell_parameters": cell_string,
-        "bss": str(not no_bss),
-        "pdb_input": str(pdboff),
-        "mtz_input": str(mtzon),
-        "nickname": str(output_dir / nickname),
-    }
-
-    if off_labels is None:
-        params["columns"] = "FPH1,SIGFPH1"  # names from scaleit output
-    else:
-        params["columns"] = off_labels  # user-provided column nanes
-
-    # if selection is not None:
-    #     params["all"] = selection  # overwrite atom selection
-
-    for key, value in params.items():
-        eff_contents = eff_contents.replace(key, value)
-
-    # either add ligands to .eff file or delete "ligands" placeholder
-    if ligands is not None:
-        ligand_string = "\n".join([f"file_name = '{l}'" for l in ligands])
-        eff_contents = eff_contents.replace("ligands", ligand_string)
-    else:
-        eff_contents = eff_contents.replace("ligands", "")
-
-    if rbr_selections is not None:
-        selection_string = "\n".join(
-            [f"rigid_body = '{sel}'" for sel in rbr_selections]
-        )
-        eff_contents = eff_contents.replace("rigid_body_sites", selection_string)
-    else:
-        eff_contents = eff_contents.replace("rigid_body_sites", "rigid_body = all")
-
-    # write out customized .eff file for use by phenix
-    with open(eff, "w") as file:
-        file.write(eff_contents)
-
-    # run refinement!
-    # print refinement output to terminal if user supplied the --verbose flag
-    subprocess.run(
-        f"phenix.refine {eff}",
-        shell=True,
-        capture_output=(not verbose),
-    )
-
-    return output_dir / nickname
-
-
 def _handle_special_positions(pdboff, output_dir):
     """
     Check if any waters happen to sit on special positions, and if so, remove them.
@@ -423,141 +277,6 @@ def _handle_special_positions(pdboff, output_dir):
     return pdboff_nospecialpositions
 
 
-def _renumber_waters(pdb):
-    """
-    Call phenix.sort_hetatms to place waters onto the nearest protein chain. This ensures that rbr selections handle waters properly
-
-    Parameters
-    ----------
-    pdb : str
-        name of pdb file
-    dir : str
-        directory in which pdb file lives
-    """
-
-    pdb_renumbered = Path(str(pdb).removesuffix(".pdb") + "_renumbered.pdb")
-
-    subprocess.run(
-        f"phenix.sort_hetatms file_name={pdb} output_file={pdb_renumbered}",
-        shell=True,
-        capture_output=True,
-    )
-
-    print(f"{time.strftime('%H:%M:%S')}: Moved waters to nearest protein chains...")
-
-    return pdb_renumbered
-
-
-def _remove_waters(
-    pdb,
-    output_dir,
-):
-    pdb_dry = pdb.name.removesuffix(".pdb") + "_dry"
-
-    subprocess.run(
-        f"phenix.pdbtools {pdb} remove='water' \
-            output.prefix='{output_dir}/' \
-            output.suffix='{pdb_dry}'",
-        shell=True,
-        capture_output=True,
-    )
-
-    return output_dir / (pdb_dry + ".pdb")
-
-
-def phaser_wrapper(
-    mtzfile,
-    pdb,
-    input_dir,
-    output_dir,
-    off_labels,
-    eff=None,
-    verbose=False,
-):
-    """
-    Handle simple phaser run from the command line
-    """
-
-    if shutil.which("phenix.phaser") is None:
-        raise OSError(
-            "Cannot find executable, phenix.phaser. Please set up your phenix environment."
-        )
-
-    if eff is None:
-        eff_contents = """
-phaser {
-  mode = ANO CCA EP_AUTO *MR_AUTO MR_FRF MR_FTF MR_PAK MR_RNP NMAXYZ SCEDS
-  hklin = mtz_input
-  labin = labels
-  model = pdb_input
-  model_identity = 100
-  component_copies = 1
-  search_copies = 1
-  chain_type = *protein dna rna
-  crystal_symmetry {
-    unit_cell = cell_parameters
-    space_group = sg
-  }
-  keywords {
-    general {
-      root = '''nickname'''
-      title = '''matchmaps_MR'''
-      mute = None
-      xyzout = True
-      xyzout_ensemble = True
-      hklout = True
-      jobs = 6
-    }
-  }
-}
-        """
-    else:
-        raise NotImplementedError("Custom phaser specifications are not yet supported")
-
-    nickname = f"{mtzfile.name.removesuffix('.mtz')}_phased_with_{pdb.name.removesuffix('.pdb')}"
-
-    similar_files = list(output_dir.glob(f"{nickname}_*"))
-    if len(similar_files) == 0:
-        nickname += "_0"
-    else:
-        nums = []
-        for s in similar_files:
-            try:
-                nums.append(int(str(s).split("_")[-1].split(".")[0]))
-            except ValueError:
-                pass
-        nickname += f"_{max(nums)+1}"
-
-    mtz = rs.read_mtz(str(mtzfile))
-    cell_string = f"{mtz.cell.a} {mtz.cell.b} {mtz.cell.c} {mtz.cell.alpha} {mtz.cell.beta} {mtz.cell.gamma}"
-    sg = mtz.spacegroup.short_name()
-
-    eff = output_dir / f"params_{nickname}.eff"
-
-    params = {
-        "sg": sg,
-        "cell_parameters": cell_string,
-        "pdb_input": str(pdb),
-        "mtz_input": str(mtzfile),
-        "nickname": str(output_dir / nickname),
-        "labels": off_labels,  # should be prepackaged as a string
-    }
-
-    for key, value in params.items():
-        eff_contents = eff_contents.replace(key, value)
-
-    with open(eff, "w") as file:
-        file.write(eff_contents)
-
-    subprocess.run(
-        f"phenix.phaser {eff}",
-        shell=True,
-        capture_output=(not verbose),
-    )
-
-    return output_dir / nickname
-
-
 def _restore_ligand_occupancy(
     pdb_to_be_restored,
     original_pdb,
@@ -1079,12 +798,17 @@ def _clean_up_files(output_dir, old_files, keep_temp_files):
     return
 
 
-def _write_script(utility, arguments, script_name):
+def _write_script(utility, arguments, script_name, phenix_version):
     from matchmaps import __version__ as version
 
+    if phenix_version is None:
+        phenix_version = _detect_phenix_version()
+
     contents = f"""#!/bin/bash
 
-# This file was produced by matchmaps version {version} on {time.strftime('%c')}
+# This file was produced on {time.strftime('%c')}
+# Using matchmaps version {version} and phenix version {phenix_version}
+#
 # The command below was originally run in the following directory:
 # {os.getcwd()}