ergo
grid_hicu.cc File Reference

Hierarchical Cubature (HiCu) grid generation. More...

#include <stdlib.h>
#include <cmath>
#include <stdio.h>
#include <errno.h>
#include <memory.h>
#include <time.h>
#include <pthread.h>
#include <stdexcept>
#include "grid_hicu.h"
#include "basisinfo.h"
#include "integrals_general.h"
#include "cubature_rules.h"
#include "utilities.h"
#include "pi.h"
#include "box_system.h"
#include "integrator.h"
#include "functionals.h"
#include "aos.h"
#include "dft_common.h"
#include "rho-mat.h"
#include "units.h"

Classes

struct  ShellSpecStructWithExtent
 
struct  DensitySpecStruct
 
struct  rhoTreeNode_
 
struct  GridGenerationParamsStruct
 
struct  compute_grid_for_box_params_struct
 
struct  ComputeGridResultValuesStruct
 
struct  compute_grid_thread_func_struct
 

Macros

#define __CVERSION__
 
#define USE_EXP_STD
 
#define USE_ERF_STD
 
#define DO_EXTRA_ERROR_CHECKING
 
#define FILE_BATCH_N   1000000
 
#define MAX_NO_OF_POINTS_PER_BATCH   100
 
#define MAX_NO_OF_SHLBLOCKS   44444
 
#define EXPONENT_DIFF_LIMIT   1e-22
 
#define DISTR_CENTER_DIST_LIMIT   1e-22
 
#define N_BATCH_JOBS   22
 
#define MAX_NO_OF_POINTS_PER_WRITE   50000
 
#define HICU_SPARSE_MATRIX_ACCESS_ROUTINE   at
 
#define solid_harmonic_s_0(x, y, z, x2, y2, z2, r2)   1
 
#define solid_harmonic_p_2(x, y, z, x2, y2, z2, r2)   x
 
#define solid_harmonic_p_0(x, y, z, x2, y2, z2, r2)   y
 
#define solid_harmonic_p_1(x, y, z, x2, y2, z2, r2)   z
 
#define solid_harmonic_d_0(x, y, z, x2, y2, z2, r2)   (x * y)
 
#define solid_harmonic_d_1(x, y, z, x2, y2, z2, r2)   (y * z)
 
#define solid_harmonic_d_2(x, y, z, x2, y2, z2, r2)   ((2 * z2 - x2 - y2) / (2 * template_blas_sqrt((ergo_real)3)))
 
#define solid_harmonic_d_3(x, y, z, x2, y2, z2, r2)   (x * z)
 
#define solid_harmonic_d_4(x, y, z, x2, y2, z2, r2)   (0.5 * (x2 - y2))
 
#define solid_harmonic_f_0(x, y, z, x2, y2, z2, r2)   ((0.5 * template_blas_sqrt(2.5) * (3 * x2 - y2) * y) / template_blas_sqrt((ergo_real)15))
 
#define solid_harmonic_f_1(x, y, z, x2, y2, z2, r2)   (x * y * z)
 
#define solid_harmonic_f_2(x, y, z, x2, y2, z2, r2)   (0.5 * template_blas_sqrt((ergo_real)1.5) * (5 * z2 - r2) * y / template_blas_sqrt((ergo_real)15))
 
#define solid_harmonic_f_3(x, y, z, x2, y2, z2, r2)   (0.5 * (5 * z2 - 3 * r2) * z / template_blas_sqrt((ergo_real)15))
 
#define solid_harmonic_f_4(x, y, z, x2, y2, z2, r2)   (0.5 * template_blas_sqrt((ergo_real)1.5) * (5 * z2 - r2) * x / template_blas_sqrt((ergo_real)15))
 
#define solid_harmonic_f_5(x, y, z, x2, y2, z2, r2)   (0.5 * (x2 - y2) * z)
 
#define solid_harmonic_f_6(x, y, z, x2, y2, z2, r2)   (0.5 * template_blas_sqrt((ergo_real)2.5) * (x2 - 3 * y2) * x / template_blas_sqrt((ergo_real)15))
 
