def numgrad(self, func, delta=1e-3):
     """
     # Numerical gradient
     """
     deltah = delta/2
     #
     # structure coefficients
     #
     structgrad = full.matrix(len(self.structs))
     #
     #
     for s in range(structgrad.shape[0]):
         self.coef[s] += deltah
         ep = func()
         self.coef[s] -= delta
         em = func()
         self.coef[s] += deltah
         structgrad[s] = (ep - em)/delta
     #
     # orbital gradient
     #
     r, c = Nod.C.shape
     orbgrad = full.matrix((r, c))
     for m in range(c):
         for t in range(r):
             Nod.C[t, m] += deltah
             ep = func()
             Nod.C[t, m] -= delta
             em = func()
             orbgrad[t, m] = (ep - em)/delta
             Nod.C[t, m] += deltah
     return (structgrad, orbgrad[:, :])

 def normhess(self):
     """
     Numerical norm Hessian
     #  N = <0|0>
     # d2<0|0> = <d20|0> + <0|d20> + 2<d0|d0> = 2<0|d20> + 2<d0|d0>
     """
     n_struct_hess = full.matrix(self.coef.shape*2)
     n_orb_structhess = full.matrix(Nod.C.shape + self.coef.shape)
     n_orb_hess = full.matrix(Nod.C.shape*2)
     #
     for s1, (str1, cstr1) in enumerate(zip(self.structs, self.coef)):
         for s2, (str2, cstr2) in enumerate(zip(self.structs, self.coef)):
             for det1, cdet1 in zip(str1.nods, str1.coef):
                 for det2, cdet2 in zip(str2.nods, str2.coef):
                     #
                     bk12 = BraKet(det1, det2)
                     #
                     n_struct_hess[s1, s2] += (cdet1*cdet2)*bk12.overlap()
                     n_orb_structhess[:, :, s1] += \
                         cdet1*cstr2*cdet2*bk12.overlap_gradient()
                     n_orb_hess += \
                         cstr1*cdet1*cstr2*cdet2*bk12.overlap_hessian()
     n_struct_hess *= 2
     n_orb_structhess *= 2
     #n_orb_hess
     return n_struct_hess, n_orb_structhess, n_orb_hess

    def energygrad(self):
        """
        #  N = <0|0>
        # dN = 2<d0|0>
        #  E = <0|H|0>/<0|0>
        # dE = 2<d0|H-E|0>/<0|0>
        """
        h_struct_grad = full.matrix(len(self.structs))
        h_orb_grad = full.matrix(Nod.C.shape)

        for s1, (str1, cstr1) in enumerate(zip(self.structs, self.coef)):
            for str2, cstr2 in zip(self.structs, self.coef):
                for det1, cdet1 in zip(str1.nods, str1.coef):
                    for det2, cdet2 in zip(str2.nods, str2.coef):
                        det12 = BraKet(det1, det2)

                        h_struct_grad[s1] += \
                            2*(cdet1*cstr2*cdet2)*\
                            det12.overlap()*\
                            det12.energy((self.h, self.h))

                        h_orb_grad += \
                            cstr1*cdet1*cstr2*cdet2*\
                            det12.energy_gradient((self.h, self.h))

        N = self.norm()
        E = self.energy()

        n_struct_grad, n_orb_grad = self.normgrad()

        structgrad = (1/N)*(h_struct_grad - E*n_struct_grad)
        orbgrad = (1/N)*(h_orb_grad - E*n_orb_grad)
        return (structgrad, orbgrad[:, :])

 def constraint(x, self):
     self.x = x
     mo = self.C[:, i]
     sg = full.matrix(self.coef.shape)
     og = full.matrix(self.C.shape)
     og[:, i] = 2*vb.Nod.S*mo
     return self.so2x(sg, og)

 def contravariant_transition_density_ao_mo(self):
     """Return contravariant density matrix in mix ao,mo basis"""
     d_am = (full.matrix(Nod.C.shape), full.matrix(Nod.C.shape))
     CK = self.K.orbitals()
     for s in (0, 1):
         if self.K(s) and self.L(s):
             d_am[s][:, self.L(s)] = CK[s]*self.transition_density[s]
     return d_am

 def contravariant_transition_density_mo_ao(self):
     """Return contravariant density matrix in mix mo,ao basis"""
     d_ma = (full.matrix(Nod.C.shape[::-1]), full.matrix(Nod.C.shape[::-1]))
     CL = self.L.orbitals()
     for s in (0, 1):
         if self.K(s) and self.L(s):
             d_ma[s][self.K(s), :] = self.transition_density[s]*CL[s].T
     return d_ma

 def fun(coef, wf):
     grad = full.matrix(coef.size)
     VBTestH2C.update_wf(coef, wf)
     mo = wf.C[:, i]
     dc2 = 2*vb.Nod.S*mo
     tmp = full.matrix(wf.C.shape)
     tmp[:, i] = dc2
     grad[wf.coef.size:] = tmp.block(*wf.blockdims).ravel(order='F')
     return grad

    def QUab(self):
        """Quadrupole moment"""
        if self._QUab is not None: return self._QUab

        D = self.D
        R = self.R
        Rc = self.Rc
        dRab = self.dRab
        Qab = self.Qab
        Dab = self.Dab

        lab = ("XXSECMOM", "XYSECMOM", "XZSECMOM", 
                           "YYSECMOM", "YZSECMOM", 
                                       "ZZSECMOM")

        xy = self.getprop(*lab)

        noa = self.noa
        QUab = full.matrix((6, noa, noa))
        rrab = full.matrix((6, noa, noa))
        rRab = full.matrix((6, noa, noa))
        RRab = full.matrix((6, noa, noa))
        Rab = self.Rab
        for a in range(noa):
            for b in range(noa):
                ij = 0
                for i in range(3):
                    for j in range(i, 3):
                        rrab[ij, a, b] = -(
                            xy[ij].subblock[a][b]&D.subblock[a][b]
                            ) 
                        rRab[ij, a, b] = Dab[i, a, b]*Rab[a,b,j]+Dab[j, a, b]*Rab[a,b,i]
                        RRab[ij, a, b] = Rab[a,b,i]*Rab[a,b,j]*Qab[a, b]
                        ij += 1
        QUab = rrab-rRab-RRab
        self._QUab = QUab
        #
        # Addition term - gauge correction summing up bonds
        #
        dQUab = full.matrix(self.QUab.shape)
        for a in range(noa):
            for b in range(noa):
                ij = 0
                for i in range(3):
                    for j in range(i, 3):
                        dQUab[ij, a, b] = dRab[a, b, i]*Dab[j, a, b] \
                                      +dRab[a, b, j]*Dab[i, a, b]
                        ij += 1
        self.dQUab = - dQUab

        return self._QUab

    def project_virtual_occupied(self, h1):
        """Rhs derivative <K|h|dL/dC(mu, m)>"""

        D_mo = self.transition_density
        delta = self.co_contravariant_transition_delta()

        K_h_dL = (full.matrix(Nod.C.shape), full.matrix(Nod.C.shape))

        CK = self.K.orbitals()
        for s in (0, 1):
            if self.K(s) and self.L(s):
                K_h_dL[s][:, self.L(s)] += \
                    delta[s]*h1[s].T*CK[s]*D_mo[s]*self.overlap()
        return K_h_dL

 def full_mo_transition_density(self):
     """
     Return mo transition density matrix in full mo basis
     """
     if self._ftd is None:
         _, mo = Nod.C.shape
         D_KL = self.transition_density
         self._ftd = (full.matrix((mo, mo)), full.matrix((mo, mo)))
         for s in (0, 1):
             if self.K(s) and self.L(s):
                 D_KL[s].scatter(
                     self._ftd[s], rows=self.K(s), columns=self.L(s)
                     )
     return self._ftd

 def test_quadrupole_allbonds(self):
     QU = full.matrix(ref.QU.shape)
     QUab = self.m.QUab
     for ab, a, b in pairs(self.m.noa):
         QU[:, ab] += QUab[:, a, b ] 
         if a != b: QU[:, ab] += QUab[:, b, a] 
     self.assert_allclose(QU, ref.QU)

