Python numpy.integer() Examples

def _add_choices(self, nchoices):
        if isinstance(nchoices, int):
            self.nchoices = nchoices
            self.choice_names = None
        elif isinstance(nchoices, list) or nchoices.__class__.__name__ == "Series" or nchoices.__class__.__name__ == "DataFrame":
            self.choice_names = np.array(nchoices).reshape(-1)
            self.nchoices = self.choice_names.shape[0]
            if np.unique(self.choice_names).shape[0] != self.choice_names.shape[0]:
                raise ValueError("Arm/choice names contain duplicates.")
        elif isinstance(nchoices, np.ndarray):
            self.choice_names = nchoices.reshape(-1)
            self.nchoices = self.choice_names.shape[0]
            if np.unique(self.choice_names).shape[0] != self.choice_names.shape[0]:
                raise ValueError("Arm/choice names contain duplicates.")
        else:
            raise ValueError("'nchoices' must be an integer or list with named arms.") 

def make_rdm12(self, fcivec, norb, nelec, link_index=None, **kwargs):
        if isinstance(nelec, (int, numpy.integer)):
            nelectrons = nelec
        else:
            nelectrons = nelec[0]+nelec[1]

        nstates = len(self.state_weights)

        # If norm != 1 then the state weights will need normalising.
        norm = sum(self.state_weights)

        two_pdm = numpy.zeros( (norb, norb, norb, norb) )

        for irdm in range(nstates):
            if self.state_weights[irdm] != 0.0:
                dm_filename = 'spinfree_TwoRDM.' + str(irdm+1)
                temp_dm = read_neci_two_pdm(self, dm_filename, norb,
                                            self.scratchDirectory)
                two_pdm += (self.state_weights[irdm]/norm)*temp_dm

        one_pdm = one_from_two_pdm(two_pdm, nelectrons)

        return one_pdm, two_pdm 

def kernel(self, h1e, eri, norb, nelec, fci_restart=None, ecore=0, **kwargs):
        if fci_restart is None:
            fci_restart = self.restart
        if isinstance(nelec, (int, numpy.integer)):
            neleca = nelec//2 + nelec%2
            nelecb = nelec - neleca
        else:
            neleca, nelecb = nelec

        write_integrals_file(h1e, eri, norb, neleca, nelecb, self, ecore)
        if self.generate_neci_input:
            write_fciqmc_config_file(self, neleca, nelecb, fci_restart)
        if self.verbose >= logger.DEBUG1:
            in_file = self.configFile
            logger.debug1(self, 'FCIQMC Input file')
            logger.debug1(self, open(in_file, 'r').read())
        execute_fciqmc(self)
        if self.verbose >= logger.DEBUG1:
            out_file = self.outputFileCurrent
            with open(out_file) as f:
                logger.debug1(self, f.read())
        rdm_energy = read_energy(self)

        return rdm_energy, None 

def make_rdm12_forSQA(self, state, norb, nelec, link_index=None, **kwargs):
        nelectrons = 0
        if isinstance(nelec, (int, numpy.integer)):
            nelectrons = nelec
        else:
            nelectrons = nelec[0] + nelec[1]

        # The 2RDMs written by "SHCIrdm::saveRDM" in DICE
        # are written as E2[i1,j2,k1,l2]
        # and stored here as E2[i1,k1,j2,l2] (for PySCF purposes)
        # This is NOT done with SQA in mind.
        twopdm = numpy.zeros((norb, norb, norb, norb))
        file2pdm = "spatialRDM.%d.%d.txt" % (state, state)
        r2RDM(twopdm, norb, os.path.join(self.scratchDirectory, file2pdm).endcode())
        twopdm = twopdm.transpose(0, 2, 1, 3)

        # (This is coherent with previous statement about indexes)
        onepdm = numpy.einsum("ijkj->ki", twopdm)
        onepdm /= nelectrons - 1
        return onepdm, twopdm 

