Skip to content
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

B-field derivatives #57

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

B-field derivatives #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4f2bf6f
added GAIO B-field derivatives for a single shell-pair
Jan 24, 2023
77e2247
B-field derivative
Jan 26, 2023
6326dc4
cleanup
Jan 26, 2023
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
35 changes: 26 additions & 9 deletions include/libecpint/ecpint.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,20 +148,37 @@ namespace libecpint {
TwoIndex<double> &values, int shiftA = 0, int shiftB = 0) const;

/**
* Computes the overall ECP integral first derivatives over the given ECP center, C, and shell pair (A | B)
* Computes the overall ECP integral first derivatives over the given ECP center, C, and shell pair (A | B)
* The results are placed in order [Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz] and are calculated so that each component
* can always be added to the relevant total derivative. E.g. if A = B, then the contribution to the total derivative
* for that coordinate on the x axis will be Ax + Bx.
* The order for each derivative matrices matches that specified in compute_shell_pair
*
* @param U - reference to the ECP
* @param shellA - the first basis shell (rows in values)
* @param shellB - the second basis shell (cols in values)
* @param results - reference to array of 9 TwoIndex arrays where the results will be stored
*/
*
* @param U - reference to the ECP
* @param shellA - the first basis shell (rows in values)
* @param shellB - the second basis shell (cols in values)
* @param results - reference to array of 9 TwoIndex arrays where the results will be stored
*/
void compute_shell_pair_derivative(
const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB,
std::array<TwoIndex<double>, 9> &results) const;

/**
* Computes the overall ECP integral first derivatives wrt an external magnetic field over the
* given ECP center, C, and GIAO shell pair (A | B).
* The results are placed in order [dECP/dBx, dECP/dBy, dECP/dBz].
* The order for each derivative matrices matches that specified in compute_shell_pair
*
* Note that the resulting integrals are imaginary, i.e., the user has to take care of skew-symmetry.
*
* @param U - reference to the ECP
* @param shellA - the first basis shell (rows in values)
* @param shellB - the second basis shell (cols in values)
* @param results - reference to array of 9 TwoIndex arrays where the results will be stored
*/
void compute_shell_pair_derivative_Bfield(
const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB,
std::array<TwoIndex<double>, 3> &results) const;

