anchor.diff:        IF(I_USE_SUPSYM.EQ.1) THEN
   *. Bring first density to form with supersymmetry blocks
   C            REFORM_RHO1_TO_GNSM(RHO1_ST,RHO1_GNSM_ST,IWAY,IREO_GNSYM_TO_TS)
  -        CALL REFORM_RHO1_TO_GNSM(WORK(KRHO1),WORK(KLMAT2),1,
  +        CALL REFORM_RHO1_TO_GNSM(dbl_mb(KRHO1),WORK(KLMAT2),1,
        &       WORK(KIREO_GNSYM_TO_TS_ACOB))
           IF(NTEST.GE.1000) THEN
             WRITE(6,*) ' Reformed density matrix '
  
positive.code:             KRHO2 = 1
            END IF
            IF(IDENSI.GE.1)
       &    CALL DENSI2(IDENSI,WORK(KRHO1),WORK(KRHO2),
       &         KPVEC1,KPVEC2,0,0,EXPS2,ISPNDEN,WORK(KSRHO1),
       &         WORK(KRHO2AA),WORK(KRHO2AB),WORK(KRHO2BB),1)
            IF(IDENSI.EQ.2) THEN
  
negative.code:         CALL AVE_SUPSYM_MAT(WORK(KLMAT2),NACOB_GS,0)
  *. And reform back to standard form of density
          ZERO = 0.0D0
          CALL SETVEC(dbl_mb(KRHO1),ZERO,NACOB**2)
          CALL REFORM_RHO1_TO_GNSM(dbl_mb(KRHO1),WORK(KLMAT2),2,
       &       WORK(KIREO_GNSYM_TO_TS_ACOB))
          IF(NTEST.GE.0000) THEN
           WRITE(6,*) ' Supersymmetry averaged density '
           CALL  WRTMAT(dbl_mb(KRHO1),NACOB,NACOB,NACOB,NACOB)
          END IF
        END IF
  *
  