def make_rdm1(civec, norb, nelec, link_index=None):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    assert(neleca == nelecb)

    if link_index is None:
        link_index = cistring.gen_linkstr_index(range(norb), neleca)
    na = cistring.num_strings(norb, neleca)
    t1 = numpy.zeros((norb,na))
    #:for str0, tab in enumerate(link_index):
    #:    for a, i, str1, sign in tab:
    #:        if a == i:
    #:            t1[i,str1] += civec[str0]
    link1 = link_index[link_index[:,:,0] == link_index[:,:,1]].reshape(na,-1,4)
    t1[link1[:,:,1],link1[:,:,2]] = civec[:,None]

    dm1 = numpy.diag(numpy.einsum('ip,p->i', t1, civec)) * 2
    return dm1 

def kernel(self, h1e, eri, norb, nelec, ci0=None, ecore=0, **kwargs):
        if isinstance(nelec, (int, numpy.integer)):
            nelecb = nelec//2
            neleca = nelec - nelecb
        else:
            neleca, nelecb = nelec
        h2e = self.absorb_h1e(h1e, eri, norb, nelec, .5)
        h2e = ao2mo.restore(1, h2e, norb)

        hdiag = self.make_hdiag(h1e, eri, norb, nelec)
        nroots = 1
        if ci0 is None:
            ci0 = self.get_init_guess(norb, nelec, nroots, hdiag)

        def hop(c):
            return self.contract_2e(h2e, c, norb, nelec)
        precond = lambda x, e, *args: x/(hdiag-e+1e-4)
        e, c = lib.davidson(hop, ci0, precond, **kwargs)
        return e+ecore, c 

def build(self, mol=None):
        if mol is None: mol = self.mol
        if mol.symmetry:
            for irname in self.irrep_nelec:
                if irname not in self.mol.irrep_name:
                    logger.warn(self, 'No irrep %s', irname)

            fix_na, fix_nb = check_irrep_nelec(mol, self.irrep_nelec, self.nelec)[:2]
            alpha_open = beta_open = False
            for ne in self.irrep_nelec.values():
                if not isinstance(ne, (int, numpy.integer)):
                    alpha_open |= ne[0] > ne[1]
                    beta_open  |= ne[0] < ne[1]

            frozen_spin = fix_na - fix_nb
            if ((alpha_open and beta_open) or
                (0 < mol.spin < frozen_spin) or (frozen_spin < 0 < mol.spin) or
                (frozen_spin < mol.spin < 0) or (mol.spin < 0 < frozen_spin)):
                raise ValueError('Low-spin configuration was found in '
                                 'the irrep_nelec input. ROHF does not '
                                 'support low-spin configuration.')
        return hf.RHF.build(self, mol) 

def make_rdm12(self, fcivec, norb, nelec, link_index=None, **kwargs):
        if isinstance(nelec, (int, numpy.integer)):
            nelectrons = nelec
        else:
            nelectrons = nelec[0]+nelec[1]

        nstates = len(self.state_weights)

        # If norm != 1 then the state weights will need normalising.
        norm = sum(self.state_weights)

        two_pdm = numpy.zeros( (norb, norb, norb, norb) )

        for irdm in range(nstates):
            if self.state_weights[irdm] != 0.0:
                dm_filename = 'spinfree_TwoRDM.' + str(irdm+1)
                temp_dm = read_neci_two_pdm(self, dm_filename, norb,
                                            self.scratchDirectory)
                two_pdm += (self.state_weights[irdm]/norm)*temp_dm

        one_pdm = one_from_two_pdm(two_pdm, nelectrons)

        return one_pdm, two_pdm 

def kernel(self, h1e, eri, norb, nelec, fci_restart=None, ecore=0, **kwargs):
        if fci_restart is None:
            fci_restart = self.restart
        if isinstance(nelec, (int, numpy.integer)):
            neleca = nelec//2 + nelec%2
            nelecb = nelec - neleca
        else:
            neleca, nelecb = nelec

        write_integrals_file(h1e, eri, norb, neleca, nelecb, self, ecore)
        if self.generate_neci_input:
            write_fciqmc_config_file(self, neleca, nelecb, fci_restart)
        if self.verbose >= logger.DEBUG1:
            in_file = self.configFile
            logger.debug1(self, 'FCIQMC Input file')
            logger.debug1(self, open(in_file, 'r').read())
        execute_fciqmc(self)
        if self.verbose >= logger.DEBUG1:
            out_file = self.outputFileCurrent
            with open(out_file) as f:
                logger.debug1(self, f.read())
        rdm_energy = read_energy(self)

        return rdm_energy, None 

