Reading MOs from GAMESSUS FMO output

S , 2012-06-09 00:30 +0300, My Th rakstīja:

Hi!

I have a GAMESS output file which kills Avogadro when loading MOs:

push back 0.013014
add MOs: 50 1052
*** glibc detected *** /usr/bin/avogadro: malloc(): memory corruption:
0x0000000002bec940 ***
When calling new in OpenQube::GAMESSUSOutput::load at
libavogadro/src/extensions/surfaces/openqube/gamessus.cpp:257

So I’m reading trough the code trying to understand where is the issue.
The relevant files seem to be gamessus.{h,cpp} and gaussianset.{h,cpp}.
I haven’t gone trough all of it yet, but the first thing I noticed in
gaussianset.h is this:

  • The GaussianSet class has a transparent data structure for storing the basis
  • sets output by many quantum mechanical codes. It has a certain hierarchy
  • where shells are built up from n primitives, in this case Gaussian Type
  • Orbitals (GTOs). Each shell has a type (S, P, D, F, etc) and is composed of
  • one or more GTOs. Each GTO has a contraction coefficient, c, and an exponent,
  • a.

It says that there is one exponent and one contraction coefficient for
each GTO and the code seems to also treating it like that, but my output
file has two contraction coefficients for L type shells (it should be
6-31G** basis set):

SHELL TYPE PRIMITIVE EXPONENT CONTRACTION COEFFICIENT(S)

O

  1   S       1          5484.6716600    0.001831074430
  1   S       2           825.2349460    0.013950172200
  1   S       3           188.0469580    0.068445078098
  1   S       4            52.9645000    0.232714335992
  1   S       5            16.8975704    0.470192897984
  1   S       6             5.7996353    0.358520852987

  2   L       7            15.5396162   -0.110777549525    0.070874268231
  2   L       8             3.5999336   -0.148026262701    0.339752839147
  2   L       9             1.0137618    1.130767015354    0.727158577316

  3   L      10             0.2700058    1.000000000000    1.000000000000

  4   D      11             0.8000000    1.000000000000

Did I read it correctly and it is a bug in Avogadro or am I missing
something? There are other things which might trow off the parser (it is
output from FMO calculation), but I want to have this clarified first.

Ok, never mind. L type is translated to SP and that checks for and reads
the extra contraction coefficient and afterwards it is expanded into two
GTOs - S and P.

The real issue is that GAMESSUSOutput::processLine() is ignoring
elements when processing “ATOMIC BASIS SET” and relies on having a basis
printed for each atom in molecule and in the right order. Now with the
output of FMO (Fragment Molecular Orbital method) calculation these
assumptions don’t hold. In this case we have only one basis printed for
each atom type, and MOs are printed for fragments, not the whole
molecule.

I have a molecule consisting of 4 elements, that gives total of 50 AOs,
but the first fragment has 14 atoms and 205 AOs (for which 192 MOs are
written). And that confuses gaussianset.cpp::GaussianSet::addMOs(), it
tries to fit 205x192=39360 elements into 50x50=2500 matrix. It looks
like Eigen::MatrixXd doesn’t check indexes for coeffRef(i, j) and it
ends up with glibc panicking.

I guess it would be relatively easy to add a check and bail out if
inconsistency is detected, but adding proper FMO output support would be
more tricky. One way would be by adding multilevel mode support for
gamessus.cpp::GAMESSUSOutput::processLine() so that m_currentMode works
as now, by reading eigenvectors, but additionally would direct them to
the right fragment. Any ideas comments?

Reinis

S , 2012-06-09 11:45 +0300, My Th rakstīja:

S , 2012-06-09 00:30 +0300, My Th rakstīja:

Hi!

I have a GAMESS output file which kills Avogadro when loading MOs:

push back 0.013014
add MOs: 50 1052
*** glibc detected *** /usr/bin/avogadro: malloc(): memory corruption:
0x0000000002bec940 ***
When calling new in OpenQube::GAMESSUSOutput::load at
libavogadro/src/extensions/surfaces/openqube/gamessus.cpp:257

The real issue is that GAMESSUSOutput::processLine() is ignoring
elements when processing “ATOMIC BASIS SET” and relies on having a basis
printed for each atom in molecule and in the right order. Now with the
output of FMO (Fragment Molecular Orbital method) calculation these
assumptions don’t hold. In this case we have only one basis printed for
each atom type, and MOs are printed for fragments, not the whole
molecule.

I have a molecule consisting of 4 elements, that gives total of 50 AOs,
but the first fragment has 14 atoms and 205 AOs (for which 192 MOs are
written). And that confuses gaussianset.cpp::GaussianSet::addMOs(), it
tries to fit 205x192=39360 elements into 50x50=2500 matrix. It looks
like Eigen::MatrixXd doesn’t check indexes for coeffRef(i, j) and it
ends up with glibc panicking.

I guess it would be relatively easy to add a check and bail out if
inconsistency is detected, but adding proper FMO output support would be
more tricky. One way would be by adding multilevel mode support for
gamessus.cpp::GAMESSUSOutput::processLine() so that m_currentMode works
as now, by reading eigenvectors, but additionally would direct them to
the right fragment. Any ideas comments?

Attached is a patch which stops loading and displaying MOs if there are
more eigenvectors than expected from basis (if there is not given basis
for each atom in the system). With this it aborts MO loading and
displays just the structures from optimization run. I’m not sure if the
checks are put in quite the right places, someone more familiar with the
code could comment on that.

Reinis

S , 2012-06-09 15:26 +0300, My Th rakstīja:

Attached is a patch which stops loading and displaying MOs if there are
more eigenvectors than expected from basis (if there is not given basis
for each atom in the system). With this it aborts MO loading and
displays just the structures from optimization run. I’m not sure if the
checks are put in quite the right places, someone more familiar with the
code could comment on that.

Interestingly emerge happily compiled it for live version of avogadro
even tough MOs is passed as const reference and I was calling
MOs.clear(). Maybe because I’m emerging with -O3 … Building manually
it complained about it.

Reinis