anchor.diff:  *. Position vectors 
            CALL REWINO(LUC)
            DO JROOT = 1, IH0ROOT
  -           IF(ICISTR.EQ.1) THEN
  -             CALL FRMDSC(WORK(KVEC1),NVAR,-1,LUC,IMZERO,IAMPACK)
  -           ELSE
  -             CALL REWINO(LUSC1)
  -             CALL COPVCD(LUC,LUSC1,WORK(KVEC1),0,LBLK)
  -           END IF
  +           CALL FRMDSC(dbl_mb(KVEC1),NVAR,-1,LUC,IMZERO,IAMPACK)
            END DO
  -         IF(ICISTR.EQ.1) THEN
  -           CALL COPVEC(WORK(KVEC1),WORK(KVEC2),NVAR)
  -           CALL DENSI2(1,WORK(KRHO1),WORK(KRHO2),
  -     &          WORK(KVEC1),WORK(KVEC2),0,0,EXPS2,
  -     &          ISPNDEN,WORK(KSRHO1),WORK(KRHO2AA),WORK(KRHO2AB),
  -     &          WORK(KRHO2BB),1 )           
  -         ELSE
  -           CALL REWINO(LUSC1)
  -           CALL COPVCD(LUSC1,LUSC2,WORK(KVEC1),1,LBLK)
  -             CALL DENSI2(1,WORK(KRHO1),WORK(KRHO2),
  -     &            WORK(KVEC1),WORK(KVEC2),LUSC1,LUSC2,EXPS2,ISPNDEN,
  -     &            WORK(KSRHO1),WORK(KRHO2AA),WORK(KRHO2AB),
  -     &            WORK(KRHO2BB),1)
  -         END IF
  +         CALL COPVEC(dbl_mb(KVEC1),dbl_mb(KVEC2),NVAR)
  +         CALL DENSI2(1,WORK(KRHO1),WORK(KRHO2),
  +     &        dbl_mb(KVEC1),dbl_mb(KVEC2),0,0,EXPS2,
  +     &        ISPNDEN,WORK(KSRHO1),WORK(KRHO2AA),WORK(KRHO2AB),
  +     &        WORK(KRHO2BB),1 )           
          END IF
   *
          IF(IEXTKOP.EQ.1) THEN
  
positive.code:         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)
  C?     WRITE(6,*) ' WORK(KPDIA_SD) before TODSC =', WORK(KPDIA_SD)
           CALL TODSC(WORK(KPDIA_SD),NVAR,NVAR,LUDIA)
         END IF
  *
         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,
       &              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,
       &   LUSC52,LUDIA,LUSC53,1,LBLK)
           CALL COPVCD(LUSC53,LUDIA,WORK(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
         END IF
  *
        ELSE
  
negative.code:             WRITE(6,*) 
  *. Vector IROOT1 on  LUSC1, IROOT2 on LUHC
  C                SKPVCD(LU,NVEC,SEGMNT,IREW,LBLK)
              CALL SKPVCD(LUC,IROOT1-1,dbl_mb(KVEC1),1,LBLK)
              CALL REWINO(LUSC1)
              CALL COPVCD(LUC,LUSC1,dbl_mb(KVEC1),0,LBLK)
  *
              CALL SKPVCD(LUC,IROOT2-1,dbl_mb(KVEC1),1,LBLK)
              CALL REWINO(LUHC)
              CALL COPVCD(LUC,LUHC,dbl_mb(KVEC1),0,LBLK)
  *  <IROOT1!E!IROOT2>
              IF(IDENSI.EQ.1) THEN
                KRHO2 = 1
  