def trans_rdm12(self, statebra, stateket, norb, nelec, link_index=None, **kwargs):
        nelectrons = 0
        if isinstance(nelec, (int, numpy.integer)):
            nelectrons = nelec
        else:
            nelectrons = nelec[0] + nelec[1]

        writeSHCIConfFile(self, nelec, True)
        executeSHCI(self)

        # The 2RDMs written by "SHCIrdm::saveRDM" in DICE
        # are written as E2[i1,j2,k1,l2]
        # and stored here as E2[i1,k1,j2,l2] (for PySCF purposes)
        # This is NOT done with SQA in mind.
        twopdm = numpy.zeros((norb, norb, norb, norb))
        file2pdm = "spatialRDM.%d.%d.txt" % (statebra, stateket)
        r2RDM(twopdm, norb, os.path.join(self.scratchDirectory, file2pdm).endcode())

        # (This is coherent with previous statement about indexes)
        onepdm = numpy.einsum("ikjj->ki", twopdm)
        onepdm /= nelectrons - 1
        return onepdm, twopdm 

def make_rdm1s(self, state, norb, nelec, link_index=None, **kwargs):
        # Ref: IJQC, 109, 3552 Eq (3)
        if isinstance(nelec, (int, numpy.integer)):
            nelecb = (nelec-self.spin) // 2
            neleca = nelec - nelecb
        else :
            neleca, nelecb = nelec

        # DO NOT call self.make_rdm12. Calling DMRGCI.make_rdm12 instead of
        # self.make_rdm12 because self.make_rdm12 may be modified
        # by state-average mcscf solver (see function mcscf.addons.state_average).
        # When calling make_rdm1s from state-average FCI solver,
        # DMRGCI.make_rdm12 ensures that the basic make_rdm12 method is called.
        # (Issue https://github.com/pyscf/pyscf/issues/335)
        dm1, dm2 = DMRGCI.make_rdm12(self, state, norb, nelec, link_index, **kwargs)
        dm1n = (2-(neleca+nelecb)/2.) * dm1 - numpy.einsum('pkkq->pq', dm2)
        dm1n *= 1./(neleca-nelecb+1)
        dm1a, dm1b = (dm1+dm1n)*.5, (dm1-dm1n)*.5
        return dm1a, dm1b 

def make_rdm12(self, state, norb, nelec, link_index=None, **kwargs):
        nelectrons = 0
        if isinstance(nelec, (int, numpy.integer)):
          nelectrons = nelec
        else:
          nelectrons = nelec[0]+nelec[1]

        # The 2RDMs written by "save_spatial_twopdm_text" in BLOCK and STACKBLOCK
        # are written as E2[i1,j2,k2,l1]
        # and stored here as E2[i1,l1,j2,k2] (for PySCF purposes)
        # This is NOT done with SQA in mind.
        twopdm = numpy.zeros( (norb, norb, norb, norb) )
        file2pdm = "spatial_twopdm.%d.%d.txt" %(state, state)
        with open(os.path.join(self.scratchDirectory, "node0", file2pdm), "r") as f:
            norb_read = int(f.readline().split()[0])
            assert(norb_read == norb)
            for line in f:
                linesp = line.split()
                i, k, l, j = [int(x) for x in linesp[:4]]
                twopdm[i,j,k,l] = 2.0 * float(linesp[4])

        # (This is coherent with previous statement about indexes)
        onepdm = numpy.einsum('ikjj->ki', twopdm)
        onepdm /= (nelectrons-1)
        return onepdm, twopdm 

def get_nocc(mp):
    if mp._nocc is not None:
        return mp._nocc
    elif mp.frozen is None:
        nocc = numpy.count_nonzero(mp.mo_occ > 0)
        assert(nocc > 0)
        return nocc
    elif isinstance(mp.frozen, (int, numpy.integer)):
        nocc = numpy.count_nonzero(mp.mo_occ > 0) - mp.frozen
        assert(nocc > 0)
        return nocc
    elif isinstance(mp.frozen[0], (int, numpy.integer)):
        occ_idx = mp.mo_occ > 0
        occ_idx[list(mp.frozen)] = False
        nocc = numpy.count_nonzero(occ_idx)
        assert(nocc > 0)
        return nocc
    else:
        raise NotImplementedError 

