Merge pull request #1282 from nwchemgit/dependabot/github_actions/actions/setup-python-7

Bump actions/setup-python from 6 to 7

.github/workflows/check_64_to_32.yml

@@ -19,7 +19,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout code
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
         with:
           show-progress:
           fetch-depth: 1

.github/workflows/docker_actions.yml

@@ -22,20 +22,20 @@ jobs:
     timeout-minutes: 420
     steps:
       - name: Checkout
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
         with:
           repository: 'nwchemgit/nwchem-dockerfiles'
           fetch-depth: 5
       - name: Setup cache
         id: setup-cache
-        uses: actions/cache@v4
+        uses: actions/cache@v6
         with:
           path: |
             ~/cache
           key: ${{ matrix.folder}}-${{ matrix.fc}}/${{ matrix.archs}}-nwchem-dockeractions-v002
       - name: Qemu
         id: qemu
-        uses: docker/setup-qemu-action@v3
+        uses: docker/setup-qemu-action@v4
         with:
           image: tonistiigi/binfmt:qemu-v8.1.5
       - name: Available platforms
@@ -43,7 +43,7 @@ jobs:
           echo ${{ steps.qemu.outputs.platforms }}
           docker images
       - name: Set up Docker Buildx
-        uses: docker/setup-buildx-action@v3
+        uses: docker/setup-buildx-action@v4
         with:
           driver-opts: |
               image=moby/buildkit:v0.18.2
@@ -62,7 +62,7 @@ jobs:
         run: |
           echo "cache_hit=N" >> $GITHUB_ENV
       - name: build_schedule
-        uses: docker/build-push-action@v6
+        uses: docker/build-push-action@v7
         with:
           push: false
           context: ${{ matrix.folder }}

.github/workflows/github_actions.yml

@@ -21,7 +21,7 @@ jobs:
       matrix:
         blas: [build_openblas]
         blas_size: [8]
-        os: [ubuntu-22.04, macos-latest, ubuntu-24.04-arm]
+        os: [ubuntu-latest, macos-latest, ubuntu-24.04-arm]
         experimental: [false]
         mpi_impl:
           - openmpi
@@ -36,6 +36,9 @@ jobs:
           - "tce"
         fc:
           - gfortran
+        exclude:
+          - os: macos-latest
+            mpi_impl: mpich
         include:
           - os: ubuntu-latest
             experimental: true
@@ -252,7 +255,7 @@ jobs:
             mpi_impl: openmpi
             armci_network: MPI-TS
             nwchem_modules: "tinyqmpw"
-            fc: gfortran-12
+            fc: gfortran-13
             cc: gcc
             use_libxc: 1
             xcode_version: 15.3
@@ -261,7 +264,7 @@ jobs:
             mpi_impl: openmpi
             armci_network: MPI-TS
             nwchem_modules: "tinyqmpw python"
-            fc: gfortran-13
+            fc: gfortran-14
             use_libxc: 1
             xcode_version: 15.2
           - os: ubuntu-22.04
@@ -301,7 +304,7 @@ jobs:
             mpi_impl: openmpi
             armci_network: MPI-PR
             nwchem_modules: "qmandpw qmd"
-            fc: gfortran-12
+            fc: gfortran-13
             xcode_version: 15.0.1
           - os: macos-latest
             experimental: true
@@ -315,16 +318,16 @@ jobs:
             experimental: true
             mpi_impl: openmpi
             armci_network: MPI-TS
-            nwchem_modules: "tinyqmpw"
-            fc: gfortran-12
-            cc: gcc-12
+            nwchem_modules: "tinyqmpw python"
+            fc: gfortran-14
+            cc: gcc-14
             xcode_version: 15.4
           - os: macos-latest
             experimental: true
             mpi_impl: openmpi
             armci_network: MPI-TS
             nwchem_modules: "tinyqmpw"
-            fc: gfortran-13
+            fc: gfortran-14
             blas: "accelerate"
             blas_size: 4
           - os: ubuntu-22.04
@@ -411,28 +414,9 @@ jobs:
             experimental: true
             mpi_impl: openmpi
             armci_network: MPI-TS
-            nwchem_modules: "nwdft driver solvation"
-            fc: gfortran
-            cc: gcc
-            blas: "brew_openblas"
-            blas_size: 4
-          - os: macos-15
-            experimental: true
-            mpi_impl: openmpi
-            armci_network: MPI-TS
-            nwchem_modules: "nwdft driver solvation"
-            fc: gfortran
-            cc: gcc
-            blas: "brew_openblas"
-            blas_size: 4
-            xcode_version: 15.0.1
-          - os: macos-latest
-            experimental: true
-            mpi_impl: mpich
-            armci_network: MPI-TS
             nwchem_modules: "all"
-            fc: gfortran-14
-            cc: gcc-14
+            fc: gfortran-15
+            cc: gcc-15
             blas: "build_openblas"
             blas_size: 8
       fail-fast: false
@@ -474,12 +458,12 @@ jobs:
           fi
         shell: bash
       - name: Checkout code
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
         with:
           show-progress:
           fetch-depth: 200
       - name: Install Python
-        uses: actions/setup-python@v6
+        uses: actions/setup-python@v7
         with:
           python-version: "3.12"
       - name: Get Instruction
@@ -488,14 +472,17 @@ jobs:
           ./travis/cpuinfo.sh
           echo "simd=$(./travis/guess_simd.sh)" >> $GITHUB_ENV
           echo "microarch=$(cat /tmp/microarch_$(id -u).txt)" >> $GITHUB_ENV
+          if [[ ! -z $USE_OPENMP ]]; then
+          echo "openmp=-openmp$USE_OPENMP" >> $GITHUB_ENV
+          fi
         shell: bash
       - name: Setup cache
         id: setup-cache
-        uses: actions/cache@v4
+        uses: actions/cache@v6
         with:
           path: |
             ~/cache
-          key: ${{ matrix.os }}-${{ matrix.mpi_impl}}-${{ matrix.fc}}-xcode${{ matrix.xcode_version}}-simint${{ matrix.use_simint}}-blas${{ matrix.blas}}-blas_size${{ matrix.blas_size}}-libxc${{matrix.use_libxc}}-tblite${{matrix.use_tblite}}-simd${{ env.simd }}${{ env.microarch}}-nwchemcache-v003
+          key: ${{ matrix.os }}-${{ matrix.mpi_impl}}-${{ matrix.fc}}-xcode${{ matrix.xcode_version}}-simint${{ matrix.use_simint}}-blas${{ matrix.blas}}-blas_size${{ matrix.blas_size}}-libxc${{matrix.use_libxc}}-tblite${{matrix.use_tblite}}-simd${{ env.openmp }}${{ env.simd }}${{ env.microarch}}-nwchemcache-v003
       - name: build environment
         run: |
           pwd
@@ -528,8 +515,10 @@ jobs:
         id: qa_test
         if: steps.compile.conclusion == 'success'
         run: |
+          if [[ $FC != 'nvfortran' && ${{ env.microarch}} != 'Zen4' ]]; then
           ./travis/run_qas.sh
+          fi
       - name: Check if QA testing has failed
-        if: ${{ failure() }} &&  steps.qa_test.outcome == 'failure'
+        if: ${{ failure() &&  steps.qa_test.outcome == 'failure' }}
         run: |
           ./travis/check_qas.sh

.github/workflows/lint_python.yml

@@ -11,8 +11,8 @@ jobs:
   lint_python:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v5
-      - uses: actions/setup-python@v6
+      - uses: actions/checkout@v7
+      - uses: actions/setup-python@v7
         with:
           python-version: '3.12.3'
       - run: pip install --upgrade pip wheel

.github/workflows/nwchem_releasetar.yaml

@@ -10,7 +10,11 @@ on:
         type: choice
         options:
           - nightly
+          - 7.2.4
           - 7.3.0
+          - 7.3.1
+          - 7.3.2
+          - 7.3.3
       branch:
         description: 'branch of NWChem repository'
         required: true
@@ -18,16 +22,16 @@ on:
         type: choice
         options:
           - master
-          - hotfix/7.2.0
+          - hotfix/release-7-2-0
           - release-7-3-0
-          - hotfix/7.3.0
+          - hotfix/release-7-3-0
 
 jobs:
   do_tar:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout code
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
         with:
           show-progress:
           fetch-depth: 200
