anchor.diff:            NKL = NK*NL
   *. Obtain block B(KSM,KGAS,LSM,LGAS)
             CALL EXTR_SYMGAS_BLK_FROM_ORBMAT
  -     &    (B,WORK(KLBBLK),KSM,KGAS,LSM,LGAS)
  +     &    (B,dbl_mb(KLBBLK),KSM,KGAS,LSM,LGAS)
   *. Fetch RHO2(ISM,JSM,KSM,LSM)
   C         GETD2(RHO2B,ISM,IGAS,JSM,JGAS,KSM,KGAS,LSM,LGAS,ISPC)
             CALL GETD2
  -     &    (WORK(KLRHO2B),ISM,IGAS,JSM,JGAS,KSM,KGAS,LSM,LGAS,1)
  +     &    (dbl_mb(KLRHO2B),ISM,IGAS,JSM,JGAS,KSM,KGAS,LSM,LGAS,1)
   * sum(kl) RHO2(ij,kl) B(kl)
             FACTORC = 0.0D0
             FACTORAB = 1.0D0
             ZERO = 0.0D0
  -          CALL SETVEC(WORK(KLVEC),ZERO,NIJ)
  -          CALL MATML7(WORK(KLVEC),WORK(KLRHO2B),WORK(KLBBLK),NIJ,1,
  -     &                NIJ,NKL,NKL,1,FACTORC,FACTORAB,0)
  +          CALL SETVEC(dbl_mb(KLVEC),ZERO,NIJ)
  +          CALL MATML7(dbl_mb(KLVEC),dbl_mb(KLRHO2B),dbl_mb(KLBBLK),NIJ,
  +     &                1,NIJ,NKL,NKL,1,FACTORC,FACTORAB,0)
             IF(NTEST.GE.1000) THEN
               WRITE(6,*) ' KLVEC: '
  -            CALL WRTMAT(WORK(KLVEC),1,NIJ,1,NIJ)
  +            CALL WRTMAT(dbl_mb(KLVEC),1,NIJ,1,NIJ)
             END IF
   *. sum(ij) (A(ij) (sum(kl) RHO2(ij,kl)B(kl) )
  -          AB2 = AB2 + INPROD(WORK(KLABLK),WORK(KLVEC),NIJ)
  +          AB2 = AB2 + INPROD(dbl_mb(KLABLK),dbl_mb(KLVEC),NIJ)
             IF(NTEST.GE.100) THEN
               WRITE(6,'(A,4I4,E20.12)') 
        &      ' Two-electron contribution after ISM, JSM, KSM, LSM = ',
  
positive.code: *
        RETURN
        END
        SUBROUTINE LSQ_FOR_ORBS(XLSQ)
  *
  * Obtain l^2 integrals for orbitals in CMOAO_ACT
  *
  *. Jeppe Olsen, HNIE, July 2, 2012
  *
        INCLUDE 'implicit.inc'
        INCLUDE 'mxpdim.inc'
        INCLUDE 'orbinp.inc'
        INCLUDE 'wrkspc-static.inc'
        INCLUDE 'glbbas.inc'
        INCLUDE 'crun.inc'
        INCLUDE 'lucinp.inc'
        INCLUDE 'multd2h.inc'
  *
        CHARACTER*8 LABEL
  *. Output
        DIMENSION XLSQ(*)
  *
        NTEST = 100
        IF(NTEST.GE.100) THEN
          WRITE(6,*) 
          WRITE(6,*) ' Output from LSQ_FOR_ORBS '
          WRITE(6,*) ' ========================='
          WRITE(6,*) 
        END IF
  *
        IDUM = 0 
        CALL MEMMAN(IDUM,IDUM,'MARK  ',IDUM,'LSQ_OB')
        LDIM_H1 = NTOOB**2
        CALL MEMMAN(KLH1,LDIM_H1,'ADDL  ',2,'H1_SCR')
        CALL MEMMAN(KLH2,LDIM_H1,'ADDL  ',2,'H2_SCR')
        CALL MEMMAN(KLH3,LDIM_H1,'ADDL  ',2,'H3_SCR')
        CALL MEMMAN(KLSCR,3*LDIM_H1,'ADDL  ',2,'SCR   ')
  *
        LEN_TOTSYM = NDIM_1EL_MAT(1,NTOOBS,NTOOBS,NSMOB,0)
        ONEM = -1.0D0
        ONE  =  1.0D0
        LABEL(1:6) = 'ANGMOM'
        CALL GET_PROP_PERMSM(LABEL,IPERMSM)
  *
  *. Lx^2:
  *
  *. Lx integrals in MO basis 
        LX_SYM = MULTD2H(IXYZSYM(2),IXYZSYM(3))
        LABEL(1:8) = 'XANGMOM '
        CALL GET_PROPINT(WORK(KLH1),LX_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &             NTOOBS,NSMOB,0,IPERMSM)
  C     MULT_H1H2(H1,IH1SM,H2,IH2SM,H12,IH12SM)
        CALL MULT_H1H2(WORK(KLH1),LX_SYM,WORK(KLH1),LX_SYM,WORK(KLH2),
       &     IH12SM)
  *. The imaginary i in Lx gives rise to a minus,so
        CALL SCALVE(WORK(KLH2),ONEM,LEN_TOTSYM)
        CALL COPVEC(WORK(KLH2),WORK(KLH3),LEN_TOTSYM)
  *
  *. Ly^2:
  *
  *. Ly integrals in MO basis 
        LY_SYM = MULTD2H(IXYZSYM(1),IXYZSYM(3))
        LABEL(1:8) = 'YANGMOM '
        CALL GET_PROPINT(WORK(KLH1),LY_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &             NTOOBS,NSMOB,0,IPERMSM)
  C     MULT_H1H2(H1,IH1SM,H2,IH2SM,H12,IH12SM)
        CALL MULT_H1H2(WORK(KLH1),LY_SYM,WORK(KLH1),LY_SYM,WORK(KLH2),
       &     IH12SM)
  *. The imaginary i in Ly gives rise to a minus,so
        CALL VECSUM(WORK(KLH3),WORK(KLH3),WORK(KLH2),ONE,ONEM,LEN_TOTSYM)
  *
  *. Lz^2:
  *
  *. Lz integrals in MO basis 
        LZ_SYM = MULTD2H(IXYZSYM(1),IXYZSYM(2))
        LABEL(1:8) = 'ZANGMOM '
        CALL GET_PROPINT(WORK(KLH1),LZ_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &             NTOOBS,NSMOB,0,IPERMSM)
  C     MULT_H1H2(H1,IH1SM,H2,IH2SM,H12,IH12SM)
        CALL MULT_H1H2(WORK(KLH1),LZ_SYM,WORK(KLH1),LZ_SYM,WORK(KLH2),
       &     IH12SM)
  *. The imaginary i in Lz gives rise to a minus,so
        CALL VECSUM(WORK(KLH3),WORK(KLH3),WORK(KLH2),ONE,ONEM,LEN_TOTSYM)
  *
        CALL COPVEC(WORK(KLH3),XLSQ,LEN_TOTSYM)
  *
        IF(NTEST.GE.100) THEN
          WRITE(6,*) ' The l^2 operator in the orbital basis '
          CALL APRBLM2(XLSQ,NTOOBS,NTOOBS,NSMOB,0)
        END IF
  *
        CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'LSQ_OB')
  *
        RETURN
        END
        SUBROUTINE GET_LML_FROM_CMOAO(ISUPSYM_ORBS,L_ORB,ML_ORB)
  *
  * Obtain L and ML values of orbitals stored in KMOAO_ACT
  * (L only if atomic symmetry)
  *
  *. Jeppe Olsen, HNIE, July 2, 2012
  *
        INCLUDE 'implicit.inc'
        INCLUDE 'mxpdim.inc'
        INCLUDE 'orbinp.inc'
        INCLUDE 'lucinp.inc'
        INCLUDE 'crun.inc'
        INCLUDE 'multd2h.inc'
        INCLUDE 'wrkspc-static.inc'
  *
        CHARACTER*8 LABEL
  *. Output
        INTEGER L_ORB(*), ML_ORB(*) 
  *
        IDUM = 0
        CALL MEMMAN(IDUM,IDUM,'MARK  ',2,'GTLMLA')
  *
        NTEST = 1000
        IF(NTEST.GE.1000) THEN
          WRITE(6,*) ' Info from GET_LML_FROM_CMOAO'
          WRITE(6,*) ' -----------------------------'
        END IF
  *
        LEN_H1 = NTOOB*2
        CALL MEMMAN(KLLZ ,LEN_H1,'ADDL  ',2,'LZ_ORB')
        CALL MEMMAN(KLLZ2,LEN_H1,'ADDL  ',2,'LZ2ORB')
        CALL MEMMAN(KLSCR,3*LEN_H1,'ADDL  ',2,'LZ2ORB')
        IF(CSUPSYM.EQ.'ATOMIC')  THEN
          CALL MEMMAN(KLL2,LEN_H1,'ADDL  ',2,'LZ2ORB')
          CALL MEMMAN(KLL2_AR,NTOOB,'ADDL  ',2,'LZ2OBA')
        END IF
  * 
        LABEL(1:6) = 'ANGMOM'
        CALL GET_PROP_PERMSM(LABEL,IPERMSM)
        LEN_TOTSYM = NDIM_1EL_MAT(1,NTOOBS,NTOOBS,NSMOB,0)
        ONE = 1.0D0
        ONEM = -1.0D0
  *
  * 1. obtain the representation of L^2 for the orbitals - for atomic supersymmetry
  * 
        IF(CSUPSYM.EQ.'ATOMIC') THEN
  *
  *.  Lx^2
  *
          LX_SYM = MULTD2H(IXYZSYM(2),IXYZSYM(3))
          LABEL(1:8) = 'XANGMOM '
          CALL GET_PROPINT(WORK(KLLZ),LX_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &               NTOOBS,NSMOB,0,IPERMSM)
  *. And Lx^2 
  C       MULT_H1H2(H1,IH1SM,H2,IH2SM,H12,IH12SM)
          CALL MULT_H1H2(WORK(KLLZ),LX_SYM,WORK(KLLZ),LX_SYM,WORK(KLLZ2),
       &       IH12SM)
  *. The imaginary i in Lx gives rise to a minus,so
          CALL SCALVE(WORK(KLLZ2),ONEM,LEN_TOTSYM)
          CALL COPVEC(WORK(KLLZ2),WORK(KLL2),LEN_TOTSYM)
  *
  *. + Ly^2
  *
          LY_SYM = MULTD2H(IXYZSYM(1),IXYZSYM(3))
          LABEL(1:8) = 'YANGMOM '
          CALL GET_PROPINT(WORK(KLLZ),LY_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &               NTOOBS,NSMOB,0,IPERMSM)
  C       MULT_H1H2(H1,IH1SM,H2,IH2SM,H12,IH12SM)
          CALL MULT_H1H2(WORK(KLLZ),LY_SYM,WORK(KLLZ),LY_SYM,WORK(KLLZ2),
       &       IH12SM)
  *. The imaginary i in Lx gives rise to a minus,so
          CALL VECSUM(WORK(KLL2),WORK(KLL2),WORK(KLLZ2),ONE,ONEM,
       &       LEN_TOTSYM)
  *
  *. + Lz^2
  *
          LZ_SYM = MULTD2H(IXYZSYM(1),IXYZSYM(2))
          LABEL(1:8) = 'ZANGMOM '
          CALL GET_PROPINT(WORK(KLLZ),LZ_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &               NTOOBS,NSMOB,0,IPERMSM)
          CALL MULT_H1H2(WORK(KLLZ),LZ_SYM,WORK(KLLZ),LZ_SYM,
       &       WORK(KLLZ2),IH12SM)
  *. The imaginary i in Lz gives rise to a minus,so
          CALL VECSUM(WORK(KLL2),WORK(KLL2),WORK(KLLZ2),ONE,ONEM,
       &       LEN_TOTSYM)
          IF(NTEST.GE.1000) THEN
            WRITE(6,*)' The L^2 matrix '
            CALL APRBLM2(WORK(KLL2),NTOOBS,NTOOBS,NSMOB,0)
          END IF
  *. Determine the L-values of the orbitals
  C DIAG_BLOC_MAT(A,ADIAG,NBLOCK,LBLOCK,ISYM)
          CALL GET_DIAG_BLOC_MAT(WORK(KLL2),WORK(KLL2_AR),NSMOB,NTOOBS,0)
          THRES = 1.0D-8
          DO IORB = 1, NTOOB
            XL = WORK(KLL2_AR-1+IORB)
            DO L = 0, MXPL
              LLP1 = L*(L+1)
              IF(ABS(XL-FLOAT(LLP1)).LE.THRES) L_ORB(IORB) = L
            END DO
          END DO
  *
          IF(NTEST.GE.100) THEN
            WRITE(6,*) ' The L-value of the orbitals '
            WRITE(6,*) ' ============================'
            CALL IWRTMA3(L_ORB,1,NTOOB,1,NTOOB)
          END IF
  *
        END IF! Atomic supersymmetry
  *
  *. The Lz and Lz2 arrays
  *
        LZ_SYM = MULTD2H(IXYZSYM(1),IXYZSYM(2))
        LABEL(1:8) = 'ZANGMOM '
        CALL GET_PROPINT(WORK(KLLZ),LZ_SYM,LABEL,WORK(KLSCR),NTOOBS,
       &             NTOOBS,NSMOB,0,IPERMSM)
        CALL MULT_H1H2(WORK(KLLZ),LZ_SYM,WORK(KLLZ),LZ_SYM,
       &     WORK(KLLZ2),IH12SM)
  *. The imaginary i in Lz gives rise to a minus,so
        CALL SCALVE(WORK(KLLZ2),ONEM,LEN_TOTSYM)
        IF(NTEST.GE.1000) THEN
          WRITE(6,*) ' The Lz and Lz^2  arrays '
          CALL APRBLM2(WORK(KLLZ),NTOOBS,NTOOBS,NSMOB,0)
          WRITE(6,*)
          CALL APRBLM2(WORK(KLLZ2),NTOOBS,NTOOBS,NSMOB,0)
        END IF
  *
  * We are now after the ML-values of the orbitals. The orbitals are
  * real spherical Harmonics and due therefore not have a well-defined ML.
  *.In the real sherical Harmonics the operator -i Lz (which we are using) have 
  * the representation
  *
  * <S(L,-M)!S(L,-M)>             <S(L,-M)!S(L,M)>           0         M
  *                                                     = 
  * <S(L,M)!S(L,-M)>              <S(L,M)!S(L,M)>           -M         0
  *
        THRES = 1.0D-10
        IJ_OFF = 1
        IORB = 0
        DO ISM = 1, NSMOB
          JSM = MULTD2H(ISM,LZ_SYM)
  *
          NI = NTOOBS(ISM)
          NJ = NTOOBS(JSM)
          DO I = 1, NI
           IORB = IORB + 1
  *.  Find largest off-diagonal element in LZ(I,*)
           JMAX = 0
           XMAX = 0.0D0
           DO J = 1, NJ
             XLZIJ = WORK(KLLZ-1+IJ_OFF - 1 + (J-1)*NI+I)
             IF(ABS(XLZIJ).GT.ABS(XMAX)) THEN
               XMAX = XLZIJ
               JMAX = J
             END IF
           END DO
  *
           ML_FOUND = 0
           DO ML = -MXPL, MXPL
             XML = FLOAT(ML)
             IF(ABS(XMAX-XML).LE.THRES) THEN
               ML_ORB(IORB) = ML
               ML_FOUND = 1
             END IF
           END DO
           IF(ML_FOUND.EQ.0) THEN
             WRITE(6,*) 
       &     ' Warning ML value of orbital could not be defined '
           END IF
          END DO ! loop over IORB
          IJ_OFF = IJ_OFF + NI*NJ
        END DO
  *
        IF(NTEST.GE.100) THEN
          WRITE(6,*) ' Pseudo ML value of orbitals '
          WRITE(6,*) ' ============================'
          CALL IWRTMA3(ML_ORB,1,NTOOB,1,NTOOB)
        END IF
  *
        CALL MEMMAN(IDUM,IDUM,'FLUSM ',2,'GTLMLA')
  *
        RETURN
        END
  
negative.code: *. Obtain block A(ISM,IGAS,JSM,JGAS)
  C       EXTR_SYMGAS_BLK_FROM_ORBMAT(A,ABLK,ISM,IGAS,JSM,JGAS)
          CALL EXTR_SYMGAS_BLK_FROM_ORBMAT
       &  (A,dbl_mb(KLABLK),ISM,IGAS,JSM,JGAS)
          DO KSM = 1, NSMOB
           LSM = MULTD2H(KSM,IBSM)
           DO KGAS = 1, NGAS
  