KContainer Class Reference

Container for k-points. More...

Inheritance diagram for KContainer:

Collaboration diagram for KContainer:

Public Types
using	ParticleLayout = PtclOnLatticeTraits::ParticleLayout
Public Types inherited from QMCTraits
enum	{ DIM = OHMMS_DIM, DIM_VGL = OHMMS_DIM + 2 }
using	QTBase = QMCTypes< OHMMS_PRECISION, DIM >
using	QTFull = QMCTypes< OHMMS_PRECISION_FULL, DIM >
using	RealType = QTBase::RealType
using	ComplexType = QTBase::ComplexType
using	ValueType = QTBase::ValueType
using	PosType = QTBase::PosType
using	GradType = QTBase::GradType
using	TensorType = QTBase::TensorType
using	IndexType = OHMMS_INDEXTYPE
	define other types More...
using	FullPrecRealType = QTFull::RealType
using	FullPrecValueType = QTFull::ValueType
using	PropertySetType = RecordNamedProperty< FullPrecRealType >
	define PropertyList_t More...
using	PtclGrpIndexes = std::vector< std::pair< int, int > >

Public Member Functions
void	updateKLists (const ParticleLayout &lattice, RealType kc, unsigned ndim, const PosType &twist=PosType(), bool useSphere=true)
	k points sorted by the |k| excluding |k|=0 More...
const auto &	get_kpts_cart_soa () const

Public Attributes
int	numk
	number of k-points More...
TinyVector< int, DIM+1 >	mmax
	maximum integer translations of reciprocal cell within kc. More...
std::vector< TinyVector< int, DIM > >	kpts
	K-vector in reduced coordinates. More...
std::vector< PosType >	kpts_cart
	K-vector in Cartesian coordinates. More...
std::vector< RealType >	ksq
	squre of kpts in Cartesian coordniates More...
std::vector< int >	minusk
	Given a k index, return index to -k. More...
std::vector< int >	kshell
	kpts which belong to the ith-shell [kshell[i], kshell[i+1]) More...

Private Member Functions
void	findApproxMMax (const ParticleLayout &lattice, unsigned ndim)
	compute approximate parallelpiped that surrounds kc More...
void	BuildKLists (const ParticleLayout &lattice, const PosType &twist, bool useSphere)
	construct the container for k-vectors More...

Private Attributes
RealType	kcutoff
	The cutoff up to which k-vectors are generated. More...
VectorSoaContainer< RealType, DIM, OffloadAllocator< RealType > >	kpts_cart_soa_
	K-vector in Cartesian coordinates in SoA layout. More...

Detailed Description

Container for k-points.

It generates a set of k-points that are unit-translations of the reciprocal-space cell. K-points are generated within a spherical cutoff set by the supercell

Definition at line 29 of file KContainer.h.

Member Typedef Documentation

ParticleLayout

using ParticleLayout = PtclOnLatticeTraits::ParticleLayout

Definition at line 37 of file KContainer.h.

Member Function Documentation

BuildKLists()

void BuildKLists ( const ParticleLayout &  lattice, const PosType &  twist, bool  useSphere  )

private

construct the container for k-vectors

Definition at line 105 of file KContainer.cpp.

References qmcplusplus::abs(), QMCTraits::DIM, qmcplusplus::dot(), KContainer::kcutoff, KContainer::kpts, KContainer::kpts_cart, KContainer::kpts_cart_soa_, KContainer::kshell, KContainer::ksq, qmcplusplus::lattice, KContainer::minusk, KContainer::mmax, and KContainer::numk.

Referenced by KContainer::updateKLists().

