Abstract : We propose a novel numerical method for fast and accurate evaluation of the exchange part of the Fock operator in the Hartree-Fock equation which is a (nonlocal) integral operator in $\mathbb{R}^3 \times \mathbb{R}^3$. Usually, this challenging computational problem is solved by analytical evaluation of two-electron integrals using the ``analytically separable'' Galerkin basis functions, like Gaussians. Instead, we employ the agglomerated ``grey-box'' numerical computation of the corresponding six-dimensional integrals in the tensor-structured format which does not require analytical separability of the basis set. The point of our method is a low-rank tensor representation of arising functions and operators on an $n\times n\times n$ Cartesian grid and the implementation of the corresponding multi-linear algebraic operations in the tensor product format. Linear scaling of the tensor operations, including the 3D convolution product, with respect to the one-dimension grid size $n$ enables computations on huge 3D Cartesian grids thus providing the required high accuracy. The presented algorithm for evaluation of the exchange operator and a recent tensor method for the computation of the Coulomb matrix are the main building blocks in the numerical solution of the Hartree-Fock equation by the tensor-structured methods. These methods provide a new tool for algebraic optimization of the Galerkin basis in the case of large molecules.