/*
 *            Copyright 2009-2019 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 *Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include <votca/xtp/aobasis.h>
#include <votca/xtp/fourcenter.h>

namespace votca {
namespace xtp {

void FCMatrix::Fill_4c_small_molecule(const AOBasis& dftbasis) {
  Index dftBasisSize = dftbasis.AOBasisSize();
  Index vectorSize = (dftBasisSize * (dftBasisSize + 1)) / 2;

  try {
    _4c_vector = Eigen::VectorXd::Zero((vectorSize * (vectorSize + 1)) / 2);
  } catch (std::bad_alloc& ba) {
    throw std::runtime_error(
        "Basisset too large for 4c calculation. Not enough RAM.");
  }
  Index shellsize = dftbasis.getNumofShells();
#pragma omp parallel for schedule(dynamic)
  for (Index i = 0; i < shellsize; ++i) {

    const AOShell& shell_3 = dftbasis.getShell(i);
    Index start_3 = shell_3.getStartIndex();
    Index NumFunc_3 = shell_3.getNumFunc();

    for (Index j = i; j < shellsize; ++j) {
      const AOShell& shell_4 = dftbasis.getShell(j);
      Index start_4 = shell_4.getStartIndex();
      Index NumFunc_4 = shell_4.getNumFunc();

      for (Index k = i; k < shellsize; ++k) {
        const AOShell& shell_1 = dftbasis.getShell(k);
        Index start_1 = shell_1.getStartIndex();
        Index NumFunc_1 = shell_1.getNumFunc();

        for (Index l = k; l < shellsize; ++l) {
          const AOShell& shell_2 = dftbasis.getShell(l);
          Index start_2 = shell_2.getStartIndex();
          Index NumFunc_2 = shell_2.getNumFunc();

          Eigen::Tensor<double, 4> block(NumFunc_1, NumFunc_2, NumFunc_3,
                                         NumFunc_4);
          block.setZero();

          bool nonzero =
              FillFourCenterRepBlock(block, shell_1, shell_2, shell_3, shell_4);

          if (nonzero) {

            for (Index i_3 = 0; i_3 < NumFunc_3; i_3++) {
              Index ind_3 = start_3 + i_3;
              Index sum_ind_3 = (ind_3 * (ind_3 + 1)) / 2;
              for (Index i_4 = 0; i_4 < NumFunc_4; i_4++) {
                Index ind_4 = start_4 + i_4;
                if (ind_3 > ind_4) {
                  continue;
                }
                Index index_34 = dftBasisSize * ind_3 - sum_ind_3 + ind_4;
                Index index_34_12_a =
                    vectorSize * index_34 - (index_34 * (index_34 + 1)) / 2;
                for (Index i_1 = 0; i_1 < NumFunc_1; i_1++) {
                  Index ind_1 = start_1 + i_1;
                  Index sum_ind_1 = (ind_1 * (ind_1 + 1)) / 2;
                  for (Index i_2 = 0; i_2 < NumFunc_2; i_2++) {
                    Index ind_2 = start_2 + i_2;
                    if (ind_1 > ind_2) {
                      continue;
                    }
                    Index index_12 = dftBasisSize * ind_1 - sum_ind_1 + ind_2;
                    if (index_34 > index_12) {
                      continue;
                    }
                    _4c_vector(index_34_12_a + index_12) =
                        block(i_1, i_2, i_3, i_4);

                  }  // i_2
                }    // i_1
              }      // i_4
            }        // i_3

          }  // end if
        }    // DFT shell_2
      }      // DFT shell_1
    }        // DFT shell_4
  }          // DFT shell_3

  return;
}  // FCMatrix_dft::Fill_4c_small_molecule

}  // namespace xtp
}  // namespace votca