Generates the plot on the convergenge of the mean field in single
       site spin hamiltonian under with N degenerate half-filled orbitals
496203
def plot_mean_field_conv(N=1, n=0.5, Uspan=np.arange(0, 3.6, 0.5)):
    

    sl = Spinon(slaves=2*N, orbitals=N, avg_particles=2*n,
                       hopping=[0.5]*2*N, orbital_e=[0]*2*N)
    hlog = solve_loop(sl, Uspan, [0.])[1]

    f, (ax1, ax2) = plt.subplots(2, sharex=True)
    for field in hlog:
        field = np.asarray(field)
        ax1.semilogy(abs(field[1:]-field[:-1]))
        ax2.plot(field)#, label = 'h, U = {}'.format(Uint))

    plt.title('Convergence of selfconsintent mean field')
    ax1.set_ylabel('$\\Delta h$')
    ax2.set_ylabel('mean field $h$')
    plt.xlabel('iterations')

    return hlog
149239
def spin_hamiltonian(self, h, l):
        """"""
        h_spin  = np.einsum('i,ijk', h[1], self.oper['O'])
        h_spin += np.einsum('i,ijk', h[0], self.oper['O_d'])
        h_spin += np.einsum('i,ijk', l, self.oper['Sz+1/2'])

        h_spin += self.oper['Hint']
        return h_spin
743362
def gaussian2d(xy, sx, sy, mx=0, my=0, rho=0, amp=1, offs=0):
    '''
    
    '''
    x,y = xy
    return offs+amp*( 
        1/(2*np.pi*sx*sy*(1-(rho**2))**0.5) *
         np.exp( (-1/(2*(1-rho**2))) *
                 (
                    ( (x-mx)**2/sx**2 )
                  + ( (y-my)**2/sy**2 )
                  - ( ( 2*rho*(x-mx)*(y-my)) / (sx*sy) )
                  )
                )
         )
424402
def term_symbols(self):
        
        L_symbols = 'SPDFGHIKLMNOQRTUVWXYZ'

        L, v_e = self.valence

        # for one electron in subshell L
        ml = list(range(-L, L + 1))
        ms = [1 / 2, -1 / 2]
        # all possible configurations of ml,ms for one e in subshell L
        ml_ms = list(product(ml, ms))

        # Number of possible configurations for r electrons in subshell L.
        n = (2 * L + 1) * 2
        # the combination of n_e electrons configurations
        # C^{n}_{n_e}
        e_config_combs = list(combinations(range(n), v_e))

        # Total ML = sum(ml1, ml2), Total MS = sum(ms1, ms2)
        TL = [sum([ml_ms[comb[e]][0] for e in range(v_e)])
              for comb in e_config_combs]
        TS = [sum([ml_ms[comb[e]][1] for e in range(v_e)])
              for comb in e_config_combs]
        comb_counter = Counter([r for r in zip(TL, TS)])

        term_symbols = []
        while sum(comb_counter.values()) > 0:
            # Start from the lowest freq combination,
            # which corresponds to largest abs(L) and smallest abs(S)
            L, S = min(comb_counter)

            J = list(np.arange(abs(L - S), abs(L) + abs(S) + 1))
            term_symbols.append([str(int(2 * (abs(S)) + 1))
                                 + L_symbols[abs(L)]
                                 + str(j) for j in J])
            # Without J
            # term_symbols.append(str(int(2 * (abs(S)) + 1)) \
            #                     + L_symbols[abs(L)])

            # Delete all configurations included in this term
            for ML in range(-L, L - 1, -1):
                for MS in np.arange(S, -S + 1, 1):
                    if (ML, MS) in comb_counter:

                        comb_counter[(ML, MS)] -= 1
                        if comb_counter[(ML, MS)] == 0:
                            del comb_counter[(ML, MS)]
        return term_symbols
732797
def get_interpolated_value(self, energy):
        
        f = {}
        for spin in self.densities.keys():
            f[spin] = get_linear_interpolated_value(self.energies,
                                                    self.densities[spin],
                                                    energy)
        return f