def get_frozen_mask(mp):
    '''Get boolean mask for the restricted reference orbitals.

    In the returned boolean (mask) array of frozen orbital indices, the
    element is False if it corresonds to the frozen orbital.
    '''
    moidx = numpy.ones(mp.mo_occ.size, dtype=numpy.bool)
    if mp._nmo is not None:
        moidx[mp._nmo:] = False
    elif mp.frozen is None:
        pass
    elif isinstance(mp.frozen, (int, numpy.integer)):
        moidx[:mp.frozen] = False
    elif len(mp.frozen) > 0:
        moidx[list(mp.frozen)] = False
    else:
        raise NotImplementedError
    return moidx 

def get_nocc(mp):
    frozen = mp.frozen
    if mp._nocc is not None:
        return mp._nocc
    elif frozen is None:
        nocca = numpy.count_nonzero(mp.mo_occ[0] > 0)
        noccb = numpy.count_nonzero(mp.mo_occ[1] > 0)
    elif isinstance(frozen, (int, numpy.integer)):
        nocca = numpy.count_nonzero(mp.mo_occ[0] > 0) - frozen
        noccb = numpy.count_nonzero(mp.mo_occ[1] > 0) - frozen
        #assert(nocca > 0 and noccb > 0)
    elif isinstance(frozen[0], (int, numpy.integer, list, numpy.ndarray)):
        if len(frozen) > 0 and isinstance(frozen[0], (int, numpy.integer)):
            # The same frozen orbital indices for alpha and beta orbitals
            frozen = [frozen, frozen]
        occidxa = mp.mo_occ[0] > 0
        occidxa[list(frozen[0])] = False
        occidxb = mp.mo_occ[1] > 0
        occidxb[list(frozen[1])] = False
        nocca = numpy.count_nonzero(occidxa)
        noccb = numpy.count_nonzero(occidxb)
    else:
        raise NotImplementedError
    return nocca, noccb 

def get_nmo(mp):
    frozen = mp.frozen
    if mp._nmo is not None:
        return mp._nmo
    elif frozen is None:
        nmoa = mp.mo_occ[0].size
        nmob = mp.mo_occ[1].size
    elif isinstance(frozen, (int, numpy.integer)):
        nmoa = mp.mo_occ[0].size - frozen
        nmob = mp.mo_occ[1].size - frozen
    elif isinstance(frozen[0], (int, numpy.integer, list, numpy.ndarray)):
        if isinstance(frozen[0], (int, numpy.integer)):
            frozen = (frozen, frozen)
        nmoa = len(mp.mo_occ[0]) - len(set(frozen[0]))
        nmob = len(mp.mo_occ[1]) - len(set(frozen[1]))
    else:
        raise NotImplementedError
    return nmoa, nmob 

def make_hdiag(h1e, g2e, ci_strs, norb, nelec):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    strsa, strsb = ci_strs
    strsa = numpy.asarray(strsa)
    strsb = numpy.asarray(strsb)
    occslista = [[i for i in range(norb) if str0 & (1<<i)] for str0 in strsa]
    occslistb = [[i for i in range(norb) if str0 & (1<<i)] for str0 in strsb]

    g2e = ao2mo.restore(1, g2e, norb)
    diagj = numpy.einsum('iijj->ij',g2e)
    diagk = numpy.einsum('ijji->ij',g2e)
    hdiag = []
    for aocc in occslista:
        for bocc in occslistb:
            e1 = h1e[aocc,aocc].sum() + h1e[bocc,bocc].sum()
            e2 = diagj[aocc][:,aocc].sum() + diagj[aocc][:,bocc].sum() \
               + diagj[bocc][:,aocc].sum() + diagj[bocc][:,bocc].sum() \
               - diagk[aocc][:,aocc].sum() - diagk[bocc][:,bocc].sum()
            hdiag.append(e1 + e2*.5)
    return numpy.array(hdiag) 

