def __init__(self,atno,x,y,z,atid=0,fx=0.0,fy=0.0,fz=0.0,vx=0.0,vy=0.0,vz=0.0):
     self.atno = atno
     self.r = array([x,y,z],'d')
     #added by Hatem Helal [email protected]
     #atom id defaults to zero so as not to break preexisting code...
     self.atid = atid
     self.f = array([fx,fy,fz],'d')
     self.vel = array([vx,vy,vz],'d')
     return

 def gamma(self):
     "Return the density gradient gamma for each point in the grid"
     if not self.do_grad_dens: return None
     gs = []
     for agr in self.atomgrids:
         gs.extend(agr.gamma())
     return array(gs)

 def make_bfgrid(self):
     "Construct a matrix with bfs in columns over the entire grid, "
     " so that R[0] is the first basis function, R[1] is the second..."
     bfs = []
     for point in self.points:
         bfs.extend(point.bfs)
     return array(bfs)

 def grad_old(self,pos):
     amp = self.amp(pos[0],pos[1],pos[2])
     alpha = self.exp
     L,M,N = self.powers
     x = pos[0]-self.origin[0]
     y = pos[1]-self.origin[1]
     z = pos[2]-self.origin[2]
     val = array([L/x - 2*x*alpha,M/y - 2*y*alpha,N/z-2*z*alpha])
     return self.norm*self.coef*val*amp

 def overlap_1(self,other):
     "Overlap derivative integral between two gaussians. THO eq. 2.14."
     return self.norm*other.norm*\
      array([overlap_1x(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin),
             overlap_1y(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin),
             overlap_1z(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin)],'d')

 def kinetic_1(self,other):
     "Kinetic derivative integral between two gaussians. THO eq. 2.14."
     return self.norm*other.norm*\
      array([kinetic_1x(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin),
             kinetic_1y(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin),
             kinetic_1z(self.exp,self.powers,self.origin,
                        other.exp,other.powers,other.origin)],'d')

 def grad_old(self,pos):
     amp = self.amp(pos[0],pos[1],pos[2])
     alpha = self._exponent
     L,M,N = self._powers
     x = pos[0]-self._origin[0]
     y = pos[1]-self._origin[1]
     z = pos[2]-self._origin[2]
     val = array([L/x - 2*x*alpha,M/y - 2*y*alpha,N/z-2*z*alpha])
     return self._normalization*self._coefficient*val*amp

def der_oneeE(a,D,bset,atoms):

    dH_dXa,dH_dYa,dH_dZa = der_Hcore_matrix(a,bset,atoms)
    
    doneE_Xa = trace2(D,dH_dXa)
    doneE_Ya = trace2(D,dH_dYa)
    doneE_Za = trace2(D,dH_dZa)
    
    #print doneE_Xa,doneE_Ya,doneE_Za
    return array([doneE_Xa,doneE_Ya,doneE_Za],'d')

def der_twoeE(a,D,bset):
    d2Ints_dXa,d2Ints_dYa,d2Ints_dZa  = der2Ints(a,bset)

    Gx,Gy,Gz = der2JmK(D,d2Ints_dXa,d2Ints_dYa,d2Ints_dZa)
    
    dtwoeE_Xa = trace2(D,Gx)
    dtwoeE_Ya = trace2(D,Gy)
    dtwoeE_Za = trace2(D,Gz)

    #print dtwoeE_Xa,dtwoeE_Ya,dtwoeE_Za
    return array([dtwoeE_Xa,dtwoeE_Ya,dtwoeE_Za],'d')

def der2Ints(a,bset):
    #modified from Ints.py -> get2ints
    """Store integrals in a long array in the form (ij|kl) (chemists
    notation. We only need i>=j, k>=l, and ij <= kl"""
    from array import array
    nbf = len(bset)
    totlen = nbf*(nbf+1)*(nbf*nbf+nbf+2)/8
    d2Ints_dXa = array('d',[0]*totlen)
    d2Ints_dYa = array('d',[0]*totlen)
    d2Ints_dZa = array('d',[0]*totlen)
    for i in xrange(nbf):
        for j in xrange(i+1):
            ij = i*(i+1)/2+j
            for k in xrange(nbf):
                for l in xrange(k+1):
                    kl = k*(k+1)/2+l
                    if ij >= kl:
                        ijkl = ijkl2intindex(i,j,k,l)
                        d2Ints_dXa[ijkl],d2Ints_dYa[ijkl],d2Ints_dZa[ijkl] =\
                                 der_Jints(a,bset[i],bset[j],bset[k],bset[l])
    return d2Ints_dXa,d2Ints_dYa,d2Ints_dZa

