Menu
▾
▴
cclib-devel
|
From: Matt W <pb...@gm...> - 2010-02-05 06:26:14
|
Hello all, I recently wrote a short bit of code for cclib to read the Hessian from a Gaussian output file. If there is interest in adding this to the project, what is the best way to submit it? Best Regards, Matt W |
|
From: Noel O'B. <bao...@gm...> - 2010-02-07 15:16:18
|
Hi Matt, Sounds good. If you could send a patch to this list, we can take a look. We already have some hessian code somewhere (I think?) so we will need to make sure that the units are standardised and so forth. If you checked our code out of subversion, "svn diff" will create the patch. Otherwise, use "diff -u oldfile newfile > a.patch". If all else fails, just attach the modified file. - Noel On 5 February 2010 06:26, Matt W <pb...@gm...> wrote: > Hello all, > I recently wrote a short bit of code for cclib to read the Hessian from a > Gaussian output file. If there is interest in adding this to the project, > what is the best way to submit it? > > Best Regards, > > Matt W > > > > ------------------------------------------------------------------------------ > The Planet: dedicated and managed hosting, cloud storage, colocation > Stay online with enterprise data centers and the best network in the > business > Choose flexible plans and management services without long-term contracts > Personal 24x7 support from experience hosting pros just a phone call away. > http://p.sf.net/sfu/theplanet-com > _______________________________________________ > cclib-devel mailing list > ccl...@li... > https://lists.sourceforge.net/lists/listinfo/cclib-devel > > |
|
From: Matt W <pb...@gm...> - 2010-02-10 03:28:30
|
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"])
|
|
From: Noel O'B. <bao...@gm...> - 2010-02-12 09:51:58
|
Looks really good. Having the tests is great. I particularly like the way you stuck a full Hessian to vibfreq conversion in there. We might move this out to the algorithms at some point to make it more widely available. I'll try to do the patch this weekend. Feel free to contribute incognito, but I'm wondering who exactly do we owe this code to. Google isn't exactly helpful in this case. :-) - Noel On 10 February 2010 03:28, Matt W <pb...@gm...> wrote: > Hi Noel, > Attached is a patch file against rev 881. I actually generated it with > mercurial, so the -p1 option is needed with patch. Molpro seems to use > atomic units (Hartrees/Bohr^2) for the Hessian, and the Molpro parser > doesn't do any conversion. Gaussian uses the same units, so I also didn't > do any conversion, although units of eV/Angstrom^2 would seem to be more > consistent with other cclib quantities. > > I noticed that the Molpro parser was reading in atomic masses but they > weren't being made public, so I changed that (in data.py) and wrote a test > (testhess.py) that compares the vib freqs calculated from the Hessian and > the masses to the freqs read in directly from the Gaussian or Molpro file. > One strange thing is that Molpro prints abundance averaged atomic weights by > default, while Gaussian prints relative atomic (isotopic) masses (carbon-12 > = 12 exactly). If you'd rather not have the additional changes, the part of > the patch against gaussianparser.py is fine on its own. > > Best Regards, > > Matt White > > On Sun, Feb 7, 2010 at 7:16 AM, Noel O'Boyle <bao...@gm...> wrote: >> >> Hi Matt, >> >> Sounds good. If you could send a patch to this list, we can take a >> look. We already have some hessian code somewhere (I think?) so we >> will need to make sure that the units are standardised and so forth. >> >> If you checked our code out of subversion, "svn diff" will create the >> patch. Otherwise, use "diff -u oldfile newfile > a.patch". If all else >> fails, just attach the modified file. >> >> - Noel >> >> On 5 February 2010 06:26, Matt W <pb...@gm...> wrote: >> > Hello all, >> > I recently wrote a short bit of code for cclib to read the Hessian >> > from a >> > Gaussian output file. If there is interest in adding this to the >> > project, >> > what is the best way to submit it? >> > >> > Best Regards, >> > >> > Matt W >> > >> > >> > >> > >> > ------------------------------------------------------------------------------ >> > The Planet: dedicated and managed hosting, cloud storage, colocation >> > Stay online with enterprise data centers and the best network in the >> > business >> > Choose flexible plans and management services without long-term >> > contracts >> > Personal 24x7 support from experience hosting pros just a phone call >> > away. >> > http://p.sf.net/sfu/theplanet-com >> > _______________________________________________ >> > cclib-devel mailing list >> > ccl...@li... >> > https://lists.sourceforge.net/lists/listinfo/cclib-devel >> > >> > > > > ------------------------------------------------------------------------------ > SOLARIS 10 is the OS for Data Centers - provides features such as DTrace, > Predictive Self Healing and Award Winning ZFS. Get Solaris 10 NOW > http://p.sf.net/sfu/solaris-dev2dev > _______________________________________________ > cclib-devel mailing list > ccl...@li... > https://lists.sourceforge.net/lists/listinfo/cclib-devel > > |
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X