@@ -69,13 +73,13 @@ jobs:
         run: |
           ls -lrt ${{ env.tmpdir_name }}/*tar*
       - name: Upload artifacts
-        uses: actions/upload-artifact@v4
+        uses: actions/upload-artifact@v7
         with:
           name: nwchem-tarball
           path: ${{ env.tmpdir_name }}/*tar*
       - if: github.event.inputs.release_version == 'nightly'
         name: Nightly
-        uses: softprops/action-gh-release@v2
+        uses: softprops/action-gh-release@v3
         with:
           prerelease: true
           name: nightly
@@ -84,7 +88,7 @@ jobs:
           fail_on_unmatched_files: true
       - if: github.event.inputs.release_version != 'nightly'
         name: Release 
-        uses: softprops/action-gh-release@v2
+        uses: softprops/action-gh-release@v3
         with:
           prerelease: true
           draft: true

QA/doqmtests.mpi

@@ -92,6 +92,8 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np dft_siosi3
 let "myexit+=$?"
+./runtests.mpi.unix procs $np adft_he2+
+let "myexit+=$?"
 ./runtests.mpi.unix procs $np adft_siosi3
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np dft_cr2
@@ -192,6 +194,8 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np tddft_ac_co
 let "myexit+=$?"
+./runtests.mpi.unix procs $np h2o_bnl2007_tddft
+let "myexit+=$?"
 ./runtests.mpi.unix procs $np dft_xdm1
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np natom_iorb
@@ -342,6 +346,10 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np cho_bp_props
 let "myexit+=$?"
+./runtests.mpi.unix procs $np dft_hyperf_zno
+let "myexit+=$?"
+./runtests.mpi.unix procs $np dft_hyperf_zno6bq
+let "myexit+=$?"
 ./runtests.mpi.unix procs $np prop_cg_nh3_b3lyp
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np localize-ibo-aa
@@ -603,13 +611,15 @@ let "myexit+=$?"
 ./runtests.mpi.unix procs $np bare_ecps
 let "myexit+=$?"
 if [[ ! -z "${BUILD_ELPA}" ]]; then
-./runtests.mpi.unix procs $np siosi3_elpa
+./runtests.mpi.unix procs $np dft_siosi3_elpa
 let "myexit+=$?"
 fi
 ./runtests.mpi.unix procs $np dft_formaldehyde_mom
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np dft_boron_imom
 let "myexit+=$?"
+./runtests.mpi.unix procs $np dft_napht_masvp
+let "myexit+=$?"
 #--- if we only want quick tests get out now!
 if [[ "$what" == "fast" ]]; then
     echo
@@ -718,6 +728,24 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np udft_aoresp
 let "myexit+=$?"
+./runtests.mpi.unix procs $np water/water_md
+let "myexit+=$?"
+./runtests.mpi.unix procs $np water/water_pme
+let "myexit+=$?"
+./runtests.mpi.unix procs $np ethanol/ethanol_md
+let "myexit+=$?"
+./runtests.mpi.unix procs $np ethanol/ethanol_ti
+let "myexit+=$?"
+./runtests.mpi.unix procs $np prep/aal
+let "myexit+=$?"
+./runtests.mpi.unix procs $np prep/bpx
+let "myexit+=$?"
+./runtests.mpi.unix procs $np prep/caa
+let "myexit+=$?"
+./runtests.mpi.unix procs $np prep/fsc
+let "myexit+=$?"
+./runtests.mpi.unix procs $np trypsin/trypsin_md
+let "myexit+=$?"
 echo 'The following  tests require few processors but a lot of time'
 ./runtests.mpi.unix procs $np ch4cl_zts
 let "myexit+=$?"
@@ -733,7 +761,7 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np H2COWat_VEMExc_CISHF
 let "myexit+=$?"
-./runtests.mpi.unix procs $np p2ta-vem
+./runtests.mpi.unix procs $np p2ta-vem-ri
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np bsse_sodft
 let "myexit+=$?"
@@ -743,6 +771,8 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np  rt_tddft_dplot
 let "myexit+=$?"
+./runtests.mpi.unix procs $np  metagga_li4h4_1_m05_0
+let "myexit+=$?"
 if [[ "$what" != "flaky" ]]; then
     echo
     echo  "the number of failed tests is" $myexit

QA/runtests.mpi.unix

@@ -243,9 +243,11 @@ do
      md_test_run=0
   fi
 #MV
-
+  sync
   if [ ! -f $TESTDIR/$INPUTFILE ]; then
-    echo failed: cannot access $TESTDIR/$INPUTFILE
+      echo failed: cannot access $TESTDIR/$INPUTFILE
+      echo directory content
+      printf "$(ls -lrt $TESTDIR) \n"
     overall_status=1
     continue
   fi
@@ -337,6 +339,7 @@ sync
         runstatus=$?
       else
           $MPIRUN $NPOPT $NPROC $NWCHEM $INPUTFILE 2> $ERRORFILE 1> $OUTPUTFILE
+	  sync
         runstatus=$?
       fi
     else

QA/tests/cosmo_h2cco2_opt/cosmo_h2cco2_opt.nw

@@ -69,6 +69,7 @@ end
 
 cosmo
  do_gasphase f
+ do_cosmo_ks
  cavity ses
  ndiv xfine
 end

QA/tests/dft_hyperf_zno6bq/dft_hyperf_zno6bq.nw

@@ -2,37 +2,54 @@ echo
 start dft_zorahyperf_zno6bq
 title dft_zorahyperf_zno6bq
 
-
+charge -10
 geometry 
-symmetry c4
+#symmetry c4
 zn 0. 0. 0.
 O  2. 0. 0
-O -2. 0. 0
+O1 -2. 0. 0
 O  0. 2. 0.
-O  0. -2. 0.
+O1  0. -2. 0.
 O  0. 0. 2.1
-O  0. 0. -2.1 
-bqZn  4.  0.  0. charge 1d-20
-bqZn -4.  0.  0. charge 1d-20
-bqZn  0.  4.  0. charge 1d-20
-bqZn  0. -4.  0. charge 1d-20
-bqZn  0.  0. -4.2 charge 1d-20
-bqZn  0.  0.  4.2 charge 1d-20
+O1  0. 0. -2.1
+bq_zn  4.  0.  0.  charge 1d-20
+bq_zn -4.  0.  0.  charge 1d-20
+bq_zn  0.  4.  0.  charge 1d-20
+bq_zn  0. -4.  0.  charge 1d-20
+bq_zn  0.  0. -4.2  charge 1d-20
+bq_zn  0.  0.  4.2  charge 1d-20
 end
-charge -10
-BASIS "ao basis" spherical 
-* library def2-svp
+
+basis bse spherical
+zn library stuttgart_rsc_1997
+O library def2-svp
 end
-set tolguess 1d-3
-charge 0
+
+ecp bse
+zn library stuttgart_rsc_1997
+bq_zn nelec 10
+bq_zn ul
+2       1.000000000            0.000000000
+bq_zn S
+2      34.150000000          399.987282000
+2      14.590000000           85.485655000
+bq_zn P
+2      39.780000000          277.148960000
+2      14.950000000           69.052205000
+bq_zn D
+2      43.800000000          -34.149349000
+2      14.980000000           -3.291831000
+end
+
+
 dft
  direct
  mult 3
  xc pbe0
- grid xfine
- tolerances tight
+ grid fine
  convergence energy 1d-8 lshift 0.
- end
+ maxiter 99
+end
 task dft
 set dft:no_prune t
 relativistic

QA/tests/dft_napht_masvp/dft_napht_masvp.nw

@@ -0,0 +1,34 @@
+ echo
+start 
+ title "naphthalene in ma-SVP basis set" 
+
+ geometry# units au  
+  C      2.4044      0.7559      0.0000
+  C      2.4328     -0.6584      0.0000
+  C      1.2672     -1.3753      0.0000
+  C      0.0142     -0.7050      0.0000
+  C     -0.0142      0.7048      0.0000
+  C      1.2108      1.4252      0.0000
+  C     -1.2672      1.3754      0.0000
+  C     -2.4328      0.6585      0.0000
+  C     -2.4043     -0.7558      0.0000
+  C     -1.2108     -1.4254      0.0000
+  H      3.3509      1.3062      0.0000
+  H      3.4006     -1.1703      0.0000
+  H      1.2810     -2.4710      0.0000
+  H      1.1803      2.5206      0.0000
+  H     -1.2808      2.4710      0.0000
+  H     -3.4008      1.1701      0.0000
+  H     -3.3508     -1.3060      0.0000
+  H     -1.1805     -2.5207      0.0000
+ end  
+ 
+ basis "ao basis" spherical
+   * library ma-svp
+ end
+
+ dft
+   xc M06-2X
+ end
+
+ task dft optimize

QA/tests/ethanol/ethanol_ti.tst

@@ -2,29 +2,29 @@ Iteration        =          0
 Energy           =    0.000E+00
 
 Iteration        =          0
-Energy           =   -8.957E+03
+Energy           =   -8.944E+03
 
 Iteration        =          0
-Energy           =   -8.961E+03
+Energy           =   -8.948E+03
 
 Iteration        =          0
-Energy           =   -8.955E+03
+Energy           =   -8.942E+03
 
 Iteration        =          0
-Energy           =   -8.957E+03
+Energy           =   -8.942E+03
 
 Iteration        =          0
-Energy           =   -8.957E+03
+Energy           =   -8.945E+03
 
 Iteration        =          0
-Energy           =   -8.956E+03
+Energy           =   -8.942E+03
 
 Iteration        =          0
-Energy           =   -8.959E+03
+Energy           =   -8.946E+03
 
 Iteration        =          0
-Energy           =   -8.964E+03
+Energy           =   -8.951E+03
 
 Iteration        =          0
-Energy           =   -8.957E+03
+Energy           =   -8.943E+03
 

QA/tests/libxc_waterdimer_bmk/libxc_waterdimer_bmk.nw

@@ -52,6 +52,4 @@ bsse
 mon first 1 2 3
 mon second 4 5 6
 end
-
-task dft optimize
-
+task dft energy

QA/tests/metagga_li4h4_1_m05_0/metagga_li4h4_1_m05_0.nw

@@ -1,14 +1,14 @@
 echo
 start metagga_dat
 geometry
- li     0.01002691    -0.94366842     1.43225550
- h     -1.37286946     0.02420936     0.52629308
- h     -0.02940706     3.13404608    -1.03053445
- h      1.37316029     0.05159532     0.52604548
- li     1.94807292    -1.21726078    -0.71290509
- li    -0.01236479     1.28755446     0.09274502
- li    -1.92253586    -1.25587660    -0.71259553
- h     -0.03460612     3.40646322    -0.30519899
+ li    1.41886094    -0.97519294    -0.00213685
+ h     0.52227920     0.01448403    -1.37641561
+ h    -0.99294465     3.13015382    -0.16071498
+ h     0.52812846     0.02093224     1.36933229
+ li   -0.72646082    -1.23567921     1.93800611
+ li    0.10745306     1.27247387    -0.00525518
+ li   -0.73084646    -1.24632666    -1.93625430
+ h    -0.26448317     3.39416572     0.20706284
 end
 charge +1
 
@@ -36,4 +36,4 @@ driver
 end
 
 task dft optimize
-task dft frequencies numerical
+task dft frequencies

QA/tests/p2ta-vem-ri/p2ta-vem-ri.nw

@@ -0,0 +1,145 @@
+# COSMO-VEM solvation results:
+# excitation spectrum data:  GS = initial state, ES = final state
+# iteration #3
+# (1)          GS equilibrium total free energy =     -2419.3240621130
+# (2)               GS polarization free energy =        -0.0469479928 (   -1.2775 eV)
+# (3)                 GSRF ES total free energy =     -2419.2518703361
+# (4)          GSRF ES polarization free energy =        -0.0462179669 (   -1.2577 eV)
+# (5)          GSRF excitation energy (3) - (1) =         0.0721917769 (    1.9644 eV)
+# (6)                  VEM ES total free energy =     -2419.2519518530
+# (7)           VEM ES polarization free energy =        -0.0462461370 (   -1.2584 eV)
+# (8)        fast polarization component of (7) =        -0.0237342086 (   -0.6458 eV)
+# (9)                 1/2 * delV * delQdyn term =        -0.0000281700 (   -0.0008 eV)
+# (10) VEM vertical excitation energy (6) - (1) =         0.0721102600 (    1.9622 eV)
+#
+echo 
+title "p2ta_VEMExc_StressTest: Porphyrin in DMSO with RI fitting"
+
+start 
+
+
+geometry units angstrom
+symmetry d2h
+ H        -12.77626055    -1.35417648     0.00000000
+ H        -12.77626055     1.35417648     0.00000000
+ C        -11.92733377    -0.67419001     0.00000000
+ C        -11.92733377     0.67419001     0.00000000
+ H        -10.91550612    -3.16929606     0.00000000
+ H        -10.91550612     3.16929606     0.00000000
+ C        -10.53054773    -1.07949669     0.00000000
+ C        -10.53054773     1.07949669     0.00000000
+ C        -10.12799741    -2.41457338     0.00000000
+ C        -10.12799741     2.41457338     0.00000000
+ N         -9.70458289     0.00000000     0.00000000
+ H         -9.05535648    -5.10630304     0.00000000
+ H         -9.05535648     5.10630304     0.00000000
+ C         -8.82495359    -2.88813257     0.00000000
+ C         -8.82495359     2.88813257     0.00000000
+ C         -8.38692961    -4.24920097     0.00000000
+ C         -8.38692961     4.24920097     0.00000000
+ H         -7.70650435    -1.08831220     0.00000000
+ H         -7.70650435     1.08831220     0.00000000
+ N         -7.69744258    -2.10961304     0.00000000
+ N         -7.69744258     2.10961304     0.00000000
+ C         -7.02157966    -4.25203838     0.00000000
+ C         -7.02157966     4.25203838     0.00000000
+ C         -6.57821439    -2.89244519     0.00000000
+ C         -6.57821439     2.89244519     0.00000000
+ H         -6.35505884    -5.11038077     0.00000000
+ H         -6.35505884     5.11038077     0.00000000
+ N         -5.65839959     0.00000000     0.00000000
+ C         -5.26649554    -2.42970347     0.00000000
+ C         -5.26649554     2.42970347     0.00000000
+ C         -4.86446289    -1.10528304     0.00000000
+ C         -4.86446289     1.10528304     0.00000000
+ H         -4.46838689    -3.17297996     0.00000000
+ H         -4.46838689     3.17297996     0.00000000
+ C         -3.45147343    -0.71676164     0.00000000
+ C         -3.45147343     0.71676164     0.00000000
+ N         -2.38151207    -1.44748897     0.00000000
+ N         -2.38151207     1.44748897     0.00000000
+ C         -1.21080050    -0.71887652     0.00000000
+ C         -1.21080050     0.71887652     0.00000000
+ H          0.00000000    -2.49209217     0.00000000
+ C          0.00000000    -1.40280773     0.00000000
+ C          0.00000000     1.40280773     0.00000000
+ H          0.00000000     2.49209217     0.00000000
+ C          1.21080050    -0.71887652     0.00000000
+ C          1.21080050     0.71887652     0.00000000
+ N          2.38151207    -1.44748897     0.00000000
+ N          2.38151207     1.44748897     0.00000000
+ C          3.45147343    -0.71676164     0.00000000
+ C          3.45147343     0.71676164     0.00000000
+ H          4.46838689    -3.17297996     0.00000000
+ H          4.46838689     3.17297996     0.00000000
+ C          4.86446289    -1.10528304     0.00000000
+ C          4.86446289     1.10528304     0.00000000
+ C          5.26649554    -2.42970347     0.00000000
+ C          5.26649554     2.42970347     0.00000000
+ N          5.65839959     0.00000000     0.00000000
+ H          6.35505884    -5.11038077     0.00000000
+ H          6.35505884     5.11038077     0.00000000
+ C          6.57821439    -2.89244519     0.00000000
+ C          6.57821439     2.89244519     0.00000000
+ C          7.02157966    -4.25203838     0.00000000
+ C          7.02157966     4.25203838     0.00000000
+ N          7.69744258    -2.10961304     0.00000000
+ N          7.69744258     2.10961304     0.00000000
+ H          7.70650435    -1.08831220     0.00000000
+ H          7.70650435     1.08831220     0.00000000
+ C          8.38692961    -4.24920097     0.00000000
+ C          8.38692961     4.24920097     0.00000000
+ C          8.82495359    -2.88813257     0.00000000
+ C          8.82495359     2.88813257     0.00000000
+ H          9.05535648    -5.10630304     0.00000000
+ H          9.05535648     5.10630304     0.00000000
+ N          9.70458289     0.00000000     0.00000000
+ C         10.12799741    -2.41457338     0.00000000
+ C         10.12799741     2.41457338     0.00000000
+ C         10.53054773    -1.07949669     0.00000000
+ C         10.53054773     1.07949669     0.00000000
+ H         10.91550612    -3.16929606     0.00000000
+ H         10.91550612     3.16929606     0.00000000
+ C         11.92733377    -0.67419001     0.00000000
+ C         11.92733377     0.67419001     0.00000000
+ H         12.77626055    -1.35417648     0.00000000
+ H         12.77626055     1.35417648     0.00000000
+end
+
+basis "ao basis" spherical
+ * library 6-31G*
+END
+
+basis "cd basis" spherical
+* library "Weigend Coulomb Fitting"
+end
+
+driver
+ maxiter 100
+ xyz p2ta-opt
+end
+
+dft
+ maxiter 100
+ xc xcamb88 1.00 lyp 1.0 hfexch 1.00
+ cam 0.33 cam_alpha 0.0 cam_beta 1.0
+ direct
+end
+
+cosmo
+ do_cosmo_vem 1
+ solvent dmso
+end
+
+tddft
+ nroots 5
+ algorithm 1
+ notriplet
+ target 1
+ civecs
+ grad
+   root 1
+   solve_thresh 1d-05
+ end
+end
+task tddft gradient

QA/tests/prop_uhf_ch3f/prop_uhf_ch3f.nw

@@ -20,12 +20,14 @@ property
 end
 
 scf
+ direct
  uhf
 end
 
 task scf property
 
 dft
+# direct
  xc b3lyp
  odft
 end

README.md

@@ -13,11 +13,13 @@
 ![Github Downloads Release 7.2.1](https://img.shields.io/github/downloads/nwchemgit/nwchem/v7.2.1-release/total)  
 ![Github Downloads Release 7.2.2](https://img.shields.io/github/downloads/nwchemgit/nwchem/v7.2.2-release/total)
 ![Github Downloads Release 7.2.3](https://img.shields.io/github/downloads/nwchemgit/nwchem/v7.2.3-release/total)
+![Github Downloads Release 7.3.0](https://img.shields.io/github/downloads/nwchemgit/nwchem/v7.3.0-release/total)
+![Github Downloads Release 7.3.1](https://img.shields.io/github/downloads/nwchemgit/nwchem/v7.3.1-release/total)
 
 
 [![All Releases tags](https://img.shields.io/github/release/nwchemgit/nwchem/all.svg)](https://github.com/nwchemgit/nwchem/releases)
-![GitHub Release Date](https://img.shields.io/github/release-date/nwchemgit/nwchem)
-[![EMSL](https://rawgit.com/nwchemgit/nwchem/master/contrib/git.nwchem/emsl_logo2.svg)](https://www.emsl.pnl.gov)
+![GitHub Release Date](https://img.shields.io/github/release-date/nwchemgit/nwchem)  
+[![EMSL](https://raw.githubusercontent.com/nwchemgit/nwchem/179015c5af4f2b1b397a7f039b85a9aa202f4139/contrib/git.nwchem/emsl_logo2.svg)](https://www.emsl.pnl.gov)
  
 ## NWChem: Open Source High-Performance Computational Chemistry
 <img alt="NWChem logo" src="https://raw.githubusercontent.com/nwchemgit/nwchem/master/contrib/git.nwchem/MS3_logo_cropped.png" align=middle width="102pt" height="117pt"/>  
@@ -32,19 +34,19 @@ Source download instructions for developers
 ```
 git clone https://github.com/nwchemgit/nwchem.git
 ```
-* 7.2.3 Release
+* 7.3.1 Release
 ```
-git clone  -b hotfix/release-7-2-0 https://github.com/nwchemgit/nwchem.git nwchem-7.2.3
+git clone  -b hotfix/release-7-3-0 https://github.com/nwchemgit/nwchem.git nwchem-7.3.1
 ```
-### Release 7.2.3
+### Release 7.3.1
 
 Release files available at the page  
-https://github.com/nwchemgit/nwchem/releases/tag/v7.2.3-release    
+https://github.com/nwchemgit/nwchem/releases/tag/v7.3.0-release    
 Tar files for download:  
-https://github.com/nwchemgit/nwchem/releases/download/v7.2.3-release/nwchem-7.2.3-release.revision-d690e065-srconly.2024-08-27.tar.bz2  
-https://github.com/nwchemgit/nwchem/releases/download/v7.2.3-release/nwchem-7.2.3-release.revision-d690e065-nonsrconly.2024-08-27.tar.bz2  
+https://github.com/nwchemgit/nwchem/releases/download/v7.3.1-release/nwchem-7.3.1-release.revision-23c3b41b-srconly.2025-11-06.tar.bz2  
+https://github.com/nwchemgit/nwchem/releases/download/v7.3.1-release/nwchem-7.3.1-release.revision-23c3b41b-nonsrconly.2025-11-06.tar.bz2  
 Release notes available at  
-https://github.com/nwchemgit/nwchem/blob/master/release.notes.7.2.3.md  
+https://github.com/nwchemgit/nwchem/blob/master/release.notes.7.3.1.md  
 
 ### Documentation
 https://nwchemgit.github.io

contrib/getfiles.nwchem

@@ -1,20 +1,26 @@
 #!/usr/bin/env bash
 #
 # $Id$
+source $NWCHEM_TOP/src/libext/libext_utils/getfiles_utils.sh
 cd $NWCHEM_TOP/src/nwpw/nwpwlib/nwpwxc/
 rm -f dftd3.tgz
 wget https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/dftd3.tgz
+#openblas
 cd $NWCHEM_TOP/src/libext/openblas
-VERSION=0.3.29
-rm -rf OpenBLAS*gz
-curl -L https://github.com/xianyi/OpenBLAS/archive/v${VERSION}.tar.gz -o OpenBLAS-${VERSION}.tar.gz
+rm -rf OpenB* openblas*
+get_openblas
+#scalapack
 cd $NWCHEM_TOP/src/libext/scalapack
-COMMIT=a23c2cdc6586c427686f6097ae66bb54ef693571
-rm -f scalapack-$COMMIT.zip
-curl -L https://github.com/Reference-ScaLAPACK/scalapack/archive/$COMMIT.zip -o scalapack-$COMMIT.zip
+rm -rf scalapa*
+get_scalapack
+#ELPA
+cd $NWCHEM_TOP/src/libext/elpa
+rm -rf elpa*
+get_elpa 
 #caching
 if [ ! -z "$1" ]; then
     rsync -av $NWCHEM_TOP/src/nwpw/nwpwlib/nwpwxc/dftd3.tgz "$1"/.
     rsync -av $NWCHEM_TOP/src/libext/openblas/OpenBLAS*gz "$1"/.
-    rsync -av $NWCHEM_TOP/src/libext/scalapack/scalapack*zip "$1"/.
+    rsync -av $NWCHEM_TOP/src/libext/scalapack/scalapack*tar*gz "$1"/.
+    rsync -av $NWCHEM_TOP/src/libext/elpa/elpa*tar.gz "$1"/.
 fi

contrib/python/nwchem_tddft_spectrum.py

@@ -34,9 +34,9 @@ from pylab import *
 def parse_nwchem_output(fname):
     import re
     from pprint import pprint
-    singlet_pat = re.compile('\s+Root\s+\d+\s+singlet')
-    triplet_pat = re.compile('\s+Root\s+\d+\s+triplet')
-    unrestricted_pat = re.compile('\s+Root\s+\d')
+    singlet_pat = re.compile(r'\s+Root\s+\d+\s+singlet')
+    triplet_pat = re.compile(r'\s+Root\s+\d+\s+triplet')
+    unrestricted_pat = re.compile(r'\s+Root\s+\d')
     singlets = []
     triplets = []
     unrestricted = []
@@ -59,7 +59,7 @@ def parse_nwchem_output(fname):
 
 def parse_element(line,f):
     import re
-    oscillator_pat = re.compile('\s+Oscillator')
+    oscillator_pat = re.compile(r'\s+Oscillator')
     words = line.split()
     try:
         energy = float(words[4])

release.notes.7.3.1.md

@@ -0,0 +1,26 @@
+NWChem Version 7.3.1 Release Notes
+==================================
+
+NWChem is now available on Github at
+https://github.com/nwchemgit/nwchem
+
+User Manual available from the NWChem website
+https://nwchemgit.github.io
+
+NWChem 7.3.1 is released as open-source under the ECL 2.0 license.
+
+NWChem 7.3.1 will be released with the latest Global Arrays Toolkit (v5.9.2).
+
+The  7.3.1 release is a maintenance release containing fixes/enhancements for the NWChem 7.3.0 tree
+
+The change log below is relative to the [7.3.0](release.notes.7.3.0.md) code base.
+
+NEW FUNCTIONALITY
+-----
+
+   N/A  
+
+GITHUB ISSUES ADDRESSED
+----
+   * [water_pme QA test fails: Illegal input to pfft3d](https://github.com/nwchemgit/nwchem/issues/1196)
+   

src/GNUmakefile

@@ -255,7 +255,7 @@ ifdef SUBDIRS
 endif
 
 .PHONY:	install
-install: all
+install: link
 	@test $${INSTALL_PREFIX?Please set the environment variable INSTALL_PREFIX}
 	@echo install_prefix is $(INSTALL_PREFIX)
 	@mkdir -p $(INSTALL_PREFIX)/bin  $(INSTALL_PREFIX)/share

src/NWints/api/exactd_mem.F

@@ -172,7 +172,7 @@ c::local
       integer bra_prim, ket_prim, ket_lim
       integer bra_gen, ket_gen
       integer max_nbf_bra,max_nbf_ket, nat
-      integer bra_fact, ket_fact
+      integer bra_fact, ket_fact, lp
       double precision factor
 c
 #include "bas_exndcf_sfn.fh"
@@ -238,7 +238,11 @@ c
             ket_gen  = 1
           endif
 c
-          memdum = 60 000 000
+c     must be larger than NPP*Lp3*ncenters (line 340 of hf1d.F)
+          lp=bra_type+ket_type+1
+          memdum= (((lp+1)*(lp+2)*(lp+3))/6*
+     *         nat*bra_prim*ket_prim)*125/100
+          memdum=max(memdum,60000000)
           nint = int_unint(brain,ibc,ketin,ikc,0,0,0,0)
           call hf1d(
      &          dbl_mb(mb_exndcf(1,1)), dbl_mb(mb_exndcf(1,1)),

src/NWints/auxil/nmcoef2c.F

@@ -20,6 +20,7 @@ c::local
       integer ip, j, index
       logical warning
       logical oldway
+      double precision factor2,summin
 c
       oldway = .false.
 c
@@ -38,7 +39,6 @@ c... number of integrals
      &       BASIS_ERR)
       nbf = (l+1)*(l+2)/2
       nint = nbf*nbf
-      if (.not.oldway) then
 *
 *...
 *... normalize for overlap of basis function in all space
@@ -68,6 +68,7 @@ c... number of integrals
 *...
 *
         sum = 0.0d00
+        call dfill(nbf,0d0,scr(lscr-nbf+1),1)
         do i = 1,np
           do j = i,np
 *debug:            write(6,*)' '
@@ -81,15 +82,20 @@ c... number of integrals
      &          xyz,z(j),cc(j),1,1,l,
      &          xyz,dumexp,dumcof,1,1,0,
      &          scr,nint,.false.,.false.,.false.,.false.,
-     &          scr(1+nint),(lscr-nint))
+     &          scr(1+nint),(lscr-nint-nbf))
+            if(oldway) then
 * find minimum diagonal element
-            call find_min_diag_nmcoef2c(scr,nbf,factor)
+               call find_min_diag_nmcoef2c(scr,nbf,factor)
 *debug:            write(6,*)' sum 1     ', sum, factor
-            if(i.ne.j) factor = factor*2d0
-            sum = sum + factor
+               if(i.ne.j) factor = factor*2d0
+               sum = sum + factor
 *debug:            write(6,*)' sum 2     ', sum
+            else
+               call sum_diag_nmcoef2c(scr(lscr-nbf+1),scr,nbf,i.eq.j)
+            endif
           enddo
         enddo
+        if(.not.oldway) sum=minval(scr(lscr-nbf+1:lscr))
 *debug:        write(6,*)' sum final ', sum
         factor = 1.0d00/sqrt(sum)
 *debug:        write(6,*)' factor    ',factor
@@ -97,38 +103,6 @@ c... number of integrals
           cc(i) = cc(i)*factor
         enddo
 c
-      else
-        call hf2 (
-     &        xyz,z,cc,np,1,l,
-     &        xyz,dumexp,dumcof,1,1,0,
-     &        xyz,z,cc,np,1,l,
-     &        xyz,dumexp,dumcof,1,1,0,
-     &        scr,nint,.false.,.false.,.false.,.false.,
-     &        scr(1+nint),(lscr-nint))
-*        
-crak:--
-c.. find minimum diagonal element
-c
-        factor = abs(scr(1))
-        index = 0
-        do i = 1,nbf
-          do j = 1,nbf
-            index = index + 1
-            if (i.eq.j) then
-              if (abs(scr(index)).gt.1.0d-05) then
-                factor = min(factor,abs(scr(index)))
-              endif
-*debug:              write(6,*)'factor ',factor,' val ',scr(index)
-            endif
-          enddo
-        enddo
-        factor = 1.0d00/sqrt(factor)
-c
-c... normalize contraction coefficients
-        do i = 1, np
-          cc(i) = cc(i)*factor
-        enddo
-      endif
 *debug:* compute sum of squares
 *debug:      factor = 0.0d00
 *debug:      do ip = 1,np
@@ -189,17 +163,30 @@ c
 *debug:      write(6,*)' nr = ',nr
       factor = abs(buf(1,1))
       negative = buf(1,1).lt.0.0d00
-      do i = 1,nr
-        do j = 1,nr
-          if (i.eq.j) then
-*debug:            write(6,10000)i,j,buf(i,j)
-            if (abs(buf(i,j)).lt.factor) then
-              factor = abs(buf(i,j))
-              negative = buf(i,j) .lt. 0.0d00
-            endif
-          endif
-        enddo
+      do i = 2,nr
+*debug:            write(6,10000)i,buf(i,i)
+         if (abs(buf(i,i)).lt.factor) then
+            factor = abs(buf(i,i))
+            negative = buf(i,i) .lt. 0.0d00
+         endif
       enddo
       if (negative) factor = -1.0d00*factor
-10000 format(' findbuffer(',i3,',',i3,') =',1pd20.10)
+      write(6,*) ' findmin factor ',factor
+10000 format(' findbuffer(',i3,') =',1pd20.10)
+      end
+      subroutine sum_diag_nmcoef2c(test,scr,nbf,ieqj)
+      implicit none
+      integer nbf
+      logical ieqj
+      double precision test(*),scr(nbf,*)
+c
+      integer k
+      double precision fact
+c
+      fact=1d0
+      if(.not.ieqj) fact=2d0
+      do k=1,nbf
+         test(k) = test(k) + fact*scr(k,k)
+      enddo
+      return
       end

src/NWints/auxil/spcartP.fh

@@ -23,7 +23,7 @@
       logical active_sp2c_lindx    ! active ma array sp2c_lindx ?
       logical active_invsp2c_lindx ! active ma array sp2c_lindx ?
       logical active_cart_norm_scale ! active ma array cart_norm_scale ?
-      logical trust_dgemm          ! is dgemm okay
+      logical trust_dgemm          ! is _dgemm_ okay
 *
       parameter (sph_cart_init_value = 1731)
 *

src/NWints/deriv/hf1d.F

@@ -16,29 +16,30 @@ C> Lachlan Belcher 14 Feb 2023
      &       O2I,KEI,NAI,canAB,dryrun,W0,maxW0)
 c $Id$
 
-      Implicit real*8 (a-h,o-z)
-      Implicit integer (i-n)
+      implicit none
       logical hlfdrv    ! [input] if true, calculate same atom as well
                         ! as different atom integrals. This is necessary
                         ! when calculating the half-derivative overlap
                         ! matrix for NADCTs
-      integer ictrA, ictrB !NB already used implicitly
+      integer ictrA, ictrB      !NB already used implicitly
+      integer NPA,NCA,La,NPB,NCB,Lb,ncenters,nint,maxw0
 
       Logical dryrun
       Logical O2I,KEI,NAI,canAB,GenCon
 
 c--> Cartesian Coordinates, Primitives & Contraction Coefficients
 
-      Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
-      Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
+      double precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
+      double precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
 
 c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents
 
-      Dimension Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
+      double precision Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
 
 c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals
 
-      Dimension bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*)
+      Double precision bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*)
+      Double Precision W0(maxW0)
 
 
       hlfdrv=.false.
@@ -66,9 +67,7 @@ C-------------------------------------------------------------------------
      &       hlfdrv)
 c $Id$
 
-      Implicit real*8 (a-h,o-z)
-      Implicit integer (i-n)
-
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
       LOGICAL smatm     ! if true, both integrals use the same atom
@@ -80,30 +79,31 @@ c $Id$
                         ! when calculating the half-derivative overlap
                         ! matrix for NADCTs
       integer ictrA, ictrB !NB already used implicitly
+      integer NPA,NCA,La,NPB,NCB,Lb,ncenters,nint,maxw0
 
       Logical dryrun
       Logical O2I,KEI,NAI,canAB,GenCon
 
 c--> Cartesian Coordinates, Primitives & Contraction Coefficients
 
-      Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
-      Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
+      Double Precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
+      Double Precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
 
 c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents
 
-      Dimension Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
+      Double Precision Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
 
 c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals
 
-      Dimension bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*)
+      Double Precision bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*)
 
 c--> Derivative Indices
 
-      Dimension inder1(3,3),inder2(6,6)
+      integer inder1(3,3),inder2(6,6)
 
 c--> Scratch Space.
 
-      Dimension W0(maxW0)
+      Double Precision W0(maxW0)
 
       Data inder1/ 1,0,0,
      &             0,1,0,
@@ -115,6 +115,12 @@ c--> Scratch Space.
      &             0,0,0,1,0,0,
      &             0,0,0,0,1,0,
      &             0,0,0,0,0,1 /
+c
+      integer i_ALPHAp,i_Ep,i_ESp,i_ff,i_IJK,i_IPAIRp,
+     I     i_left,i_P,i_PC,i_pf,i_R0,i_R0C,i_right,
+     I     i_Rj,i_RS,i_Ti,i_top,i_VA,i_VB,i_VP,i_VR,
+     I     Li, Lp,Lp3, lprod, MaxMem, MXD,
+     I     n, NCP, nd, nder, npp
 c
 c Compute gradient of the overlap, kinetic energy, and nuclear attraction 
 c integrals for 2 shells of contracted Gaussians functions.

src/NWints/dk/GNUmakefile

@@ -4,7 +4,7 @@
 include ../../config/makefile.h
      OBJ_OPTIMIZE = int_1edk_ga.o \
                     dk_input.o dkdiv_energy.o int_pxvp_ga.o \
-                    int_1edkso_ga.o
+                    int_1edkso_ga.o dk_utils.o
 
 ifdef XLF10
     OBJ = bas_create_uncont.o
@@ -16,5 +16,10 @@ endif
 
 USES_BLAS = int_1edk_ga.F  int_1edkso_ga.F
 
+_PDSYEVR := $(shell [ -e ${NWCHEM_TOP}/src/tools/install/lib/libga.a ] &&  nm ${NWCHEM_TOP}/src/tools/install/lib/libga.a|egrep T |egrep  ga_pdsyevr |awk '/ga_pdsyevr/ {print "Y"}')
+ifeq ($(_PDSYEVR),Y)
+ LIB_DEFINES += -DUSE_PDSYEVR
+endif
+
 include ../../config/makelib.h
 

src/NWints/dk/dk_utils.F

@@ -0,0 +1,23 @@
+      subroutine ga_serialcopy(g_inpg,g_inworld)
+      implicit none
+#include "mafdecls.fh"
+#include "errquit.fh"
+#include "global.fh"
+      integer g_inpg,g_inworld
+c
+      integer itype,dim
+      integer k_a,l_a,i,j
+c
+      if(ga_nodeid().eq.0) then
+         call ga_inquire(g_inpg, itype, dim, dim)
+         if (.not. ma_push_get(MT_DBL,dim*dim,'MA temp',l_a,k_a))
+     $        call errquit('serialc: ma_push_get failed MA temp',
+     Z        0       ,     MA_ERR)
+         call ga_get(g_inpg,1,dim,1,dim,dbl_mb(k_a),dim)
+         call ga_put(g_inworld,1,dim,1,dim,dbl_mb(k_a),dim)
+         if (.not. ma_pop_stack(l_a))
+     $        call errquit('serialc: ma_pop_stack l_a',0, MA_ERR)
+      endif
+      call ga_sync()
+      return
+      end

src/NWints/dk/int_1edk_ga.F

@@ -39,6 +39,8 @@ c:tex-\noindent Author: W. A. de Jong
 c:tex-
 c:tex-{\it Syntax:}
 c:tex-\begin{verbatim}
+#define PGROUP
+#define PGROUP_MAXPROC 2
       subroutine int_1edk_ga(basis_in, g_in, type, kskel)
 c:tex-\end{verbatim}
       implicit none
@@ -122,9 +124,13 @@ c
       integer ZpxVp_ga          ! Unmodified pxVp for DK3
       integer tmp3g_a           ! Temporary g_a matrix #3 handle for
                                 ! full DK3
+      integer max2,list2(PGROUP_MAXPROC),p_handle,g_inpg,g_inworld
+      integer comm_world
+      logical status
       integer info
       integer ga_cholesky,ga_llt_i
       external ga_cholesky,ga_llt_i
+      integer itype,dim
 c
 c     cau = 137.0360000d0         ! Speed of light by Bernd Hess
 c
@@ -135,6 +141,25 @@ c
 c
 c     Define the dimensions of the arrays
 c
+      call ga_sync()
+#ifdef PGROUP
+#if defined(SCALAPACK)
+      call slexit
+      call slexit2
+#endif
+      max2=min(ga_nnodes(),PGROUP_MAXPROC)
+      do iloop=1,PGROUP_MAXPROC
+         list2(iloop)=iloop-1
+      enddo
+      p_handle=ga_pgroup_create(list2,max2)
+      call ga_inquire(g_in, itype, dim, dim)
+      comm_world=ga_pgroup_get_world()
+      if (.not. ga_create(MT_DBL,dim,dim,'GA tmp',1,1,g_inworld))
+     $     call errquit('int_1edk_ga: ga_create failed GA temp',0,
+     &     GA_ERR)
+      if(ga_nodeid().lt.max2) then
+      call ga_pgroup_set_default(p_handle)
+#endif
       if (.not. bas_numbf(dk_fitbas,nbf))
      $   call errquit('int_1edk_ga: bas_numbf failed',0, BASIS_ERR)
 c
@@ -182,7 +207,9 @@ c
       if (.not. ma_push_get(MT_DBL,nbf,'psquared',p2m_a,p2Im_a))
      $     call errquit('int_1edk_ga: ma_push_get failed psquared',0,
      &       MA_ERR)
-#ifdef SCALAPACK
+#if defined(SCALAPACK)
+!      call slexit
+!      call slexit2
 c     pdsygv is unstable since it uses the inaccurate pdsyevx
 cold         call ga_pdsygv(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a))
 c     S=chol*chol_t
@@ -190,7 +217,7 @@ c     S=chol*chol_t
 c     inv(a)
       info = ga_llt_i('L',tempg_a,-1)
 c     U=chol(inv(a))_t
-      info = ga_cholesky('L',tempg_a,-1)
+      info = ga_cholesky('L',tempg_a)
 
 c     f_prime=U_t*f*inv(chol)*U
       call ga_dgemm('t','n',nbf,nbf,nbf,1.0d0,tempg_a,Tg_a,0.0d0,tmpg_a)
@@ -201,7 +228,9 @@ c     diag f_prime c_prime
 c     c=inv(chol)_t*c_prime
       call ga_dgemm('n','n',nbf,nbf,nbf,1.0d0,tempg_a,tmpg_a,
      Z     0.0d0,ug_a)
-#elif defined(PEIGS)
+      call slexit
+      call slexit2
+#elif defined(PEIGS) && !defined(PGROUP)
       call ga_diag(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a))
 #else
       call ga_diag_seq(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a))
@@ -1744,6 +1773,7 @@ c
 c     S_-1 is inverse overlap in the uncontracted basis
 c     Scontr_uncon is the mixed basis overlap
 c
+
       if (.not. bas_numbf(basis_in,nbf_in))
      $   call errquit('int_1edk_ga: bas_numbf failed',0, BASIS_ERR)
       if (.not. ga_create(MT_DBL,nbf,nbf_in,'Mixed overlap',1,1,
@@ -1765,7 +1795,11 @@ c     Create overlap matrix dk_fitbas, invert to S^-1 and multiply with Mixg_a
 c
       call ga_zero(tempg_a)
       call int_1e_oldga(dk_fitbas,dk_fitbas,tempg_a,'overlap',.false.)
+#ifdef PGROUP
+      call ga_matpow_serial(tempg_a, -1.0d0, 0.0d0)
+#else
       call ga_matpow(tempg_a, -1.0d0, 0.0d0)
+#endif
       call ga_dgemm('n','n',nbf,nbf_in,nbf,1.0d0,tempg_a,Mixg_a,
      $              0.0d0,tmpg_a)
       call ga_copy(tmpg_a,Mixg_a)
@@ -1782,9 +1816,18 @@ c
       call ga_zero(tempg_a)
       call ga_dgemm('t','n',nbf_in,nbf,nbf,1.0d0,Mixg_a,
      $              DKg_a,0.0d0,tempg_a)
+      call ga_zero(DKg_a)
+#ifdef PGROUP
+      if (.not. ga_create(MT_DBL,dim,dim,'GA tmp',1,1,g_inpg))
+     $     call errquit('int_1edk_ga: ga_create failed GA temp',0,
+     &     GA_ERR)
+      call ga_zero(g_inpg)
+#else
+      g_inpg=g_in
+#endif
       call ga_dgemm('n','n',nbf_in,nbf_in,nbf,1.0d0,tempg_a,Mixg_a,
-     $              1.0d0,g_in)
-      call ga_symmetrize(g_in)
+     $              1.0d0,g_inpg)
+      call ga_symmetrize(g_inpg)
 c
 c     Destroy transformation matrices Mixg_a, tempg_a and DKg_a
 c
@@ -1795,6 +1838,22 @@ c
       if (.not. ga_destroy(DKg_a))
      $     call errquit('int_1edk_ga: ga_destroy DKg_a',0, GA_ERR)
 c
+#ifdef PGROUP
+            call ga_copy(g_inpg,g_inworld)
+      endif
+      call ga_pgroup_set_default(comm_world)
+      call ga_sync()
+c      call ga_serialcopy(g_inpg,g_inworld)
+      if(ga_nodeid().lt.max2) then
+         if (.not. ga_destroy(g_inpg))
+     $        call errquit('int_1edk_ga: ga_destroy g_inpg',0, GA_ERR)
+      endif
+      call ga_add(1d0,g_inworld,1d0,g_in, g_in)
+      if (.not. ga_destroy(g_inworld))
+     $     call errquit('int_1edk_ga: ga_destroy g_inworld',0, GA_ERR)
+      if(.not.ga_pgroup_destroy(p_handle))
+     $     call errquit('int_1edk_ga: pgroupdestroy',0, GA_ERR)
+#endif
       return
 c
       end

src/NWints/simint/libsimint_source/build_simint.sh

@@ -127,7 +127,8 @@ fi
 CMAKE_VER_MAJ=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f1)
 CMAKE_VER_MIN=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f2)
 echo CMAKE_VER is ${CMAKE_VER_MAJ} ${CMAKE_VER_MIN}
-if ((CMAKE_VER_MAJ < 3)) || (((CMAKE_VER_MAJ == 3) && (CMAKE_VER_MIN < 24))); then
+echo CMAKE_VERSION_REQUIRED is ${CMAKE_VER_REQ_MAJ} ${CMAKE_VER_REQ_MIN}
+if ((CMAKE_VER_MAJ < ${CMAKE_VER_REQ_MAJ})) || (((CMAKE_VER_MAJ == ${CMAKE_VER_REQ_MAJ}) && (CMAKE_VER_MIN < ${CMAKE_VER_REQ_MIN}))); then
     cmake_instdir=../../../libext/libext_utils/
     get_cmake_release $cmake_instdir
     echo PATH is $PATH
@@ -222,7 +223,7 @@ if [[ ${FC_EXTRA} == gfortran  || ${FC_EXTRA} == flang || ${FC_EXTRA} == armflan
 fi
 elif  [ ${FC} == xlf ] || [ ${FC} == xlf_r ] || [ ${FC} == xlf90 ]|| [ ${FC} == xlf90_r ]; then
     Fortran_FLAGS=" -qintsize=8 -qextname -qpreprocess"
-elif  [[ ${FC} == ifort || (${FC} == ftn && ${PE_ENV} == INTEL) ]]; then
+elif  [[ ${FC} == ifort || ${FC} == ifx || (${FC} == ftn && ${PE_ENV} == INTEL) ]]; then
     Fortran_FLAGS="-i8 -fpp"
 elif  [ ${FC} == ftn ]  && [ ${PE_ENV} == CRAY  ]; then
     Fortran_FLAGS=" -ffree -s integer64 -e F "

src/NWints/simint/module/GNUmakefile

@@ -11,7 +11,6 @@ ifdef SIMINT_GRADIENT
       LIB_DEFINES +=  -DSIMINT_GRADIENT
 endif
 
-        
        HEADERS = 
 
        LIBRARY = libnwints.a

src/NWints/simint/module/nwcsim_module.F

@@ -7,11 +7,12 @@ c      use ISO_FORTRAN_ENV
 #define SIMINT_SCREEN_SCHWARZ      1
 #define SIMINT_SCREEN_FASTSCHWARZ  2
 C     shells 
-      type(c_simint_shell), target :: smnt_sh(nw_max_shells,3)
+      integer, parameter :: max_bases=4
+      type(c_simint_shell), target :: smnt_sh(nw_max_shells,max_bases)
 c dummy shell for 3c- 
-      type(c_simint_shell), target :: zero_sh(nw_max_shells,3)
+      type(c_simint_shell), target :: zero_sh(nw_max_shells,max_bases)
       type(c_simint_multi_shellpair), target :: bra_msh, ket_msh
-      integer :: nwcsim_noshell(3), nwcsim_bas(3)
+      integer :: nwcsim_noshell(max_bases), nwcsim_bas(max_bases)
       integer :: nwcsim_nbas
       logical :: nwcsim_initialized = .false.
       double precision :: smnt_screen_tol 

src/NWints/simint/source/nwcsim_facef90.F

@@ -15,18 +15,18 @@
 #include "rtdb.fh"
 #include "stdio.fh"
 #include "nwcsim.fh"
-c      use iso_c_binding
-      integer rtdb,nbas,bases(3)
+c     use iso_c_binding
+      integer rtdb,nbas,bases(max_bases)
       integer :: num_der
 c      type(c_simint_multi_shellpair), target :: msh
-      double precision :: alpha(3), coef(3)
       integer :: ibasis,basis,bas
       logical cart_2_sphe
       integer :: geom, natoms
       integer :: max_der, max_ang
       integer :: mxmx
-      integer :: iat,icset,ic1,ic2, nprim
-      integer :: l, ncontr, isphere
+      integer :: iat,icset,ic1,ic2
+      integer :: nprim, l
+      integer :: ncontr, isphere
       character*16 tag 
       integer l_coef,k_coef,l_exp,k_exp
       double precision coord(3)
@@ -47,9 +47,9 @@ c     fetch stuff from rtdb
 c reset screen_tol when screening is off
       if(smnt_screen_method.eq.SIMINT_SCREEN_NONE)smnt_screen_tol=0d0
 
-c     stick to nbas>=2 for now
-         if(nbas.gt.2) call errquit(
-     C        ' simint interface not ready for no basis gt 2',0,0)
+         if(nbas.gt.max_bases) call errquit(
+     C     ' simint interface not tested yet for no basis gt ',
+     M     max_bases, 0)
 c     init
       max_der=0
 #ifdef SIMINT_GRADIENT
@@ -168,7 +168,7 @@ c      iszb_2e4c=isz_2e4c
       implicit none
 #include "errquit.fh"
       integer nbas     ! [in]
-      integer bases(3) ! [in]
+      integer bases(max_bases) ! [in]
 c
       integer icsh,ibasis,bas
       character*14 pname

src/basis/basP.fh

@@ -234,6 +234,7 @@ c.. incore structures
       integer bas_nassoc
       logical bsactive 
       logical bas_spherical
+      logical bas_upcase
       logical bas_any_gc
       logical bas_any_sp_shell
 c
@@ -259,6 +260,7 @@ c......... above for each basis is stored ... below is in-core
      &       bas_spherical(nbasis_bsmx),
      &       bas_any_gc(nbasis_bsmx),
      &       bas_any_sp_shell(nbasis_bsmx),
+     L	     bas_upcase(nbasis_bsmx),
      &       bas_norm_id(nbasis_bsmx),
      &       angular_bs(nbasis_bsmx),
      &       nbfmax_bs(nbasis_bsmx),

src/basis/bas_input.F

@@ -12,7 +12,7 @@ c $Id$
 c
 c   basis [<name>] [library [<standard set>] [file <filename>] \
 c         [spherical|cartesian] [segment||nosegment] [print|noprint]\
-c         [ecpset ecp_name] [soset so_name] [rel] [bse]
+c         [ecpset ecp_name] [soset so_name] [rel] [bse] [uppercase]
 c
 c     tag library [tag in library] <standard set> [file <filename>]
 c     tag <shell type>
@@ -53,23 +53,24 @@ c::passed
       integer rtdb              ! [input] handle to database
 c::local
       integer nopt
-      parameter (nopt = 13)
+      parameter (nopt = 14)
       character*10 opts(nopt)
       character*255 test, name, filename, standard
       character*255 ecpname, soname
       character*256 mydir
       character*1000 errmsg
-      logical status, ospherical, osegment, oprint
+      logical status, ospherical, osegment, oprint, oupcase
       logical o_add_ecpname, o_add_soname
       logical oHas_Star
       logical oIs_ecp, oIs_so, oIs_rel
       integer ind, basis
-      logical bas_add_ucnt_init, bas_set_spherical, bas_add_ucnt_tidy
-      external bas_add_ucnt_init, bas_set_spherical
+      logical bas_add_ucnt_init, bas_set_spherical, bas_add_ucnt_tidy,
+     B     bas_set_upcase
+      external bas_add_ucnt_init, bas_set_spherical,bas_set_upcase
       data opts /
      &    'spherical', 'cartesian', 'segment', 'nosegment', 'library',
      &    'file',      'print',     'noprint', 'ecpset',    'soset',
-     &    'version', 'rel','bse'/
+     &    'version', 'rel','bse', 'uppercase'/
 c
       ecpname = '                                                '
       soname = '                                                '
@@ -105,6 +106,7 @@ c
       o_add_soname = .false.  ! Default is to not associate a so name
       oIs_rel = .false.    ! Default is nonrelativistic
       oHas_Star = .false.  ! Default, not using star-tags 
+      oupcase = .false. ! Default is do not upper case ouput
 c
  10   if (inp_a(test)) then
 c
@@ -122,7 +124,7 @@ c
          endif
 c
          goto (100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
-     &         1100,1200,1300) ind
+     &         1100,1200,1300,1400) ind
          goto 10000
 c
 c     spherical
@@ -212,6 +214,13 @@ c
 01300    continue
          mydir='libraries.bse'
          goto 10
+c
+c
+c uppercase bs_tags
+c
+01400    continue
+         oupcase=.true.
+         goto 10
 c
       endif
 c
@@ -270,6 +279,9 @@ c
      &       BASIS_ERR)
       if (.not. bas_set_spherical(basis, ospherical))
      &      call errquit
+     &      (' bas_set_spherical failed ',911, BASIS_ERR)
+      if (.not. bas_set_upcase(basis, oupcase))
+     &      call errquit
      &      (' bas_set_spherical failed ',911, BASIS_ERR)
       
       call bas_input_body(basis, osegment, oIs_rel, oHas_Star)

src/basis/basis.F

@@ -330,6 +330,7 @@ c
       bas_spherical(basis)    = .false.
       bas_any_gc(basis)       = .false.
       bas_any_sp_shell(basis) = .false.
+      bas_upcase              = .false.
 c
       h_tmp = exndcf(H_exndcf,basis)
       if (h_tmp .ne. -1) then
@@ -549,6 +550,8 @@ c::-local
           endif
         enddo
         dum_tag = bs_tags(i_tag,basis)
+        if (bas_upcase(basis)) call inp_ucase(dum_tag)
+        
         tmp1 = inp_strlen(bs_stdname(i_tag,basis))
         if (tmp1 .lt. (30-1)) then
           tmp2 = (30-tmp1)/2
@@ -663,6 +666,7 @@ c:: local
 *. . . . . . . . . . . ! Room for tag+space+(+element+) = 16+1+1+16+1
       character*35 buffer  
       character*12 polynomial
+      character*16 case_tag
 c
 #include "bas_exndcf.fh"
 #include "ecpso_sfnP.fh"
@@ -723,8 +727,11 @@ c
      $        atn)) then
             len_tag = inp_strlen(bs_tags(i,basis))
             len_ele = inp_strlen(element)
+            case_tag=bs_tags(i,basis)
+            if (bas_upcase(basis)) call inp_ucase(case_tag)
+
             write(buffer,'(a,'' ('',a,'')'')')
-     $           bs_tags(i,basis)(1:len_tag), element(1:len_ele)
+     $           case_tag(1:len_tag), element(1:len_ele)
          else
             buffer = bs_tags(i,basis)
          endif

src/basis/basisP.F

@@ -513,3 +513,24 @@ c
         endif
       endif
       end
+      logical function bas_set_upcase(basisin, oupcase)
+      implicit none
+#include "nwc_const.fh"
+#include "basP.fh"
+#include "basdeclsP.fh"
+c
+      integer basisin    ! [input] basis set handle
+      logical oupcase ! [input] logical for upcase setting
+c
+      integer basis
+c
+      bas_set_upcase = .false.
+c
+      basis = basisin + BASIS_HANDLE_OFFSET
+      if (oupcase) then
+        bas_upcase(basis) = .true.
+      else
+        bas_upcase(basis) = .false.
+      endif
+      bas_set_upcase = .true.
+      end

src/basis/libraries/qmmm_zhang_3-21g_ecp

@@ -1,3 +1,8 @@
+# QMMM Zhang 3-21G ecp
+# Elements                             References
+# --------                             ----------
+# F      :  Y. Zhang, T.-S. Lee, W. Yang, J. Chem. Phys.  110, 46–54 (1999)
+# 
 ecp "F_qmmm zhang 3-21g ecp"
   F nelec  2 # down to carbon
   F ul

src/basis/libraries/qmmm_zhang_6-31gs_ecp

@@ -1,3 +1,8 @@
+# QMMM Zhang 6-31G* ecp
+# Elements                             References
+# --------                             ----------
+# F      :  Y. Zhang, T.-S. Lee, W. Yang, J. Chem. Phys.  110, 46–54 (1999)
+# 
 ecp "F_qmmm zhang 6-31G* ecp"
   F nelec  2 # down to carbon
   F ul

src/bq/bq_data.F

@@ -664,7 +664,7 @@ c      local variables
      &      0, MA_ERR)
 
        if (oprint) then
-         call util_print_centered(6,
+         call util_print_centered(luout,
      >      "Bq Structure Information (Angstroms)",
      >              36, .true.)
       
@@ -674,7 +674,7 @@ c      local variables
 c        == tally up bq charges ==
          bq_charge_total = 0.d0
          do i=1,bq_ncent(handle)
-           write(6,FMT=9000)
+           write(luout,FMT=9000)
 c     >           i,(dbl_mb(i_c+3*(i-1)+k-1),k=1,3),
      >           i,(dbl_mb(i_c+3*(i-1)+k-1)*0.529177249d00,k=1,3),
      >           dbl_mb(i_q+i-1)

src/config/32_to_64

@@ -10,19 +10,8 @@ if [[ -z "${NWCHEM_TOP}" ]]; then
     MYPWD=`pwd`
     NWCHEM_TOP=`echo ${MYPWD}/${DIRUTIL} | sed -e 's/\/src.*//' `
 fi
