HybridRepSetReader< SA > Class Template Reference

General HybridRepSetReader to handle any unitcell. More...

Inheritance diagram for HybridRepSetReader< SA >:

Collaboration diagram for HybridRepSetReader< SA >:

Public Member Functions
	HybridRepSetReader (EinsplineSetBuilder *e)
std::unique_ptr< SPOSet >	create_spline_set (const std::string &my_name, int spin, const BandInfoGroup &bandgroup) override
	create the actual spline sets More...
void	initialize_hybridrep_atomic_centers (SA &bspline) const
	initialize basic parameters of atomic orbitals More...
void	create_atomic_centers_Gspace (const Vector< std::complex< double >> &cG, Communicate &band_group_comm, const int iorb, const std::complex< double > &rotate_phase, SA &bspline) const
	initialize construct atomic orbital radial functions from plane waves More...
void	initialize_hybrid_pio_gather (const int spin, const BandInfoGroup &bandgroup, SA &bspline) const
	transforming planewave orbitals to 3D B-spline orbitals and 1D B-spline radial orbitals in real space. More...
Public Member Functions inherited from BsplineReader
	BsplineReader (EinsplineSetBuilder *e)
virtual	~BsplineReader ()
std::string	getSplineDumpFileName (const BandInfoGroup &bandgroup) const
template<typename GT , typename BCT >
bool	set_grid (const TinyVector< int, 3 > &halfg, GT *xyz_grid, BCT *xyz_bc) const
	read gvectors and set the mesh, and prepare for einspline More...
template<typename SPE >
void	check_twists (SPE &bspline, const BandInfoGroup &bandgroup) const
	initialize twist-related data for N orbitals More...
std::string	psi_g_path (int ti, int spin, int ib) const
	return the path name in hdf5 More...
std::string	psi_r_path (int ti, int spin, int ib) const
	return the path name in hdf5 More...
void	setCommon (xmlNodePtr cur)
	setting common parameters More...
std::unique_ptr< SPOSet >	create_spline_set (int spin, xmlNodePtr cur, SPOSetInputInfo &input_info)
	create the spline after one of the kind is created More...
std::unique_ptr< SPOSet >	create_spline_set (int spin, xmlNodePtr cur)
	create the spline set More...
void	setCheckNorm (bool new_checknorm)
	Set the checkNorm variable. More...
void	setRotate (bool new_rotate)
	Set the orbital rotation flag. More...
void	initialize_spo2band (int spin, const std::vector< BandInfo > &bigspace, SPOSetInfo &sposet, std::vector< int > &band2spo)
	build index tables to map a state to band with k-point folidng More...

Private Types
using	SplineReader = SplineSetReader< typename SA::SplineBase >
using	DataType = typename SA::DataType

Private Attributes
SplineReader	spline_reader_

Additional Inherited Members
Public Attributes inherited from BsplineReader
EinsplineSetBuilder *	mybuilder
	pointer to the EinsplineSetBuilder More...
Communicate *	myComm
	communicator More...
bool	checkNorm
	mesh size check the norm of orbitals More...
bool	saveSplineCoefs
	save spline coefficients to storage More...
bool	rotate
	apply orbital rotations More...
std::vector< std::vector< int > >	spo2band
	map from spo index to band index More...

Detailed Description

template<typename SA>
class qmcplusplus::HybridRepSetReader< SA >

General HybridRepSetReader to handle any unitcell.

Definition at line 30 of file HybridRepCenterOrbitals.h.

Member Typedef Documentation

DataType

using DataType = typename SA::DataType

private

Definition at line 47 of file HybridRepSetReader.h.

SplineReader

using SplineReader = SplineSetReader<typename SA::SplineBase>

private

Definition at line 46 of file HybridRepSetReader.h.

Constructor & Destructor Documentation

HybridRepSetReader()

HybridRepSetReader

(

EinsplineSetBuilder *

e

)