def vb_transform(dens, delta, **kwargs):
    filename = kwargs.get('filename', '/tmp/AOTWOINT')
    a, m = dens[0].shape
    H_uvmn = matrix((a, a, m, m))
    for ig, g in list_integrals(filename):
        p, q, r, s = ig
        s, r, q, p = p-1, q-1, r-1, s-1
        if p == q: g *= 0.5
        if r == s: g *= 0.5
        if (p, q) == (r, s): g *= 0.5

        for d1, D1 in zip(dens, delta):
            for d2, D2 in zip(dens, delta):
                H_uvmn += D1[:, q].x(D2[:, s].x(d1[p, :].x(d2[r, :]*g)))
                H_uvmn += D1[:, p].x(D2[:, s].x(d1[q, :].x(d2[r, :]*g)))
                H_uvmn += D1[:, q].x(D2[:, r].x(d1[p, :].x(d2[s, :]*g)))
                H_uvmn += D1[:, p].x(D2[:, r].x(d1[q, :].x(d2[s, :]*g)))
                H_uvmn += D1[:, s].x(D2[:, q].x(d1[r, :].x(d2[p, :]*g)))
                H_uvmn += D1[:, r].x(D2[:, q].x(d1[s, :].x(d2[p, :]*g)))
                H_uvmn += D1[:, s].x(D2[:, p].x(d1[r, :].x(d2[q, :]*g)))
                H_uvmn += D1[:, r].x(D2[:, p].x(d1[s, :].x(d2[q, :]*g)))

            H_uvmn -= D1[:, s].x(D1[:, q].x(d1[p, :].x(d1[r, :]*g)))
            H_uvmn -= D1[:, s].x(D1[:, p].x(d1[q, :].x(d1[r, :]*g)))
            H_uvmn -= D1[:, r].x(D1[:, q].x(d1[p, :].x(d1[s, :]*g)))
            H_uvmn -= D1[:, r].x(D1[:, p].x(d1[q, :].x(d1[s, :]*g)))
            H_uvmn -= D1[:, q].x(D1[:, s].x(d1[r, :].x(d1[p, :]*g)))
            H_uvmn -= D1[:, q].x(D1[:, r].x(d1[s, :].x(d1[p, :]*g)))
            H_uvmn -= D1[:, p].x(D1[:, s].x(d1[r, :].x(d1[q, :]*g)))
            H_uvmn -= D1[:, p].x(D1[:, r].x(d1[s, :].x(d1[q, :]*g)))
        
    return H_uvmn.transpose(0, 2, 1, 3)

