390. NWChem: "Fix collapse/expand in xc_nucder_gen" when using actlist and frequency calc

I ended up with the error in the title when running a job today -- the job does a calc on a complex system, where I'm only interested in optimising and doing normal mode analysis on part of the system. I (am supposed to) achieve this using actlist.

Googling shows the following two posts:
Question: http://www.emsl.pnl.gov/docs/nwchem/nwchem-support/2004/03/0010.Problem_with_Vibrational_frequencies
Answer: http://www.emsl.pnl.gov/docs/nwchem/nwchem-support/2004/03/0012.RE:_Problem_with_Vibrational_frequencies

There doesn't seem to be a follow up between then and now (nine years later).

Note that this is completely different from this post which seems to be due to the OP not putting the actlist in deck A2 (i.e. s/he computes the frequencies for different systems in A and B)

The issue:
The smallest job I can use to trigger this is the following:

scratch_dir /scratch
Title "trigger"

Start  trigger

echo

charge 0
geometry autosym units angstrom
 O     -1.79757     -0.236404     0.264744
 H     -2.77094     -0.198938     0.271156
 H     -1.44498     0.380299     -0.412669
 O     -0.0458786     1.77247     -0.0605460
 H     -0.317322     1.32556     0.781197
 H     0.894700     1.49387     -0.199468
end

ecce_print ecce.out

basis "ao basis" cartesian print
  H library "6-31G"
  O library "6-31G"
END

dft
  mult 1
  direct
  XC b3lyp
  grid fine
  iterations 999 
  mulliken
end

task dft energy
set geometry:actlist 1:3 
task dft freq

Here's the error message:

HESSIAN: the one electron contributions are done in       0.1s

 Fix collapse/expand in xc_nucder_gen        0
 ------------------------------------------------------------------------
 ------------------------------------------------------------------------
  current input line : 
     0:   
 ------------------------------------------------------------------------
 ------------------------------------------------------------------------
 This error has not yet been assigned to a category
 ------------------------------------------------------------------------
 For more information see the NWChem manual at http://www.nwchem-sw.org/index.php/NWChem_Documentation

I can't trigger it with just one water and using an actlist of 1:3 -- seems like I need inactive atoms. However, I've run similar calcs in the past without issue.

The solution:
Do a numerical instead of analytical frequency calculation.

This works:

scratch_dir /scratch
Title "trigger"

Start  trigger

echo

charge 0
geometry autosym units angstrom
 O     -1.79757     -0.236404     0.264744
 H     -2.77094     -0.198938     0.271156
 H     -1.44498     0.380299     -0.412669
 O     -0.0458786     1.77247     -0.0605460
 H     -0.317322     1.32556     0.781197
 H     0.894700     1.49387     -0.199468
end

ecce_print ecce.out

basis "ao basis" cartesian print
  H library "6-31G"
  O library "6-31G"
END

dft
  mult 1
  direct
  XC b3lyp
  grid fine
  iterations 999 
  mulliken
end

task dft energy
set geometry:actlist 1:3 
task dft freq numerical

389. Patches for NWChem 6.1.1 on Debian Wheezy/Testing

There are a couple of issues with the current version of NWChem (27th of June 2012):
* PSPW is broken when NWChem is compiled with/run on systems with gcc 4.7 (here)
* Python support requires patching to include -lz -lssl (here)
* for GabEdit to work more detail needs to be printed (here)

---

To fix all those issues in one go, do the following:

1. Copy the text at the end of the post, and paste it into a file, e.g. diff.patch.

2. Put the patch file in NWCHEM_TOP (i.e. the root of the source code) e.g. /opt/nwchem/nwchem-6.1.1
If you are patching a previously compiled version of nwchem then do

