From ad55d487e6927fbca12895d35d90e75e09bc7e07 Mon Sep 17 00:00:00 2001
From: edoapra <edoardo.apra@gmail.com>
Date: Fri, 22 Nov 2024 17:14:07 -0800
Subject: [PATCH] fix for https://github.com/nwchemgit/nwchem/issues/1049

---
 QA/doqmtests.mpi              |  2 ++
 src/nwdft/util/dft_genutils.F | 16 ++++++++++++++--
 src/task/task_bsse.F          | 13 ++++++++++++-
 3 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/QA/doqmtests.mpi b/QA/doqmtests.mpi
index 93cae666b6..1fce9f11c3 100755
--- a/QA/doqmtests.mpi
+++ b/QA/doqmtests.mpi
@@ -728,6 +728,8 @@ let "myexit+=$?"
 let "myexit+=$?"
 ./runtests.mpi.unix procs $np p2ta-vem
 let "myexit+=$?"
+./runtests.mpi.unix procs $np bsse_sodft
+let "myexit+=$?"
 if [[ "$what" != "flaky" ]]; then
     echo
     echo  "the number of failed tests is" $myexit
diff --git a/src/nwdft/util/dft_genutils.F b/src/nwdft/util/dft_genutils.F
index f84b8011e9..d8c1837d5c 100644
--- a/src/nwdft/util/dft_genutils.F
+++ b/src/nwdft/util/dft_genutils.F
@@ -42,7 +42,7 @@ C
      $     pname//" only valid for double precision matricies", 0, 0)
       
       if (n1 .ne. n2)
-     $     call errquit (pname//"n1 must equal n2")
+     $     call errquit (pname//"n1 must equal n2",0,0)
 C
 C     Do all work on processor 0.
 C     
@@ -131,7 +131,7 @@ C
      $     pname//" only valid for double precision  matricies", 0, 0)
 C      
       if (n1 .ne. n2)
-     $     call errquit (pname//"n1 must equal n2")
+     $     call errquit (pname//"n1 must equal n2",0,0)
 C
 C     Do all work on processor 0.
 C     
@@ -238,3 +238,15 @@ c
       return
       end
 c $Id$
+      subroutine dft_forceatomic(rtdb)
+      implicit none
+      integer rtdb
+#include "cdft.fh"
+#include "rtdb.fh"
+#include "errquit.fh"
+      movecs_in='atomic'
+      if (.not. rtdb_cput(rtdb, 'dft:input vectors', 1, movecs_in))
+     $    call errquit(' dft_forceatomic: rtdb_cput',0,RTDB_ERR)
+      write(6,*) ' force mos in ',movecs_in
+      return
+      end
diff --git a/src/task/task_bsse.F b/src/task/task_bsse.F
index 073eb140c8..4fe66d7843 100644
--- a/src/task/task_bsse.F
+++ b/src/task/task_bsse.F
@@ -434,6 +434,7 @@ c
 c:    multiplicity
 c:      density methods
       if ( theory(1:lentheo).eq.'dft' .or.
+     $     theory(1:lentheo).eq.'sodft'.or.
      $     theory(1:lentheo).eq.'tddft') then
         if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, mult))
      $    call errquit('bsse_param: rtdb_put of mult failed',
@@ -441,6 +442,7 @@ c:      density methods
 c:      wavefuntion methods
 
       elseif( theory(1:lentheo).ne.'dft' .and.
+     $        theory(1:lentheo).ne.'sodft'.and.
      $        theory(1:lentheo).ne.'tddft') then
         if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, mult-1))
      $    call errquit('bsse_param: rtdb_put of nopen failed',
@@ -905,6 +907,7 @@ c: take supermolecule total charge
 c: multiplicity
 c:    density methods
         if ( theory(1:lentheo).eq.'dft' .or.
+     $       theory(1:lentheo).eq.'sodft'.or.
      $       theory(1:lentheo).eq.'tddft') then
           if (.not. rtdb_get(rtdb, 'dft:mult', mt_int, 1, m_spr))
      $      call errquit('bsse_energy: rtdb_get of mult failed',
@@ -912,6 +915,7 @@ c:    density methods
 c:   wavefuntion methods
 
         elseif( theory(1:lentheo).ne.'dft' .and.
+     $          theory(1:lentheo).ne.'sodft'.and.
      $          theory(1:lentheo).ne.'tddft') then
           if (.not. rtdb_get(rtdb,'scf:nopen', MT_INT, 1, m_spr))
      $      call errquit('bsse_energy: rtdb_get of nopen failed',
@@ -921,7 +925,11 @@ c:   wavefuntion methods
 
 c
 c: name of the original movecs
-        if(theory(1:lentheo).eq.'dft' .or.
+        if(theory(1:lentheo).eq.'sodft') then
+          vec_dbo= 'dft:output vectors'
+          vec_dbi= 'atomic'
+          call dft_forceatomic(rtdb)
+        elseif(theory(1:lentheo).eq.'dft' .or.
      $     theory(1:lentheo).eq.'tddft') then
           vec_dbo= 'dft:output vectors'
           vec_dbi= 'dft:input vectors'
@@ -966,6 +974,7 @@ c
         call bsse_param(rtdb, mmon(j), qmon(j), j_mon_name, input(i),
      $                  theory)
 c
+        if(theory(1:lentheo).eq.'sodft') call dft_forceatomic(rtdb)
 c:evaluate energy
         if (.not. task_energy_doit(rtdb,theory, mon_energy(i)))
      $    call 
@@ -1028,6 +1037,7 @@ c: multiplicity
 c:    density methods
 
         if ( theory(1:lentheo).eq.'dft' .or.
+     $        theory(1:lentheo).eq.'sodft'.or.
      $       theory(1:lentheo).eq.'tddft') then
 
           if (.not. rtdb_put(rtdb, 'dft:mult', mt_int, 1, m_spr))
@@ -1037,6 +1047,7 @@ c:    density methods
 c:   wavefuntion methods
 
         elseif( theory(1:lentheo).ne.'dft' .and.
+     $         theory(1:lentheo).ne.'sodft'.and.
      $          theory(1:lentheo).ne.'tddft') then
 
           if (.not. rtdb_put(rtdb,'scf:nopen', MT_INT, 1, m_spr))