Skip to content

Commit

Permalink
Modularize Libint2 constructor (#3037)
Browse files Browse the repository at this point in the history 
* Modularize the construction of Libint2 shell data

* n_prim_per_shell_ was never initialized
  • Loading branch information
@susilehtola
susilehtola committed Aug 28, 2023
1 parent b5f5dea commit d58c518
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 25 deletions.
58 changes: 37 additions & 21 deletions psi4/src/psi4/libmints/basisset.cc
View file
Expand Up @@ -118,12 +118,11 @@ BasisSet::BasisSet() {
name_ = "(Empty Basis Set)";
key_ = "(Empty Basis Set)";
target_ = "(Empty Basis Set)";
shells_[0] = GaussianShell(Gaussian, 0, nprimitive_, uoriginal_coefficients_.data(), ucoefficients_.data(), uerd_coefficients_.data(),
uexponents_.data(), GaussianType(0), 0, xyz_.data(), 0);
shells_[0] = GaussianShell(Gaussian, 0, nprimitive_, uoriginal_coefficients_.data(), ucoefficients_.data(),
uerd_coefficients_.data(), uexponents_.data(), GaussianType(0), 0, xyz_.data(), 0);
}

BasisSet::~BasisSet() {
}
BasisSet::~BasisSet() {}

std::shared_ptr<BasisSet> BasisSet::build(std::shared_ptr<Molecule> /*molecule*/,
const std::vector<ShellInfo> & /*shells*/) {
Expand Down Expand Up @@ -213,7 +212,7 @@ int BasisSet::n_frozen_core(const std::string &depth, SharedMolecule mol) {
// If center is a post-lanthanide or a post-actinide in Nth period,
// freeze all 14 of its (N-2)f electrons too
if (current_shell > 5) {
if ((Z + ECP - delta) >= 18 ) delta += 14;
if ((Z + ECP - delta) >= 18) delta += 14;
}
// If this center has an ECP, some electrons are already frozen
if (ECP > 0) delta -= ECP;
Expand Down Expand Up @@ -262,7 +261,9 @@ int BasisSet::n_frozen_core(const std::string &depth, SharedMolecule mol) {
int current_shell = atom_to_period(Z + ECP);
int delta = period_to_full_shell(std::max(current_shell - req_shell, 0));
// If this center has an ECP, some electrons are already frozen
if (delta < ECP) throw PSIEXCEPTION("ECP on atom freezes more electrons than requested by choosing a previous shell.");
if (delta < ECP)
throw PSIEXCEPTION(
"ECP on atom freezes more electrons than requested by choosing a previous shell.");
if (ECP > 0) delta -= ECP;
// Keep track of current valence electrons
mol_valence = mol_valence + Z - delta;
Expand Down Expand Up @@ -766,10 +767,10 @@ BasisSet::BasisSet(const std::string &basistype, SharedMolecule mol,
}

shell_first_ao_ = std::vector<int>(n_shells_);
shell_first_exponent_ = std::vector<int>(n_shells_);
shell_first_basis_function_ = std::vector<int>(n_shells_);
shells_ = std::vector<GaussianShell>(n_shells_);
ecp_shells_ = std::vector<GaussianShell>(n_ecp_shells_);
l2_shells_.resize(n_shells_);
ao_to_shell_ = std::vector<int>(nao_);
function_to_shell_ = std::vector<int>(nbf_);
function_center_ = std::vector<int>(nbf_);
Expand Down Expand Up @@ -823,12 +824,8 @@ BasisSet::BasisSet(const std::string &basistype, SharedMolecule mol,
GaussianShell(shelltype, am, shell_nprim, &uoriginal_coefficients_[ustart + atom_nprim],
&ucoefficients_[ustart + atom_nprim], &uerd_coefficients_[ustart + atom_nprim],
&uexponents_[ustart + atom_nprim], puream, n, &xyz_.data()[3 * n], bf_count);
auto l2c = libint2::svector<double>(&uoriginal_coefficients_[ustart + atom_nprim],
&uoriginal_coefficients_[ustart + atom_nprim + shell_nprim]);
auto l2e = libint2::svector<double>(&uexponents_[ustart + atom_nprim],
&uexponents_[ustart + atom_nprim + shell_nprim]);
l2_shells_[shell_count] =
libint2::Shell{l2e, {{am, static_cast<bool>(puream), l2c}}, {{xyz[0], xyz[1], xyz[2]}}};
shell_first_exponent_[shell_count] = ustart + atom_nprim;
n_prim_per_shell_[shell_count] = shell_nprim;
} else {
throw PSIEXCEPTION("Unexpected shell type in BasisSet constructor!");
}
Expand All @@ -846,6 +843,8 @@ BasisSet::BasisSet(const std::string &basistype, SharedMolecule mol,
throw PSIEXCEPTION("Problem with nprimitive in basis set construction!");
}
}
// Update the libint2 shell data
update_l2_shells();

/*
* Now loop over ECPs and finalize metadata
Expand All @@ -872,9 +871,9 @@ BasisSet::BasisSet(const std::string &basistype, SharedMolecule mol,
max_ecp_am_ = max_ecp_am_ > std::abs(am) ? max_ecp_am_ : std::abs(am);
ecp_shell_center_[ecp_shell_count] = n;
if (shelltype == ECPType1 || shelltype == ECPType2) {
ecp_shells_[ecp_shell_count] =
GaussianShell(shelltype, am, ecp_shell_nprim, &uecpcoefficients_[ustart + atom_nprim],
&uecpexponents_[ustart + atom_nprim], &uecpns_[ustart + atom_nprim], n, &xyz_.data()[3 * n]);
ecp_shells_[ecp_shell_count] = GaussianShell(
shelltype, am, ecp_shell_nprim, &uecpcoefficients_[ustart + atom_nprim],
&uecpexponents_[ustart + atom_nprim], &uecpns_[ustart + atom_nprim], n, &xyz_.data()[3 * n]);
} else {
throw PSIEXCEPTION("Unknown ECP shell type in BasisSet constructor!");
}
Expand All @@ -888,6 +887,21 @@ BasisSet::BasisSet(const std::string &basistype, SharedMolecule mol,
}
}

void BasisSet::update_l2_shells() {
l2_shells_.resize(n_shells_);
for (auto ishell = 0; ishell < n_shells_; ishell++) {
auto am = shells_[ishell].am();
auto center_index = shell_to_center(ishell);
Vector3 xyz = molecule_->xyz(center_index);

auto offset = shell_first_exponent_[ishell];
auto nprim = n_prim_per_shell_[ishell];
auto l2c = libint2::svector<double>(&uoriginal_coefficients_[offset], &uoriginal_coefficients_[offset + nprim]);
auto l2e = libint2::svector<double>(&uexponents_[offset], &uexponents_[offset + nprim]);
l2_shells_[ishell] = libint2::Shell{l2e, {{am, puream_, l2c}}, {{xyz[0], xyz[1], xyz[2]}}};
}
}

std::string BasisSet::make_filename(const std::string &name) {
// Modify the name of the basis set to generate a filename: STO-3G -> sto-3g
std::string basisname = name;
Expand Down Expand Up @@ -1207,7 +1221,7 @@ void BasisSet::move_atom(int atom, const Vector3 &trans) {
}
}

void BasisSet::compute_phi(double* phi_ao, double x, double y, double z) {
void BasisSet::compute_phi(double *phi_ao, double x, double y, double z) {
zero_arr(phi_ao, nbf());

int ao = 0;
Expand All @@ -1233,8 +1247,9 @@ void BasisSet::compute_phi(double* phi_ao, double x, double y, double z) {

for (int l = 0; l < INT_NCART(am); l++) {
Vector3 &components = exp_ao[am][l];
cart_buffer[l] += pow(dx, static_cast<double>(components[0])) * pow(dy, static_cast<double>(components[1])) *
pow(dz, static_cast<double>(components[2])) * cexpr;
cart_buffer[l] += pow(dx, static_cast<double>(components[0])) *
pow(dy, static_cast<double>(components[1])) *
pow(dz, static_cast<double>(components[2])) * cexpr;
}

for (int ind = 0; ind < s_transform.n(); ind++) {
Expand All @@ -1248,8 +1263,9 @@ void BasisSet::compute_phi(double* phi_ao, double x, double y, double z) {
} else {
for (int l = 0; l < INT_NCART(am); l++) {
Vector3 &components = exp_ao[am][l];
phi_ao[ao + l] += pow(dx, static_cast<double>(components[0])) * pow(dy, static_cast<double>(components[1])) *
pow(dz, static_cast<double>(components[2])) * cexpr;
phi_ao[ao + l] += pow(dx, static_cast<double>(components[0])) *
pow(dy, static_cast<double>(components[1])) *
pow(dz, static_cast<double>(components[2])) * cexpr;
}
}

Expand Down
13 changes: 9 additions & 4 deletions psi4/src/psi4/libmints/basisset.h
View file
Expand Up @@ -48,7 +48,7 @@ PRAGMA_WARNING_IGNORE_DEPRECATED_DECLARATIONS
PRAGMA_WARNING_POP

namespace libint2 {
struct Shell;
struct Shell;
}
namespace psi {

Expand Down Expand Up @@ -81,7 +81,7 @@ class PSI_API BasisSet {
std::vector<GaussianShell> shells_;
//! Array of ECP shells
std::vector<GaussianShell> ecp_shells_;
//! Array of Libint2 shells
//! Array of Libint2 shells; updated from shells with update_l2_shells()
std::vector<libint2::Shell> l2_shells_;

//! The number of core electrons for each atom type
Expand Down Expand Up @@ -128,6 +128,8 @@ class PSI_API BasisSet {
std::vector<int> n_prim_per_shell_;
/// The first (Cartesian) atomic orbital in each shell
std::vector<int> shell_first_ao_;
/// First exponent for i:th shell
std::vector<int> shell_first_exponent_;
/// The first (Cartesian / spherical) basis function in each shell
std::vector<int> shell_first_basis_function_;
/// Shell number to atomic center.
Expand Down Expand Up @@ -166,6 +168,9 @@ class PSI_API BasisSet {
/// The flattened list of Cartesian coordinates for each atom
std::vector<double> xyz_;

/// Update Libint2 shells
void update_l2_shells();

public:
BasisSet();

Expand Down Expand Up @@ -400,8 +405,8 @@ class PSI_API BasisSet {
void move_atom(int atom, const Vector3 &trans);
// Returns the values of the basis functions at a point
void compute_phi(double *phi_ao, double x, double y, double z);
private:

private:
/// Helper functions for frozen core to reduce LOC
int atom_to_period(int Z);
int period_to_full_shell(int p);
Expand Down

0 comments on commit d58c518

Please sign in to comment.