#define MAX_DEPTH   888
 
#define MAX_NO_OF_TEST_POINTS   1000
 

Typedefs

typedef struct rhoTreeNode_ rhoTreeNode
 
typedef std::vector< std::vector< std::vector< int > > > tripleVectorOfInt
 
typedef real coor3DPtr[3]
 

Functions

static void print_box (BoxStruct *box)
 
static void get_distribution_box (BoxStruct *box, DistributionSpecStruct *distr, real targetRhoError)
 
static void get_shell_box (BoxStruct *box, ShellSpecStructWithExtent *shell)
 
static real compute_value_at_point (DensitySpecStruct *density, int noOfNonzeroShells, int *nonZeroShellsIndexList, int noOfNonzeroBasFuncs, int *nonZeroBasFuncsIndexList, const real *localFullDensityMatrix, real(*coor)[3], real *workList)
 
static real compute_integral_from_points (const BasisInfoStruct &bis, DensitySpecStruct *density, int noOfNonzeroShells, int *nonZeroShellsIndexList, int noOfNonzeroBasFuncs, int *nonZeroBasFuncsIndexList, const real *localFullDensityMatrix, int nPoints, real(*coor)[3], real *weight, real *workList, real &minValue, real &maxValue, real &maxAbsValue)
 
template<class Treal >
Treal hicuErf (Treal a)
 
template<>
float hicuErf (float a)
 
template<>
double hicuErf (double a)
 
static real to_power (real x, int n)
 
static real compute_1d_gaussian_integral_recursive (real a, real b, int n, real alpha)
 
static real compute_1d_gaussian_integral (real a, real b, int n, real alpha)
 
static real compute_integral_over_box (DistributionSpecStruct *distr, BoxStruct *box)
 
static int get_rhotree_indexes_for_box (int *resultList, int resultListMaxCount, const rhoTreeNode *node, const BoxStruct *inputBoxPtr)
 
static void callbackGga (DftIntegratorBl *grid, int bllen, real &energy)
 
static void callbackLda (DftIntegratorBl *grid, int bllen, real &energy)
 
static void integrate_density_and_energy (const BasisInfoStruct &bis, DensitySpecStruct *density, DftIntegratorBl *integrator, real &electrons, real &energy, int noOfGridPoints, real(*coor)[3], real *weight, real *dmagao)
 
static int compute_grid_for_box (compute_grid_for_box_params_struct *params, int maxlen, real(*coor)[3], real *weight, BoxStruct *box, real analyticalIntegralValue, real *workList, ComputeGridResultValuesStruct &resultValues, bool resolutionIsOk)
 
static rhoTreeNodeBuildRhoTreeBranch (int noOfDistributionsTot, DistributionSpecStruct *rho_alt_1, ShellSpecStructWithExtent *rho_alt_2, int distrIndexListN, int *distrIndexList, real targetRhoError)
 
static rhoTreeNodeBuildRhoTree (int noOfDistributions, DistributionSpecStruct *rho_alt_1, ShellSpecStructWithExtent *rho_alt_2, real targetRhoError)
 
static void free_rho_tree_memory (rhoTreeNode *rootNode)
 
static int round_real (real x)
 
static void getSubBox (const BoxStruct &startBox, BoxStruct &subBox, int Nx, int Ny, int Nz, int ix, int iy, int iz)
 
static void * compute_grid_thread_func (void *arg)
 
static int compute_grid (const BasisInfoStruct &bis, DensitySpecStruct *density, const GridGenerationParamsStruct &gridGenerationParams, real boxdist, real startBoxSizeDebug, const char *gridFileName, int noOfThreads, bool generateSparsePatternOnly, Dft::SparsePattern *sparsePattern)
 
static int do_merge_sort_distrs (int n, DistributionSpecStruct *list, DistributionSpecStruct *workList)
 
static int compute_extent_for_shells (ShellSpecStructWithExtent *shellList, const BasisInfoStruct &bis, real targetRhoError)
 
