Added Hatem's Analytic Force Module

git-svn-id: https://pyquante.svn.sourceforge.net/svnroot/pyquante/trunk@89 64417113-1622-0410-aef8-ef15d1a3721e

PyQuante/Atom.py

@@ -23,13 +23,17 @@
 #  whatever you store you get back.
 
 class Atom:
-    def __init__(self,atno,x,y,z):
+    def __init__(self,atno,x,y,z,atid=0,fx=0.0,fy=0.0,fz=0.0):
         self.atno = atno
         self.r = array([x,y,z],'d')
+        #added by Hatem Helal hhh23@cam.ac.uk
+        #atom id defaults to zero so as not to break preexisting code...
+        self.atid = atid
+        self.forces = array([fx,fy,fz],'d')
         return
 
-    def __repr__(self): return "Atom %2d (%6.3f,%6.3f,%6.3f)" % \
-        (self.atno,self.r[0],self.r[1],self.r[2])
+    def __repr__(self): return "Atom ID: %d Atomic Num: %2d (%6.3f,%6.3f,%6.3f)" % \
+        (self.atid,self.atno,self.r[0],self.r[1],self.r[2])
     def __getitem__(self, i):
         return self.r[i]
     def mass(self): return mass[self.atno]
@@ -52,6 +56,8 @@ def get_nel_mindo(self): return self.Z
 
     def update_coords(self,xyz): self.r = array(xyz)
     def update_from_atuple(self,(atno,xyz)): self.update_coords(xyz)
+
+    def set_force(self,fxfyfz): self.forces = array(fxfyfz)
 
 def test():
     at1 = Atom(1,0,0,0)

PyQuante/CGBF.py

@@ -25,19 +25,23 @@
 
 class CGBF:
     "Class for a contracted Gaussian basis function"
-    def __init__(self,origin,powers=(0,0,0)):
+    def __init__(self,origin,powers=(0,0,0),atid=0):
         self._origin = tuple([float(i) for i in origin])
         self._powers = powers
         self._normalization = 1.
         self._prims = []
         self._pnorms = []
         self._pexps = []
         self._pcoefs = []
+        #added by Hatem H Helal hhh23@cam.ac.uk
+        #stores atom id number for easy method to identify which atom
+        #a particular bf is centered on
+        self.atid = atid
         return
 
     def __repr__(self):
