How to use the qcelemental.molparse function in qcelemental

# r'\s+(?:-+)\s*' +
        # r'^\s+' + r'Z-matrix   Atomic            Coordinates (in bohr)' + r'\s*' +
        r'^\s+' + r'Symbol    Number           X              Y              Z' + r'\s*' +
        r'^\s+(?:-+)\s*' +
        r'((?:\s+[A-Z]+\s+[0-9]+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s*\n)+)' +
        r'^\s+(?:-+)\s*',
        outtext, re.MULTILINE)
    if mobj:
        print('matched geom')
        molxyz = '%d bohr\n\n' % len(mobj.group(1).splitlines())
        for line in mobj.group(1).splitlines():
            lline = line.split()
            molxyz += '%s %16s %16s %16s\n' % (lline[0], lline[-3], lline[-2], lline[-1])
        # Rather a dinky Molecule as no ghost, charge, or multiplicity
        psivar_coord = Molecule(validate=False,
                                **qcel.molparse.to_schema(qcel.molparse.from_string(molxyz,
                                                                                    dtype='xyz+',
                                                                                    fix_com=True,
                                                                                    fix_orientation=True)["qm"],
                                                          dtype=2))

    # Process atom geometry
    mobj = re.search(r'^\s+' + r'@GETXYZ-I,     1 atoms read from ZMAT.' + r'\s*$', outtext, re.MULTILINE)
    mobj2 = re.search(
        r'^([A-Z]+)#1' + r'\s+' + NUMBER + r'\s+' + NUMBER + r'\s+' + NUMBER + r'\s*$',
        outtext, re.MULTILINE)
    if mobj and mobj2:
        print('matched atom2')  # unsavory for when atom never printed except for basis file
        # Dinky Molecule
        molxyz = '1 bohr\n\n%s 0.0 0.0 0.0\n' % (mobj2.group(1))
        psivar_coord = Molecule(validate=False,
                                **qcel.molparse.to_schema(qcel.molparse.from_string(molxyz,

molxyz += '%s %16s %16s %16s\n' % (lline[-5], lline[-3], lline[-2], lline[-1])
                    # Jiyoung was collecting charge (-4)? see if this is ok for ghosts
                    # Tag    ,    X,        Y,        Z
                psivar_coord = Molecule(validate=False,
                                        **qcel.molparse.to_schema(qcel.molparse.from_string(
                                            molxyz, dtype='xyz+', fix_com=True, fix_orientation=True)["qm"],
                                                                  dtype=2))

            else:  # unit = a.u.
                molxyz = '%d au\n%d %d tag\n' % (len(mobj.group(2).splitlines()), out_charge, out_mult)
                for line in mobj.group(2).splitlines():
                    lline = line.split()
                    molxyz += '%s %16s %16s %16s\n' % (int(float(lline[-4])), lline[-3], lline[-2], lline[-1])
                    # Tag    ,    X,        Y,        Z
                psivar_coord = Molecule(validate=False,
                                        **qcel.molparse.to_schema(qcel.molparse.from_string(
                                            molxyz, dtype='xyz+', fix_com=True, fix_orientation=True)["qm"],
                                                                  dtype=2))

        # Process gradient
        mobj = re.search(
            r'^\s+' + r'.*' + r'ENERGY GRADIENTS' + r'\s*' + r'\s+' + r'\n' + r'^\s+' +
            r'atom               coordinates                        gradient' + r'\s*' + r'^\s+' +
            r'x          y          z           x          y          z' + r'\s*' +
            r'((?:\s+([1-9][0-9]*)+\s+([A-Z][a-x]*)+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s+[-+]?\d+\.\d+\s*\n)+)'
            + r'\s*$', outtext, re.MULTILINE)

        if mobj:
            logger.debug('matched molgrad')
            atoms = []
            psivar_grad = []
            for line in mobj.group(1).splitlines():

"""
    grd = grd.splitlines()
    Nat = int(grd[0].split()[0])
    molxyz = f"{Nat} bohr\n\n"

    grad = []
    for at in range(Nat):
        mline = grd[at + 1].split()
        el = "GH" if int(float(mline[0])) == 0 else qcel.periodictable.to_E(int(float(mline[0])))
        molxyz += "%s %16s %16s %16s\n" % (el, mline[-3], mline[-2], mline[-1])
        lline = grd[at + 1 + Nat].split()
        grad.append([float(lline[-3]), float(lline[-2]), float(lline[-1])])
    mol = Molecule(
        validate=False,
        **qcel.molparse.to_schema(
            qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
        ),
    )

    return mol, grad

"""
    grd = grd.splitlines()
    Nat = int(grd[0].split()[0])
    molxyz = f"{Nat} bohr\n\n"

    grad = []
    for at in range(Nat):
        mline = grd[at + 1].split()
        el = "GH" if int(float(mline[0])) == 0 else qcel.periodictable.to_E(int(float(mline[0])))
        molxyz += "%s %16s %16s %16s\n" % (el, mline[-3], mline[-2], mline[-1])
        lline = grd[at + 1 + Nat].split()
        grad.append([float(lline[-3]), float(lline[-2]), float(lline[-1])])
    mol = Molecule(
        validate=False,
        **qcel.molparse.to_schema(
            qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
        ),
    )

    return mol, grad

if mtd.startswith("mp2d-"):
            mtd = mtd[5:]

        if input_model.driver.derivative_int() > 1:
            raise InputError(f"""MP2D valid driver options are 'energy' and 'gradient', not {input_model.driver}""")

        # temp until actual options object
        input_model.extras["info"] = empirical_dispersion_resources.from_arrays(
            name_hint=mtd,
            level_hint=input_model.keywords.get("level_hint", None),
            param_tweaks=input_model.keywords.get("params_tweaks", None),
            dashcoeff_supplement=input_model.keywords.get("dashcoeff_supplement", None),
        )

        # Need 'real' field later and that's only guaranteed for molrec
        molrec = qcel.molparse.from_schema(input_model.molecule.dict())
        xyz = qcel.molparse.to_string(molrec, dtype="xyz", units="Angstrom", ghost_format="")
        infiles = {"mp2d_geometry": xyz}
        # jobrec['molecule']['real'] = molrec['real']

        # env = {
        #    'HOME': os.environ.get('HOME'),
        #    'PATH': os.environ.get('PATH'),
        #    #'PATH': os.pathsep.join([os.path.abspath(x) for x in os.environ.get('PSIPATH', '').split(os.pathsep) if x != '']) + \
        #    #        os.pathsep + os.environ.get('PATH'),
        #    #'LD_LIBRARY_PATH': os.environ.get('LD_LIBRARY_PATH'),
        # }

        command = ["mp2d", "mp2d_geometry"]
        command.extend(
            """--TT_a1={a1} --TT_a2={a2} --rcut={rcut} --w={w} --s8={s8}""".format(
                **input_model.extras["info"]["dashparams"]

#    #'PATH': os.pathsep.join([os.path.abspath(x) for x in os.environ.get('PSIPATH', '').split(os.pathsep) if x != '']) + \
        #    #        os.pathsep + os.environ.get('PATH'),
        #    #'LD_LIBRARY_PATH': os.environ.get('LD_LIBRARY_PATH'),
        # }

        command = ["dftd3", "dftd3_geometry.xyz"]
        if input_model.driver == "gradient":
            command.append("-grad")
        if input_model.extras["info"]["dashlevel"] == "atmgr":
            command.append("-abc")

        infiles = {
            ".dftd3par.local": dftd3_coeff_formatter(
                input_model.extras["info"]["dashlevel"], input_model.extras["info"]["dashparams"]
            ),
            "dftd3_geometry.xyz": qcel.molparse.to_string(molrec, dtype="xyz", units="Angstrom", ghost_format=""),
        }

        return {
            "command": command,
            "infiles": infiles,
            "outfiles": ["dftd3_gradient", "dftd3_abc_gradient"],
            "scratch_directory": config.scratch_directory,
            "input_result": input_model.copy(deep=True),
            "blocking_files": [os.path.join(pathlib.Path.home(), ".dftd3par." + socket.gethostname())],
        }

elez=celez[fr],
                    units='Bohr',
                    fix_com=True,
                    fix_orientation=True) for fr in frag_pattern
            ]
            if isinstance(self, Molecule):
                ret_mols = [Molecule.from_dict(molrec) for molrec in molrecs]
            else:
                from psi4 import core
                ret_mols = [core.Molecule.from_dict(molrec) for molrec in molrecs]
            outputs.append(ret_mols)

        if return_molecule:
            dcontig = qcel.molparse.contiguize_from_fragment_pattern(
                frag_pattern, geom=cgeom, elez=celez, elem=celem, mass=cmass)
            molrec = qcel.molparse.from_arrays(
                geom=dcontig['geom'],
                mass=dcontig['mass'],
                elem=dcontig['elem'],
                elez=dcontig['elez'],
                units='Bohr',
                molecular_charge=self.molecular_charge(),
                # molecular_multiplicity may not be conservable upon fragmentation
                #   potentially could do two passes and try to preserve it
                fix_com=self.com_fixed(),
                fix_orientation=self.orientation_fixed(),
                fix_symmetry=(None if self.symmetry_from_input() == '' else self.symmetry_from_input()),
                fragment_separators=dcontig['fragment_separators'])
            if isinstance(self, Molecule):
                ret_mol = Molecule.from_dict(molrec)
            else:
                from psi4 import core

missing_enabled_return_efp='none',
                 missing_enabled_return='error',
                 tooclose=0.1,
                 zero_ghost_fragments=False,
                 nonphysical=False,
                 mtol=1.e-3,
                 verbose=1):
        """Initialize Molecule object from LibmintsMolecule"""
        super(Molecule, self).__init__()

        if molinit is not None or geom is not None:
            if isinstance(molinit, dict):
                molrec = molinit

            elif isinstance(molinit, str):
                compound_molrec = qcel.molparse.from_string(
                    molstr=molinit,
                    dtype=dtype,
                    name=name,
                    fix_com=fix_com,
                    fix_orientation=fix_orientation,
                    fix_symmetry=fix_symmetry,
                    return_processed=False,
                    enable_qm=enable_qm,
                    enable_efp=enable_efp,
                    missing_enabled_return_qm=missing_enabled_return_qm,
                    missing_enabled_return_efp=missing_enabled_return_efp,
                    verbose=verbose)
                molrec = compound_molrec['qm']

            elif molinit is None and geom is not None:
                molrec = qcel.molparse.from_arrays(