static int get_product_distrs (const BasisInfoStruct &bis, const Dft::Matrix &dmat, real targetRhoError, DistributionSpecStruct *rho, int maxCount)
 
static void get_shell_list_with_extents (const BasisInfoStruct &bis, int maxCountShellList, ShellSpecStructWithExtent *shellList, real targetRhoError)
 
static int get_density (const BasisInfoStruct &bis, DistributionSpecStruct *rho, int maxCountRho, real targetRhoError, int nbast, const Dft::Matrix &dmat, BasisFuncStruct *basisFuncList)
 
int hicu_grid_generate (const char *grid_file_name, const BasisInfoStruct &bis, ergo_real maxError, ergo_real boxSize, ergo_real startBoxSizeDebug, int use_error_per_volume, int do_double_checking, int compare_to_refined, int use_energy_criterion, int use_energy_criterion_only, int do_variation_checking, const Dft::Matrix *dmat, Dft::SparsePattern *sparsePattern, int nThreads, bool generateSparsePatternOnly)
 
void grid_generate_sparse_pattern (const BasisInfoStruct &bis, ergo_real maxError, ergo_real boxSize, ergo_real startBoxSizeDebug, Dft::SparsePattern &sparsePattern)
 

Variables

static const int CUBATURE_RULE = 3
 
static const int CUBATURE_RULE_2 = 6
 
const real COORD_DIFF_FOR_SAMEPOINT_CRITERION = 1.0e-11
 
const real DISTR_PRODUCT_THRESHOLD = 1e-12
 
const real RELATIVE_DENSITY_VARIATION_LIMIT = 0.5
 
const real DISTR_COEFF_CUTOFF_VALUE = 1e-12
 
const real TARGET_RHO_ERROR_FACTOR = 1e-4
 
pthread_mutex_t global_main_hicu_mutex = PTHREAD_MUTEX_INITIALIZER
 
const int HICU_GRID_PLOT_RESOLUTION = 50
 

Detailed Description

Hierarchical Cubature (HiCu) grid generation.

Author
: Elias Rudberg responsible

Macro Definition Documentation

◆ __CVERSION__

#define __CVERSION__

◆ DISTR_CENTER_DIST_LIMIT

#define DISTR_CENTER_DIST_LIMIT   1e-22

◆ DO_EXTRA_ERROR_CHECKING

#define DO_EXTRA_ERROR_CHECKING

◆ EXPONENT_DIFF_LIMIT

#define EXPONENT_DIFF_LIMIT   1e-22

◆ FILE_BATCH_N

#define FILE_BATCH_N   1000000

◆ HICU_SPARSE_MATRIX_ACCESS_ROUTINE

#define HICU_SPARSE_MATRIX_ACCESS_ROUTINE   at

◆ MAX_DEPTH

#define MAX_DEPTH   888

◆ MAX_NO_OF_POINTS_PER_BATCH

#define MAX_NO_OF_POINTS_PER_BATCH   100

◆ MAX_NO_OF_POINTS_PER_WRITE

#define MAX_NO_OF_POINTS_PER_WRITE   50000

◆ MAX_NO_OF_SHLBLOCKS

#define MAX_NO_OF_SHLBLOCKS   44444

◆ MAX_NO_OF_TEST_POINTS

#define MAX_NO_OF_TEST_POINTS   1000

◆ N_BATCH_JOBS

#define N_BATCH_JOBS   22

◆ solid_harmonic_d_0

#define solid_harmonic_d_0 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (x * y)

◆ solid_harmonic_d_1

#define solid_harmonic_d_1 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (y * z)

◆ solid_harmonic_d_2

#define solid_harmonic_d_2 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    ((2 * z2 - x2 - y2) / (2 * template_blas_sqrt((ergo_real)3)))

◆ solid_harmonic_d_3

#define solid_harmonic_d_3 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (x * z)

◆ solid_harmonic_d_4

#define solid_harmonic_d_4 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * (x2 - y2))

◆ solid_harmonic_f_0

