def numeric_forces(atoms,D=None,**opts):
    "Compute numerical forces on atoms"
    # D is ignored here.
    dx = opts.get('dx',1e-6)
    sym = opts.get('sym',True)
    return_energy = opts.get('return_energy',False)
    nat = len(atoms)
    Forces = zeros((nat,3),'d')
    E0 = scf(atoms)
    for iat in range(nat):
        for idir in range(3):
            dr = zeros(3,'d')
            dr[idir] = dx
            atoms[iat].translate(dr)
            Ep = scf(atoms)
            atoms[iat].translate(-dr)
            if sym:
                atoms[iat].translate(-dr)
                Em = scf(atoms)
                atoms[iat].translate(dr)
                Forces[iat,idir] = 0.5*(Ep-Em)/dx
            else:
                Forces[iat,idir] = (Ep-E0)/dx
    if return_energy: return E0,Forces
    return Forces

def numeric_forces(atoms,D=None,**kwargs):
    "Compute numerical forces on atoms"
    # D is ignored here.
    dx = kwargs.get('dx',settings.NumericForceDx)
    sym = kwargs.get('sym',settings.NumericForceSym)
    return_energy = kwargs.get('return_energy')
    nat = len(atoms)
    Forces = zeros((nat,3),'d')
    E0 = scf(atoms)
    for iat in xrange(nat):
        for idir in xrange(3):
            dr = zeros(3,'d')
            dr[idir] = dx
            atoms[iat].translate(dr)
            Ep = scf(atoms)
            atoms[iat].translate(-dr)
            if sym:
                atoms[iat].translate(-dr)
                Em = scf(atoms)
                atoms[iat].translate(dr)
                Forces[iat,idir] = 0.5*(Ep-Em)/dx
            else:
                Forces[iat,idir] = (Ep-E0)/dx
    if return_energy: return E0,Forces
    return Forces

 def zero_density(self):
     """Initialize the density matrices for later. If nothing else,
     this assures that I don't have to worry about an undefined
     gamma matrix later on."""
     self.density = zeros((self.ng,2),'d')
     self.grada = zeros((self.ng,3),'d')
     self.gradb = zeros((self.ng,3),'d')
     self.gamma = zeros((self.ng,3),'d')

def TransformInts(Ints,orbs):
    """O(N^5) 4-index transformation of the two-electron integrals. Not as
    efficient as it could be, since it inflates to the full rectangular
    matrices rather than keeping them compressed. But at least it gets the
    correct result."""

    from time import time

    t0 = time()

    nbf,nmo = orbs.shape
    totlen = nmo*(nmo+1)*(nmo*nmo+nmo+2)/8

    temp = zeros((nbf,nbf,nbf,nmo),'d')
    tempvec = zeros(nbf,'d')
    temp2 = zeros((nbf,nbf,nmo,nmo),'d')

    mos = range(nmo) # preform so we don't form inside loops
    bfs = range(nbf)

    # Start with (mu,nu|sigma,eta)
    # Unpack aoints and transform eta -> l
    for mu in bfs:
        for nu in bfs:
            for sigma in bfs:
                for l in mos:
                    for eta in bfs:
                        tempvec[eta] = Ints[ijkl2intindex(mu,nu,sigma,eta)]
                    temp[mu,nu,sigma,l] = dot(orbs[:,l],tempvec)

    # Transform sigma -> k
    for mu in bfs:
        for nu in bfs:
            for l in mos:
                for k in mos:
                    temp2[mu,nu,k,l] = dot(orbs[:,k],temp[mu,nu,:,l])

    # Transform nu -> j
    for mu in bfs:
        for k in mos:
            for l in mos:
                for j in mos:
                    temp[mu,j,k,l] = dot(orbs[:,j],temp2[mu,:,k,l])

    # Transform mu -> i and repack integrals:
    MOInts = zeros(totlen,'d')
    for i in mos:
        for j in xrange(i+1):
            ij = i*(i+1)/2+j
            for k in mos:
                for l in xrange(k+1):
                    kl = k*(k+1)/2+l
                    if ij >= kl:
                        ijkl = ijkl2intindex(i,j,k,l)
                        MOInts[ijkl] = dot(orbs[:,i],temp[:,j,k,l])

    del temp,temp2,tempvec #force garbage collection now
    return MOInts