-  if [ -x /msrc/apps/bin/perl ]; then
-      perlexe=/msrc/apps/bin/perl
-  elif [ -x /usr/local/bin/perl ]; then
-      perlexe=/usr/local/bin/perl
-  elif [ -x /usr/local/gnu/bin/perl ]; then
-      perlexe=/usr/local/gnu/bin/perl
-  elif [ -x /usr/gnu/bin/perl ]; then
-      perlexe=/usr/gnu/bin/perl
-  else
-# assume perl is in your path
-      perlexe=perl
-  fi
-   perlscript=${NWCHEM_TOP}/src/config/32_to_64.pl
+pythonexe=python
+pythonscript=${NWCHEM_TOP}/src/config/32_to_64.py
 
 
  JOB_LIMIT=6
@@ -31,16 +20,18 @@ fi
 
 for file in "$@"
 do
-
+if [ -f "$file" ]; then
    echo converting "$file" to 64-bit integers for BLAS/LAPACK
-      ($perlexe $perlscript $file) &	
+      ($pythonexe $pythonscript $file) &	
 
    let njob++
    if [ $njob -gt $JOB_LIMIT ]; then
      wait
      let njob=0
    fi
-
+else
+    echo "!!!!file missing: " "$file"
+fi
 done
 
 if [ $njob -gt  0 ] ; then

src/config/32_to_64.py

@@ -0,0 +1,218 @@
+#!/usr/bin/env python3
+"""
+Python script to do transliteration from "single" values to "double" values
+Usage: python 32_to_64.py file1.f [file2.f ...]
+Original Perl script:
+    Written:  3/14/97
+    By:       Ricky A. Kendall
+    High Performance Computational Chemistry Group
+    Theory Modeling and Simulation Program [2]
+
+Converted from Perl to Python using Claude (claude-4-5-opus-4-5-20251101-v1)
+Optimizations applied:
+    - Early exits to skip unnecessary processing
+    - Reduced redundant regex passes
+    - Combined line-type checks
+    - Precompiled regex patterns
+    - Lookahead for correct underscore boundary matching
+"""
+
+import os
+import sys
+import re
+
+debug = False
+
+
+def copy_case(from_str, to_str):
+    """Take case from 'from_str' and apply it to 'to_str' and return that new string"""
+    result = []
+    for i, char in enumerate(to_str):
+        if i < len(from_str):
+            if from_str[i].isupper():
+                result.append(char.upper())
+            elif from_str[i].islower():
+                result.append(char.lower())
+            else:
+                result.append(char)
+        else:
+            result.append(char)
+    return ''.join(result)
+
+
+class ConversionPattern:
+    """Holds precompiled regex patterns for a single from/to conversion pair"""
+    
+    def __init__(self, from_str, to_str):
+        self.from_str = from_str
+        self.to_str = to_str
+        self.from_str_upper = from_str.upper()
+        # Precompile all regex patterns
+        escaped = re.escape(from_str)
+        
+        # Pattern for finding the from_str (case insensitive)
+        self.search_pattern = re.compile(escaped, re.IGNORECASE)
+        
+        # Pattern for Fortran continuation lines (5 spaces + non-space)
+        self.fortran_pattern = re.compile(
+            r'([ ]{5}.)' + escaped + r'(?=_|\W)',
+            re.IGNORECASE
+        )
+
+        # Pattern for general substitution
+        # Use capturing group for preceding \W (consumed, restored via group(1))
+        # Use lookahead (?=_|\W) so underscore after pattern is not consumed [2]
+        # This correctly handles: void ycopy_() -> void dcopy_()
+        self.general_pattern = re.compile(
+            r'(\W)' + escaped + r'(?=_|\W)',
+            re.IGNORECASE
+        )
+
+
+def load_data_file():
+    """Load the from/to conversion pairs from the data file and precompile patterns"""
+    patterns = []
+    
+    data_path = os.path.dirname(os.path.abspath(__file__))
+    data_path = os.path.join(data_path, "data.64_to_32")
+    
+    if debug:
+        print(f"Data path: {data_path}")
+    
+    try:
+        with open(data_path, 'r') as data_file:
+            for line in data_file:
+                if line and not line.startswith('#'):
+                    tokens = line.split()
+                    if len(tokens) >= 2:
+                        # For 32_to_64: tokens[1] is "from", tokens[0] is "to" [2]
+                        patterns.append(ConversionPattern(tokens[1], tokens[0]))
+    except IOError:
+        sys.exit(f"Unable to open: {data_path}")
+    
+    if debug:
+        print(f"Loaded {len(patterns)} conversion patterns")
+    
+    return patterns
+
+
+def process_line(line, patterns, line_upper):
+    """
+    Process a single line with all conversion patterns.
+    Uses early exits and optimized regex passes [2]
+    """
+    # Early exit: Skip comment lines (Fortran style) or empty lines [2]
+    if not line or line[0] in ('c', 'C', '*') or line.strip() == '':
+        return line
+
+    # Early exit: Quick check if any pattern might match using uppercase comparison
+    has_potential_match = any(conv.from_str_upper in line_upper for conv in patterns)
+    if not has_potential_match:
+        return line
+
+    # Determine line type once (not for each pattern) [2]
+    is_fortran_continuation = (
+        line.startswith('     ') and len(line) > 5 and not line[5].isspace()
+    )
+    # Matches lines starting with space or tab [2]
+    is_general_line = len(line) > 0 and line[0] in ' \t'
+    # Matches any line starting with space or non-whitespace char [2]
+    # ^[ \S] in Perl matches virtually every non-empty line
+    # including C declarations like "void ycopy_()"
+    is_declaration_line = len(line) > 0 and (line[0] == ' ' or not line[0].isspace())
+
+    # Process each conversion pattern
+    for conv in patterns:
+        # Early exit: Skip if pattern not in line (case-insensitive quick check)
+        if conv.from_str_upper not in line_upper:
+            continue
+
+        match = conv.search_pattern.search(line)
+        if not match:
+            continue
+
+        # Found a match - compute replacement
+        froom = match.start()
+        toot = copy_case(
+            line[froom:froom + len(conv.from_str)],
+            conv.to_str
+        )
+
+        # Apply Fortran continuation substitution [2]
+        if is_fortran_continuation:
+            line = conv.fortran_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Apply general substitution for tab/space lines [2]
+        if is_general_line:
+            line = conv.general_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Apply declaration substitution for C-style lines [2]
+        # Handles cases like "void ycopy_()" where line starts with non-space
+        if is_declaration_line:
+            line = conv.general_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Update line_upper after substitution for subsequent patterns
+        line_upper = line.upper()
+
+    return line
+
+
+def process_file(file_path, patterns):
+    """Process a single file and perform transliteration"""
+    pid = os.getpid()
+    orgfile = file_path
+    filebak = f"{file_path}.{pid}"
+    
+    if debug:
+        print(f"Processing: {file_path}")
+        print(f"Backup file: {filebak}")
+    
+    # Rename original to backup [2]
+    os.rename(file_path, filebak)
+    
+    try:
+        with open(filebak, 'r') as fh_in:
+            with open(orgfile, 'w') as fh_out:
+                for line in fh_in:
+                    line_upper = line.upper()
+                    processed_line = process_line(line, patterns, line_upper)
+                    fh_out.write(processed_line)
+    except IOError as e:
+        sys.exit(f"Can't open file: {e}")
+    
+    # Remove backup [2]
+    os.unlink(filebak)
+
+
+def main():
+    if debug:
+        print(f"Arguments: {sys.argv[1:]}")
+    
+    # Load patterns once (precompiled)
+    patterns = load_data_file()
+    
+    if len(patterns) == 0:
+        sys.exit("Fatal sngl2dbl error: No conversion patterns loaded")
+    
+    files = sys.argv[1:]
+
+    if len(files) == 0:
+        print("Usage: python 32_to_64.py file1.f [file2.f ...]")
+        sys.exit(1)
+
+    # Process each file sequentially
+    for file_path in files:
+        process_file(file_path, patterns)
+
+
+if __name__ == "__main__":
+    main()

