From JMS to semileptonic Fierz basis for Class V.
    `ddll` should be of the form 'sbl_enu_tau', 'dbl_munu_e' etc.
11237
def JMS_to_Fierz_lep(C, ddll):
    """"""
    if ddll[:2] == 'uc':
        s = uflav[ddll[0]]
        b = uflav[ddll[1]]
        q = 'u'
    else:
        s = dflav[ddll[0]]
        b = dflav[ddll[1]]
        q = 'd'
    l = lflav[ddll[4:ddll.find('n')]]
    lp = lflav[ddll[ddll.find('_',5)+1:len(ddll)]]
    ind = ddll.replace('l_','').replace('nu_','')
    return {
        'F' + ind + '9' : C["V" + q + "eLR"][s, b, l, lp] / 2
                          + C["Ve" + q + "LL"][l, lp, s, b] / 2,
        'F' + ind + '10' : C["V" + q + "eLR"][s, b, l, lp] / 2
                           - C["Ve" + q + "LL"][l, lp, s, b] / 2,
        'F' + ind + 'S' : C["Se" + q + "RL"][lp, l, b, s].conj() / 2
                          + C["Se" + q + "RR"][l, lp, s, b] / 2,
        'F' + ind + 'P' : - C["Se" + q + "RL"][lp, l, b, s].conj() / 2
                          + C["Se" + q + "RR"][l, lp, s, b] / 2,
        'F' + ind + 'T' : C["Te" + q + "RR"][l, lp, s, b] / 2
                          + C["Te" + q + "RR"][lp, l, b, s].conj() / 2,
        'F' + ind + 'T5' : C["Te" + q + "RR"][l, lp, s, b] / 2
                           - C["Te" + q + "RR"][lp, l, b, s].conj() / 2,
        'F' + ind + '9p' : C["Ve" + q + "LR"][l, lp, s, b] / 2
                           + C["Ve" + q + "RR"][l, lp, s, b] / 2,
        'F' + ind + '10p' : -C["Ve" + q + "LR"][l, lp, s, b] / 2
                            + C["Ve" + q + "RR"][l, lp, s, b] / 2,
        'F' + ind + 'Sp' : C["Se" + q + "RL"][l, lp, s, b] / 2
                           + C["Se" + q + "RR"][lp, l, b, s].conj() / 2,
        'F' + ind + 'Pp' : C["Se" + q + "RL"][l, lp, s, b] / 2
                            - C["Se" + q + "RR"][lp, l, b, s].conj() / 2,
        }