581528
def multiplicity(self):
        
        if self._multiplicity is None:
            # generate multiplicity based on space group symmetry operations performed on defect coordinates
            try:
                d_structure = create_saturated_interstitial_structure(self)
            except ValueError:
                logger.debug('WARNING! Multiplicity was not able to be calculated adequately '
                             'for interstitials...setting this to 1 and skipping for now...')
                return 1

            sga = SpacegroupAnalyzer(d_structure)
            periodic_struc = sga.get_symmetrized_structure()
            poss_deflist = sorted(
                periodic_struc.get_sites_in_sphere(self.site.coords, 2, include_index=True),
                key=lambda x: x[1])
            defindex = poss_deflist[0][2]

            equivalent_sites = periodic_struc.find_equivalent_sites(periodic_struc[defindex])
            return len(equivalent_sites)

        else:
            return self._multiplicity
1134771
def smeftpar(scale, C, basis):
    """"""
    # start with a zero dict and update it with the input values
    MW = p['m_W']
    # MZ = p['m_Z']
    GF = p['GF']
    Mh = p['m_h']
    vb = sqrt(1 / sqrt(2) / GF)
    v = vb  # TODO
    _d = vMh2_to_m2Lambda(v=v, Mh2=Mh**2, C=C)
    m2 = _d['m2'].real
    Lambda = _d['Lambda'].real
    gsbar = sqrt(4 * pi * p['alpha_s'])
    gs = (1 - C['phiG'] * (v**2)) * gsbar
    gbar = 2 * MW / v
    g = gbar * (1 - C['phiW'] * (v**2))
    ebar = sqrt(4 * pi * p['alpha_e'])
    gp = get_gpbar(ebar, gbar, v, C)
    c = {}
    c['m2'] = m2
    c['Lambda'] = Lambda
    c['g'] = g
    c['gp'] = gp
    c['gs'] = gs
    K = ckmutil.ckm.ckm_tree(p['Vus'], p['Vub'], p['Vcb'], p['delta'])
    if basis == 'Warsaw':
        Mu = K.conj().T @ np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = np.diag([p['m_d'], p['m_s'], p['m_b']])
    elif basis == 'Warsaw up':
        Mu = np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = K @ np.diag([p['m_d'], p['m_s'], p['m_b']])
    else:
        raise ValueError("Basis '{}' not supported".format(basis))
    Me = np.diag([p['m_e'], p['m_mu'], p['m_tau']])
    c['Gd'] = Md / (v / sqrt(2)) + C['dphi'] * (v**2) / 2
    c['Gu'] = Mu / (v / sqrt(2)) + C['uphi'] * (v**2) / 2
    c['Ge'] = Me / (v / sqrt(2)) + C['ephi'] * (v**2) / 2
    return c
689442
def mass_2d(self, R, Rs, rho0, r_core):

        

        x = R * Rs ** -1
        p = Rs * r_core ** -1
        gx = self._G(x, p)

        m_2d = 2 * np.pi * rho0 * Rs ** 3 * gx

        return m_2d
1280216
def evolve(self, rho: Density) -> Density:
        """"""
        N = rho.qubit_nb
        qubits = rho.qubits

        indices = list([qubits.index(q) for q in self.qubits]) + \
            list([qubits.index(q) + N for q in self.qubits])

        tensor = bk.tensormul(self.tensor, rho.tensor, indices)
        return Density(tensor, qubits, rho.memory)
1089703
def pressision_try(orbitals, U, beta, step):
    """"""
    mu, lam = main(orbitals, U, beta, step)
    mu2, lam2 = linspace(0, U*orbitals, step), zeros(step)
    for i in range(99):
        lam2[i+1] = fsolve(restriction, lam2[i], (mu2[i+1], orbitals, U, beta))
    plot(mu2, 2*orbitals*fermi_dist(-(mu2+lam2), beta), label='Test guess')
    legend(loc=0)
19857
def set_filling(self, populations):
        """"""
        populations = np.asarray(populations)
#
#        self.param['orbital_e'] -= bethe_findfill_zeroT( \
#                                        self.param['avg_particles'],
#                                        self.param['orbital_e'],
#                                        self.param['hopping'])
        efermi = - bethe_find_crystalfield(
                            populations, self.param['hopping'])
        self.param['populations'] = populations
#        fermion_avg(efermi, self.param['hopping'], 'ocupation')
        self.param['ekin'] = fermion_avg(efermi, self.param['hopping'], 'ekin')

        spin_gen_op(self.oper, estimate_gauge(populations))