Definition at line 135 of file HybridRepSetReader.cpp.

                                                                 : BsplineReader(e), spline_reader_(e)
{}

Member Function Documentation

create_atomic_centers_Gspace()

void create_atomic_centers_Gspace ( const Vector< std::complex< double >> &  cG, Communicate &  band_group_comm, const int  iorb, const std::complex< double > &  rotate_phase, SA &  bspline  ) const

inline

initialize construct atomic orbital radial functions from plane waves

Definition at line 340 of file HybridRepSetReader.cpp.

References qmcplusplus::app_log(), Communicate::bcast(), bspline(), Gvectors< ST, LT >::calc_jlm_G(), Gvectors< ST, LT >::calc_phase_shift(), Gvectors< ST, LT >::calc_Ylm_G(), qmcplusplus::create(), FairDivideLow(), Communicate::getGroupID(), Communicate::isGroupLeader(), omptarget::min(), Gvectors< ST, LT >::NumGvecs, omp_get_max_threads(), omp_get_num_threads(), omp_get_thread_num(), Communicate::reduce_in_place(), Communicate::size(), and qmcplusplus::Ylm().

{
  auto& mybuilder       = spline_reader_.mybuilder;
  double rotate_phase_r = rotate_phase.real();
  double rotate_phase_i = rotate_phase.imag();

  band_group_comm.bcast(rotate_phase_r);
  band_group_comm.bcast(rotate_phase_i);
  //distribute G-vectors over processor groups
  const int Ngvecs      = mybuilder->Gvecs[0].size();
  const int Nprocs      = band_group_comm.size();
  const int Ngvecgroups = std::min(Ngvecs, Nprocs);
  Communicate gvec_group_comm(band_group_comm, Ngvecgroups);
  std::vector<int> gvec_groups(Ngvecgroups + 1, 0);
  FairDivideLow(Ngvecs, Ngvecgroups, gvec_groups);
  const int gvec_first = gvec_groups[gvec_group_comm.getGroupID()];
  const int gvec_last  = gvec_groups[gvec_group_comm.getGroupID() + 1];

  // prepare Gvecs Ylm(G)
  using UnitCellType = typename EinsplineSetBuilder::UnitCellType;
  Gvectors<double, UnitCellType> Gvecs(mybuilder->Gvecs[0], mybuilder->PrimCell, bspline.HalfG, gvec_first, gvec_last);
  // if(band_group_comm.isGroupLeader()) std::cout << "print band=" << iorb << " KE=" << Gvecs.evaluate_KE(cG) << std::endl;

  std::vector<AtomicOrbitals<DataType>>& centers = bspline.AtomicCenters;
  app_log() << "Transforming band " << iorb << " on Rank 0" << std::endl;
  // collect atomic centers by group
  std::vector<int> uniq_species;
  for (int center_idx = 0; center_idx < centers.size(); center_idx++)
  {
    auto& ACInfo         = mybuilder->AtomicCentersInfo;
    const int my_GroupID = ACInfo.GroupID[center_idx];
    int found_idx        = -1;
    for (size_t idx = 0; idx < uniq_species.size(); idx++)
      if (my_GroupID == uniq_species[idx])
      {
        found_idx = idx;
        break;
      }
    if (found_idx < 0)
      uniq_species.push_back(my_GroupID);
  }
  // construct group list
  std::vector<std::vector<int>> group_list(uniq_species.size());
  for (int center_idx = 0; center_idx < centers.size(); center_idx++)
  {
    auto& ACInfo         = mybuilder->AtomicCentersInfo;
    const int my_GroupID = ACInfo.GroupID[center_idx];
    for (size_t idx = 0; idx < uniq_species.size(); idx++)
      if (my_GroupID == uniq_species[idx])
      {
        group_list[idx].push_back(center_idx);
        break;
      }
  }

  for (int group_idx = 0; group_idx < group_list.size(); group_idx++)
  {
    const auto& mygroup        = group_list[group_idx];
    const double spline_radius = centers[mygroup[0]].getSplineRadius();
    const int spline_npoints   = centers[mygroup[0]].getSplineNpoints();
    const int lmax             = centers[mygroup[0]].getLmax();
    const double delta         = spline_radius / static_cast<double>(spline_npoints - 1);
    const int lm_tot           = (lmax + 1) * (lmax + 1);
    const size_t natoms        = mygroup.size();
    const int policy           = lm_tot > natoms ? 0 : 1;

    std::vector<std::complex<double>> i_power(lm_tot);
    // rotate phase is introduced here.
    std::complex<double> i_temp(rotate_phase_r, rotate_phase_i);
    for (size_t l = 0; l <= lmax; l++)
    {
      for (size_t lm = l * l; lm < (l + 1) * (l + 1); lm++)
        i_power[lm] = i_temp;
      i_temp *= std::complex<double>(0.0, 1.0);
    }

    std::vector<Matrix<double>> all_vals(natoms);
    std::vector<std::vector<aligned_vector<double>>> vals_local(spline_npoints * omp_get_max_threads());
    VectorSoaContainer<double, 3> myRSoA(natoms);
    for (size_t idx = 0; idx < natoms; idx++)
    {
      all_vals[idx].resize(spline_npoints, lm_tot * 2);
      all_vals[idx] = 0.0;
      myRSoA(idx)   = centers[mygroup[idx]].getCenterPos();
    }

#pragma omp parallel
    {
      const size_t tid = omp_get_thread_num();
      const size_t nt  = omp_get_num_threads();

      for (int ip = 0; ip < spline_npoints; ip++)
      {
        const size_t ip_idx = tid * spline_npoints + ip;
        if (policy == 1)
        {
          vals_local[ip_idx].resize(lm_tot * 2);
          for (size_t lm = 0; lm < lm_tot * 2; lm++)
          {
            auto& vals = vals_local[ip_idx][lm];
            vals.resize(natoms);
            std::fill(vals.begin(), vals.end(), 0.0);
          }
        }
        else
        {
          vals_local[ip_idx].resize(natoms * 2);
          for (size_t iat = 0; iat < natoms * 2; iat++)
          {
            auto& vals = vals_local[ip_idx][iat];
            vals.resize(lm_tot);
            std::fill(vals.begin(), vals.end(), 0.0);
          }
        }
      }

      const size_t size_pw_tile = 32;
      const size_t num_pw_tiles = (Gvecs.NumGvecs + size_pw_tile - 1) / size_pw_tile;
      aligned_vector<double> j_lm_G(lm_tot, 0.0);
      std::vector<aligned_vector<double>> phase_shift_r(size_pw_tile);
      std::vector<aligned_vector<double>> phase_shift_i(size_pw_tile);
      std::vector<aligned_vector<double>> YlmG(size_pw_tile);
      for (size_t ig = 0; ig < size_pw_tile; ig++)
      {
        phase_shift_r[ig].resize(natoms);
        phase_shift_i[ig].resize(natoms);
        YlmG[ig].resize(lm_tot);
      }
      SoaSphericalTensor<double> Ylm(lmax);

#pragma omp for
      for (size_t tile_id = 0; tile_id < num_pw_tiles; tile_id++)
      {
        const size_t ig_first = tile_id * size_pw_tile;
        const size_t ig_last  = std::min((tile_id + 1) * size_pw_tile, Gvecs.NumGvecs);
        for (size_t ig = ig_first; ig < ig_last; ig++)
        {
          const size_t ig_local = ig - ig_first;
          // calculate phase shift for all the centers of this group
          Gvecs.calc_phase_shift(myRSoA, ig, phase_shift_r[ig_local], phase_shift_i[ig_local]);
          Gvecs.calc_Ylm_G(ig, Ylm, YlmG[ig_local]);
        }

        for (int ip = 0; ip < spline_npoints; ip++)
        {
          double r            = delta * static_cast<double>(ip);
          const size_t ip_idx = tid * spline_npoints + ip;

          for (size_t ig = ig_first; ig < ig_last; ig++)
          {
            const size_t ig_local = ig - ig_first;
            // calculate spherical bessel function
            Gvecs.calc_jlm_G(lmax, r, ig, j_lm_G);
            for (size_t lm = 0; lm < lm_tot; lm++)
              j_lm_G[lm] *= YlmG[ig_local][lm];

            const double cG_r = cG[ig + gvec_first].real();
            const double cG_i = cG[ig + gvec_first].imag();
            if (policy == 1)
            {
              for (size_t lm = 0; lm < lm_tot; lm++)
              {
                double* restrict vals_r         = vals_local[ip_idx][lm * 2].data();
                double* restrict vals_i         = vals_local[ip_idx][lm * 2 + 1].data();
                const double* restrict ps_r_ptr = phase_shift_r[ig_local].data();
                const double* restrict ps_i_ptr = phase_shift_i[ig_local].data();
                double cG_j_r                   = cG_r * j_lm_G[lm];
                double cG_j_i                   = cG_i * j_lm_G[lm];
#pragma omp simd aligned(vals_r, vals_i, ps_r_ptr, ps_i_ptr : QMC_SIMD_ALIGNMENT)
                for (size_t idx = 0; idx < natoms; idx++)
                {
                  const double ps_r = ps_r_ptr[idx];
                  const double ps_i = ps_i_ptr[idx];
                  vals_r[idx] += cG_j_r * ps_r - cG_j_i * ps_i;
                  vals_i[idx] += cG_j_i * ps_r + cG_j_r * ps_i;
                }
              }
            }
            else
            {
              for (size_t idx = 0; idx < natoms; idx++)
              {
                double* restrict vals_r           = vals_local[ip_idx][idx * 2].data();
                double* restrict vals_i           = vals_local[ip_idx][idx * 2 + 1].data();
                const double* restrict j_lm_G_ptr = j_lm_G.data();
                double cG_ps_r = cG_r * phase_shift_r[ig_local][idx] - cG_i * phase_shift_i[ig_local][idx];
                double cG_ps_i = cG_i * phase_shift_r[ig_local][idx] + cG_r * phase_shift_i[ig_local][idx];
#pragma omp simd aligned(vals_r, vals_i, j_lm_G_ptr : QMC_SIMD_ALIGNMENT)
                for (size_t lm = 0; lm < lm_tot; lm++)
                {
                  const double jlm = j_lm_G_ptr[lm];
                  vals_r[lm] += cG_ps_r * jlm;
                  vals_i[lm] += cG_ps_i * jlm;
                }
              }
            }
          }
        }
      }

#pragma omp for collapse(2)
      for (int ip = 0; ip < spline_npoints; ip++)
        for (size_t idx = 0; idx < natoms; idx++)
        {
          double* vals = all_vals[idx][ip];
          for (size_t tid = 0; tid < nt; tid++)
            for (size_t lm = 0; lm < lm_tot; lm++)
            {
              double vals_th_r, vals_th_i;
              const size_t ip_idx = tid * spline_npoints + ip;
              if (policy == 1)
              {
                vals_th_r = vals_local[ip_idx][lm * 2][idx];
                vals_th_i = vals_local[ip_idx][lm * 2 + 1][idx];
              }
              else
              {
                vals_th_r = vals_local[ip_idx][idx * 2][lm];
                vals_th_i = vals_local[ip_idx][idx * 2 + 1][lm];
              }
              const double real_tmp = 4.0 * M_PI * i_power[lm].real();
              const double imag_tmp = 4.0 * M_PI * i_power[lm].imag();
              vals[lm] += vals_th_r * real_tmp - vals_th_i * imag_tmp;
              vals[lm + lm_tot] += vals_th_i * real_tmp + vals_th_r * imag_tmp;
            }
        }
    }
    //app_log() << "Building band " << iorb << " at center " << center_idx << std::endl;

    for (size_t idx = 0; idx < natoms; idx++)
    {
      // reduce all_vals
      band_group_comm.reduce_in_place(all_vals[idx].data(), all_vals[idx].size());
      if (!band_group_comm.isGroupLeader())
        continue;
#pragma omp parallel for
      for (int lm = 0; lm < lm_tot; lm++)
      {
        auto& mycenter = centers[mygroup[idx]];
        aligned_vector<double> splineData_r(spline_npoints);
        UBspline_1d_d* atomic_spline_r = nullptr;
        for (size_t ip = 0; ip < spline_npoints; ip++)
          splineData_r[ip] = all_vals[idx][ip][lm];
        atomic_spline_r = einspline::create(atomic_spline_r, 0.0, spline_radius, spline_npoints, splineData_r.data(),
                                            ((lm == 0) || (lm > 3)));
        if (!bspline.isComplex())
        {
          mycenter.set_spline(atomic_spline_r, lm, iorb);
          einspline::destroy(atomic_spline_r);
        }
        else
        {
          aligned_vector<double> splineData_i(spline_npoints);
          UBspline_1d_d* atomic_spline_i = nullptr;
          for (size_t ip = 0; ip < spline_npoints; ip++)
            splineData_i[ip] = all_vals[idx][ip][lm + lm_tot];
          atomic_spline_i = einspline::create(atomic_spline_i, 0.0, spline_radius, spline_npoints, splineData_i.data(),
                                              ((lm == 0) || (lm > 3)));
          mycenter.set_spline(atomic_spline_r, lm, iorb * 2);
          mycenter.set_spline(atomic_spline_i, lm, iorb * 2 + 1);
          einspline::destroy(atomic_spline_r);
          einspline::destroy(atomic_spline_i);
        }
      }
    }
  }
}