{
  TinyVector<int, DIM + 1> TempActualMax;
  TinyVector<int, DIM> kvec;
  TinyVector<RealType, DIM> kvec_cart;
  RealType modk2;
  std::vector<TinyVector<int, DIM>> kpts_tmp;
  std::vector<PosType> kpts_cart_tmp;
  std::vector<RealType> ksq_tmp;
  // reserve the space for memory efficiency
  if (useSphere)
  {
    const RealType kcut2 = kcutoff * kcutoff;
    //Loop over guesses for valid k-points.
    for (int i = -mmax[0]; i <= mmax[0]; i++)
    {
      kvec[0] = i;
      for (int j = -mmax[1]; j <= mmax[1]; j++)
      {
        kvec[1] = j;
        for (int k = -mmax[2]; k <= mmax[2]; k++)
        {
          kvec[2] = k;
          //Do not include k=0 in evaluations.
          if (i == 0 && j == 0 && k == 0)
            continue;
          //Convert kvec to Cartesian
          kvec_cart = lattice.k_cart(kvec + twist);
          //Find modk
          modk2 = dot(kvec_cart, kvec_cart);
          if (modk2 > kcut2)
            continue; //Inside cutoff?
          //This k-point should be added to the list
          kpts_tmp.push_back(kvec);
          kpts_cart_tmp.push_back(kvec_cart);
          ksq_tmp.push_back(modk2);
          //Update record of the allowed maximum translation.
          for (int idim = 0; idim < 3; idim++)
            if (std::abs(kvec[idim]) > TempActualMax[idim])
              TempActualMax[idim] = std::abs(kvec[idim]);
        }
      }
    }
  }
  else
  {
    // Loop over all k-points in the parallelpiped and add them to kcontainer
    // note layout is for interfacing with fft, so for each dimension, the
    // positive indexes come first then the negative indexes backwards
    // e.g.    0, 1, .... mmax, -mmax+1, -mmax+2, ... -1
    const int idimsize = mmax[0] * 2;
    const int jdimsize = mmax[1] * 2;
    const int kdimsize = mmax[2] * 2;
    for (int i = 0; i < idimsize; i++)
    {
      kvec[0] = i;
      if (kvec[0] > mmax[0])
        kvec[0] -= idimsize;
      for (int j = 0; j < jdimsize; j++)
      {
        kvec[1] = j;
        if (kvec[1] > mmax[1])
          kvec[1] -= jdimsize;
        for (int k = 0; k < kdimsize; k++)
        {
          kvec[2] = k;
          if (kvec[2] > mmax[2])
            kvec[2] -= kdimsize;
          // get cartesian location and modk2
          kvec_cart = lattice.k_cart(kvec);
          modk2     = dot(kvec_cart, kvec_cart);
          // add k-point to lists
          kpts_tmp.push_back(kvec);
          kpts_cart_tmp.push_back(kvec_cart);
          ksq_tmp.push_back(modk2);
        }
      }
    }
    // set allowed maximum translation
    TempActualMax[0] = mmax[0];
    TempActualMax[1] = mmax[1];
    TempActualMax[2] = mmax[2];
  }

  //Update a record of the number of k vectors
  numk = kpts_tmp.size();
  std::map<int64_t, std::vector<int>*> kpts_sorted;
  //create the map: use simple integer with resolution of 0.00000001 in ksq
  for (int ik = 0; ik < numk; ik++)
  {
    //This is a workaround for ewald bug (Issue #2105).  Basically, 1e-7 is the resolution of |k|^2 for doubles,
    //so we jack up the tolerance to match that.
    const int64_t k_ind = static_cast<int64_t>(ksq_tmp[ik] * 10000000);
    auto it(kpts_sorted.find(k_ind));
    if (it == kpts_sorted.end())
    {
      std::vector<int>* newSet = new std::vector<int>;
      kpts_sorted[k_ind]       = newSet;
      newSet->push_back(ik);
    }
    else
    {
      (*it).second->push_back(ik);
    }
  }
  std::map<int64_t, std::vector<int>*>::iterator it(kpts_sorted.begin());
  kpts.resize(numk);
  kpts_cart.resize(numk);
  kpts_cart_soa_.resize(numk);
  ksq.resize(numk);
  kshell.resize(kpts_sorted.size() + 1, 0);
  int ok = 0, ish = 0;
  while (it != kpts_sorted.end())
  {
    std::vector<int>::iterator vit((*it).second->begin());
    while (vit != (*it).second->end())
    {
      int ik             = (*vit);
      kpts[ok]           = kpts_tmp[ik];
      kpts_cart[ok]      = kpts_cart_tmp[ik];
      kpts_cart_soa_(ok) = kpts_cart_tmp[ik];
      ksq[ok]            = ksq_tmp[ik];
      ++vit;
      ++ok;
    }
    kshell[ish + 1] = kshell[ish] + (*it).second->size();
    ++it;
    ++ish;
  }
  kpts_cart_soa_.updateTo();
  it = kpts_sorted.begin();
  std::map<int64_t, std::vector<int>*>::iterator e_it(kpts_sorted.end());
  while (it != e_it)
  {
    delete it->second;
    it++;
  }
  //Finished searching k-points. Copy list of maximum translations.
  mmax[DIM] = 0;
  for (int idim = 0; idim < DIM; idim++)
  {
    mmax[idim] = TempActualMax[idim];
    mmax[DIM]  = std::max(mmax[idim], mmax[DIM]);
    //if(mmax[idim] > mmax[DIM]) mmax[DIM] = mmax[idim];
  }
  //Now fill the array that returns the index of -k when given the index of k.
  minusk.resize(numk);

  //Assigns a unique hash value to each kpoint.
  auto getHashOfVec = [](const auto& inpv, int hashparam) -> int64_t {
    int64_t hash = 0; // this will cause integral promotion below
    for (int i = 0; i < inpv.Size; ++i)
      hash += inpv[i] + hash * hashparam;
    return hash;
  };

  // Create a map from the hash value for each k vector to the index
  std::map<int64_t, int> hashToIndex;
  for (int ki = 0; ki < numk; ki++)
  {
    hashToIndex[getHashOfVec(kpts[ki], numk)] = ki;
  }
  // Use the map to find the index of -k from the index of k
  for (int ki = 0; ki < numk; ki++)
  {
    minusk[ki] = hashToIndex[getHashOfVec(-1 * kpts[ki], numk)];
  }
}