def TransformIntsMP2(Ints,orbs,nclosed):
    """\
    O(N^5) 4-index transformation of the two-electron integrals.
    Only transform the ones needed for MP2, which reduces the
    scaling to O(nN^4), where n are the occs (<<N).
    """

    from time import time

    t0 = time()

    nbf,nmo = orbs.shape
    totlen = nmo*(nmo+1)*(nmo*nmo+nmo+2)/8

    occs = range(nclosed)
    mos = range(nmo)
    bfs = range(nbf)

    # Start with (mu,nu|sigma,eta)
    # Unpack aoints and transform sigma -> b
    temp = zeros((nbf,nbf,nclosed,nbf),'d')
    tempvec = zeros(nbf,'d')
    for mu in bfs:
        for nu in bfs:
            for eta in bfs:
                for b in occs:
                    for sigma in bfs:
                        tempvec[sigma] = Ints[ijkl2intindex(mu,nu,sigma,eta)]
                    temp[mu,nu,b,eta] = dot(orbs[:,b],tempvec)

    temp2 = zeros((nclosed,nbf,nclosed,nbf),'d')
    for nu in bfs:
        for eta in bfs:
            for b in occs:
                for a in occs:
                    temp2[a,nu,b,eta] = dot(orbs[:,a],temp[:,nu,b,eta])

    temp = zeros((nclosed,nbf,nclosed,nmo),'d')
    for a in occs:
        for nu in bfs:
            for b in occs:
                for j in mos:
                    temp[a,nu,b,j] = dot(orbs[:,j],temp2[a,nu,b,:])

    # Transform mu -> i and repack integrals:
    MOInts = zeros(totlen,'d')
    for a in occs:
        for j in mos:
            for b in occs:
                for i in mos:
                    aibj = ijkl2intindex(a,i,b,j)
                    MOInts[aibj] = dot(orbs[:,i],temp[a,:,b,j])

    #print "Integral transform time = ",time()-t0
    del temp,temp2,tempvec #force garbage collection now
    return MOInts

def EXXC1(dens,gamma):
    "AEM's EXX compatible correlation #1 (note: no spin). AEM June 2006."
    npts = len(dens)
    ec = zeros(npts,'d')
    vc = zeros(npts,'d')
    for i in range(npts):
        rho = float(dens[i]) # Density 
        gam = gamma[i]
        ecpnt,dnedrho,dnedgamma = c1(rho,gam)
        ec[i] = ecpnt
        vc[i] = dnedrho   # more derivatives need to be added
    return ec,vc

def CPBE(dens,gamma):
    "PBE Correlation Functional"
    npts = len(dens)
    ec = zeros(npts,'d')
    vc = zeros(npts,'d')
    for i in range(npts):
        rho = 0.5*float(dens[i]) # Density of the alpha spin
        gam = 0.25*gamma[i]
        ecab,vca,vcb = cpbe(rho,rho,gam,gam,gam)
        ec[i] = ecab
        vc[i] = vca
    return ec,vc

def XPBE(dens,gamma):
    "PBE Exchange Functional"
    npts = len(dens)
    assert len(gamma) == npts
    ex = zeros(npts,'d')
    vx = zeros(npts,'d')
    for i in range(npts):
        rho = 0.5*float(dens[i])  # Density of the alpha spin
        gam = 0.25*gamma[i]
        exa,vxa = xpbe(rho,gam)
        ex[i] = 2*exa
        vx[i] = vxa
    return ex,vx

 def set_bf_amps(self,bfs):
     x,y,z,w = self.xyzw()
     nbf = len(bfs)
     self.bfs = zeros(nbf,'d')
     for i in xrange(nbf):
         self.bfs[i] = bfs[i].amp(x,y,z)
     # This *if* statement is potentially slow. If it becomes
     #  a factor, pull it up to AtomicGrids and have two
     #  explicit cases here.
     if self.do_grad_dens:
         self.bfgrads = zeros((nbf,3),'d')
         for i in xrange(nbf):
             self.bfgrads[i,:] = bfs[i].grad(x,y,z)
     return

def get1ints(bfs,atoms):
    "Form the overlap S and h=t+vN one-electron Hamiltonian matrices"
    nbf = len(bfs)
    S = zeros((nbf,nbf),'d')
    h = zeros((nbf,nbf),'d')

    for i in xrange(nbf):
        bfi = bfs[i]
        for j in xrange(nbf):
            bfj = bfs[j]
            S[i,j] = bfi.overlap(bfj)
            h[i,j] = bfi.kinetic(bfj)
            for atom in atoms:
                h[i,j] = h[i,j] + atom.atno*bfi.nuclear(bfj,atom.pos())
    return S,h