#define solid_harmonic_f_0 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    ((0.5 * template_blas_sqrt(2.5) * (3 * x2 - y2) * y) / template_blas_sqrt((ergo_real)15))

◆ solid_harmonic_f_1

#define solid_harmonic_f_1 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (x * y * z)

◆ solid_harmonic_f_2

#define solid_harmonic_f_2 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * template_blas_sqrt((ergo_real)1.5) * (5 * z2 - r2) * y / template_blas_sqrt((ergo_real)15))

◆ solid_harmonic_f_3

#define solid_harmonic_f_3 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * (5 * z2 - 3 * r2) * z / template_blas_sqrt((ergo_real)15))

◆ solid_harmonic_f_4

#define solid_harmonic_f_4 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * template_blas_sqrt((ergo_real)1.5) * (5 * z2 - r2) * x / template_blas_sqrt((ergo_real)15))

◆ solid_harmonic_f_5

#define solid_harmonic_f_5 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * (x2 - y2) * z)

◆ solid_harmonic_f_6

#define solid_harmonic_f_6 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    (0.5 * template_blas_sqrt((ergo_real)2.5) * (x2 - 3 * y2) * x / template_blas_sqrt((ergo_real)15))

◆ solid_harmonic_p_0

#define solid_harmonic_p_0 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    y

◆ solid_harmonic_p_1

#define solid_harmonic_p_1 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    z

◆ solid_harmonic_p_2

#define solid_harmonic_p_2 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    x

◆ solid_harmonic_s_0

#define solid_harmonic_s_0 (   x,
  y,
  z,
  x2,
  y2,
  z2,
  r2 
)    1

◆ USE_ERF_STD

#define USE_ERF_STD

◆ USE_EXP_STD

#define USE_EXP_STD

Typedef Documentation

◆ coor3DPtr

typedef real coor3DPtr[3]

◆ rhoTreeNode

typedef struct rhoTreeNode_ rhoTreeNode

◆ tripleVectorOfInt

typedef std::vector< std::vector< std::vector<int> > > tripleVectorOfInt

Function Documentation

◆ BuildRhoTree()

static rhoTreeNode* BuildRhoTree ( int  noOfDistributions,
DistributionSpecStruct rho_alt_1,
ShellSpecStructWithExtent rho_alt_2,
real  targetRhoError 
)
static

◆ BuildRhoTreeBranch()

◆ callbackGga()

◆ callbackLda()

◆ compute_1d_gaussian_integral()

static real compute_1d_gaussian_integral ( real  a,
real  b,
int  n,
real  alpha 
)
static

◆ compute_1d_gaussian_integral_recursive()

static real compute_1d_gaussian_integral_recursive ( real  a,
real  b,
int  n,
real  alpha 
)
static

◆ compute_extent_for_shells()

◆ compute_grid()

◆ compute_grid_for_box()

static int compute_grid_for_box ( compute_grid_for_box_params_struct params,
int  maxlen,
real(*)  coor[3],
real weight,
BoxStruct box,
real  analyticalIntegralValue,
real workList,
ComputeGridResultValuesStruct resultValues,
bool  resolutionIsOk 
)
static

