Skip to content

Commit

Permalink
McMurchie-Davidson integrals for angular momentum (psi4#2483)
Browse files Browse the repository at this point in the history
  • Loading branch information
@maxscheurer @zachglick
maxscheurer authored and zachglick committed Apr 18, 2022
1 parent f501736 commit dfba9f1
Show file tree
Hide file tree
Showing 10 changed files with 326 additions and 234 deletions.
8 changes: 4 additions & 4 deletions .azure-pipelines/azure-pipelines-linux.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,13 @@ jobs:
strategy:
maxParallel: 8
matrix:
gcc_6:
CXX_COMPILER: 'g++-6'
gcc_7:
CXX_COMPILER: 'g++-7'
PYTHON_VER: '3.8'
C_COMPILER: 'gcc-6'
C_COMPILER: 'gcc-7'
F_COMPILER: 'gfortran'
BUILD_TYPE: 'release'
APT_INSTALL: 'gfortran gcc-6 g++-6'
APT_INSTALL: 'gfortran gcc-7 g++-7'
vmImage: 'ubuntu-18.04'

gcc_latest:
Expand Down
20 changes: 9 additions & 11 deletions cmake/custom_cxxstandard.cmake
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,20 @@
cmake_policy(PUSH)
cmake_policy(SET CMP0057 NEW) # support IN_LISTS

# We require C++14 support from the compiler and standard library.
list(APPEND _allowed_cxx_standards 14 17)
# We require C++17 support from the compiler and standard library.
list(APPEND _allowed_cxx_standards 17)
if(NOT psi4_CXX_STANDARD IN_LIST _allowed_cxx_standards)
message(FATAL_ERROR "Psi4 requires C++14 at least")
message(FATAL_ERROR "Psi4 requires C++17 at least")
endif()

if (CMAKE_CXX_COMPILER_ID MATCHES GNU)
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 6.0)
message(FATAL_ERROR "GCC version must be at least 6.0!")
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 7.0)
message(FATAL_ERROR "GCC version must be at least 7.0!")
endif()

elseif (CMAKE_CXX_COMPILER_ID MATCHES Intel)
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS "17.0.0")
message(FATAL_ERROR "ICPC version must be at least 2017.0.0 to work with pybind11 2.1!") # v1.2
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS "19.0.0")
message(FATAL_ERROR "ICPC version must be at least 2019.0.0 to work with pybind11 2.1!") # v1.2
endif()

set(_testfl ${CMAKE_BINARY_DIR}/test_gcc_version.cc)
Expand Down Expand Up @@ -55,14 +55,12 @@ elseif (CMAKE_CXX_COMPILER_ID MATCHES AppleClang)
endif()

elseif (CMAKE_CXX_COMPILER_ID MATCHES Clang)
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.6)
message(FATAL_ERROR "CLANG version must be at least 3.6!")
if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 5.0)
message(FATAL_ERROR "CLANG version must be at least 5.0!")
endif()