def _trans2(fcivec, norb, nelec):
    if isinstance(nelec, (int, numpy.integer)):
        neleca = nelecb = nelec//2
    else:
        neleca, nelecb = nelec
    link_indexa = fci.cistring.gen_linkstr_index(range(norb), neleca)
    link_indexb = fci.cistring.gen_linkstr_index(range(norb), nelecb)
    na, nlinka = link_indexa.shape[:2]
    nb, nlinkb = link_indexb.shape[:2]
    fcivec = fcivec.reshape(na,nb)
    t1 = _trans1(fcivec, norb, nelec)
    t2 = numpy.zeros((na,nb,norb,norb,norb,norb))
    for str0, tab in enumerate(link_indexa):
        for a, i, str1, sign in tab:
            t2[str1,:,a,i] += sign * t1[str0]
    for k in range(na):
        for str0, tab in enumerate(link_indexb):
            for a, i, str1, sign in tab:
                t2[k,str1,a,i] += sign * t1[k,str0]
    return t2 

def _trans1(fcivec, norb, nelec):
    if isinstance(nelec, (int, numpy.integer)):
        neleca = nelecb = nelec//2
    else:
        neleca, nelecb = nelec
    link_indexa = fci.cistring.gen_linkstr_index(range(norb), neleca)
    link_indexb = fci.cistring.gen_linkstr_index(range(norb), nelecb)
    na, nlinka = link_indexa.shape[:2]
    nb, nlinkb = link_indexb.shape[:2]
    fcivec = fcivec.reshape(na,nb)
    t1 = numpy.zeros((na,nb,norb,norb))
    for str0, tab in enumerate(link_indexa):
        for a, i, str1, sign in tab:
            t1[str1,:,a,i] += sign * fcivec[str0]
    for k in range(na):
        for str0, tab in enumerate(link_indexb):
            for a, i, str1, sign in tab:
                t1[k,str1,a,i] += sign * fcivec[k,str0]
    return t1

#
# NOTE: this rdm3 is defined as
# rdm3(p,q,r,s,t,u) = <p^+ q r^+ s t^+ u> 

def contract_1e(f1e, fcivec, norb, nelec):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    link_indexa = cistring.gen_linkstr_index(range(norb), neleca)
    link_indexb = cistring.gen_linkstr_index(range(norb), nelecb)
    na = cistring.num_strings(norb, neleca)
    nb = cistring.num_strings(norb, nelecb)
    ci0 = fcivec.reshape(na,nb)
    t1 = numpy.zeros((norb,norb,na,nb))
    for str0, tab in enumerate(link_indexa):
        for a, i, str1, sign in tab:
            t1[a,i,str1] += sign * ci0[str0]
    for str0, tab in enumerate(link_indexb):
        for a, i, str1, sign in tab:
            t1[a,i,:,str1] += sign * ci0[:,str0]
    fcinew = numpy.dot(f1e.reshape(-1), t1.reshape(-1,na*nb))
    return fcinew.reshape(fcivec.shape) 

def make_hdiag(h1e, eri, norb, nelec, opt=None):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec

    occslista = cistring._gen_occslst(range(norb), neleca)
    occslistb = cistring._gen_occslst(range(norb), nelecb)
    eri = ao2mo.restore(1, eri, norb)
    diagj = numpy.einsum('iijj->ij', eri)
    diagk = numpy.einsum('ijji->ij', eri)
    hdiag = []
    for aocc in occslista:
        for bocc in occslistb:
            e1 = h1e[aocc,aocc].sum() + h1e[bocc,bocc].sum()
            e2 = diagj[aocc][:,aocc].sum() + diagj[aocc][:,bocc].sum() \
               + diagj[bocc][:,aocc].sum() + diagj[bocc][:,bocc].sum() \
               - diagk[aocc][:,aocc].sum() - diagk[bocc][:,bocc].sum()
            hdiag.append(e1 + e2*.5)
    return numpy.array(hdiag) 

