def _load(self, **kwargs):
        """ Initialize OrcaHess Hessian object from .hess file

        Searches indicated file for data blocks within the .hess file.  The
        geometry, Hessian block, frequencies, and normal modes must be present
        and will be retrieved; other blocks will be imported if present, or
        ignored if absent (stored as |None| in the resulting object).
        If malformed/inaccurate data is found in any
        block that is present, some flavor of :exc:`~opan.error.HessError`
        will be raised.

        Parameters
        ----------
        path
            |str| --
            `kwargs` parameter specifying complete path to the
            .hess file to be read.

        Raises
        ------
        ~opan.error.HessError
            (various typecodes) If indicated Hessian file
            is malformed in some fashion.

        ~excoptions.KeyError
            If invalid atomic symbol appears in .hess file.

        ~exceptions.IOError
            If the indicated file does not exist or cannot be read.


        """

        # Local method(s)
        def parse_multiblock(hesstext, p_block, p_sec, p_line, num_ats,
                                                            blockname, tc):
            """ Helper function for importing blocks with multiple sections.

            Parsing of data spanning multiple sections of columns is somewhat
            involved. This function encapsulates the process for cleaner
            code.  The structure depends critically on several formatting
            features of |orca| .hess files.

            Note the search groups that must be present in the `p_block` and
            `p_line` Regex patterns.

            Parameters
            ----------
            hesstest
                |str| --
                Complete text of the .hess file

            p_block
                |re.compile| pattern --
                Retrieves the **entirety** of the relevant block
                Required groups:
                    dim     : overall dimension of the data block
                    block   : contents of the block, including column headers

            p_sec
                |re.compile| pattern --
                Retrieves each section of data columns

            p_line
                |re.compile| pattern --
                Retrieves individual lines of a section
                Required groups:
                    row     : Row index into the final data block
                    e#      : Column index into the data section (# in range(6))

            num_ats
                |int| --
                Number of atoms in the geometry

            blockname
                |str| --
                Brief text description of the block being imported, if needed
                for error reporting purposes

            tc
                :class:`~opan.error.HessError` typecode --
                Type of error to be thrown, if required

            Returns
            -------
            workmtx
                3N x 3N |npfloat_| --
                Assembled array for storage

            """

            # Pull the block
            m_block = p_block.search(hesstext)

            # Confirm the anticipated matrix size matches 3*num_ats
            if not int(m_block.group("dim")) == 3*num_ats:
                raise HessError(tc, blockname + " dimension " +
                        "specification mismatched with geometry",
                        "HESS File: " + hess_path)
            ## end if

            # Initialize the working matrix
            workmtx = np.zeros((3*num_ats, 3*num_ats), dtype=np.float_)

            # Initialize the column offset for populating the matrix
            col_offset = 0

            # Initialize the counter for the number of row sections imported
            rows_counter = 0

            # Loop through the subsections of the matrix
            for m_sec in p_sec.finditer(m_block.group("block")):
                # Loop through each entry
                for m_line in p_line.finditer(m_sec.group(0)):
                    # Store the row; expect base zero so no adjustment needed.
                    rowval = scast(m_line.group("row"), np.int_)

                    # Loop to fill the row
                    for i in range(6):
                        # Store the element
                        val = scast(m_line.group('e{0}'.format(i)), np.float_)

                        # Only store to matrix if a value actually retrieved.
                        #  This protects against the final three-column section
                        #  in blocks for systems with an odd number of atoms.
                        if not np.isnan(val):
                            workmtx[rowval, col_offset + i] = val
                        ## end if
                    ## next i

                    # Increment the row-read counter
                    rows_counter += 1

                ## next m_hess_line

                # Last thing is to increment the offset by six. Don't have to
                #  worry about offsetting for a section only three columns wide
                #  because that SHOULD only ever occur at the last section
                col_offset += 6

            ## next m_hess_sec

            # Check to ensure that the column offset is high enough to have
            #  fully populated the matrix based on the number of atoms. This is
            #  a check against malformation of the HESS file, resulting in a
            #  reduced number of sections being retrieved.
            if not col_offset >= (3 * num_ats):
                raise HessError(tc, "Insufficient number of " + blockname +
                        " sections found",
                        "HESS File: " + hess_path)
            ## end if

            # Additional cross-check on number of rows imported
            if (rows_counter % (3*num_ats)) != 0:
                # Not the right number of rows; complain
                raise HessError(tc,
                            blockname + " row count mismatch",
                            "HESS File: " + hess_path)
            ## end if

            # Return the working matrix
            return workmtx

        ## end def parse_multiblock

        # Imports
        import re, numpy as np
        from .error import HessError
        from .utils import safe_cast as scast
        from .const import atom_num, atom_sym, PRM, DEF

        # Check if instantiated; complain if so
        if 'hess_path' in dir(self):
            raise HessError(HessError.OVERWRITE,
                    "Cannot overwrite contents of existing OrcaHess", "")
        ## end if

        # Retrieve the file target and store
        hess_path = kwargs['path']
        self.hess_path = hess_path

        # Open file, read contents, close stream
        with open(hess_path,'rU') as in_fl:
            self.in_str = in_fl.read()
        ## end with

        # Store the source string
        srcstr = "HESS File: {0}".format(hess_path)

        # Check to ensure all required data blocks are found
        if not OrcaHess.Pat.at_block.search(self.in_str):
            raise HessError(HessError.AT_BLOCK,
                    "Atom specification block not found", srcstr)
        ## end if
        if not OrcaHess.Pat.hess_block.search(self.in_str):
            raise HessError(HessError.HESS_BLOCK,
                    "Hessian block not found", srcstr)
        ## end if
        if not OrcaHess.Pat.freq_block.search(self.in_str):
            raise HessError(HessError.FREQ_BLOCK,
                    "Frequencies block (cm**-1 units) not found", srcstr)
        ## end if
        if not OrcaHess.Pat.modes_block.search(self.in_str):
            raise HessError(HessError.MODES_BLOCK,
                    "Normal modes block not found", srcstr)
        ## end if
        if not OrcaHess.Pat.energy.search(self.in_str):
            raise HessError(HessError.ENERGY,
                    "Energy value not found", srcstr)
        ## end if
        if not OrcaHess.Pat.temp.search(self.in_str):
            raise HessError(HessError.TEMP,
                    "'Actual temperature' value not found", srcstr)
        ## end if


        #=== Geometry spec block ===#
        # Store the block
        m_work = self.Pat.at_block.search(self.in_str)

        # Bring in the number of atoms
        self.num_ats = np.int_(m_work.group("num"))

        # Initialize the vectors of atomic symbols, masses, and the
        #  geometry
        self.atom_syms = []
        self.atom_masses = []
        self.geom = []

        # Parse the list of atoms
        for m in OrcaHess.Pat.at_line.finditer(m_work.group("block")):
            # Parse the element symbol or atomic number
            try:
                # See if it casts as an int
                num = scast(m.group('el'), np.int_)
            except (ValueError, TypeError):
                # Nope, must be letters. Check if valid by poking it into the
                #  dict. If not valid, should raise another error.
                num = atom_num[m.group('el').upper()]
            ## end try

            # If this point is reached, either it's cast ok to an int,
            #  or the string has been passed through the atom_num dict and
            #  is a valid element symbol. SO, convert back to atom_sym and
            #  store - this will also check for an invalid int entry
            self.atom_syms.append(atom_sym[num])

            # Parse the atomic weight; should be a simple append?
            self.atom_masses.append(scast(m.group("mass"), np.float_))

            # Build geometry as a 1-D list
            self.geom.extend([scast(m.group(gp), np.float_) for
                                                gp in ['c1', 'c2', 'c3']])
        ## next m

        # Double-check that the number of atoms retrieved matches the
        #  number indicated in the HESS file. Checks of atomic masses
        #  and geometry length are presumably redundant.
        if not len(self.atom_syms) == self.num_ats:
            raise HessError(HessError.AT_BLOCK, "Atomic symbol dimension " +
                    "mismatch with HESS atom specification", srcstr)
        ## end if

        # Convert geometry to numpy storage form
        self.geom = np.array(self.geom, dtype=np.float_)

        #=== Hessian and modes ===#
        # Pull the Hessian
        self.hess = parse_multiblock(self.in_str, self.Pat.hess_block,
                self.Pat.hess_sec, self.Pat.hess_line, self.num_ats,
                "Hessian", HessError.HESS_BLOCK)

        # Pull the modes
        self.modes = parse_multiblock(self.in_str, self.Pat.modes_block,
                self.Pat.modes_sec, self.Pat.modes_line, self.num_ats,
                "modes", HessError.MODES_BLOCK)

        # Extra check of 'dim' vs 'dim2' on modes
        m_work = self.Pat.modes_block.search(self.in_str)
        if scast(m_work.group("dim"), np.int_) != \
                            scast(m_work.group("dim2"), np.int_):
            raise HessError(HessError.MODES_BLOCK,
                    "Normal mode block dimension specification mismatch",
                    srcstr)
        ## end if

        #=== Frequencies ===#
        # Pull the block
        m_work = self.Pat.freq_block.search(self.in_str)

        # Check that number of frequencies indicated in the block matches
        #  that expected from the number of atoms
        if 3*self.num_ats != scast(m_work.group("num"), np.int_):
            raise HessError(HessError.FREQ_BLOCK,
                    "Count in frequencies block != 3 * number of atoms",
                    srcstr)
        ## end if

        # Retrieve the frequencies
        self.freqs = np.array(
                [scast(m.group("freq"), np.float_) for m in
                self.Pat.freq_line.finditer(m_work.group("block"))])

        # Proofread for proper size
        if not self.freqs.shape[0] == 3*self.num_ats:
            raise HessError(HessError.FREQ_BLOCK,
                    "Number of frequencies != 3 * number of atoms", srcstr)
        ## end if


        #=== Pull the single values that should always be present ===#
        # Store the reported energy
        self.energy = scast(self.Pat.energy.search(self.in_str).group("en"),
                                                                    np.float_)

        # Store the reported 'actual temperature'
        self.temp = scast(self.Pat.temp.search(self.in_str).group("temp"),
                                                                    np.float_)

        #=== Dipole derivatives ===#
        # Check if block found. Store None if not; otherwise import
        m_work = self.Pat.dipder_block.search(self.in_str)
        if m_work is None:
            self.dipders = None
        else:
            # Check that number of derivatives rows indicated in the block
            #  matches that expected from the number of atoms
            if 3*self.num_ats != scast(m_work.group("dim"), np.int_):
                raise HessError(HessError.DIPDER_BLOCK,
                        "Count in dipole derivatives block != 3 * # of atoms",
                        srcstr)
            ## end if

            # Retrieve the derivatives
            self.dipders = np.array(
                    [[scast(m.group("e{0}".format(i)), np.float_)
                        for i in range(3)]
                        for m in self.Pat.dipder_line.finditer(
                                                    m_work.group("block")) ])

            # Proofread for proper size. Don't have to proofread the width of
            #  three, since any row not containing three numerical values will
            #  result in the block getting truncated, per the Regex constrution.
            if not self.dipders.shape[0] == 3*self.num_ats:
                raise HessError(HessError.DIPDER_BLOCK,
                        "Number of dipole derivative rows != \
                        3 * number of atoms", srcstr)
            ## end if

            # If max-absolute element is too big, overwrite with None
            if np.max(np.abs(self.dipders)) > PRM.MAX_SANE_DIPDER:
                self.dipders = None
            ## end if
        ## end if


        #=== IR Spectrum ===#
        # If dipole derivs absent or munged, or if block missing, then skip
        m_work = self.Pat.ir_block.search(self.in_str)
        if self.dipders is None or m_work is None:
            self.ir_comps = None
            self.ir_mags = None
        else:
            # Complain if number of stated modes mismatches expectation
            if 3*self.num_ats != np.int_(m_work.group("dim")):
                raise HessError(HessError.IR_BLOCK,
                        "Count in IR spectrum block != 3 * # of atoms",
                        srcstr)
            ## end if

            # Pull the blocks
            self.ir_comps = np.array(
                    [[scast(m.group("e{0}".format(i)), np.float_)
                        for i in range(3)]
                        for m in self.Pat.ir_line.finditer(
                                                    m_work.group("block")) ])
            self.ir_mags = np.array(
                    [scast(m.group("mag"), np.float_)
                        for m in self.Pat.ir_line.finditer(
                                                    m_work.group("block")) ])

            # Confirm length of ir_mags conforms. Shouldn't need to check both,
            #  since they both rely equally on Pat.ir_line.finditer.
            if 3*self.num_ats != self.ir_mags.shape[0]:
                raise HessError(HessError.IR_BLOCK,
                        "Number of IR spectrum rows != \
                        3 * number of atoms", srcstr)
            ## end if

            # Confirm match of all frequencies with those reported separately
            if not np.allclose(
                    self.freqs,
                    np.array([np.float_(m.group('freq'))
                                for m in self.Pat.ir_line.finditer(
                                m_work.group("block") )]),
                    rtol=0,
                    atol=DEF.HESS_IR_MATCH_TOL):
                raise HessError(HessError.IR_BLOCK,
                        "Frequency mismatch between freq and IR blocks",
                        srcstr)
            ## end if
        ## end if


        #=== Polarizability Derivatives ===#
        # If block is missing, skip it
        m_work = self.Pat.polder_block.search(self.in_str)
        if m_work is None:
            self.polders = None
        else:
            # Check that number of derivatives rows indicated in the block
            #  matches that expected from the number of atoms
            if 3*self.num_ats != np.int_(m_work.group("dim")):
                raise HessError(HessError.POLDER_BLOCK,
                        "Count in polarizability derivatives block \
                        != 3 * # of atoms", srcstr)
            ## end if

            # Retrieve the derivatives
            self.polders = np.array(
                    [[scast(m.group("e{0}".format(i)), np.float_)
                        for i in range(6)]
                        for m in self.Pat.polder_line.finditer(
                                                    m_work.group("block")) ])

            # Proofread for proper size. Don't have to proofread the width of
            #  six, since any row not containing six numerical values will
            #  result in the block getting truncated.
            if not self.polders.shape[0] == 3*self.num_ats:
                raise HessError(HessError.POLDER_BLOCK,
                        "Number of polarizability derivative rows != \
                        3 * number of atoms", srcstr)
            ## end if
        ## end if


        #=== Raman Spectrum ===#
        # If polarizability derivs absent or munged, or if block missing,
        #  then skip
        m_work = self.Pat.raman_block.search(self.in_str)
        if self.polders is None or m_work is None:
            self.raman_acts = None
            self.raman_depols = None
        else:
            # Complain if number of stated modes mismatches expectation
            if 3*self.num_ats != np.int_(m_work.group("dim")):
                raise HessError(HessError.RAMAN_BLOCK,
                        "Count in Raman spectrum block != 3 * # of atoms",
                        srcstr)
            ## end if

            # Pull the blocks
            self.raman_acts = np.array(
                    [np.float_(m.group("act"))
                        for m in self.Pat.raman_line.finditer(
                                                    m_work.group("block")) ])
            self.raman_depols = np.array(
                    [np.float_(m.group("depol"))
                        for m in self.Pat.raman_line.finditer(
                                                    m_work.group("block")) ])

            # Confirm length of raman_acts conforms. Shouldn't need to check
            #  both, since they both rely equally on Pat.raman_line.finditer.
            if 3*self.num_ats != self.raman_acts.shape[0]:
                raise HessError(HessError.RAMAN_BLOCK,
                        "Number of Raman spectrum rows != \
                        3 * number of atoms", srcstr)
            ## end if

            # Confirm match of all frequencies with those reported separately
            if not np.allclose(
                    self.freqs,
                    np.array([np.float_(m.group('freq'))
                                for m in self.Pat.raman_line.finditer(
                                m_work.group("block")) ]),
                    rtol=0,
                    atol=DEF.HESS_IR_MATCH_TOL):
                raise HessError(HessError.RAMAN_BLOCK,
                        "Frequency mismatch between freq and Raman blocks",
                        srcstr)
            ## end if
        ## end if


        #=== Job list ===#
        # Check if block found. Store None if not; otherwise import
        m_work = self.Pat.jobs_block.search(self.in_str)
        if m_work is None:
            self.joblist = None
        else:
            # Check that number of joblist rows indicated in the block
            #  matches that expected from the number of atoms
            if 3*self.num_ats != np.int_(m_work.group("dim")):
                raise HessError(HessError.JOB_BLOCK,
                        "Count in job list block != 3 * # of atoms", srcstr)
            ## end if

            # Retrieve the job list
            self.joblist = np.array(
                    [[scast(m.group("e{0}".format(i)), np.float_)
                        for i in range(3)]
                        for m in self.Pat.jobs_line.finditer(
                                                    m_work.group("block")) ])

            # Proofread for proper size. Don't have to proofread the width of
            #  three, since any row not containing three numerical values will
            #  result in the block getting truncated.
            if not self.joblist.shape[0] == self.num_ats:
                raise HessError(HessError.JOB_BLOCK,
                        "Number of job list rows != number of atoms", srcstr)
            ## end if

            # Convert to boolean
            self.joblist = np.equal(self.joblist, np.ones_like(self.joblist))

        ## end if


        #=== Mass-weighted Hessian -- Eigenvalues ===#
        # Pull the block; continue only if block is present
        m_work = self.Pat.eigvals_block.search(self.in_str)
        if m_work is None:
            self.mwh_eigvals = None
        else:
            # Check that number of eigenvalues indicated in the block matches
            #  that expected from the number of atoms
            if 3*self.num_ats != np.int_(m_work.group("dim")):
                raise HessError(HessError.EIGVAL_BLOCK,
                        "Count in MWH eigenvalues block != 3 * number of atoms",
                        srcstr)
            ## end if

            # Retrieve the eigenvalues
            self.mwh_eigvals = np.array(
                    [np.float_(m.group("eig")) for m in
                    self.Pat.eigvals_line.finditer(m_work.group("block"))])

            # Proofread for proper size
            if not self.mwh_eigvals.shape[0] == 3*self.num_ats:
                raise HessError(HessError.EIGVAL_BLOCK,
                        "Number of MWH eigenvalues != 3 * number of atoms",
                        srcstr)
            ## end if

        ## end if


        #=== Mass-weighted Hessian -- Eigenvalues ===#
        # See if the block is there; import if so
        m_work = self.Pat.eigvecs_block.search(self.in_str)
        if m_work is None:
            self.mwh_eigvecs = None
        else:
            # Pull the eigenvectors
            self.mwh_eigvecs = \
                    parse_multiblock(self.in_str, self.Pat.eigvecs_block,
                    self.Pat.eigvecs_sec, self.Pat.eigvecs_line, self.num_ats,
                    "MWH eigenvectors", HessError.EIGVEC_BLOCK)

            # Extra check of 'dim' vs 'dim2' on modes
            if scast(m_work.group("dim"), np.int_) != \
                                scast(m_work.group("dim2"), np.int_):
                raise HessError(HessError.EIGVEC_BLOCK,
                        "MWH eigenvectors dimension specification mismatch",
                        srcstr)