patch -p0 < diff.patch
export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=`pwd`
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all python"
export PYTHONVERSION=2.7
export PYTHONHOME=/usr
export BLASOPT="-L/opt/openblas/lib -lopenblas"
export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export MPI_LOC=/usr/lib/openmpi/lib
export MPI_INCLUDE=/usr/lib/openmpi/include
export LIBRARY_PATH=$LIBRARY_PATH:/usr/lib/openmpi/lib:/opt/openblas/lib
export LIBMPI="-lmpi -lopen-rte -lopen-pal -ldl -lmpi_f77 -lpthread"
export FC=gfortran
cd $NWCHEM_TOP/src/ddscf
make
cd $NWCHEM_TOP/src/nwdft/scf_dft
make
cd $NWCHEM_TOP/src/mcscf
make
cd $NWCHEM_TOP/src
make link
cd $NWCHEM_TOP/contrib
./getmem.nwchem

If it's a freshly extracted source, otherwise look at http://verahill.blogspot.com.au/2012/09/briefly-compiling-nwchem-611-with.html

---

Patch:

diff -rupN src.original/config/makefile.h src/config/makefile.h
--- src.original/config/makefile.h 2013-04-15 12:41:45.016853322 +1000
+++ src/config/makefile.h 2013-04-15 12:38:44.933319544 +1000
@@ -1169,7 +1169,7 @@ endif
         FOPTIONS   = -Wextra #-Wunused #-ffast-math
         FOPTIMIZE  = -O2 -ffast-math -Wuninitialized 
         DEFINES  += -DGFORTRAN
-        _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6./ {print "Y";exit};/4.7./ {print "Y";exit};{print "N"}')
+        _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6/ {print "Y";exit};/4.7/ {print "Y";exit};{print "N"}')
         ifeq ($(_GCC46),Y) 
         DEFINES  += -DGCC46
     endif
@@ -1298,7 +1298,7 @@ endif
          FVECTORIZE=-O3 -ffast-math -mtune=native -mfpmath=sse -msse3 -ftree-vectorize -ftree-vectorizer-verbose=1   -fprefetch-loop-arrays  -funroll-all-loops 
 #         FOPTIMIZE=-O1
 #         FVECTORIZE=-O1
-       _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6./ {print "Y";exit};/4.7./ {print "Y";exit};{print "N"}')
+       _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6/ {print "Y";exit};/4.7/ {print "Y";exit};{print "N"}')
        ifeq ($(_GCC46),Y) 
          DEFINES  += -DGCC46
        endif
@@ -1890,7 +1890,7 @@ endif
           FOPTIONS +=  -ff2c -fno-second-underscore
         endif
         DEFINES  += -DCHKUNDFLW -DGCC4
-        _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6./ {print "Y";exit};/4.7./ {print "Y";exit};{print "N"}')
+        _GCC46= $(shell gfortran -dumpversion  2>&1|awk ' /4.6/ {print "Y";exit};/4.7/ {print "Y";exit};{print "N"}')
         ifeq ($(_GCC46),Y) 
           DEFINES  += -DGCC46
         endif
@@ -1954,7 +1954,7 @@ endif
 
      ifeq ($(BUILDING_PYTHON),python)
 #   EXTRA_LIBS += -ltk -ltcl -L/usr/X11R6/lib -lX11 -ldl
-     EXTRA_LIBS +=    -lnwcutil  -lpthread -lutil -ldl
+     EXTRA_LIBS +=    -lnwcutil  -lpthread -lutil -ldl -lssl -lz
   LDOPTIONS = -Wl,--export-dynamic 
      endif
 ifeq ($(NWCHEM_TARGET),CATAMOUNT)
diff -rupN src.original/ddscf/movecs_pr_anal.F src/ddscf/movecs_pr_anal.F
--- src.original/ddscf/movecs_pr_anal.F 2013-04-15 12:41:45.036852381 +1000
+++ src/ddscf/movecs_pr_anal.F 2013-04-15 12:23:28.100409225 +1000
@@ -195,7 +195,7 @@ c
  22         format(1x,2('  Bfn.  Coefficient  Atom+Function  ',5x))
             write(LuOut,23)
  23         format(1x,2(' ----- ------------  ---------------',5x))