src/config/64_to_32

@@ -10,39 +10,35 @@ if [[ -z "${NWCHEM_TOP}" ]]; then
     MYPWD=`pwd`
     NWCHEM_TOP=`echo ${MYPWD}/${DIRUTIL} | sed -e 's/\/src.*//' `
 fi
-  if [ -x /msrc/apps/bin/perl ]; then
-      perlexe=/msrc/apps/bin/perl
-  elif [ -x /usr/local/bin/perl ]; then
-      perlexe=/usr/local/bin/perl
-  elif [ -x /usr/local/gnu/bin/perl ]; then
-      perlexe=/usr/local/gnu/bin/perl
-  elif [ -x /usr/gnu/bin/perl ]; then
-      perlexe=/usr/gnu/bin/perl
-  else
-# assume perl is in your path
-      perlexe=perl
-  fi
-   perlscript=${NWCHEM_TOP}/src/config/64_to_32.pl
-
+pythonexe=python
+pythonscript=${NWCHEM_TOP}/src/config/64_to_32.py
+filecheck="/tmp/myfail.txt"
 
+rm -f $filecheck
  JOB_LIMIT=6
 
  njob=0
-
+ nfail=0
 for file in "$@"
 do
 if [ -f "$file" ]; then
    echo converting "$file" to 32-bit integers for BLAS/LAPACK
-      ($perlexe $perlscript $file) &	
+   ($pythonexe $pythonscript $file ; if [ $? == 123 ]; then echo $file> $filecheck;fi ) &
 
    let njob++
    if [ $njob -gt $JOB_LIMIT ]; then
      wait
      let njob=0
    fi
+else
+    echo "!!!!file missing: " "$file"
 fi
 done
-
 if [ $njob -gt  0 ] ; then
     wait
 fi
+if [  -f $filecheck ]; then
+    echo 'Error: DOS file '$(cat $filecheck)
+    exit 1
+fi
+    

src/config/64_to_32.py

@@ -0,0 +1,238 @@
+#!/usr/bin/env python3
+"""
+Python script to do transliteration from "double" values to "8wrap" values
+Usage: python 64_to_32.py file1.f [file2.f ...]
+
+Original Perl script:
+    Written:  3/14/97
+    By:       Ricky A. Kendall
+
+Converted from Perl to Python using Claude (claude-4-5-opus-4-5-20251101-v1)
+Optimizations applied:
+    - Early exits to skip unnecessary processing
+    - Reduced redundant regex passes
+    - Combined line-type checks
+    - Precompiled regex patterns
+    - Lookahead for correct underscore boundary matching
+"""
+
+import os
+import sys
+import re
+import subprocess # nosec B404
+
+debug = False
+
+
+def copy_case(from_str, to_str):
+    """Take case from 'from_str' and apply it to 'to_str' and return that new string"""
+    result = []
+    for i, char in enumerate(to_str):
+        if i < len(from_str):
+            if from_str[i].isupper():
+                result.append(char.upper())
+            elif from_str[i].islower():
+                result.append(char.lower())
+            else:
+                result.append(char)
+        else:
+            result.append(char)
+    return ''.join(result)
+
+
+class ConversionPattern:
+    """Holds precompiled regex patterns for a single from/to conversion pair"""
+    
+    def __init__(self, from_str, to_str):
+        self.from_str = from_str
+        self.to_str = to_str
+        self.from_str_upper = from_str.upper()
+        
+        # Precompile all regex patterns
+        escaped = re.escape(from_str)
+        
+        # Pattern for finding the from_str (case insensitive)
+        self.search_pattern = re.compile(escaped, re.IGNORECASE)
+        
+        # Pattern for Fortran continuation lines (5 spaces + non-space)
+        self.fortran_pattern = re.compile(
+            r'([ ]{5}.)' + escaped + r'(?=_\(|\W)',
+            re.IGNORECASE
+        )
+
+        # Pattern for general substitution
+        # Use capturing group for preceding \W (consumed, restored via group(1))
+        # Lookahead (?=_\(|\W) matches underscore only when followed by ( [1]
+        # This correctly handles:
+        #   DCOPY_()    -> YCOPY_()    (matches - underscore followed by paren)
+        #   DCOPY_OMP() -> DCOPY_OMP() (no match - underscore followed by letter)
+        #   dcopy(      -> ycopy(      (matches - non-word char after)
+        self.general_pattern = re.compile(
+            r'(\W)' + escaped + r'(?=_\(|\W)',
+            re.IGNORECASE
+        )
+
+
+def load_data_file():
+    """Load the from/to conversion pairs from the data file and precompile patterns"""
+    patterns = []
+    
+    data_path = os.path.dirname(os.path.abspath(__file__))
+    data_path = os.path.join(data_path, "data.64_to_32")
+    
+    if debug:
+        print(f"Data path: {data_path}")
+
+    try:
+        with open(data_path, 'r') as data_file:
+            for line in data_file:
+                if line and not line.startswith('#'):
+                    tokens = line.split()
+                    if len(tokens) >= 2:
+                        # For 64_to_32: tokens[0] is "from", tokens[1] is "to" [1]
+                        patterns.append(ConversionPattern(tokens[0], tokens[1]))
+    except IOError:
+        sys.exit(f"Unable to open: {data_path}")
+    
+    if debug:
+        print(f"Loaded {len(patterns)} conversion patterns")
+    
+    return patterns
+
+
+# Precompile line-type detection patterns (used for every line)
+FORTRAN_CONTINUATION = re.compile(r'^[ ]{5}[^\s]')
+
+
+def process_line(line, patterns, line_upper):
+    """
+    Process a single line with all conversion patterns.
+    Uses early exits and optimized regex passes.
+    """
+    # Early exit: Skip comment lines (Fortran style) or empty lines [1]
+    if not line or line[0] in ('c', 'C', '*') or line.strip() == '':
+        return line
+    
+    # Early exit: Quick check if any pattern might match using uppercase comparison
+    # This avoids regex overhead for lines with no potential matches
+    has_potential_match = False
+    for conv in patterns:
+        if conv.from_str_upper in line_upper:
+            has_potential_match = True
+            break
+    
+    if not has_potential_match:
+        return line
+    
+    # Determine line type once (not for each pattern) [1]
+    is_fortran_continuation = line.startswith('     ') and len(line) > 5 and not line[5].isspace()
+    is_general_line = len(line) > 0 and line[0] in ' \t'
+    # Matches any line starting with space or non-whitespace char [1]
+    # ^[ \S] in Perl matches virtually every non-empty line
+    # including C declarations like "void dcopy_()"
+    is_declaration_line = len(line) > 0 and (line[0] == ' ' or not line[0].isspace())
+
+    # Process each conversion pattern
+    for conv in patterns:
+        # Early exit: Skip if pattern not in line (case-insensitive quick check)
+        if conv.from_str_upper not in line_upper:
+            continue
+
+        match = conv.search_pattern.search(line)
+        if not match:
+            continue
+        
+        # Found a match - compute replacement
+        froom = match.start()
+        toot = copy_case(
+            line[froom:froom + len(conv.from_str)],
+            conv.to_str
+        )
+
+        # Apply Fortran continuation substitution [1]
+        if is_fortran_continuation:
+            line = conv.fortran_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Apply general substitution for tab/space lines [1]
+        if is_general_line:
+            line = conv.general_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Apply declaration substitution for C-style lines [1]
+        # Handles cases like "void dcopy_()" where line starts with non-space
+        if is_declaration_line:
+            line = conv.general_pattern.sub(
+                lambda m: m.group(1) + toot,
+                line
+            )
+
+        # Update line_upper after substitution for subsequent patterns
+        line_upper = line.upper()
+    
+    return line
+
+
+def process_file(file_path, patterns):
+    """Process a single file and perform transliteration"""
+    pid = os.getpid()
+    orgfile = file_path
+    filebak = f"{file_path}.{pid}"
+    
+    if debug:
+        print(f"Processing: {file_path}")
+        print(f"Backup file: {filebak}")
+    
+    # Runs the Linux 'file' command
+    
+    result = subprocess.run(['file', file_path], shell=False, capture_output=True, text=True) # nosec
+
+    # DOS files output contains "with CRLF line terminators"
+    if "CRLF" in result.stdout:
+        print(f"ERROR: DOS file with CRLF {file_path}", file=sys.stderr)
+        sys.exit(123)  
+    
+    # Rename original to backup [1]
+    os.rename(file_path, filebak)
+    
+    try:
+        with open(filebak, 'r') as fh_in:
+            with open(orgfile, 'w') as fh_out:
+                for line in fh_in:
+                    line_upper = line.upper()
+                    processed_line = process_line(line, patterns, line_upper)
+                    fh_out.write(processed_line)
+    except IOError as e:
+        sys.exit(f"Can't open file: {e}")
+    
+    # Remove backup [1]
+    os.unlink(filebak)
+
+
+def main():
+    if debug:
+        print(f"Arguments: {sys.argv[1:]}")
+    
+    # Load patterns once (precompiled)
+    patterns = load_data_file()
+    
+    if len(patterns) == 0:
+        sys.exit("Fatal dbl2sngl error: No conversion patterns loaded")
+    
+    files = sys.argv[1:]
+    
+    if len(files) == 0:
+        print("Usage: python 64_to_32.py file1.f [file2.f ...]")
+        sys.exit(1)
+    
+    # Process each file sequentially
+    for file_path in files:
+        process_file(file_path, patterns)
+
+
+if __name__ == "__main__":
+    main()

src/config/makefile.h

@@ -20,7 +20,7 @@
 # For development tree 
 #RELEASE := 
 # For current release tree
-RELEASE := 7.3.0
+RELEASE := 7.3.1
 
 #
 ifndef NWCHEM_TOP
@@ -287,6 +287,14 @@ endif
 # their header files are needed for dependency analysis of
 # other NWChem modules
 ifdef USE_INTERNALBLAS
+    ifneq ($(or $(BLASOPT),$(LAPACK_LIB),$(BLAS_LIB)),)
+        $(info     )
+        $(info You must unset)
+        $(info BLASOPT ,LAPACK_LIB and BLAS_LIB)
+        $(info when using USE_INTERNALBLAS )
+        $(info )
+        $(error )
+    endif
     NW_CORE_SUBDIRS += blas lapack
 endif
 ifdef USE_LIBXC
@@ -305,10 +313,10 @@ ifdef BUILD_OPENBLAS
     NW_CORE_SUBDIRS += libext
 
     #bail out if BLASOPT or LAPACK_LIB or BLAS_LIB are defined by user
-    ifneq ($(or $(BLASOPT),$(LAPACK_LIB),$(BLAS_LIB)),)
+    ifneq ($(or $(BLASOPT),$(LAPACK_LIB),$(BLAS_LIB),$(USE_INTERNALBLAS)),)
         $(info     )
         $(info You must unset)
-        $(info BLASOPT ,LAPACK_LIB and BLAS_LIB)
+        $(info USE_INTERNALBLAS, BLASOPT ,LAPACK_LIB and BLAS_LIB)
         $(info when using BUILD_OPENBLAS )
         $(info )
         $(error )
@@ -376,9 +384,23 @@ ifdef BUILD_ELPA
 #   endif
 
     ifndef SCALAPACK_SIZE
-        SCALAPACK_SIZE=8
+        SCALAPACK_SIZE = 8
+    endif
+    ELPA := -L$(NWCHEM_TOP)/src/libext/lib -lnwc_elpa -I$(NWCHEM_TOP)/src/libext/include/elpa/modules
+    ifdef ELPA_NVIDIA
+       ifdef CUDA_ROOT
+          ELPA := $(ELPA) -L$(CUDA_ROOT)/lib64
+#          ELPA := $(ELPA) $(shell PKG_CONFIG_PATH=$(NWCHEM_TOP)/src/libext/lib/pkgconfig pkg-config --libs-only-l elpa | sed -e 's/-lelpa//g')
+	  ELPA := $(ELPA) -lcusolver -lcudart -lcublasLt -lcublas
+          ELPA := $(ELPA) -lstdc++
+       else
+          $(info )
+          $(info please set the env. variable CUDA_ROOT)
+          $(info when compiling ELPA for Nvidia GPUs)
+          $(info )
+          $(error )
+       endif
     endif
-    ELPA=-L$(NWCHEM_TOP)/src/libext/lib -lnwc_elpa -I$(NWCHEM_TOP)/src/libext/include/elpa/modules
 endif
 
 
@@ -1357,7 +1379,11 @@ ifeq ($(TARGET),MACX64)
 
         ifdef USE_OPENMP
             FOPTIONS  += -fopenmp
-            LDOPTIONS += -fopenmp
+            ifneq (,$(wildcard ${HOMEBREW_PREFIX}/opt/libomp/lib/libomp.a))
+                LDOPTIONS += -L$(HOMEBREW_PREFIX)/opt/libomp/lib -lomp
+            else
+                LDOPTIONS += -fopenmp
+            endif
         endif
 
         ifeq ($(USE_FLANG),1)

src/config/makelib.h

@@ -264,7 +264,7 @@ endif
 endif # QUICKBUILD
 ifdef CONVERT_ALL
 64_to_32:
-	$(CNFDIR)/64_to_32 *.F *.f *.c *.f90
+	$(CNFDIR)/64_to_32 *.F *.f *.c *.f90 *fh
 ifdef SUBDIRS
 	$(MAKESUBDIRS)
 endif

src/ddscf/localize.F

@@ -227,6 +227,7 @@ c
 #include "geom.fh"
 #include "bas.fh"
 #include "util.fh"
+#include "stdio.fh"
 c
 c     Localize the nloc orbitals in iloc() by mixing with each other
 c
@@ -251,10 +252,10 @@ c     locpow (localization sum exponent) should be 2 or 4
       gamma_tol = 1d-10
 c
       if (ga_nodeid() .eq. 0) then
-        write(6,9001)
+        write(luout,9001)
  9001   format(/10x,' iter   Max. delocal   Mean delocal    Converge'/
      &          10x,' ----   ------------   ------------    --------')
-         call util_flush(6)
+         call util_flush(luout)
       end if
 c
       dprev = 0.0d0
@@ -281,7 +282,7 @@ c
                end do
                qs = qs + qas**2
             end do
-*            write(6,*) ' ds ', s, 1.0d0/qs
+*            write(luout,*) ' ds ', s, 1.0d0/qs
             dmax = max(dmax, 1.0d0/qs)
             d = d + 1.0d0/qs
          end do
@@ -298,9 +299,9 @@ c
          call ga_dgop(2, d , 1, '+')
 c
          if (ga_nodeid() .eq. 0) then
-           write(6,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)')
+           write(luout,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)')
      &       iter, dmax, d/dble(nloc), gamma_max
-           call util_flush(6)
+           call util_flush(luout)
          end if
          call ga_sync
 c
@@ -318,7 +319,7 @@ c
                if (tt .le. nloc) then
                   s = iloc(ss)
                   t = iloc(tt)
-*                  write(6,*) nloc, neven, set, pair, ss, tt, s, t
+*                  write(luout,*) nloc, neven, set, pair, ss, tt, s, t
                   call ga_get(g_c,  1, nbf, s, s, c(1,1), 1)
                   call ga_get(g_sc, 1, nbf, s, s,sc(1,1), 1)
                   call ga_get(g_c,  1, nbf, t, t, c(1,2), 1)
@@ -408,6 +409,7 @@ c
 #include "geom.fh"
 #include "bas.fh"
 #include "util.fh"
+#include "stdio.fh"
 c
 c     Localize the nloc orbitals in iloc(*) by mixing with each other
 c
@@ -441,17 +443,17 @@ c
       u21 = util_random(12345)
 c
       if (ga_nodeid() .eq. 0) then
-        write(6,9001)
+        write(luout,9001)
  9001   format(/10x,' iter   Max. dipole2   Mean dipole2    Converge'/
      &          10x,' ----   ------------   ------------    --------')
-        call util_flush(6)
+        call util_flush(luout)
       end if
 c
       dprev = 0.0d0
       gamma_max = 0.0d0
       call ga_distribution(g_uc(1), ga_nodeid(),
      I          ilo, ihi, jlo, jhi)
-      do iter = 1, 100
+      do iter = 1, 900
          call ga_sync
          nrot = 0
 c
@@ -496,9 +498,9 @@ c
          call ga_dgop(2, d , 1, '+')
 c
          if (ga_nodeid() .eq. 0) then
-           write(6,'(10x, i5, 2f17.8, 1p,2d12.2)')
+           write(luout,'(10x, i5, 2f17.8, 1p,2d12.2)')
      &       iter, dmax, d/dble(nloc), gamma_max
-           call util_flush(6)
+           call util_flush(luout)
          end if
          call ga_sync
 c
@@ -521,7 +523,7 @@ cold            do pair = 1+ga_nodeid(), neven/2, ga_nnodes()
 c     if g_uxc.s belongs to ga_nodeid, do
                   if(s.ge.jlo.and.s.le.jhi) then
                   t = iloc(tt)
-*                  write(6,*) nloc, neven, set, pair, ss, tt, s, t
+*                  write(luout,*) nloc, neven, set, pair, ss, tt, s, t
                   call ga_get(g_c,ilo, ihi, s, s, c(1,1), 1)
                   call ga_get(g_c,ilo, ihi, t, t, c(1,2), 1)
                   do x = 1, 4
@@ -547,7 +549,7 @@ c
                      ast = ast + u12*u12 - 0.25d0*(u1-u2)**2
                      bst = bst + u12*(u1 - u2)
                      if (abs(u12-u21)/max(1.0d0,abs(u12)).gt.1d-8) then
-                        write(6,*) ' U12, U21 ', u12, u21
+                        write(luout,*) ' U12, U21 ', u12, u21
                         call errquit('bad u12', 0, UNKNOWN_ERR)
                      endif
                   end do
@@ -584,13 +586,16 @@ c
             call ga_sync
          end do ! do set
       end do ! do iter
+      if(ga_nodeid().eq.0) write(luout,*)
+     W     ' WARNING: Foster-Boys solver not converged'
 c
  1000 continue
 c
 c     Analyze localization of each mo
 c
       if (ga_nodeid() .eq. 0) then
-         write(6,*)
+         write(luout,*)
+         write(luout,*) ' Orbital Localization Analysis by atom'
          do ss = 1, nloc
             s = iloc(ss)
             call ga_get(g_c,  1, nbf, s, s, c(1,1), 1)
@@ -621,10 +626,10 @@ c
                   end if
                end do
             end do
-            write(6,9002) s, (list(a), pop(a), a=1,nlist)
+            write(luout,9002) s, (list(a), pop(a), a=1,nlist)
  9002         format(i5, 100(2x,i4,'(',f5.2,')'))
          end do
-         call util_flush(6)
+         call util_flush(luout)
       end if
 c
       call ga_sync

src/ddscf/rohf_diis.F

@@ -13,6 +13,7 @@ C     $Id$
 #include "util.fh"
 #include "cscf.fh"
 #include "cosmo.fh"
+#include "stdio.fh"
 c     
 c     Solve the ROHF equations using DIIS
 c     
@@ -69,7 +70,7 @@ c
 c
       odbug=.false.
       if(odbug) then
-         write(6,*) 'in rohf_diis_solve ... with node id =',
+         write(luout,*) 'in rohf_diis_solve ... with node id =',
      $              ga_nodeid()
       endif
 c
@@ -81,8 +82,8 @@ c   Print info
 c
       nodezero = ga_nodeid().eq.0
       if (nodezero.and. oprint_parm) then
-         write(6,1) gnorm_tol, maxiter, tol2e, mxdiisbas
-         write(6,111)
+         write(luout,1) gnorm_tol, maxiter, tol2e, mxdiisbas
+         write(luout,111)
  1       format(//,
      $        1x,'----------------------------------------------',/
      $        1x,'               DIIS ROHF',//,
@@ -92,7 +93,7 @@ c
      $        1x,'DIIS subspace             :',9x,i4,/)
  111     format(
      $        1x,'----------------------------------------------',/)
-         call util_flush(6)
+         call util_flush(luout)
       endif
 c     
 c     Allocate global arrays.  
@@ -174,21 +175,22 @@ c
         call rohf_errvec( g_fock, g_movecs, g_over, g_xorth,
      $                    g_tmp, g_tmp2, g_err )
         call ga_maxelt( g_err, diis_emax )
+        call ga_sync()
 c
 c   Print energy
 c
         if(iter.eq.1) then
            if(cosmo_on.and.cosmo_phase.eq.1) then
-              if(nodezero) write(6,909)
+              if(nodezero) write(luout,909)
            else
-              if(nodezero) write(6,910)
+              if(nodezero) write(luout,910)
            endif
         endif
  909    format(2x,'   COSMO gas phase')
  910    format(2x,'   COSMO solvation phase')
 c
         if(nodezero) then
-           write(6,911) iter, ndiis, energy, gnorm, diis_emax
+           write(luout,911) iter, ndiis, energy, gnorm, diis_emax
         endif
  911    format(5x,i5,i3,5x,f20.12,5x,2e10.2)
 c
@@ -216,7 +218,7 @@ c
 c     End SCF minimisation
 c     
  333  continue
-      if (nodezero) call util_flush(6)
+      if (nodezero) call util_flush(luout)
 c     
 c     *ALWAYS* return canonical MOs ... do not do aufbau
 c     
@@ -257,19 +259,19 @@ c
          egas      =energy
          cosmo_phase = 2
          if(odbug) then
-            write(6,*) 'gas_phase done, do sol-phase now ...'
+            write(luout,*) 'gas_phase done, do sol-phase now ...'
          endif
          call ga_sync()
          go to 1000
       elseif(cosmo_on.and.cosmo_phase.eq.2) then
          esol      =energy
          if(nodezero) then
-            write(6,*) '                 COSMO solvation results'
-            write(6,*) '                 -----------------------'
-            write(6,*) 
-            write(6,912) egas
-            write(6,913) esol
-            write(6,914) (egas-esol),(egas-esol)*627.509451d+00
+            write(luout,*) '                 COSMO solvation results'
+            write(luout,*) '                 -----------------------'
+            write(luout,*) 
+            write(luout,912) egas
+            write(luout,913) esol
+            write(luout,914) (egas-esol),(egas-esol)*627.509451d+00
  912     format('                 gas phase energy = ',f20.10)
  913     format('                 sol phase energy = ',f20.10)
  914     format(' (electrostatic) solvation energy = ',f20.10,
@@ -293,7 +295,7 @@ c
       if (.not. ga_destroy(g_xorth)) call errquit
      $     ('rohf_diis_solve: ga_destroy xorth', 0, GA_ERR)
 c     
-      if (nodezero) call util_flush(6)
+      if (nodezero) call util_flush(luout)
       call ga_sync()
 c     
       rohf_diis_solve = converged
@@ -444,24 +446,32 @@ c
       integer g_xorth
 c
       integer l_e, k_e, i
+      integer l_scr,k_scr
       double precision xx
 c
       if (.not.ma_push_get(MT_DBL, nbf, 'overlap', l_e, k_e ))
      $     call errquit('rohf_ao_orthonorm: cannot allocate over',0,
-     &       MA_ERR)
+     &     MA_ERR)
+      if (.not.ma_push_get(MT_DBL, nbf, 'scratch', l_scr, k_scr ))
+     $     call errquit('rohf_ao_orthonorm: cannot allocate over',0,
+     &     MA_ERR)
       call util_diag_std( g_over, g_xorth, dbl_mb(k_e))
 c               1/2
 c   X   = U  / s
 c    ij    ij   j
 c
       do i=ga_nodeid()+1,nbf,ga_nnodes()
