Skip to content

Commit

Permalink
E(30)exch-ind in SAPT2+3 without the S^2 approximation (psi4#2314)
Browse files Browse the repository at this point in the history 
* First implementation of nonapproximated E(30)exch-ind, based on Jonathan Waldrop's Psi4NumPy code.

* Reduced the number of JK builds for nonexpanded E(30)exch-ind

* Code cleanup

* More code cleanup

* Revert the inclusion of E(30)exch-ind,resp(S^All) in totals - back to regular S^2 in totals

* Added a sapt-exch-ind30-inf test.

* Further cleanup

* Furtherer cleanup

* Added a short documentation of the new method to the manual.

* Removed the unused threading variable

* Small fixes in the documentation

Co-authored-by: kjp0013 <kjp0013@c20-login01.cm.cluster>
  • Loading branch information
@konpat @zachglick
2 people authored and zachglick committed Apr 18, 2022
1 parent d156ea2 commit caffbed
Show file tree
Hide file tree
Showing 12 changed files with 3,346 additions and 24 deletions.
13 changes: 13 additions & 0 deletions doc/sphinxman/source/bibliography.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -568,3 +568,16 @@ Bibliography
T. H. Thompson and C. Ochsenfeld
*J. Chem. Phys.* **147**, 144101 (2017).
doi: 10.1063/1.4994190
.. [Smith:2020:184108]
"PSI4 1.4: Open-source software for high-throughput quantum chemistry",
D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O'Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer III, A. Yu. Sokolov, K. Patkowski, A. E. DePrince III, U. Bozkaya, R. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill
*J. Chem. Phys.* **152**, 184108 (2020).
doi: 10.1063/5.0006002
.. [Waldrop:2021:024103]
"Nonapproximated third-order exchange induction energy in symmetry-adapted perturbation theory",
J. M. Waldrop and K. Patkowski
*J. Chem. Phys.* **154**, 024103 (2021).
doi: 10.1063/1.4994190
24 changes: 18 additions & 6 deletions doc/sphinxman/source/sapt.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -917,28 +917,40 @@ following publications: [Patkowski:2018:164110]_

.. _`sec:saptinf`:

Second-Order Exchange Terms without Single-Exchange Approximation
Higher-Order Exchange Terms without Single-Exchange Approximation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Recently, the SAPT second-order exchange terms have been derived without
the :math:`S^{2}` approximation in the works [Schaffer:2012:1235]_ and
[Schaffer:2013:2570]_. These new terms can be computed with the following
Recently, several SAPT higher-order exchange terms have been derived without
the :math:`S^{2}` approximation: :math:`E_{exch-ind}^{(20)}` [Schaffer:2012:1235]_,
:math:`E_{exch-disp}^{(20)}` [Schaffer:2013:2570]_, and :math:`E_{exch-ind}^{(30)}`
[Waldrop:2021:024103]_. The second-order terms can be computed with the following
settings::

set SAPT_DFT_FUNCTIONAL HF
set DO_IND_EXCH_SINF true # calculate Exch-Ind20 (S^inf)
set SAPT_DFT_MP2_DISP_ALG fisapt
set DO_DISP_EXCH_SINF true # calculate Exch-Disp20 (S^inf)
energy('sapt(dft)')
and the third-order exchange-induction term is computed as follows::

set DO_IND30_EXCH_SINF true # calculate Exch-Ind30 (S^inf)
energy('sapt2+3')
These calculations are performed with the atomic orbital and
density-fitting scheme of [J. M. Waldrop et al., to be published].
density-fitting scheme described in the Supplementary Material to
[Smith:2020:184108]_ for the second-order terms and in [Waldrop:2021:024103]_
for the third-order exchange induction. The coupled (response) version of the
exchange-induction corrections are also calculated, exactly for
:math:`E_{exch-ind,resp}^{(20)}` and by scaling the uncoupled term for
:math:`E_{exch-ind,resp}^{(30)}`.

S^inf Keywords
~~~~~~~~~~~~~~

.. include:: autodir_options_c/sapt__do_ind_exch_sinf.rst
.. include:: autodir_options_c/sapt__do_disp_exch_sinf.rst
.. include:: autodir_options_c/sapt__do_ind30_exch_sinf.rst

.. _`sec:saptd`:

Expand Down
257 changes: 257 additions & 0 deletions psi4/src/psi4/libsapt_solver/exch-ind30.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,10 @@
#include "psi4/libpsio/psio.hpp"
#include "psi4/libqt/qt.h"
#include "psi4/libpsi4util/PsiOutStream.h"
#include "psi4/libfock/jk.h"
#include "psi4/libmints/basisset.h"
#include "psi4/libmints/matrix.h"
#include "psi4/libpsi4util/process.h"

namespace psi {
namespace sapt {
Expand Down Expand Up @@ -314,5 +318,258 @@ double SAPT2p3::exch_ind30_3(double **sBS) {

return (energy);
}

// Compute the nonapproximated third-order exchange-induction //
// Konrad Patkowski, based on Jonathan Waldrop's Psi4NumPy code //
// September 2021 //
void SAPT2p3::sinf_e30ind() {

if (print_) {
outfile->Printf(" ==> Nonapproximated third-order induction <==\n\n");
}

// => Sizing <= //

int nn = basisset_->nbf();
int na = noccA_;
int nb = noccB_;
int nr = nvirA_;
int ns = nvirB_;

auto uAR = Matrix("Ind30 uAR Amplitudes", na, nr);
auto uBS = Matrix("Ind30 uBS Amplitudes", nb, ns);
uAR.load(psio_, PSIF_SAPT_AMPS, Matrix::SaveType::Full);
uBS.load(psio_, PSIF_SAPT_AMPS, Matrix::SaveType::Full);

// => Intermolecular overlap matrix and inverse <= //
auto Sab = linalg::triplet(CoccA_, Smat_, CoccB_, true, false, false);

double** Sabp = Sab->pointer();
auto D = Matrix("D", na + nb, na + nb);
D.identity();
double** Dp = D.pointer();
for (int a = 0; a < na; a++) {
for (int b = 0; b < nb; b++) {
Dp[a][b + na] = Dp[b + na][a] = Sabp[a][b];
}
}
D.power(-1.0, 1.0E-12);

Dimension zero(1);
Dimension nadim(1);
Dimension nabdim(1);
zero[0] = 0;
nadim[0] = na;
nabdim[0] = na + nb;
auto Daa = D.get_block(Slice(zero, nadim), Slice(zero, nadim));
auto Dab = D.get_block(Slice(zero, nadim), Slice(nadim, nabdim));
auto Dbb = D.get_block(Slice(nadim, nabdim), Slice(nadim, nabdim));

// => New Stuff <= //
// Start with T's
auto Sbr = linalg::triplet(CoccB_, Smat_, CvirA_, true, false, false);
auto Sas = linalg::triplet(CoccA_, Smat_, CvirB_, true, false, false);
auto Tar = linalg::doublet(Dab, Sbr, false, false);
auto Tbr = linalg::doublet(Dbb, Sbr, false, false);
auto Tas = linalg::doublet(Daa, Sas, false, false);
auto Tbs = linalg::doublet(Dab, Sas, true, false);

// C1's and C2's from D's and C's.
// C1's are C times D diagonal blocks.
// C2's are times off-diagonal blocks.
auto C1a = linalg::doublet(CoccA_, Daa, false, false);
auto C1b = linalg::doublet(CoccB_, Dbb, false, false);
auto C2a = linalg::doublet(CoccB_, Dab, false, true);
auto C2b = linalg::doublet(CoccA_, Dab, false, false);

// Coeffs for all occupied
std::vector<std::shared_ptr<Matrix> > hold_these;
hold_these.push_back(CoccA_);
hold_these.push_back(CoccB_);

auto Cocc0AB = linalg::horzcat(hold_these);
hold_these.clear();

// Half transform D_ia and D_ib for JK
auto D_Ni_a = std::make_shared<Matrix>("D_Ni_a", nn, na + nb);
auto D_Ni_b = std::make_shared<Matrix>("D_Ni_b", nn, na + nb);

C_DGEMM('N', 'N', nn, na + nb, na, 1.0, CoccA_->pointer()[0], na, &Dp[0][0], na + nb, 0.0,
D_Ni_a->pointer()[0], na + nb);
C_DGEMM('N', 'N', nn, na + nb, nb, 1.0, CoccB_->pointer()[0], nb, &Dp[na][0], na + nb, 0.0,
D_Ni_b->pointer()[0], na + nb);

// Global JK object
std::shared_ptr<JK> jk_;

// Make JK's
jk_ = JK::build_JK(basisset_, get_basisset("DF_BASIS_SCF"), options_, false, mem_);
jk_->set_memory(mem_);
jk_->set_do_J(true);
jk_->set_do_K(true);
jk_->initialize();
jk_->print_header();

auto& Cl = jk_->C_left();
auto& Cr = jk_->C_right();
const auto& J = jk_->J();
const auto& K = jk_->K();

Cl.clear();
Cr.clear();
Cl.push_back(Cocc0AB);
Cr.push_back(D_Ni_a);
Cl.push_back(Cocc0AB);
Cr.push_back(D_Ni_b);
jk_->compute();

auto J_D_ia = J[0];
auto K_D_ia = K[0];

auto J_D_ib = J[1];
auto K_D_ib = K[1];

// Finish D_ia and D_ib transformation to make tilded C's
auto D_ia = linalg::doublet(Cocc0AB, D_Ni_a, false, true);
auto D_ib = linalg::doublet(Cocc0AB, D_Ni_b, false, true);

auto Ct_Kr = linalg::triplet(D_ib, Smat_, CvirA_, false, false, false);
Ct_Kr->scale(-1);
Ct_Kr->add(CvirA_);
auto Ct_Ks = linalg::triplet(D_ia, Smat_, CvirB_, false, false, false);
Ct_Ks->scale(-1);
Ct_Ks->add(CvirB_);

// Make omega cores
std::shared_ptr<Matrix> AJK(J_D_ia->clone());
AJK->scale(2);
AJK->add(VAmat_);
AJK->subtract(K_D_ia);

auto AJK_ar = linalg::triplet(C2a, AJK, Ct_Kr, true, false, false);
auto AJK_as = linalg::triplet(C2a, AJK, Ct_Ks, true, false, false);
auto AJK_br = linalg::triplet(C1b, AJK, Ct_Kr, true, false, false);
auto AJK_bs = linalg::triplet(C1b, AJK, Ct_Ks, true, false, false);

std::shared_ptr<Matrix> BJK(J_D_ib->clone());
BJK->scale(2);
BJK->add(VBmat_);
BJK->subtract(K_D_ib);

auto BJK_ar = linalg::triplet(C1a, BJK, Ct_Kr, true, false, false);
auto BJK_as = linalg::triplet(C1a, BJK, Ct_Ks, true, false, false);
auto BJK_br = linalg::triplet(C2b, BJK, Ct_Kr, true, false, false);
auto BJK_bs = linalg::triplet(C2b, BJK, Ct_Ks, true, false, false);

// Finish omega terms
Matrix omega_ar(AJK_ar->clone());
omega_ar.add(BJK_ar);
omega_ar.scale(2);

Matrix omega_as(AJK_as->clone());
omega_as.add(BJK_as);
omega_as.scale(2);

Matrix omega_br(AJK_br->clone());
omega_br.add(BJK_br);
omega_br.scale(2);

Matrix omega_bs(AJK_bs->clone());
omega_bs.add(BJK_bs);
omega_bs.scale(2);

Cocc0AB.reset();
D_Ni_a.reset();
D_Ni_b.reset();
J_D_ia.reset();
K_D_ia.reset();
J_D_ib.reset();
K_D_ib.reset();
AJK.reset();
BJK.reset();
AJK_ar.reset();
AJK_as.reset();
AJK_br.reset();
AJK_bs.reset();
BJK_ar.reset();
BJK_as.reset();
BJK_br.reset();
BJK_bs.reset();

//Konrad - all the stuff above still needed but I adjusted the scale factors in \omega's
//we don't need the DF stuff that was below

double CompleteInd30 = 0.0;

auto sAR = std::make_shared<Matrix>("sAR", na, nr);
double **sARp = sAR->pointer();

for (int a = 0; a < na; a++) {
for (int r = 0; r < nr; r++) {
sARp[a][r] = wBAR_[a][r] / (evalsA_[a] - evalsA_[na+r]);
}
}

auto sBS = std::make_shared<Matrix>("sBS", nb, ns);
double **sBSp = sBS->pointer();

for (int b = 0; b < nb; b++) {
for (int s = 0; s < ns; s++) {
sBSp[b][s] = wABS_[b][s] / (evalsB_[b] - evalsB_[nb+s]);
}
}

auto STS_br = linalg::triplet(sBS, Tas, sAR, false, true, false);
auto STS_as = linalg::triplet(sAR, Tbr, sBS, false, true, false);
auto STS_ar = linalg::triplet(sAR, Tar, sAR, false, true, false);
auto STS_bs = linalg::triplet(sBS, Tbs, sBS, false, true, false);

CompleteInd30 += uAR.vector_dot(omega_ar);
CompleteInd30 += uBS.vector_dot(omega_bs);
CompleteInd30 -= STS_br->vector_dot(omega_br);
CompleteInd30 -= STS_as->vector_dot(omega_as);
CompleteInd30 -= STS_ar->vector_dot(omega_ar);
CompleteInd30 -= STS_bs->vector_dot(omega_bs);

//All Omega-dependent contributions have been completed, now Xi-dependent (J/K) contributions

auto SCt_Na = linalg::doublet(Ct_Kr, sAR, false, true);
auto SCt_Nb = linalg::doublet(Ct_Ks, sBS, false, true);

auto preXiAA = linalg::doublet(SCt_Na, Daa, false, false);
auto preXiBA = linalg::doublet(SCt_Na, Dab, false, false);
auto preXiAB = linalg::doublet(SCt_Nb, Dab, false, true);
auto preXiBB = linalg::doublet(SCt_Nb, Dbb, false, false);

auto XiAA = linalg::doublet(CoccA_, preXiAA, false, true);
auto XiAB = linalg::doublet(CoccA_, preXiAB, false, true);
auto XiBA = linalg::doublet(CoccB_, preXiBA, false, true);
auto XiBB = linalg::doublet(CoccB_, preXiBB, false, true);

preXiBB->add(preXiBA); //enough to calculate JK for the sum of the two

Cl.clear();
Cr.clear();
Cl.push_back(CoccB_);
Cr.push_back(preXiBB);
jk_->compute();

std::shared_ptr<Matrix> J_XiBB = J[0];
std::shared_ptr<Matrix> K_XiBB = K[0]->transpose();

CompleteInd30 += 4.0 * XiAA->vector_dot(J_XiBB);
CompleteInd30 -= 2.0 * XiAA->vector_dot(K_XiBB);
CompleteInd30 += 4.0 * XiAB->vector_dot(J_XiBB);
CompleteInd30 -= 2.0 * XiAB->vector_dot(K_XiBB);

e_exch_ind30_sinf_ = CompleteInd30 - e_ind30_;

Process::environment.globals["SAPT EXCH-IND30(S^INF) ENERGY"] = e_exch_ind30_sinf_;
if (print_) {
outfile->Printf(" Exch-Ind30 (S^inf) = %18.12lf [Eh]\n", e_exch_ind30_sinf_);
outfile->Printf("\n");
}
}

} // namespace sapt
} // namespace psi
24 changes: 15 additions & 9 deletions psi4/src/psi4/libsapt_solver/sapt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,10 +217,10 @@ void SAPT::initialize(SharedWavefunction MonomerA, SharedWavefunction MonomerB)
auto intfact = std::make_shared<IntegralFactory>(basisset_, basisset_, basisset_, basisset_);

std::shared_ptr<OneBodyAOInt> Sint(intfact->ao_overlap());
auto Smat = std::make_shared<Matrix>(fact->create_matrix("Overlap"));
Sint->compute(Smat);
Smat_ = std::make_shared<Matrix>(fact->create_matrix("Overlap"));
Sint->compute(Smat_);

double **sIJ = Smat->pointer();
double **sIJ = Smat_->pointer();
double **sAJ = block_matrix(nmoA_, nso_);
sAB_ = block_matrix(nmoA_, nmoB_);

Expand All @@ -245,8 +245,8 @@ void SAPT::initialize(SharedWavefunction MonomerA, SharedWavefunction MonomerB)
}
}
potA->set_charge_field(ZxyzA);
auto VAmat = std::make_shared<Matrix>(fact->create_matrix("Nuclear Attraction (Monomer A)"));
potA->compute(VAmat);
VAmat_ = std::make_shared<Matrix>(fact->create_matrix("Nuclear Attraction (Monomer A)"));
potA->compute(VAmat_);