def LYP(dens,gamma):
    """Transformed version of LYP. See 'Results obtained with correlation
    energy density functionals of Becke and Lee, Yang, and Parr.' Miehlich,
    Savin, Stoll and Preuss. CPL 157, 200 (1989).
    """
    npts = len(dens)
    ec = zeros(npts,'d')
    vc = zeros(npts,'d')
    for i in range(npts):
        rho = 0.5*float(dens[i]) # Density of the alpha spin
        gam = 0.25*gamma[i]
        ecab,vca,vcb = clyp(rho,rho,gam,gam,gam)
        ec[i] = ecab
        vc[i] = vca
    return ec,vc

 def add_basis(self,bfs):
     # Compute the amplitudes of the basis functions over the grid
     self.nbf = len(bfs)
     self.bfgrid = zeros((self.ng,self.nbf),'d')
     for ibf in xrange(self.nbf):
         for ig in xrange(self.ng):
             x,y,z,w = self.xyzw[ig,:]
             self.bfgrid[ig,ibf] = bfs[ibf].amp(x,y,z)
     if self.do_grad_dens:
         self.bfgrads = zeros((self.ng,self.nbf,3),'d')
         for ibf in xrange(self.nbf):
             for ig in xrange(self.ng):
                 x,y,z,w = self.xyzw[ig,:]
                 self.bfgrads[ig,ibf,:] = bfs[ibf].grad(x,y,z)
     return

def TransformInts(Ints,orbs1,orbs2, nocc):

    nbf,nmo = orbs1.shape
    totlen = nmo*nmo*nmo*nmo

    occs = range(nocc)
    mos = range(nmo)
    bfs = range(nbf)

    # Start with (mu,nu|sigma,eta)
    # Unpack aoints and transform sigma -> b
    # Here sigma, b, nu, j are of first same spin group,
    #  others of second same spin group
    temp = zeros((nbf,nbf,nbf,nbf),'d')
    tempvec = zeros(nbf,'d')
    for mu in bfs:
        for nu in bfs:
            for eta in bfs:
                for b in bfs:
                    for sigma in bfs:
                        tempvec[sigma] = Ints[ijkl2intindex(mu,nu,sigma,eta)]
                    temp[mu,nu,b,eta] = dot(orbs1[b,:],tempvec)

    temp2 = zeros((nbf,nbf,nbf,nbf),'d')
    for nu in bfs:
        for eta in bfs:
            for b in bfs:
                for a in bfs:
                    temp2[a,nu,b,eta] = dot(orbs2[a,:],temp[:,nu,b,eta])

    temp = zeros((nbf,nbf,nbf,nbf),'d')
    for a in bfs:
        for nu in bfs:
            for b in bfs:
                for j in bfs:
                    temp[a,nu,b,j] = dot(orbs1[j,:],temp2[a,nu,b,:])

    # Transform mu -> i and repack integrals:
    MOInts = zeros(totlen,'d')
    for a in bfs:
        for j in bfs:
            for b in bfs:
                for i in bfs:
                    aibj = ijkl2intindex(a,i,b,j)
                    MOInts[aibj] = dot(orbs2[i,:],temp[a,:,b,j])

    del temp,temp2,tempvec #force garbage collection now
    return MOInts, nbf