def getTA1B(bfsA,bfsB):
    """get the derivatives of T_ij between two molecules A and B.
The derivative refers to molecule A basis functions"""
    nbfA = len(bfsA)
    nbfB = len(bfsB)
    T = zeros((nbfA,nbfB,3),'d')
    
    for i in xrange(nbfA):
        bfi = bfsA[i]
        for j in xrange(nbfB):
            bfj = bfsB[j]
            T[i,j,:] = array( bfi.kinetic_1(bfj) , 'd' )
    return T

 def __init__(self,x,y,z,w=1.0,**opts):
     self.do_grad_dens = opts.get('do_grad_dens',False)
     self._x = float(x)
     self._y = float(y)
     self._z = float(z)
     self._w = float(w)
     self.xyz = array((self._x,self._y,self._z),'d')
     self._r = sqrt(self._x*self._x+self._y*self._y+self._z*self._z)
     self._gamma = None
     self._dens = 0
     self._dens0 = None
     self._grad = None
     self.bfs = []
     return

 def allbfs(self):
     "Construct a matrix with bfs in columns over the entire grid, "
     " so that R[0] is the first basis function, R[1] is the second..."
     bfs = []
     for agr in self.atomgrids:
         bfs.extend(agr.allbfs())
     bfs = array(bfs)
     npts = self.npts()
     nbf,nrem = divmod(len(bfs),npts)
     if nrem != 0: raise Exception("Remainder in divmod allbfs")
     nbf2 = self.nbf()
     if nbf != nbf2: raise Exception("Wrong # bfns %d %d" % (nbf,nbf2))
     bfs = reshape(bfs,(npts,nbf))
     return bfs

 def __init__(self,x,y,z,w=1.0,**kwargs):
     self.do_grad_dens = kwargs.get('do_grad_dens',settings.DFTDensityGradient)
     self._x = float(x)
     self._y = float(y)
     self._z = float(z)
     self._w = float(w)
     self.xyz = array((self._x,self._y,self._z),'d')
     self._r = sqrt(self._x*self._x+self._y*self._y+self._z*self._z)
     self._gamma = None
     self._dens = 0
     self._dens0 = None
     self._grad = None
     self.bfs = []
     return

 def grad(self,x,y,z):
     alpha = self.exp
     I,J,K = self.powers
     C = self.norm*self.coef
     x0,y0,z0 = self.origin
     fx = pow(x-x0,I)*exp(-alpha*pow(x-x0,2))
     fy = pow(y-y0,J)*exp(-alpha*pow(y-y0,2))
     fz = pow(z-z0,K)*exp(-alpha*pow(z-z0,2))
     gx = -2*alpha*(x-x0)*fx
     gy = -2*alpha*(y-y0)*fy
     gz = -2*alpha*(z-z0)*fz
     if I > 0: gx += pow(x-x0,I-1)*exp(-alpha*pow(x-x0,2))
     if J > 0: gy += pow(y-y0,J-1)*exp(-alpha*pow(y-y0,2))
     if K > 0: gz += pow(z-z0,K-1)*exp(-alpha*pow(z-z0,2))
     return array([C*gx*fy*fz,C*fx*gy*fz,C*fx*fy*gz])

 def grad(self,x,y,z):
     alpha = self._exponent
     I,J,K = self._powers
     C = self._normalization*self._coefficient
     x0,y0,z0 = self._origin
     fx = pow(x-x0,I)*exp(-alpha*pow(x-x0,2))
     fy = pow(y-y0,J)*exp(-alpha*pow(y-y0,2))
     fz = pow(z-z0,K)*exp(-alpha*pow(z-z0,2))
     gx = -2*alpha*(x-x0)*fx
     gy = -2*alpha*(y-y0)*fy
     gz = -2*alpha*(z-z0)*fz
     if I > 0: gx += pow(x-x0,I-1)*exp(-alpha*pow(x-x0,2))
     if J > 0: gy += pow(y-y0,J-1)*exp(-alpha*pow(y-y0,2))
     if K > 0: gz += pow(z-z0,K-1)*exp(-alpha*pow(z-z0,2))
     return array([C*gx*fy*fz,C*fx*gy*fz,C*fx*fy*gz])

