Python lib.dot函数代码示例

 def h_op(x):
     x = self.unpack_uniq_var(x)
     norb = x.shape[0]
     #:hx = numpy.einsum('qp,xjj,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('qp,xqq,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('jk,xkk,xkp->pj', x, dip, dip)
     #:hx+= numpy.einsum('jk,xpp,xkp->pj', x, dip, dip)
     #:hx+= numpy.einsum('qj,xjq,xjp->pj', x, dip, dip)
     #:hx+= numpy.einsum('pk,xjp,xkp->pj', x, dip, dip)
     #:hx-= numpy.einsum('qp,xpp,xjq->pj', x, dip, dip) * 2
     #:hx-= numpy.einsum('qp,xjp,xpq->pj', x, dip, dip) * 2
     #:hx+= numpy.einsum('qj,xjp,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('pk,xkp,xjp->pj', x, dip, dip)
     #:hx-= numpy.einsum('jk,xjj,xkp->pj', x, dip, dip) * 2
     #:hx-= numpy.einsum('jk,xkj,xjp->pj', x, dip, dip) * 2
     #:return -self.pack_uniq_var(hx)
     #:hx = numpy.einsum('iq,qp->pi', g0, x)
     hx = lib.dot(x.T, g0.T).conj()
     #:hx+= numpy.einsum('qi,xiq,xip->pi', x, dip, dip) * 2
     hx+= numpy.einsum('xip,xi->pi', dip, numpy.einsum('qi,xiq->xi', x, dip)) * 2
     #:hx-= numpy.einsum('qp,xpp,xiq->pi', x, dip, dip) * 2
     hx-= numpy.einsum('xpp,xip->pi', dip,
                       lib.dot(dip.reshape(-1,norb), x).reshape(3,norb,norb)) * 2
     #:hx-= numpy.einsum('qp,xip,xpq->pi', x, dip, dip) * 2
     hx-= numpy.einsum('xip,xp->pi', dip, numpy.einsum('qp,xpq->xp', x, dip)) * 2
     return -self.pack_uniq_var(hx-hx.conj().T)

def general(mydf, mo_coeffs, kpts=None, compact=True):
    if isinstance(mo_coeffs, numpy.ndarray) and mo_coeffs.ndim == 2:
        mo_coeffs = (mo_coeffs,) * 4

    eri = mydf.get_eri(kpts)

####################
# gamma point, the integral is real and with s4 symmetry
    if eri.dtype == numpy.float64:
        return ao2mo.general(eri, mo_coeffs, compact=compact)
    else:
        mokl, klslice = ao2mo.incore._conc_mos(mo_coeffs[2], mo_coeffs[3],
                                               False)[2:]
        if mokl.dtype == numpy.float64:
            mokl = mokl + 0j
        nao = mo_coeffs[0].shape[0]
        nmoi = mo_coeffs[0].shape[1]
        nmoj = mo_coeffs[1].shape[1]
        nmok = mo_coeffs[2].shape[1]
        nmol = mo_coeffs[3].shape[1]
        moi = numpy.asarray(mo_coeffs[0], order='F')
        moj = numpy.asarray(mo_coeffs[1], order='F')
        tao = [0]
        ao_loc = None
        pqkl = _ao2mo.r_e2(eri.reshape(-1,nao**2), mokl, klslice, tao, ao_loc, aosym='s1')
        pqkl = pqkl.reshape(nao,nao,nmok*nmol)
        pjkl = numpy.empty((nao,nmoj,nmok*nmol), dtype=numpy.complex128)
        for i in range(nao):
            lib.dot(moj.T, pqkl[i], 1, pjkl[i], 0)
        pqkl = None
        eri_mo = lib.dot(moi.T.conj(), pjkl.reshape(nao,-1))
        return eri_mo.reshape(nmoi*nmoj,-1)

def get_j_kpts(mf, cell, dm_kpts, kpts, kpt_band=None):
    coords = gen_grid.gen_uniform_grids(cell)
    nkpts = len(kpts)
    ngs = len(coords)
    dm_kpts = np.asarray(dm_kpts)
    nao = dm_kpts.shape[-1]

    ni = numint._KNumInt(kpts)
    aoR_kpts = ni.eval_ao(cell, coords, kpts)
    if kpt_band is not None:
        aoR_kband = numint.eval_ao(cell, coords, kpt_band)

    dms = dm_kpts.reshape(-1,nkpts,nao,nao)
    nset = dms.shape[0]

    vjR = [get_vjR(cell, dms[i], aoR_kpts) for i in range(nset)]
    if kpt_band is not None:
        vj_kpts = [cell.vol/ngs * lib.dot(aoR_kband.T.conj()*vjR[i], aoR_kband)
                   for i in range(nset)]
        if dm_kpts.ndim == 3:  # One set of dm_kpts for KRHF
            vj_kpts = vj_kpts[0]
        return lib.asarray(vj_kpts)
    else:
        vj_kpts = []
        for i in range(nset):
            vj = [cell.vol/ngs * lib.dot(aoR_k.T.conj()*vjR[i], aoR_k)
                  for aoR_k in aoR_kpts]
            vj_kpts.append(lib.asarray(vj))
        return lib.asarray(vj_kpts).reshape(dm_kpts.shape)

