diff --git a/include/libecpint/ecpint.hpp b/include/libecpint/ecpint.hpp index 27967c6..77246cc 100644 --- a/include/libecpint/ecpint.hpp +++ b/include/libecpint/ecpint.hpp @@ -148,20 +148,37 @@ namespace libecpint { TwoIndex &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, 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, 3> &results) const; /** * Computes the overall ECP integral second derivatives over the given ECP center, C, and shell pair (A | B) @@ -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) @@ -185,7 +202,7 @@ namespace libecpint { * Worker function to calculate the derivative of the integral with respect to A. * This is given as = l_q* - 2*mu* * 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) diff --git a/include/libecpint/gshell.hpp b/include/libecpint/gshell.hpp index 08e73d1..86b3acf 100644 --- a/include/libecpint/gshell.hpp +++ b/include/libecpint/gshell.hpp @@ -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; } diff --git a/src/lib/ecpint.cpp b/src/lib/ecpint.cpp index 7737f56..340d090 100644 --- a/src/lib/ecpint.cpp +++ b/src/lib/ecpint.cpp @@ -697,6 +697,86 @@ namespace libecpint { } } + void ECPIntegral::compute_shell_pair_derivative_Bfield( + const ECP &U, const GaussianShell &shellA, const GaussianShell &shellB, + std::array, 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 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, 45> &results) const {