-        s = "<cgbf origin=\"(%f,%f,%f)\" powers=\"(%d,%d,%d)\">\n" % \
-            (self._origin[0],self._origin[1],self._origin[2],
+        s = "<cgbf atomid=%d origin=\"(%f,%f,%f)\" powers=\"(%d,%d,%d)\">\n" % \
+            (self.atid,self._origin[0],self._origin[1],self._origin[2],
              self._powers[0],self._powers[1],self._powers[2])
         for prim in self._prims:
             s = s + prim.prim_str(self.norm())

PyQuante/LA2.py

@@ -144,3 +144,17 @@ def mkdens_spinavg(c,nclosed,nopen):
        shell orbitals and *nopen* open shell orbitals"""
     return mkdens(c,0,nclosed) + 0.5*mkdens(c,nclosed,nclosed+nopen)
 
+#added by Hatem H Helal 18.07.2007
+#RPM: Consider using diag() for this
+def diagonal_mat(vector):
+    #takes a vector with N components and returns an NXN matrix whose diagonal elements
+    #are the components of the vector
+    len = vector.shape[0]
+
+    matrix = zeros((len,len),'d')
+
+    i=0
+    for i in range(len):
+        matrix[i][i]=vector[i]
+
+    return matrix    

PyQuante/Molecule.py

@@ -73,15 +73,16 @@ def translate(self,pos):
 
     def add_atom(self,atom): self.atoms.append(atom)
 
-    def add_atuple(self,atno,xyz):
+    def add_atuple(self,atno,xyz,atid):
         if self.units != 'bohr': xyz = toBohr(xyz[0],xyz[1],xyz[2])
         if type(atno) == type(''): atno = sym2no[atno]
-        self.atoms.append(Atom(atno,xyz[0],xyz[1],xyz[2]))
+        self.atoms.append(Atom(atno,xyz[0],xyz[1],xyz[2],atid))
 
     def add_atuples(self,atoms):
         "Add a list of (atno,(x,y,z)) tuples to the atom list"
         from Atom import Atom
-        for atno,xyz in atoms: self.add_atuple(atno,xyz)
+        id=0
+        for atno,xyz in atoms: self.add_atuple(atno,xyz,id); id+=1
         return
 
     def atuples(self):
@@ -140,7 +141,7 @@ def get_closedopen(self):
             print "Incompatible number of electrons and spin multiplicity"
             print "nel = ",nel
             print "multiplicity = ",multiplicity
-            raise Exception("Incompatible number of electrons and spin multiplicity")
+            raise Exception("Incompatible # electrons and multiplicity")
 
         nopen = multiplicity-1
         nclosed,ierr = divmod(nel-nopen,2)
@@ -191,7 +192,7 @@ def ParseXYZLines(name,xyz_lines,**opts):
         sym = cleansym(words[0])
         xyz = map(float,words[1:4])
         atoms.append((sym,xyz))
-    return Molecule(name,atoms,units="Angstrom",**opts)
+    return Molecule(name,atoms,**opts)
 
 def ParseXYZFile(fname,**opts):
     from PyQuante.Util import parseline

PyQuante/NumWrap.py

@@ -22,6 +22,9 @@
 # As of 2/12/2007 we have updated PyQuante to the numpy convention,
 # where eigenvectors are kept in columns.
 if use_numpy:
+    from numpy import array2string #added by Hatem Helal hhh23@cam.ac.uk
+    from numpy import abs #added by Hatem Helal hhh23@cam.ac.uk
+
     from numpy import array,zeros,concatenate,dot,ravel,arange
     from numpy import arcsinh,diagonal,identity,choose,transpose
     from numpy import reshape,take,trace
@@ -34,7 +37,7 @@
     from numpy.linalg import solve
     from numpy.linalg import inv
     from numpy.linalg import cholesky
-
+    
     # still need to kill these two, which are used by Optimize:
     import numpy.oldnumeric.mlab as MLab
     from numpy.oldnumeric import NewAxis

PyQuante/Wavefunction.py

@@ -0,0 +1,170 @@
+#!/usr/bin/env python
+"""\
+ Wavefunction.py: trial module defining an orbital class
+ The idea is to make handling computations over orbitals slightly easier
+
+ Started by Hatem H Helal hhh23@cam.ac.uk on 17/07/2007
+
+ This program is part of the PyQuante quantum chemistry program suite.
+
+ Copyright (c) 2004, Richard P. Muller. All Rights Reserved. 
+
+ PyQuante version 1.2 and later is covered by the modified BSD
+ license. Please see the file LICENSE that is part of this
+ distribution. 
+"""
+
+from NumWrap import matrixmultiply,transpose
+from LA2 import diagonal_mat
+class Wavefunction:
+    """\
+    Data class to hold all relevant wavefunction data in PyQuante
+
+    Example of how to use this in a restricted HF (or DFT or ....) calculation:
+    Wavefunction(orbs=orbital_coef_matrix,orbe=orbital_eigenvalues, \
+                nclosed=closed_shell_orbitals, nopen=open_shell_orbitals) 
+
+    The unrestricted analog:
+    Wavefunction(orbs_a=spin_up_coef_matrix, orbs_b=spin_down_coef_matrix, \
+            orbe_a=spin_up_orbital_eigenvalues, orbe_b=spin_down_orbital_eigenvalues, \
+            nalpha=spin_up_electrons, nbeta=spin_down_electrons)
+    
+    This class also allows for a container for orbital occupation numbers which
+    could be useful for caclulations involving non-standard occupations, excitations,
+    or thermal smearing through a Fermi-Dirac function.
+
+    Options:       Description
+    --------       -----------
+Restricted
+    orbs            Orbital coefficient matrix, columns are eigenstates of 
+                    the Fock matrix arranged in ascending order by eigenvalue
+
+    orbe            Orbital eigenvalues in list form
+    nclosed         Number of closed shell orbitals
+    nopen           Number of open shell orbitals
+    occs            Orbital occupation numbers, useful for non-standard occupations
+
+Unrestricted 
+    orbs_a          Spin up orbital coefficient matrix, columns are the eigenstates
+                    of the spin up Fock matrix arranged in ascending order by eigenvalue
+    orbs_b          Spin down analog of orbs_a
+    
+    orbe_a          Spin up orbital eigenvalues in list form
+    orbe_b          Spin down analog of orbe_a
+
+    nalpha          Number of spin up electrons
+    nbeta           Number of spin down electrons
+
+    occs_a          Spin up occupation numbers
+    occs_b          Spin down analog of occs_a
+    """
+    def __init__(self,**opts):
+        #orbital coefs, energies, and occupations for restricted calculations
+        self.orbs    = opts.get('orbs',None)
+        self.orbe    = opts.get('orbe',None)
+        self.nclosed = opts.get('nclosed',None)
+        self.nopen   = opts.get('nopen',None)
+        self.occs    = opts.get('occs',None)
+
+        #orbital coefs, energies, and occupations for unrestricted calculations        
+        self.orbs_a = opts.get('orbs_a',None)
+        self.orbs_b = opts.get('orbs_b',None)
+
+        self.orbe_a = opts.get('orbe_a',None)
+        self.orbe_b = opts.get('orbe_b',None)
+
+        nalpha      = opts.get('nalpha',None)
+        nbeta       = opts.get('nbeta',None)
+
+        self.occs_a = opts.get('occs_a',None)
+        self.occs_b = opts.get('occs_b',None)
+
+        self.restricted   = opts.get('restricted',None)
+        self.unrestricted = opts.get('unrestricted',None)
+
+        return
+
+    def update_wf(self,**opts):
+        #use this function to update coef and eigenvalues within SCF loop
+        #restricted
+        self.orbs    = opts.get(orbs,None)
+        self.orbe    = opts.get(orbe,None)
+
+        #unrestricted
+        self.orbs_a = opts.get(orbs_a,None)
+        self.orbs_b = opts.get(orbs_b,None)
+
+        self.orbe_a = opts.get(orbe_a,None)
+        self.orbe_b = opts.get(orbe_b,None)
+
+        return
+
+    def set_occs(self,**opts):
+        #sets orbital occupation arrays
+        self.occs   = opts.get(occs,None)
+        self.occs_a = opts.get(occs_a,None)
+        self.occs_b = opts.get(occs_b,None)
+
+        return
+
+
+    def mkdens(self):
+        #form density matrix D for restricted wavefuntions
+        #for unrestricted returns Da and Db matrices
+        if self.restricted:
+            #makes spin avg density matrix for open shell systems
+            D = self.mk_R_Dmat(0,self.nclosed) #
+            if self.nopen>=1:
+                D += 0.5*self.mk_R_Dmat(0,self.nclosed+self.nopen)
+            return D
+        if self.unrestricted:
+            Da = self.mk_A_Dmat(0,self.nalpha)
+            Db = self.mk_B_Dmat(0,self.nbeta)
+            return Da,Db
+
+    def mk_R_Dmat(self,nstart,nstop):
+        "Form a density matrix C*Ct given eigenvectors C[nstart:nstop,:]"
+        d = self.orbs[:,nstart:nstop]
+        Dmat = matrixmultiply(d,transpose(d))
+        return Dmat
+
+    def mk_A_Dmat(self,nstart,nstop):
+        "Form a density matrix C*Ct given eigenvectors C[nstart:nstop,:]"
+        d = self.orbs_a[:,nstart:nstop]
+        Dmat = matrixmultiply(d,transpose(d))
+        return Dmat
+
+    def mk_B_Dmat(self,nstart,nstop):
+        "Form a density matrix C*Ct given eigenvectors C[nstart:nstop,:]"
+        d = self.orbs_b[:,nstart:nstop]
+        Dmat = matrixmultiply(d,transpose(d))
+        return Dmat
+
+    def mkQmatrix(self):
+        #computes Q matrix which is basically density matrix weighted by the orbital eigenvalues
+        if self.restricted:
+            occ_orbe = self.orbe[0:self.nclosed]
+            occ_orbs = self.orbs[:,0:self.nclosed]
+            Qmat = matrixmultiply( occ_orbs, matrixmultiply( diagonal_mat(occ_orbe), transpose(occ_orbs) ) )
+            return Qmat
+
+        if self.unrestricted:
+            occ_orbs_A = self.orbs_a[:,0:self.nalpha]
+            occ_orbs_B = self.orbs_b[:,0:self.nbeta]
+            occ_orbe_A = self.orbe_a[0:self.nalpha]
+            occ_orbe_B = self.orbe_b[0:self.nbeta]
+
+            Qa = matrixmultiply( occ_orbs_A, matrixmultiply( diagonal_mat(occ_orbe_A), transpose(occ_orbs_A) ) )
+            Qb = matrixmultiply( occ_orbs_B, matrixmultiply( diagonal_mat(occ_orbe_B), transpose(occ_orbs_B) ) )
+            return Qa,Qb
+
+    def mk_auger_dens(self):
+        "Forms a density matrix from a coef matrix c and occupations in occ"
+        #count how many states we were given
+        nstates = self.occc.shape[0]
+        D = 0.0
+        for i in range(nstates):
+            D += self.occs[i]*dot( self.orbs[:,i:i+1], transpose(self.orbs[:,i:i+1]))
+        #pad_out(D)
+        return D
+