findApproxMMax()

void findApproxMMax ( const ParticleLayout &  lattice, unsigned  ndim  )

private

compute approximate parallelpiped that surrounds kc

Parameters

lattice

supercell

Definition at line 44 of file KContainer.cpp.

References qmcplusplus::app_log(), QMCTraits::DIM, qmcplusplus::dot(), qmcplusplus::floor(), LRCoulombSingleton::isQuasi2D(), KContainer::kcutoff, qmcplusplus::lattice, KContainer::mmax, and qmcplusplus::sqrt().

Referenced by KContainer::updateKLists().

{
  //Estimate the size of the parallelpiped that encompasses a sphere of kcutoff.
  //mmax is stored as integer translations of the reciprocal cell vectors.
  //Does not require an orthorhombic cell.
  /* Old method.
  //2pi is not included in lattice.b
  Matrix<RealType> mmat;
  mmat.resize(3,3);
  for(int j=0;j<3;j++)
    for(int i=0;i<3;i++){
      mmat[i][j] = 0.0;
      for(int k=0;k<3;k++)
  mmat[i][j] = mmat[i][j] + 4.0*M_PI*M_PI*lattice.b(k)[i]*lattice.b(j)[k];
    }

  TinyVector<RealType,3> x,temp;
  RealType tempr;
  for(int idim=0;idim<3;idim++){
    int i = ((idim)%3);
    int j = ((idim+1)%3);
    int k = ((idim+2)%3);

    x[i] = 1.0;
    x[j] = (mmat[j][k]*mmat[k][i] - mmat[k][k]*mmat[i][j]);
    x[j]/= (mmat[j][j]*mmat[k][k] - mmat[j][k]*mmat[j][k]);
    x[k] = -(mmat[k][i] + mmat[j][k]*x[j])/mmat[k][k];

    for(i=0;i<3;i++){
      temp[i] = 0.0;
  for(j=0;j<3;j++)
    temp[i] += mmat[i][j]*x[j];
    }

    tempr = dot(x,temp);
    mmax[idim] = static_cast<int>(sqrt(4.0*kcut2/tempr)) + 1;
  }
  */
  // see rmm, Electronic Structure, p. 85 for details
  for (int i = 0; i < DIM; i++)
    mmax[i] = static_cast<int>(std::floor(std::sqrt(dot(lattice.a(i), lattice.a(i))) * kcutoff / (2 * M_PI))) + 1;

  mmax[DIM] = mmax[0];
  for (int i = 1; i < DIM; ++i)
    mmax[DIM] = std::max(mmax[i], mmax[DIM]);

  //overwrite the non-periodic directon to be zero
  if (LRCoulombSingleton::isQuasi2D())
  {
    app_log() << "  No kspace sum perpendicular to slab " << std::endl;
    mmax[2] = 0;
  }
  if (ndim < 3)
  {
    app_log() << "  No kspace sum along z " << std::endl;
    mmax[2] = 0;
  }
  if (ndim < 2)
    mmax[1] = 0;
}

get_kpts_cart_soa()

const auto& get_kpts_cart_soa ( ) const

inline

Definition at line 82 of file KContainer.h.

References KContainer::kpts_cart_soa_.

{ return kpts_cart_soa_; }

updateKLists()

void updateKLists	(	const ParticleLayout &	lattice,
		RealType	kc,
		unsigned	ndim,
		const PosType &	twist = PosType(),
		bool	useSphere = true
	)

k points sorted by the |k| excluding |k|=0

The first for |k| The second for a map to the full index. The size of the second is the degeneracy. update k-vectors

Parameters

sc	supercell
kc	cutoff radius in the K
twist	shifts the center of the grid of k-vectors
useSphere	if true, use the |K|

Definition at line 24 of file KContainer.cpp.

References APP_ABORT, qmcplusplus::app_log(), KContainer::BuildKLists(), KContainer::findApproxMMax(), KContainer::kcutoff, KContainer::kpts, KContainer::kshell, and qmcplusplus::lattice.

Referenced by FreeOrbitalBuilder::createSPOSetFromXML(), SimulationCell::resetLRBox(), LRBreakup< BreakupBasis >::SetupKVecs(), and qmcplusplus::TEST_CASE().

