manage_dft_calculations.py

"""
Implementation of an automated workflow for catalyst calculations.

In particular, this script implements the following workflow:
1. For a given set of substrate and dopant species, GPAW input files and
   FireWorks workflows for a single-atom alloy (SAA) surface are generated (via
   `autocat`).
2. The SAA surface calculation is submitted via the FireWorks launch interface.
3. For successfully-completed SAA surface calculations, a single adsorbate is
   placed on the surface (at pre-specified sites), and corresponding GPAW input
   files and FireWorks workflows are generated.
4. The SAA surface + adsorbate calculation is submitted via the FireWorks launch
   interface.
5. For successfully-completed SAA surface + adsorbate calculations, the DFT
   output is parsed, binding energies calculated, and the data is converted into
   a Physical Information File (PIF) and added to an ASE DB.
"""

import os
import glob
import json
import logging
from typing import List
from typing import Dict
from typing import Tuple

import numpy as np
import matplotlib.pyplot as plt

import ase.db
from ase import io
from ase.data import reference_states
from ase.data import atomic_numbers
from ase.visualize.plot import plot_atoms

from pypif import pif
from pypif.obj import Property
from pypif.obj import FileReference

from fireworks.core.launchpad import LaunchPad
from fireworks.core.firework import Firework
from fireworks.core.fworker import FWorker
from fireworks.utilities.fw_serializers import load_object_from_file
from fireworks.queue.queue_launcher import launch_rocket_to_queue
from fireworks.user_objects.firetasks.script_task import ScriptTask

from dfttopif import directory_to_pif
from dftinputgen.gpaw import GPAWInputGenerator

from autocat.saa import generate_saa_structures
from autocat.adsorption import generate_adsorbed_structures


# Module level variables initialized on every invocation of the script
# Correspond to various FireWorks-related objects used to manage batch jobs
FW_CONFIG_DIR = os.path.expanduser("~/.fireworks")

LPAD_YAML = os.path.join(FW_CONFIG_DIR, "my_launchpad.yaml")
LPAD = LaunchPad.from_file(LPAD_YAML)

FWORKER_YAML = os.path.join(FW_CONFIG_DIR, "my_fworker.yaml")
FWORKER = FWorker.from_file(FWORKER_YAML)

QADAPTER_YAML = os.path.join(FW_CONFIG_DIR, "my_qadapter.yaml")
QADAPTER = load_object_from_file(QADAPTER_YAML)

LOG_FILE = os.path.join(FW_CONFIG_DIR, "manage.log")
logging.basicConfig(filename=LOG_FILE, level=logging.DEBUG)

# json with the reference energies (and corrections) used to calculate adsorbate
# binding energies for each catalyst system
reference_energies_json = os.path.join(
    os.path.dirname(__file__), "reference_energies.json"
)
with open(reference_energies_json, "r") as fr:
    REFERENCE_ENERGIES = json.load(fr)


def _recover_offline_runs():
    """Recovers all current "offline" fireworks in the database."""
    failed_fws = []
    recovered_fws = []
    # mongodb query: find(query, projection) [projection = what data is to be returned]
    for orun in LPAD.offline_runs.find(
        {"completed": False, "deprecated": False}, {"launch_id": 1, "fw_id": 1}
    ):
        count = LPAD.launches.count(
            {"launch_id": orun["launch_id"], "fworker.name": FWORKER.name}
        )
        if count == 0:
            continue
        fw = LPAD.recover_offline(orun["launch_id"])
        if fw:
            failed_fws.append(orun["fw_id"])
        else:
            recovered_fws.append(orun["fw_id"])
    msg = f"{len(recovered_fws)} job(s) SUCCESSFULLY recovered: {recovered_fws}"
    logging.info(msg)
    if failed_fws:
        msg = f"FAILED to recover offline fw_ids: {failed_fws}"
        logging.info(msg)


def _detect_lostruns():
    """Detects any fireworks in the database with a running status that have been lost"""
    lost_launch_ids, lost_fw_ids, inconsistent_fw_ids = LPAD.detect_lostruns(
        fizzle=True,
        expiration_secs=604800,
        # TODO: launch_query to restrict host
    )
    msg = f"{len(lost_fw_ids)} lost FWs detected and set to FIZZLED: {lost_fw_ids}"
    logging.info(msg)


