How to use the qcelemental.constants function in qcelemental

def mol_tester(lbl, molstr, pg, sigma, refgeomang, isbohr=False, iso=False):
    symmol = qcdb.Molecule(molstr)
    if iso is not False:
        if isinstance(iso, int):
            iso = [iso]
        for at in iso:
            # mass needn't make sense for element, just breaking the symmetry
            symmol.set_mass(at, 2.014)
    symmol.update_geometry()
    symmol.axis_representation()
    assert compare_strings(pg, symmol.get_full_point_group(), pg + " point group: " + lbl)
    assert compare_integers(sigma, symmol.rotational_symmetry_number(), pg + " sigma")
    if isbohr:
        geom_now = symmol.full_geometry()
    else:
        geom_now = qcdb.mscale(symmol.full_geometry(), qcel.constants.bohr2angstroms)
    if refgeomang:
        assert compare_matrices(refgeomang, geom_now, 6, pg + " orientation")

def test_kabsch_identity():
    oco10 = qcel.molparse.from_string(soco10)
    oco12 = qcel.molparse.from_string(soco10)

    oco10_geom_au = oco10["qm"]["geom"].reshape((-1, 3)) / qcel.constants.bohr2angstroms
    oco12_geom_au = oco12["qm"]["geom"].reshape((-1, 3)) / qcel.constants.bohr2angstroms

    rmsd, rot, shift = qcel.molutil.kabsch_align(oco10_geom_au, oco12_geom_au)

    assert compare_values(0.0, rmsd, "identical")
    assert compare_values(np.identity(3), rot, "identity rotation matrix")
    assert compare_values(np.zeros(3), shift, "identical COM")

driver = query.pop("driver")
                qname = query.pop("name")
                data = self.get_records(
                    query.pop("method").upper(), include=["return_result"], merge=True, subset=subset, **query
                )
                new_data[qname] = data["return_result"]
                units[qname] = au_units[driver]
                query["name"] = qname
        else:
            for query in new_queries:
                query["native"] = True
            new_data, units = self._view.get_values(new_queries, subset)

        for query in new_queries:
            qname = query["name"]
            new_data[qname] *= constants.conversion_factor(units[qname], self.units)
            self._column_metadata[qname].update({"native": True, "units": self.units})

        self._update_cache(new_data)
        return self.df.loc[subset, names]

print(df.dtypes)
            print(df.head())

            df = df.merge(records, how="left", on="molecule")
            df.set_index("index", inplace=True)
            df.drop("molecule", axis=1, inplace=True)
            ret.append(df)

        if len(ret) == 1:
            retdf = ret[0]
        else:
            retdf = ret[0]
            for df in ret[1:]:
                retdf += df

        retdf[retdf.select_dtypes(include=['number']).columns] *= constants.conversion_factor('hartree', self.units)

        return retdf

P = np.identity(3 * nat)
    for irt in TRspace:
        P -= np.outer(irt, irt)
    text.append(mat_symm_info(P, lbl='total projector') + f' ({nrt})')

    # mass-weight & solve
    sqrtmmm = np.repeat(np.sqrt(mass), 3)
    sqrtmmminv = np.divide(1.0, sqrtmmm)
    mwhess = np.einsum('i,ij,j->ij', sqrtmmminv, nmwhess, sqrtmmminv)
    text.append(mat_symm_info(mwhess, lbl='mass-weighted Hessian') + ' (0)')

    pre_force_constant_au = np.linalg.eigvalsh(mwhess)

    idx = np.argsort(pre_force_constant_au)
    pre_force_constant_au = pre_force_constant_au[idx]
    uconv_cm_1 = (np.sqrt(qcel.constants.na * qcel.constants.hartree2J * 1.0e19) /
                  (2 * np.pi * qcel.constants.c * qcel.constants.bohr2angstroms))
    pre_frequency_cm_1 = np.lib.scimath.sqrt(pre_force_constant_au) * uconv_cm_1

    pre_lowfreq = np.where(np.real(pre_frequency_cm_1) < 100.0)[0]
    pre_lowfreq = np.append(pre_lowfreq, np.arange(nrt_expected))  # catch at least nrt modes
    for lf in set(pre_lowfreq):
        vlf = pre_frequency_cm_1[lf]
        if vlf.imag > vlf.real:
            text.append('  pre-proj  low-frequency mode: {:9.4f}i [cm^-1]'.format(vlf.real, vlf.imag))
        else:
            text.append('  pre-proj  low-frequency mode: {:9.4f}  [cm^-1]'.format(vlf.real, ''))
    text.append('  pre-proj  all modes:' + str(_format_omega(pre_frequency_cm_1, 4)))

    # project & solve
    mwhess_proj = np.dot(P.T, mwhess).dot(P)
    text.append(mat_symm_info(mwhess_proj, lbl='projected mass-weighted Hessian') + f' ({nrt})')