def convert_to_uccsd(mycc):
    from pyscf import scf
    from pyscf.cc import uccsd, gccsd
    if isinstance(mycc, uccsd.UCCSD):
        return mycc
    elif isinstance(mycc, gccsd.GCCSD):
        raise NotImplementedError

    mf = scf.addons.convert_to_uhf(mycc._scf)
    ucc = uccsd.UCCSD(mf)
    assert(mycc._nocc is None)
    assert(mycc._nmo is None)
    ucc.__dict__.update(mycc.__dict__)
    ucc._scf = mf
    ucc.mo_coeff = mf.mo_coeff
    ucc.mo_occ = mf.mo_occ
    if not (mycc.frozen is None or isinstance(mycc.frozen, (int, numpy.integer))):
        raise NotImplementedError
    ucc.t1, ucc.t2 = uccsd.amplitudes_from_rccsd(mycc.t1, mycc.t2)
    return ucc 

def convert_to_gccsd(mycc):
    from pyscf import scf
    from pyscf.cc import gccsd
    if isinstance(mycc, gccsd.GCCSD):
        return mycc

    mf = scf.addons.convert_to_ghf(mycc._scf)
    gcc = gccsd.GCCSD(mf)
    assert(mycc._nocc is None)
    assert(mycc._nmo is None)
    gcc.__dict__.update(mycc.__dict__)
    gcc._scf = mf
    gcc.mo_coeff = mf.mo_coeff
    gcc.mo_occ = mf.mo_occ
    if isinstance(mycc.frozen, (int, numpy.integer)):
        gcc.frozen = mycc.frozen * 2
    elif not (mycc.frozen is None or mycc.frozen == 0):
        raise NotImplementedError
    gcc.t1 = spatial2spin(mycc.t1, mf.mo_coeff.orbspin)
    gcc.t2 = spatial2spin(mycc.t2, mf.mo_coeff.orbspin)
    return gcc 

def convert_to_gcisd(myci):
    from pyscf.ci import gcisd
    if isinstance(myci, gcisd.GCISD):
        return myci

    mf = scf.addons.convert_to_ghf(myci._scf)
    gci = gcisd.GCISD(mf)
    assert(myci._nocc is None)
    assert(myci._nmo is None)
    gci.__dict__.update(myci.__dict__)
    gci._scf = mf
    gci.mo_coeff = mf.mo_coeff
    gci.mo_occ = mf.mo_occ
    if isinstance(myci.frozen, (int, np.integer)):
        gci.frozen = myci.frozen * 2
    else:
        raise NotImplementedError
    gci.ci = gcisd.from_rcisdvec(myci.ci, myci.nocc, mf.mo_coeff.orbspin)
    return gci 