def is_calc_running(calc_dir: str = None, name_check: str = None) -> bool:
    """
    Checks if the Firework in the current directory is still running on a
    FireWorker. Optionally, checks that the name of the Firework is consistent
    with the current calculation directory.
    """
    fw_json = os.path.join(calc_dir, "FW.json")
    fw_id = Firework.from_file(fw_json).fw_id
    fw = LPAD.get_fw_by_id(fw_id)
    if name_check is not None:
        assert name_check in fw.name
    if fw.state.lower() in ["reserved", "running"]:
        return True
    else:
        return False


def is_calc_success(output_traj_file: str = None, max_force: float = 0.05) -> bool:
    """Checks if the GPAW calculation has converged (force-based convergence)."""
    output_traj = io.read(output_traj_file)
    if "forces" in output_traj.calc.results:
        force_vectors = output_traj.get_forces()
        force_mags = np.linalg.norm(force_vectors, axis=1)
        if max(force_mags) < max_force:
            return True
    return False


def get_calc_status(
    output_traj_path: str = "output.traj", name_check: str = None
) -> Tuple[str, str]:
    """
    Returns the status [running/restart/complete] of a calculation corresponding
    to the GPAW output file specified. Optionally, checks for consistency in the
    names of the current calculation directory and the location of the previous
    calculation.
    """
    # Check if output files are available
    output_traj_files = glob.glob(output_traj_path)
    if not output_traj_files or os.path.getsize(output_traj_files[0]) == 0:
        # If no output files are available (or is empty), check if calculation is running
        fw_json_path = os.path.join(os.path.dirname(output_traj_path), "FW.json")
        fw_json_files = glob.glob(fw_json_path)
        if fw_json_files:
            assert len(fw_json_files) == 1
            fw_json_file = fw_json_files[0]
            calc_dir = os.path.dirname(fw_json_file)
            calc_running = is_calc_running(calc_dir=calc_dir, name_check=name_check)
            if calc_running:
                return "running", calc_dir
            else:
                # restarts if crashes on first iteration
                return "restart", calc_dir
        else:
            return "start", None
    else:
        assert len(output_traj_files) == 1
        output_traj_file = output_traj_files[0]
        calc_dir = os.path.dirname(output_traj_file)
        calc_success = is_calc_success(output_traj_file=output_traj_file)
        if calc_success:
            return "complete", calc_dir
        calc_running = is_calc_running(calc_dir=calc_dir, name_check=name_check)
        if calc_running:
            return "running", calc_dir
        else:
            return "restart", calc_dir


def restart_calc(calc_dir: str = None):
    """
    Resubmits the GPAW calculation in the specified location to the batch
    scheduler (in the Fireworks "reserved" mode).
    """
    fw_json = os.path.join(calc_dir, "FW.json")
    fw_id = Firework.from_file(fw_json).fw_id
    fw = LPAD.get_fw_by_id(fw_id)
    launch = fw.launches[-1]
    # make sure doesn't try and submit jobs from other clusters
    assert launch.fworker.name == FWORKER.name
    rerun_id = LPAD.rerun_fw(fw_id)
    logging.info(f"Rerun fw_ids: {rerun_id}")
    logging.info(f"Launch id: {launch.launch_id}")
    logging.info(f"Launch dir: {launch.launch_dir}")
    # resubmit job
    launch_rocket_to_queue(
        LPAD, FWORKER, QADAPTER, launcher_dir=launch.launch_dir, reserve=True
    )


def start_sub_calc(
    sub: str = None,
    dop: str = None,
    facet: str = None,
    basedir: str = ".",
):
    """
    Submits a surface calculation (from scratch) for the substrate/dopant SAA
    system (for the specified surface facet) in the specified location.
    """
    msg = f"Setting up {sub}/{dop}/{facet} substrate calculation..."
    logging.info(msg)

    _facet = facet.strip("fcc").strip("bcc").strip("hcp")

    saa_dict = generate_saa_structures(
        host_species=[sub],
        dopant_species=[dop],
        facets={sub: [_facet]},
        write_location=basedir,
        write_to_disk=True,
        default_lat_param_lib="beefvdw_fd",
        n_fixed_layers=2,
    )

    # Write GPAW input script
    sub_dir = os.path.dirname(saa_dict[sub][dop][facet]["traj_file_path"])
    gig = GPAWInputGenerator(
        crystal_structure=saa_dict[sub][dop][facet]["structure"],
        write_location=sub_dir,
        calculation_presets="surface_relax",
        custom_sett_dict={"h": 0.18},
    )
    gig.write_input_files()

    # submit job from scratch
    ntasks = QADAPTER.get("ntasks")
    ft = ScriptTask({"script": f"mpirun -np {ntasks} gpaw python surface_relax_in.py"})
    fw = Firework(ft, name=f"{sub_dir} Substrate Relax")

    LPAD.add_wf(fw)
    launch_rocket_to_queue(LPAD, FWORKER, QADAPTER, launcher_dir=sub_dir, reserve=True)
    msg = "Added the {sub_dir} calculation to the Fireworks launchpad"
    logging.info(msg)


