anchor.diff:          END IF
          END IF ! determinant diagonal should be calculated
   *
  -       IF(NOCSF.EQ.0.AND.IH0_CSF.EQ.1) THEN
  -*. Average of determinant diagonal in each conf
  -C?       WRITE(6,*) ' IH0_CSF = ', IH0_CSF
  -         CALL CSDIAG(WORK(KPVEC2),WORK(KPDIA_SD),
  -     &        NCONF_PER_OPEN(1,ISM),MAXOP,ISM,
  -     &        int_mb(KSDREO_I(ISM)),NPCMCNF,NPCSCNF,IPRCSF,
  -     &        ICNFBAT,NCOCCLS_ACT,dbl_mb(KCIOCCLS_ACT),
  -     &        int_mb(KCLBT),LUDIA,LUSC52)
  -C        CSDIAG(CSFDIA,DETDIA,NCNFTP,MAXOP,ISM,
  -C    &         ICTSDT,NPCMCNF,NPCSCNF,IPRCSF,
  -C    &         ICNFBAT,NOCCLS_ACT,IOCCLS_ACT,
  -C    &         LBLOCK,LUDIA_DET,LUDIA_CSF)
  -*
  -         IF(ICNFBAT.EQ.1) THEN
  -           CALL COPVEC(WORK(KPVEC2),WORK(KPDIA_SD),NVAR)
  -         ELSE
  -           CALL COPVCD(LUSC52,LUDIA,WORK(KPVEC2),1,-1)
  -         END IF
  -       END IF
  -*
  -       IF(NOCSF.EQ.0.AND.IH0_CSF.GE.2) THEN
  -*. Use exact diagonal or blocks
  -C?       WRITE(6,*) ' Ecore before GET_CSF.... ', ECORE
  -C?       WRITE(6,*) ' KPDIA_SD = ', KPDIA_SD
  -         CALL GET_CSF_H_PRECOND(NCONF_PER_OPEN(1,ISM),
  -     &        int_mb(KICONF_OCC(ISM)),WORK(KPDIA_SD),ECORE,
  -     &        LUDIA,NCOCCLS_ACT,dbl_mb(KCIOCCLS_ACT),ISM)
  -C        GET_CSF_H_PRECOND(NCONF_FOR_OPEN,ICONF_OCC,H0,ECORE,
  -C    &   LUDIA,NOCCLS_ACT,IOCCLS_ACT),ISYM)
  -       END IF
  -*
  -       IF(ICISTR.EQ.1) THEN
  -         CALL REWINO(LUDIA)
  +       CALL REWINO(LUDIA)
   C?     WRITE(6,*) ' WORK(KPDIA_SD) before TODSC =', WORK(KPDIA_SD)
  -         CALL TODSC(WORK(KPDIA_SD),NVAR,NVAR,LUDIA)
  -       END IF
  +       CALL TODSC(KPDIA_SD,NVAR,NVAR,LUDIA)
   *
          IF(IIUSEH0P.EQ.1) THEN
   *. Is pt not reprogrammed for CSF's and ICISTR = 1
   *. Diagonal with F
            CALL SWAPVE(WORK(KFI),WORK(KINT1O),NINT1)
  -         CALL GASDIAT(WORK(KVEC1),LUSC52,SHIFT,ICISTR,1,
  +         CALL GASDIAT(dbl_mb(KVEC1),LUSC52,SHIFT,ICISTR,1,
        &              dbl_mb(KCBLTP),NBLOCK,int_mb(KCIBT),IUSE_EXP)
            CALL SWAPVE(WORK(KFI),WORK(KINT1O),NINT1)
   *. diag of (1-Lambda) F + Lambda H
            FAC1 = 1.0D0 - XLAMBDA
            FAC2 = XLAMBDA
  -         CALL VECSMD(WORK(KVEC1),WORK(KVEC2),FAC1,FAC2,
  +         CALL VECSMD(dbl_mb(KVEC1),dbl_mb(KVEC2),FAC1,FAC2,
        &   LUSC52,LUDIA,LUSC53,1,LBLK)
  -         CALL COPVCD(LUSC53,LUDIA,WORK(KVEC1),1,LBLK)
  +         CALL COPVCD(LUSC53,LUDIA,dbl_mb(KVEC1),1,LBLK)
          END IF !IIUSEH0P.EQ.1
   *
          IF(IPRCIX.GE.2) WRITE(6,*) ' Diagonal constructed  '
   *
          IF(IPRNT.GE.100) THEN
           WRITE(6,*) ' Constructed diagonal '
  -        IF(ICISTR.GE.2) THEN
  -          CALL WRTVCD(WORK(KVEC1),LUDIA,1,LBLK)
  -        ELSE
  -          CALL WRTMAT(WORK(KPDIA_SD),1,NVAR,1,NVAR)
  -        END IF
  +        CALL WRTMAT(KPDIA_SD,1,NVAR,1,NVAR)
          END IF
   *
         ELSE
  
positive.code:      &  ' **************************************************'
          WRITE(6,*)
          END IF
          IF(ICISTR.EQ.1) THEN
            CALL FRMDSC(WORK(KPVEC1),NVAR,-1,LUC,IMZERO,IAMPACK)
            CALL COPVEC(WORK(KPVEC1),WORK(KPVEC2),NVAR)
  *
  C?        WRITE(6,*) ' Testing: CI vector read in '
  C?        CALL WRTMAT(WORK(KPVEC1),1,NVAR,1,NVAR)
  *
            IF(IDENSI.EQ.1) THEN
              KRHO2 = 1
            END IF
            IF(IDENSI.GE.1)
       &    CALL DENSI2(IDENSI,WORK(KRHO1),WORK(KRHO2),
       &         WORK(KPVEC1),WORK(KPVEC2),0,0,EXPS2,ISPNDEN,WORK(KSRHO1),
       &         WORK(KRHO2AA),WORK(KRHO2AB),WORK(KRHO2BB),1)
            IF(IDENSI.EQ.2) THEN
  *. Test calculation of energy
              CALL EN_FROM_DENS(ENERGY,2,0)
              WRITE(6,*) ' Energy from density ', ENERGY
            END IF
  *
  C?        IF(I_DO_LZ2.EQ.1) THEN
  C?            CALL EXP_LZ2(KPVEC1,KPVEC2,RLZEFF,RL2EFF,0)
  C?            WORK(KLZEXP-1+JROOT) = RLZEFF
  C?            WORK(KL2EXP-1+JROOT) = RL2EFF
  C?        END IF
  *
            IF(IPRNT.GT.2) THEN
             IF(NOCSF.EQ.1) THEN
              CALL GASANA(WORK(KPVEC1),NBLOCK,int_mb(KCIBT),
       &                  int_mb(KCBLTP),LUC,ICISTR)
             ELSE
              THRES = 0.1D0
              MAXTRM = 100
              IOUT = 6
  C                ANACSF(CIVEC,ICONF_OCC,NCONF_FOR_OPEN,IPROCS,THRES,MAXTRM,IOUT)
              IF(ICNFBAT.EQ.1) THEN
                CALL ANACSF(WORK(KPVEC1),int_mb(KICONF_OCC(ISM)),
       &                    NCONF_PER_OPEN(1,ISM),int_mb(KCFTP),
       &                    THRES,MAXTRM,IOUT)
              ELSE
              THRES = 0.1D0
              MAXTRM = 100
                IOUT = 6
                CALL ANACSF2(LUC,NCOCCLS_ACT,dbl_mb(KCIOCCLS_ACT),ISM,
       &             WORK(KPVEC1),int_mb(KICONF_OCC(ISM)),
       &             NCONF_PER_OPEN(1,ISM),int_mb(KCFTP),THRES,
       &             MAXTRM,IOUT)
  C      ANACSF2(LUC,NOCCLS_SPC,IOCCLS_SPC,
  C     &           CIVEC,ICONF_OCC,NCONF_FOR_OPEN,IPROCS,THRES,
  C     &           MAXTRM,IOUT)
                END IF ! CNF batch switch
             END IF! CSF switch
            END IF ! IPRNT > 2
          ELSE ! ICISTR switch
            IF(IDENSI.GE.1.OR.NROOT.GT.1) THEN
               CALL REWINO(LUSC1)
               CALL COPVCD(LUC,LUSC1,WORK(KVEC1),0,LBLK)
               CALL COPVCD(LUSC1,LUHC,WORK(KVEC1),1,LBLK)
            END IF
            IF(IDENSI.GE.1) 
       &    CALL DENSI2(IDENSI,WORK(KRHO1),WORK(KRHO2),
       &          WORK(KVEC1),WORK(KVEC2),LUSC1,LUHC,EXPS2,
       &          ISPNDEN,WORK(KSRHO1),WORK(KRHO2AA),WORK(KRHO2AB),
       &          WORK(KRHO2BB),1 ) 
            IF(IPRNT.GT.1.AND.IDENSI.EQ.2) 
       &    WRITE(6,*) ' Expectation value of S2', EXPS2
            IF(IDENSI.EQ.2) THEN
  *. Test calculation of energy
              CALL EN_FROM_DENS(ENERGY,2,0)
              WRITE(6,*) ' Energy from density ', ENERGY
            END IF
  *
  C?        IF(I_DO_LZ2.EQ.1) THEN
  C?            IF(NROOT.EQ.1) THEN
  C?              LUCEFF = LUC
  C?            ELSE
  C?              LUCEFF = LUSC1
  C?            END IF
  C?            WRITE(6,*) ' EXP_LZ2 will be called '
  C?            CALL EXP_LZ2(KPVEC1,KPVEC2,RLZEFF,RL2EFF,
  C?   &              LUCEFF)
  C?            WORK(KLZEXP-1+JROOT) = RLZEFF
  C?            WORK(KL2EXP-1+JROOT) = RL2EFF
  C?        END IF
  *
            IF(ISPC.EQ.1.AND.JROOT.EQ.IH0ROOT) THEN 
  *. Calculate inactive Fock matrix for reference root in first space 
               CALL COPVEC(WORK(KINT1),WORK(KFIZ),NINT1)  
               CALL FIFAM(WORK(KFIZ))
            END IF
  *
            IF(IPRNT.GT.2) THEN
             IF(NROOT.GT.1) THEN
              CALL REWINO(LUSC1)
              IF(NOCSF.EQ.1) THEN
                CALL GASANA(WORK(KVEC1),NBLOCK,int_mb(KCIBT),
       &                    int_mb(KCBLTP),LUSC1,ICISTR)
              ELSE
                IOUT = 6
                THRES = 0.1D0
                MAXTRM = 100
                CALL ANACSF2(LUC,NCOCCLS_ACT,dbl_mb(KCIOCCLS_ACT),ISM,
       &             WORK(KPVEC1),int_mb(KICONF_OCC(ISM)),
       &             NCONF_PER_OPEN(1,ISM),int_mb(KCFTP),THRES,
       &             MAXTRM,IOUT)
              END IF
             ELSE 
              CALL REWINO(LUC)
              IF(NOCSF.EQ.1) THEN
                CALL GASANA(WORK(KVEC1),NBLOCK,int_mb(KCIBT),
       &                    int_mb(KCBLTP),LUC,ICISTR)
              ELSE
                IOUT = 6
                THRES = 0.1D0
                MAXTRM = 100
                CALL ANACSF2(LUC,NCOCCLS_ACT,dbl_mb(KCIOCCLS_ACT),ISM,
       &             WORK(KPVEC1),int_mb(KICONF_OCC(ISM)),
       &             NCONF_PER_OPEN(1,ISM),int_mb(KCFTP),THRES,
       &             MAXTRM,IOUT)
              END IF !NOCSF switch
             END IF !Nroot switch
            END IF !IPRNT .gt. 2
  C?        WRITE(6,*) ' Back from GASANA ' 
  *
            IF(IPRNCIV.EQ.1) THEN 
               WRITE(6,*)
               WRITE(6,*) ' ======================'
               WRITE(6,*) ' Complete CI expansion '
               WRITE(6,*) ' ======================'
               WRITE(6,*)
               IF(NROOT.EQ.1) THEN
                 CALL WRTVCD(WORK(KVEC1),LUC  ,1,LBLK) 
               ELSE
                 CALL WRTVCD(WORK(KVEC1),LUSC1,1,LBLK) 
               END IF
            END IF
  *         ^ End if print of CI vector
          END IF
  *       ^ End of ICISTR switch
        END DO
  *     ^ End of loop over roots
  *
  
negative.code:             IPRDEN = 0
              XDUM = 0
              CALL DENSI2(1     ,WORK(KRHO1),WORK(KRHO2),
       &           dbl_mb(KVEC1),dbl_mb(KVEC2),LUSC1,LUHC,EXPS2,
       &           0,XDUM,XDUM,XDUM,XDUM,1)
              IPRDEN = IPRDEN_SAVE
              IF(IROOT1.EQ.IROOT2) THEN
  