FreemanTheMaverick
commented
1 year ago Now I shorten the snippet so that it becomes easier to read. It is supposed to compute all second-order overlap derivatives and store them in a 2D array of matrices (i.e., a dimension of 3Natoms×3Natoms×Nbasis×Nbasis.
#include <Eigen/Dense>
#include <libint2.hpp>
#include <iostream>
#define EigenMatrix Eigen::MatrixXd
#define EigenZero Eigen::MatrixXd::Zero
int main(){
// Initialization of a water molecule
std::vector<libint2::Atom> libint2atoms(3);
libint2atoms[0].atomic_number=8;
libint2atoms[0].x=0.000000;
libint2atoms[0].y=0.000000;
libint2atoms[0].z=0.000000;
libint2atoms[1].atomic_number=1;
libint2atoms[1].x=0.000000;
libint2atoms[1].y=-1.430523;
libint2atoms[1].z=1.109269;
libint2atoms[2].atomic_number=1;
libint2atoms[2].x=0.000000;
libint2atoms[2].y=1.430523;
libint2atoms[2].z=1.109269;
libint2::BasisSet obs("sto-3g",libint2atoms);
int nbasis=0;
for (const auto& shell:obs) nbasis+=shell.size();
// Initialization of libint2
EigenMatrix ovl_hess[3*3][3*3];
for (int i=0;i<3*3;i++)
for (int j=0;j<3*3;j++)
ovl_hess[i][j]=EigenZero(nbasis,nbasis);
int atomlist[]={0,0};
libint2::initialize();
libint2::Engine ovl_engine(libint2::Operator::overlap,obs.max_nprim(),obs.max_l(),2);
const auto & ovl_buf=ovl_engine.results();
auto shell2bf=obs.shell2bf();
auto shell2atom=obs.shell2atom(libint2atoms);
for (short int s1=0;s1!=(short int)obs.size();s1++){ // Looping over unique shell pairs
atomlist[0]=shell2atom[s1];
const short int bf1_first=shell2bf[s1];
const short int n1=obs[s1].size();
for (short int s2=0;s2<=s1;s2++){
atomlist[1]=shell2atom[s2];
const short int bf2_first=shell2bf[s2];
const short int n2=obs[s2].size();
// Overlap.
ovl_engine.compute(obs[s1],obs[s2]);
if (ovl_buf[0]){
for (short int f1=0,f12=0;f1!=n1;f1++){
const short int bf1=bf1_first+f1;
for (short int f2=0;f2!=n2;f2++,f12++){
const short int bf2=bf2_first+f2;
if (bf2<=bf1){ // Looping over unique basis function pairs
int xpert; // The x^th coordinate of hessian
int ypert; // The y^th coordinate of hessian
for (int p=0,ptqs=0;p<2;p++) for (int t=0;t<3;t++){
xpert=3*atomlist[p]+t;
for (int q=p;q<2;q++) for (int s=((q==p)?t:0);s<3;s++,ptqs++){ // Looping over unique coordinate pairs
ypert=3*atomlist[q]+s;
ovl_hess[xpert][ypert](bf1,bf2)=ovl_hess[xpert][ypert](bf2,bf1)=ovl_hess[ypert][xpert](bf1,bf2)=ovl_hess[ypert][xpert](bf2,bf1)=ovl_buf[ptqs][f12];
}
}
}
}
}
}
}
}
libint2::finalize(); // Printing the overlap derivatives
std::cout.setf(std::ios::fixed);
std::cout<<std::setprecision(6);
for (int xpert=0;xpert<3*3;xpert++)
for (int ypert=xpert;ypert<3*3;ypert++){
std::cout<<xpert<<" "<<ypert<<std::endl;
std::cout<<ovl_hess[xpert][ypert]<<std::endl;
}
return 0;
}However, the problem still exists. Several diagonal or near-diagonal elements are wrong in the matrices whose two perturbations are about the same atom (i.e., d^2/dXdY where X and Y are coordinates of the same atom). Some incorrect results are presented as follows along with the corresponding correct ones given by hartree-fock++.cc.
X=0, Y=0, x coordinate of the first atom. A few diagonal elements (blocks) are wrong.
My code
-19.335467 0.112007 0.000000 0.000000 0.000000 -0.037385 -0.037385
0.112007 -0.538752 0.000000 0.000000 0.000000 -0.191744 -0.191744
0.000000 0.000000 -3.034477 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 -1.011492 0.000000 0.175216 -0.175216
0.000000 0.000000 0.000000 0.000000 -1.011492 -0.135867 -0.135867
-0.037385 -0.191744 0.000000 0.175216 -0.135867 0.000000 0.000000
-0.037385 -0.191744 0.000000 -0.175216 -0.135867 0.000000 0.000000
hartree-fock++.cc
-0.000000 0.000000 0.000000 0.000000 0.000000 -0.037385 -0.037385
0.000000 0.000000 0.000000 0.000000 0.000000 -0.191744 -0.191744
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.175216 -0.175216
0.000000 0.000000 0.000000 0.000000 -0.000000 -0.135868 -0.135868
-0.037385 -0.191744 0.000000 0.175216 -0.135868 0.000000 0.000000
-0.037385 -0.191744 0.000000 -0.175216 -0.135868 0.000000 0.000000X=0, Y=1, x and y of the first atom. A pair of symmetric off-diagonal elements in the middle are wrong.
My code
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 -1.011492 0.000000 0.175216 -0.175216
0.000000 0.000000 -1.011492 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.175216 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -0.175216 0.000000 0.000000 0.000000 0.000000
hartree-fock++.cc
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.175216 -0.175216
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.175216 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -0.175216 0.000000 0.000000 0.000000 0.000000X=8, Y=8, z of the last atom. The last diagonal element is wrong.
My code
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.006654
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.072069
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.014628
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.283078
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.066909
0.006654 -0.072069 0.000000 -0.014628 -0.283078 -0.066909 -0.506688
hartree-fock++.cc
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.006654
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.072069
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.014628
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.283078
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.066909
0.006654 -0.072069 0.000000 -0.014628 -0.283078 -0.066909 0.000000I must have mistaken the meaning of what is in the integral engine buffer, especially of the diagonal elements of hessian. I would appreciate it if anyone explains this to me.
evaleev
Hi. I want to implement molecular hessian im my code. Since there are usually too many second-order derivatives of electron integrals to save in memory, in my code I computed each integral derivative's contribution to molecular hessian (different from the example
hartree-fock++.ccwhich saves all derivatives in memory). However, the resultant diagonal hessian elements are all incorrect while the off-diagonal ones are correct. The minimal codes are put as follows. They seem long but are actually quite basic for a libint2 user.The results are
while the true answer should be as given by
hartree-fock++.ccInitially, I thought the comments beginning from the 1134th line of
hartree-fock++.ccmight be helpful. It says that the diagonal elements of a second-order derivative (i.e., d2 (s1|s2) / dX2 ) should be treated differently from the off-diagonal ones (i.e., d2 (s1|s2) / dXdY ). However, it appears to me that I don't need to do anything different becauselibint2already gives the final derivative (e.g. d2 (s1|s2) / dX2 ) instead of its components (e.g. (d2 s1 / dX2 | s2), (d s1 / dX | d s2 / dX) and (s1| d2 s2 / dX2 ) ).Thanks for any help!