Generates the generic spin operator in z basis for a system of
        N=particles and for the selected spin index name. where index=0..N-1
        The gauge term sets the behavoir for a system away from half-filling
1145
def spin_gen(particles, index, gauge=1):
    """"""
    mat = np.zeros((2**particles, 2**particles))

    flipper = 2**index
    for i in range(2**particles):
        ispin = btest(i, index)
        if ispin == 1:
            mat[i ^ flipper, i] = 1
        else:
            mat[i ^ flipper, i] = gauge
    return mat
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
1004152
def add_spin_by_element(self, spins):
        
        for site in self.sites:
            new_sp = {}
            for sp, occu in site.species.items():
                sym = sp.symbol
                oxi_state = getattr(sp, "oxi_state", None)
                new_sp[Specie(sym, oxidation_state=oxi_state,
                              properties={'spin': spins.get(str(sp), spins.get(sym, None))})] = occu
            site.species = new_sp
331279
def primary_xi(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
    
    spinx = primary_spin(mass1, mass2, spin1x, spin2x)
    spiny = primary_spin(mass1, mass2, spin1y, spin2y)
    return chi_perp_from_spinx_spiny(spinx, spiny)
739802
def _generate_a_indptr(num_states, s_indices, out):
    
    idx = 0
    out[0] = 0
    for s in range(num_states-1):
        while(s_indices[idx] == s):
            idx += 1
        out[s+1] = idx
    out[num_states] = len(s_indices)
554853
public float[] t4(float[] z, int k, int M) {
    float[] result = new float[M];
    float[] w = new float[z.length];

    for (int i = 0; i < z.length; i++) {
      w[i] = (float) 2.0 * (i + 1);
    }

    for (int i = 1; i <= M - 1; i++) {
      int head = (i - 1) * k / (M - 1) + 1;
      int tail = i * k / (M - 1);
      float[] subZ = subVector(z, head - 1, tail - 1);
      float[] subW = subVector(w, head - 1, tail - 1);

      result[i - 1] = (new Transformations()).rSum(subZ, subW);
    }

    int head = k + 1 - 1;
    int tail = z.length - 1;
    float[] subZ = subVector(z, head, tail);
    float[] subW = subVector(w, head, tail);
    result[M - 1] = (new Transformations()).rSum(subZ, subW);

    return result;
  }
821664
def yield_pair_energies(self, index1, index2):
        """"""
        d_1 = 1/self.distances[index1, index2]
        if self.charges is not None:
            c1 = self.charges[index1]
            c2 = self.charges[index2]
            yield c1*c2*d_1, 1
        if self.dipoles is not None:
            d_3 = d_1**3
            d_5 = d_1**5
            delta = self.deltas[index1, index2]
            p1 = self.dipoles[index1]
            p2 = self.dipoles[index2]
            yield d_3*np.dot(p1, p2), 1
            yield -3*d_5, np.dot(p1, delta)*np.dot(delta, p2)
            if self.charges is not None:
                yield c1*d_3, np.dot(p2, delta)
                yield c2*d_3, np.dot(p1, -delta)
1041330
def spinflipandhop(slaves):
    """"""

    Sdw = [csr_matrix(spin_gen(slaves, i, 0)) for i in range(slaves)]
    Sup = [mat.T for mat in Sdw]

    sfh = np.zeros_like(Sup[0])
    orbitals = slaves//2
    for n in range(orbitals):
        for m in range(n+1, orbitals):
            sfh += Sup[2*n  ] * Sdw[2*n + 1] * Sup[2*m + 1] * Sdw[2*m  ]
            sfh += Sup[2*n+1] * Sdw[2*n    ] * Sup[2*m    ] * Sdw[2*m+1]

            sfh += Sup[2*n] * Sup[2*n + 1] * Sdw[2*m] * Sdw[2*m+1]
            sfh += Sup[2*m] * Sup[2*m + 1] * Sdw[2*n] * Sdw[2*n+1]

    return sfh
413121
def operator(self):
        
        # is projection the zero operator?
        if self.V.shape[1] == 0:
            N = self.V.shape[0]
            return ZeroLinearOperator((N, N))
        return self._get_operator(self.apply, self.apply_adj)
973040
def get_bond_symmetry(site_symmetry,
                      lattice,
                      positions,
                      atom_center,
                      atom_disp,
                      symprec=1e-5):
    
    bond_sym = []
    pos = positions
    for rot in site_symmetry:
        rot_pos = (np.dot(pos[atom_disp] - pos[atom_center], rot.T) +
                   pos[atom_center])
        diff = pos[atom_disp] - rot_pos
        diff -= np.rint(diff)
        dist = np.linalg.norm(np.dot(lattice, diff))
        if dist < symprec:
            bond_sym.append(rot)

    return np.array(bond_sym)
1052669
def as_spinor_array(a):
    
    a = np.atleast_1d(a)
    assert a.dtype == np.dtype(np.quaternion)
    # I'm not sure why it has to be so complicated, but all of these steps
    # appear to be necessary in this case.
    return a.view(np.float).reshape(a.shape + (4,))[..., [0, 3, 2, 1]].ravel().view(np.complex).reshape(a.shape + (2,))
353855
def from_Z(z: int):
        
        for sym, data in _pt_data.items():
            if data["Atomic no"] == z:
                return Element(sym)
        raise ValueError("No element with this atomic number %s" % z)
858124
def Jzjmcoeff(ls, m, shift) -> sympy.Expr:
    
    '''
    assert isinstance(ls, SpinSpace)
    n = ls.dimension
    s = sympify(n - 1) / 2
    assert n == int(2 * s + 1)
    if isinstance(m, str):
        return ls.basis.index(m) - s
    elif isinstance(m, int):
        if shift:
            assert 0 <= m < n
            return m - s
    else:
        return sympify(m)
734044
def pauli_product(*elements: Pauli) -> Pauli:
    """"""
    result_terms = []

    for terms in product(*elements):
        coeff = reduce(mul, [term[1] for term in terms])
        ops = (term[0] for term in terms)
        out = []
        key = itemgetter(0)
        for qubit, qops in groupby(heapq.merge(*ops, key=key), key=key):
            res = next(qops)[1]  # Operator: X Y Z
            for op in qops:
                pair = res + op[1]
                res, rescoeff = PAULI_PROD[pair]
                coeff *= rescoeff
            if res != 'I':
                out.append((qubit, res))

        p = Pauli(((tuple(out), coeff),))
        result_terms.append(p)

    return pauli_sum(*result_terms)
1047762
def spin_sz(self):
        """"""
        return conversions.secondary_spin(self.mass1, self.mass2, self.spin1z,
                                        self.spin2z)
645024
def follow_cf(save, Uspan, target_cf, nup, n_tot=5.0, slsp=None):
    

    if slsp == None:
        slsp = Spinon(slaves=6, orbitals=3, avg_particles=n_tot,
                       hopping=[0.5]*6, populations = np.asarray([n_tot]*6)/6)

    zet, lam, mu, mean_f = [], [], [], []

    for co in Uspan:
        print('U=', co, 'del=', target_cf)
        res=root(targetpop, nup[-1],(co,target_cf,slsp, n_tot))
        print(res.x)
        if res.x>nup[-1]: break

        nup.append(res.x)
        slsp.param['populations']=population_distri(nup[-1])
        mean_f.append(slsp.mean_field())
        zet.append(slsp.quasiparticle_weight())
        lam.append(slsp.param['lambda'])
        mu.append(orbital_energies(slsp.param, zet[-1]))

#    plt.plot(np.asarray(zet)[:,0], label='d={}, zl'.format(str(target_cf)))
#    plt.plot(np.asarray(zet)[:,5], label='d={}, zh'.format(str(target_cf)))

    case = save.createGroup('cf={}'.format(target_cf))
    varis = st.setgroup(case)
    st.storegroup(varis, Uspan[:len(zet)], zet, lam, mu, nup[1:],target_cf,mean_f)
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)