References compute_grid_for_box_params_struct::bis, GridGenerationParamsStruct::compareToRefined, compute_integral_from_points(), compute_integral_over_box(), CUBATURE_RULE, CUBATURE_RULE_2, compute_grid_for_box_params_struct::density, compute_grid_for_box_params_struct::dftIntegrator, DensitySpecStruct::distrList, compute_grid_for_box_params_struct::dmagao, GridGenerationParamsStruct::doDoubleChecking, GridGenerationParamsStruct::doVariationChecking, ComputeGridResultValuesStruct::estimatedIntegralErrorDensity, ComputeGridResultValuesStruct::estimatedIntegralErrorEnergy, compute_grid_for_box_params_struct::gridGenerationParams, integrate_density_and_energy(), compute_grid_for_box_params_struct::localFullDensityMatrix, BoxStruct_::max, MAX_NO_OF_TEST_POINTS, GridGenerationParamsStruct::maxerrorPerBox, BoxStruct_::min, NO_OF_DIMENSIONS, compute_grid_for_box_params_struct::nonZeroBasisFuncIndexList, compute_grid_for_box_params_struct::nonZeroShellsIndexList, DensitySpecStruct::noOfDistributions, compute_grid_for_box_params_struct::noOfNonzeroBasisFuncs, compute_grid_for_box_params_struct::noOfNonzeroShells, RELATIVE_DENSITY_VARIATION_LIMIT, GridGenerationParamsStruct::targetRhoError, template_blas_fabs(), ComputeGridResultValuesStruct::totalIntegralResultAnalytical, ComputeGridResultValuesStruct::totalIntegralResultEnergy, ComputeGridResultValuesStruct::totalIntegralResultNumerical, use_cubature_rule(), GridGenerationParamsStruct::useEnergyCriterion, GridGenerationParamsStruct::useEnergyCriterionOnly, and GridGenerationParamsStruct::useErrorPerVolume.

Referenced by compute_grid_thread_func().

◆ compute_grid_thread_func()

static void* compute_grid_thread_func ( void *  arg)
static

References Dft::Matrix::at(), DftIntegratorBl_::bas_bl_cnt, DftIntegratorBl_::basblocks, compute_grid_thread_func_struct::bis, compute_grid_for_box(), compute_integral_over_box(), compute_grid_thread_func_struct::counterArrForPlot, compute_grid_thread_func_struct::currJobNumber, compute_grid_for_box_params_struct::density, compute_grid_thread_func_struct::density, dft_integrator_bl_free(), dft_integrator_bl_new(), DFT_MAX_BLLEN, compute_grid_for_box_params_struct::dftIntegrator, DensitySpecStruct::distrList, compute_grid_for_box_params_struct::dmagao, DensitySpecStruct::dmat, do_output(), ergoShellsToOrbs(), FILE_BATCH_N, compute_grid_thread_func_struct::fileMutex, compute_grid_thread_func_struct::generateSparsePatternOnly, get_rhotree_indexes_for_box(), getSubBox(), compute_grid_thread_func_struct::gridFile, compute_grid_for_box_params_struct::gridGenerationParams, compute_grid_thread_func_struct::gridGenerationParams, HICU_GRID_PLOT_RESOLUTION, compute_grid_thread_func_struct::jobMutex, compute_grid_for_box_params_struct::listShlblocks_otherformat, compute_grid_for_box_params_struct::localFullDensityMatrix, LOG_AREA_DFT, LOG_CAT_ERROR, LOG_CAT_INFO, BoxStruct_::max, MAX_NO_OF_POINTS_PER_BATCH, MAX_NO_OF_POINTS_PER_WRITE, MAX_NO_OF_SHLBLOCKS, compute_grid_thread_func_struct::maxNoOfRelevantDistrsPerBox, BoxStruct_::min, N_BATCH_JOBS, DensitySpecStruct::nbast, compute_grid_for_box_params_struct::nonZeroBasisFuncIndexList, compute_grid_for_box_params_struct::nonZeroShellsIndexList, DensitySpecStruct::noOfDistributions, compute_grid_for_box_params_struct::noOfNonzeroBasisFuncs, compute_grid_for_box_params_struct::noOfNonzeroShells, compute_grid_thread_func_struct::noOfPoints, DensitySpecStruct::noOfShells, compute_grid_thread_func_struct::noOfWrittenBatches, compute_grid_for_box_params_struct::nShlblocks, compute_grid_thread_func_struct::Nx, compute_grid_thread_func_struct::Ny, compute_grid_thread_func_struct::Nz, compute_grid_thread_func_struct::resultCode, compute_grid_thread_func_struct::resultValues, compute_grid_thread_func_struct::rhoTreeRootNode, compute_grid_thread_func_struct::rhoTreeRootNodeShells, ShellSpecStructWithExtent::s, selected_func, DensitySpecStruct::shellList, ShellSpecStruct::shellType, DftIntegratorBl_::shl_bl_cnt, DftIntegratorBl_::shlblocks, compute_grid_thread_func_struct::sparsePattern, compute_grid_thread_func_struct::startBox, ShellSpecStruct::startIndexInMatrix, and compute_grid_thread_func_struct::threadNo.