def to_fcivec(cisdvec, nelec, orbspin, frozen=None):
    assert(numpy.count_nonzero(orbspin == 0) ==
           numpy.count_nonzero(orbspin == 1))
    norb = len(orbspin)
    frozen_mask = numpy.zeros(norb, dtype=bool)
    if frozen is None:
        pass
    elif isinstance(frozen, (int, numpy.integer)):
        frozen_mask[:frozen] = True
    else:
        frozen_mask[frozen] = True
    frozen = (numpy.where(frozen_mask[orbspin == 0])[0],
              numpy.where(frozen_mask[orbspin == 1])[0])
    nelec = (numpy.count_nonzero(orbspin[:nelec] == 0),
             numpy.count_nonzero(orbspin[:nelec] == 1))
    orbspin = orbspin[~frozen_mask]
    nmo = len(orbspin)
    nocc = numpy.count_nonzero(~frozen_mask[:sum(nelec)])
    ucisdvec = to_ucisdvec(cisdvec, nmo, nocc, orbspin)
    return ucisd.to_fcivec(ucisdvec, norb//2, nelec, frozen) 

def from_fcivec(ci0, nelec, orbspin, frozen=None):
    if not (frozen is None or frozen == 0):
        raise NotImplementedError

    assert(numpy.count_nonzero(orbspin == 0) ==
           numpy.count_nonzero(orbspin == 1))
    norb = len(orbspin)
    frozen_mask = numpy.zeros(norb, dtype=bool)
    if frozen is None:
        pass
    elif isinstance(frozen, (int, numpy.integer)):
        frozen_mask[:frozen] = True
    else:
        frozen_mask[frozen] = True
    #frozen = (numpy.where(frozen_mask[orbspin == 0])[0],
    #          numpy.where(frozen_mask[orbspin == 1])[0])
    nelec = (numpy.count_nonzero(orbspin[:nelec] == 0),
             numpy.count_nonzero(orbspin[:nelec] == 1))
    ucisdvec = ucisd.from_fcivec(ci0, norb//2, nelec, frozen)
    nocc = numpy.count_nonzero(~frozen_mask[:sum(nelec)])
    return from_ucisdvec(ucisdvec, nocc, orbspin[~frozen_mask]) 

def make_rdm12(self, state, norb, nelec, link_index=None, **kwargs):
        nelectrons = 0
        if isinstance(nelec, (int, numpy.integer)):
          nelectrons = nelec
        else:
          nelectrons = nelec[0]+nelec[1]

        # The 2RDMs written by "save_spatial_twopdm_text" in BLOCK and STACKBLOCK
        # are written as E2[i1,j2,k2,l1]
        # and stored here as E2[i1,l1,j2,k2] (for PySCF purposes)
        # This is NOT done with SQA in mind.
        twopdm = numpy.zeros( (norb, norb, norb, norb) )
        file2pdm = "spatial_twopdm.%d.%d.txt" %(state, state)
        with open(os.path.join(self.scratchDirectory, "node0", file2pdm), "r") as f:
            norb_read = int(f.readline().split()[0])
            assert(norb_read == norb)
            for line in f:
                linesp = line.split()
                i, k, l, j = [int(x) for x in linesp[:4]]
                twopdm[i,j,k,l] = 2.0 * float(linesp[4])

        # (This is coherent with previous statement about indexes)
        onepdm = numpy.einsum('ikjj->ki', twopdm)
        onepdm /= (nelectrons-1)
        return onepdm, twopdm 

def parse_type(self, hdat, dataobj=None):
        dtype = super(MGHImageType, self).parse_type(hdat, dataobj=dataobj)
        if   np.issubdtype(dtype, np.floating): dtype = np.float32
        elif np.issubdtype(dtype, np.int8):     dtype = np.int8
        elif np.issubdtype(dtype, np.int16):    dtype = np.int16
        elif np.issubdtype(dtype, np.integer):  dtype = np.int32
        else: raise ValueError('Could not deduce appropriate MGH type for dtype %s' % dtype)
        return dtype 

def save_freesurfer_annot(filename, obj, index=None):
    '''
    save_freesurfer_annot(filename, prop) saves the given integer property prop to the given
      filename as a FreeSurfer annotation file.

    The optional argument index specifies how the colortab and names of the property labels should
    be handles. By default this is None, in which case names are generated using the formatter
    'label%d' and color values are generated using the label_colors function. If index is not None,
    then it is coerced to an index using the to_label_index() function, and the label index is
    used to create the colortable and names.
    '''
    # first parse the index
    if index is None:
        if len(obj) == 2 and pimms.is_vector(obj[0]): (obj, index) = obj
        else: index = label_index(obj)
    index = to_label_index(index)
    # okay, let's get the data in the right format
    (u,ris) = np.unique(obj, return_inverse=True)
    es = [index[l] for l in u]
    clrs = np.round([np.asarray(e.color)*255 for e in es]).astype('int')
    clrs[:,3] = 255 - clrs[:,3] # alpha -> transparency
    nms  = [e.name for e in es]
    lbls = [e.id for e in es]
    fsio.write_annot(filename, ris, clrs, nms, fill_ctab=True)
    return filename

# A few annot labels we can just save: 

def __modify_schema__(cls, field_schema: Dict[str, Any]) -> None:
        dt = cls._dtype
        if dt is int or np.issubdtype(dt, np.integer):
            items = {"type": "number", "multipleOf": 1.0}
        elif dt is float or np.issubdtype(dt, np.floating):
            items = {"type": "number"}
        elif dt is str or np.issubdtype(dt, np.string_):
            items = {"type": "string"}
        elif dt is bool or np.issubdtype(dt, np.bool_):
            items = {"type": "boolean"}
        field_schema.update(type="array", items=items)