def gamma1_intermediates(mycc, t1, t2, l1, l2):
    nocc, nvir = t1.shape
    doo =-numpy.einsum('ja,ia->ij', l1, t1)
    dvv = numpy.einsum('ia,ib->ab', l1, t1)
    dvo = l1.T
    xtv = numpy.einsum('ie,me->im', t1, l1)
    dov = t1 - numpy.einsum('im,ma->ia', xtv, t1)
    #:doo -= numpy.einsum('jkab,ikab->ij', l2, theta)
    #:dvv += numpy.einsum('jica,jicb->ab', l2, theta)
    #:xt1  = numpy.einsum('mnef,inef->mi', l2, make_theta(t2))
    #:xt2  = numpy.einsum('mnaf,mnef->ea', l2, make_theta(t2))
    #:dov += numpy.einsum('imae,me->ia', make_theta(t2), l1)
    #:dov -= numpy.einsum('ma,ie,me->ia', t1, t1, l1)
    #:dov -= numpy.einsum('mi,ma->ia', xt1, t1)
    #:dov -= numpy.einsum('ie,ae->ia', t1, xt2)
    max_memory = mycc.max_memory - lib.current_memory()[0]
    unit = nocc*nvir**2
    blksize = max(ccsd.BLKMIN, int(max_memory*.95e6/8/unit))
    for p0, p1 in prange(0, nocc, blksize):
        theta = make_theta(t2[p0:p1])
        doo[p0:p1] -= lib.dot(theta.reshape(p1-p0,-1), l2.reshape(nocc,-1).T)
        dov[p0:p1] += numpy.einsum('imae,me->ia', theta, l1)
        xt1 = lib.dot(l2.reshape(nocc,-1), theta.reshape(p1-p0,-1).T)
        dov[p0:p1] -= numpy.einsum('mi,ma->ia', xt1, t1)
        xt2 = lib.dot(theta.reshape(-1,nvir).T, l2[p0:p1].reshape(-1,nvir))
        dov -= numpy.einsum('ie,ae->ia', t1, xt2)
        dvv += lib.dot(l2[p0:p1].reshape(-1,nvir).T, theta.reshape(-1,nvir))
    return doo, dov, dvo, dvv

 def get_eri(self):
     nao = self.mol.nao_nr()
     nao_pair = nao * (nao+1) // 2
     ao_eri = numpy.zeros((nao_pair,nao_pair))
     for eri1 in self.loop():
         lib.dot(eri1.T, eri1, 1, ao_eri, 1)
     return ao2mo.restore(8, ao_eri, nao)

 def add_k_(vk, dm, pqk, coulG, buf=None):
     nG = pqk.shape[-1]
     pqk *= numpy.sqrt(coulG)
     ipk = lib.dot(dm, pqk.reshape(nao,-1)).reshape(nao,nao,-1)
     pik = numpy.ndarray((nao,nao,nG), buffer=buf)
     pik[:] = ipk.transpose(1,0,2)
     vk += lib.dot(pqk.reshape(nao,-1), pik.reshape(nao,-1).T)
     return vk

 def cost_function(self, u=None):
     if u is None: u = numpy.eye(self.mo_coeff.shape[1])
     mo_coeff = lib.dot(self.mo_coeff, u)
     dip = dipole_integral(self.mol, mo_coeff)
     r2 = self.mol.intor_symmetric('int1e_r2')
     r2 = numpy.einsum('pi,pi->', mo_coeff, lib.dot(r2, mo_coeff))
     val = r2 - numpy.einsum('xii,xii->', dip, dip) * 2
     return val

