Skip to content

Add Hirshfeld partitioning to weights generation #209

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
awvwgk wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add Hirshfeld partitioning to weights generation #209

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions include/gauxc/c/enums.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,10 @@ enum GauXC_AtomicGridSizeDefault {
*/
enum GauXC_XCWeightAlg {
GauXC_XCWeightAlg_NOTPARTITIONED, ///< Not partitioned
GauXC_XCWeightAlg_Becke, ///< The original Becke weighting scheme
GauXC_XCWeightAlg_SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme
GauXC_XCWeightAlg_LKO ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme
GauXC_XCWeightAlg_Becke, ///< The original Becke weighting scheme
GauXC_XCWeightAlg_SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme
GauXC_XCWeightAlg_LKO, ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme
GauXC_XCWeightAlg_Hirshfeld ///< Hirshfeld pro-atom density weighting scheme
};

/**
Expand Down
3 changes: 2 additions & 1 deletion include/gauxc/enums.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,8 @@ enum class XCWeightAlg {
NOTPARTITIONED, ///< Not partitioned
Becke, ///< The original Becke weighting scheme
SSF, ///< The Stratmann-Scuseria-Frisch weighting scheme
LKO ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme
LKO, ///< The Lauqua-Kuessman-Ochsenfeld weighting scheme
Hirshfeld ///< Hirshfeld pro-atom density weighting scheme
};

/**
Expand Down
236 changes: 236 additions & 0 deletions src/xc_integrator/local_work_driver/host/reference/weights.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@

#include <gauxc/molgrid/defaults.hpp>

#include <cmath>

namespace GauXC {

// Reference Becke weights impl
Expand Down Expand Up @@ -385,6 +387,74 @@ void reference_lko_weights_host(

}

void reference_hirshfeld_weights_host(
const Molecule& mol,
const MolMeta& meta,
task_iterator task_begin,
task_iterator task_end
) {

(void)meta;

constexpr double pi = 3.141592653589793238462643383279502884;

const size_t ntasks = std::distance(task_begin,task_end);
const size_t natoms = mol.natoms();

struct hirshfeld_atom_data {
double x;
double y;
double z;
double alpha;
Comment on lines +404 to +408
double scale;
};

std::vector<hirshfeld_atom_data> atom_data;
atom_data.reserve(natoms);

for( const auto& atom : mol ) {
const double radius = default_atomic_radius(atom.Z);
const double alpha = 0.5 / (radius * radius);
const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5);
atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } );
}

#pragma omp parallel
{

#pragma omp for
for( size_t iT = 0; iT < ntasks; ++iT )
for( size_t i = 0; i < (task_begin+iT)->points.size(); ++i ) {

auto& task = *(task_begin+iT);
auto& weight = task.weights[i];
const auto& point = task.points[i];
const auto parent = static_cast<size_t>(task.iParent);

double sum = 0.;
double parent_density = 0.;
for(size_t iA = 0; iA < natoms; iA++) {
const auto& atom = atom_data[iA];

const double da_x = point[0] - atom.x;
const double da_y = point[1] - atom.y;
const double da_z = point[2] - atom.z;
const double r2 = da_x*da_x + da_y*da_y + da_z*da_z;
const double density = atom.scale * std::exp(-atom.alpha * r2);

sum += density;
if( iA == parent ) parent_density = density;

}

weight *= (sum > 0.) ? parent_density / sum : 0.;

} // Collapsed loop over tasks and points

} // OMP context

}


/**
* 1st derivative which expects weight_deri to be preallocated as (ngrid*natoms*3)
Expand Down Expand Up @@ -682,6 +752,81 @@ void reference_ssf_weights_1st_derivative_host(
}
}

void reference_hirshfeld_weights_1st_derivative_host(
const Molecule& mol,
const MolMeta& meta,
const XCTask& task,
double* weight_deri
){

(void)meta;

constexpr double pi = 3.141592653589793238462643383279502884;

const size_t natoms = mol.natoms();

struct hirshfeld_atom_data {
double x;
double y;
double z;
double alpha;
double scale;
};
Comment on lines +768 to +774

std::vector<hirshfeld_atom_data> atom_data;
atom_data.reserve(natoms);
for( const auto& atom : mol ) {
const double radius = default_atomic_radius(atom.Z);
const double alpha = 0.5 / (radius * radius);
const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5);
atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } );
}

std::vector<double> density(natoms);

for( size_t i = 0; i < task.points.size(); ++i ) {
auto* weight_deri_ith = weight_deri + 3*natoms*i;
std::fill(weight_deri_ith, weight_deri_ith + 3*natoms, 0.);

const auto& point = task.points[i];
const double weight = task.weights[i];
if( std::abs(weight) < 1.e-14 ) continue;

const size_t iParent = task.iParent;

double sum = 0.;
for( size_t iA = 0; iA < natoms; ++iA ) {
const auto& atom = atom_data[iA];
const double da_x = point[0] - atom.x;
const double da_y = point[1] - atom.y;
const double da_z = point[2] - atom.z;
const double r2 = da_x*da_x + da_y*da_y + da_z*da_z;
density[iA] = atom.scale * std::exp(-atom.alpha * r2);
sum += density[iA];
}

if( sum <= 0. ) continue;

auto* weight_deri_iParent = weight_deri_ith + 3*iParent;
for( size_t iB = 0; iB < natoms; ++iB ) {
if( iB == iParent ) continue;

const auto& atom = atom_data[iB];
const double scaled_density = density[iB] / sum;
const double coef = -2. * atom.alpha * scaled_density * weight;

auto* weight_deri_iB = weight_deri_ith + 3*iB;
weight_deri_iB[0] = coef * (point[0] - atom.x);
weight_deri_iB[1] = coef * (point[1] - atom.y);
weight_deri_iB[2] = coef * (point[2] - atom.z);

weight_deri_iParent[0] -= weight_deri_iB[0];
weight_deri_iParent[1] -= weight_deri_iB[1];
weight_deri_iParent[2] -= weight_deri_iB[2];
}
}
}



/**
Expand Down Expand Up @@ -983,6 +1128,97 @@ void reference_ssf_weights_1std_contraction_host(

}

void reference_hirshfeld_weights_1std_contraction_host(
const Molecule& mol,
const MolMeta& meta,
const XCTask& task,
const double* w_times_f,
double* exc_grad_w
){

(void)meta;

constexpr double pi = 3.141592653589793238462643383279502884;
constexpr double w_times_f_thresh = 1.e-14;

const size_t natoms = mol.natoms();

struct hirshfeld_atom_data {
double x;
double y;
double z;
double alpha;
double scale;
};
Comment on lines +1146 to +1152

std::vector<hirshfeld_atom_data> atom_data;
atom_data.reserve(natoms);
for( const auto& atom : mol ) {
const double radius = default_atomic_radius(atom.Z);
const double alpha = 0.5 / (radius * radius);
const double scale = atom.Z.get() * std::pow(alpha / pi, 1.5);
atom_data.push_back( hirshfeld_atom_data{ atom.x, atom.y, atom.z, alpha, scale } );
}

std::vector<double> density(natoms);

for( size_t i = 0; i < task.points.size(); ++i ) {
const double w_times_f_i = w_times_f[i];
if( std::abs(w_times_f_i) < w_times_f_thresh ) continue;

const auto& point = task.points[i];
const size_t iParent = task.iParent;

double sum = 0.;
for( size_t iA = 0; iA < natoms; ++iA ) {
const auto& atom = atom_data[iA];
const double da_x = point[0] - atom.x;
const double da_y = point[1] - atom.y;
const double da_z = point[2] - atom.z;
const double r2 = da_x*da_x + da_y*da_y + da_z*da_z;
density[iA] = atom.scale * std::exp(-atom.alpha * r2);
sum += density[iA];
}

if( sum <= 0. ) continue;

double parent_x = 0.;
double parent_y = 0.;
double parent_z = 0.;

for( size_t iB = 0; iB < natoms; ++iB ) {
if( iB == iParent ) continue;

const auto& atom = atom_data[iB];
const double scaled_density = density[iB] / sum;
const double coef = -2. * atom.alpha * scaled_density * w_times_f_i;

const double grad_x = coef * (point[0] - atom.x);
const double grad_y = coef * (point[1] - atom.y);
const double grad_z = coef * (point[2] - atom.z);

#pragma omp atomic
exc_grad_w[3*iB + 0] += grad_x;
#pragma omp atomic
exc_grad_w[3*iB + 1] += grad_y;
#pragma omp atomic
exc_grad_w[3*iB + 2] += grad_z;

parent_x -= grad_x;
parent_y -= grad_y;
parent_z -= grad_z;
}

#pragma omp atomic
exc_grad_w[3*iParent + 0] += parent_x;
#pragma omp atomic
exc_grad_w[3*iParent + 1] += parent_y;
#pragma omp atomic
exc_grad_w[3*iParent + 2] += parent_z;
}

}



}
22 changes: 22 additions & 0 deletions src/xc_integrator/local_work_driver/host/reference/weights.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,13 @@ void reference_lko_weights_host(
task_iterator task_end
);

void reference_hirshfeld_weights_host(
const Molecule& mol,
const MolMeta& meta,
task_iterator task_begin,
task_iterator task_end
);

void reference_becke_weights_1st_derivative_host(
const Molecule& mol,
const MolMeta& meta,
Expand All @@ -50,6 +57,13 @@ void reference_ssf_weights_1st_derivative_host(
double* weight_deri
);

void reference_hirshfeld_weights_1st_derivative_host(
const Molecule& mol,
const MolMeta& meta,
const XCTask& task,
double* weight_deri
);

// Becke weights 1st derivative contracted with integrator
void reference_becke_weights_1std_contraction_host(
const Molecule& mol,
Expand All @@ -68,4 +82,12 @@ void reference_ssf_weights_1std_contraction_host(
double* exc_grad_w
);

void reference_hirshfeld_weights_1std_contraction_host(
const Molecule& mol,
const MolMeta& meta,
const XCTask& task,
const double* w_times_f,
double* exc_grad_w
);

}
6 changes: 6 additions & 0 deletions src/xc_integrator/local_work_driver/host/reference_local_host_work_driver.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ namespace GauXC {
case XCWeightAlg::LKO:
reference_lko_weights_host( mol, meta, task_begin, task_end );
break;
case XCWeightAlg::Hirshfeld:
reference_hirshfeld_weights_host( mol, meta, task_begin, task_end );
break;
default:
GAUXC_GENERIC_EXCEPTION("Weight Alg Not Supported");
}
Expand All @@ -66,6 +69,9 @@ namespace GauXC {
case XCWeightAlg::SSF:
reference_ssf_weights_1std_contraction_host( mol, meta, task, w_times_f, exc_grad_w );
break;
case XCWeightAlg::Hirshfeld:
reference_hirshfeld_weights_1std_contraction_host( mol, meta, task, w_times_f, exc_grad_w );
break;
default:
GAUXC_GENERIC_EXCEPTION("Weight Alg Not Supported");
}
Expand Down
Loading
Loading