SPM.cpp

Go to the documentation of this file.

/* 
 * @BEGIN LICENSE
 *
 * Copyright (C) 2014-2015  Ward Poelmans
 *
 * This file is part of v2DM-DOCI.
 * 
 * v2DM-DOCI is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Foobar is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
 *
 * @END LICENSE
 */

#include <algorithm>
#include "include.h"
#include "lapack.h"

#include "OptIndex.h"
#include "Irreps.h"

using namespace doci2DM;

SPM::SPM(int L, int N): BlockVector(1)
{
   this->L = L;
   this->N = N;

   setDim(0,L,2);
}

SPM::SPM(const TPM &tpm): BlockVector(1)
{
   this->L = tpm.gL();
   this->N = tpm.gN();

   setDim(0,L,2);

   bar3(1.0/(N-1.0), tpm);
}

int SPM::gN() const
{
   return N;
}

int SPM::gL() const
{
   return L;
}

int SPM::gn() const
{
   return 2*L;
}

double SPM::GetElement(int a, int b) const
{
   if(a!=b)
      return 0;

   return (*this)(0,a%L);
}

void SPM::bar(double scal, const TPM &tpm)
{
   for(int i=0;i<L;i++)
      (*this)(0,i) = scal * tpm(0,i,i);

//   int incx = L+1;
//   int incy = 1;
//   dcopy_(&L, tpm.getMatrix(0).gMatrix(), &incx, (*this)[0].gVector(), &incy);
}

void SPM::bar2(double scal, const TPM &tpm)
{
   for(int i=0;i<L;i++)
   {
      (*this)(0,i) = 0;

      for(int a=0;a<L;a++)
         (*this)(0,i) += tpm.getDiag(i,a);

      (*this)(0,i) *= scal;
   }
}

void SPM::bar3(double scal, const TPM &tpm)
{
   for(int i=0;i<L;i++)
   {
      (*this)(0,i) = tpm.getMatrix(0)(i,i);

      for(int a=0;a<L;a++)
         (*this)(0,i) += 2 * tpm.getDiag(i,a);

      (*this)(0,i) *= scal;
   }
}


namespace doci2DM {
std::ostream &operator<<(std::ostream &output,SPM &spm)
{
   output << "SPM (2x): " << std::endl;

   for(int i=0;i<spm.gL();i++)
      output << i << "\t|\t" << i << "  " << i << "\t\t" << spm(0,i) << std::endl;

   return output;
}
}

void SPM::bar(double scal, const PHM &phm)
{
   for(int i=0;i<L;i++)
   {
      (*this)(0,i) = 0;

      for(int a=0;a<L;a++)
         (*this)(0,i) += phm(i,a,i,a)+phm(i,a+L,i,a+L);

      (*this)(0,i) *= scal;
   }
}

void SPM::PrintSorted() const
{
   std::vector<std::pair<int, double>> spm_elems;
   spm_elems.reserve(L);

   for(int i=0;i<L;i++)
      spm_elems.push_back(std::make_pair(i, (*this)(0,i)));

   std::sort(spm_elems.begin(), spm_elems.end(),
         [](const std::pair<int,double> & a, const std::pair<int,double> & b) -> bool
         {
         return a.second < b.second;
         });

   for(auto& elem: spm_elems)
      std::cout << elem.first << "\t|\t" << elem.first << "  " << elem.first << "\t\t" << elem.second << std::endl;
}

std::vector<double> SPM::Particlesperirrep(const simanneal::OptIndex &index, bool print) const
{
   CheMPS2::Irreps symgroup(index.getNgroup());
   const auto orbtoirrep = index.get_irrep_each_orbital();

   std::vector<double> results(symgroup.getNumberOfIrreps(), 0);

   for(int a=0;a<L;a++)
      results[orbtoirrep[a]] += 2*(*this)(0,a);

   if(print)
   {
      std::cout << "Group: " << symgroup.getGroupName() << std::endl;
      for(int i=0;i<symgroup.getNumberOfIrreps();i++)
         std::cout << symgroup.getIrrepName(i) << ":\t" << results[i] << std::endl;
   }

   return results;
}

/*  vim: set ts=3 sw=3 expandtab :*/