/**
* Computes the overall ECP integral second derivatives over the given ECP center, C, and shell pair (A | B)
Expand All @@ -171,7 +188,7 @@ namespace libecpint {
* especially in the instance where A=B. It's recommended to look at the compute_second_derivatives interface in api.cpp for how
* to handle this.
* The order for each derivative matrices matches that specified in compute_shell_pair
*
*
* @param U - reference to the ECP
* @param shellA - the first basis shell (rows in values)
* @param shellB - the second basis shell (cols in values)
Expand All @@ -185,7 +202,7 @@ namespace libecpint {
* Worker function to calculate the derivative of the integral <A | C | B> with respect to A.
* This is given as <d_q A(l_q) | C | B> = l_q*<A(l_q-1) | C | B> - 2*mu*<A(l_q + 1) | C | B>
* where l_q is the angular momentum component of A in the q coordinate, and mu is the exponent of A.
*
*
* @param U - reference to the ECP
* @param shellA - the first basis shell (rows in values)
* @param shellB - the second basis shell (cols in values)
Expand Down
2 changes: 1 addition & 1 deletion include/libecpint/gshell.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ namespace libecpint {
void addPrim(double exp, double c);

/// @return the number of primitives
int nprimitive() const { return exps.size(); }
int nprimitive() const { return (int)exps.size(); }

/// @return the number of cartesian basis functions in a shell with this angular momentum
int ncartesian() const { return ((l+1)*(l+2))/2; }
Expand Down
80 changes: 80 additions & 0 deletions src/lib/ecpint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -697,6 +697,86 @@ namespace libecpint {
}
}

void ECPIntegral::compute_shell_pair_derivative_Bfield(
const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB,
std::array<TwoIndex<double>, 3> &results) const {
// First we check centres
double A[3], B[3], C[3], AmC[3], BmC[3];
for (int i = 0; i < 3; i++) {
A[i] = shellA.center()[i];
B[i] = shellB.center()[i];
C[i] = U.center()[i];
AmC[i] = A[i] - C[i];
BmC[i] = B[i] - C[i];
}

int lqnA = shellA.am();
int lqnB = shellB.am();

int ncartA = (lqnA+1) * (lqnA+2) / 2;
int ncartB = (lqnB+1) * (lqnB+2) / 2;

results[0].assign(ncartA, ncartB, 0.0);
results[1].assign(ncartA, ncartB, 0.0);
results[2].assign(ncartA, ncartB, 0.0);

if ((std::abs(A[0] - B[0]) + std::abs(A[1] - B[1]) + std::abs(A[2] - B[2])) < 1e-5){ // same centers, integrals are zero...
return;
}

int ncartAp1 = (lqnA+2) * (lqnA+3) / 2;
int ncartBp1 = (lqnB+2) * (lqnB+3) / 2;

double AxC[3];
AxC[0] = A[1]*C[2] - A[2]*C[1];
AxC[1] = A[2]*C[0] - A[0]*C[2];
AxC[2] = A[0]*C[1] - A[1]*C[0];
double BxC[3];
BxC[0] = B[1]*C[2] - B[2]*C[1];
BxC[1] = B[2]*C[0] - B[0]*C[2];
BxC[2] = B[0]*C[1] - B[1]*C[0];

// Calculate shell derivatives
TwoIndex<double> ABints, ApBints, ABpints;
compute_shell_pair(U,shellA,shellB,ABints); // regular ints
compute_shell_pair(U,shellA,shellB,ApBints,1,0); // A-shell increased l-qn by 1
compute_shell_pair(U,shellA,shellB,ABpints,0,1); // B-shell increased l-qn by 1

int iB = 0;
for(int LxB=lqnB;LxB>=0;LxB--){
for(int LyB=lqnB-LxB;LyB>=0;LyB--){
const int LzB = lqnB-LxB-LyB;
int iA = 0;
for(int LxA=lqnA;LxA>=0;LxA--){
for(int LyA=lqnA-LxA;LyA>=0;LyA--){
const int LzA = lqnA-LxA-LyA;

const int pos0 = ncartB*iA + iB;

const int posAx = ncartB*iA + iB; // LxA+1
const int posAy = ncartB*(LzA + ((LyA+LzA+2)*(LyA+LzA+1))/2) + iB; // LyA+1
const int posAz = posAy + ncartB; // LzA+1

const int posBx = ncartBp1*iA + iB; // LxB+1
const int posBy = ncartBp1*iA + (LzB + ((LyB+LzB+2)*(LyB+LzB+1))/2); // LyB+1
const int posBz = posBy + 1; // LzB+1

results[0].data[pos0] = AxC[0] * ABints.data[pos0] + AmC[1] * ApBints.data[posAz] - AmC[2] * ApBints.data[posAy];
results[1].data[pos0] = AxC[1] * ABints.data[pos0] + AmC[2] * ApBints.data[posAx] - AmC[0] * ApBints.data[posAz];
results[2].data[pos0] = AxC[2] * ABints.data[pos0] + AmC[0] * ApBints.data[posAy] - AmC[1] * ApBints.data[posAx];
results[0].data[pos0] -= BxC[0] * ABints.data[pos0] + BmC[1] * ABpints.data[posBz] - BmC[2] * ABpints.data[posBy];
results[1].data[pos0] -= BxC[1] * ABints.data[pos0] + BmC[2] * ABpints.data[posBx] - BmC[0] * ABpints.data[posBz];
results[2].data[pos0] -= BxC[2] * ABints.data[pos0] + BmC[0] * ABpints.data[posBy] - BmC[1] * ABpints.data[posBx];

iA++;
}
}
iB++;
}
}
}


void ECPIntegral::compute_shell_pair_second_derivative(
const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB,
std::array<TwoIndex<double>, 45> &results) const {
Expand Down