Hartree-Fock的C++代码的SCF迭代问题

最爱喵星人

我想重复一下crawdad 练习里hf能量的计算，但是以下代码不能得到正确的hf结果，能量一直在震荡。请教大家是什么原因呢？

到第一步重建fock矩阵的时候还能重复别人结果，但是迭代第一步的能量就错了。所以应该是迭代第一步密度矩阵哪里出了问题。但是当用Hcore做初猜时fock矩阵和密度矩阵的建立还都是正确的。所以不知道哪里出错了。

谢谢。


//-------------------------------------------------
//form guess fock matrix
//------------------------------------------------
  CMatrix F0,ST,C0;//Fock matrix, transpose of S(-1/2), MO coeff
  F0.NewMat(row,col);
  ST.NewMat(row,col);
  C0.NewMat(row,col);
  Fock_calc(F0,inv_sqrt_overlap,Hcore,C0);  //计算 S^(-1/2)*Hcore*S^(-1/2), 对角化，对角化后的fock矩阵存到F0，本征向量矩阵存到C0

//----------------------
//initial density
//------------------------
  CMatrix D;//density
  D.NewMat(row,col);

  int occ=5;//10 electrons for water, 5 orbitals
  Density_calc(D,C0,occ); //计算密度矩阵， 见最后

//------------------
//initial SCF energy
//------------------
  double Eele=0.0;
  for(int mu=0; mu<row; mu++){
    for(int nu=0; nu<col; nu++){
       Eele+=D[mu][nu]*(Hcore[mu][nu]+Hcore_copy[mu][nu]);
    }
  }

  double Etotal=0.0;
  double ENN=  13.497304462036480;
  Etotal=Eele+ENN;
  msg.Header(outfile,"SCF iterations");
  cout<<"Iter          "<<"E(elec)          "<<"Etot"<<endl;
  cout<<"00  "<<fixed<<setw(20)<<setprecision(15)<<Eele<<"     "<<Etotal<<endl; //初始能量能够重复别人结果


  double deltaE=1e-8;
  double RMSD=1e-8;
  double tmpE, tmpD=0.0;
  int iter=0;

  double Eold=0.0;
  double Ecurrent=0.0;


  while (iter<200)
{
  iter++;
//------------------------------------
//construct new Fock matrix
//------------------------------------
  for(int mu=0; mu<row; mu++){
   for(int nu=0; nu<col; nu++){
     Fock[mu][nu]=Hcore[mu][nu];
     for(int lemda=0; lemda<row; lemda++){
       for(int sigma=0; sigma<col; sigma++){
          Fock[mu][nu]+=Dnew[lemda][sigma]*(2.0*I(mu,nu,lemda,sigma)-I(mu,lemda,nu,sigma)); //迭代第一步还能重复fock矩阵
       }//sigma
      }//lemda
    }//nu
   }//mu

  Fnew_unDiag.CopyMat(Fock);  //将未对角化的fock矩阵保存到Fnew_unDiag
//--------------------
//diagnolize Fock
//new density matrix
//--------------------
  Fnew.Init();
  Cnew.Init();
  Dnew.Init();
  Fock_calc(Fnew,inv_sqrt_overlap,Fock,Cnew);
  Density_calc(Dnew,Cnew,occ);

//---------------------
//new scf energy
//---------------------
Ecurrent=0.0;
  for(int mu=0; mu<row; mu++){
    for(int nu=0; nu<col; nu++){
       Ecurrent+=Dnew[mu][nu]*(Hcore[mu][nu]+Fnew_unDiag[mu][nu]); //能量不正确
    }
  }


//convergence check
  tmpE=Ecurrent-Eold;
  tmpD=0.0;
  for(int mu=0; mu<occ; mu++){
    for(int nu=0; nu<occ; nu++){
       tmpD+=pow( (Dnew[mu][nu]-D[mu][nu]),2.0);
       D[mu][nu]=Dnew[mu][nu];
    }
  }
  tmpD=sqrt(tmpD);

  if(tmpE<deltaE and tmpD<RMSD)
    {
      msg.Header(outfile,"SCF converged");
      break;
    }

  Eold=Ecurrent;
  Etotal=Ecurrent+ENN;
//  cout<<"Iter          "<<"E(elec)          "<<"Etot"<<endl;
  cout<<iter  <<"   "<<Ecurrent<<"     "<<Etotal<<"    "<<tmpE<<"  "<<tmpD<<endl;

}//while

最爱喵星人

//密度矩阵计算
  double tmp;
  for(int i=0; i<C.rows; i++){
    for(int j=0; j<C.cols; j++){
       tmp=0.0;
       for(int m=0; m<occ; m++){  //10 electrons for water, 5 orbitals
         tmp+=C[i][m]*C[j][m];
//          tmp+=C[m][i]*C[m][j];

      }
       D[i][j]=tmp;
    }
  }

sobereva

注意帖子标题应尽可能明确体现帖子内容、不容易引起歧义（见http://bbs.keinsci.com/thread-9348-1-1.html），我已把你的标题“hf cpp迭代问题”改了，望以后注意

Warm_Cloud

检查一下Fock_calc吧，怀疑你没有把波函数投影回去。

最爱喵星人

Warm_Cloud 发表于 2018-12-14 19:41
检查一下Fock_calc吧，怀疑你没有把波函数投影回去。

谢谢。我看你在别的帖子里好像也写过这个。如果方便的话可不可以把每步迭代生成的MO矩阵，密度矩阵和fock矩阵打印出来，我好比较一下。

最爱喵星人

Warm_Cloud 发表于 2018-12-14 19:41
检查一下Fock_calc吧，怀疑你没有把波函数投影回去。

可能我的对角化矩阵有问题。我再看一下。谢谢。

granvia

平方用pow()函数来算。这编程的例子也太不地道了吧

Warm_Cloud

参考一下这个吧http://bbs.keinsci.com/thread-6139-1-1.html