730872
def _JMS_to_Fierz_III_IV_V(C, qqqq):
    """"""
    #case dduu
    classIII = ['sbuc', 'sbcu', 'dbuc', 'dbcu', 'dsuc', 'dscu']
    classVdduu = ['sbuu' , 'dbuu', 'dsuu', 'sbcc' , 'dbcc', 'dscc']
    if qqqq in classIII + classVdduu:
        f1 = dflav[qqqq[0]]
        f2 = dflav[qqqq[1]]
        f3 = uflav[qqqq[2]]
        f4 = uflav[qqqq[3]]
        return {
            'F' + qqqq + '1' : C["V1udLL"][f3, f4, f1, f2]
                                - C["V8udLL"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '2' : C["V8udLL"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '3' : C["V1duLR"][f1, f2, f3, f4]
                                - C["V8duLR"][f1, f2, f3, f4] / (2 * Nc),
            'F' + qqqq + '4' : C["V8duLR"][f1, f2, f3, f4] / 2,
            'F' + qqqq + '5' : C["S1udRR"][f3, f4, f1, f2]
                                - C["S8udduRR"][f3, f2, f1, f4] / 4
                                - C["S8udRR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '6' : -C["S1udduRR"][f3, f2, f1, f4] / 2
                                + C["S8udduRR"][f3, f2, f1, f4] /(4 * Nc)
                                + C["S8udRR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '7' : -C["V8udduLR"][f4, f1, f2, f3].conj(),
            'F' + qqqq + '8' : -2 * C["V1udduLR"][f4, f1, f2, f3].conj()
                                + C["V8udduLR"][f4, f1, f2, f3].conj() / Nc,
            'F' + qqqq + '9' : -C["S8udduRR"][f3, f2, f1, f4] / 16,
            'F' + qqqq + '10' : -C["S1udduRR"][f3, f2, f1, f4] / 8
                                + C["S8udduRR"][f3, f2, f1, f4] / (16 * Nc),
            'F' + qqqq + '1p' : C["V1udRR"][f3, f4, f1, f2]
                                - C["V8udRR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '2p' : C["V8udRR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '3p' : C["V1udLR"][f3, f4, f1, f2]
                                - C["V8udLR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '4p' : C["V8udLR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '5p' : C["S1udRR"][f4, f3, f2, f1].conj() -
                                C["S8udduRR"][f4, f1, f2, f3].conj() / 4
                                - C["S8udRR"][f4, f3, f2, f1].conj() / (2 * Nc),
            'F' + qqqq + '6p' : -C["S1udduRR"][f4, f1, f2, f3].conj() / 2 +
                                C["S8udduRR"][f4, f1, f2, f3].conj()/(4 * Nc)
                                + C["S8udRR"][f4, f3, f2, f1].conj() / 2,
            'F' + qqqq + '7p' : -C["V8udduLR"][f3, f2, f1, f4],
            'F' + qqqq + '8p' : - 2 * C["V1udduLR"][f3, f2, f1, f4]
                                + C["V8udduLR"][f3, f2, f1, f4] / Nc,
            'F' + qqqq + '9p' : -C["S8udduRR"][f4, f1, f2, f3].conj() / 16,
            'F' + qqqq + '10p' : -C["S1udduRR"][f4, f1, f2, f3].conj() / 8
                                + C["S8udduRR"][f4, f1, f2, f3].conj() / 16 / Nc
                            }
    classVuudd = ['ucdd', 'ucss', 'ucbb']
    if qqqq in classVuudd:
        f3 = uflav[qqqq[0]]
        f4 = uflav[qqqq[1]]
        f1 = dflav[qqqq[2]]
        f2 = dflav[qqqq[3]]
        return {
            'F' + qqqq + '1' : C["V1udLL"][f3, f4, f1, f2]
                                - C["V8udLL"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '2' : C["V8udLL"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '3p' : C["V1duLR"][f1, f2, f3, f4]
                                - C["V8duLR"][f1, f2, f3, f4] / (2 * Nc),
            'F' + qqqq + '4p' : C["V8duLR"][f1, f2, f3, f4] / 2,
            'F' + qqqq + '5' : C["S1udRR"][f3, f4, f1, f2]
                                - C["S8udduRR"][f3, f2, f1, f4] / 4
                                - C["S8udRR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '6' : -C["S1udduRR"][f3, f2, f1, f4] / 2
                                + C["S8udduRR"][f3, f2, f1, f4] /(4 * Nc)
                                + C["S8udRR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '7p' : -C["V8udduLR"][f4, f1, f2, f3].conj(),
            'F' + qqqq + '8p' : -2 * C["V1udduLR"][f4, f1, f2, f3].conj()
                                + C["V8udduLR"][f4, f1, f2, f3].conj() / Nc,
            'F' + qqqq + '9' : -C["S8udduRR"][f3, f2, f1, f4] / 16,
            'F' + qqqq + '10' : -C["S1udduRR"][f3, f2, f1, f4] / 8
                                + C["S8udduRR"][f3, f2, f1, f4] / (16 * Nc),
            'F' + qqqq + '1p' : C["V1udRR"][f3, f4, f1, f2]
                                - C["V8udRR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '2p' : C["V8udRR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '3' : C["V1udLR"][f3, f4, f1, f2]
                                - C["V8udLR"][f3, f4, f1, f2] / (2 * Nc),
            'F' + qqqq + '4' : C["V8udLR"][f3, f4, f1, f2] / 2,
            'F' + qqqq + '5p' : C["S1udRR"][f4, f3, f2, f1].conj() -
                                C["S8udduRR"][f4, f1, f2, f3].conj() / 4
                                - C["S8udRR"][f4, f3, f2, f1].conj() / (2 * Nc),
            'F' + qqqq + '6p' : -C["S1udduRR"][f4, f1, f2, f3].conj() / 2 +
                                C["S8udduRR"][f4, f1, f2, f3].conj()/(4 * Nc)
                                + C["S8udRR"][f4, f3, f2, f1].conj() / 2,
            'F' + qqqq + '7' : -C["V8udduLR"][f3, f2, f1, f4],
            'F' + qqqq + '8' : - 2 * C["V1udduLR"][f3, f2, f1, f4]
                                + C["V8udduLR"][f3, f2, f1, f4] / Nc,
            'F' + qqqq + '9p' : -C["S8udduRR"][f4, f1, f2, f3].conj() / 16,
            'F' + qqqq + '10p' : -C["S1udduRR"][f4, f1, f2, f3].conj() / 8
                                + C["S8udduRR"][f4, f1, f2, f3].conj() / 16 / Nc
                            }
    #case dddd
    classIV = ['sbsd', 'dbds', 'bsbd']
    classVdddd = ['sbss', 'dbdd', 'dsdd', 'sbbb', 'dbbb', 'dsss']
    classVddddind = ['sbdd', 'dsbb', 'dbss']
    if qqqq in classIV + classVdddd + classVddddind:
        f1 = dflav[qqqq[0]]
        f2 = dflav[qqqq[1]]
        f3 = dflav[qqqq[2]]
        f4 = dflav[qqqq[3]]
        return {
                'F'+ qqqq +'1' : C["VddLL"][f3, f4, f1, f2],
                 'F'+ qqqq +'2' : C["VddLL"][f1, f4, f3, f2],
                 'F'+ qqqq +'3' : C["V1ddLR"][f1, f2, f3, f4]
                                  - C["V8ddLR"][f1, f2, f3, f4]/(2 * Nc),
                 'F'+ qqqq +'4' : C["V8ddLR"][f1, f2, f3, f4] / 2,
                 'F'+ qqqq +'5' : C["S1ddRR"][f3, f4, f1, f2]
                                  - C["S8ddRR"][f3, f2, f1,f4] / 4
                                  - C["S8ddRR"][f3, f4, f1, f2] / (2 * Nc),
                 'F'+ qqqq +'6' : -C["S1ddRR"][f1, f4, f3, f2] / 2
                                  + C["S8ddRR"][f3, f2, f1, f4] / (4 * Nc)
                                  + C["S8ddRR"][f3, f4, f1, f2] / 2,
                 'F'+ qqqq +'7' : -C["V8ddLR"][f1, f4, f3, f2],
                 'F'+ qqqq +'8' : -2 * C["V1ddLR"][f1, f4, f3, f2]
                                  + C["V8ddLR"][f1, f4, f3, f2] / Nc,
                 'F'+ qqqq +'9' : -C["S8ddRR"][f3, f2, f1, f4] / 16,
                 'F'+ qqqq +'10' : -C["S1ddRR"][f1, f4, f3, f2] / 8
                                   + C["S8ddRR"][f3, f2, f1, f4] / (16 * Nc),
                 'F'+ qqqq +'1p' : C["VddRR"][f3, f4, f1, f2],
                 'F'+ qqqq +'2p' : C["VddRR"][f1, f4, f3, f2],
                 'F'+ qqqq +'3p' : C["V1ddLR"][f3, f4, f1, f2]
                                   - C["V8ddLR"][f3, f4, f1,f2] / (2 * Nc),
                 'F'+ qqqq +'4p' : C["V8ddLR"][f3, f4, f1, f2] / 2,
                 'F'+ qqqq +'5p' : C["S1ddRR"][f4, f3, f2, f1].conj() -
                                   C["S8ddRR"][f4, f1, f2, f3].conj() / 4
                                   -C["S8ddRR"][f4, f3, f2, f1].conj() / 2 / Nc,
                 'F'+ qqqq +'6p' : -C["S1ddRR"][f4, f1, f2, f3].conj() / 2 +
                                    C["S8ddRR"][f4, f1, f2, f3].conj() / 4 / Nc
                                    + C["S8ddRR"][f4, f3, f2, f1].conj() / 2,
                 'F'+ qqqq +'7p' : -C["V8ddLR"][f3, f2, f1, f4],
                 'F'+ qqqq +'8p' : -2 * C["V1ddLR"][f3, f2, f1, f4]
                                    + C["V8ddLR"][f3, f2, f1, f4] / Nc,
                 'F'+ qqqq +'9p' : -C["S8ddRR"][f4, f1, f2, f3].conj() / 16,
                 'F'+ qqqq +'10p' : -C["S1ddRR"][f4, f1, f2, f3].conj() / 8 +
                                    C["S8ddRR"][f4, f1, f2, f3].conj() / 16 / Nc
                                    }
    #case uuuu
    classVuuuu = ['ucuu', 'cucc', 'cuuu', 'uccc']
    if qqqq in classVuuuu:
        f1 = uflav[qqqq[0]]
        f2 = uflav[qqqq[1]]
        f3 = uflav[qqqq[2]]
        f4 = uflav[qqqq[3]]
        return {
                'F' + qqqq + '1' : C["VuuLL"][f3, f4, f1, f2],
                'F' + qqqq + '2' : C["VuuLL"][f1, f4, f3, f2],
                'F' + qqqq + '3' : C["V1uuLR"][f1, f2, f3, f4]
                                   - C["V8uuLR"][f1, f2, f3, f4] / (2 * Nc),
                'F' + qqqq + '4' : C["V8uuLR"][f1, f2, f3, f4] / 2,
                'F' + qqqq + '5' : C["S1uuRR"][f3, f4, f1, f2]
                                    - C["S8uuRR"][f3, f2, f1, f4] / 4
                                    - C["S8uuRR"][f3, f4, f1, f2] / (2 * Nc),
                'F' + qqqq + '6' : -C["S1uuRR"][f1, f4, f3, f2] / 2
                                    + C["S8uuRR"][f3, f2, f1, f4] / (4 * Nc)
                                    + C["S8uuRR"][f3, f4, f1, f2] / 2,
                'F' + qqqq + '7' : -C["V8uuLR"][f1, f4, f3, f2],
                'F' + qqqq + '8' : -2 * C["V1uuLR"][f1, f4, f3, f2]
                                    + C["V8uuLR"][f1, f4, f3, f2] / Nc,
                'F' + qqqq + '9' : -C["S8uuRR"][f3, f2, f1, f4] / 16,
                'F' + qqqq + '10' : -C["S1uuRR"][f1, f4, f3, f2] / 8
                                    + C["S8uuRR"][f3, f2, f1, f4] / (16 * Nc),
                'F'+ qqqq + '1p': C["VuuRR"][f3, f4, f1, f2],
                'F' + qqqq + '2p': C["VuuRR"][f1, f3, f4, f2],
                'F' + qqqq + '3p' : C["V1uuLR"][f3, f4, f1, f2]
                                    - C["V8uuLR"][f3, f4, f1,f2] / (2 * Nc),
                'F' + qqqq + '4p' : C["V8uuLR"][f3, f4, f1, f2] / 2,
                'F' + qqqq + '5p' : C["S1uuRR"][f4, f3, f2, f1].conj() -
                                    C["S8uuRR"][f4, f1, f2, f3].conj() / 4 -
                                    C["S8uuRR"][f4, f3, f2, f1].conj() / 2 / Nc,
                'F' + qqqq + '6p' : -C["S1uuRR"][f4, f1, f2, f3].conj() / 2 +
                                    C["S8uuRR"][f4, f1, f2, f3].conj() / 4 / Nc
                                    + C["S8uuRR"][f4, f3, f2, f1].conj() / 2,
                'F' + qqqq + '7p' : -C["V8uuLR"][f3, f2, f1, f4],
                'F' + qqqq + '8p' : -2 * C["V1uuLR"][f3, f2, f1, f4]
                                    + C["V8uuLR"][f3, f2, f1, f4] / Nc,
                'F' + qqqq + '9p' : -C["S8uuRR"][f4, f1, f2, f3].conj() / 16,
                'F' + qqqq + '10p' : -C["S1uuRR"][f4, f1, f2, f3].conj() / 8 +
                                    C["S8uuRR"][f4, f1, f2, f3].conj() / 16 / Nc
                                    }
    else:
        raise ValueError("Case not implemented: {}".format(qqqq))
725003
def as_dict(self):
        
        d = {"vasp_version": self.vasp_version,
             "has_vasp_completed": self.converged,
             "nsites": len(self.final_structure)}
        comp = self.final_structure.composition
        d["unit_cell_formula"] = comp.as_dict()
        d["reduced_cell_formula"] = Composition(comp.reduced_formula).as_dict()
        d["pretty_formula"] = comp.reduced_formula
        symbols = [s.split()[1] for s in self.potcar_symbols]
        symbols = [re.split(r"_", s)[0] for s in symbols]
        d["is_hubbard"] = self.is_hubbard
        d["hubbards"] = self.hubbards

        unique_symbols = sorted(list(set(self.atomic_symbols)))
        d["elements"] = unique_symbols
        d["nelements"] = len(unique_symbols)

        d["run_type"] = self.run_type

        vin = {"incar": {k: v for k, v in self.incar.items()},
               "crystal": self.initial_structure.as_dict(),
               "kpoints": self.kpoints.as_dict()}
        actual_kpts = [{"abc": list(self.actual_kpoints[i]),
                        "weight": self.actual_kpoints_weights[i]}
                       for i in range(len(self.actual_kpoints))]
        vin["kpoints"]["actual_points"] = actual_kpts
        vin["nkpoints"] = len(actual_kpts)
        vin["potcar"] = [s.split(" ")[1] for s in self.potcar_symbols]
        vin["potcar_spec"] = self.potcar_spec
        vin["potcar_type"] = [s.split(" ")[0] for s in self.potcar_symbols]
        vin["parameters"] = {k: v for k, v in self.parameters.items()}
        vin["lattice_rec"] = self.lattice_rec.as_dict()
        d["input"] = vin

        nsites = len(self.final_structure)

        try:
            vout = {"ionic_steps": self.ionic_steps,
                    "final_energy": self.final_energy,
                    "final_energy_per_atom": self.final_energy / nsites,
                    "crystal": self.final_structure.as_dict(),
                    "efermi": self.efermi}
        except (ArithmeticError, TypeError):
            vout = {"ionic_steps": self.ionic_steps,
                    "final_energy": self.final_energy,
                    "final_energy_per_atom": None,
                    "crystal": self.final_structure.as_dict(),
                    "efermi": self.efermi}

        if self.eigenvalues:
            eigen = {str(spin): v.tolist()
                     for spin, v in self.eigenvalues.items()}
            vout["eigenvalues"] = eigen
            (gap, cbm, vbm, is_direct) = self.eigenvalue_band_properties
            vout.update(dict(bandgap=gap, cbm=cbm, vbm=vbm,
                             is_gap_direct=is_direct))

            if self.projected_eigenvalues:
                vout['projected_eigenvalues'] = {
                    str(spin): v.tolist()
                    for spin, v in self.projected_eigenvalues.items()}

        vout['epsilon_static'] = self.epsilon_static
        vout['epsilon_static_wolfe'] = self.epsilon_static_wolfe
        vout['epsilon_ionic'] = self.epsilon_ionic
        d['output'] = vout
        return jsanitize(d, strict=True)
456387
def sm_dict2lha(d):
    """"""
    blocks = OrderedDict([
        ('GAUGE', {'values': [[1, d['g'].real], [2, d['gp'].real], [3, d['gs'].real]]}),
        ('SCALAR', {'values': [[1, d['Lambda'].real], [2, d['m2'].real]]}),
        ('GU', {'values': matrix2lha(d['Gu'].real)}),
        ('IMGU', {'values': matrix2lha(d['Gu'].imag)}),
        ('GD', {'values': matrix2lha(d['Gd'].real)}),
        ('IMGD', {'values': matrix2lha(d['Gd'].imag)}),
        ('GE', {'values': matrix2lha(d['Ge'].real)}),
        ('IMGE', {'values': matrix2lha(d['Ge'].imag)}),
        ('THETA', {'values': [[1, d['Theta'].real], [2, d['Thetap'].real], [3, d['Thetas'].real]]}),
        ])
    return {'BLOCK': blocks}
313524
def cublasDgemv(handle, trans, m, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasDgemv_v2(handle,
                                       _CUBLAS_OP[trans], m, n,
                                       ctypes.byref(ctypes.c_double(alpha)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(ctypes.c_double(beta)),
                                       int(y), incy)
    cublasCheckStatus(status)
1150953
def _Bern_to_JMS_II(C, udlnu):
    """"""
    u = uflav[udlnu[0]]
    d = dflav[udlnu[1]]
    l = lflav[udlnu[4:udlnu.find('n')]]
    lp = lflav[udlnu[udlnu.find('_',5)+1:len(udlnu)]]
    ind = udlnu[0]+udlnu[1]+udlnu[4:udlnu.find('n')]+udlnu[udlnu.find('_',5)+1
                                                                    :len(udlnu)]
    return {
        "VnueduLL_{}{}{}{}".format(lp + 1, l + 1, d + 1, u + 1): C['1' + ind].conjugate(),
        "SnueduRL_{}{}{}{}".format(lp + 1, l + 1, d + 1, u + 1): C['5' + ind].conjugate(),
        "VnueduLR_{}{}{}{}".format(lp + 1, l + 1, d + 1, u + 1): C['1p' + ind].conjugate(),
        "SnueduRR_{}{}{}{}".format(lp + 1, l + 1, d + 1, u + 1): C['5p' + ind].conjugate(),
        "TnueduRR_{}{}{}{}".format(lp + 1, l + 1, d + 1, u + 1): C['7p' + ind].conjugate()
        }
778291
def calc_freefree_snu_ujy(ne, t, width, elongation, dist, ghz):
    
    hz = ghz * 1e9
    eta = calc_freefree_eta(ne, t, hz)
    kappa = calc_freefree_kappa(ne, t, hz)
    snu = calc_snu(eta, kappa, width, elongation, dist)
    ujy = snu * cgs.jypercgs * 1e6
    return ujy
273081
def vbem_aBl(self):
        
        vbem_aBl = np.zeros((self.T, self.num_states))
        ids, dds, eds = self.init_dynamics_distns, self.dynamics_distns, self.emission_distns

        for k, (id, dd) in enumerate(zip(ids, dds)):
            vbem_aBl[0, k] = expected_gaussian_logprob(id.mu, id.sigma, self.E_init_stats)
            vbem_aBl[:-1, k] += expected_regression_log_prob(dd, self.E_dynamics_stats)

        if self.single_emission:
            ed = self.emission_distns[0]
            vbem_aBl += expected_regression_log_prob(ed, self.E_emission_stats)[:,None]
        else:
            for k, ed in enumerate(self.emission_distns):
                vbem_aBl[:, k] += expected_regression_log_prob(ed, self.E_emission_stats)

        vbem_aBl[np.isnan(vbem_aBl).any(1)] = 0.
        return vbem_aBl
282039
def read_turbomole(basis_lines, fname):
    '''
    '''

    skipchars = '*#$'
    basis_lines = [l for l in basis_lines if l and not l[0] in skipchars]

    bs_data = create_skel('component')

    i = 0
    while i < len(basis_lines):
        line = basis_lines[i]
        elementsym = line.split()[0]

        element_Z = lut.element_Z_from_sym(elementsym)
        element_Z = str(element_Z)

        if not element_Z in bs_data['elements']:
            bs_data['elements'][element_Z] = {}

        element_data = bs_data['elements'][element_Z]

        if "ecp" in line.lower():
            if not 'ecp_potentials' in element_data:
                element_data['ecp_potentials'] = []

            i += 1
            line = basis_lines[i]

            lsplt = line.split('=')
            maxam = int(lsplt[2])
            n_elec = int(lsplt[1].split()[0])

            amlist = [maxam]
            amlist.extend(list(range(0, maxam)))

            i += 1
            for shell_am in amlist:
                shell_am2 = lut.amchar_to_int(basis_lines[i][0])[0]
                if shell_am2 != shell_am:
                    raise RuntimeError("AM not in expected order?")

                i += 1

                ecp_shell = {
                    'ecp_type': 'scalar_ecp',
                    'angular_momentum': [shell_am],
                }
                ecp_exponents = []
                ecp_rexponents = []
                ecp_coefficients = []

                while i < len(basis_lines) and basis_lines[i][0].isalpha() is False:
                    lsplt = basis_lines[i].split()
                    ecp_exponents.append(lsplt[2])
                    ecp_rexponents.append(int(lsplt[1]))
                    ecp_coefficients.append(lsplt[0])
                    i += 1

                ecp_shell['r_exponents'] = ecp_rexponents
                ecp_shell['gaussian_exponents'] = ecp_exponents
                ecp_shell['coefficients'] = [ecp_coefficients]
                element_data['ecp_potentials'].append(ecp_shell)

            element_data['ecp_electrons'] = n_elec

        else:
            if not 'electron_shells' in element_data:
                element_data['electron_shells'] = []

            i += 1
            while i < len(basis_lines) and basis_lines[i][0].isalpha() == False:
                lsplt = basis_lines[i].split()
                shell_am = lut.amchar_to_int(lsplt[1])
                nprim = int(lsplt[0])

                if max(shell_am) <= 1:
                    func_type = 'gto'
                else:
                    func_type = 'gto_spherical'

                shell = {
                    'function_type': func_type,
                    'region': '',
                    'angular_momentum': shell_am
                }

                exponents = []
                coefficients = []

                i += 1
                for j in range(nprim):
                    line = basis_lines[i].replace('D', 'E')
                    line = line.replace('d', 'E')
                    lsplt = line.split()
                    exponents.append(lsplt[0])
                    coefficients.append(lsplt[1:])
                    i += 1

                shell['exponents'] = exponents

                # We need to transpose the coefficient matrix
                # (we store a matrix with primitives being the column index and
                # general contraction being the row index)
                shell['coefficients'] = list(map(list, zip(*coefficients)))

                element_data['electron_shells'].append(shell)

    return bs_data
1208923
def warsaw_to_warsawmass(C, parameters=None, sectors=None):
    
    p = default_parameters.copy()
    if parameters is not None:
        # if parameters are passed in, overwrite the default values
        p.update(parameters)
    # start out with a 1:1 copy
    C_out = C.copy()
    # rotate left-handed up-type quark fields in uL-uR operator WCs
    C_rotate_u = ['uphi', 'uG', 'uW', 'uB']
    for name in C_rotate_u:
        _array = smeft_toarray(name, C)
        V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
        UuL = V.conj().T
        _array = UuL.conj().T @ _array
        _dict = smeft_fromarray(name, _array)
        C_out.update(_dict)
    # diagonalize dimension-5 Weinberg operator
    _array = smeft_toarray('llphiphi', C)
    _array = np.diag(ckmutil.diag.msvd(_array)[1])
    _dict = smeft_fromarray('llphiphi', _array)
    C_out.update(_dict)
    return C_out
241475
def d2Ibr_dV2(Ybr, V, lam):
    
    nb = len(V)
    diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = spdiag(mul(-(Ybr.T * lam), V))
    Hva = -1j * Haa * diaginvVm
    Hav = Hva
    Hvv = spmatrix([], [], [], (nb, nb))

    return Haa, Hav, Hva, Hvv
489107
def blt(f: List[SYM], x: List[SYM]) -> Dict[str, Any]:
    
    J = ca.jacobian(f, x)
    nblock, rowperm, colperm, rowblock, colblock, coarserow, coarsecol = J.sparsity().btf()
    return {
        'J': J,
        'nblock': nblock,
        'rowperm': rowperm,
        'colperm': colperm,
        'rowblock': rowblock,
        'colblock': colblock,
        'coarserow': coarserow,
        'coarsecol': coarsecol
    }
251937
def Fierz_to_EOS_chrom(C, dd, parameters):
    """"""
    p = parameters
    V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
    Vtb = V[2,2]
    Vts = V[2,1]
    e = sqrt(4 * pi * parameters['alpha_e'])
    gs = sqrt(4 * pi * parameters['alpha_s'])
    mb = parameters['m_b']
    dic = {"b->s::c7": 16 * pi**2 / mb / e * C["F7gamma" + dd],
            "b->s::c7'": 16 * pi**2 / mb / e * C["F7pgamma" + dd],
            "b->s::c8": 16 * pi**2 / mb / gs * C["F8g" + dd],
            "b->s::c8'": 16 * pi**2 / mb / gs * C["F8pg" + dd]
                }
    prefactor = sqrt(2)/p['GF']/Vtb/Vts.conj()/4
    return {k: prefactor * v for k,v in dic.items()}
236469
def Flavio_to_Fierz_nunu(C, ddll, parameters, norm_gf=True):
    """"""
    p = parameters
    V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
    if ddll[:2] == 'sb':
        xi = V[2, 2] * V[2, 1].conj()
    elif ddll[:2] == 'db':
        xi = V[2, 2] * V[2, 0].conj()
    elif ddll[:2] == 'ds':
        xi = V[2, 1] * V[2, 0].conj()
    else:
        raise ValueError("Unexpected flavours: {}".format(ddll[:2]))
    q1, q2 = ddll[:2]
    l1 = ddll[4:ddll.find('n')]
    l2 = ddll[ddll.find('_', 5) + 1:]
    ind = q1 + q2 + l1 + l2
    # flavio has indices within currents inverted
    indnu = q2 + q1 + 'nu' + l2 + 'nu' + l1
    e = sqrt(4* pi * parameters['alpha_e'])
    dic = {
        'F' + ind + 'nu': C["CL_" + indnu] / ((8 * pi**2) / e**2),
        'F' + ind + 'nup': C["CR_" + indnu] / ((8 * pi**2) / e**2),
    }
    if norm_gf:
        prefactor = sqrt(2)/p['GF']/xi/4
    else:
        prefactor = 1 / xi
    return {k: v / prefactor for k, v in dic.items()}
920716
def _Fierz_to_JMS_III_IV_V(Fqqqq, qqqq):
    """"""
    F = Fqqqq.copy()
    #case dduu
    classIII = ['sbuc', 'sbcu', 'dbuc', 'dbcu', 'dsuc', 'dscu']
    classVdduu = ['sbuu' , 'dbuu', 'dsuu', 'sbcc' , 'dbcc', 'dscc']
    if qqqq in classIII + classVdduu:
        f1 = str(dflav[qqqq[0]] + 1)
        f2 = str(dflav[qqqq[1]] + 1)
        f3 = str(uflav[qqqq[2]] + 1)
        f4 = str(uflav[qqqq[3]] + 1)
        d = {'V1udLL_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1'] + F['F' + qqqq + '2'] / Nc,
            'V8udLL_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '2'],
            'V1duLR_' + f1 + f2 + f3 + f4: F['F' + qqqq + '3'] + F['F' + qqqq + '4'] / Nc,
            'V8duLR_' + f1 + f2 + f3 + f4: 2 * F['F' + qqqq + '4'],
            'S1udRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '5'] + F['F' + qqqq + '6'] / Nc - 4 * F['F' + qqqq + '9'] - (4 * F['F' + qqqq + '10']) / Nc,
            'S8udRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '6'] - 8 * F['F' + qqqq + '10'],
            'S1udduRR_' + f3 + f2 + f1 + f4: -((8 * F['F' + qqqq + '9']) / Nc) - 8 * F['F' + qqqq + '10'],
            'V8udduLR_' + f4 + f1 + f2 + f3: -F['F' + qqqq + '7'].conjugate(),
            'V1udduLR_' + f4 + f1 + f2 + f3: -(F['F' + qqqq + '7'].conjugate() / (2 * Nc)) - F['F' + qqqq + '8'].conjugate() / 2,
            'S8udduRR_' + f3 + f2 + f1 + f4: -16 * F['F' + qqqq + '9'],
            'V1udRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1p'] + F['F' + qqqq + '2p'] / Nc,
            'V8udRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '2p'],
            'V1udLR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '3p'] + F['F' + qqqq + '4p'] / Nc,
            'V8udLR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '4p'],
            'S1udRR_' + f4 + f3 + f2 + f1: F['F' + qqqq + '5p'].conjugate() + F['F' + qqqq + '6p'].conjugate() / Nc - 4 * F['F' + qqqq + '9p'].conjugate() - (4 * F['F' + qqqq + '10p'].conjugate()) /  Nc,
            'S8udRR_' + f4 + f3 + f2 + f1: 2 * F['F' + qqqq + '6p'].conjugate() - 8 * F['F' + qqqq + '10p'].conjugate(),
            'S1udduRR_' + f4 + f1 + f2 + f3: -((8 * F['F' + qqqq + '9p'].conjugate()) / Nc) - 8 * F['F' + qqqq + '10p'].conjugate(),
            'V8udduLR_' + f3 + f2 + f1 + f4: -F['F' + qqqq + '7p'],
            'V1udduLR_' + f3 + f2 + f1 + f4: -(F['F' + qqqq + '7p'] / (2 * Nc)) - F['F' + qqqq + '8p'] / 2,
            'S8udduRR_' + f4 + f1 + f2 + f3: -16 * F['F' + qqqq + '9p'].conjugate(),
            }
        return symmetrize_JMS_dict(d)
    #case uudd
    classVuudd = ['ucdd', 'ucss','ucbb']
    if qqqq in classVuudd:
        f3 = str(uflav[qqqq[0]] + 1)
        f4 = str(uflav[qqqq[1]] + 1)
        f1 = str(dflav[qqqq[2]] + 1)
        f2 = str(dflav[qqqq[3]] + 1)
        d = {'V1udLL_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1'] + F['F' + qqqq + '2'] / Nc,
            'V8udLL_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '2'],
            'V1duLR_' + f1 + f2 + f3 + f4: F['F' + qqqq + '3p'] + F['F' + qqqq + '4p'] / Nc,
            'V8duLR_' + f1 + f2 + f3 + f4: 2 * F['F' + qqqq + '4p'],
            'S1udRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '5'] + F['F' + qqqq + '6'] / Nc - 4 * F['F' + qqqq + '9'] - (4 * F['F' + qqqq + '10']) / Nc,
            'S8udRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '6'] - 8 * F['F' + qqqq + '10'],
            'S1udduRR_' + f3 + f2 + f1 + f4: -((8 * F['F' + qqqq + '9']) / Nc) - 8 * F['F' + qqqq + '10'],
            'V8udduLR_' + f4 + f1 + f2 + f3: -F['F' + qqqq + '7p'].conjugate(),
            'V1udduLR_' + f4 + f1 + f2 + f3: -(F['F' + qqqq + '7p'].conjugate() / (2 * Nc)) - F['F' + qqqq + '8p'].conjugate() / 2,
            'S8udduRR_' + f3 + f2 + f1 + f4: -16 * F['F' + qqqq + '9'],
            'V1udRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1p'] + F['F' + qqqq + '2p'] / Nc,
            'V8udRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '2p'],
            'V1udLR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '3'] + F['F' + qqqq + '4'] / Nc,
            'V8udLR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '4'],
            'S1udRR_' + f4 + f3 + f2 + f1: F['F' + qqqq + '5p'].conjugate() + F['F' + qqqq + '6p'].conjugate() / Nc - 4 * F['F' + qqqq + '9p'].conjugate() - (4 * F['F' + qqqq + '10p'].conjugate()) /  Nc,
            'S8udRR_' + f4 + f3 + f2 + f1: 2 * F['F' + qqqq + '6p'].conjugate() - 8 * F['F' + qqqq + '10p'].conjugate(),
            'S1udduRR_' + f4 + f1 + f2 + f3: -((8 * F['F' + qqqq + '9p'].conjugate()) / Nc) - 8 * F['F' + qqqq + '10p'].conjugate(),
            'V8udduLR_' + f3 + f2 + f1 + f4: -F['F' + qqqq + '7'],
            'V1udduLR_' + f3 + f2 + f1 + f4: -(F['F' + qqqq + '7'] / (2 * Nc)) - F['F' + qqqq + '8'] / 2,
            'S8udduRR_' + f4 + f1 + f2 + f3: -16 * F['F' + qqqq + '9p'].conjugate(),
            }
        return symmetrize_JMS_dict(d)
    #case dddd
    classIV = ['sbsd', 'dbds', 'bsbd']
    classVdddd = ['sbss', 'dbdd', 'dsdd', 'sbbb', 'dbbb', 'dsss']
    classVddddind = ['sbdd', 'dsbb', 'dbss']
    classVuuuu = ['ucuu', 'cucc', 'uccc', 'cuuu']
    if qqqq in classVdddd + classIV + classVuuuu:
        # if 2nd and 4th or 1st and 3rd fields are the same, Fierz can be used
        # to express the even coeffs in terms of the odd ones
        for key in F:
            # to make sure we're not screwing things up, check that none
            # of the even WCs is actually present
            assert int(key[5:].replace('p', '')) % 2 == 1, "Unexpected key in Fierz basis: " + key
        for p in ['', 'p']:
            if qqqq in ['sbbb', 'dbbb', 'dsss', 'uccc']:
                F['F' + qqqq + '2' + p] = F['F' + qqqq + '1' + p]
                F['F' + qqqq + '4' + p] = -1 / 2 * F['F' + qqqq + '7' + p]
                F['F' + qqqq + '6' + p] = -1 / 2 * F['F' + qqqq + '5' + p] - 6 * F['F' + qqqq + '9' + p]
                F['F' + qqqq + '8' + p] = -2 * F['F' + qqqq + '3' + p]
                F['F' + qqqq + '10' + p] = -1 / 8 * F['F' + qqqq + '5' + p] + 1 / 2 * F['F' + qqqq + '9' + p]
            elif qqqq in ['sbss', 'dbdd', 'dsdd', 'sbsd', 'dbds', 'bsbd', 'ucuu']:
                notp = 'p' if p == '' else ''
                F['F' + qqqq + '2' + p] = F['F' + qqqq + '1' + p]
                F['F' + qqqq + '4' + p] = -1 / 2 * F['F' + qqqq + '7' + notp]
                F['F' + qqqq + '6' + notp] = -1 / 2 * F['F' + qqqq + '5' + notp] - 6 * F['F' + qqqq + '9' + notp]
                F['F' + qqqq + '8' + notp] = -2 * F['F' + qqqq + '3' + p]
                F['F' + qqqq + '10' + notp] = -1 / 8 * F['F' + qqqq + '5' + notp] + 1 / 2 * F['F' + qqqq + '9' + notp]
    if qqqq in classIV + classVdddd + classVddddind:
        f1 = str(dflav[qqqq[0]] + 1)
        f2 = str(dflav[qqqq[1]] + 1)
        f3 = str(dflav[qqqq[2]] + 1)
        f4 = str(dflav[qqqq[3]] + 1)
        d = {
        'VddLL_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1'],
        'VddLL_' + f1 + f4 + f3 + f2: F['F' + qqqq + '2'],
        'V1ddLR_' + f1 + f2 + f3 + f4: F['F' + qqqq + '3'] + F['F' + qqqq + '4'] / Nc,
        'V8ddLR_' + f1 + f2 + f3 + f4: 2 * F['F' + qqqq + '4'],
        'S1ddRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '5'] + F['F' + qqqq + '6'] / Nc - 4 * F['F' + qqqq + '9'] - (4 * F['F' + qqqq + '10']) / Nc,
        'S8ddRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '6'] - 8 * F['F' + qqqq + '10'],
        'V8ddLR_' + f1 + f4 + f3 + f2: -F['F' + qqqq + '7'],
        'V1ddLR_' + f1 + f4 + f3 + f2: -(F['F' + qqqq + '7'] / (2 * Nc)) - F['F' + qqqq + '8'] /  2,
        'S1ddRR_' + f1 + f4 + f3 + f2: -((8 * F['F' + qqqq + '9']) / Nc) - 8 * F['F' + qqqq + '10'],
        'S8ddRR_' + f3 + f2 + f1 + f4: -16 * F['F' + qqqq + '9'],
        'VddRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1p'],
        'VddRR_' + f1 + f4 + f3 + f2: F['F' + qqqq + '2p'],
        'V1ddLR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '3p'] + F['F' + qqqq + '4p'] / Nc,
        'V8ddLR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '4p'],
        'S1ddRR_' + f4 + f3 + f2 + f1: F['F' + qqqq + '5p'].conjugate() + F['F' + qqqq + '6p'].conjugate() / Nc - 4 * F['F' + qqqq + '9p'].conjugate() - (4 * F['F' + qqqq + '10p'].conjugate()) /  Nc,
        'S8ddRR_' + f4 + f3 + f2 + f1: 2 * F['F' + qqqq + '6p'].conjugate() - 8 * F['F' + qqqq + '10p'].conjugate(),
        'V8ddLR_' + f3 + f2 + f1 + f4: -F['F' + qqqq + '7p'],
        'V1ddLR_' + f3 + f2 + f1 + f4: -(F['F' + qqqq + '7p'] / (2 * Nc)) - F['F' + qqqq + '8p'] / 2,
        'S1ddRR_' + f4 + f1 + f2 + f3: -((8 * F['F' + qqqq + '9p'].conjugate()) / Nc) - 8 * F['F' + qqqq + '10p'].conjugate(),
        'S8ddRR_' + f4 + f1 + f2 + f3: -16 * F['F' + qqqq + '9p'].conjugate(),
        }
        return symmetrize_JMS_dict(d)
    #case uuuu
    if qqqq in classVuuuu:
        f1 = str(uflav[qqqq[0]] + 1)
        f2 = str(uflav[qqqq[1]] + 1)
        f3 = str(uflav[qqqq[2]] + 1)
        f4 = str(uflav[qqqq[3]] + 1)
        d = {
        'VuuLL_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1'],
        'VuuLL_' + f1 + f4 + f3 + f2: F['F' + qqqq + '2'],
        'V1uuLR_' + f1 + f2 + f3 + f4: F['F' + qqqq + '3'] + F['F' + qqqq + '4'] / Nc,
        'V8uuLR_' + f1 + f2 + f3 + f4: 2 * F['F' + qqqq + '4'],
        'S1uuRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '5'] + F['F' + qqqq + '6'] / Nc - 4 * F['F' + qqqq + '9'] - (4 * F['F' + qqqq + '10']) / Nc,
        'S8uuRR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '6'] - 8 * F['F' + qqqq + '10'],
        'V8uuLR_' + f1 + f4 + f3 + f2: -F['F' + qqqq + '7'],
        'V1uuLR_' + f1 + f4 + f3 + f2: -(F['F' + qqqq + '7'] / (2 * Nc)) - F['F' + qqqq + '8'] /  2,
        'S1uuRR_' + f1 + f4 + f3 + f2: -((8 * F['F' + qqqq + '9']) / Nc) - 8 * F['F' + qqqq + '10'],
        'S8uuRR_' + f3 + f2 + f1 + f4: -16 * F['F' + qqqq + '9'],
        'VuuRR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '1p'],
        'VuuRR_' + f1 + f4 + f3 + f2: F['F' + qqqq + '2p'],
        'V1uuLR_' + f3 + f4 + f1 + f2: F['F' + qqqq + '3p'] + F['F' + qqqq + '4p'] / Nc,
        'V8uuLR_' + f3 + f4 + f1 + f2: 2 * F['F' + qqqq + '4p'],
        'S1uuRR_' + f4 + f3 + f2 + f1: F['F' + qqqq + '5p'].conjugate() + F['F' + qqqq + '6p'].conjugate() / Nc - 4 * F['F' + qqqq + '9p'].conjugate() - (4 * F['F' + qqqq + '10p'].conjugate()) /  Nc,
        'S8uuRR_' + f4 + f3 + f2 + f1: 2 * F['F' + qqqq + '6p'].conjugate() - 8 * F['F' + qqqq + '10p'].conjugate(),
        'V8uuLR_' + f3 + f2 + f1 + f4: -F['F' + qqqq + '7p'],
        'V1uuLR_' + f3 + f2 + f1 + f4: -(F['F' + qqqq + '7p'] / (2 * Nc)) - F['F' + qqqq + '8p'] / 2,
        'S1uuRR_' + f4 + f1 + f2 + f3: -((8 * F['F' + qqqq + '9p'].conjugate()) / Nc) - 8 * F['F' + qqqq + '10p'].conjugate(),
        'S8uuRR_' + f4 + f1 + f2 + f3: -16 * F['F' + qqqq + '9p'].conjugate()
        }
        return symmetrize_JMS_dict(d)
    raise ValueError("Case not implemented: {}".format(qqqq))
92524
def read_g94(basis_lines, fname):
    '''
    '''

    skipchars = '!'
    basis_lines = [l for l in basis_lines if l and not l[0] in skipchars]

    bs_data = create_skel('component')

    i = 0
    if basis_lines[0] == '****':
        i += 1  # skip initial ****

    while i < len(basis_lines):
        line = basis_lines[i]
        elementsym = line.split()[0]

        # Some gaussian files have a dash before the element
        if elementsym[0] == '-':
            elementsym = elementsym[1:]

        element_Z = lut.element_Z_from_sym(elementsym)
        element_Z = str(element_Z)

        if not element_Z in bs_data['elements']:
            bs_data['elements'][element_Z] = {}

        element_data = bs_data['elements'][element_Z]

        i += 1

        # Try to guess if this is an ecp
        # Electron basis almost always end in 1.0 (scale factor)
        # ECP lines would end in an integer, so isdecimal() = true for ecp
        if basis_lines[i].split()[-1].isdecimal():
            if not 'ecp_potentials' in element_data:
                element_data['ecp_potentials'] = []

            lsplt = basis_lines[i].split()
            maxam = int(lsplt[1])
            n_elec = int(lsplt[2])
            element_data['ecp_electrons'] = n_elec

            # Highest AM first, then the rest in order
            am_list = list(range(maxam + 1))
            am_list.insert(0, am_list.pop())

            i += 1

            for j in range(maxam + 1):
                i += 1  # Skip the 'title' block - unused according to gaussian docs
                n_entries = int(basis_lines[i])
                i += 1  # Skip title block

                shell_am = am_list[j]
                ecp_shell = {'angular_momentum': [shell_am], 'ecp_type': 'scalar_ecp'}
                rexponents = []
                gexponents = []
                coefficients = []

                for k in range(n_entries):
                    lsplt = basis_lines[i].split()
                    rexponents.append(int(lsplt[0]))
                    gexponents.append(lsplt[1])
                    coefficients.append(lsplt[2:])
                    i += 1

                ecp_shell['r_exponents'] = rexponents
                ecp_shell['gaussian_exponents'] = gexponents

                # We need to transpose the coefficient matrix
                # (we store a matrix with primitives being the column index and
                # general contraction being the row index)
                ecp_shell['coefficients'] = list(map(list, zip(*coefficients)))

                element_data['ecp_potentials'].append(ecp_shell)
        else:
            if not 'electron_shells' in element_data:
                element_data['electron_shells'] = []

            while basis_lines[i] != '****':
                lsplt = basis_lines[i].split()
                shell_am = lut.amchar_to_int(lsplt[0])
                nprim = int(lsplt[1])

                if max(shell_am) <= 1:
                    func_type = 'gto'
                else:
                    func_type = 'gto_spherical'

                shell = {
                    'function_type': func_type,
                    'region': '',
                    'angular_momentum': shell_am
                }

                exponents = []
                coefficients = []

                i += 1
                for j in range(nprim):
                    line = basis_lines[i].replace('D', 'E')
                    line = line.replace('d', 'E')
                    lsplt = line.split()
                    exponents.append(lsplt[0])
                    coefficients.append(lsplt[1:])
                    i += 1

                shell['exponents'] = exponents

                # We need to transpose the coefficient matrix
                # (we store a matrix with primitives being the column index and
                # general contraction being the row index)
                shell['coefficients'] = list(map(list, zip(*coefficients)))

                element_data['electron_shells'].append(shell)

            i += 1

    return bs_data
536581
def flavor_rotation(C_in, Uq, Uu, Ud, Ul, Ue, sm_parameters=True):
    """"""
    C = {}
    if sm_parameters:
        # shift of theta terms, see 0907.4763
        C['Thetas'] = C_in['Thetas'] - 2*argdet(Uq) + argdet(Uu) + argdet(Ud)
        C['Theta'] = C_in['Theta'] - 3*argdet(Uq) - argdet(Ul)
        C['Thetap'] = ( C_in['Thetap'] - (1/6)*argdet(Uq) + (4/3)*argdet(Uu) + (1/3)*argdet(Ud)
                                       - (1/2)*argdet(Ul) + argdet(Ue))
        # nothing to do for scalar SM parameters
        for k in ['g', 'gp', 'gs', 'Lambda', 'm2']:
            C[k] = C_in[k]
        C['Ge'] = Ul.conj().T @ C_in['Ge'] @ Ue
        C['Gu'] = Uq.conj().T @ C_in['Gu'] @ Uu
        C['Gd'] = Uq.conj().T @ C_in['Gd'] @ Ud
    # nothing to do for purely bosonic operators
    for k in WC_keys_0f:
        C[k] = C_in[k]
    # see 1704.03888 table 4 (but staying SU(2) invariant here)
    # LR
    for k in ['ephi', 'eW', 'eB']:
        C[k] = Ul.conj().T @ C_in[k] @ Ue
    for k in ['uphi', 'uW', 'uB', 'uG']:
        C[k] = Uq.conj().T @ C_in[k] @ Uu
    for k in ['dphi', 'dW', 'dB', 'dG']:
        C[k] = Uq.conj().T @ C_in[k] @ Ud
    # LL
    for k in ['phil1', 'phil3']:
        C[k] = Ul.conj().T @ C_in[k] @ Ul
    for k in ['phiq1', 'phiq3']:
        C[k] = Uq.conj().T @ C_in[k] @ Uq
    C['llphiphi'] = Ul.T @ C_in['llphiphi'] @ Ul
    # RR
    C['phie'] = Ue.conj().T @ C_in['phie'] @ Ue
    C['phiu'] = Uu.conj().T @ C_in['phiu'] @ Uu
    C['phid'] = Ud.conj().T @ C_in['phid'] @ Ud
    C['phiud'] = Uu.conj().T @ C_in['phiud'] @ Ud
    # 4-fermion
    C['ll'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ul, Ul.conj(), Ul.conj(), C_in['ll'])
    C['ee'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Ue, Ue.conj(), Ue.conj(), C_in['ee'])
    C['le'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ue, Ul.conj(), Ue.conj(), C_in['le'])
    C['qq1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uq, Uq.conj(), Uq.conj(), C_in['qq1'])
    C['qq3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uq, Uq.conj(), Uq.conj(), C_in['qq3'])
    C['dd'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ud, Ud, Ud.conj(), Ud.conj(), C_in['dd'])
    C['uu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Uu, Uu.conj(), Uu.conj(), C_in['uu'])
    C['ud1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uu.conj(), Ud.conj(), C_in['ud1'])
    C['ud8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uu.conj(), Ud.conj(), C_in['ud8'])
    C['qu1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uu, Uq.conj(), Uu.conj(), C_in['qu1'])
    C['qu8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uu, Uq.conj(), Uu.conj(), C_in['qu8'])
    C['qd1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ud, Uq.conj(), Ud.conj(), C_in['qd1'])
    C['qd8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ud, Uq.conj(), Ud.conj(), C_in['qd8'])
    C['quqd1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uq.conj(), Uq.conj(), C_in['quqd1'])
    C['quqd8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uq.conj(), Uq.conj(), C_in['quqd8'])
    C['lq1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uq, Ul.conj(), Uq.conj(), C_in['lq1'])
    C['lq3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uq, Ul.conj(), Uq.conj(), C_in['lq3'])
    C['ld'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ud, Ul.conj(), Ud.conj(), C_in['ld'])
    C['lu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uu, Ul.conj(), Uu.conj(), C_in['lu'])
    C['qe'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ue, Uq.conj(), Ue.conj(), C_in['qe'])
    C['ed'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Ud, Ue.conj(), Ud.conj(), C_in['ed'])
    C['eu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ue.conj(), Uu.conj(), C_in['eu'])
    C['ledq'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uq, Ul.conj(), Ud.conj(), C_in['ledq'])
    C['lequ1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ul.conj(), Uq.conj(), C_in['lequ1'])
    C['lequ3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ul.conj(), Uq.conj(), C_in['lequ3'])
    C['duql'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ul, Ud, Uq, C_in['duql'])
    C['qque'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ue, Uq, Uu, C_in['qque'])
    C['qqql'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ul, Uq, Uq, C_in['qqql'])
    C['duue'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ue, Ud, Uu, C_in['duue'])
    return C