-            do klo = 0, min(n-1,9), 2
+            do klo = 0, min(n-1,199), 2
                khi = min(klo+1,n-1)
                write(LuOut,2) (
      $              int_mb(k_list+k)+1, 
diff -rupN src.original/ddscf/rohf.F src/ddscf/rohf.F
--- src.original/ddscf/rohf.F 2013-04-15 12:41:45.036852381 +1000
+++ src/ddscf/rohf.F 2013-04-15 12:23:28.100409225 +1000
@@ -153,7 +153,7 @@ c
             ilo = 1
             ihi = nmo
          endif
-         call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs, 
+         call movecs_print_anal(basis, ilo, ihi, 0.01d0, g_movecs, 
      $        'ROHF Final Molecular Orbital Analysis', 
      $        .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
      $        .true., dbl_mb(k_occ))
diff -rupN src.original/ddscf/scf_vec_guess.F src/ddscf/scf_vec_guess.F
--- src.original/ddscf/scf_vec_guess.F 2013-04-15 12:41:45.036852381 +1000
+++ src/ddscf/scf_vec_guess.F 2013-04-15 12:23:28.100409225 +1000
@@ -505,19 +505,19 @@ c
          nprint = min(nclosed+nopen+30,nmo)
          if (scftype.eq.'RHF' .or. scftype.eq.'ROHF') then
             call movecs_print_anal(basis, 1,
-     &           nprint, 0.15d0, g_movecs, 
+     &           nprint, 0.01d0, g_movecs, 
      &           'ROHF Initial Molecular Orbital Analysis', 
      &           .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
      &           .true., dbl_mb(k_occ))
          else
             nprint = min(nalpha+20,nmo)
             call movecs_print_anal(basis, max(1,nbeta-20),
-     &           nprint, 0.15d0, g_movecs, 
+     &           nprint, 0.01d0, g_movecs, 
      &           'UHF Initial Alpha Molecular Orbital Analysis', 
      &           .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
      &           .true., dbl_mb(k_occ))
             call movecs_print_anal(basis, max(1,nbeta-20),
-     &           nprint, 0.15d0, g_movecs(2), 
+     &           nprint, 0.01d0, g_movecs(2), 
      &           'UHF Initial Beta Molecular Orbital Analysis', 
      &           .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
      &           .true., dbl_mb(k_occ+nbf))
diff -rupN src.original/ddscf/uhf.F src/ddscf/uhf.F
--- src.original/ddscf/uhf.F 2013-04-15 12:41:45.036852381 +1000
+++ src/ddscf/uhf.F 2013-04-15 12:23:28.096409414 +1000
@@ -144,11 +144,11 @@ C
          enddo
          ihi = max(ihi-1,1)
  9611    continue
-         call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs, 
+         call movecs_print_anal(basis, ilo, ihi, 0.01d0, g_movecs, 
      $        'UHF Final Alpha Molecular Orbital Analysis', 
      $        .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
      $        .true., dbl_mb(k_occ))
-         call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs(2), 
+         call movecs_print_anal(basis, ilo, ihi, 0.01d0, g_movecs(2), 
      $        'UHF Final Beta Molecular Orbital Analysis', 
      $        .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
      $        .true., dbl_mb(k_occ+nbf))
diff -rupN src.original/mcscf/mcscf.F src/mcscf/mcscf.F
--- src.original/mcscf/mcscf.F 2013-04-15 12:41:45.000854073 +1000
+++ src/mcscf/mcscf.F 2013-04-15 12:23:23.748613695 +1000
@@ -680,7 +680,7 @@ c
       if (util_print('final vectors analysis', print_default))
      $     call movecs_print_anal(basis, 
      $     max(1,nclosed-10), min(nbf,nclosed+nact+10),
-     $     0.15d0, g_movecs, 'Analysis of MCSCF natural orbitals',
+     $     0.01d0, g_movecs, 'Analysis of MCSCF natural orbitals',
      $     .true., dbl_mb(k_evals), .true., int_mb(k_sym), 
      $     .true., dbl_mb(k_occ))
 c     
diff -rupN src.original/nwdft/scf_dft/dft_mxspin_ovlp.F src/nwdft/scf_dft/dft_mxspin_ovlp.F
--- src.original/nwdft/scf_dft/dft_mxspin_ovlp.F 2013-04-15 12:41:45.604825677 +1000
+++ src/nwdft/scf_dft/dft_mxspin_ovlp.F 2013-04-15 12:23:28.228403211 +1000
@@ -184,14 +184,14 @@ c
       call ga_sync()
 c
       call movecs_print_anal(basis,int_mb(k_non),int_mb(k_non)
-     & ,0.15d0,g_alpha,'Alpha Orbitals without Beta Partners',
+     & ,0.01d0,g_alpha,'Alpha Orbitals without Beta Partners',
      &   .false., 0.0 ,.false., 0 , .false., 0 )
 c
       if (nct.GE.2) then
       do i = 2,nct
       ind = int_mb(k_non+i-1)
       call movecs_print_anal(basis,ind,ind
-     & ,0.15d0,g_alpha,' ',
+     & ,0.01d0,g_alpha,' ',
      &   .false., 0.0 ,.false., 0 , .false., 0 )
       enddo
       endif
@@ -350,7 +350,7 @@ c      endif
 c      endif
 c 9990 format(/,18x,'THERE ARE',i3,1x,'UN-PARTNERED ALPHA ORBITALS')
 c
-       call movecs_print_anal(basis, 1, nalp, 0.15d0, g_ualpha,
+       call movecs_print_anal(basis, 1, nalp, 0.01d0, g_ualpha,
      & 'Alpha Orb. w/o Beta Partners (after maxim. alpha/beta overlap)',
      &   .false., 0.0 ,.false., 0 , .false., 0 )
 c
diff -rupN src.original/nwdft/scf_dft/dft_scf.F src/nwdft/scf_dft/dft_scf.F
--- src.original/nwdft/scf_dft/dft_scf.F 2013-04-15 12:41:45.608825490 +1000
+++ src/nwdft/scf_dft/dft_scf.F 2013-04-15 12:23:28.228403211 +1000
@@ -1733,7 +1733,7 @@ c
             else
                blob='DFT Final Beta Molecular Orbital Analysis' 
             endif
-            call movecs_print_anal(ao_bas_han, ilo, ihi, 0.15d0, 
+            call movecs_print_anal(ao_bas_han, ilo, ihi, 0.01d0, 
      &           g_movecs(ispin), 
      &           blob, 
      &           .true., dbl_mb(k_eval(ispin)), oadapt, 
diff -rupN src.original/nwdft/scf_dft_cg/dft_cg_solve.F src/nwdft/scf_dft_cg/dft_cg_solve.F
--- src.original/nwdft/scf_dft_cg/dft_cg_solve.F 2013-04-15 12:41:45.612825303 +1000
+++ src/nwdft/scf_dft_cg/dft_cg_solve.F 2013-04-15 12:23:28.220403588 +1000
@@ -164,7 +164,7 @@ c
             blob = 'DFT Final Beta Molecular Orbital Analysis'
           endif
           call movecs_fix_phase(g_movecs(ispin))
-          call movecs_print_anal(basis, ilo, ihi, 0.15d0,
+          call movecs_print_anal(basis, ilo, ihi, 0.01d0,
      &         g_movecs(ispin),blob,
      &         .true., dbl_mb(k_eval+(ispin-1)*nbf),
      &         oadapt, int_mb(k_irs+(ispin-1)*nbf),

388. NWChem, PSPW and Fortran runtime error

I'm following this post: http://www.emsl.pnl.gov/docs/nwchem/nwchem-support/2012/12/0024.Re:_NWCHEM_errors_in_running_nwchem-6.1.

The issue:
If, on debian testing/wheezy 64 bit with gcc 4.7,  you run the following using NWChem 6.1.1

Title "boric acid"
Start  boric

echo

charge 0

geometry autosym units angstrom
 B     0.00000     0.00000     0.00000
 O     -4.93432e-17     1.55000     0.00000
 H     1.06537     1.92667     0.00000
 O     1.34234     -0.775000     0.00000
 H     1.13586     -1.88597     0.00000
 O     -1.34234     -0.775000     0.00000
 H     -2.20123     -0.0406974     0.00000
end

ecce_print ecce.out

nwpw
  mult 1
  xc pbe96
  cutoff 90.0
  np_dimensions -1  -1  
  tolerances 1e-7  1e-7
  car-parrinello
    nose-hoover 1.000000e+03 2.981500e+02 1.000000e+03 2.981500e+02
    time_step 0.500000e+00
    fake_mass 1.000000e+02
    loop 10 100 
    scaling 1.000000e+00 1.000000e+00
  end 
end

task pspw energy
task pspw car-parrinello

You'll end up with

****************************************************
          *                                                  *
          *               NWPW PSPW Calculation              *
          *                                                  *
          *  [ (Grassman/Stiefel manifold implementation) ]  *
          *                                                  *
          *      [ NorthWest Chemistry implementation ]      *
          *                                                  *
          *            version #5.10   06/12/02              *
          *                                                  *
          *    This code was developed by Eric J. Bylaska,   *
          *    and was based upon algorithms and code        *
          *    developed by the group of Prof. John H. Weare *
          *                                                  *
          ****************************************************
     >>>  JOB STARTED       AT Mon Apr 15 10:05:42 2013  <<<
          ================ input data ========================
  library name resolved from: compiled reference
  NWCHEM_NWPW_LIBRARY set to: 
Generating 1d pseudopotential for B
At line 649 of file psp_generator_input.F (unit = 99, file = './junk.inp')
Fortran runtime error: Sequential READ or WRITE not allowed after EOF marker, possibly use REWIND or BACKSPACE
--------------------------------------------------------------------------
mpirun has exited due to process rank 0 with PID 28100 on
node neon exiting without calling "finalize". This may
have caused other processes in the application to be
terminated by signals sent by mpirun (as reported here).
-------------------------------------------------------------------------

The fix
/opt/nwchem/nwchem-6.1.1-src is the NWCHEM_TOP of my particular nwchem source. See e.g. here for build instructions in general and not that e.g. BLASOPT settings etc. correspond to what I need on my system. Modify as necessary.

cd /opt/nwchem/nwchem-6.1.1-src/
wget http://www.nwchem-sw.org/images/Makefile.h.gcc46.patch.gz
gunzip Makefile.h.gcc46.patch.gz
patch -p0 < Makefile.h.gcc46.patch
cd src/nwpw/
touch `egrep -l GCC46 */*/*F`

export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=/opt/nwchem/nwchem-6.1.1-src
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all"
export PYTHONHOME=/usr
export BLASOPT="-L/opt/acml/acml5.2.0/gfortran64_fma4_int64/lib -lacml"
export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export MPI_LOC=/usr/lib/openmpi/lib
export MPI_INCLUDE=/usr/lib/openmpi/include
export LIBRARY_PATH="$LIBRARY_PATH:/usr/lib/openmpi/lib:/opt/acml/acml5.2.0/gfortran64_fma4_int64/lib"
export LIBMPI="-lmpi -lopen-rte -lopen-pal -ldl -lmpi_f77 -lpthread"
export FC=gfortran
make 
cd ../
make link
./getmem.nwchem

Because we just relink instead of recompiling from scratch, recompiling is fairly quick (2-3 minutes in total)

Testing

It now continues as it should:

================ input data ========================
  library name resolved from: compiled reference
  NWCHEM_NWPW_LIBRARY set to: 
Generating 1d pseudopotential for B   

and everything is good from there-on.


NOTE: you need to delete the output files from any previous run, or you might end up with errors during the Car-Parrinello portion (e.g.

At line 170 of file MOTION.F (unit = 19, file = './boric.ion_motion')
Fortran runtime error: Sequential READ or WRITE not allowed after EOF marker, possibly use REWIND or BACKSPACE

)

Stay tuned...
What remains for me is to start collecting the patches that are needed for nwchem to work properly (from my perspective -- I only have to support myself and my debian boxes, not a heterogeneous set of users with equally heterogeneous hardware, so my job is easier than that of the nwchem devs) so that installing and upgrading nwchem becomes less painful (yes, I've been patching by hand up till now...).

The issues that need to be patched are:
* including libz and libssl on debian when building with python support (e.g. step 1 in this post)
* changing the details in the output to support gabEdit (e.g. this post)
* patching for GCC 4.7 (i.e. what you've just read)

Update: The post is here now: http://verahill.blogspot.com.au/2013/04/389-patches-for-nwchem-611-on-debian.html