auto potB = std::shared_ptr<PotentialInt>(dynamic_cast<PotentialInt *>(intfact->ao_potential()));
auto ZxyzB = std::make_shared<Matrix>("Charges B (Z,x,y,z)", natomsB_, 4);
Expand All @@ -264,8 +264,8 @@ void SAPT::initialize(SharedWavefunction MonomerA, SharedWavefunction MonomerB)
}
}
potB->set_charge_field(ZxyzB);
auto VBmat = std::make_shared<Matrix>(fact->create_matrix("Nuclear Attraction (Monomer B)"));
potB->compute(VBmat);
VBmat_ = std::make_shared<Matrix>(fact->create_matrix("Nuclear Attraction (Monomer B)"));
potB->compute(VBmat_);

double **vIB = block_matrix(nso_, nmoB_);
double **vAJ = block_matrix(nmoA_, nso_);
Expand All @@ -274,20 +274,26 @@ void SAPT::initialize(SharedWavefunction MonomerA, SharedWavefunction MonomerB)
vBAA_ = block_matrix(nmoA_, nmoA_);
vBAB_ = block_matrix(nmoA_, nmoB_);

double **vIJ = VAmat->pointer();
double **vIJ = VAmat_->pointer();

C_DGEMM('N', 'N', nso_, nmoB_, nso_, 1.0, vIJ[0], nso_, CB_[0], nmoB_, 0.0, vIB[0], nmoB_);
C_DGEMM('T', 'N', nmoA_, nmoB_, nso_, 1.0, CA_[0], nmoA_, vIB[0], nmoB_, 0.0, vAAB_[0], nmoB_);
C_DGEMM('T', 'N', nmoB_, nmoB_, nso_, 1.0, CB_[0], nmoB_, vIB[0], nmoB_, 0.0, vABB_[0], nmoB_);

vIJ = VBmat->pointer();
vIJ = VBmat_->pointer();

C_DGEMM('T', 'N', nmoA_, nso_, nso_, 1.0, CA_[0], nmoA_, vIJ[0], nso_, 0.0, vAJ[0], nso_);
C_DGEMM('N', 'N', nmoA_, nmoA_, nso_, 1.0, vAJ[0], nso_, CA_[0], nmoA_, 0.0, vBAA_[0], nmoA_);
C_DGEMM('N', 'N', nmoA_, nmoB_, nso_, 1.0, vAJ[0], nso_, CB_[0], nmoB_, 0.0, vBAB_[0], nmoB_);

free_block(vIB);
free_block(vAJ);

//Konrad: extra stuff for nonapproximate E(30)ex-ind in AOs
CoccA_ = MonomerA->Ca_subset("AO", "OCC");
CoccB_ = MonomerB->Ca_subset("AO", "OCC");
CvirA_ = MonomerA->Ca_subset("AO", "VIR");
CvirB_ = MonomerB->Ca_subset("AO", "VIR");
}

void SAPT::get_denom() {
Expand Down
Loading

0 comments on commit caffbed

Please sign in to comment.