create_spline_set()

std::unique_ptr< SPOSet > create_spline_set ( const std::string &  my_name, int  spin, const BandInfoGroup &  bandgroup  )

overridevirtual

create the actual spline sets

Implements BsplineReader.

Definition at line 139 of file HybridRepSetReader.cpp.

References qmcplusplus::app_log(), bspline(), hdf_archive::close(), hdf_archive::create(), Timer::elapsed(), hdf_archive::open(), and hdf_archive::write().

{
  auto bspline = std::make_unique<SA>(my_name);
  app_log() << "  ClassName = " << bspline->getClassName() << std::endl;
  // set info for Hybrid
  initialize_hybridrep_atomic_centers(*bspline);
  bool foundspline = spline_reader_.createSplineDataSpaceLookforDumpFile(bandgroup, *bspline);
  typename SA::HYBRIDBASE& hybrid_center_orbs = *bspline;
  hybrid_center_orbs.resizeStorage(bspline->myV.size());
  if (foundspline)
  {
    Timer now;
    hdf_archive h5f(myComm);
    const auto splinefile = getSplineDumpFileName(bandgroup);
    h5f.open(splinefile, H5F_ACC_RDONLY);
    foundspline = bspline->read_splines(h5f);
    if (foundspline)
      app_log() << "  Successfully restored 3D B-spline coefficients from " << splinefile << ". The reading time is "
                << now.elapsed() << " sec." << std::endl;
  }

  if (!foundspline)
  {
    hybrid_center_orbs.flush_zero();
    initialize_hybrid_pio_gather(spin, bandgroup, *bspline);

    if (saveSplineCoefs && myComm->rank() == 0)
    {
      Timer now;
      const std::string splinefile(getSplineDumpFileName(bandgroup));
      hdf_archive h5f;
      h5f.create(splinefile);
      std::string classname = bspline->getClassName();
      h5f.write(classname, "class_name");
      int sizeD = sizeof(DataType);
      h5f.write(sizeD, "sizeof");
      bspline->write_splines(h5f);
      h5f.close();
      app_log() << "  Stored spline coefficients in " << splinefile << " for potential reuse. The writing time is "
                << now.elapsed() << " sec." << std::endl;
    }
  }

  {
    Timer now;
    bspline->bcast_tables(myComm);
    app_log() << "  Time to bcast the table = " << now.elapsed() << std::endl;
  }
  return bspline;
}