162158
def coordination_number_delta_E( self ):
        
        initial_site_neighbours = [ s for s in self.initial_site.p_neighbours if s.is_occupied ] # excludes final site, since this is always unoccupied
        final_site_neighbours = [ s for s in self.final_site.p_neighbours if s.is_occupied and s is not self.initial_site ] # excludes initial site
        initial_cn_occupation_energy = ( self.initial_site.cn_occupation_energy() + 
            sum( [ site.cn_occupation_energy() for site in initial_site_neighbours ] ) +
            sum( [ site.cn_occupation_energy() for site in final_site_neighbours ] ) )
        final_cn_occupation_energy = ( self.final_site.cn_occupation_energy( delta_occupation = { self.initial_site.label : -1 } ) +
            sum( [ site.cn_occupation_energy( delta_occupation = { self.initial_site.label : -1 } ) for site in initial_site_neighbours ] ) +
            sum( [ site.cn_occupation_energy( delta_occupation = { self.final_site.label : +1 } ) for site in final_site_neighbours ] ) )
        return ( final_cn_occupation_energy - initial_cn_occupation_energy )
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
1124921
def repeat(nx=1, ny=1, nz=1)
      
      raise "Not a periodic system." if self.lattice_vectors.nil?
      
      u = self.copy
      v1 = self.lattice_vectors[0]
      v2 = self.lattice_vectors[1]
      v3 = self.lattice_vectors[2]

      nx_sign = (0 < nx) ? 1 : -1
      ny_sign = (0 < ny) ? 1 : -1
      nz_sign = (0 < nz) ? 1 : -1
      
      new_atoms = []
      nx.to_i.abs.times do |i|        
        ny.to_i.abs.times do |j|
          nz.to_i.abs.times do |k|
            new_atoms << self.displace(nx_sign*i*v1[0] + ny_sign*j*v2[0] + nz_sign*k*v3[0], 
                                       nx_sign*i*v1[1] + ny_sign*j*v2[1] + nz_sign*k*v3[1], 
                                       nx_sign*i*v1[2] + ny_sign*j*v2[2] + nz_sign*k*v3[2]).atoms
          end
        end
      end

      u.atoms = new_atoms.flatten
      u.lattice_vectors = [Vector[nx.abs*v1[0], nx.abs*v1[1], nx.abs*v1[2]], 
                           Vector[ny.abs*v2[0], ny.abs*v2[1], ny.abs*v2[2]], 
                           Vector[nz.abs*v3[0], nz.abs*v3[1], nz.abs*v3[2]]]
      # u.make_bonds 
      return u     
    end
781387
def get_projection_on_elements(self, structure):
        
        dico = {}
        for spin in self.data.keys():
            dico[spin] = [[defaultdict(float)
                           for i in range(self.nkpoints)]
                          for j in range(self.nbands)]

        for iat in range(self.nions):
            name = structure.species[iat].symbol
            for spin, d in self.data.items():
                for k, b in itertools.product(range(self.nkpoints),
                                              range(self.nbands)):
                    dico[spin][b][k][name] = np.sum(d[k, b, iat, :])

        return dico
1140865
def _compute_mean_on_rock(self, C, mag, rrup, rvol, hypo_depth, CN, CR,
                              f4HW):
        

        # Define subduction flag (page 23)
        # SI=1 for subduction interface, 0 otherwise
        # DS=1 for subduction intraslab, 0 otherwise
        SI = 0
        DS = 1

        lnSA_AB = (
            # line 1 and 2 of equation 2
            C['c11'] + (C['c12y'] + (C['c15'] - C['c17']) * C['c19y']) *
            (mag - 6) +
            # line 3
            C['c13y'] * (10 - mag) ** 3 +
            # line 4
            C['c17'] * np.log(rrup + C['c18y'] * np.exp(C['c19y'] * mag)) +
            # line 5
            C['c20'] * hypo_depth + C['c24'] * SI +
            # line 6
            C['c46'] * rvol * (1 - DS)
        )

        return lnSA_AB
1003465
def infos_on_basis_set(self):
        
        o = []
        o.append("=========================================")
        o.append("Reading basis set:")
        o.append("")
        o.append(" Basis set for {} atom ".format(str(self.filename)))
        o.append(" Maximum angular momentum = {}".format(self.data['lmax']))
        o.append(" Number of atomics orbitals = {}".format(self.data['n_nlo']))
        o.append(" Number of nlm orbitals = {}".format(self.data['n_nlmo']))
        o.append("=========================================")

        return str(0)