Permalink
Browse files

Merge pull request #1138 from PeterKraus/cubeprop_improvements_01

Cubeprop improvements
  • Loading branch information...
dgasmith committed Sep 6, 2018
2 parents ba459ab + 66711d7 commit 2f23252e612f644c28d702347154282b6117fdf8
@@ -394,3 +394,11 @@ Bibliography
.. [Patkowski:2018:164110]
Konrad Patkowski, Piotr S. \.Zuchowski, Daniel G. A. Smith
*J. Chem. Phys.* **148**, 164110 (2018).
.. [Morell:2005:205]
Christophe Morell, Andr\ |e_acute|\ Grand, and Alejandro Toro-Labb\ |e_acute|
*J. Phys. Chem. A* **109**, 205 (2005).
.. [Martinez-Araya:2015:451]
Jorge Ignacio Mart\ |i_acute|\ nez-Araya
*J. Math. Chem.* **53**, 451 (2015).
@@ -38,8 +38,8 @@
Generation of Cube Files |w---w| :py:func:`~psi4.cubeprop`
==========================================================
.. codeauthor:: Robert M. Parrish and Francesco A. Evangelista
.. sectionauthor:: Francesco A. Evangelista
.. codeauthor:: Robert M. Parrish, Francesco A. Evangelista and Peter Kraus
.. sectionauthor:: Francesco A. Evangelista and Peter Kraus
.. autofunction:: psi4.cubeprop(wfn)
@@ -109,6 +109,13 @@ ORBITALS [Default if |globals__cubeprop_tasks| is not specified]
Produces cube representations of the molecular orbitals
:math:`\psi_q(\mathbf{r})`. Orbitals are sorted according to increasing
orbital energy ignoring symmetry.
FRONTIER_ORBITALS
Produces cube representations of the frontier molecular orbitals. For closed shell
species, the highest occupied (HOMO) and the lowest unoccupied (LUMO) alpha orbitals (ie.
:math:`\psi_{\alpha}(\mathbf{r})`) are printed, while for open shell species a total
of :math:`(4 + M_s)` orbitals are printed (:math:`\alpha` and :math:`\beta`
spin for both lowest virtual (LVMO) and highest doubly occupied
orbitals (DOMO), along with all :math:`\alpha` singly occupied (SOMO) orbitals).
DENSITY
This task can be used to obtain the alpha and beta electron densities,
:math:`\rho_\alpha(\mathbf{r})` and :math:`\rho_\beta(\mathbf{r})`, together
@@ -122,6 +129,14 @@ BASIS_FUNCTIONS
ESP
Calculates the total (nuclear + electronic) electrostatic potential
:math:`V(\mathbf{r})`.
DUAL_DESCRIPTOR
Calculates the dual descriptor from frontier orbitals:
:math:`f^2(\mathbf{r})=\rho_{\mathrm{LUMO}}(\mathbf{r})-\rho_{\mathrm{HOMO}}(\mathbf{r})`.
The dual descriptor is a good measure of nucleophilicity and electrophilicity,
containing information essentially equivalent to both Fukui functions combined.
More details on the dual descriptor itself can be found in [Morell:2005:205]_,
while the current implementation is described in [Martinez-Araya:2015:451]_.
This feature is currently only supported for closed shell systems.
.. note:: The ``ESP`` task requires the user to specify a density-fitting basis
via the |scf__df_basis_scf| keyword.
@@ -598,7 +598,7 @@ void CubicScalarGrid::compute_density(std::shared_ptr<Matrix> D, const std::stri
double density_percent = 100.0 * options_.get_double("CUBEPROP_ISOCONTOUR_THRESHOLD");
std::stringstream comment;
comment << " [e/a0^3]. Isocontour range for " << density_percent << "% of the density: (" << isocontour_range.first
<< "," << isocontour_range.second << ")";
<< "," << isocontour_range.second << ")." << ecp_header();
// Write to disk
write_gen_file(v, name, type, comment.str());
delete[] v;
@@ -649,6 +649,34 @@ void CubicScalarGrid::compute_orbitals(std::shared_ptr<Matrix> C, const std::vec
}
free_block(v);
}
void CubicScalarGrid::compute_difference(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::string& label, const bool square, const std::string& type) {
auto C2 = std::make_shared<Matrix>(primary_->nbf(), indices.size());
double** Cp = C->pointer();
double** C2p = C2->pointer();
for (int k = 0; k < indices.size(); k++) {
C_DCOPY(primary_->nbf(), &Cp[0][indices[k]], C->colspi()[0], &C2p[0][k], C2->colspi()[0]);
}
auto v_t = std::make_shared<Matrix>(indices.size(), npoints_);
double** v_tp = v_t->pointer();
auto v = std::make_shared<Vector>(npoints_);
double* vp = v->pointer();
add_orbitals(&v_tp[0], C2);
for (int i = 0; i < npoints_; i++) {
if (square) {
v->set(0, i, (v_t->get(0,i) - v_t->get(1,i))*(v_t->get(0,i) + v_t->get(1,i)));
} else {
v->set(0, i, (v_t->get(0,i) - v_t->get(1,i)));
}
}
std::pair<double, double> isocontour_range = compute_isocontour_range(&vp[0], 2.0);
double density_percent = 100.0 * options_.get_double("CUBEPROP_ISOCONTOUR_THRESHOLD");
std::stringstream comment;
comment << ". Isocontour range for " << density_percent << "% of the density: (" << isocontour_range.first
<< "," << isocontour_range.second << ")";
// Write to disk
write_gen_file(&vp[0], label, type, comment.str());
}
void CubicScalarGrid::compute_LOL(std::shared_ptr<Matrix> D, const std::string& name, const std::string& type) {
double* v = new double[npoints_];
memset(v, '\0', npoints_ * sizeof(double));
@@ -695,4 +723,22 @@ std::pair<double, double> CubicScalarGrid::compute_isocontour_range(double* v2,
}
return std::make_pair(positive_isocontour, negative_isocontour);
}
std::string CubicScalarGrid::ecp_header() {
std::stringstream ecp_head;
ecp_head.str("");
if (primary_->has_ECP()) {
ecp_head << " Total core: " << primary_->n_ecp_core() << " [e] from 1-indexed atoms (";
std::stringstream ecp_atoms;
std::stringstream ecp_ncore;
for (int A = 0; A < mol_->natom(); A++) {
if (primary_->n_ecp_core(mol_->label(A)) > 0) {
ecp_atoms << A+1 << "[" << mol_->symbol(A) << "], ";
ecp_ncore << primary_->n_ecp_core(mol_->label(A)) << ", ";
}
}
ecp_head << ecp_atoms.str().substr(0,ecp_atoms.str().length()-2) << ") electrons ("
<< ecp_ncore.str().substr(0,ecp_ncore.str().length()-2) << ").";
}
return ecp_head.str();
}
}
@@ -181,6 +181,10 @@ class CubicScalarGrid {
void compute_orbitals(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::vector<std::string>& labels, const std::string& name,
const std::string& type = "CUBE");
/// Compute a set of orbital-type properties and drop files corresponding to name, index, symmetry label, and type
void compute_difference(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::string& label, bool square = false, const std::string& type = "CUBE");
/// Compute a LOL-type property and drop a file corresponding to name and type
void compute_LOL(std::shared_ptr<Matrix> D, const std::string& name, const std::string& type = "CUBE");
/// Compute an ELF-type property and drop a file corresponding to name and type (TODO: this seems very unstable)
@@ -189,6 +193,9 @@ class CubicScalarGrid {
/// Compute the isocountour range that capture a given fraction of a property. Exponent is used
/// to properly compute the density. E.g. for orbitals exponent = 2, for densities exponent = 1
std::pair<double, double> compute_isocontour_range(double* v2, double exponent);
/// A helper function to construct ECP comment in the cubefile header.
std::string ecp_header();
};
} // End namespace
@@ -70,14 +70,18 @@ CubeProperties::CubeProperties(SharedWavefunction wfn) : options_(Process::envir
int nirrep = wfn->nirrep();
Dimension nmopi = wfn->nmopi();
nalpha_ = 0;
// Gather orbital information
for (int h = 0; h < nirrep; h++) {
nalpha_ += wfn->nalphapi()[h];
for (int i = 0; i < (int)nmopi[h]; i++) {
info_a_.push_back(std::tuple<double, int, int>(wfn->epsilon_a()->get(h, i), i, h));
}
}
std::sort(info_a_.begin(), info_a_.end(), std::less<std::tuple<double, int, int> >()); // Sort as in wfn
nbeta_ = 0;
for (int h = 0; h < nirrep; h++) {
nbeta_ += wfn->nbetapi()[h];
for (int i = 0; i < (int)nmopi[h]; i++) {
info_b_.push_back(std::tuple<double, int, int>(wfn->epsilon_b()->get(h, i), i, h));
}
@@ -145,15 +149,68 @@ void CubeProperties::raw_compute_properties() {
for (size_t ind = 0; ind < indsa0.size(); ++ind) {
int i = std::get<1>(info_a_[indsa0[ind]]);
int h = std::get<2>(info_a_[indsa0[ind]]);
labelsa.push_back(to_string(i + 1) + "-" + ct.gamma(h).symbol());
labelsa.push_back(std::to_string(i + 1) + "-" + ct.gamma(h).symbol());
}
for (size_t ind = 0; ind < indsb0.size(); ++ind) {
int i = std::get<1>(info_b_[indsb0[ind]]);
int h = std::get<2>(info_b_[indsb0[ind]]);
labelsb.push_back(to_string(i + 1) + "-" + ct.gamma(h).symbol());
labelsb.push_back(std::to_string(i + 1) + "-" + ct.gamma(h).symbol());
}
if (indsa0.size()) compute_orbitals(Ca_, indsa0, labelsa, "Psi_a");
if (indsb0.size()) compute_orbitals(Cb_, indsb0, labelsb, "Psi_b");
} else if (task == "FRONTIER_ORBITALS") {
std::vector<int> indsa0;
std::vector<int> indsb0;
std::vector<std::string> labelsa;
std::vector<std::string> labelsb;
CharacterTable ct = basisset_->molecule()->point_group()->char_table();
if (nalpha_ == nbeta_) {
indsa0.push_back(nalpha_-1);
labelsa.push_back(std::to_string(std::get<1>(info_a_[nalpha_-1]) + 1) + "-" +
ct.gamma(std::get<2>(info_a_[nalpha_-1])).symbol() + "_HOMO");
indsa0.push_back(nalpha_);
labelsa.push_back(std::to_string(std::get<1>(info_a_[nalpha_]) + 1) + "-" +
ct.gamma(std::get<2>(info_a_[nalpha_])).symbol() + "_LUMO");
} else {
int orb_index = nalpha_;
indsa0.push_back(orb_index);
labelsa.push_back(std::to_string(std::get<1>(info_a_[orb_index]) + 1) + "-" +
ct.gamma(std::get<2>(info_a_[orb_index])).symbol() + "_LVMO");
indsb0.push_back(orb_index);
labelsb.push_back(std::to_string(std::get<1>(info_b_[orb_index]) + 1) + "-" +
ct.gamma(std::get<2>(info_b_[orb_index])).symbol() + "_LVMO");
for (int i = 1; i <= nalpha_ - nbeta_; i++) {
orb_index = nalpha_ - i;
indsa0.push_back(orb_index);
labelsa.push_back(std::to_string(std::get<1>(info_a_[orb_index]) + 1) + "-" +
ct.gamma(std::get<2>(info_a_[orb_index])).symbol() + "_SOMO");
}
orb_index = nbeta_ - 1;
indsa0.push_back(orb_index);
labelsa.push_back(std::to_string(std::get<1>(info_a_[orb_index]) + 1) + "-" +
ct.gamma(std::get<2>(info_a_[orb_index])).symbol() + "_DOMO");
orb_index = nbeta_ - 1;
indsb0.push_back(orb_index);
labelsb.push_back(std::to_string(std::get<1>(info_b_[orb_index]) + 1) + "-" +
ct.gamma(std::get<2>(info_b_[orb_index])).symbol() + "_DOMO");
}
if (indsa0.size()) compute_orbitals(Ca_, indsa0, labelsa, "Psi_a");
if (indsb0.size()) compute_orbitals(Cb_, indsb0, labelsb, "Psi_b");
} else if (task == "DUAL_DESCRIPTOR") {
// Calculates the dual descriptor from frontier molecular orbitals.
// The dual descriptor is a good measure of electro-/nucleophilicity:
// f^2(r) = rho_lumo(r) - rho_homo(r)
// See 10.1021/jp046577a and 10.1007/s10910-014-0437-7
std::vector<int> indsa0;
if (nalpha_ != nbeta_) {
throw PSIEXCEPTION(task + "is not implemented for open-shell systems");
} else {
indsa0.push_back(nalpha_);
indsa0.push_back(nalpha_-1);
std::stringstream ss;
ss << "DUAL_" << (nalpha_+1) << "_LUMO-" << nalpha_ << "_HOMO";
compute_difference(Ca_, indsa0, ss.str(), true);
}
} else if (task == "BASIS_FUNCTIONS") {
std::vector<int> inds0;
if (options_["CUBEPROP_BASIS_FUNCTIONS"].size() == 0) {
@@ -188,6 +245,10 @@ void CubeProperties::compute_orbitals(std::shared_ptr<Matrix> C, const std::vect
const std::vector<std::string>& labels, const std::string& key) {
grid_->compute_orbitals(C, indices, labels, key);
}
void CubeProperties::compute_difference(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::string& label, bool square) {
grid_->compute_difference(C, indices, label, square);
}
void CubeProperties::compute_basis_functions(const std::vector<int>& indices, const std::string& key) {
grid_->compute_basis_functions(indices, key);
}
@@ -64,6 +64,9 @@ class PSI_API CubeProperties {
std::vector<std::tuple<double, int, int> > info_b_;
/// Auxiliary Basis Set if Any
std::shared_ptr<BasisSet> auxiliary_;
/// Total number of alpha/beta electrons
int nalpha_;
int nbeta_;
// => Computers <= //
@@ -100,6 +103,9 @@ class PSI_API CubeProperties {
/// Compute an orbital task (key_N.cube, for 0-based indices of C)
void compute_orbitals(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::vector<std::string>& labels, const std::string& key);
/// Compute a difference task between two indices of matrix C
void compute_difference(std::shared_ptr<Matrix> C, const std::vector<int>& indices,
const std::string& label, bool square);
/// Compute a basis function task (key_N.cube, for 0-based indices of basisset_)
void compute_basis_functions(const std::vector<int>& indices, const std::string& key);
/// Compute a LOL grid task (key.cube)
View
@@ -207,12 +207,14 @@ int read_options(const std::string &name, Options & options, bool suppress_print
options.add_str_i("CUBEPROP_FILEPATH", ".");
/*- Properties to compute. Valid tasks include:
``DENSITY`` - Da, Db, Dt, Ds
``ESP`` - Dt, ESP
``ORBITALS`` - Psi_a_N, Psi_b_N
``BASIS_FUNCTIONS`` - Phi_N
``LOL`` - LOLa, LOLb
``ELF`` - ELFa, ELFb
``DENSITY`` - Da, Db, Dt, Ds;
``ESP`` - Dt, ESP;
``ORBITALS`` - Psi_a_N, Psi_b_N;
``BASIS_FUNCTIONS`` - Phi_N;
``LOL`` - LOLa, LOLb;
``ELF`` - ELFa, ELFb;
``FRONTIER_ORBITALS`` - Psi_a_N_HOMO + Psi_a_N_LUMO;
``DUAL_DESCRIPTOR`` - DUAL_N_HOMO-M_LUMO.
-*/
options.add("CUBEPROP_TASKS", new ArrayType());
/*- List of orbital indices for which cube files are generated (1-based,
View
@@ -45,7 +45,7 @@ foreach(test_name adc1 adc2 casscf-fzc-sp casscf-semi casscf-sa-sp ao-casscf-sp
cc9 cc9a cdomp2-1 cdomp2-2 cepa0-grad1 cepa0-grad2 cepa1
cepa2 cepa3 cepa4 cepa-module ci-multi cisd-h2o+-0 cisd-h2o+-1
cisd-h2o+-2 cisd-h2o-clpse cisd-opt-fd cisd-sp cisd-sp-2
ci-property cubeprop decontract dcft-grad1 dcft-grad2
ci-property cubeprop cubeprop-frontier decontract dcft-grad1 dcft-grad2
dcft-grad3 dcft-grad4 dcft1 dcft2 dcft3 dcft4 dcft5 dcft6
dcft7 dcft8 dcft9 ao-dfcasscf-sp dfcasscf-sa-sp dfcasscf-fzc-sp dfcasscf-sp
dfccd1 dfccdl1 dfccd-grad1 dfccsd1 dfccsdl1 dfccsd-grad1 dfccsd-t-grad1
@@ -0,0 +1,8 @@
include(TestingMacros)
file(COPY O1_HOMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
file(COPY O1_LUMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
file(COPY O3_B2g_SOMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
file(COPY O3_B3g_SOMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
file(COPY O3_DOMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
file(COPY O3_LVMO.cube.ref DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
add_regression_test(cubeprop-frontier "psi;cubeprop")
Oops, something went wrong.

0 comments on commit 2f23252

Please sign in to comment.