{
  kcutoff = kc;
  if (kcutoff <= 0.0)
  {
    APP_ABORT("  Illegal cutoff for KContainer");
  }
  findApproxMMax(lattice, ndim);
  BuildKLists(lattice, twist, useSphere);

  app_log() << "  KContainer initialised with cutoff " << kcutoff << std::endl;
  app_log() << "   # of K-shell  = " << kshell.size() << std::endl;
  app_log() << "   # of K points = " << kpts.size() << std::endl;
  app_log() << std::endl;
}

Member Data Documentation

kcutoff

RealType kcutoff

private

The cutoff up to which k-vectors are generated.

Definition at line 33 of file KContainer.h.

Referenced by KContainer::BuildKLists(), KContainer::findApproxMMax(), and KContainer::updateKLists().

kpts

std::vector<TinyVector<int, DIM> > kpts

K-vector in reduced coordinates.

Definition at line 50 of file KContainer.h.

Referenced by KContainer::BuildKLists(), FreeOrbitalBuilder::createSPOSetFromXML(), QMCFiniteSize::initialize(), qmcplusplus::TEST_CASE(), and KContainer::updateKLists().

kpts_cart

std::vector<PosType> kpts_cart

K-vector in Cartesian coordinates.

Definition at line 53 of file KContainer.h.

Referenced by KContainer::BuildKLists(), StructFact::computeRhok(), FreeOrbitalBuilder::createSPOSetFromXML(), EwaldHandler3D::filldFk_dk(), LRHandlerSRCoulomb< Func, BreakupBasis >::filldFk_dk(), EwaldHandler2D::fillFk(), EwaldHandlerQuasi2D::fillFk(), EwaldHandler3D::fillFk(), LRRPABFeeHandlerTemp< Func, BreakupBasis >::fillFk(), LRRPAHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkg(), EwaldHandler3D::fillYkgstrain(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkgstrain(), and qmcplusplus::TEST_CASE().

kpts_cart_soa_

VectorSoaContainer<RealType, DIM, OffloadAllocator<RealType> > kpts_cart_soa_

private

K-vector in Cartesian coordinates in SoA layout.

Definition at line 94 of file KContainer.h.

Referenced by KContainer::BuildKLists(), and KContainer::get_kpts_cart_soa().

kshell

std::vector<int> kshell

kpts which belong to the ith-shell [kshell[i], kshell[i+1])

Definition at line 61 of file KContainer.h.

Referenced by KContainer::BuildKLists(), QMCFiniteSize::calcPotentialDiscrete(), QMCFiniteSize::calcPotentialInt(), EwaldHandlerQuasi2D::fillFk(), EwaldHandler2D::fillFk(), EwaldHandler3D::fillFk(), LRRPABFeeHandlerTemp< Func, BreakupBasis >::fillFk(), LRRPAHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkg(), EwaldHandler3D::fillYkgstrain(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkgstrain(), QMCFiniteSize::printSkRawSphAvg(), LRBreakup< BreakupBasis >::SetupKVecs(), and KContainer::updateKLists().

ksq

std::vector<RealType> ksq

squre of kpts in Cartesian coordniates

Definition at line 56 of file KContainer.h.

Referenced by KContainer::BuildKLists(), EwaldHandler3D::filldFk_dk(), LRHandlerSRCoulomb< Func, BreakupBasis >::filldFk_dk(), EwaldHandler2D::fillFk(), EwaldHandlerQuasi2D::fillFk(), EwaldHandler3D::fillFk(), LRRPABFeeHandlerTemp< Func, BreakupBasis >::fillFk(), LRRPAHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerTemp< Func, BreakupBasis >::fillFk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYk(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkg(), EwaldHandler3D::fillYkgstrain(), LRHandlerSRCoulomb< Func, BreakupBasis >::fillYkgstrain(), QMCFiniteSize::printSkRawSphAvg(), and LRBreakup< BreakupBasis >::SetupKVecs().

minusk

std::vector<int> minusk

Given a k index, return index to -k.

Definition at line 59 of file KContainer.h.

Referenced by KContainer::BuildKLists(), and FreeOrbitalBuilder::createSPOSetFromXML().

mmax

TinyVector<int, DIM + 1> mmax

maximum integer translations of reciprocal cell within kc.

Last index is max. of first dimension+1

Definition at line 46 of file KContainer.h.

Referenced by KContainer::BuildKLists(), and KContainer::findApproxMMax().

numk

int numk

number of k-points

Definition at line 40 of file KContainer.h.

Referenced by KContainer::BuildKLists(), StructFact::computeRhok(), LRBreakup< BreakupBasis >::SetupKVecs(), and qmcplusplus::TEST_CASE().

---

The documentation for this class was generated from the following files:

• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/Particle/LongRange/KContainer.h
• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/Particle/LongRange/KContainer.cpp