def start_ads_calc(
    sub: str = None,
    ads: str = None,
    substrate_dir: str = "substrate",
):
    """
    Submits a surface calculation (from scratch) for the substrate/dopant SAA
    system (for the specified surface facet) + adsorbate placed on a suitable
    surface site, in the specified location.
    """
    msg = f"Starting the {ads} adsorbate calculation for {substrate_dir}"
    logging.info(msg)

    sub_out_traj = os.path.join(substrate_dir, "output.traj")
    sub_struct = io.read(sub_out_traj)

    # use specified sites
    bv = reference_states[atomic_numbers[sub]]["symmetry"]
    if bv == "fcc":
        s = (sub_struct[15].x, sub_struct[15].y)
    elif bv == "bcc":
        s = (sub_struct[24].x, sub_struct[24].y)
    elif bv == "hcp":
        x = sub_struct[32].x
        y = (sub_struct[32].y + sub_struct[29].y) / 2
        s = (x, y)

    sites = {"sa": [s]}

    ads_base_dir = os.path.dirname(substrate_dir)
    ads_dict = generate_adsorbed_structures(
        surface=sub_struct,
        adsorbates=[ads],
        adsorption_sites=sites,
        use_all_sites=False,
        write_location=ads_base_dir,
        write_to_disk=True,
    )

    for a in ads_dict:
        for typ in ads_dict[a]:
            for loc in ads_dict[a][typ]:
                ads_dir = os.path.dirname(ads_dict[a][typ][loc]["traj_file_path"])
                gig = GPAWInputGenerator(
                    crystal_structure=ads_dict[a][typ][loc]["structure"],
                    write_location=ads_dir,
                    calculation_presets="surface_relax",
                    custom_sett_dict={"h": 0.18},
                )
                gig.write_input_files()

                # submit job from scratch
                ntasks = QADAPTER.get("ntasks")
                ft = ScriptTask(
                    {"script": f"mpirun -np {ntasks} gpaw python surface_relax_in.py"}
                )
                fw = Firework(ft, name=f"{ads_dir} Substrate Relax")

                LPAD.add_wf(fw)
                launch_rocket_to_queue(
                    LPAD, FWORKER, QADAPTER, launcher_dir=ads_dir, reserve=True
                )
                msg = f"Added the {ads_dir} calculation to the Fireworks launchpad"
                logging.info(msg)


def _write_crystal_structure(traj_file: str = None):
    """
    Reads a crystal structure from the input ASE trajectory file, and writes
    its ASE/matplotlib-generated schematic to a PNG file.
    """
    if traj_file is None:
        return
    atoms = io.read(traj_file)

    fig = plt.figure(figsize=(5, 5))
    ax = fig.add_subplot(111)
    ax.set_axis_off()
    plot_atoms(atoms, ax, radii=0.4, rotation=("90x,45y,0z"))

    traj_basename = os.path.splitext(os.path.basename(traj_file))[0]
    image_file = os.path.join(os.path.dirname(traj_file), f"{traj_basename}.png")
    fig.savefig(image_file, bbox_inches="tight", dpi=100)