def transform_ci_for_orbital_rotation(ci, norb, nelec, u):
    '''Transform CI coefficients to the representation in new one-particle basis.
    Solving CI problem for Hamiltonian h1, h2 defined in old basis,
    CI_old = fci.kernel(h1, h2, ...)
    Given orbital rotation u, the CI problem can be either solved by
    transforming the Hamiltonian, or transforming the coefficients.
    CI_new = fci.kernel(u^T*h1*u, ...) = transform_ci_for_orbital_rotation(CI_old, u)

    Args:
        u : 2D array or a list of 2D array
            the orbital rotation to transform the old one-particle basis to new
            one-particle basis
    '''
    neleca, nelecb = _unpack(nelec)
    strsa = numpy.asarray(cistring.gen_strings4orblist(range(norb), neleca))
    strsb = numpy.asarray(cistring.gen_strings4orblist(range(norb), nelecb))
    one_particle_strs = numpy.asarray([1<<i for i in range(norb)])
    na = len(strsa)
    nb = len(strsb)

    if isinstance(u, numpy.ndarray) and u.ndim == 2:
        ua = ub = u
    else:
        ua, ub = u

    # Unitary transformation array trans_ci is the overlap between two sets of CI basis.
    occ_masks = (strsa[:,None] & one_particle_strs) != 0
    trans_ci_a = numpy.zeros((na,na))
    #for i in range(na): # for old basis
    #    for j in range(na):
    #        uij = u[occ_masks[i]][:,occ_masks[j]]
    #        trans_ci_a[i,j] = numpy.linalg.det(uij)
    occ_idx_all_strs = numpy.where(occ_masks)[1]
    for i in range(na):
        ui = ua[occ_masks[i]].T.copy()
        minors = numpy.take(ui, occ_idx_all_strs, axis=0).reshape(na,neleca,neleca)
        trans_ci_a[i,:] = numpy.linalg.det(minors)

    if neleca == nelecb and numpy.allclose(ua, ub):
        trans_ci_b = trans_ci_a
    else:
        occ_masks = (strsb[:,None] & one_particle_strs) != 0
        trans_ci_b = numpy.zeros((nb,nb))
        #for i in range(nb):
        #    for j in range(nb):
        #        uij = u[occ_masks[i]][:,occ_masks[j]]
        #        trans_ci_b[i,j] = numpy.linalg.det(uij)
        occ_idx_all_strs = numpy.where(occ_masks)[1]
        for i in range(nb):
            ui = ub[occ_masks[i]].T.copy()
            minors = numpy.take(ui, occ_idx_all_strs, axis=0).reshape(nb,nelecb,nelecb)
            trans_ci_b[i,:] = numpy.linalg.det(minors)

    # Transform old basis to new basis for all alpha-electron excitations
    ci = lib.dot(trans_ci_a.T, ci.reshape(na,nb))
    # Transform old basis to new basis for all beta-electron excitations
    ci = lib.dot(ci.reshape(na,nb), trans_ci_b)
    return ci

 def h_op(x):
     x = self.unpack_uniq_var(x)
     norb = x.shape[0]
     hx = lib.dot(x.T, g0.T)
     hx+= numpy.einsum('xip,xi->pi', pop, numpy.einsum('qi,xiq->xi', x, pop)) * 2
     hx-= numpy.einsum('xpp,xip->pi', pop,
                       lib.dot(pop.reshape(-1,norb), x).reshape(-1,norb,norb)) * 2
     hx-= numpy.einsum('xip,xp->pi', pop, numpy.einsum('qp,xpq->xp', x, pop)) * 2
     return -self.pack_uniq_var(hx-hx.T)

