Error in NICS computation - hnd giaox: geom tag to element failed
Just Got Here
7:12:59 PM PDT - Mon, Jun 23rd 2014
Hi all,

I am trying to compute NICS (nucleus-independent chemical shifts). I am encountering an error when trying to run a test included with NWChem, benzene_nics.nw:

  echo

  start benzene_nics

  title "benzene nics"
 
  geometry noautoz units angstrom nocenter
    symmetry c1
    C   0.98348719 0.98348719 0.00000000
    C  -1.34346849 0.35998130 0.00000000
    C   0.35998130 -1.34346849 0.00000000
    C  -0.98348719 -0.98348719 0.00000000
    C   1.34346849 -0.35998130 0.00000000
    C  -0.35998130 1.34346849 0.00000000
    H   1.74853940 1.74853940 0.00000000
    H  -2.38854924 0.64000984 0.00000000
    H   0.64000984 -2.38854924 0.00000000
    H  -1.74853940 -1.74853940 0.00000000
    H   2.38854924 -0.64000984 0.00000000
    H  -0.64000984 2.38854924 0.00000000
    bqH  0.00000000  0.0000000 0.00000000 
  end

  basis "ao basis" spherical
   c library 6-31g*
   h library 6-31g*
   bqH library H 6-31g*
  end

  dft
   xc hfexch
  end

  property
   shielding 1 13
  end
  task dft property


The error is get is the following:


          -----------------------------------------
          Chemical Shielding Tensors (GIAO, in ppm)
          -----------------------------------------

                                NWChem CPHF Module
                                ------------------

<edited to avoid error on forum software>

 Calc. par tensor-> nonrel
 ------------------------------------------------------------------------
 hnd_giaox: geom_tag_to_element failed                   0
 ------------------------------------------------------------------------
 ------------------------------------------------------------------------
  current input line : 
    37: task dft property
 ------------------------------------------------------------------------
 ------------------------------------------------------------------------
 There is an error related to the specified geometry
 ------------------------------------------------------------------------
                                                                                                                                                                                                                                                           
0:hnd_giaox: geom_tag_to_element failed:Received an Error in Communication
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI COMMUNICATOR 3 DUP FROM 0 
with errorcode -1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------


I saw the same error reported in a previous forum post (http://nwchemgit.github.io/Special_AWCforum/st/id594/NMR_NICS.html), which pointed to a solution by using a development version from August 2012. I used the latest development version, dated June 09, 2014.

I compiled NWChem on Linux (Ubuntu 14.04) with gfortran/gcc 4.8.2.

I will appreciate any pointers for help!

Regards,
Pablo

Gets Around
2:39:30 AM PDT - Tue, Jun 24th 2014
I don't have a solution but I have done a little digging. The code involved is under src/property. There appears to have been some copy-paste programming at work because hnd_giaox.F, hnd_giaox_zora.F, and hnd_spinspin.F all contain the line

errquit('hnd_giaox: geom_tag_to_element failed',0,GEOM_ERR)


After changing the error messages to be non-identical and recompiling, I see the problem with the benzene_nics test calculation actually appears to come from hnd_giaox_zora.F. This is a bit weird as I don't know why a non-relativistic calculation would end up in a file for ZORA.

Forum Regular
7:19:18 PM PDT - Tue, Jun 24th 2014
Indeed the error occurs in hnd_giaox_zora. The problem is that the code checks whether a tag for a center in the geometry corresponds to an element. However, it is valid for a center to be a Bq center (which is not an atom), but the code did not check whether the center is possibly a Bq center and chose to fail instead. Essentially the test in question in hnd_giaox_zora.F should match the one in hnd_giaox.F, including the test for a Bq center. I.e. the code in hnd_giaox_zora.F that reads 

 #include "case.fh"

 #include "zora.fh"



          if (.not. geom_tag_to_element(tag, symbol, element, atn)) call

      &       errquit('hnd_giaox: geom_tag_to_element failed',0,GEOM_ERR)

should read

 #include "case.fh"

 #include "inp.fh"

 #include "zora.fh"



          if (.not. geom_tag_to_element(tag, symbol, element, atn)) then

            if (.not. inp_compare(0,tag(1:2),'bq')) call ! check for bq as a fall back

      &       errquit('hnd_giaox_zora: geom_tag_to_element failed',

      &               0,GEOM_ERR)

Please let me know if you need help implementing that fix.

Best wishes, Huub

Gets Around
9:19:52 PM PDT - Tue, Jun 24th 2014
The calculation completes with the new code. The new code produces the same total shielding tensor as in the old code that generated the QA output, but the paramagnetic and diamagnetic tensors are different. I don't know if the change represents an improvement or a regression.

Old output, from QA/tests/benzene_nics/benzene_nics.out:
      Atom:    1  bq
        Diamagnetic
     26.0935      0.0000      0.0000
      0.0000     26.0935      0.0000
      0.0000      0.0000     39.7821

        Paramagnetic
    -17.7957      0.0000      0.0000
      0.0000    -17.7957      0.0000
      0.0000      0.0000    -22.4191

        Total Shielding Tensor
      8.2978      0.0000      0.0000
      0.0000      8.2978      0.0000
      0.0000      0.0000     17.3630

           isotropic =      11.3195
          anisotropy =       9.0652

          Principal Components and Axis System
                 1           2           3
               17.3630      8.2978      8.2978

      1         0.0000      0.7071      0.7071
      2         0.0000      0.7071     -0.7071
      3         1.0000      0.0000      0.0000


New output, from NWChem compiled today with changes by Huub:
      Atom:    1  bq
        Diamagnetic
     26.0935      0.0000     -0.0000
      0.0000     26.0935      0.0000
     -0.0000      0.0000    -29.8045

        Paramagnetic
    -17.7957      0.0000     -0.0000
      0.0000    -17.7957     -0.0000
     -0.0000     -0.0000     47.1675

        Total Shielding Tensor
      8.2978      0.0000     -0.0000
      0.0000      8.2978     -0.0000
     -0.0000     -0.0000     17.3630

           isotropic =      11.3195
          anisotropy =       9.0652

          Principal Components and Axis System
                 1           2           3
               17.3630      8.2978      8.2978

      1         0.0000      0.7071      0.7071
      2         0.0000      0.7071     -0.7071
      3         1.0000      0.0000      0.0000

Forum Regular
4:32:54 PM PDT - Wed, Jun 25th 2014
Talking about this with our experts I learned that the separation of the total shielding tensor into a diamagnetic and paramagnetic term is to some degree arbitrary. Only the total shielding tensor, the isotropic, and the anisotropic quantities have physical significance. As you can see those terms come out the same. The definition of diamagnetic and paramagnetic terms, however, has recently changed. I believe that change was related to something to do with the gauge invariance (but I could be wrong on that point).

Best wishes, Huub


Forum >> NWChem's corner >> Running NWChem