Referenced by compute_grid().

◆ compute_integral_from_points()

static real compute_integral_from_points ( const BasisInfoStruct bis,
DensitySpecStruct density,
int  noOfNonzeroShells,
int *  nonZeroShellsIndexList,
int  noOfNonzeroBasFuncs,
int *  nonZeroBasFuncsIndexList,
const real localFullDensityMatrix,
int  nPoints,
real(*)  coor[3],
real weight,
real workList,
real minValue,
real maxValue,
real maxAbsValue 
)
static

◆ compute_integral_over_box()

◆ compute_value_at_point()

◆ do_merge_sort_distrs()

static int do_merge_sort_distrs ( int  n,
DistributionSpecStruct list,
DistributionSpecStruct workList 
)
static

◆ free_rho_tree_memory()

static void free_rho_tree_memory ( rhoTreeNode rootNode)
static

References rhoTreeNode_::child1, and rhoTreeNode_::child2.

Referenced by compute_grid().

◆ get_density()

◆ get_distribution_box()

◆ get_product_distrs()

◆ get_rhotree_indexes_for_box()

static int get_rhotree_indexes_for_box ( int *  resultList,
int  resultListMaxCount,
const rhoTreeNode node,
const BoxStruct inputBoxPtr 
)
static

◆ get_shell_box()

◆ get_shell_list_with_extents()

static void get_shell_list_with_extents ( const BasisInfoStruct bis,
int  maxCountShellList,
ShellSpecStructWithExtent shellList,
real  targetRhoError 
)
static

◆ getSubBox()

static void getSubBox ( const BoxStruct startBox,
BoxStruct subBox,
int  Nx,
int  Ny,
int  Nz,
int  ix,
int  iy,
int  iz 
)
static

◆ grid_generate_sparse_pattern()

◆ hicu_grid_generate()

◆ hicuErf() [1/3]

template<>
double hicuErf ( double  a)

◆ hicuErf() [2/3]

template<>
float hicuErf ( float  a)

◆ hicuErf() [3/3]

template<class Treal >
Treal hicuErf ( Treal  a)

◆ integrate_density_and_energy()

◆ print_box()

static void print_box ( BoxStruct box)
static

◆ round_real()

static int round_real ( real  x)
static

Referenced by compute_grid().

◆ to_power()

static real to_power ( real  x,
int  n 
)
static

Variable Documentation

◆ COORD_DIFF_FOR_SAMEPOINT_CRITERION

const real COORD_DIFF_FOR_SAMEPOINT_CRITERION = 1.0e-11

Referenced by BuildRhoTreeBranch().

◆ CUBATURE_RULE

const int CUBATURE_RULE = 3
static

Referenced by compute_grid_for_box().

◆ CUBATURE_RULE_2

const int CUBATURE_RULE_2 = 6
static

Referenced by compute_grid_for_box().

◆ DISTR_COEFF_CUTOFF_VALUE

const real DISTR_COEFF_CUTOFF_VALUE = 1e-12

Referenced by get_density(), and get_product_distrs().

◆ DISTR_PRODUCT_THRESHOLD

const real DISTR_PRODUCT_THRESHOLD = 1e-12

Referenced by get_product_distrs().

◆ global_main_hicu_mutex

pthread_mutex_t global_main_hicu_mutex = PTHREAD_MUTEX_INITIALIZER

◆ HICU_GRID_PLOT_RESOLUTION

const int HICU_GRID_PLOT_RESOLUTION = 50

◆ RELATIVE_DENSITY_VARIATION_LIMIT

const real RELATIVE_DENSITY_VARIATION_LIMIT = 0.5

Referenced by compute_grid_for_box().

◆ TARGET_RHO_ERROR_FACTOR

const real TARGET_RHO_ERROR_FACTOR = 1e-4