ci_configuration2.h

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
//
// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
//////////////////////////////////////////////////////////////////////////////////////


#ifndef QMCPLUSPLUS_CI_CONFIGURATION2_H
#define QMCPLUSPLUS_CI_CONFIGURATION2_H

#include <algorithm>
#include <iostream>
#include <stdexcept>

namespace qmcplusplus
{
// Defines a ci configuration based of the occupied MOs
struct ci_configuration2
{
  // occupied orbitals
  std::vector<size_t> occup;

  ci_configuration2() {}

  ci_configuration2(std::vector<size_t>& v) : occup(v.begin(), v.end()) {}
  ci_configuration2(const ci_configuration2& c) : occup(c.occup) {}

  ~ci_configuration2() {}

  inline bool operator==(const ci_configuration2& c) const
  {
    if (occup.size() != c.occup.size())
      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");

    for (int i = 0; i < occup.size(); i++)
      if (occup[i] != c.occup[i])
        return false;
    return true;
  }

  /*
  * The routine determines the number excitations that
  * transform the given configuration (argument) into (*this).
  */
  size_t calculateNumOfExcitations(const ci_configuration2& c) const
  {
    if (occup.size() != c.occup.size())
      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");

    size_t n = 0;
    for (size_t i = 0; i < occup.size(); i++)
      if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end())
        n++;
    return n;
  }

  /*
   * The routine determines the excitations that
   * transform the given configuration (argument) into (*this).
   * The sign of the permutation is returned.
   *    - n: number of excitations
   *    - pos: positions of MOs in c that need to be changed (not the label of the MO)
   *    - ocp: label of the MO being removed (needed to calculate matrix elements)
   *    - uno: label of the MO that replaces ocp[i] (should be the same as the position in array, as opposed to ocp)
   *
   */
  double calculateExcitations(const ci_configuration2& c,
                              size_t& n,
                              std::vector<size_t>& pos,
                              std::vector<size_t>& ocp,
                              std::vector<size_t>& uno) const
  {
    if (occup.size() != c.occup.size())
      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");

    n = 0;
    for (size_t i = 0; i < occup.size(); i++)
      if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end())
      {
        pos[n]   = i;
        ocp[n++] = occup[i];
      }

    if (n == 0)
      return 1.0;

    size_t cnt = 0;
    for (size_t i = 0; i < c.occup.size(); i++)
      if (std::find(occup.begin(), occup.end(), c.occup[i]) == occup.end())
        uno[cnt++] = c.occup[i];

    if (cnt != n)
      throw std::runtime_error(" Error #1 in ci_configuration2::calculateExcitations() \n");

    double res = 1.0;
    // this is needed because ci coefficients are given wrt standard ordering,
    // but by defining the determinant through excitations from a reference might change
    // the parity
    // inefficient but easy, is there a sort routine in STL that gives me the parity too???
    auto ref0(occup);
    for (size_t i = 0; i < n; i++)
      ref0[pos[i]] = uno[i];

    for (size_t i = 0; i < ref0.size(); i++)
      for (size_t k = i + 1; k < ref0.size(); k++)
        if (ref0[i] > ref0[k])
        {
          size_t q = ref0[i];
          ref0[i]  = ref0[k];
          ref0[k]  = q;
          res *= -1.0;
        }
        else if (ref0[i] == ref0[k])
          throw std::runtime_error(" Error #2 in ci_configuration2::calculateExcitations() \n");

    return res;
  }
};

inline std::ostream& operator<<(std::ostream& out, const ci_configuration2& c)
{
  out << "ci ci_configuration2: ";
  for (size_t i = 0; i < c.occup.size(); i++)
    out << c.occup[i] << " ";
  out << std::endl;
  return out;
};

} // namespace qmcplusplus
#endif