anchor.diff:  
         deallocate(V)
   
  -      allocate( V(a_start:a_end, b_start:b_end, c_start:c_end) )
  -      V = 0.d0
  -      do b=b_start,b_end
  -        call DGEMM('N','N', (a_end-a_start+1), (c_end-c_start+1),    &
  -            (k_end-k_start+1), 1.d0,                                 &
  -            T(a_start,k_start,b), size(T,1),                         &
  -            matrix_B(k_start,k_start), size(matrix_B,1), 1.d0,       &
  -            V(a_start,c_start,b), size(V,1) )
  +      do b=b_start,d
  +        call DGEMM('N','N', (b-a_start+1), (c_end-c_start+1),    &
  +            (k_end-k_start+1), matrix_B(l, d),                   &
  +            T(a_start,k_start,b), size(T,1),                     &
  +            matrix_B(k_start,k_start), size(matrix_B,1), 1.d0,   &
  +            U(a_start,c_start,b), size(U,1) )
         enddo
   
  -      U = U + V*matrix_B(l, d)
  -      deallocate(T,V)
  +      deallocate(T)
   
       enddo
   
       idx = 0_8
  -    do c=c_start,c_end
  -      do b=b_start,b_end
  -        do a=a_start,a_end
  -!    do c=c_start,c_end
  -!      do b=b_start,d
  -!        do a=a_start,min(b,c)
  +    do b=b_start,d
  +      do c=c_start,c_end
  +        do a=a_start,min(b,c)
             if (dabs(U(a,c,b)) < 1.d-15) then
               cycle
             endif
  
positive.code:             call bielec_integrals_index(i,j,k,l,idx)
              call map_get(map_a,idx,tmp)
              T(i, k,j) = tmp
              T(k, i,j) = tmp
            enddo
          enddo
        enddo
  
  
        call DGEMM('T','N', (a_end-a_start+1),                         &
            (k_end-k_start+1)*(j_end-j_start+1),                       &
  
negative.code:           V(a_start,k_start,j_start), size(V, 1) )
  
        deallocate(T)
        allocate( T(a_start:a_end, k_start:k_end, b_start:d) )
  
        call DGEMM('N','N', (a_end-a_start+1)*(k_end-k_start+1),       &
                (d-b_start+1),                                     &
                (j_end-j_start+1), 1.d0,                               &
                V(a_start,k_start,j_start), size(V,1)*size(V,2),       &
                matrix_B(j_start,b_start), size(matrix_B,1),0.d0,      &
  