-        xx = 1.d0/sqrt(dbl_mb(k_e+i-1))
-        call ga_dscal_patch( g_xorth, 1, nbf, i, i, xx )
+         xx = 1.d0/sqrt(dbl_mb(k_e+i-1))
+c     call ga_dscal_patch( g_xorth, 1, nbf, i, i, xx )
+        call ga_get( g_xorth, 1, nbf, i, i, dbl_mb(k_scr), nbf )
+        call dscal(nbf,xx,dbl_mb(k_scr),1)
+        call ga_put( g_xorth, 1, nbf, i, i, dbl_mb(k_scr), nbf )
       enddo
+      call ga_sync()
 c
 c
 c
-      if (.not.ma_pop_stack(l_e))
+      if (.not.ma_chop_stack(l_e))
      $     call errquit('rohf_ao_orthonorm: cannot pop stack',0,
      &       MA_ERR)
       return

src/ddscf/rohf_fock.F

@@ -60,7 +60,7 @@ c
       integer g_cos1e
       dimension g_cos1e(3)
       double precision ecosmo
-      logical osome
+      logical do_x2c
 c
       logical odebug
       integer itype_wfn
@@ -160,13 +160,16 @@ c
          endif
         call ga_dadd(1.d0,g_hcore,1.d0,g_zora_Kinetic(1),g_hcore) ! zora kinetic
       endif
+c     == X2C options ==
+      if (rtdb_get(rtdb,'x2c:on',mt_log,1,do_x2c)) then
+            if(do_x2c) call
+     A           errquit('X2C available only for DFT ',0,CAPMIS_ERR)
+      endif
       call int_1e_ga(basis, basis, g_hcore, 'potential', oskel) ! potential
       if (util_print('rohf_fock',print_debug)) call ga_print(g_hcore)
 c
 c     cosmo charges, potential and energy contribution
 c
-      osome=.false.
-      osome=osome.and.ga_nodeid().eq.0
       if(cosmo_on.and.cosmo_phase.eq.2) then
        cosmo_file = "cosmo.xyz"
        g_dens = ga_create_atom_blocked(geom,basis,'density matrix')

src/ddscf/scf.F

@@ -14,6 +14,7 @@ C$Id$
 #include "cscf.fh"
 #include "case.fh"
 #include "frozemb.fh"
+#include "stdio.fh"
 c
 c     ROHF module.
 c
@@ -85,13 +86,14 @@ c
 c     Print info
 c 
       if (ga_nodeid().eq.0 .and. oprint) then
-         call util_print_centered(6, 'NWChem SCF Module', 40, .true.)
-         write(6,*)
-         write(6,*)
+         call util_print_centered(luout, 'NWChem SCF Module', 40,
+     $                            .true.)
+         write(luout,*)
+         write(luout,*)
          if (title .ne. ' ') then
-            call util_print_centered(6, title, 40, .false.)
-            write(6,*)
-            write(6,*)
+            call util_print_centered(luout, title, 40, .false.)
+            write(luout,*)
+            write(luout,*)
          endif
 c
          if(.not. geom_systype_get(geom,itype))
@@ -102,7 +104,7 @@ c
          if (.not. geom_ncent(geom, natoms)) call errquit
      $        ('scf: geom_ncent failed', 0, GEOM_ERR)
 c
-         write(6,1) trans(1:inp_strlen(trans)), nbf
+         write(luout,1) trans(1:inp_strlen(trans)), nbf
  1       format(/
      $        '  ao basis        = "',a,'"'/
      $        '  functions       = ', i5)
@@ -111,23 +113,23 @@ c
             if (.not. bas_name(riscf_basis, name, ri_trans))
      $           call errquit('scf: bas_name?', 0,
      &       BASIS_ERR)
-            write(6,11) ri_trans(1:inp_strlen(ri_trans)), nff
+            write(luout,11) ri_trans(1:inp_strlen(ri_trans)), nff
  11         format(
      $           '  ri basis        = "',a,'"'/
      $           '  ri functions    = ', i5)
          endif
-         write(6,12) natoms
+         write(luout,12) natoms
  12      format('  atoms           = ', i5)
          if (scftype .eq. 'UHF') then
-            write(6,121) nalpha, nbeta
+            write(luout,121) nalpha, nbeta
  121        format('  alpha electrons = ', i5/
      $           '  beta  electrons = ', i5)
          else
-            write(6,122) nclosed, nopen
+            write(luout,122) nclosed, nopen
  122        format('  closed shells   = ', i5/
      $           '  open shells     = ', i5)
          endif
-         write(6,123)
+         write(luout,123)
      $        charge, scftype,
      $        movecs_in(1:inp_strlen(movecs_in)),
      $        movecs_out(1:inp_strlen(movecs_out)), oskel, oadapt
@@ -138,14 +140,14 @@ c
      $        '  output vectors  = ', a/
      $        '  use symmetry    = ', l1/
      $        '  symmetry adapt  = ', l1)
-         if (olock) write(6,124) olock
+         if (olock) write(luout,124) olock
  124     format('  lock orbitals   = ', l1)
          if (util_print('basis summary', print_default)) then
             if (.not.bas_summary_print(basis))
      &           call errquit('scf: basis summary print failed',911,
      &       BASIS_ERR)
          else
-            write(6,*)
+            write(luout,*)
          endif
          if (oadapt .and. util_print('char table', print_high)) then
             call sym_print_char_table(geom)
@@ -153,7 +155,7 @@ c
          if (oadapt .and. util_print('basis sym', print_default)) then
             call sym_bas_irreps(basis, .true., nbf_per_ir)
          endif
-         call util_flush(6)
+         call util_flush(luout)
       endif
 c
 c     Print out other info as requested by user
@@ -162,17 +164,17 @@ c
          if (util_print('geometry',print_high)) then
             if (.not. geom_print(geom))
      $        call errquit('scf: geom_print', 0, GEOM_ERR)
-            write(6,*)
+            write(luout,*)
          endif
          if (util_print('symmetry', print_debug)) then
             call sym_print_all(geom, .true., .true., .true., 
      $           .true., .true.)
-            write(6,*)
+            write(luout,*)
          endif
          if (util_print('basis', print_high)) then
             if (.not. bas_print(basis)) 
      $           call errquit('scf: bas_print', 0, BASIS_ERR)
-            write(6,*)
+            write(luout,*)
          endif
          if (util_print('basis labels',print_high) .and.
      $        ga_nodeid().eq.0) then
@@ -181,9 +183,9 @@ c
          if (util_print('geombas', print_debug)) then
             if (.not. gbs_map_print(basis)) 
      $           call errquit('scf:gmp', 0, BASIS_ERR)
-            write(6,*)
+            write(luout,*)
          endif
-         call util_flush(6)
+         call util_flush(luout)
       endif
 c
 c     Get and check the Coulomb Attenuation Method (CAM) parameters
@@ -208,10 +210,10 @@ c
      $           call errquit('scf: failed getting converged energy',0,
      &       RTDB_ERR)
             if (ga_nodeid().eq.0 .and. oprint) then
-               write(6,1101) energy
+               write(luout,1101) energy
  1101          format(/'  The SCF is already converged '//,
      $              '         Total SCF energy =', f20.12/)
-               call util_flush(6)
+               call util_flush(luout)
             endif
             goto 3131
          endif
@@ -225,14 +227,14 @@ c
      &       UNKNOWN_ERR)
         if (ga_nodeid().eq.0 .and. oprint) then
           if (nriscf.eq.1) then
-            call util_print_centered(6, 'RI hessian', 40, .true.)
+            call util_print_centered(luout, 'RI hessian', 40, .true.)
           else if (nriscf.eq.2) then
-            call util_print_centered(6, 'RISCF', 40, .true.)
+            call util_print_centered(luout, 'RISCF', 40, .true.)
           else if (nriscf.eq.3) then
-            call util_print_centered(6, 'preconverge with RISCF', 40, 
-     $           .true.)
+            call util_print_centered(luout, 'preconverge with RISCF',
+     $           40, .true.)
           endif
-          write(6,*)
+          write(luout,*)
         endif
       endif
 c
@@ -248,7 +250,7 @@ c
 c
 c     Form intial guess vectors
 c
-      if (ga_nodeid().eq.0 .and. oprint) write(6,711) util_wallsec()
+      if (ga_nodeid().eq.0 .and. oprint) write(luout,711) util_wallsec()
  711  format(/' Forming initial guess at ',f9.1,'s'/)
 c      call scf_vectors_guess(rtdb)
 c
@@ -266,7 +268,7 @@ c     Here branch to the appropriate wavefunction type ... returns
 c     true if happily converged, false otherwise (e.g., if restart
 c     is needed or problems encountered).
 c
-      if (ga_nodeid().eq.0 .and. oprint) write(6,712) util_wallsec()
+      if (ga_nodeid().eq.0 .and. oprint) write(luout,712) util_wallsec()
  712  format(/' Starting SCF solution at ',f9.1,'s'/)
 c
 c     frozen embedding 
@@ -289,9 +291,9 @@ c     reconverge with exact energy/gradient if requested
 c
       if (nriscf.eq.3) then
         if (ga_nodeid().eq.0 .and. oprint) then