def dft_output_to_pif(
    calc_dir: str = None,
    plot_crystal_structure: bool = True,
    system_names: List[str] = None,
    system_tags: List[str] = None,
):
    """
    Parses the DFT output files and writes the parsed information as a
    serialized PIF object in "pif.json", in each of the input directories.
    If specified, writes a schematic of the crystal structure to a PNG file
    in each directory as well.
    Any input names, tags, and properties are appended to those of the
    ChemicalSystem object in the PIF.
    All raw DFT output files, parsed PIFs, crystal structure schematics are
    referenced in the PIFs w.r.t. the local (relative) directory structure.
    """
    if calc_dir is None:
        return

    logging.info(f'  Parsing calculation in "{calc_dir}"...')
    # Parse DFT output files to generate a PIF (`ChemicalSystem` object)
    chem_system = directory_to_pif(calc_dir)

    # Append any additional names
    if system_names is not None:
        logging.info(f'    Appending names "{system_names}" to the system...')
        if chem_system.names:
            chem_system.names.extend(system_names)
        else:
            chem_system.names = system_names

    # Append any additional tags specified
    if system_tags is not None:
        logging.info(f'    Appending tags "{system_tags}" to the system...')
        if chem_system.tags:
            chem_system.tags.extend(system_tags)
        else:
            chem_system.tags = system_tags

    # Attach DFT standard output files to the PIF
    stdout_files = filter(
        lambda x: os.path.splitext(x)[-1] == ".txt", os.listdir(calc_dir)
    )
    file_refs = [
        FileReference(relative_path=os.path.join(calc_dir, f)) for f in stdout_files
    ]
    chem_system.properties.append(Property(name="DFT stdout", files=file_refs))

    # Write crystal structure figures using ASE, if specified
    if plot_crystal_structure:
        traj_files = filter(
            lambda x: x in ["input.traj", "output.traj"], os.listdir(calc_dir)
        )
        for traj_file in traj_files:
            _write_crystal_structure(os.path.join(calc_dir, traj_file))

        # Attach the crystal structure image files to the PIF
        cs_image_files = filter(
            lambda x: os.path.splitext(x)[-1] == ".png", os.listdir(calc_dir)
        )
        file_refs = [
            FileReference(
                relative_path=os.path.join(calc_dir, f), mime_type="image/png"
            )
            for f in cs_image_files
        ]
        chem_system.properties.append(
            Property(name="Crystal structure schematic", files=file_refs)
        )

    # Write the PIF to a JSON file in the specified directory
    json_path = os.path.join(calc_dir, "pif.json")
    with open(json_path, "w") as fw:
        pif.dump(chem_system, fw, indent=2)
    logging.info(f"    Parsed data written to {json_path}.")


def _append_properties_to_pif(
    pif_json: str = "pif.json", properties: List[pif.Property] = None
):
    if not properties:
        return
    with open(pif_json, "r") as fr:
        system = pif.load(fr)
    for prop in properties:
        system.properties.append(prop)
    with open(pif_json, "w") as fw:
        pif.dump(system, fw, indent=2)


def _get_total_energy(system: pif.System = None) -> float:
    prop = filter(lambda x: "Total Energy" in x.name, system.properties)
    prop = list(prop)[0]
    return prop.scalars[0].value


def calculate_binding_energy(
    adsorbate: str = "N",
    adsorbate_dir: str = "adsorbate",
    substrate_dir: str = "substrate",
) -> pif.Property:
    """
    Calculates and returns binding energy of the specified adsorbate on the
    specified substrate as a pif.Property object, using:
    Binding energy = E_tot(a+s) + corrections(a+s) - E_tot(s) + corrections(a+s) - E_tot(a) + corrections(a)

    For example, for N*,
    corrections(a+s) = zpve(N*), corrections(a) = - 0.5*[zpve(N2)] + 0.5*[TdS(N2)]
    """
    ads_pif = os.path.join(adsorbate_dir, "pif.json")
    with open(ads_pif, "r") as fr:
        ads = pif.load(fr)
    sub_pif = os.path.join(substrate_dir, "pif.json")
    with open(sub_pif, "r") as fr:
        sub = pif.load(fr)
    etot_ads = _get_total_energy(ads)
    etot_sub = _get_total_energy(sub)
    etot_mol = REFERENCE_ENERGIES[adsorbate]["total_energy:mol"]
    ecor_ads = REFERENCE_ENERGIES[adsorbate].get("correction:ads", 0)
    ecor_sub = REFERENCE_ENERGIES[adsorbate].get("correction:sub", 0)
    ecor_mol = REFERENCE_ENERGIES[adsorbate].get("correction:mol", 0)
    be = etot_ads + ecor_ads - etot_sub + ecor_sub - etot_mol + ecor_mol
    return pif.Property(name="Binding energy", scalars=[pif.Scalar(value=be)])