initialize_hybrid_pio_gather()

void initialize_hybrid_pio_gather	(	const int	spin,
		const BandInfoGroup &	bandgroup,
		SA &	bspline
	)		const

transforming planewave orbitals to 3D B-spline orbitals and 1D B-spline radial orbitals in real space.

Parameters

spin	orbital dataset spin index
bandgroup	band info
bspline	the spline object being worked on

Definition at line 613 of file HybridRepSetReader.cpp.

References qmcplusplus::app_log(), Communicate::bcast(), bspline(), Timer::elapsed(), FairDivideLow(), Communicate::getGroupID(), Communicate::getGroupLeaderComm(), BandInfoGroup::getNumDistinctOrbitals(), Communicate::isGroupLeader(), omptarget::min(), and BandInfoGroup::myBands.

{
  auto& mybuilder = spline_reader_.mybuilder;
  //distribute bands over processor groups
  int Nbands            = bandgroup.getNumDistinctOrbitals();
  const int Nprocs      = myComm->size();
  const int Nbandgroups = std::min(Nbands, Nprocs);
  Communicate band_group_comm(*myComm, Nbandgroups);
  std::vector<int> band_groups(Nbandgroups + 1, 0);
  FairDivideLow(Nbands, Nbandgroups, band_groups);
  int iorb_first = band_groups[band_group_comm.getGroupID()];
  int iorb_last  = band_groups[band_group_comm.getGroupID() + 1];

  app_log() << "Start transforming plane waves to 3D B-splines and atomic radial orbital 1D B-splines." << std::endl;
  OneSplineOrbData oneband(mybuilder->MeshSize, bspline.HalfG, bspline.isComplex());
  hdf_archive h5f(&band_group_comm, false);
  Vector<std::complex<double>> cG(mybuilder->Gvecs[0].size());
  const std::vector<BandInfo>& cur_bands = bandgroup.myBands;
  if (band_group_comm.isGroupLeader())
    h5f.open(mybuilder->H5FileName, H5F_ACC_RDONLY);
  for (int iorb = iorb_first; iorb < iorb_last; iorb++)
  {
    if (band_group_comm.isGroupLeader())
    {
      const auto& cur_band = cur_bands[bspline.BandIndexMap[iorb]];
      const int ti         = cur_band.TwistIndex;
      spline_reader_.readOneOrbitalCoefs(psi_g_path(ti, spin, cur_band.BandIndex), h5f, cG);
      oneband.fft_spline(cG, mybuilder->Gvecs[0], mybuilder->primcell_kpoints[ti], rotate);
      bspline.set_spline(&oneband.get_spline_r(), &oneband.get_spline_i(), cur_band.TwistIndex, iorb, 0);
    }
    band_group_comm.bcast(cG);
    create_atomic_centers_Gspace(cG, band_group_comm, iorb, oneband.getRotatePhase(), bspline);
  }

  myComm->barrier();
  if (band_group_comm.isGroupLeader())
  {
    Timer now;
    bspline.gather_tables(band_group_comm.getGroupLeaderComm());
    app_log() << "  Time to gather the table = " << now.elapsed() << std::endl;
  }
}