def get_F0_old(atoms):
    "Form the zero-iteration (density matrix independent) Fock matrix"
    nbf = get_nbf(atoms)
    nat = len(atoms)
    
    F0 = zeros((nbf,nbf),'d')

    ibf = 0 # bf number of the first bfn on iat
    for iat in range(nat):
        atomi = atoms[iat]
        for i in range(atomi.nbf):
            bfi = atomi.basis[i]
            F0[ibf+i,ibf+i] = bfi.u
            
            jbf = 0
            for jat in range(nat):
                atomj = atoms[jat]
                if iat != jat:
                    gammaij = get_gamma(atomi,atomj)
                    betaij = get_beta0(atomi.atno,atomj.atno)
                    F0[ibf+i,ibf+i] -= gammaij*atomj.Z
                    for j in range(atomj.nbf):
                        bfj = atomj.basis[j]
                        Sij = bfi.cgbf.overlap(bfj.cgbf)
                        #Sij = mopac_overlap(bfi,bfj)
                        IPij = bfi.ip+bfj.ip
                        F0[ibf+i,jbf+j] = betaij*IPij*Sij
                        F0[jbf+j,ibf+i] = F0[ibf+i,jbf+j]
                jbf += atomj.nbf
        ibf += atomi.nbf
    return F0

 def Ffunc_num(cnew):
     E0 = Efunc(cnew)
     F = zeros(3*nat,'d')
     ei = zeros(3*nat,'d')
     dx = 1e-7
     for i in range(nat):
         for j in range(3):
             ei[3*i+j] = 1.0
             E1 = Efunc(cnew+ei*dx)
             ei[3*i+j] = 0.0
             F[3*i+j] = (E1-E0)/dx
     if verbose_level > 0:
         print "MINDO3 gradient calculation requested:"
         print atoms
         print Hf
     return F

def get_F2_open(atoms,Da,Db):
    "Two-center corrections to the core fock matrix"
    nbf = get_nbf(atoms)
    nat = len(atoms)
    
    F2 = zeros((nbf,nbf),'d')

    ibf = 0 # bf number of the first bfn on iat
    for iat in range(nat):
        atomi = atoms[iat]
        jbf = 0
        for jat in range(nat):
            atomj = atoms[jat]
            if iat != jat:
                gammaij = get_gamma(atomi,atomj)
                for i in range(atomi.nbf):
                    for j in range(atomj.nbf):
                        pija = Da[ibf+i,jbf+j] 
                        pijb = Db[ibf+i,jbf+j] 
                        pij = pija+pijb
                        qja = Da[jbf+j,jbf+j]
                        qjb = Db[jbf+j,jbf+j]
                        qj = qja+qjb
                        qia = Da[ibf+i,ibf+i]
                        qib = Db[ibf+i,ibf+i]
                        qi = qia+qib
                        F2[ibf+i,jbf+j] -= 0.25*pij*gammaij
                        F2[jbf+j,ibf+i] = F2[ibf+i,jbf+j]
                        # The following 0.5 is a kludge
                        F2[ibf+i,ibf+i] += 0.5*qj*gammaij
                        F2[jbf+j,jbf+j] += 0.5*qi*gammaij
            jbf += atomj.nbf
        ibf += atomi.nbf
    return F2

def get_F1(atoms,D):
    "One-center corrections to the core fock matrix"
    nbf = get_nbf(atoms)
    nat = len(atoms)
    
    F1 = zeros((nbf,nbf),'d')

    ibf = 0 # bf number of the first bfn on iat
    for iat in range(nat):
        atomi = atoms[iat]
        for i in range(atomi.nbf):
            bfi = atomi.basis[i]
            gii = get_g(bfi,bfi)
            qi =  D[ibf+i,ibf+i]
            F1[ibf+i,ibf+i] = 0.5*qi*gii
            
            for j in range(atomi.nbf):  # ij on same atom
                if j != i:
                    bfj = atomi.basis[j]
                    qj = D[ibf+j,ibf+j]
                    gij = get_g(bfi,bfj)
                    pij = D[ibf+i,ibf+j]
                    hij = get_h(bfi,bfj)
                    # the following 0.5 is something of a kludge to match
                    #  the mopac results.
                    F1[ibf+i,ibf+i] += qj*gij - 0.5*qj*hij
                    F1[ibf+i,ibf+j] += 0.5*pij*(3*hij-gij)
        ibf += atomi.nbf
    return F1

 def bfs(self,i):
     "Return a vector of the product of two basis functions "
     " over the entire grid"
     bfs = zeros(len(self.points),'d')
     for j in xrange(len(self.points)):
         bfs[j] = self.points[j].bfs[i]
     return bfs

 def get_gamma(self):
     "Return the density gradient gamma for each point in the grid"
     if not self.do_grad_dens: return None
     gamma = zeros(len(self.points),'d')
     for i in xrange(len(self.points)):
         gamma[i] = self.points[i].get_gamma()
     return gamma

 def grad(self):
     pts = self._points
     npts = len(pts)
     gr = zeros((npts,3),'d')
     for i in range(npts):
         gr[i,:] = pts[i].grad()
     return gr