def gamma1_intermediates(mycc, t1, t2, l1, l2):
    nocc, nvir = t1.shape
    goo = -numpy.einsum('ja,ia->ij', l1, t1)
    gvv = numpy.einsum('ia,ib->ab', l1, t1)
    #:goo -= numpy.einsum('jkab,ikab->ij', l2, theta)
    #:gvv += numpy.einsum('jica,jicb->ab', l2, theta)
    theta = make_theta(t2)
    goo -= lib.dot(theta.reshape(nocc,-1), l2.reshape(nocc,-1).T)
    gvv += lib.dot(l2.reshape(-1,nvir).T, theta.reshape(-1,nvir))
    return goo, gvv

 def short_range_k(dm, Lpq, j3c):
     Lpq = Lpq.reshape(-1,nao)
     j3c = j3c.reshape(-1,nao)
     iLp = lib.dot(dm, Lpq.T).reshape(-1,nao)
     vk  = lib.dot(j3c.T, iLp)
     if hermi:
         vk = vk + vk.T
     else:
         iLp = lib.dot(dm.T, Lpq.T).reshape(-1,nao)
         vk += lib.dot(iLp.T, j3c)
     return vk

 def shift_grids(r):
     r_frac = lib.dot(r - origin, b.T)
     # Grids on the boundary (r_frac == +/-0.5) of the new cell may lead to
     # unbalanced contributions to the dipole moment. Exclude them from the
     # dipole and quadrupole
     r_frac[r_frac== 0.5] = 0
     r_frac[r_frac==-0.5] = 0
     r_frac[r_frac > 0.5] -= 1
     r_frac[r_frac <-0.5] += 1
     r = lib.dot(r_frac, a)
     return r

 def trans(aoiR, aojR, fac=1):
     if id(aoiR) == id(aojR) and ao2mo.incore.iden_coeffs(mo_coeffs[0], mo_coeffs[1]):
         moiR = mojR = numpy.asarray(lib.dot(mo_coeffs[0].T,aoiR.T), order='C')
     else:
         moiR = numpy.asarray(lib.dot(mo_coeffs[0].T, aoiR.T), order='C')
         mojR = numpy.asarray(lib.dot(mo_coeffs[1].T, aojR.T), order='C')
     mo_pairs_G = numpy.empty((nmoi,nmoj,ngs), dtype=numpy.complex128)
     for i in range(nmoi):
         mo_pairs_G[i] = tools.fft(fac * moiR[i].conj() * mojR, mydf.gs)
     mo_pairs_G = mo_pairs_G.reshape(-1,ngs).T
     return mo_pairs_G

 def trans(aoi, aoj, fac=1):
     if id(aoi) == id(aoj) and iden_coeffs(mo_coeffs[0], mo_coeffs[1]):
         moi = moj = numpy.asarray(lib.dot(mo_coeffs[0].T,aoi.T), order='C')
     else:
         moi = numpy.asarray(lib.dot(mo_coeffs[0].T, aoi.T), order='C')
         moj = numpy.asarray(lib.dot(mo_coeffs[1].T, aoj.T), order='C')
     mo_pairs_G = numpy.empty((nmoi,nmoj,ngrids), dtype=numpy.complex128)
     for i in range(nmoi):
         mo_pairs_G[i] = tools.fft(fac * moi[i].conj() * moj, mydf.mesh)
     mo_pairs_G = mo_pairs_G.reshape(-1,ngrids).T
     return mo_pairs_G

def general(mydf, mo_coeffs, kpts=None,
            compact=getattr(__config__, 'pbc_df_ao2mo_general_compact', True)):
    '''General MO integral transformation'''
    from pyscf.pbc.df.df_ao2mo import warn_pbc2d_eri
    warn_pbc2d_eri(mydf)
    cell = mydf.cell
    kptijkl = _format_kpts(kpts)
    kpti, kptj, kptk, kptl = kptijkl
    if isinstance(mo_coeffs, numpy.ndarray) and mo_coeffs.ndim == 2:
        mo_coeffs = (mo_coeffs,) * 4
    mo_coeffs = [numpy.asarray(mo, order='F') for mo in mo_coeffs]
    if not _iskconserv(cell, kptijkl):
        lib.logger.warn(cell, 'fft_ao2mo: momentum conservation not found in '
                        'the given k-points %s', kptijkl)
        return numpy.zeros([mo.shape[1] for mo in mo_coeffs])

    allreal = not any(numpy.iscomplexobj(mo) for mo in mo_coeffs)
    q = kptj - kpti
    coulG = tools.get_coulG(cell, q, mesh=mydf.mesh)
    coords = cell.gen_uniform_grids(mydf.mesh)
    max_memory = mydf.max_memory - lib.current_memory()[0]

    if gamma_point(kptijkl) and allreal:
        ao = mydf._numint.eval_ao(cell, coords, kpti)[0]
        if ((iden_coeffs(mo_coeffs[0], mo_coeffs[1]) and
             iden_coeffs(mo_coeffs[0], mo_coeffs[2]) and
             iden_coeffs(mo_coeffs[0], mo_coeffs[3]))):
            moiT = mojT = numpy.asarray(lib.dot(mo_coeffs[0].T,ao.T), order='C')
            ao = None
            max_memory = max_memory - moiT.nbytes*1e-6
            eri = _contract_compact(mydf, (moiT,mojT), coulG, max_memory=max_memory)
            if not compact:
                nmo = moiT.shape[0]
                eri = ao2mo.restore(1, eri, nmo).reshape(nmo**2,nmo**2)
        else:
            mos = [numpy.asarray(lib.dot(c.T, ao.T), order='C') for c in mo_coeffs]
            ao = None
            fac = numpy.array(1.)
            max_memory = max_memory - sum([x.nbytes for x in mos])*1e-6
            eri = _contract_plain(mydf, mos, coulG, fac, max_memory=max_memory).real
        return eri

    else:
        aos = mydf._numint.eval_ao(cell, coords, kptijkl)
        mos = [numpy.asarray(lib.dot(c.T, aos[i].T), order='C')
               for i,c in enumerate(mo_coeffs)]
        aos = None
        fac = numpy.exp(-1j * numpy.dot(coords, q))
        max_memory = max_memory - sum([x.nbytes for x in mos])*1e-6
        eri = _contract_plain(mydf, mos, coulG, fac, max_memory=max_memory)
        return eri