Energy = 1e+30

[mchatzikos]

This is another edge case, focusing on the 2s orbital of He-like Zn.

The attached zip file contains the script (and isodata file) to reproduce the problem. It uses an MR that consists of the ground state and the 2s orbital without excitations. rmcdhf_csfg_mpi reports the total energy of the 0+ levels as 1e+30, viz:

Level  1    Energy =  1.000000000000D+30    Weight =  2.00000D-01
Level  2    Energy =  1.000000000000D+30    Weight =  2.00000D-01
Level  1    Energy = -5.640448871352D+02    Weight =  6.00000D-01

If I exclude the ground state, things look normal:

Level  1    Energy = -5.636745059204D+02    Weight =  2.50000D-01
Level  1    Energy = -5.640443034698D+02    Weight =  7.50000D-01

I've added some printing statements, namely:

diff --git a/srcg/appl/rmcdhf90_csfg_mpi/maneig_csfg.f90 b/srcg/appl/rmcdhf90_csfg_mpi/maneig_csfg.f90
index 1afbc8b..c0e166b 100644
--- a/srcg/appl/rmcdhf90_csfg_mpi/maneig_csfg.f90
+++ b/srcg/appl/rmcdhf90_csfg_mpi/maneig_csfg.f90
@@ -227,6 +227,14 @@
 !
       NEND = NCF*NVEX
       DO J = 1, NEVBLK(JBLOCK)
+         WRITE(*,*) 'J=', J, 'NEVBLK(JBLOCK)=', NEVBLK(JBLOCK), &
+                 'NEND=', NEND, 'NCMINPAT=', NCMINPAT, &
+                 'J + NCMINPAT=', J + NCMINPAT, &
+                 'ICCMIN(J + NCMINPAT)=', ICCMIN(J + NCMINPAT), &
+                 'NEND + ICCMIN(J + NCMINPAT)=', &
+                 NEND + ICCMIN(J + NCMINPAT), &
+                 'WORK(NEND + ICCMIN(J + NCMINPAT))=', &
+                 WORK(NEND + ICCMIN(J + NCMINPAT))
          EVAL(NCMINPAT+J) = WORK(NEND + ICCMIN(J + NCMINPAT))
          CALL DCOPY (NCF, WORK(NCF*(ICCMIN(J + NCMINPAT) - 1) + 1), 1, EVEC(&
             NEVECPAT+NCF*(J-1)+1), 1)
diff --git a/srcg/appl/rmcdhf90_csfg_mpi/newcompi.f90 b/srcg/appl/rmcdhf90_csfg_mpi/newcompi.f90
index 214c5ea..2409b4f 100644
--- a/srcg/appl/rmcdhf90_csfg_mpi/newcompi.f90
+++ b/srcg/appl/rmcdhf90_csfg_mpi/newcompi.f90
@@ -113,6 +113,8 @@
       DO JALL = 1, NCMIN
          NB = IDXBLK(JALL)                       ! Block number of this state
          EE = EAVBLK(NB) + EVAL(JALL)
+         WRITE(*,*) 'JALL=', JALL, 'NB=', NB, 'EAVBLK(NB)=', EAVBLK(NB), &
+                 'EVAL(JALL)=', EVAL(JALL), 'EE=', EE
          IF (MYID .EQ. 0) then
             WRITE (*, 301) ICCMIN(JALL), EE, WT(JALL)
             IF (LDBPG(5)) THEN

and it appears that the problem lies in the EVAL, and ultimately in the WORK arrays:

 J=           1 NEVBLK(JBLOCK)=           2 NEND=           4 NCMINPAT=           0 J + NCMINPAT=           1 ICCMIN(J + NCMINPAT)=           1 NEND + ICCMIN(J + NCMINPAT)=           5 WORK(NEND + ICCMIN(J + NCMINPAT))=   1.0000000000000000E+030
 J=           2 NEVBLK(JBLOCK)=           2 NEND=           4 NCMINPAT=           0 J + NCMINPAT=           2 ICCMIN(J + NCMINPAT)=           2 NEND + ICCMIN(J + NCMINPAT)=           6 WORK(NEND + ICCMIN(J + NCMINPAT))=   1.0000000000000000E+030

...

 JALL=           1 NB=           1 EAVBLK(NB)=  -727.34013468349042      EVAL(JALL)=   1.0000000000000000E+030 EE=   1.0000000000000000E+030
Level  1    Energy =  1.000000000000D+30    Weight =  2.00000D-01
 JALL=           2 NB=           1 EAVBLK(NB)=  -727.34013468349042      EVAL(JALL)=   1.0000000000000000E+030 EE=   1.0000000000000000E+030
Level  2    Energy =  1.000000000000D+30    Weight =  2.00000D-01
 JALL=           3 NB=           2 EAVBLK(NB)=  -564.04488713524256      EVAL(JALL)=   0.0000000000000000      EE=  -564.04488713524256     
Level  1    Energy = -5.640448871352D+02    Weight =  6.00000D-01

He-like_Zn_2s.zip

[maj-mas]

I've also encountered this behaviour, but in a more complicated calculation. Did you find a fix @mchatzikos ?

[maj-mas]

When I try to run this example, I get different behaviour, namely a crash due to an illegal parameter being passed to one of the diag routines of LAPACK:

** On entry to DSPEVX parameter number  9 had an illegal value

This is related to what I talked about over in #2.

Inspecting rlabel.out, we can see that only the 1s subshell is in the labelling space:

           1   = nonsym
  1s

However, in the second J-Block, there is only one CSF which is not a labelling CSF because the 2s subshell is occupied:

Core subshells:

Peel subshells:
  1s   2s
CSF(s):
  1s ( 2)

         0+
  1s ( 1)  2s ( 1)
      1/2      1/2
                  0+
 *
  1s ( 1)  2s ( 1)   <------------- !!
      1/2      1/2
                  1+

graspg cannot deal with this as I explained in the other issue. If you change your call to rcsfggenerate_csfg as such:

rcsfggenerate_csfg <<EOF &> rcsfggenerate_csfg.mr.STDOUTERR
0
1s(2,*)
1s(1,*)2s(1,*)

1s,2s
0,2
0
n
EOF

the calculation finishes just fine with no changes to rcsfg.out but the 2s orbital now added to rlabel.out. note that the first list you generate with rcsfggenerate_csfg defines the labelling space while any further lists add orbitals to the correlation space. one needs to take care to avoid adding states to the multireference whose orbitals are not in the labelling space. in this regard, the usage of rcsfgenerate and rcsfggenerate_csfg is subtly different.

I hope this helps @mchatzikos ! :)