Re: [cclib-devel] Hessian code for gaussianparser.py

[Matt W <pb...@gm...>]

diff -r 881155faef88 -r 7456e8b3840b src/cclib/parser/data.py
--- a/src/cclib/parser/data.py
+++ b/src/cclib/parser/data.py
@@ -69,7 +69,7 @@
         """
 
         # Names of all supported attributes.
-        self._attrlist = ['aonames', 'aooverlaps', 'atombasis',
+        self._attrlist = ['amass', 'aonames', 'aooverlaps', 'atombasis',
                           'atomcoords', 'atomnos',
                           'ccenergies', 'charge', 'coreelectrons',
                           'etenergies', 'etoscs', 'etrotats', 'etsecs', 'etsyms',
@@ -82,7 +82,8 @@
                           'vibdisps', 'vibfreqs', 'vibirs', 'vibramans', 'vibsyms']
 
         # The expected types for all supported attributes.
-        self._attrtypes = { "aonames":        list,
+        self._attrtypes = { "amass":          numpy.ndarray,
+                            "aonames":        list,
                             "aooverlaps":     numpy.ndarray,
                             "atombasis":      list,
                             "atomcoords":     numpy.ndarray,
diff -r 881155faef88 -r 7456e8b3840b src/cclib/parser/gaussianparser.py
--- a/src/cclib/parser/gaussianparser.py
+++ b/src/cclib/parser/gaussianparser.py
@@ -969,6 +969,58 @@
         if line[1:7] == "ONIOM:":
             self.oniom = True
 
+        # Read in entire archive section, then extract Hessian if present
+        # The lower-triangular part of the Hessian is returned flattened.
+        # To convert to a 3N x 3N rank 2 matrix:
+        #   hess = numpy.zeros([3N,3N]); hess[numpy.tril_indices(3N)] = data.hessian
+        # The archive section is the only place the Hessian is found unless
+        #  additional print ("p") was specified in the Gaussian input file
+        # Under Windows, Gaussian uses "|" instead of "\" as separator
+        #  in archive section
+        if line[1:5] == "1|1|" or line[1:5] == "1\\1\\":
+
+            self.updateprogress(inputfile, "Archive Section", self.fupdate)            
+            arcsep = line[4];
+            arcstr = ""
+            # archive section is terminated by "||@", but just look for 
+            #  blank line following
+            while line.strip():
+                arcstr += line.strip()
+                line = inputfile.next()
+            # sections separated by "||"
+            arcsects = arcstr.split(arcsep + arcsep)
+            
+            # Section immediately preceeding Hessian ends with "NImag=<int>"
+            # Note that output file may contain multiple archive sections (e.g.,
+            #  for an opt + freq job, Hessian is in 2nd one).
+            for ii in range(len(arcsects)):
+                if arcsects[ii].find("NImag=") > -1:
+                    self.hessian = map(float, arcsects[ii+1].split(","))
+                       
+        # Atomic masses: In "Isotopes and Nuclear Properties:" section if
+        #  additional print ("p") is on, otherwise only present if freq job
+        #  (in Thermochemistry section).
+        # The "Isotopes and Nuclear Properties:" section may appear more than
+        #  once (e.g opt + freq job); this code will yield values from the final one
+        if line.find("Isotopes and Nuclear Properties") > -1:
+            self.amass = []
+
+        if line[1:8] == "AtmWgt=":
+            self.amass += map(float, line.split()[1:])
+        
+        # Get mass from Thermochemistry section if we haven't found it yet
+        # looking for "Atom  1 has atomic number x and mass y"; space 
+        #  between "Atom" and "1" is different in G03 and G09
+        if line.find("1 has atomic number") > -1 and not hasattr(self, "amass"):
+            self.amass = []
+            broken = line.split()
+            # no blank line following, so use len(split()) to detect end
+            while len(broken) == 9:
+                self.amass.append(float(broken[8]))
+                line = inputfile.next()
+                broken = line.split()
+
+
 if __name__ == "__main__":
     import doctest, gaussianparser
     doctest.testmod(gaussianparser, verbose=False)
diff -r 881155faef88 -r 7456e8b3840b test/testdata
--- a/test/testdata
+++ b/test/testdata
@@ -100,3 +100,7 @@
 vib       Molpro      MolproIRTest          basicMolpro2006       dvb_ir.out            dvb_ir.log
 vib       ORCA        OrcaIRTest            basicORCA2.6          dvb_ir.out
 vib       ORCA        OrcaRamanTest         basicORCA2.6          dvb_raman.out
+
+hess      Gaussian    Gaussian03HessTest    basicGaussian03       dvb_ir.out
+hess      Gaussian    Gaussian09HessTest    basicGaussian09       dvb_raman.log
+hess      Molpro      MolproHessTest        basicMolpro2006       dvb_ir.out            dvb_ir.log
diff -r 881155faef88 -r 7456e8b3840b test/testhess.py
--- /dev/null
+++ b/test/testhess.py
@@ -0,0 +1,67 @@
+__revision__ = "$Revision$"
+
+import numpy
+import bettertest
+
+
+class GenericHessTest(bettertest.TestCase):
+    """Generic Hessian unittest."""
+
+    def modes_from_hess(self, flathess, mass, units=5140.4872066):
+        """ Calculate normal mode freqs (in 1/cm) and vectors from hessian.
+        Default units=5140.4872066 assumes Hessian in Hartree/Bohr^2 and mass 
+        in amu (12C := 12 amu) """
+  
+        # get full Hessian as square matrix
+        N = 3*len(mass);  # number of modes
+        hess = numpy.zeros((N,N))
+        # trilidx = numpy.where(numpy.tril(numpy.ones((hdim,hdim),int),0) != 0)  # NumPy 1.3
+        hess[numpy.tril_indices(N)] = flathess
+        # fill in upper triangle
+        hess = hess.transpose()
+        hess[numpy.tril_indices(N)] = flathess
+      
+        # create mass weighting matrix - repeat for x,y,z 
+        invrootmass = numpy.repeat(1/numpy.sqrt(numpy.array(mass)), 3)
+        # for weighing hessian
+        ivmassouter = numpy.outer(invrootmass, invrootmass)
+        # we need a column vector for scaling normal modes all at once
+        invrootmass = numpy.reshape(invrootmass, (-1,1))
+  
+        # Perform mass-weighting and diagonalize to get normal modes
+        freq, modes = numpy.linalg.eigh(hess * ivmassouter)
+        # convert freqs to 1/cm
+        freq = numpy.sign(freq) * numpy.sqrt(numpy.abs(freq)) * units
+        
+        # scale and renormalize eigenvectors
+        modes = modes * invrootmass
+        modes = numpy.transpose(modes/numpy.sqrt(numpy.sum(modes**2, 0)))
+        # after reshape, indicies are (<mode number>, <atom number>, X, Y, or Z) 
+        modes = numpy.reshape(modes, (N, N/3, 3))
+        # The LAPACK routine used by numpy's eigh() always returns eigenvalues in 
+        #  ascending order, so we can safely assume the first 6 modes are the zero
+        #  freq modes. This only works for a PES minimum, not a transition state
+        return freq[6:], modes[6:]
+
+    def testfreqval(self):
+        """ Do the vibfreqs read directly match those calculated from the Hessian within +/-0.1/cm? """
+        freqs, modes = self.modes_from_hess(self.data.hessian, self.data.amass)
+        self.assertEqual(len(self.data.vibfreqs), len(freqs))
+        for delta in (self.data.vibfreqs - freqs):
+            self.assertInside(delta, 0, 0.1)
+      
+
+class Gaussian03HessTest(GenericHessTest):
+    """Gaussian 03 Hessian unittest."""
+
+class Gaussian09HessTest(GenericHessTest):
+    """Gaussian 09 Hessian unittest."""
+
+class MolproHessTest(GenericHessTest):
+    """Molpro Hessian unittest."""
+
+
+if __name__=="__main__":
+
+    from testall import testall
+    testall(modules=["hess"])