initialize_hybridrep_atomic_centers()

void initialize_hybridrep_atomic_centers

(

SA &

bspline

)

const

initialize basic parameters of atomic orbitals

Definition at line 193 of file HybridRepSetReader.cpp.

References OhmmsAttributeSet::add(), APP_ABORT, qmcplusplus::app_error(), qmcplusplus::app_log(), bspline(), qmcplusplus::ceil(), qmcplusplus::COSCOS, qmcplusplus::Units::charge::e, qmcplusplus::LEKS2018, qmcplusplus::LINEAR, omptarget::min(), OhmmsAttributeSet::put(), and AtomicOrbitals< ST >::set_info().

{
  auto& mybuilder = spline_reader_.mybuilder;

  OhmmsAttributeSet a;
  std::string scheme_name("Consistent");
  std::string s_function_name("LEKS2018");
  a.add(scheme_name, "smoothing_scheme");
  a.add(s_function_name, "smoothing_function");
  a.put(mybuilder->XMLRoot);
  // assign smooth_scheme
  if (scheme_name == "Consistent")
    bspline.smooth_scheme = SA::smoothing_schemes::CONSISTENT;
  else if (scheme_name == "SmoothAll")
    bspline.smooth_scheme = SA::smoothing_schemes::SMOOTHALL;
  else if (scheme_name == "SmoothPartial")
    bspline.smooth_scheme = SA::smoothing_schemes::SMOOTHPARTIAL;
  else
    APP_ABORT("initialize_hybridrep_atomic_centers wrong smoothing_scheme name! Only allows Consistent, SmoothAll or "
              "SmoothPartial.");

  // assign smooth_function
  if (s_function_name == "LEKS2018")
    bspline.smooth_func_id = smoothing_functions::LEKS2018;
  else if (s_function_name == "coscos")
    bspline.smooth_func_id = smoothing_functions::COSCOS;
  else if (s_function_name == "linear")
    bspline.smooth_func_id = smoothing_functions::LINEAR;
  else
    APP_ABORT(
        "initialize_hybridrep_atomic_centers wrong smoothing_function name! Only allows LEKS2018, coscos or linear.");
  app_log() << "Hybrid orbital representation uses " << scheme_name << " smoothing scheme and " << s_function_name
            << " smoothing function." << std::endl;

  bspline.set_info(*(mybuilder->SourcePtcl), mybuilder->TargetPtcl, mybuilder->Super2Prim);
  auto& centers = bspline.AtomicCenters;
  auto& ACInfo  = mybuilder->AtomicCentersInfo;
  // load atomic center info only when it is not initialized
  if (centers.size() == 0)
  {
    bool success = true;
    app_log() << "Reading atomic center info for hybrid representation" << std::endl;
    for (int center_idx = 0; center_idx < ACInfo.Ncenters; center_idx++)
    {
      const int my_GroupID = ACInfo.GroupID[center_idx];
      if (ACInfo.cutoff[center_idx] < 0)
      {
        app_error() << "Hybrid orbital representation needs parameter 'cutoff_radius' for atom " << center_idx
                    << std::endl;
        success = false;
      }

      if (ACInfo.inner_cutoff[center_idx] < 0)
      {
        const double inner_cutoff = std::max(ACInfo.cutoff[center_idx] - 0.3, 0.0);
        app_log() << "Hybrid orbital representation setting 'inner_cutoff' to " << inner_cutoff << " for group "
                  << my_GroupID << " as atom " << center_idx << std::endl;
        // overwrite the inner_cutoff of all the atoms of the same species
        for (int id = 0; id < ACInfo.Ncenters; id++)
          if (my_GroupID == ACInfo.GroupID[id])
            ACInfo.inner_cutoff[id] = inner_cutoff;
      }
      else if (ACInfo.inner_cutoff[center_idx] > ACInfo.cutoff[center_idx])
      {
        app_error() << "Hybrid orbital representation 'inner_cutoff' must be smaller than 'spline_radius' for atom "
                    << center_idx << std::endl;
        success = false;
      }

      if (ACInfo.cutoff[center_idx] > 0)
      {
        if (ACInfo.lmax[center_idx] < 0)
        {
          app_error() << "Hybrid orbital representation needs parameter 'lmax' for atom " << center_idx << std::endl;
          success = false;
        }

        if (ACInfo.spline_radius[center_idx] < 0 && ACInfo.spline_npoints[center_idx] < 0)
        {
          app_log() << "Parameters 'spline_radius' and 'spline_npoints' for group " << my_GroupID << " as atom "
                    << center_idx << " are not specified." << std::endl;
          const double delta     = std::min(0.02, ACInfo.cutoff[center_idx] / 4.0);
          const int n_grid_point = std::ceil((ACInfo.cutoff[center_idx] + 1e-4) / delta) + 3;
          for (int id = 0; id < ACInfo.Ncenters; id++)
            if (my_GroupID == ACInfo.GroupID[id])
            {
              ACInfo.spline_npoints[id] = n_grid_point;
              ACInfo.spline_radius[id]  = (n_grid_point - 1) * delta;
            }
          app_log() << "  Based on default grid point distance " << delta << std::endl;
          app_log() << "  Setting 'spline_npoints' to " << ACInfo.spline_npoints[center_idx] << std::endl;
          app_log() << "  Setting 'spline_radius' to " << ACInfo.spline_radius[center_idx] << std::endl;
        }
        else
        {
          if (ACInfo.spline_radius[center_idx] < 0)
          {
            app_error() << "Hybrid orbital representation needs parameter 'spline_radius' for atom " << center_idx
                        << std::endl;
            success = false;
          }

          if (ACInfo.spline_npoints[center_idx] < 0)
          {
            app_error() << "Hybrid orbital representation needs parameter 'spline_npoints' for atom " << center_idx
                        << std::endl;
            success = false;
          }
        }

        // check maximally allowed cutoff_radius
        double max_allowed_cutoff = ACInfo.spline_radius[center_idx] -
            2.0 * ACInfo.spline_radius[center_idx] / (ACInfo.spline_npoints[center_idx] - 1);
        if (success && ACInfo.cutoff[center_idx] > max_allowed_cutoff)
        {
          app_error() << "Hybrid orbital representation requires cutoff_radius<=" << max_allowed_cutoff
                      << " calculated by spline_radius-2*spline_radius/(spline_npoints-1) for atom " << center_idx
                      << std::endl;
          success = false;
        }
      }
      else
      {
        // no atomic regions for this atom type
        ACInfo.spline_radius[center_idx]  = 0.0;
        ACInfo.spline_npoints[center_idx] = 0;
        ACInfo.lmax[center_idx]           = 0;
      }
    }

    if (!success)
      spline_reader_.myComm->barrier_and_abort(
          "initialize_hybridrep_atomic_centers Failed to initialize atomic centers "
          "in hybrid orbital representation!");

    for (int center_idx = 0; center_idx < ACInfo.Ncenters; center_idx++)
    {
      AtomicOrbitals<DataType> oneCenter(ACInfo.lmax[center_idx]);
      oneCenter.set_info(ACInfo.ion_pos[center_idx], ACInfo.cutoff[center_idx], ACInfo.inner_cutoff[center_idx],
                         ACInfo.spline_radius[center_idx], ACInfo.non_overlapping_radius[center_idx],
                         ACInfo.spline_npoints[center_idx]);
      centers.push_back(oneCenter);
    }
  }
}

Member Data Documentation

spline_reader_

SplineReader spline_reader_

private

Definition at line 48 of file HybridRepSetReader.h.

---

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

• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCWaveFunctions/BsplineFactory/HybridRepCenterOrbitals.h
• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCWaveFunctions/BsplineFactory/HybridRepSetReader.h
• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCWaveFunctions/BsplineFactory/HybridRepSetReader.cpp