anchor.diff:  
       end function get_one_e_int
   
  -    pure function get_two_e_int(phi1, phi2, phi3, phi4) result(two_e_int)
  +    elemental function get_two_e_int(phi1, phi2, phi3, phi4) result(two_e_int)
   
           real(dp) :: two_e_int
           integer, intent(in) :: phi1, phi2, phi3, phi4
  
positive.code:         if (all(delta_k == 0)) then
              conserved = .true.
          else
              forall (i=1:ndim) delta_kc(i) = sum(delta_k*lattice(i,:))
              conserved = all(mod(delta_kc, 1.0_dp) < depsilon)
          end if
  
      end function momentum_conserved
  
negative.code: 
      end function get_two_e_int
  
      elemental function momentum_conserved(i, j, k, l) result(conserved)
  
          logical :: conserved
          integer, intent(in) :: i, j, k, l
  