def vb_transform2(Dma, Dam, Delta1, Delta2, **kwargs):
    filename = kwargs.get('filename', '/tmp/AOTWOINT')
    H_umvn = matrix(Dam[0].shape + Dam[0].shape)
    for ig, g in list_integrals(filename):
        p, q, r, s = ig
        s, r, q, p = p-1, q-1, r-1, s-1
        if p == q: g *= 0.5
        if r == s: g *= 0.5
        if (p, q) == (r, s): g *= 0.5

        for d1a, d2a, D1a, D2a in zip(Dma, Dam, Delta1, Delta2):
            for d1b, d2b, D1b, D2b in zip(Dma, Dam, Delta1, Delta2):
               H_umvn += D1a[p,:].x(d1a[:, q].x(D2b[:, s].x(d2b[r, :]*g)))
               H_umvn += D1a[q,:].x(d1a[:, p].x(D2b[:, s].x(d2b[r, :]*g)))
               H_umvn += D1a[p,:].x(d1a[:, q].x(D2b[:, r].x(d2b[s, :]*g)))
               H_umvn += D1a[q,:].x(d1a[:, p].x(D2b[:, r].x(d2b[s, :]*g)))

               H_umvn += D1a[r,:].x(d1a[:, s].x(D2b[:, q].x(d2b[p, :]*g)))
               H_umvn += D1a[r,:].x(d1a[:, s].x(D2b[:, p].x(d2b[q, :]*g)))
               H_umvn += D1a[s,:].x(d1a[:, r].x(D2b[:, q].x(d2b[p, :]*g)))
               H_umvn += D1a[s,:].x(d1a[:, r].x(D2b[:, p].x(d2b[q, :]*g)))

            H_umvn -= D1a[p, :].x(d1a[:, s].x(D2a[:, q].x(d2a[r, :]*g)))
            H_umvn -= D1a[q, :].x(d1a[:, s].x(D2a[:, p].x(d2a[r, :]*g)))
            H_umvn -= D1a[p, :].x(d1a[:, r].x(D2a[:, q].x(d2a[s, :]*g)))
            H_umvn -= D1a[q, :].x(d1a[:, r].x(D2a[:, p].x(d2a[s, :]*g)))

            H_umvn -= D1a[r, :].x(d1a[:, q].x(D2a[:, s].x(d2a[p, :]*g)))
            H_umvn -= D1a[r, :].x(d1a[:, p].x(D2a[:, s].x(d2a[q, :]*g)))
            H_umvn -= D1a[s, :].x(d1a[:, q].x(D2a[:, r].x(d2a[p, :]*g)))
            H_umvn -= D1a[s, :].x(d1a[:, p].x(D2a[:, r].x(d2a[q, :]*g)))

    return H_umvn

    def Bab(self):
        """Localized hyperpolariziabilities"""
        if self._Bab is not None: return self._Bab

        D2k = self.D2k
        Rab = self.Rab
        d2Qa = self.d2Qa
        x = self.x

        labs = ('XDIPLEN ', 'YDIPLEN ', 'ZDIPLEN ')
        qlabs = [labs[i] + labs[j] for i in range(3) for j in range(i,3)]
        Bab = full.matrix( (self.nfreqs, 3, 6, self.noa, self.noa) )


        #correction term for shifting origin from O to Rab
        for i, li in enumerate(labs):
            for jk,ljk in enumerate(qlabs):
                for a in range(self.noa):
                    for b in range(self.noa):
                        for iw, w in enumerate(self.freqs):
                            Bab[iw, i, jk, a, b] = (
                                -x[i].subblock[a][b] & D2k [(ljk, w, w)].subblock[a][b]
                                )
                    for iw in self.rfreqs:
                        Bab[iw, i, jk, a, a] -= d2Qa[iw, a, jk]*Rab[a, a, i]

        self._Bab = Bab
        return self._Bab

 def test_dipole_allbonds(self):
     D = full.matrix(ref.D.shape)
     Dab = self.m.Dab
     for ab, a, b in pairs(self.m.noa):
         D[:, ab] += Dab[:, a, b ] 
         if a != b: D[:, ab] += Dab[:, b, a] 
     self.assert_allclose(D, ref.D)

