Python numpy.integer() Examples

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 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_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 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 _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 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 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 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 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 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 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 _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 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 _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 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_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 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_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 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_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 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 _response_dm1(mycc, Xvo, eris=None):
    nvir, nocc = Xvo.shape
    nmo = nocc + nvir
    with_frozen = not ((mycc.frozen is None)
                       or (isinstance(mycc.frozen, (int, numpy.integer)) and mycc.frozen == 0)
                       or (len(mycc.frozen) == 0))
    if eris is None or with_frozen:
        mo_energy = mycc._scf.mo_energy
        mo_occ = mycc.mo_occ
        mo_coeff = mycc.mo_coeff
        def fvind(x):
            x = x.reshape(Xvo.shape)
            dm = reduce(numpy.dot, (mo_coeff[:,nocc:], x, mo_coeff[:,:nocc].T))
            v = mycc._scf.get_veff(mycc.mol, dm + dm.T)
            v = reduce(numpy.dot, (mo_coeff[:,nocc:].T, v, mo_coeff[:,:nocc]))
            return v * 2
    else:
        mo_energy = eris.mo_energy
        mo_occ = numpy.zeros_like(mo_energy)
        mo_occ[:nocc] = 2
        ovvo = numpy.empty((nocc,nvir,nvir,nocc))
        for i in range(nocc):
            ovvo[i] = eris.ovvo[i]
            ovvo[i] = ovvo[i] * 4 - ovvo[i].transpose(1,0,2)
            ovvo[i]-= eris.oovv[i].transpose(2,1,0)
        def fvind(x):
            return numpy.einsum('iabj,bj->ai', ovvo, x.reshape(Xvo.shape))
    dvo = cphf.solve(fvind, mo_energy, mo_occ, Xvo, max_cycle=30)[0]
    dm1 = numpy.zeros((nmo,nmo))
    dm1[nocc:,:nocc] = dvo
    dm1[:nocc,nocc:] = dvo.T
    return dm1 

def make_rdm12(self, fakewfn_by_rdm2, ncas, nelec, **kwargs):
        if not isinstance(nelec, (int, numpy.integer)):
            nelec = sum(nelec)
# CheMPS2 uses physics notation
        rdm2 = fakewfn_by_rdm2.transpose(0,2,1,3)
        rdm1 = numpy.einsum('ijkk->ij', rdm2) / (nelec-1)
        return rdm1, rdm2 

def make_rdm1(civec_strs, norb, nelec):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    ci_coeff, ci_strs = civec_strs
    strsa, strsb = ci_strs
    strsa = numpy.asarray(strsa)
    strsb = numpy.asarray(strsb)

    cd_indexa = cre_des_linkstr(strsa, norb, neleca)
    cd_indexb = cre_des_linkstr(strsb, norb, nelecb)
    na = len(strsa)
    nb = len(strsb)

    fcivec = ci_coeff.reshape(na,nb)
    rdm1 = numpy.zeros((norb,norb))
    for str1, tab in enumerate(cd_indexa):
        for a, i, str0, sign in tab:
            if a >= 0:
                rdm1[a,i] += sign * numpy.dot(fcivec[str1], fcivec[str0])

    for str1, tab in enumerate(cd_indexb):
        for a, i, str0, sign in tab:
            if a >= 0:
                rdm1[a,i] += sign * numpy.dot(fcivec[:,str1], fcivec[:,str0])
    return rdm1.T

# dm_pq,rs = <|p^+ q r^+ s|>