计算化学公社

标题: Hartree-Fock的C++代码的SCF迭代问题 [打印本页]

作者
Author: 最爱喵星人 时间: 2018-12-13 16:33
标题: 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

作者
Author: 最爱喵星人 时间: 2018-12-13 16:57
//密度矩阵计算
  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;
}
  }

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

作者
Author: Warm_Cloud 时间: 2018-12-14 19:41
检查一下Fock_calc吧，怀疑你没有把波函数投影回去。

作者
Author: 最爱喵星人 时间: 2018-12-14 21:15

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

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

作者
Author: 最爱喵星人 时间: 2018-12-14 21:21
本帖最后由最爱喵星人于 2018-12-14 22:19 编辑

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

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

作者
Author: granvia 时间: 2018-12-14 21:41
平方用pow()函数来算。这编程的例子也太不地道了吧

作者
Author: Warm_Cloud 时间: 2018-12-15 16:36
参考一下这个吧http://bbs.keinsci.com/thread-6139-1-1.html

欢迎光临计算化学公社 (http://bbs.keinsci.com/)