def Dao(K, L):
    """
    # Return intermediate normalized ao transition density matrix given
    # determinants K and L
    # as [Dalpha, Dbeta]
    """

    if abs(K*L) < SINGULAR_OVERLAP_THRESHOLD:
        raise SingularOverlapError

    CK = K.orbitals()
    CL = L.orbitals()
    #
    D = []

    for s in range(2):
        if CK[s] is None or CL[s] is None:
            #
            # if None orbitals set atomic density to zero matrix
            #
            D.append(full.matrix(Nod.S.shape))
        else:
            SLK = CL[s].T*Nod.S*CK[s]
            try:
                D.append(CK[s]*(CL[s].T/SLK))
            except LinAlgError:
                raise SingularOverlapError

    return D

 def constraint_norm_grad(coef, wf):
     grad = full.matrix(coef.size)
     VBTestH2C.update_wf(coef, wf)
     sg, og = wf.normgrad()
     grad[:wf.coef.size] = sg
     grad[wf.coef.size:] = og.block(*wf.blockdims).ravel(order='F')
     return grad

 def wf_gradient(coef, wf):
     grad = full.matrix(coef.size)
     VBTestH2C.update_wf(coef, wf)
     sg, og = wf.energygrad()
     grad[:wf.coef.size] = sg
     grad[wf.coef.size:] = og.block(*wf.blockdims).ravel(order='F')
     return grad

    def project_occupied_virtual(self, h1):
        """Lhs derivative <dK/dC(mu,m)|h|L>"""

        D_mo = self.transition_density
        delta = self.contra_covariant_transition_delta()

        dK_h_La = full.matrix(Nod.C.shape[::-1])
        dK_h_Lb = full.matrix(Nod.C.shape[::-1])

        CL = self.L.orbitals()
        if self.L(0):
            dK_h_La[self.K(0), :] += \
                D_mo[0]*CL[0].T*h1[0].T*delta[0]*self.overlap()
        if self.L(1):
            dK_h_Lb[self.K(1), :] += \
                D_mo[1]*CL[1].T*h1[1].T*delta[1]*self.overlap()

        return dK_h_La.T, dK_h_Lb.T

 def grad(x, self):
     _x = self.x
     self.x = x
     p = self.p
     dL = full.matrix(len(x))
     dL[:len(p)] = self.gp(p) - sum(l*c['jac'](p) for l, c in zip(self.l, self.cp))
     dL[len(p):] = [-c['fun'](p) for c in self.cp]
     self.x = _x
     return dL