np.array
            1 by 3 of rotational constants or moments of inertia in units of `return_units`.

        Notes
        -----
        This used to return a list with inf values as None.

        """
        evals, evecs = diagonalize3x3symmat(self.inertia_tensor())
        evals = sorted(evals)
        evals = np.asarray(evals)

        im_amuA = qcel.constants.bohr2angstroms * qcel.constants.bohr2angstroms
        im_ghz = qcel.constants.h * qcel.constants.na * 1e14 / (8 * math.pi * math.pi * qcel.constants.bohr2angstroms * qcel.constants.bohr2angstroms)
        im_mhz = im_ghz * 1000.
        im_cm = im_ghz * 1.e7 / qcel.constants.c

        rc_moi = {}
        rc_moi['u a0^2'] = evals
        rc_moi['u A^2'] = evals * im_amuA
        with np.errstate(divide='ignore'):
            rc_moi['GHz'] = im_ghz / evals
            rc_moi['MHz'] = im_mhz / evals
            rc_moi['cm^-1'] = im_cm / evals

        fmt = """  {:12}    {a:3} {:16.8f}    {b:3} {:16.8f}    {c:3} {:16.8f}\n"""
        text = "        Moments of Inertia and Rotational Constants\n\n"
        text += fmt.format('[u a0^2]', a='I_A', b='I_B', c='I_C', *rc_moi['u a0^2'])
        text += fmt.format('[u A^2]', a='I_A', b='I_B', c='I_C', *rc_moi['u A^2'])
        text += fmt.format('[GHz]', a='A', b='B', c='C', *rc_moi['GHz'])
        text += fmt.format('[MHz]', a='A', b='B', c='C', *rc_moi['MHz'])
        text += fmt.format('[cm^-1]', a='A', b='B', c='C', *rc_moi['cm^-1'])

N = 0
    for i in range(self.natom()):
        if self.Z(i):
            N += 1

    # header
    text = '%s\n' % (self.name())
    text += '  Generated by xyz2mol\n\n'
    text += '%3i%3i  0  0  0  0  0  0  0  0999 V2000\n' % (N, len(bonds))

    # coordinates
    for i in range(self.natom()):
        x = self.x(i) * qcel.constants.bohr2angstroms
        y = self.y(i) * qcel.constants.bohr2angstroms
        z = self.z(i) * qcel.constants.bohr2angstroms
        if self.Z(i):
            text += ' %9.4f %9.4f %9.4f %-2s  0  0  0  0  0\n' % (x, y, z, self.symbol(i))

    # bonds
    for p in range(len(bonds)):
        text += '%3i%3i%3i' % (bonds[p][0] + 1, bonds[p][1] + 1, bonds[p][2])
        text += '  0  0  0\n'

    text += 'M  END\n'
    return text

#        for item in oriCoord:
#            print('       %16.8f %16.8f %16.8f' % (item[0], item[1], item[2]))
#
#    print '    <<<   [1] C4-GRD-GRAD   >>>'
#    if grdGrad is not None:
#        for item in grdGrad:
#            print('       %16.8f %16.8f %16.8f' % (item[0], item[1], item[2]))
#    print '    <<<   [2] C4-ORI-GRAD   >>>'
#    if oriGrad is not None:
#        for item in oriGrad:
#            print('       %16.8f %16.8f %16.8f' % (item[0], item[1], item[2]))

    retMol = None if p4Mol else grdMol

    if oriDip is not None:
        outPsivar['CURRENT DIPOLE X'] = str(oriDip[0] * qcel.constants.dipmom_au2debye)
        outPsivar['CURRENT DIPOLE Y'] = str(oriDip[1] * qcel.constants.dipmom_au2debye)
        outPsivar['CURRENT DIPOLE Z'] = str(oriDip[2] * qcel.constants.dipmom_au2debye)

    if oriGrad is not None:
        retGrad = oriGrad
    elif grdGrad is not None:
        retGrad = grdGrad
    else:
        retGrad = None

    if oriHess is not None:
        retHess = oriHess
    else:
        retHess = None

    return outPsivar, retGrad, retMol