def der_twoeE_uhf(a,Da,Db,bset):
    d2Ints_dXa,d2Ints_dYa,d2Ints_dZa  = der2Ints(a,bset)

    dJax,dJay,dJaz = derJ(Da,d2Ints_dXa,d2Ints_dYa,d2Ints_dZa)
    dJbx,dJby,dJbz = derJ(Db,d2Ints_dXa,d2Ints_dYa,d2Ints_dZa)
    
    dKax,dKay,dKaz = derK(Da,d2Ints_dXa,d2Ints_dYa,d2Ints_dZa)
    dKbx,dKby,dKbz = derK(Db,d2Ints_dXa,d2Ints_dYa,d2Ints_dZa)

    Dab = Da + Db
    
    dtwoeE_Xa = 0.5*trace2(Dab,dJax+dJbx) - 0.5*(trace2(Da,dKax)+trace2(Db,dKbx))
    dtwoeE_Ya = 0.5*trace2(Dab,dJay+dJby) - 0.5*(trace2(Da,dKay)+trace2(Db,dKby))
    dtwoeE_Za = 0.5*trace2(Dab,dJaz+dJbz) - 0.5*(trace2(Da,dKaz)+trace2(Db,dKbz))
    
    return array([dtwoeE_Xa,dtwoeE_Ya,dtwoeE_Za],'d')

def set_boltzmann_velocities(atoms,T):
    from random import gauss,randint
    Eavg = Rgas*T/2000 # kT/2 per degree of freedom (kJ/mol)

    vels = []
    for atom in atoms:
        m = atom.mass()
        vavg = sqrt(2*Eavg*fconst/m)
        
        stdev = 0.01 #I'm setting the std dev wrong here
        atom.v = array((pow(-1,randint(0,1))*gauss(vavg,stdev),
                        pow(-1,randint(0,1))*gauss(vavg,stdev),
                        pow(-1,randint(0,1))*gauss(vavg,stdev)))

    subtract_com_velocity(atoms)
    rescale_velocities(atoms,T)
    return

def der_dmat(a,Qmat,bset):
    """
    Looking at Szabo's equation C.12 you can see that this term results 
    from adding up the terms that involve derivatives of the density matrix
    elements P_mu,nu.  Following Szabo we compute this term as
      sum_mu,nu Q_mu,nu dS_mu,nu / dRa
    where the Q matrix is essentially a density matrix that is weighted by 
    the orbital eigenvalues.  Computing the derivative of the overlap integrals
    is done in the der_overlap_matrix function later in this module.
    """
    dS_dXa,dS_dYa,dS_dZa = der_overlap_matrix(a,bset)
    
    dDmat_Xa = trace2(Qmat,dS_dXa)
    dDmat_Ya = trace2(Qmat,dS_dYa)
    dDmat_Za = trace2(Qmat,dS_dZa)
    
    #print dDmat_Xa,dDmat_Ya,dDmat_Za
    return array([dDmat_Xa,dDmat_Ya,dDmat_Za],'d')

def chi_rk4(tend):
    # Same thing as chi_integrate, but uses a fourth-order Runge-Kutta
    #  scheme
    t = 0
    X = 1.0
    h = 1e-4 # 1e-5 is not small enough and still takes forever!
    v = -1.58807
    ts = []
    Xs = []
    while t < 1:
        ts.append(t)
        Xs.append(X)
        A = array((X,v))
        k1 = h*F(A,t)
        k2 = h*F(A+k1/2,t+h/2)
        k3 = h*F(A+k2/2,t+h/2)
        k4 = h*F(A+k3,t+h)
        A += (k1+2*k2+2*k3+k4)/6
        X,v = A
        t += h
        if X < 0: break
    return ts,Xs