elseif (CMAKE_CXX_COMPILER_ID MATCHES MSVC)
# As for MSVC 14.0, it is not possible to set anything below C++14
# FIXME Remove following line when we switch to C++14
set(psi4_CXX_STANDARD 14)
if(MSVC_TOOLSET_VERSION LESS 140)
message(FATAL_ERROR "MSVC toolset version must be at least 14.0!")
endif()
Expand Down
1 change: 1 addition & 0 deletions psi4/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -288,6 +288,7 @@ install(DIRECTORY ../tests/pytests/
PATTERN "conftest.py"
PATTERN "standard_suite_runner.py"
PATTERN "tdscf_reference_data.json"
PATTERN "oei_reference_data.json"
PATTERN "utils.py")
install(DIRECTORY ../tests/pytests/test_fchk_writer
../tests/pytests/test_molden_writer
Expand Down
1 change: 1 addition & 0 deletions psi4/src/psi4/libmints/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ list(APPEND sources
gaussquad.cc
ecpint.cc
orthog.cc
mcmurchiedavidson.cc
)

# l2intf is a listing of all the files that include libint2's engine.h include. they must
Expand Down
269 changes: 67 additions & 202 deletions psi4/src/psi4/libmints/angularmomentum.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,8 @@
* @END LICENSE
*/

#include "psi4/libciomr/libciomr.h"
#include "psi4/libmints/angularmomentum.h"
#include "psi4/libmints/basisset.h"
#include "psi4/libmints/molecule.h"
#include "psi4/libmints/integral.h"

#include <memory>
Expand All @@ -38,230 +36,97 @@
#include <libint2/shell.h>

using namespace psi;
using namespace mdintegrals;

// Initialize overlap_recur_ to +1 basis set angular momentum, +1 on each center is sufficient
// to compute the dipole derivatives
AngularMomentumInt::AngularMomentumInt(std::vector<SphericalTransform>& spherical_transforms,
std::shared_ptr<BasisSet> bs1, std::shared_ptr<BasisSet> bs2, int nderiv)
: OneBodyAOInt(spherical_transforms, bs1, bs2, nderiv), overlap_recur_(bs1->max_am() + 1, bs2->max_am() + 1) {
int maxam1 = bs1_->max_am();
int maxam2 = bs2_->max_am();
: OneBodyAOInt(spherical_transforms, bs1, bs2, nderiv), MDHelper(bs1->max_am(), bs2->max_am()) {
int maxnao1 = INT_NCART(maxam1_);
int maxnao2 = INT_NCART(maxam2_);

int maxnao1 = (maxam1 + 1) * (maxam1 + 2) / 2;
int maxnao2 = (maxam2 + 1) * (maxam2 + 2) / 2;

// Increase buffer size to handle x, y, and z components
if (deriv_ == 0) {
// one chunk for each Cartesian component
buffer_ = new double[3 * maxnao1 * maxnao2];
set_chunks(3);
} else {
throw PSIEXCEPTION("AngularMomentumInt does not provide derivatives");
throw PSIEXCEPTION("AngularMomentumInt does not provide derivatives.");
}
buffers_.resize(nchunk_);
}

AngularMomentumInt::~AngularMomentumInt() { delete[] buffer_; }


void AngularMomentumInt::compute_pair(const libint2::Shell &s1, const libint2::Shell &s2) {
int ao12;
void AngularMomentumInt::compute_pair(const libint2::Shell& s1, const libint2::Shell& s2) {
// Computes the angular momentum integrals for a pair of shells using eq 9.3.33
// equation numbers from Molecular Electronic-Structure Theory (10.1002/9781119019572)
int am1 = s1.contr[0].l;
int am2 = s2.contr[0].l;
int nprim1 = s1.nprim();
int nprim2 = s2.nprim();
double A[3], B[3];
double Sxy1, Sxy2, Sxz1, Sxz2;
double Syx1, Syx2, Syz1, Syz2;
double Szx1, Szx2, Szy1, Szy2;
double S0x1, S0x2, S0y1, S0y2, S0z1, S0z2;
double muxy1, muxy2, muxz1, muxz2;
double muyx1, muyx2, muyz1, muyz2;
double muzx1, muzx2, muzy1, muzy2;

A[0] = s1.O[0];
A[1] = s1.O[1];
A[2] = s1.O[2];
B[0] = s2.O[0];
B[1] = s2.O[1];
B[2] = s2.O[2];
int am = am1 + am2;

int ydisp = INT_NCART(am1) * INT_NCART(am2);
int zdisp = ydisp + INT_NCART(am1) * INT_NCART(am2);
const auto& comps_am1 = am_comps_[am1];
const auto& comps_am2 = am_comps_[am2];

// compute intermediates
double AB2 = 0.0;
AB2 += (A[0] - B[0]) * (A[0] - B[0]);
AB2 += (A[1] - B[1]) * (A[1] - B[1]);
AB2 += (A[2] - B[2]) * (A[2] - B[2]);
auto A = s1.O;
auto B = s2.O;
// TODO: settable origin? (not supported previously)
Point C = {0.0, 0.0, 0.0};
int nprim1 = s1.nprim();
int nprim2 = s2.nprim();

memset(buffer_, 0, 3 * INT_NCART(am1) * INT_NCART(am2) * sizeof(double));
int dim1 = INT_NCART(am1);
int dim2 = INT_NCART(am2);
int size = dim1 * dim2;
memset(buffer_, 0, 3 * size * sizeof(double));

double** x = overlap_recur_.x();
double** y = overlap_recur_.y();
double** z = overlap_recur_.z();
// dimensions of the E matrix
int edim2 = am2 + 2;
int edim3 = (am1 + 1) + edim2;

int ao12 = 0;
for (int p1 = 0; p1 < nprim1; ++p1) {
double a1 = s1.alpha[p1];
double c1 = s1.contr[0].coeff[p1];
double a = s1.alpha[p1];
double ca = s1.contr[0].coeff[p1];
for (int p2 = 0; p2 < nprim2; ++p2) {
double a2 = s2.alpha[p2];
double c2 = s2.contr[0].coeff[p2];
double gamma = a1 + a2;
double oog = 1.0 / gamma;

double PA[3], PB[3];
double P[3];

P[0] = (a1 * A[0] + a2 * B[0]) * oog;
P[1] = (a1 * A[1] + a2 * B[1]) * oog;
P[2] = (a1 * A[2] + a2 * B[2]) * oog;
PA[0] = P[0] - A[0];
PA[1] = P[1] - A[1];
PA[2] = P[2] - A[2];
PB[0] = P[0] - B[0];
PB[1] = P[1] - B[1];
PB[2] = P[2] - B[2];

double over_pf = exp(-a1 * a2 * AB2 * oog) * sqrt(M_PI * oog) * M_PI * oog * c1 * c2;

// Do recursion
overlap_recur_.compute(PA, PB, gamma, am1 + 1, am2 + 1);

double b = s2.alpha[p2];
double cb = s2.contr[0].coeff[p2];
double p = a + b;
Point P{(a * A[0] + b * B[0]) / p, (a * A[1] + b * B[1]) / p, (a * A[2] + b * B[2]) / p};
auto PC = point_diff(P, C);
double prefac = ca * cb * std::pow(M_PI / p, 1.5);
fill_E_matrix(am1, am2 + 1, P, A, B, a, b, Ex, Ey, Ez);
ao12 = 0;
for (int ii = 0; ii <= am1; ii++) {
int l1 = am1 - ii;
for (int jj = 0; jj <= ii; jj++) {
int m1 = ii - jj;
int n1 = jj;
/*--- create all am components of sj ---*/
for (int kk = 0; kk <= am2; kk++) {
int l2 = am2 - kk;
for (int ll = 0; ll <= kk; ll++) {
int m2 = kk - ll;
int n2 = ll;

double x00 = x[l1][l2], y00 = y[m1][m2], z00 = z[n1][n2];
double x10 = x[l1 + 1][l2], y10 = y[m1 + 1][m2], z10 = z[n1 + 1][n2];
double x01 = x[l1][l2 + 1], y01 = y[m1][m2 + 1], z01 = z[n1][n2 + 1];

Sxy1 = Sxy2 = Sxz1 = Sxz2 = 0.0;

//
// Overlaps
//

// (a+1x|b+1y)
Sxy1 = x10 * y01 * z00 * over_pf;
// (a+1x|b-1y)
if (m2) Sxy2 = x10 * y[m1][m2 - 1] * z00 * over_pf;
// (a+1x|b+1z)
Sxz1 = x10 * y00 * z01 * over_pf;
// (a+1x|b-1z)
if (n2) Sxz2 = x10 * y00 * z[n1][n2 - 1] * over_pf;

// outfile->Printf( "Sxy1 %f Sxy2 %f Sxz1 %f Sxz2 %f\n", Sxy1,
// Sxy2, Sxz1, Sxz2);

Syx1 = Syx2 = Syz1 = Syz2 = 0.0;

// (a+1y|b+1x)
Syx1 = x01 * y10 * z00 * over_pf;
// (a+1y|b-1x)
if (l2) Syx2 = x[l1][l2 - 1] * y10 * z00 * over_pf;
// (a+1y|b+1z)
Syz1 = x00 * y10 * z01 * over_pf;
// (a+1y|b-1z)
if (n2) Syz2 = x00 * y10 * z[n1][n2 - 1] * over_pf;

// outfile->Printf( "Syx1 %f Syx2 %f Syz1 %f Syz2 %f\n", Syx1,
// Syx2, Syz1, Syz2);

Szx1 = Szx2 = Szy1 = Szy2 = 0.0;

// (a+1z|b+1x)
Szx1 = x01 * y00 * z10 * over_pf;
// (a+1z|b-1x)
if (l2) Szx2 = x[l1][l2 - 1] * y00 * z10 * over_pf;
// (a+1z|b+1y)
Szy1 = x00 * y01 * z10 * over_pf;
// (a+1z|b-1y)
if (m2) Szy2 = x00 * y[m1][m2 - 1] * z10 * over_pf;

// outfile->Printf( "Szx1 %f Szx2 %f Szy1 %f Szy2 %f\n", Szx1,
// Szx2, Szy1, Szy2);

S0x1 = S0x2 = S0y1 = S0y2 = S0z1 = S0z2 = 0.0;

// (a|b+1x)
S0x1 = x01 * y00 * z00 * over_pf;
// (a|b-1x)
if (l2) S0x2 = x[l1][l2 - 1] * y00 * z00 * over_pf;
// (a|b+1y)
S0y1 = x00 * y01 * z00 * over_pf;
// (a|b-1y)
if (m2) S0y2 = x00 * y[m1][m2 - 1] * z00 * over_pf;
// (a|b+1z)
S0z1 = x00 * y00 * z01 * over_pf;
// (a|b-1z)
if (n2) S0z2 = x00 * y00 * z[n1][n2 - 1] * over_pf;

// outfile->Printf( "S0x1 %f S0x2 %f S0y1 %f S0y2 %f S0z1
// S0z2\n", S0x1, S0x2, S0y1, S0y2, S0z1, S0z2);

//
// Moment integrals
//

muxy1 = muxy2 = muxz1 = muxz2 = 0.0;
// (a|x|b+1y)
muxy1 = Sxy1 + (A[0] - origin_[0]) * S0y1;
// (a|x|b+1z)
muxz1 = Sxz1 + (A[0] - origin_[0]) * S0z1;
// (a|x|b-1y)
muxy2 = Sxy2 + (A[0] - origin_[0]) * S0y2;
// (a|x|b-1z)
muxz2 = Sxz2 + (A[0] - origin_[0]) * S0z2;

muyx1 = muyx2 = muyz1 = muyz2 = 0.0;
// (a|y|b+1x)
muyx1 = Syx1 + (A[1] - origin_[1]) * S0x1;
// (a|y|b+1z)
muyz1 = Syz1 + (A[1] - origin_[1]) * S0z1;
// (a|y|b-1x)
muyx2 = Syx2 + (A[1] - origin_[1]) * S0x2;
// (a|y|b-1z)
muyz2 = Syz2 + (A[1] - origin_[1]) * S0z2;

muzx1 = muzx2 = muzy1 = muzy2 = 0.0;
// (a|z|b+1x)
muzx1 = Szx1 + (A[2] - origin_[2]) * S0x1;
// (a|z|b+1y)
muzy1 = Szy1 + (A[2] - origin_[2]) * S0y1;
// (a|z|b+1x)
muzx2 = Szx2 + (A[2] - origin_[2]) * S0x2;
// (a|z|b+1y)
muzy2 = Szy2 + (A[2] - origin_[2]) * S0y2;

/* (a|Lx|b) = 2 a2 * (a|(y-Cy)|b+1z) - B.z * (a|(y-Cy)|b-1z)
- 2 a2 * (a|(z-Cz)|b+1y) + B.y * (a|(z-Cz)|b-1y) */
double Lx = (2.0 * a2 * muyz1 - n2 * muyz2 - 2.0 * a2 * muzy1 + m2 * muzy2) /* * over_pf*/;

/* (a|Ly|b) = 2 a2 * (a|(z-Cz)|b+1x) - B.x * (a|(z-Cz)|b-1x)
- 2 a2 * (a|(x-Cx)|b+1z) + B.z * (a|(x-Cx)|b-1z) */
double Ly = (2.0 * a2 * muzx1 - l2 * muzx2 - 2.0 * a2 * muxz1 + n2 * muxz2) /* * over_pf*/;

/* (a|Lz|b) = 2 a2 * (a|(x-Cx)|b+1y) - B.y * (a|(x-Cx)|b-1y)
- 2 a2 * (a|(y-Cy)|b+1x) + B.x * (a|(y-Cy)|b-1x) */
double Lz = (2.0 * a2 * muxy1 - m2 * muxy2 - 2.0 * a2 * muyx1 + l2 * muyx2) /* * over_pf*/;

// Electrons have a negative charge
buffer_[ao12] += (Lx);
buffer_[ao12 + ydisp] += (Ly);
buffer_[ao12 + zdisp] += (Lz);

ao12++;
}
for (const auto& [l1, m1, n1, index1] : comps_am1) {
for (const auto& [l2, m2, n2, index2] : comps_am2) {
int xidx = address_3d(l1, l2, 0, edim2, edim3);
int yidx = address_3d(m1, m2, 0, edim2, edim3);
int zidx = address_3d(n1, n2, 0, edim2, edim3);
double Sx0 = Ex[xidx];
double Sy0 = Ey[yidx];
double Sz0 = Ez[zidx];
// Multipole integral, eq 9.5.43
double Sx1 = Ex[xidx + 1] + PC[0] * Ex[xidx];
double Sy1 = Ey[yidx + 1] + PC[1] * Ey[yidx];
double Sz1 = Ez[zidx + 1] + PC[2] * Ez[zidx];
// Use eq 9.3.30 to define derivative w.r.t. ket center
double Dx1 = -2.0 * b * Ex[address_3d(l1, l2 + 1, 0, edim2, edim3)];
double Dy1 = -2.0 * b * Ey[address_3d(m1, m2 + 1, 0, edim2, edim3)];
double Dz1 = -2.0 * b * Ez[address_3d(n1, n2 + 1, 0, edim2, edim3)];
if (l2 > 0) {
Dx1 += l2 * Ex[address_3d(l1, l2 - 1, 0, edim2, edim3)];
}
if (m2 > 0) {
Dy1 += m2 * Ey[address_3d(m1, m2 - 1, 0, edim2, edim3)];
}
if (n2 > 0) {
Dz1 += n2 * Ez[address_3d(n1, n2 - 1, 0, edim2, edim3)];
}
// Lx
buffer_[ao12 + size * 0] += (Sx0 * Dy1 * Sz1 - Sx0 * Sy1 * Dz1) * prefac;
// Ly
buffer_[ao12 + size * 1] += (Sx1 * Sy0 * Dz1 - Dx1 * Sy0 * Sz1) * prefac;
// Lz
buffer_[ao12 + size * 2] += (Dx1 * Sy1 * Sz0 - Sx1 * Dy1 * Sz0) * prefac;
++ao12;
}
}
}
Expand All @@ -270,4 +135,4 @@ void AngularMomentumInt::compute_pair(const libint2::Shell &s1, const libint2::S
buffers_[0] = buffer_ + 0 * s1.size() * s2.size();
buffers_[1] = buffer_ + 1 * s1.size() * s2.size();
buffers_[2] = buffer_ + 2 * s1.size() * s2.size();
}
}
Loading

0 comments on commit dfba9f1

Please sign in to comment.