-          call util_print_centered(6,
+          call util_print_centered(luout,
      $         'switching to exact energy/gradient', 40, .true.)
-          write(6,*)
+          write(luout,*)
         endif
         nriscf = 1
         if (scftype .eq. 'UHF') then

src/ddscf/uhf.F

@@ -12,6 +12,7 @@ C$Id$
 #include "cscf.fh"
 #include "geom.fh"
 #include "bas.fh"
+#include "stdio.fh"
 c
 c     Compute UHF wavefunction .
 c
@@ -80,16 +81,16 @@ c
 c
       if (ga_nodeid().eq.0 .and. (oprint .or. .not.converged)) then
          if (.not. converged) then
-            write(6,*)
-            call util_print_centered(6,
+            write(luout,*)
+            call util_print_centered(luout,
      $           'Calculation failed to converge', 20, .true.)
-            write(6,*)
+            write(luout,*)
          end if
-         write(6,2) scftype, energy, eone, etwo, enrep  
+         write(luout,2) scftype, energy, eone, etwo, enrep
          if (abs(ecosmo).gt.0.0d0) then
-           write(6,3) ecosmo
+           write(luout,3) ecosmo
          endif
-         write(6,4) sz, sz*(sz+1), s2, uhf_time
+         write(luout,4) sz, sz*(sz+1), s2, uhf_time
  2       format(//
      $        '       Final ',a4,' results '/
      $        '       ------------------ '//
@@ -104,7 +105,7 @@ c
      $        '                 Sz(Sz+1) =', f12.4/
      $        '                      S^2 =', f12.4//
      $        '        Time for solution =', f9.1,'s'//)
-         call util_flush(6)
+         call util_flush(luout)
          call ecce_print1('total energy', mt_dbl, energy, 1)
          call ecce_print1('one-electron energy', mt_dbl, eone, 1)
          call ecce_print1('two-electron energy', mt_dbl, etwo, 1)
@@ -134,14 +135,14 @@ C
       endif
       if (ga_nodeid() .eq. 0) then
          if (util_print('final evals', print_default)) then
-            call util_print_centered(6,'Final alpha eigenvalues',
+            call util_print_centered(luout,'Final alpha eigenvalues',
      $           20,.true.)
             call output(dbl_mb(k_eval), 1, nprint, 1, 1, nmo, 1, 1)
-            write(6,*)
-            call util_print_centered(6,'Final beta eigenvalues',
+            write(luout,*)
+            call util_print_centered(luout,'Final beta eigenvalues',
      $           20,.true.)
             call output(dbl_mb(k_eval+nbf), 1, nprint, 1, 1, nmo, 1, 1)
-            call util_flush(6)
+            call util_flush(luout)
          end if
       endif
       if (util_print('final vectors analysis', print_default)) then
@@ -164,10 +165,10 @@ C
       endif
       if (ga_nodeid() .eq. 0) then
          if (util_print('final vectors', print_debug)) then
-            write(6,*)
-            call util_print_centered(6,'Final MO vectors',40,.true.)
-            write(6,*)
-            call util_flush(6)
+            write(luout,*)
+            call util_print_centered(luout,'Final MO vectors',40,.true.)
+            write(luout,*)
+            call util_flush(luout)
          end if
       end if
       if (util_print('final vectors', print_debug)) then
@@ -240,6 +241,7 @@ C     $Id$
 #include "cscfps.fh"
 #include "util.fh"
 #include "cscf.fh"
+#include "stdio.fh"
 c     
 c     Solve the UHF equations using a hybrid NR/PCG method
 c     
@@ -282,7 +284,7 @@ c
 C      external uhf_hessv, uhf_hessv_precond
 c     
       
-      if (ga_nodeid().eq.0) call util_flush(6)
+      if (ga_nodeid().eq.0) call util_flush(luout)
 c     
 c     Allocate global arrays.  Gradient, search direction
 c     and work space for the PCG routine.
@@ -309,9 +311,9 @@ c
       tlastwrite = util_wallsec()
 c     
       if (ga_nodeid().eq.0.and. oprint_parm) then
-         write(6,1) gnorm_tol, maxiter, tol2e
-         if (ouser_changed_conv) write(6,11) shifts, nr_gswitch
-         write(6,111)
+         write(luout,1) gnorm_tol, maxiter, tol2e
+         if (ouser_changed_conv) write(luout,11) shifts, nr_gswitch
+         write(luout,111)
  1       format(//,
      $        1x,'----------------------------------------------',/
      $        1x,'        Quadratically convergent UHF',//,
@@ -328,7 +330,7 @@ c
      $        1x,'NR  enabled at maxg       :',9x,f10.3)
  111     format(
      $        1x,'----------------------------------------------',/)
-         call util_flush(6)
+         call util_flush(luout)
       end if
 c
       ododiag = .true.
@@ -388,11 +390,11 @@ c     go back to top of loop
 c     
                if (oprint_vecs) then
                   if (ga_nodeid() .eq. 0) then
-                     write(6,*)
-                     call util_print_centered(6,
+                     write(luout,*)
+                     call util_print_centered(luout,
      $                    'Intermediate MO vectors',40,.true.)
-                     write(6,*)
-                     call util_flush(6)
+                     write(luout,*)
+                     call util_flush(luout)
                   end if
                   call ga_print(g_movecs)
                   call ga_print(g_movecs(2))
@@ -413,7 +415,7 @@ c
 c     
 c     End SCF minimisation
 c     
-      if (ga_nodeid().eq.0) call util_flush(6)
+      if (ga_nodeid().eq.0) call util_flush(luout)
 c
       if (.not.noscf) then  ! bypass for noscf
 c     
@@ -442,7 +444,7 @@ c
       if (.not. ga_destroy(g_work)) call errquit
      $     ('uhf_nr_solve: ga_destroy work', 0, GA_ERR)
 c     
-      if (ga_nodeid().eq.0) call util_flush(6)
+      if (ga_nodeid().eq.0) call util_flush(luout)
       call ga_sync()
 c     
       uhf_nr_solve = converged
@@ -697,6 +699,7 @@ cc AJL/Begin/SPIN ECPs
       external bas_get_ecp_handle
       logical ecp_get_high_chan
       external ecp_get_high_chan
+      logical do_x2c
 cc AJL/End
 c     
 c ... jochen 05/20: added an option to print contributions to
@@ -800,10 +803,10 @@ c
       call ga_maxelt(g_b_coul, errmaxb)
       if (max(errmaxa,errmaxb).gt.1d-4) then
          if (ga_nodeid().eq.0) then
-            write(6,77) errmaxa,errmaxb
+            write(luout,77) errmaxa,errmaxb
  77         format(' Warning: spatial symmetry breaking in UHF: ',
      $           1p,2d9.2)
-            call util_flush(6)
+            call util_flush(luout)
          endif
       endif
 c
@@ -937,9 +940,9 @@ c       Coulomb, alpha spin MOs
           dbl_mb(k_eps(1)+imo-1) = dbl_mb(k_eps(1)+imo-1)
      &      + dbl_mb(k_diag+imo-1)
         end do
-        if (master) write(6,*) 'V(C) alfa MOs'
+        if (master) write(luout,*) 'V(C) alfa MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_eps(1)+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_eps(1)+imo-1)
         end do
 c
 c       Coulomb, beta spin MOs
@@ -960,9 +963,9 @@ c       Coulomb, beta spin MOs
           dbl_mb(k_eps(2)+imo-1) = dbl_mb(k_eps(2)+imo-1)
      &      + dbl_mb(k_diag+imo-1)
         end do
-        if (master) write(6,*) 'V(C) beta MOs'
+        if (master) write(luout,*) 'V(C) beta MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_eps(2)+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_eps(2)+imo-1)
         end do
 c        
 c       XC, alfa spin MOs
@@ -971,9 +974,9 @@ c       XC, alfa spin MOs
         call two_index_transf(g_a_exch, g_vecs(1), g_vecs(1),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'V(XC) alfa MOs'
+        if (master) write(luout,*) 'V(XC) alfa MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, -dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, -dbl_mb(k_diag+imo-1)
         end do
         do imo = 1,nbf
           dbl_mb(k_eps(1)+imo-1) = dbl_mb(k_eps(1)+imo-1)
@@ -986,9 +989,9 @@ c       XC, beta spin MOs
         call two_index_transf(g_b_exch, g_vecs(2), g_vecs(2),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'V(XC) beta MOs'
+        if (master) write(luout,*) 'V(XC) beta MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, -dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, -dbl_mb(k_diag+imo-1)
         end do
         do imo = 1,nbf
           dbl_mb(k_eps(2)+imo-1) = dbl_mb(k_eps(2)+imo-1)
@@ -997,9 +1000,9 @@ c       XC, beta spin MOs
       end if ! epsana
 c
       if (odebug .and. ga_nodeid().eq.0) then
-         write(6,*) ' coulomb energies', e_a_coul, e_b_coul
-         write(6,*) ' exchang energies', e_a_exch, e_b_exch
-         call util_flush(6)
+         write(luout,*) ' coulomb energies', e_a_coul, e_b_coul
+         write(luout,*) ' exchang energies', e_a_exch, e_b_exch
+         call util_flush(luout)
       endif
       if (odebug) then
          call ga_print(g_a_coul)
@@ -1054,6 +1057,11 @@ c
          endif
         call ga_dadd(1.d0,g_a_hcore,1.d0,g_zora_Kinetic(1),g_a_hcore) ! zora kinetic
       endif
+c     == X2C options ==
+      if (rtdb_get(rtdb,'x2c:on',mt_log,1,do_x2c)) then
+            if(do_x2c) call
+     A           errquit('X2C available only for DFT ',0,CAPMIS_ERR)
+      endif
 
       if (epsana) then
 c       MO energy analysis: print T contribs, but do not yet
@@ -1065,9 +1073,9 @@ c       T, alfa spin MOs
         call two_index_transf(g_a_hcore, g_vecs(1), g_vecs(1),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'T alfa MOs'
+        if (master) write(luout,*) 'T alfa MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_diag+imo-1)
         end do
 c       T, beta spin MOs
         call ga_zero(g_tmp(1))
@@ -1075,9 +1083,9 @@ c       T, beta spin MOs
         call two_index_transf(g_a_hcore, g_vecs(2), g_vecs(2),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'T beta MOs'
+        if (master) write(luout,*) 'T beta MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_diag+imo-1)
         end do
       end if ! epsana
 cc If spin polarised ECP, split g_hcore
@@ -1103,9 +1111,9 @@ c       alfa spin MOs
         call two_index_transf(g_a_hcore, g_vecs(1), g_vecs(1),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'T+V(nuc) alfa MOs'
+        if (master) write(luout,*) 'T+V(nuc) alfa MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_diag+imo-1)
           dbl_mb(k_eps(1)+imo-1) = dbl_mb(k_eps(1)+imo-1)
      &      + dbl_mb(k_diag+imo-1)
         end do
@@ -1116,16 +1124,16 @@ c       beta spin MOs
         call two_index_transf(g_b_hcore, g_vecs(2), g_vecs(2),
      &    g_tmp(1), g_tmp(2))
         call ga_get_diagonal(g_tmp(2),dbl_mb(k_diag))
-        if (master) write(6,*) 'T+V(nuc) beta MOs'
+        if (master) write(luout,*) 'T+V(nuc) beta MOs'
         do imo=1,nbf
-          if (master) write (6,*) imo, dbl_mb(k_diag+imo-1)
+          if (master) write (luout,*) imo, dbl_mb(k_diag+imo-1)
           dbl_mb(k_eps(2)+imo-1) = dbl_mb(k_eps(2)+imo-1)
      &      + dbl_mb(k_diag+imo-1)
         end do
-        if (master) write (6,'(//1x,a)')
+        if (master) write (luout,'(//1x,a)')
      &    'spin MO energies with T + V(nuc) + V(C) + V(XC)'
         do imo = 1,nbf
-          if (master) write (6,*) imo,
+          if (master) write (luout,*) imo,
      &       dbl_mb(k_eps(1)+imo-1) ,  dbl_mb(k_eps(2)+imo-1)
         end do
       end if ! epsana      
@@ -1139,7 +1147,7 @@ c     check if we are doing dft (e.g. cphf case)
      &        INPUT_ERR)
          if (theory.eq.'hyb') theory = 'dft'
          if(theory.ne.'dft') then
-            if(master) write(6,*) ' UHF theory ',theory
+            if(master) write(luout,*) ' UHF theory ',theory
             call errquit('uhf_energy not compatible with COSMO',0,
      C           CAPMIS_ERR)
          endif
@@ -1237,10 +1245,10 @@ c ... jochen: full MO energies, with all additional contributions
          call ga_get_diagonal(cuhf_g_falpha,dbl_mb(k_eps(1)))
          call ga_get_diagonal(cuhf_g_fbeta,dbl_mb(k_eps(2)))
          
-         if (master) write (6,'(//1x,a)')
+         if (master) write (luout,'(//1x,a)')
      &     'spin MO energies with all contributions (COSMO, etc.)'
          do imo = 1,nbf
-           if (master) write (6,*) imo,
+           if (master) write (luout,*) imo,
      &       dbl_mb(k_eps(1)+imo-1) ,  dbl_mb(k_eps(2)+imo-1)
          end do
 
@@ -1257,7 +1265,7 @@ c        shorter variable names
          na = nalpha
          nb = nbeta
 
-         if (master) write(6,'(/1x,a/)')
+         if (master) write(luout,'(/1x,a/)')
      &     'Some ERIs, Mulliken notation [1 1* | 2 2*]'
 
 c        ------------------------------------------------------         
@@ -1282,13 +1290,13 @@ c        Fock 2e matrix from alfa HOMO-1 density
          
          call ga_get(f(2), na-1, na-1, na-1, na-1, rtemp,1)
  1001    format(1x,"[",i3,"a,",i3,"a |",i3,"a,",i3,"a]=",f15.8)
-         if (master) write(6,1001) na-1,na-1,na-1,na-1,rtemp
+         if (master) write(luout,1001) na-1,na-1,na-1,na-1,rtemp
 
          call ga_get(f(2), na, na, na, na, rtemp,1)
-         if (master) write(6,1001) na-1,na-1,na,na,rtemp
+         if (master) write(luout,1001) na-1,na-1,na,na,rtemp
          
          call ga_get(f(2), na+1, na+1, na+1, na+1, rtemp,1)
-         if (master) write(6,1001) na-1,na-1,na+1,na+1,rtemp
+         if (master) write(luout,1001) na-1,na-1,na+1,na+1,rtemp
 
          call two_index_transf(f(1), g_vecs(2), g_vecs(2),
      $     d(2), f(2))
@@ -1296,10 +1304,10 @@ c        Fock 2e matrix from alfa HOMO-1 density
          call ga_sync()
          
          call ga_get(f(2), nb, nb, nb, nb, rtemp,1)
-         if (master) write(6,1001) na-1,na-1,nb,nb,rtemp
+         if (master) write(luout,1001) na-1,na-1,nb,nb,rtemp
          
          call ga_get(f(2), nb+1, nb+1, nb+1, nb+1, rtemp,1)
-         if (master) write(6,1001) na-1,na-1,nb+1,nb+1,rtemp            
+         if (master) write(luout,1001) na-1,na-1,nb+1,nb+1,rtemp
 
 c        ------------------------------------------------------         
 c        Fock 2e matrix from alfa HOMO density
@@ -1322,13 +1330,13 @@ c        Fock 2e matrix from alfa HOMO density
          call ga_sync()
          
          call ga_get(f(2), na-1, na-1, na-1, na-1, rtemp,1)
-         if (master) write(6,1001) na,na,na-1,na-1,rtemp
+         if (master) write(luout,1001) na,na,na-1,na-1,rtemp
          
          call ga_get(f(2), na, na, na, na, rtemp,1)
-         if (master) write(6,1001) na,na,na,na,rtemp
+         if (master) write(luout,1001) na,na,na,na,rtemp
          
          call ga_get(f(2), na+1, na+1, na+1, na+1, rtemp,1)
-         if (master) write(6,1001) na,na,na+1,na+1,rtemp
+         if (master) write(luout,1001) na,na,na+1,na+1,rtemp
 
          call two_index_transf(f(1), g_vecs(2), g_vecs(2),
      $     d(2), f(2))
@@ -1336,10 +1344,10 @@ c        Fock 2e matrix from alfa HOMO density
          call ga_sync()
          
          call ga_get(f(2), nb, nb, nb, nb, rtemp,1)
-         if (master) write(6,1001) na,na,nb,nb,rtemp
+         if (master) write(luout,1001) na,na,nb,nb,rtemp
          
          call ga_get(f(2), nb+1, nb+1, nb+1, nb+1, rtemp,1)
-         if (master) write(6,1001) na,na,nb+1,nb+1,rtemp         
+         if (master) write(luout,1001) na,na,nb+1,nb+1,rtemp         
 
          
 c        ------------------------------------------------------
@@ -1363,13 +1371,13 @@ c        Fock 2e matrix from alfa LUMO density
          call ga_sync()
          
          call ga_get(f(2), na-1, na-1, na-1, na-1, rtemp,1)
-         if (master) write(6,1001) na+1,na+1,na-1,na-1,rtemp         
+         if (master) write(luout,1001) na+1,na+1,na-1,na-1,rtemp         
 
          call ga_get(f(2), na, na, na, na, rtemp,1)
-         if (master) write(6,1001) na+1,na+1,na,na,rtemp         
+         if (master) write(luout,1001) na+1,na+1,na,na,rtemp         
 
          call ga_get(f(2), na+1, na+1, na+1, na+1, rtemp,1)
-         if (master) write(6,1001) na+1,na+1,na+1,na+1,rtemp         
+         if (master) write(luout,1001) na+1,na+1,na+1,na+1,rtemp         
 
          call two_index_transf(f(1), g_vecs(2), g_vecs(2),
      $     d(2), f(2))
@@ -1377,10 +1385,10 @@ c        Fock 2e matrix from alfa LUMO density
          call ga_sync()
          
          call ga_get(f(2), nb, nb, nb, nb, rtemp,1)
-         if (master) write(6,1001) na+1,na+1,nb,nb,rtemp
+         if (master) write(luout,1001) na+1,na+1,nb,nb,rtemp
          
          call ga_get(f(2), nb+1, nb+1, nb+1, nb+1, rtemp,1)
-         if (master) write(6,1001) na+1,na+1,nb+1,nb+1,rtemp             
+         if (master) write(luout,1001) na+1,na+1,nb+1,nb+1,rtemp             
 
 
 c        ------------------------------------------------------
@@ -1410,11 +1418,11 @@ c        Fock 2e matrix from SYMMETRIZED alfa HOMO-HOMO-1 product
 
          rtemp = zero
          call ga_get(f(2), na-1, na-1, na, na, rtemp,1)
-         if (master) write(6,1001) na,na-1,na,na-1,rtemp
+         if (master) write(luout,1001) na,na-1,na,na-1,rtemp
          
          rtemp = zero
          call ga_get(f(2), na, na, na-1, na-1, rtemp,1)
-         if (master) write(6,1001) na-1,na,na-1,na,rtemp
+         if (master) write(luout,1001) na-1,na,na-1,na,rtemp
 
 
 c        ------------------------------------------------------
@@ -1442,11 +1450,11 @@ c        Fock 2e matrix from SYMMETRIZED alfa LUMO-HOMO-1 product
          
          rtemp = zero
          call ga_get(f(2), na-1, na-1, na+1, na+1, rtemp,1)
-         if (master) write(6,1001) na+1,na-1,na+1,na-1,rtemp         
+         if (master) write(luout,1001) na+1,na-1,na+1,na-1,rtemp         
 
          rtemp = zero
          call ga_get(f(2), na+1, na+1, na-1, na-1, rtemp,1)
-         if (master) write(6,1001) na-1,na+1,na-1,na+1,rtemp  
+         if (master) write(luout,1001) na-1,na+1,na-1,na+1,rtemp  
          
 c        ------------------------------------------------------
 c        Fock 2e matrix from SYMMETRIZED alfa LUMO-HOMO product
@@ -1475,11 +1483,11 @@ c        Fock 2e matrix from SYMMETRIZED alfa LUMO-HOMO product
 
          rtemp = zero
          call ga_get(f(2), na, na, na+1, na+1, rtemp,1)
-         if (master) write(6,1001) na+1,na,na+1,na,rtemp
+         if (master) write(luout,1001) na+1,na,na+1,na,rtemp
 
          rtemp = zero
          call ga_get(f(2), na+1, na+1, na, na, rtemp,1)
-         if (master) write(6,1001) na,na+1,na,na+1,rtemp
+         if (master) write(luout,1001) na,na+1,na,na+1,rtemp
 
 
 c        ------------------------------------------------------
@@ -1509,11 +1517,11 @@ c        Fock 2e matrix from SYMMETRIZED beta LUMO-HOMO product
 
          rtemp = zero
          call ga_get(f(2), nb, nb, nb+1, nb+1, rtemp,1)
-         if (master) write(6,1001) nb+1,nb,nb+1,nb,rtemp
+         if (master) write(luout,1001) nb+1,nb,nb+1,nb,rtemp
 
          rtemp = zero
          call ga_get(f(2), nb+1, nb+1, nb, nb, rtemp,1)
-         if (master) write(6,1001) nb,nb+1,nb,nb+1,rtemp
+         if (master) write(luout,1001) nb,nb+1,nb,nb+1,rtemp
          
          
 c        ------------------------------------------------------
@@ -1613,7 +1621,7 @@ c     DIM energy term
       end if
 c
       if (odebug .and. ga_nodeid().eq.0) then
-         write(6,*) ' eone, etwo, enrep, energy ',
+         write(luout,*) ' eone, etwo, enrep, energy ',
      $        eone, etwo, enrep, energy
       endif
 c
@@ -1660,6 +1668,7 @@ c
 #include "util.fh"
 #include "mafdecls.fh"
 #include "rtdb.fh"
+#include "stdio.fh"
       integer rtdb
       integer g_grad
       integer g_work
@@ -1732,10 +1741,10 @@ c
       if (lshift .lt. min_shift) then
          lshift = min_shift + 2.0d0
          if (ga_nodeid().eq.0 .and. oprint_conv) then
-            write(6,3131) lshift
+            write(luout,3131) lshift
  3131       format('  Setting level-shift to ', f6.2,
      $           ' to force positive preconditioner')
-            call util_flush(6)
+            call util_flush(luout)
          end if
       end if
 c
@@ -1766,9 +1775,9 @@ c
                odisable_nr = .true.
                maxiter = maxiter + 10
                if (ga_nodeid() .eq. 0 .and. oprint_conv) then
-                  write(6,22) maxiter
+                  write(luout,22) maxiter
  22               format(/' Disabled NR: increased maxiter to ',i3/)
-                  call util_flush(6)
+                  call util_flush(luout)
                endif
                goto 30
             else
@@ -1780,9 +1789,9 @@ c
                endif
             endif
             if (ga_nodeid() .eq. 0 .and. oprint_conv) then
-               write(6,2) lshift
+               write(luout,2) lshift
  2             format(' Increased level shift to ', f8.2)
-               call util_flush(6)
+               call util_flush(luout)
             endif
             goto 20
          endif
@@ -1945,6 +1954,7 @@ C$Id$
 #include "cscfps.fh"
 #include "cscf.fh"
 #include "cuhf.fh"
+#include "stdio.fh"
       logical oaufbau
       logical oprint
 c
@@ -2038,18 +2048,18 @@ c
       end do
 c
       if (oprint .and. ga_nodeid().eq.0) then
-         write(6,*)
-         write(6,*)
-         call util_print_centered(6, 'Alpha-spin eigenvalues',
+         write(luout,*)
+         write(luout,*)
+         call util_print_centered(luout, 'Alpha-spin eigenvalues',
      $        20, .true.)
          call output(dbl_mb(k_eval), 1, min(nalpha+5,nmo), 
      $        1, 1, nmo, 1, 1)
-         write(6,*)
-         call util_print_centered(6, 'Beta-spin eigenvalues',
+         write(luout,*)
+         call util_print_centered(luout, 'Beta-spin eigenvalues',
      $        20, .true.)
          call output(dbl_mb(k_eval+nbf), 1, min(nalpha+5,nmo), 
      $        1, 1, nmo, 1, 1)
-         call util_flush(6)
+         call util_flush(luout)
       end if
 c
       if (.not. ga_destroy(g_u))
@@ -2124,6 +2134,7 @@ c
 #include "rtdb.fh"
 #include "mafdecls.fh"
 #include "bas.fh"
+#include "stdio.fh"
       integer rtdb
 c
 c     Analyze the UHF wavefunction.  Optionally print the mulliken
@@ -2178,31 +2189,31 @@ c
 c
       if (oprintmulliken) then
          if (ga_nodeid() .eq. 0) then
-            write(6,*)
-            call util_print_centered(6,
+            write(luout,*)
+            call util_print_centered(luout,
      $           'Mulliken analysis of the total density', 20,.true.)
          endif
          call mull_pop(geom, basis, g_dens, g_over, 'total')
 c
          if (ga_nodeid() .eq. 0) then
-            write(6,*)
-            call util_print_centered(6,
+            write(luout,*)
+            call util_print_centered(luout,
      $           'Mulliken analysis of the alpha density', 20,.true.)
          endif
          call mull_pop(geom, basis, g_adens, g_over, 'alpha')
 c     
          if (nbeta .gt. 0) then
             if (ga_nodeid() .eq. 0) then
-               write(6,*)
-               call util_print_centered(6,
+               write(luout,*)
+               call util_print_centered(luout,
      $              'Mulliken analysis of the beta density', 
      $              20,.true.)
             endif
             call mull_pop(geom, basis, g_bdens, g_over, 'beta')
 c
             if (ga_nodeid() .eq. 0) then
-               write(6,*)
-               call util_print_centered(6,
+               write(luout,*)
+               call util_print_centered(luout,
      $              'Mulliken analysis of the spin density', 
      $              20,.true.)
             endif
@@ -2235,14 +2246,14 @@ c
       enddo
 c
       if (oprintmultipole .and. ga_nodeid().eq.0) then
-         write(6,*)
-         call util_print_centered(6,
+         write(luout,*)
+         call util_print_centered(luout,
      $     'Multipole analysis of the density wrt the origin',
      $      30, .true.)
-         write(6,*)
-         write(6,*) '    L   x y z        total         alpha',
+         write(luout,*)
+         write(luout,*) '    L   x y z        total         alpha',
      $                 '         beta         nuclear'
-         write(6,*) '    -   - - -        -----         -----',
+         write(luout,*) '    -   - - -        -----         -----',
      $                 '         ----         -------'
          ind = 0
          do ltotal = 0, lmax
@@ -2250,13 +2261,13 @@ c
                do l = ltotal-k, 0, -1
                   m = ltotal - k - l
                   ind = ind + 1
-                  write(6,12) ltotal, k, l, m, totalmoments(ind),
+                  write(luout,12) ltotal, k, l, m, totalmoments(ind),
      $                 alphamoments(ind), betamoments(ind),
      $                 nuclmoments(ind)
  12               format(4x,i2,2x,3i2,4f14.6)
                enddo
             enddo
-            write(6,*)
+            write(luout,*)
          enddo
       endif
 c

src/ddscf/vectors.F

@@ -55,7 +55,7 @@ c
          read(unitno, err=1001, end=2001) nsets
          read(unitno, err=1001, end=2001) nbf
          if (nsets .gt. ldnmo) then
-            write(6,*) ' movecs_read_header: ldnmo too small ',
+            write(luout,*) ' movecs_read_header: ldnmo too small ',
      $           nsets, ldnmo
             close(unitno, err=1002)
             goto 10
@@ -80,28 +80,28 @@ c
 c
       return
 c
- 1000 write(6,*) ' movecs_read_header: failed to open ',
+ 1000 write(luout,*) ' movecs_read_header: failed to open ',
      $     filename(1:inp_strlen(filename)),
      A ' IERR = ', ioserr
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 1001 write(6,*) ' movecs_read_header: failing reading from ',
+ 1001 write(luout,*) ' movecs_read_header: failing reading from ',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       close(unitno, err=1002)
       goto 10
 c
- 2001 write(6,*) ' movecs_read_header:eof: failing reading from '
+ 2001 write(luout,*) ' movecs_read_header:eof: failing reading from '
      &    , filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       close(unitno, err=1002)
       goto 10
 c
- 1002 write(6,*) ' movecs_read_header: failed to close',
+ 1002 write(luout,*) ' movecs_read_header: failed to close',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
@@ -199,27 +199,27 @@ c
       movecs_read = ok .eq. 1
       if (ga_nodeid() .eq. 0 .and. movecs_read .and.
      $     util_print('vectors i/o', print_high)) then
-         write(6,22) filename(1:inp_strlen(filename))
+         write(luout,22) filename(1:inp_strlen(filename))
  22      format(/' Read molecular orbitals from ',a/)
          call util_flush(luout)
       endif
       if (oscfps) call pstat_off(ps_vecio)
       return
 c
- 1000 write(6,*) ' movecs_read: failed to open ',
+ 1000 write(luout,*) ' movecs_read: failed to open ',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 1001 write(6,*) ' movecs_read: failing reading from ',
+ 1001 write(luout,*) ' movecs_read: failing reading from ',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       close(unitno,err=1002)
       goto 10
 c
- 1002 write(6,*) ' movecs_read: failed to close',
+ 1002 write(luout,*) ' movecs_read: failed to close',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
@@ -409,7 +409,7 @@ c get nuclear repulsion energy
 c
       if (ga_nodeid() .eq. 0) then
          if (nsets .gt. 2) then
-            write(6,*) ' movecs_write: nsets > 2 ', nsets
+            write(luout,*) ' movecs_write: nsets > 2 ', nsets
             goto 10
          endif
          open(unitno, status='unknown', form='unformatted',
@@ -470,7 +470,7 @@ c
       movecs_write = ok .eq. 1
       if (ga_nodeid() .eq. 0 .and.
      $     util_print('vectors i/o', print_high)) then
-         write(6,22) filename(1:inp_strlen(filename))
+         write(luout,22) filename(1:inp_strlen(filename))
  22      format(/' Wrote molecular orbitals to ',a/)
          call util_flush(luout)
       endif
@@ -478,20 +478,20 @@ c
       if (oscfps) call pstat_off(ps_vecio)
       return
 c
- 1000 write(6,*) ' movecs_write: failed to open ',
+ 1000 write(luout,*) ' movecs_write: failed to open ',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 1001 write(6,*) ' movecs_write: failing writing to ',
+ 1001 write(luout,*) ' movecs_write: failing writing to ',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
       close(unitno,err=1002)
       goto 10
 c
- 1002 write(6,*) ' movecs_write: failed to close',
+ 1002 write(luout,*) ' movecs_write: failed to close',
      $     filename(1:inp_strlen(filename))
       call util_flush(luout)
       ok = 0
@@ -574,7 +574,7 @@ c
          call ga_sync()
          call ga_inquire(g_vecs, type, dim1, dim2)
          if (ga_nodeid() .eq. 0) then
-            if (oprint) write(6,*)
+            if (oprint) write(luout,*)
             if (.not. ma_push_get(mt_dbl, dim1, 'swap', l_veci, k_veci))
      $           call errquit('movecs_swap: no scratch space', dim1,
      &       MA_ERR)
@@ -584,7 +584,7 @@ c
             do pair = 1, nelem, 2
                i = int_mb(index+pair-1)
                j = int_mb(index+pair  )
-               if (oprint) write(6,1) text(1:inp_strlen(text)), i, j
+               if (oprint) write(luout,1) text(1:inp_strlen(text)), i, j
  1             format(' Swapping ',a,' orbitals ', 2i5)
                if (i.lt.0 .or. i.gt.dim2) call errquit
      $              ('movecs_swap: invalid vector ', i, INPUT_ERR)
@@ -1064,13 +1064,13 @@ c
          endif
          goto 10
 c
- 1000    write(6,*) ' file_read_ga_info: failed to open ',
+ 1000    write(luout,*) ' file_read_ga_info: failed to open ',
      $        fname(1:inp_strlen(fname))
          call util_flush(luout)
          ok = 0
          goto 10
 c     
- 1001    write(6,*) ' file_read_ga_info: failing reading from ',
+ 1001    write(luout,*) ' file_read_ga_info: failing reading from ',
      $        fname(1:inp_strlen(fname))
          call util_flush(luout)
          ok = 0
@@ -1086,7 +1086,7 @@ c
          endif
          goto 10
 c     
- 1002    write(6,*) ' file_read_ga_info: failed to close',
+ 1002    write(luout,*) ' file_read_ga_info: failed to close',
      $        fname(1:inp_strlen(fname))
          call util_flush(luout)
          ok = 0
@@ -1213,19 +1213,19 @@ c
       file_write_ga_patch = ok .eq. 1
       if (ga_nodeid() .eq. 0 .and.
      $     util_print('vectors i/o', print_high)) then
-         write(6,22) gtitle(1:len1), fname(1:inp_strlen(fname))
+         write(luout,22) gtitle(1:len1), fname(1:inp_strlen(fname))
  22      format(/' Wrote ',a,' to ',a/)
          call util_flush(luout)
       endif
       return
 c
- 1000 write(6,*) ' file_write_ga: failed to open ',
+ 1000 write(luout,*) ' file_write_ga: failed to open ',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 1001 write(6,*) ' file_write_ga: failing writing to ',
+ 1001 write(luout,*) ' file_write_ga: failing writing to ',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
@@ -1241,7 +1241,7 @@ c
       endif
       goto 10
 c
- 1002 write(6,*) ' file_write_ga: failed to close',
+ 1002 write(luout,*) ' file_write_ga: failed to close',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
@@ -1360,19 +1360,19 @@ c
       file_read_ga_patch = ok .eq. 1
       if (ga_nodeid() .eq. 0 .and.
      $     util_print('vectors i/o', print_high)) then
-         write(6,22) gtitle(1:len1), fname(1:inp_strlen(fname))
+         write(luout,22) gtitle(1:len1), fname(1:inp_strlen(fname))
  22      format(/' Read ',a,' from ',a/)
          call util_flush(luout)
       endif
       return
 c
- 1000 write(6,*) ' file_read_ga: failed to open ',
+ 1000 write(luout,*) ' file_read_ga: failed to open ',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 1001 write(6,*) ' file_read_ga: failing reading from ',
+ 1001 write(luout,*) ' file_read_ga: failing reading from ',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
@@ -1388,13 +1388,13 @@ c
       endif
       goto 10
 c
- 1002 write(6,*) ' file_read_ga: failed to close',
+ 1002 write(luout,*) ' file_read_ga: failed to close',
      $     fname(1:inp_strlen(fname))
       call util_flush(luout)
       ok = 0
       goto 10
 c
- 2000 write(6,*) ' file_read_ga: GA and file ',
+ 2000 write(luout,*) ' file_read_ga: GA and file ',
      $           'contents mismatch dimensions',
      $           nrow, ncol, grow, gcol
       call util_flush(luout)

src/ddscf/vectors_inp.F

@@ -5,6 +5,7 @@ C     $Id$
 #include "inp.fh"
 #include "rtdb.fh"
 #include "mafdecls.fh"
+#include "stdio.fh"
       integer rtdb
       character*(*) module
 c     
@@ -252,13 +253,13 @@ c
 c     
       return
 c     
-10000 write(6,10001)
+10000 write(luout,10001)
 10001 format(/' vectors [[input] filename|atomic|hcore] '
      $     /'         [output  filename] [lock]'
      $     /'         [swap [alpha|beta] pair_list]'
      $     /'         [reorder pair_list]'
      $     /'         [rotate geometry filename]'/)
-      call util_flush(6)
+      call util_flush(luout)
       call errquit('vectors_input: invalid format', 0, INPUT_ERR)
 c     
       end

src/driver/coptopt.fh

@@ -73,7 +73,7 @@ c
      $     ostress2,            ! include lattice gradients
      $     ocheckgrad,          ! verify the gradient
      $     diagonal_hessian,    ! tells the code to use a diagonal hessian
-     $     ofast_diag      ,    ! use dsyev instead of jacobi eigensolver
+     $     ofast_diag      ,    ! use _dsyev_ instead of jacobi eigensolver
      $     oactive(max_cent)    ! active atom flags
 c
       character*8

src/fft/pfft1.0/pfft/wrkspc2.F

@@ -56,17 +56,14 @@ c
       integer nelebk
       INTEGER           NAPROC, ME, N1MAX, N2MAX, I
       INTEGER           LSTCNT
-      integer           max
-      INTEGER           MXNPRC, MXMYND
       integer m2, mmax,nmax, iii, nele
-      intrinsic         max
 c     
       EXTERNAL          LSTCNT
       EXTERNAL          MXNPRC, MXMYND
-*-----------------------------------------------------------------------
-c     
+      INTEGER           MXNPRC, MXMYND
       NAPROC = MXNPRC()
       ME     = MXMYND()
+ 
       INFO = 0
 c     
 c     computes the large number of panels that is stored on map1(1:n3) and

src/fft/pfftwrap/iessl.F

@@ -7,12 +7,13 @@ c************************************************************************
 #include "errquit.fh"
       INTEGER icontext
       integer nprow, npcol, myrow, mycol
-      integer nnodes, nodeid, nproc, me
+#include "nwfft3d.fh"
+      integer nproc, me
       integer icontext1, icontext2
       common /xparms/nprow,npcol,myrow,mycol
 
-      nproc = nnodes()
-      me = nodeid()
+      me = util_mpicommrank()
+      nproc = util_mpicommsize()
       nprow = 1
       npcol = nproc
 c obtain system default and map to a row

src/fft/pfftwrap/nwfft3d.F

@@ -68,8 +68,9 @@ c *** local variables
 c       integer isize
 c       parameter (isize=128*128*128)
 c       integer iwork(isize)
-c	double precision work(isize)
-       integer nodeid,nnodes,me,nproc
+c     double precision work(isize)
+#include "nwfft3d.fh"
+       integer me,nproc
 #ifdef ESSLFFT
       data iESSL /0/
       save iESSL
@@ -82,8 +83,12 @@ c	double precision work(isize)
       data computedMap /0/
       save computedMap
 *
-      me = nodeid()
-      nproc = nnodes()
+      me = util_mpicommrank()
+      nproc = util_mpicommsize()
+      if(nproc.eq.0) then
+         write(0,*) ' error me nproc ',me,nproc
+         call errquit(" error ",0,0)
+      endif
       isign=idir
       if(isign.lt.-1) isign=-1
 

src/fft/pfftwrap/nwfft3d.fh

@@ -0,0 +1,3 @@
+        integer util_mpicommsize,util_mpicommrank
+        external util_mpicommsize,util_mpicommrank
+	

src/gwmol/gw_cdgw_graph.F

@@ -127,7 +127,7 @@
           rms3(:) = 1d100
 
           !fixed contribution
-          const = Eref(llmo+iqp-1,isp) + Sigma_X -  SigmaOld
+          const = Eold(llmo+iqp-1,isp) + Sigma_X -  SigmaOld
 
           ! we only need a guess
           do iter=1,8
@@ -275,7 +275,7 @@
 
               ! Finish up building the residual
               resmat(ipoints,iqp) = resmat(ipoints,iqp) +
-     $          Eref(llmo+iqp-1,isp) + sigmax(llmo+iqp-1,isp) -
+     $          Eold(llmo+iqp-1,isp) + sigmax(llmo+iqp-1,isp) -
      $          sigma(llmo+iqp-1,isp)
 
               ! If residual is small, solution has been found

src/gwmol/gw_cdgw_newton.F

@@ -121,7 +121,8 @@
           call ga_zero(pars%g_sols)
           bracket = .false.
           rms6 = 1d100
-          const = Eref(llmo+iqp-1,isp) + Sigma_X - sigmaold
+c          const = Eref(llmo+iqp-1,isp) + Sigma_X - sigmaold
+          const = Eold(llmo+iqp-1,isp) + Sigma_X - sigmaold
           lower = -1d100
           upper = 1d100
 c        

src/libext/elpa/GNUmakefile

@@ -23,7 +23,7 @@ LIB_TARGETS += elpa
 include ../../config/makelib.h
 
 clean:
-	@rm -rf  elpa* ../lib/lib*elpa* ../include/elpa*
+	@rm -rf  elpa* ../lib/lib*elpa* ../include/elpa* ../lib/pkgconfig/elpa.pc
 
 
 

src/libext/elpa/build_elpa.sh

@@ -1,31 +1,28 @@
 #!/usr/bin/env bash
+dump_build_env(){
+echo %%%%%%% debug make failures
+env|grep -i pkg
+echo %%%%%%% debug make failures
+grep install- Makefile |grep -v uni ||true
+echo %%%%%%% debug make failures
+grep -dskip PKG *
+echo %%%%%%% debug make failures
+grep pkg Makefile ||true
+echo %%%%%%% debug make failures
+cat Makefile || true
+echo %%%%%%% debug make failures
+cat  config.log || true
+echo %%%%%%% debug make failures
+}
 #set -v
 arch=`uname -m`
-#SHORTVERSION=2021.11.001
-SHORTVERSION=2025.01.002
-#VERSION=new_release_${SHORTVERSION}
-VERSION=${SHORTVERSION}
+source ../libext_utils/getfiles_utils.sh
+get_elpa
+if [ $? -ne 0 ]; then echo "elpa download failed" ; exit 1 ;  fi
+# version set in ../libext_utils/getfiles_utils.sh
 echo mpif90 is `which mpif90`
-if [ -f  elpa-${VERSION}.tar.gz ]; then
-    echo "using existing"  elpa-${VERSION}.tar.gz
-else
-    rm -rf elpa*
-#    echo    curl -L https://github.com/marekandreas/elpa/archive/refs/tags/${VERSION}.tar.gz -o elpa-${VERSION}.tar.gz
-#    curl -L https://github.com/marekandreas/elpa/archive/refs/tags/${VERSION}.tar.gz -o elpa-${VERSION}.tar.gz
-#    echo curl https://gitlab.mpcdf.mpg.de/elpa/elpa/-/archive/${VERSION}/elpa-${VERSION}.tar.gz -o elpa-${VERSION}.tar.gz
-#    curl -L https://gitlab.mpcdf.mpg.de/elpa/elpa/-/archive/${VERSION}/elpa-${VERSION}.tar.gz -o elpa-${VERSION}.tar.gz
-echo    curl -LJO https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/${VERSION}/elpa-${VERSION}.tar.gz
-    tries=1 ; until [ "$tries" -ge 6 ] ; do
-		  if [ "$tries" -gt 1 ]; then sleep 9; echo attempt no.  $tries ; fi
-    curl -LJO https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/${VERSION}/elpa-${VERSION}.tar.gz
-    # check tar.gz integrity
-    echo check tar.gz integrity
-		  gzip -t elpa-${VERSION}.tar.gz >&  /dev/null
-		  if [ $? -eq 0 ]; then echo "download successful"; break ;  fi
-		  tries=$((tries+1)) ;  done
-fi
-tar xzf elpa-${VERSION}.tar.gz
-ln -sf elpa-${VERSION} elpa
+tar xzf elpa-*.tar.gz && rm elpa-*.tar.gz
+ln -sf elpa-* elpa
 cd elpa
 rm -f check_thread_affinity.patch
 wget https://raw.githubusercontent.com/conda-forge/elpa-feedstock/main/recipe/check_thread_affinity.patch
@@ -40,8 +37,6 @@ UNAME_S=$(uname -s)
 if [[ ${UNAME_S} == Linux ]]; then
     export ARFLAGS=rU
 fi
-MYCFLAGS+=" -Wno-error=implicit-function-declaration "
-MYCFLAGS+=" -Wno-error=format "
 if [[ ${UNAME_S} == Darwin ]]; then
     MYLINK+=" -Wl,-no_compact_unwind"
 if [[  -z "$HOMEBREW_PREFIX" ]]; then
@@ -78,6 +73,12 @@ else
 #	CFLAGS+="-I`${NWCHEM_TOP}/src/tools/guess-mpidefs --mpi_include`"
     fi
 fi
+GOTNVC=$( "$MPICC" -dM -E - </dev/null 2> /dev/null |grep __NVCOMPILER\  |cut -d " " -f 3)
+echo GOTNVC $GOTNVC
+if [[ ${GOTNVC} != 1 ]]; then
+    MYCFLAGS+=" -Wno-error=implicit-function-declaration "
+    MYCFLAGS+=" -Wno-error=format "
+fi
 if [[  -z "${FC}" ]]; then
     if [[ ! -z ${PE_ENV} ]]; then
 	FC=ftn
@@ -188,7 +189,7 @@ fi
 echo I_MPI_F90 is "$I_MPI_F90"
 
 if [[  -z "${FORCETARGET}" ]]; then
-FORCETARGET="-disable-sse -disable-sse-assembly --disable-avx --disable-avx2  --disable-avx512  "
+FORCETARGET="--enable-generic-kernels -disable-sse-kernels -disable-sse-assembly-kernels --disable-avx-kernels --disable-avx2-kernels  --disable-avx512-kernels  "
 fi #FORCETARGET
 if [[ ${CC} == icx ]] ; then
     MYCFLAGS+=" -xhost "
@@ -217,7 +218,7 @@ fi
 if [[ ${CC} == ifort ]] ; then
     MYFCFLAGS+=" -O3 -xhost "
 elif [[ ${FC} == nvfortran ]]  || [[ ${PE_ENV} == NVIDIA ]] ; then
-    MYCFLAGS+=" -tp native"
+    MYFCFLAGS+=" -tp native"
 elif [[ ${FC_EXTRA} == gfortran ]] ; then
     MYFCFLAGS+=" -O3 -g -mtune=native -march=native "
 #    MYFCFLAGS+=" -Wno-lto-type-mismatch "
@@ -225,16 +226,16 @@ elif [[ ${FC_EXTRA} == gfortran ]] ; then
 fi    
     if [[ "${GOTAVX}" == "Y" ]]; then
 	echo "using AVX instructions"
-	FORCETARGET=" --disable-sse-assembly --enable-avx --disable-avx2  --disable-avx512  "
+	FORCETARGET=" --disable-sse-assembly-kernels --enable-avx-kernels --disable-avx2-kernels  --disable-avx512-kernels  "
     fi
     if [[ "${GOTAVX2}" == "Y" ]]; then
 	echo "using AVX2 instructions"
-	FORCETARGET=" --enable-sse-assembly --enable-avx --enable-avx2  --disable-avx512  "
+	FORCETARGET=" --enable-sse-assembly-kernels --enable-avx-kernels --enable-avx2-kernels  --disable-avx512-kernels  "
 #	CFLAGS+=" -mmmx -msse -msse2 -msse3 -mssse3 -msse4.1 -msse4.2 -maes -mavx -mfma -mavx2 "
     fi
     if [[ "${GOTAVX512}" == "Y" ]]; then
 	echo "using AVX512 instructions"
-	FORCETARGET=" --disable-sse-assembly --enable-avx --enable-avx2  --enable-avx512 "
+	FORCETARGET=" --disable-sse-assembly-kernels --enable-avx-kernels --enable-avx2-kernels  --enable-avx512-kernels "
     fi
 fi #USE_HWOPT
 if [[ `${CC} -dM -E - < /dev/null 2> /dev/null | grep -c GNU` > 0 ]] ; then
@@ -262,6 +263,45 @@ fi
 if [[ !  -z "${BLASOPT}" ]]; then
     MYLINK+=" ${BLASOPT} "
 fi
+if [[ ! -z "${USE_OPENMP}" ]]; then
+    export OMPFLAGS="--enable-openmp"
+    if [ -f ${HOMEBREW_PREFIX}/opt/libomp/lib/libomp.a ]
+    then
+        MYFCFLAGS+="-I${HOMEBREW_PREFIX}/opt/libomp/include"
+        MYCFLAGS+="-I${HOMEBREW_PREFIX}/opt/libomp/include"
+        MYLINK+=" -L${HOMEBREW_PREFIX}/opt/libomp/lib -lomp"
+    fi
+fi
+if [[ ! -z "${ELPA_NVIDIA}" ]]; then
+# check if we have nvcc in the PATH    
+    if ! [ -x "$(command -v nvcc)" ]; then
+	echo
+	echo nvcc not found
+	echo nvcc is required for building Elpa with Nvidia GPUs
+	echo
+	exit 1
+    fi
+    GPUFLAGS=--enable-nvidia-gpu-kernels
+    if [[ ! -z "${GPU_ARCH}" ]]; then
+	export NVCC_APPEND_FLAGS=-arch=${GPU_ARCH}
+	GPUFLAGS+=" --with-NVIDIA-GPU-compute-capability=${GPU_ARCH} "
+#	if [[ "${GPU_ARCH}" == sm_80 ]]; then
+#	    GPUFLAGS+=" --enable-nvidia-sm80-gpu-kernels "
+#	fi
+    fi
+#    GPUFLAGS+=" --with-default-real-kernel=nvidia_gpu "
+    if [[ ! -z "${CUDA_ROOT}" ]]; then
+	GPUFLAGS+=" --with-cuda-path=${CUDA_ROOT} "
+    else
+	echo " "
+	echo "Please specify location of CUDA installation"
+	echo "by setting the env. variable CUDA_ROOT"
+	echo " "
+	exit 1
+	
+    fi
+    MYLDFLAGS+=-lstdc++
+fi
 echo MYFCFLAGS is $MYFCFLAGS
 echo MYCFLAGS is $MYCFLAGS
 echo 64ints is $sixty4_int
@@ -273,17 +313,21 @@ export LIBS="${MYLINK}"
 export    FC=$MPIF90
 export CC=$MPICC
 echo FC is $MPIF90 CC is $MPICC CXX is $MPICXX
+     set -x
  FC=$MPIF90 CC=$MPICC CXX=$MPICXX ../configure \
     $sixty4_int \
   CFLAGS="$MYCFLAGS" \
   FCFLAGS="$MYFCFLAGS" \
+  LDFLAGS="$MYLDFLAGS" \
   --enable-option-checking=fatal \
- --disable-openmp  \
  --disable-dependency-tracking \
  --disable-shared --enable-static  \
  --disable-c-tests \
+ --disable-detect-mpi-launcher \
      ${FORCETARGET} \
---prefix=${NWCHEM_TOP}/src/libext
+     ${GPUFLAGS} \
+     ${OMPFLAGS} \
+     --prefix=${NWCHEM_TOP}/src/libext  ||  { echo config libs failure; dump_build_env ; exit 1; }
 unset FORCETARGET
 unset LIBS
 unset FCFLAGS
@@ -293,7 +337,11 @@ unset SCALAPACK_LDFLAGS
 echo mpif90 is `which mpif90`
 echo MPIF90 is "$MPIF90"
 if [[ "$USE_MANUALCPP" == 1 ]]; then echo @@@@ MANUALCPP @@@; fi
-make FC=$MPIF90 CC=$MPICC CXX=$MPICXX -j4
+make FC=$MPIF90 CC=$MPICC CXX=$MPICXX install-libLTLIBRARIES -j4 || { echo make libs failure; dump_build_env ; exit 1; }
+make FC=$MPIF90 CC=$MPICC CXX=$MPICXX install-nobase_elpa_includeHEADERS || { echo make headers1 failure; dump_build_env ; exit 1; }
+make FC=$MPIF90 CC=$MPICC CXX=$MPICXX install-nobase_nodist_elpa_includeHEADERS || { echo make headers2 failure; dump_build_env ; exit 1; }
+make FC=$MPIF90 CC=$MPICC CXX=$MPICXX install-pkgconfigDATA || { echo make installpkg failure; dump_build_env ; exit 1; }
+make FC=$MPIF90 CC=$MPICC CXX=$MPICXX install-binPROGRAMS || { echo make instalbins failure; dump_build_env ; exit 1; }
 if [[ "$?" != "0" ]]; then
     echo " "
     echo "Elpa compilation failed"
@@ -302,7 +350,12 @@ if [[ "$?" != "0" ]]; then
     cat config.log
     exit 1
 fi
-    make V=0   install
+#    make V=0   install
 
-cp ${NWCHEM_TOP}/src/libext/lib/libelpa.a  ${NWCHEM_TOP}/src/libext/lib/libnwc_elpa.a
-cp -r ${NWCHEM_TOP}/src/libext/include/elpa-${SHORTVERSION}  ${NWCHEM_TOP}/src/libext/include/elpa
+if [[ ! -z "${USE_OPENMP}" ]]; then
+    ompsuffix="_openmp"
+fi
+cp ${NWCHEM_TOP}/src/libext/lib/libelpa${ompsuffix}.a  ${NWCHEM_TOP}/src/libext/lib/libnwc_elpa.a
+${NWCHEM_TOP}/src/libext/bin/elpa2_print_kernels${ompsuffix}
+cp -r ${NWCHEM_TOP}/src/libext/include/elpa${ompsuffix}-*  ${NWCHEM_TOP}/src/libext/include/elpa
+mkdir -p ${NWCHEM_TOP}/src/libext/lib/pkgconfig

src/libext/libext_utils/cmake.sh

@@ -1,7 +1,11 @@
+export CMAKE_VERSION_REQUIRED=3.26.0
+CMAKE_VER_REQ_MAJ=$(echo ${CMAKE_VERSION_REQUIRED}|cut -d . -f 1)
+CMAKE_VER_REQ_MIN=$(echo ${CMAKE_VERSION_REQUIRED}|cut -d . -f 2)
+
 get_cmake_release(){
     UNAME_S=$(uname -s)
     CPU=$(uname -m)
-    CMAKE_VER=3.26.0
+    CMAKE_VER=${CMAKE_VERSION_REQUIRED}
     orgdir=`pwd`
     cmake_instdir=$1
     echo "Parameter #1 is $1"
@@ -47,13 +51,15 @@ get_cmake_release(){
 }
 
 get_cmake_master(){
-    CMAKE_COMMIT=09dd52c9d2684e933a3e013abc4f6848cb1befbf
+    CMAKE_COMMIT=v3.26.6
     if [[ -f "cmake-$CMAKE_COMMIT.zip" ]]; then
 	echo "using existing"  "cmake-$CMAKE_COMMIT.zip" >> /tmp/cmake.log
     else
-	curl -L https://gitlab.kitware.com/cmake/cmake/-/archive/$CMAKE_COMMIT.zip -o cmake-$CMAKE_COMMIT.zip
+	curl -L https://github.com/Kitware/CMake/archive/refs/tags/$CMAKE_COMMIT.zip -o cmake-$CMAKE_COMMIT.zip
     fi
+    rm -rf CMake*
     unzip -n -q cmake-$CMAKE_COMMIT.zip
+    ln -sf CMake* cmake-$CMAKE_COMMIT
     mkdir -p  cmake-$CMAKE_COMMIT/build
     cd cmake-$CMAKE_COMMIT/build
     if [[ -x "$(command -v cmake)" ]]; then

src/libext/libext_utils/getfiles_utils.sh

@@ -0,0 +1,95 @@
+
+get_elpa(){
+    if [[ $# -eq 0 ]] ; then
+	elpa_shortv=2025.06.002
+    else
+	elpa_shortv=$1
+    fi
+    echo "Parameter #1 is $1"
+#    VERSION=new_release_${elpa_shortv}
+    VERSION=${elpa_shortv}
+    echo ELPA VERSION is $VERSION
+    if [ -f  elpa-${VERSION}.tar.gz ]; then
+	echo "using existing"  elpa-${VERSION}.tar.gz
+	return 0
+    else
+	rm -rf elpa*
+	ELPA_URL=("https://elpa.mpcdf.mpg.de/software/tarball-archive/Releases/${VERSION}/elpa-${VERSION}.tar.gz" \
+		      "https://web.archive.org/web/20260214233634/https://gitlab.mpcdf.mpg.de/elpa/elpa/-/archive/new_release_2025.06.002/elpa-new_release_2025.06.002.tar.gz")
+	for url in  "${ELPA_URL[@]}"
+	do
+	    echo using $url
+	    tries=1
+	    until [ "$tries" -ge 3 ]
+	    do
+		if [ "$tries" -gt 1 ]; then echo sleeping for 9s ;sleep 9; echo attempt no.  $tries ; fi
+		curl -L --progress-bar  $url -o elpa-${VERSION}.tar.gz
+		#echo checking gz integrity
+		if [ -f elpa-${VERSION}.tar.gz ]; then
+		    gzip -t elpa-${VERSION}.tar.gz >&  /dev/null
+		    if [ $? -eq 0 ]; then break ;  fi
+		fi
+		tries=$((tries+1))
+	    done
+	    gzip -t elpa-${VERSION}.tar.gz >&  /dev/null
+	    if [ $? -eq 0 ]; then return 0 ;  fi
+	done
+    fi
+    return 1
+}
+get_scalapack(){
+    if [[ $# -eq 0 ]] ; then
+	#	version=2.1.0
+ 	COMMIT=a23c2cdc6586c427686f6097ae66bb54ef693571
+	#COMMIT=bd1768b91262b4cdc7dd5f87b373b9b18eda4636
+	#COMMIT=b935167ca4d244735abc04a3cd4f6d56699702a0
+    else
+	COMMIT=$1
+    fi
+    echo "Parameter #1 is $1"
+    echo Scalapack commit is $COMMIT
+    rm -rf scalapack 
+    if [[ -f "scalapack-$COMMIT.tar.gz" ]]; then
+	echo "using existing"  "scalapack-$COMMIT.tar.gz"
+    else
+	echo "downloading"  "scalapack-$COMMIT.tar.gz"
+	rm -f scalapack-$COMMIT.tar.gz
+	tries=1
+	until [ "$tries" -ge 6 ]
+	do
+	    if [ "$tries" -gt 1 ]; then sleep 9; echo attempt no.  $tries ; fi
+	    curl -L https://github.com/Reference-ScaLAPACK/scalapack/archive/$COMMIT.tar.gz -o scalapack-$COMMIT.tar.gz
+	    # check tar.gz integrity
+	    gzip -t scalapack-$COMMIT.tar.gz >&  /dev/null
+	    if [ $? -eq 0 ]; then return 0 ;  fi
+	    tries=$((tries+1)) ;  done
+    fi
+    return 1
+}
+
+get_openblas()
+{
+    if [[ $# -eq 0 ]] ; then
+	VERSION=0.3.29
+    else
+	VERSION=$1
+    fi
+    echo "Parameter #1 is $1"
+    echo OpenBLAS VERSION is $VERSION
+    if [ -f  OpenBLAS-${VERSION}.tar.gz ]; then
+	echo "using existing"  OpenBLAS-${VERSION}.tar.gz
+    else
+	rm -rf OpenBLAS* openblas*
+	tries=1
+	until [ "$tries" -ge 6 ]
+	do
+	    if [ "$tries" -gt 1 ]; then sleep 9; echo attempt no.  $tries ; fi
+	    curl -L https://github.com/OpenMathLib/OpenBLAS/archive/v${VERSION}.tar.gz -o OpenBLAS-${VERSION}.tar.gz ;
+	    # check tar.gz integrity
+	    gzip -t OpenBLAS-${VERSION}.tar.gz >&  /dev/null
+	    if [ $? -eq 0 ]; then return 0 ;  fi
+	    tries=$((tries+1)) ;  done
+    fi
+    if [ $? -ne 0 ]; then echo  "openBLAS tarball not ready"; rm -f OpenBLAS-${VERSION}.tar.gz; return 1 ; fi
+    return 0
+}

src/libext/libxc/build_libxc.sh

@@ -78,7 +78,8 @@ fi
 CMAKE_VER_MAJ=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f1)
 CMAKE_VER_MIN=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f2)
 echo CMAKE_VER is ${CMAKE_VER_MAJ} ${CMAKE_VER_MIN}
-if ((CMAKE_VER_MAJ < 3)) || (((CMAKE_VER_MAJ == 3) && (CMAKE_VER_MIN < 24))); then
+echo CMAKE_VERSION_REQUIRED is ${CMAKE_VER_REQ_MAJ} ${CMAKE_VER_REQ_MIN}
+if ((CMAKE_VER_MAJ < ${CMAKE_VER_REQ_MAJ})) || (((CMAKE_VER_MAJ == ${CMAKE_VER_REQ_MAJ}) && (CMAKE_VER_MIN < ${CMAKE_VER_REQ_MIN}))); then
     cmake_instdir=../libext_utils
     get_cmake_release  $cmake_instdir
     echo PATH is $PATH

src/libext/openblas/build_openblas.sh

@@ -1,23 +1,9 @@
 #!/usr/bin/env bash
 #set -v
 arch=`uname -m`
+source ../libext_utils/getfiles_utils.sh
 VERSION=0.3.29
-#COMMIT=974acb39ff86121a5a94be4853f58bd728b56b81
-BRANCH=develop
-if [ -f  OpenBLAS-${VERSION}.tar.gz ]; then
-    echo "using existing"  OpenBLAS-${VERSION}.tar.gz
-else
-    rm -rf OpenBLAS*
-    tries=1 ; until [ "$tries" -ge 6 ] ; do
-		  if [ "$tries" -gt 1 ]; then sleep 9; echo attempt no.  $tries ; fi
-		  curl -L https://github.com/xianyi/OpenBLAS/archive/v${VERSION}.tar.gz -o OpenBLAS-${VERSION}.tar.gz ;
-		  # check tar.gz integrity
-		  gzip -t OpenBLAS-${VERSION}.tar.gz >&  /dev/null
-		  if [ $? -eq 0 ]; then break ;  fi
-		  tries=$((tries+1)) ;  done
-fi
-gzip -t OpenBLAS-${VERSION}.tar.gz >&  /dev/null
-if [ $? -ne 0 ]; then echo  "openBLAS tarball not ready"; rm -f OpenBLAS-${VERSION}.tar.gz; exit 1 ; fi
+get_openblas $VERSION
 tar xzf OpenBLAS-${VERSION}.tar.gz
 ln -sf OpenBLAS-${VERSION} OpenBLAS
 cd OpenBLAS
@@ -98,7 +84,8 @@ else
     if [[ -n "${USE_DYNAMIC_ARCH}" ]] || [[ "${USE_HWOPT}" == "n" ]]; then
 	if [[ "$arch" == "x86_64" ]]; then
 	    echo   "not cross compiling, therefore using DYNAMIC_ARCH "
-	    FORCETARGET+="DYNAMIC_ARCH=1 DYNAMIC_OLDER=1"
+	    FORCETARGET+="DYNAMIC_ARCH=1"
+#	    FORCETARGET+="DYNAMIC_ARCH=1 DYNAMIC_OLDER=1"
 	fi
     fi
 fi
@@ -249,7 +236,7 @@ if [[  ! -z "${USE_OPENMP}" ]]; then
 fi
 GOTFREEBSD=$(uname -o 2>&1|awk ' /FreeBSD/ {print "1";exit}')
 MYMAKE=make
-MAKEJ="MAKE_NB_JOBS=2"
+MAKEJ="MAKE_NB_JOBS=4"
 MAKE_MAJOR=$(make --version 2>& 1|head -1| cut -d " " -f 3 |cut -d .  -f 1)
 MAKE_MINOR=$(make --version 2>& 1|head -1| cut -d " " -f 3 |cut -d .  -f 2)
 if [[ ${MAKE_MAJOR} -ge 4 ]] && [[ ${MAKE_MINOR} -ge 4 ]]; then

src/libext/scalapack/build_scalapa.sh

@@ -10,7 +10,7 @@ source ../libext_utils/cmake.sh
 cd $myscalapwd
 
 if [[ -z "${MPIF90}" ]]; then
-if [[ "$FC" = "ftn"  ]] ; then
+if [[ "$FC" = "_ftn"  ]] ; then
     MPIF90="ftn"
     MPICC="cc"
 else
@@ -70,7 +70,7 @@ if [[ ! -z "${BUILD_MPICH}" ]]; then
     fi
     echo LDFLAGS for hwloc is $LDFLAGS
 fi
-if [[ "$FC" = "ftn"  ]] || [[ ! -z "$USE_CMAKE_MASTER" ]] ; then
+if [[ ! -z "$USE_CMAKE_MASTER" ]] ; then
     get_cmake_master
 else
 if [[ -z "${CMAKE}" ]]; then
@@ -95,7 +95,8 @@ fi
 CMAKE_VER_MAJ=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f1)
 CMAKE_VER_MIN=$(${CMAKE} --version|cut -d " " -f 3|head -1|cut -d. -f2)
 echo CMAKE_VER is ${CMAKE_VER_MAJ} ${CMAKE_VER_MIN}
-if ((CMAKE_VER_MAJ < 3)) || (((CMAKE_VER_MAJ == 3) && (CMAKE_VER_MIN < 24))); then
+echo CMAKE_VERSION_REQUIRED is ${CMAKE_VER_REQ_MAJ} ${CMAKE_VER_REQ_MIN}
+if ((CMAKE_VER_MAJ < ${CMAKE_VER_REQ_MAJ})) || (((CMAKE_VER_MAJ == ${CMAKE_VER_REQ_MAJ}) && (CMAKE_VER_MIN < ${CMAKE_VER_REQ_MIN}))); then
     cmake_instdir=../libext_utils
     get_cmake_release $cmake_instdir
     echo PATH is $PATH
@@ -138,31 +139,11 @@ fi
 if [[ ! -z "$BUILD_OPENBLAS"   ]] ; then
     BLASOPT="-L`pwd`/../lib -lnwc_openblas"
 fi
-#git clone https://github.com/scibuilder/scalapack.git
-#svn co --non-interactive --trust-server-cert https://icl.utk.edu/svn/scalapack-dev/scalapack/trunk/ scalapack
-VERSION=2.1.0
-#curl -L https://github.com/Reference-ScaLAPACK/scalapack/archive/v${VERSION}.tar.gz -o scalapack.tgz
-#COMMIT=bc6cad585362aa58e05186bb85d4b619080c45a9
-#COMMIT=ea5d20668a6b8bbee645b7ffe44623c623969d33
-COMMIT=5bad7487f496c811192334640ce4d3fc5f88144b
-COMMIT=782e739f8eb0e7f4d51ad7dd23fc1d03dc99d240
+source ../libext_utils/getfiles_utils.sh
 COMMIT=a23c2cdc6586c427686f6097ae66bb54ef693571
 #COMMIT=bd1768b91262b4cdc7dd5f87b373b9b18eda4636
 #COMMIT=b935167ca4d244735abc04a3cd4f6d56699702a0
-rm -rf scalapack 
-if [[ -f "scalapack-$COMMIT.tar.gz" ]]; then
-    echo "using existing"  "scalapack-$COMMIT.tar.gz"
-else
-    echo "downloading"  "scalapack-$COMMIT.tar.gz"
-    rm -f scalapack-$COMMIT.tar.gz
-    tries=1 ; until [ "$tries" -ge 6 ] ; do
-		  if [ "$tries" -gt 1 ]; then sleep 9; echo attempt no.  $tries ; fi
-		  curl -L https://github.com/Reference-ScaLAPACK/scalapack/archive/$COMMIT.tar.gz -o scalapack-$COMMIT.tar.gz
-		  # check tar.gz integrity
-		  gzip -t scalapack-$COMMIT.tar.gz >&  /dev/null
-		  if [ $? -eq 0 ]; then break ;  fi
-		  tries=$((tries+1)) ;  done
-fi
+get_scalapack $COMMIT
 tar xzf scalapack-$COMMIT.tar.gz
 ln -sf scalapack-*$COMMIT scalapack
 #ln -sf scalapack-${VERSION} scalapack
@@ -240,10 +221,18 @@ if [[  -z "$MPICH_FC"   ]] ; then
     export MPICH_FC="$FC"
 fi
 echo MPICH_FC is "$MPICH_FC"
+if [[  -z "$OMPI_FC"   ]] ; then
+    export OMPI_FC="$FC"
+fi
+echo OMPI_FC is "$OMPI_FC"
 if [[  -z "$MPICH_CC"   ]] ; then
     export MPICH_CC="$CC"
 fi
 echo MPICH_CC is "$MPICH_CC"
+if [[  -z "$OMPI_CC"   ]] ; then
+    export OMPI_CC="$CC"
+fi
+echo OMPI_CC is "$OMPI_CC"
 echo $(${MPICC} -show)
 #Intel MPI
 if [[  -z "$I_MPI_F90"   ]] ; then
@@ -268,10 +257,13 @@ if [[  -z "$PE_ENV"   ]] ; then
 fi
 #fix for clang 12 error in implicit-function-declaration
 GOTCLANG=$( "$MPICC" -dM -E - </dev/null 2> /dev/null |grep __clang__|head -1|cut -c19)
+GOTNVC=$( "$MPICC" -dM -E - </dev/null 2> /dev/null |grep __NVCOMPILER\  |cut -d " " -f 3)
 if [[ `${MPICC} -dM -E - < /dev/null 2> /dev/null | grep -c GNU` > 0 ]] ; then
+    if [[ ${GOTNVC} != 1 ]]; then
     let GCCVERSIONGT12=$(expr `${MPICC} -dumpversion | cut -f1 -d.` \> 12)
     let GCCVERSIONGT13=$(expr `${MPICC} -dumpversion | cut -f1 -d.` \> 13)
     let GCCVERSIONGT14=$(expr `${MPICC} -dumpversion | cut -f1 -d.` \> 14)
+    fi
 fi
 
 if [[ ${GOTCLANG} == "1" ]] || [[ ${GCCVERSIONGT12} == "1" ]]  ; then

src/mp2_grad/mp2_make_tuhf.F

@@ -50,6 +50,7 @@ c
       parameter (eps=1d-14)
 c
       logical otdebug
+      logical dumpabij
       logical do_mirra,do_mirrb,do_mirr
       integer g_ta_mirr,g_tb_mirr,g_ta_org,g_tb_org
       logical util_mirrmat,tuhfreplicated
@@ -58,12 +59,14 @@ c
       external util_mirrmem
 c
       integer l_ia_uv, k_ia_uv, l_tmp, k_tmp, l_ia_jb, k_ia_jb
+      character*255 dumpabij_fn
       double precision ddot
       external ddot
 c
 #include "bitops.fh"
 c
       otdebug = util_print('tdebug',print_debug)
+      dumpabij = util_print('dumpabij',print_never)
 c     
       eaa=0.0d0
       eab=0.0d0
@@ -133,6 +136,11 @@ c
       endif
       tunitptra=1
       tunitptrb=1
+      if(dumpabij) then
+         call util_file_name('dumpabij',.false.,.true.,dumpabij_fn)
+         write(6,*)' dumpabij filename', dumpabij_fn
+         open(90, status='unknown', form='formatted', file=dumpabij_fn)
+      endif
       do i=oseg_lo,oseg_hi 
          symi=irs_a(i)
          call ga_zero(g_t_a)
@@ -260,11 +268,18 @@ c
      $                    -eval_a(i)-eval_a(j)
                      t=dbl_mb(counta+k_ia_jb+j-sym_lo_oa(symj))/denom
                      if(abs(t).lt.eps)t=0.0d0
+                  if(dumpabij.and.abs(t).gt.0d0)
+     D                    write(90,9012)
+     A            nbf*(nbf*(nbf*(i-1)+(j-1))+a_a-1)+b,
+     I                    i, j, a_a, b,
+     N                    -0.25d0*t*t*denom*denom,
+     D                    denom
 c
                      eaa=eaa+t*t*denom
                      dbl_mb(counta+k_ia_jb+j-sym_lo_oa(symj))=t
                   enddo
                   counta=counta+(sym_hi_oa(symj)-sym_lo_oa(symj)+1)
+ 9012               format(1x,i12,4i6,2(1pe15.7))
                endif
             enddo
             if(task.eq.'gradient' .and. counta.gt.0)then
@@ -294,6 +309,12 @@ c
      $                    -eval_a(i)-eval_b(j)
                      t=dbl_mb(countb+k_ia_jb+j-sym_lo_ob(symj))/denom
                      if(abs(t).lt.eps)t=0.0d0
+                     if(dumpabij.and.abs(t).gt.0d0)
+     D                    write(90,9012)
+     N                    nbf*(nbf*(nbf*(i-1)+(j-1))+a_a-1)+b,
+     I                    i, j, a_a, b,
+     N                    -0.5d0*t*t*denom*denom,
+     D                    denom
                      eab=eab+t*t*denom
                      dbl_mb(countb+k_ia_jb+j-sym_lo_ob(symj))=t
                   enddo

src/nwc_columbus/aoints/GNUmakefile

@@ -3,8 +3,10 @@
 
 
 include ../../config/makefile.h
-OBJ_OPTIMIZE =      nwc_sym_mod.o \
-     		    dft_aoints.o       wrt_dft_aoints.o int_1e_sifs.o \
+MAKEFLAGS = -j 1 --no-print-directory
+OBJ =      nwc_sym_mod.o 
+OBJ_OPTIMIZE =       \
+		    dft_aoints.o       wrt_dft_aoints.o int_1e_sifs.o \
                     int_so_sifs.o      int_2e_sifs.o    int_2e_sifs_a.o \
                     sifs_2e_task.o     int_2e_sifs_b.o  nadct_trace.o \
                     rdhcid.o           egrad_trace.o    rd1mat.o sif2ga.o \

src/nwchem.F

@@ -30,7 +30,7 @@ c
 c  $Id$
 c
 c  ======================================================================================================
-C> \mainpage Northwest Computational Chemistry Package (NWChem) 7.3.0
+C> \mainpage Northwest Computational Chemistry Package (NWChem) 7.3.1
 C>
 C> NWChem is an open-source computational chemistry package distributed under the terms of
 C> the Educational Community License (ECL) 2.0
@@ -143,7 +143,9 @@ c     Initialize local memory allocator & global array tools
 C
       call ga_initialize_ltd(ma_sizeof(mt_dbl,global,mt_byte))
       ! this must happen after GA and before MA!
+#if defined(USE_CUDA_AFFINITY) || defined(USE_OPENACC_AFFINITY)
       call util_setup_gpu_affinity()
+#endif
       if ( ga_uses_ma() ) then
         if (.not.ma_init(mt_dbl, stack, heap+global))
      &      call errquit('nwchem.F: ma_init failed (ga_uses_ma=T)',911,
@@ -171,6 +173,7 @@ C
  99   format(2x,'NWChem w/ OpenMP: maximum threads = ',i4)
       do_printaff=.false.
 #if defined(USE_OPENMP)
+      call util_setmaxthreads()
       !$omp parallel
       !$omp master
       if (ga_nodeid().eq.0) write(luout,99) omp_get_max_threads()
@@ -472,7 +475,7 @@ c
       write(LuOut,*)
       write(LuOut,*)
       call util_print_centered(LuOut,
-     $     'Northwest Computational Chemistry Package (NWChem) 7.3.0',
+     $     'Northwest Computational Chemistry Package (NWChem) 7.3.1',
      $     40, .true.)
       write(LuOut,*)
       write(LuOut,*)
@@ -568,7 +571,7 @@ C
 C  Print version information to ecce.out file
 C
       compiled(inp_strlen(compiled)+1:inp_strlen(compiled)+25) = 
-     $         ' Version 7.3.0'
+     $         ' Version 7.3.1'
       call ecce_print_version(compiled(1:inp_strlen(compiled)))
 C
       If ( Inp_StrLen(BatchID) .gt. 0) then

src/nwdft/dftgrad/dft_gradients.F

@@ -622,9 +622,7 @@ c
      &        call errquit('dft_gradients: could not alloc j hessian',
      &        1, MA_ERR)
          call dfill(9*nat*nat, 0.0d0, dbl_mb(k_hess), 1)
-         if (bgj_print() .gt. 0)
-     &        write(LuOut,*)'*** In dft_gradients: calling j_hessian'
-         call j_hessian(iga_dens, log_mb(k_act), nactive,
+         call j_hessian(rtdb,iga_dens, log_mb(k_act), nactive,
      &        dbl_mb(k_hess))
          status = MA_pop_stack(l_hess)
          if (.not.status) call
@@ -661,8 +659,6 @@ c!!!         call int_app_set_no_texas(rtdb)
 c
 c        Allocate and initialize temp GA's for RHS
 c
-         if (bgj_print() .gt. 0)
-     &        write(*,*)'*** j cpks rhs test: nactive =',nactive
          if (nat.gt.100)
      &        call errquit('dft_gradients: dimension error in test',0,
      &       UNKNOWN_ERR)
@@ -680,9 +676,7 @@ c
             endif
          enddo
 
-         if (bgj_print() .gt. 0)
-     &        write(LuOut,*)'*** In dft_gradients: calling j_cpks_rhs'
-         call j_cpks_rhs(iga_dens, log_mb(k_act), nactive, g_rhs)
+         call j_cpks_rhs(rtdb,iga_dens, log_mb(k_act), nactive, g_rhs)
 
          do i = 1, nat
             if (log_mb(k_act+i-1)) then

src/nwdft/dftgrad/dftg_getxc.F

@@ -95,7 +95,7 @@ c
       endif
 c
 c        Compute via direct numerical quadrature.
-         if(.not.grid_written) then
+         if(.not.grid_written.and.(.not.largenode)) then
 c
 c        check if we can get the grid from a file
 c         
@@ -131,8 +131,6 @@ c
      &         call errquit('dftg_getxc: could not allocate xc hessian',
      &                      1, MA_ERR)
             call dfill(9*natom*natom, 0.0d0, dbl_mb(k_hess), 1)
-            if (bgj_print() .gt. 0)
-     &           write(LuOut,*)'*** In dftg_getxc: calling xc_hessian'
             call xc_hessian(geom,rtdb,nbf_ao,ao_bas_han,
      .           iga_dens, oactive, nactive, dbl_mb(k_hess),
      .           scftype)
@@ -152,8 +150,6 @@ c
 c
 c           Allocate and initialize temp GA's for RHS
 c
-            if (bgj_print() .gt. 0)
-     &           write(*,*)'*** xc cpks rhs test: nactive =',nactive
             if (ipol*3*natom.gt.100)
      &           call errquit('dftg_getxc: dimension error in test',0,
      &       INPUT_ERR)
@@ -193,8 +189,6 @@ c
              endif
             enddo
 
-            if (bgj_print() .gt. 0)
-     &           write(LuOut,*)'*** In dftg_getxc: calling xc_cpks_rhs'
             call xc_cpks_rhs(geom,rtdb,nbf_ao,ao_bas_han,
      .           iga_dens, oactive, nactive, g_rhs,
      .           scftype)

src/nwdft/dftgrad/dftg_griddo.F

@@ -136,9 +136,6 @@ c
       Ec = 0.d0
       Ex = 0.d0
       oprint = util_print('quadrature', print_high)
-      if (bgj_print() .gt. 0)
-     &    write(*,*)'*** dftg_gridv0a: logicals ',
-     &    do_force,do_hess,do_cpks_r
 c
 c     Determine the XC energy functionals.
 c

src/nwdft/dftgrad/dftg_gridv0a.F

@@ -169,9 +169,6 @@ c
       me = ga_nodeid()
       Ec = 0.d0
       Ex = 0.d0
-      if (bgj_print() .gt. 0)
-     &    write(*,*)'*** dftg_gridv0a: logicals ',
-     &    do_force,do_hess,do_cpks_r
 c
 c     Determine the XC energy functionals.
 c

src/nwdft/dftgrad/j_nucder_gen.F

@@ -7,41 +7,43 @@ c     BGJ - 9/98
 c
 c     $Id$
 c
-      subroutine j_hessian(g_dens, oactive, nactive, hess)
+      subroutine j_hessian(rtdb,g_dens, oactive, nactive, hess)
 c
 c     Wrapper routine for J hessian without fitting
 c
       implicit none
 c
+      integer rtdb
       integer g_dens            ! [input] GA handle for total density
       integer nactive           ! [input] number of active atoms
       logical oactive(nactive)  ! [input] key to active atoms
       double precision hess(*)  ! [output] J hessian matrix
 c
-      call j_deriv_gen(g_dens, oactive, nactive, hess, 0, 2)
+      call j_deriv_gen(rtdb,g_dens, oactive, nactive, hess, 0, 2)
 c
       return
       end
 c
 c
-      subroutine j_cpks_rhs(g_dens, oactive, nactive, g_rhs)
+      subroutine j_cpks_rhs(rtdb,g_dens, oactive, nactive, g_rhs)
 c
 c     Wrapper routine for J CPKS RHS without fitting
 c
       implicit none
 c
+      integer rtdb
       integer g_dens            ! [input] handle to total density
       integer nactive           ! [input] number of active atoms
       logical oactive(nactive)  ! [input] key to active atoms
       integer g_rhs(*)          ! [output] J CPKS RHS matrices
 c
-      call j_deriv_gen(g_dens, oactive, nactive, 0d0, g_rhs, 3)
+      call j_deriv_gen(rtdb,g_dens, oactive, nactive, 0d0, g_rhs, 3)
 c
       return
       end
 c
 c
-      subroutine j_deriv_gen(g_dens, oactive, nactive, hess, g_rhs,
+      subroutine j_deriv_gen(rtdb,g_dens, oactive, nactive, hess, g_rhs,
      &                       calc_type)
 c
 c     $Id$
@@ -65,6 +67,7 @@ c
       external ga_create_atom_blocked
 c!!! End BGJ temp test code
 c
+      integer rtdb
       integer g_dens            ! [input] handle to total density
       integer nactive           ! [input] number of active atoms
       logical oactive(nactive)  ! [input] key to active atoms
@@ -78,7 +81,6 @@ c
 c     Local declarations
 c
       logical do_hess, do_cpks_r
-      integer rtdb
       integer natoms
       integer l_d_ija, k_d_ija, l_d_kla, k_d_kla, ldim1, ldim2,
      &        l_db, k_db, l_scr, k_scr, lscr, l_buf, k_buf, lbuf,
@@ -92,16 +94,11 @@ c !!! Start BGJ temp test code
       integer l_ti_db, k_ti_db
 c !!! End BGJ temp test code
 c
-      if (bgj_print() .gt. 0)
-     &     write(*,*)'--------Entered j_deriv_gen-------------'
       do_hess = calc_type .eq. 2
       do_cpks_r = calc_type .eq. 3
       if (.not. (do_hess .or. do_cpks_r))
      &     call errquit('j_deriv_gen: illegal calculation type',0,
      &       INPUT_ERR)
-      if (bgj_print() .gt. 0)
-     &     write(*,*)'logicals ', do_hess, do_cpks_r
-      rtdb = bgj_get_rtdb_handle()
       if (.not. geom_ncent(geom, natoms))
      &   call errquit('j_deriv_gen: geom_ncent failed',1, GEOM_ERR)
 #if 0
@@ -147,8 +144,6 @@ c     of dftg_cdfit_gen (allocation is set up properly in j_nucder_alloc)
 c     The s34 array is used for the fitting coefficients in the
 c     fitted case (allocation is set up properly in j_nucder_alloc)
 c
-         if (bgj_print().gt.0) write(*,*)
-     &      ': what about oactive for fitted J?'
          call dftg_cdfit_gen(geom,AO_bas_han, CD_bas_han,
      &                       nbf_cd, natoms, tol2e,
      &                       dbl_mb(k_scr), lscr, dbl_mb(k_buf), lbuf,
@@ -243,70 +238,6 @@ c
          enddo
       endif
 c
-c !!! BGJ test !!!
-c All the code from here to the next !!! BGJ test !!! can be removed once
-c the code is stable
-c
-      if (bgj_print() .gt. 0) then
-         if (do_hess) then
-c
-            write(LuOut,*) 'j_deriv_gen: j hessian'
-            call output(hess,1,3*natoms,1,3*natoms,3*natoms,3*natoms,1)
-c
-c     TI check
-c
-            if (.not. MA_push_get(MT_DBL, natoms*natoms*9, 'ti_hess',
-     &           l_ti_hess, k_ti_hess))
-     &           call errquit('j_nucder_gen: could not alloc ti_hess',1,
-     &       MA_ERR)
-            call dcopy(natoms*natoms*9, hess, 1, dbl_mb(k_ti_hess), 1)
-            call hessti(dbl_mb(k_ti_hess),natoms)
-            if (.not. MA_pop_stack(l_ti_hess))
-     &           call errquit('j_nucder_gen: could not pop ti_hess',1,
-     &       MA_ERR)
-c
-         else if (do_cpks_r) then
-c
-            do i = 1, natoms
-               if (oactive(i)) then
-                  do j = 1, 3
-                     write(*,*)'*** j cpks rhs matrix',j,i
-                     call ga_print(g_rhs(j,i))
-                  enddo
-               endif
-            enddo
-c
-c     TI check
-c
-            g_ti(1) = ga_create_atom_blocked (geom, ao_bas_han,
-     &           'TI check x')
-            g_ti(2) = ga_create_atom_blocked (geom, ao_bas_han,
-     &           'TI check y')
-            g_ti(3) = ga_create_atom_blocked (geom, ao_bas_han,
-     &           'TI check z')
-            do j = 1, 3
-               call ga_zero(g_ti(j))
-            enddo
-            do i = 1, natoms
-               if (oactive(i)) then
-                  do j = 1, 3
-                     call ga_add(1d0, g_rhs(j,i), 1d0, g_ti(j), g_ti(j))
-                  enddo
-               endif
-            enddo
-            do j = 1, 3
-               write(*,*)'*** TI check',j
-               call ga_print(g_ti(j))
-               if (.not.ga_destroy(g_ti(j))) then
-                  call errquit('j_deriv_gen: could not destroy TI',j,
-     &       GA_ERR)
-               endif
-            enddo
-c
-         endif
-         write(*,*)'--------Leaving j_deriv_gen-------------'
-      endif
-c !!! BGJ test !!!
 c
       return
       end

src/nwdft/grid/grid_rtrunc.F

@@ -246,8 +246,10 @@ c
       if (ilo .eq. 0) then
          dexpo=sqrt(alpha)*r
          if(dexpo.gt.40d0) then
+            valexp=0d0
             value=0d0
          else
+            valexp=exp(-alpha*r*r)
 #if defined(WIN32) ||defined(LINUX)
             value = 0.5d0*sqrt(4.0d0*atan(1.0d0)/alpha)*
      $           derfc(dexpo)
@@ -256,13 +258,14 @@ c
      $           erfc(dexpo)
 #endif
          endif
-      endif
-      if(alpha*r*r.gt.500d0) then
-         valexp=0d0
-         value =0d0
       else
-         valexp=exp(-alpha*r*r)
-         value =valexp/(2.0d0*alpha)
+        if(alpha*r*r.gt.500d0) then
+           valexp=0d0
+           value =0d0
+        else
+           valexp=exp(-alpha*r*r)
+           value =valexp/(2.0d0*alpha)
+        endif
       endif
 c
       do i = ilo+2,k,2

src/nwdft/input_dft/dft_rdinput.F

@@ -1158,7 +1158,7 @@ c
         CDFIT=bas_rtdb_load(rtdb, geom, CD_bas_han, 'cd basis')
       endif
       if (CDFIT)then
-         call int_init(rtdb, 1, cd_bas_han)         
+         call int_init(rtdb,2,(/ao_bas_han,cd_bas_han/))
          if (.not. int_norm_2c(rtdb,CD_bas_han))
      &      call errquit('rdinput: int_norm_2c failed', 208, BASIS_ERR)
          call int_terminate()

src/nwdft/lr_tddft/tddft_analysis.F

@@ -1882,7 +1882,7 @@ c
       if (.not.rtdb_put(rtdb,'tddft:energy',mt_dbl,1,energy))
      1  call errquit('tddft_analysis: failed to put tddft energy',0,
      &       RTDB_ERR)
-      if(ga_nodeid().eq.0) write(6,*) ' stored tddft:energy ',energy
+      if(ga_nodeid().eq.0) write(luout,*) ' stored tddft:energy ',energy
 c     storing all calculated excited state energies on RTDB
       if(.not.rtdb_put(rtdb,'tddft:energy-all',mt_dbl,nroots,
      $                 apbval(1:nroots)))