def add_pif_to_db(
    saa_key: str = None,
    adsorbate: bool = None,
    basedir: str = ".",
    calc_dir: str = None,
    db_name: str = "dft_data.db",
):
    """
    Adds pif containing DFT data to an `ase.db`
    """
    db_path = os.path.join(basedir, db_name)
    dft_data_db = ase.db.connect(db_path, type="json")

    json_path = os.path.join(calc_dir, "pif.json")
    with open(json_path, "r") as f:
        pif_data = json.load(f)

    dft_data_db.write(atoms=None, data=pif_data, saa_key=saa_key, adsorbate=adsorbate)


def manage_calculations(
    saa_systems: List[str] = None,
    facets: Dict[str, str] = None,
    adsorbates: List[str] = None,
    basedir: str = ".",
    db_name: str = "dft_data.db",
):
    """
    Manage DFT calculations for the specified SAA catalyst systems.

    Args:

        saa_systems: Strings of the form "[substrate]-[dopant]", each
            representing an SAA catalyst.
            E.g., "Cu-Fe" where Cu is the substrate, and Fe is the dopant atom.
        facets: Dictionary with the surface facet to be considered for each SAA
            system (specifically, for the substrate element of the SAA).
            E.g., {"Cu": "111", "Ru": "0001"}
        adsorbates: Symbols of adsorbates to be calculated.
            E.g., ["N", "NH2"]
        basedir: The base directory with respect to which previous calculations
            exist and new calculations must be set up.
        db_name: Name of the ASE DB to store the obtained DFT data

    Returns:
        None. Writes the status of calculations of all the SAA systems specified
        as a JSON file to disk in the current working directory.

    """
    if saa_systems is None:
        return
    if adsorbates is None:
        adsorbates = []

    substrate_species = [saa.split("-")[0] for saa in saa_systems]

    # recover "offline" (Fireworks terminology) calculations
    logging.info("Recovering offline runs...")
    _recover_offline_runs()

    logging.info("Detecting lost runs...")
    _detect_lostruns()

    # default to only using 111 for fcc, 100 for bcc, 0001 for hcp
    facet_lookup = {"fcc": "fcc111", "bcc": "bcc110", "hcp": "hcp0001"}
    cs_library = {
        sp: reference_states[atomic_numbers[sp]].get("symmetry")
        for sp in substrate_species
    }
    _facets = {sp: facet_lookup[cs_library[sp]] for sp in substrate_species}
    if not facets:
        logging.info(f"Facets not provided. Using: {_facets}")
        facets = _facets
    else:
        facets = facets.update(**_facets)

    # check that all substrate species have a facet specified
    assert set(substrate_species) == set(facets.keys())

    # check that the specified basedir exists
    if basedir == "." or basedir == "./":
        basedir = ""
    else:
        assert os.path.isdir(basedir)

    # status of system will be one of:
    # (1) substrate-[start/restart/running/complete]
    # (2) adsorbate-[start/restart/running/complete]
    # (3) dft-complete
    status_dict = {}
    for saa_system in saa_systems:
        ssp, dsp = saa_system.split("-")
        status_dict[saa_system] = {}
        facet = facets[ssp]

        # SUBSTRATE CALCULATION CHECK
        substrate_dir = os.path.join(ssp, dsp, facet, "substrate")
        output_traj_path = os.path.join(basedir, substrate_dir, "output.traj")

        msg = f"Checking for {output_traj_path}..."
        logging.debug(msg)
        status, calc_dir = get_calc_status(
            output_traj_path=output_traj_path, name_check=substrate_dir
        )
        status_dict[saa_system]["substrate"] = status

        # appropriate next step according to status
        if status == "running":
            msg = f"{substrate_dir} calculation is still running"
            logging.info(msg)
            continue
        elif status == "restart":
            msg = f"{substrate_dir} calculation needs a restart"
            logging.info(msg)
            restart_calc(calc_dir=calc_dir)
            continue
        elif status == "start":
            msg = f"{substrate_dir} calculation needs to be set up"
            logging.info(msg)
            start_sub_calc(sub=ssp, dop=dsp, facet=facet, basedir=basedir)
            continue

        # check (redundant?) that the substrate calculation was successful
        assert status == "complete"

        # Check if system information already has already been added to the ASE DB
        pif_json = os.path.join(calc_dir, "pif.json")
        if not os.path.isfile(pif_json):
            msg = f"Parsing the {calc_dir} calculation into a PIF"
            logging.info(msg)
            system_names = [" ".join(substrate_dir.split("/"))]
            system_tags = [
                f"substrate:{ssp}",
                f"dopant:{dsp}",
                f"facet:{facet}",
                "substrate",
            ]
            dft_output_to_pif(
                calc_dir=calc_dir,
                system_names=system_names,
                system_tags=system_tags,
            )
            msg = f"Adding the {calc_dir} calculation to the db"
            logging.info(msg)
            add_pif_to_db(
                saa_key=saa_system,
                adsorbate="clean",
                basedir=basedir,
                calc_dir=calc_dir,
                db_name=db_name,
            )

        # ADSORBATE CALCULATIONS CHECK
        status_dict[saa_system]["adsorbates"] = {}
        for adsorbate in adsorbates:
            adsorbate_dir = os.path.join(ssp, dsp, facet, "adsorbates", adsorbate)
            output_traj_path = os.path.join(
                basedir, adsorbate_dir, "*", "*", "output.traj"
            )

            msg = f"Checking for {output_traj_path}..."
            logging.debug(msg)
            status, calc_dir = get_calc_status(
                output_traj_path=output_traj_path, name_check=adsorbate_dir
            )
            status_dict[saa_system]["adsorbates"][adsorbate] = status

            # appropriate next step according to status
            if status == "running":
                msg = f"{adsorbate_dir} calculation is still running"
                logging.info(msg)
                continue
            elif status == "restart":
                msg = f"{adsorbate_dir} calculation needs a restart"
                logging.info(msg)
                restart_calc(calc_dir=calc_dir)
                continue
            elif status == "start":
                msg = f"{adsorbate_dir} calculation needs to be set up"
                logging.info(msg)
                start_ads_calc(
                    sub=ssp,
                    ads=adsorbate,
                    substrate_dir=substrate_dir,
                )
                continue

            # check (redundant?) that the adsorbate calculation was successful
            assert status == "complete"

            # Check if system information has already been added to the ASE DB
            pif_json = os.path.join(calc_dir, "pif.json")
            if not os.path.isfile(pif_json):
                msg = f"Parsing the {calc_dir} calculation into a PIF"
                logging.info(msg)
                system_names = [" ".join(adsorbate_dir.split("/"))]
                site, loc = calc_dir.split("/")[-2:]
                system_tags = [
                    f"substrate:{ssp}",
                    f"dopant:{dsp}",
                    f"facet:{facet}",
                    f"adsorbate:{adsorbate}",
                    f"site:{site}",
                    f"loc:{loc}",
                ]
                dft_output_to_pif(
                    calc_dir=calc_dir,
                    system_names=system_names,
                    system_tags=system_tags,
                )
                binding_energy = calculate_binding_energy(
                    adsorbate=adsorbate,
                    adsorbate_dir=calc_dir,
                    substrate_dir=substrate_dir,
                )
                _append_properties_to_pif(
                    pif_json=pif_json, properties=[binding_energy]
                )
                msg = f"Adding the {calc_dir} calculation to the db"
                logging.info(msg)
                add_pif_to_db(
                    saa_key=saa_system,
                    adsorbate=adsorbate,
                    basedir=basedir,
                    calc_dir=calc_dir,
                    db_name=db_name,
                )

    # write the status of all systems to disk
    with open("calculation_status.json", "w") as fw:
        json.dump(status_dict, fw, indent=2)


if __name__ == "__main__":
    """
    Example usage of the calculation manager.
    Use a custom `saa_candidate_systems.json` (a list of SAA systems, see the signature
    of `manage_calculations`), uncomment the lines below, with suitable values for the
    base directory and the list of adsorbates.
    """
    # basedir = "."
    # adsorbates = ["N"]
    # saa_candidates_json = os.path.join(".", "saa_candidate_systems.json")
    # with open(saa_candidates_json, "r") as fr:
    #     saa_systems = json.load(fr)
    # manage_calculations(saa_systems=saa_systems, adsorbates=adsorbates, basedir=basedir)