量化编程（二）:代码优化

nagami · 发表于 Post on 2015-6-5 11:07:58

本帖最后由 nagami 于 2015-12-15 12:22 编辑

继续上一帖几个未完成的方案：
   （一）将双电子积分的PRISM-CCTTT算法移植到单电子积分，提前收缩计算
   （二）Boys函数的低浮点数高精度单点计算
   （三）DIIS外推加速技术
一些难度更高的未实现的问题：
   （四）分子对称性的考虑
   （五）特征值求解的迭代算法
下面逐个简单总结性讨论下，并附上计算效果图：
（一）单电子积分的PRISM-CCTTT算法
      这个方案以Gill的文章《An Efficient Algorithmfor the Generation of Two-Electron Repulsion Integrals》为依据展开话题，
小卒的双电子积分就是程序化这篇文章的算法，关于单电子积分的PRISM-CCTTT算法的公式目前都已推导过，并实现了程序化，
效果比原先提升了一个档次，这也是显而易见的，而且结构也更加清晰。可以说，重叠积分、动量积分、动能积分、偶极矩积分、
核吸引能积分、双电子排斥积分，都可以通过同一种模式完成高效递推计算。对于Pople型基组，效果还是不错的。对于月份基组，
大量的非收缩高角动量基函数未必有效，但可以让程序进行智能选择，不一定取CCTTT路径，改用其它O(1)更小的路径，易或者
采用Rys数值积分，通过这种混合方式进行处理也不错。
      一个数据对比：上一帖的单电子积分采用Obara-Saika递推，NH3分子的HF/STO-3G水平计算单电子积分，所花费的时间为60ms，
而新算法则只需10ms。
   对该算法做一个简短介绍，首先递推公式见附件图片，分子积分计算的基本任务非常明确，快速高精度计算(a|O|b),(ab|O|cd),
这里O是算符。常见的算符O有两种形式，一种是函数式算符，比如重叠积分，偶极矩积分，核吸引能积分、双电子排斥积分，etc；
另一种是微分算符，比如动量积分、动能积分，etc。不同的算符类型，HRR公式将会不同，同时注意到HRR公式使得基函数a，b，
c，d的地位是非对称的，下面小卒会再回到这个话题，去讨论该算法的一些加速小技巧。
   算法的基本结构呈如下形式，无论单电子积分还是双电子积分都是如此：
   1.HRR递推公式
   2.HGTO递推公式
   3.[r]递推公式
   4.初值计算公式
   在推导公式的时候，你会发现计算了重叠积分，等于完成了动量积分、动能积分、偶极矩积分的初值计算公式，和[r]递推公式，
因为不需要额外再附加计算。HGTO递推公式是通过Hermite Gaussian作为辅助计算而引入的，首先是算法本身计算Hermite Gaussian
非常高效，同时CartesianGaussian可以通过对Hermite Gaussian降阶来获得，当然千方百计修改算法，还有其它方面的妙用。
   下面考虑计算加速的小技巧：对角动量做排序。为何会加速呢？下面以双电子积分(ab|cd)阐述
首先注意到算法的流程：(0)->(m0|00)->(m0|n0)->(ab|n0)->(ab|cd),计算路径本身是固化的，不会改变
我们具体以(ps|ss),(sp|ss),(ss|ps),(ss|sp)为例说明，根据分子积分的对称性，它们都是一样的值,(p:p型gaussian基函数，s:s型gaussian型基函数)
(ps|ss)算法的计算流程变为：(0)->(10|00)->完成输出
(sp|ss)算法的计算流程变为：(0)->(10|00)->(01|00)->完成输出
(ss|ps)算法的计算流程变为：(0)->空白调用->(00|10)->完成输出
(ss|sp)算法的计算流程变为：(0)->空白调用->(00|10)->(00|01)->完成输出
从上可见，同样的积分值(ps|ss)是最快的，为此必须在分子对称性的基础上把高角动量往前面排，就是如下的做法排序，以(12|34)为例,l为角动量
Shell[1].l > Shell[2].l
Shell[3].l > Shell[4].l
Shell[1].l + Shell[2].l  > Shell[3].l + Shell[4].l
这种略微的改变会让你的算法变得快速很多，单电子积分类似的做法。
（二）Boys函数计算
在我们讨论组对此有过详细的讨论，@卡开发发 ,Boys函数计算目前有下面几种公式：
1）Taylor级数展开（交错级数）
2）正项级数展开（含负指数函数）
3）渐近式展开（含根号）
4）连分式展开（含负指数函数）
5）修正有理分式展开（含根号）
6）数值积分
上面几种算法的组合构成了Boys函数计算的基本方法，括号里面的函数是影响速度的关键，首先负指数函数比根号运算慢2倍，依据小卒的编译器计算。
下面谈谈上面几种展开的个人评价：
1）优点：对任意阶有效，不包含指数函数和根号运算，缺点：交错收敛很慢，只适合于x特别小的情况，难以应用
2）优点：收敛效果较好，对x的取值可以适中，一般小值部分用此法。缺点：包含负指数函数计算，级数本身固有的一些缺陷。
3）优点：公式是幂型，适合多极子展开，高效收缩，大x很适合。缺点：不同的阶，截断距离变化大，含根号运算，修正项还含指数运算
4）优点：收敛性较好，对任意阶有效，精度可变。缺点：含负指数运算。
5）优点：低浮点运算，单精度以上分子和分母只需4-5阶多项式运算。缺点：含根号，不同的阶，需要不同的系数配置，没通用公式，精度无法变化。
6）优点：程序化方便，精度以不同数值积分方法调整。缺点：运算量高
目前，小卒的组合是3）+5），采用此种Boys函数计算算法和新的单电子积分算法，单电子积分计算比原先快了2倍。可以看图片数据
（三）DIIS外推加速
以下面的main程序为例子，很清晰的表述

int main(){
clock_t start,endt;
string file = "nh3.txt";
auto atoms = read_dotxyz(file);//读入几何信息
string file2 = "sto-3g.txt";
auto shells = read_basis(atoms,file2);//读入与原子配对的基组信息
cout<<"Molecular : NH3"<<endl;
cout<<"Level/Basis : HF/STO-3G"<<endl;
start = clock();
Molecular<Shell,Atom> mole(shells,atoms);//初始化分子积分程序，构建RHF模型
auto& S = mole.Overlap_Matrix();//重叠矩阵
auto& H = mole.Core_Hamilton_Matrix(); //Core Hamilton矩阵
SCF_Roothaan scf(S,mole.ndocc); //Roothaan不动点迭代求解
auto D = scf.SolveDensityMatrix(H); //由Core Hamilton矩阵初猜密度矩阵
auto n = mole.nbasis;
const auto maxiter = 20;//最大迭代次数
auto iter = 0;//迭代次数
auto rmsd = 0.0;//密度矩阵的Frobenius范数误差
auto ediff = 0.0;//能量误差
auto ehf = 0.0;//电子总能量
auto enuc = mole.Nuclear_Repulsion_Energy();//分子核间排斥能
auto m = 7;//DIIS-m算法
MatDoub D_err(n,n); //相邻密度误差
vector<MatDoub> vector_Derr(m);
vector<MatDoub> vector_D(m);
for(auto i = 0;i < m;i++){ //分配内存
vector_Derr[i].resize(n,n);
vector_D[i].resize(n,n);
}
cout<<fixed<<setprecision(13);
for(auto k = 0;k < maxiter;k++){
auto ehf_last = ehf;
auto D_last = D;
auto F = mole.Fock_Matrix(D);//由密度矩阵构建Fock矩阵
F += H;
D = scf.SolveDensityMatrix(F); //求解HF方程的广义本征值问题，并更新密度矩阵
ehf = 0.0;
for(auto i = 0; i < n;i++) //RHF能量计算
for(auto j = 0; j < n;j++)
ehf += D[i][j] * (H[i][j] + F[i][j]);
ediff = ehf - ehf_last; //能量逼近误差
D_err = D - D_last;
rmsd = sqrt(scf.trace(D_err,D_err)); //密度矩阵逼近误差
vector_D[k % m] = D_last; //保存最新前m序列
vector_Derr[k % m] = D_err; //保存外推矩阵，简单的误差D_err外推
if( k > 0 && k % m == 0 ){ //更新数据之后，开始DIIS递推
MatDoub A(m + 1,m + 1,-1.0);
A[m][m] = 0.0;
VecDoub b(m + 1,0.0);
b[m] = -1.0;
for(auto i = 0 ;i < m;i++)
for(auto j = 0;j < i + 1;j++){
A[i][j] = scf.trace(vector_Derr[i],vector_Derr[j]); //一种矩阵的内积
if(i != j) A[j][i] = A[i][j];
}
LUdcmp lu(A);
lu.solve(b,b);//求解组合系数
for(auto i = 0; i < n;i++)
for(auto j = 0; j < n;j++){
D[i][j] = 0.0;
for(auto l = 0; l < m ;l++)//获得外推后的密度矩阵
D[i][j] += b[l] * vector_D[l][i][j];
}
}
endt = clock();
if (k== 0)
cout<<"\n\nIter E(tot) Delta(E) RMS(D) Time(s)\n";
cout<<setw(3)<<k + 1<<setw(18)<<ehf + enuc<<setw(18)<<ediff<<setw(18)<<rmsd<<setw(18)<<(endt-start)/1000.0<<endl;
}
cin.get();
return 0;
}

复制代码

   上面的程序，就是在原来的基础上多加了一个DIIS的简单外推，不仅对计算时间没什么影响，计算效率还提高了一个档次，可以看下图片，对比下效果，
更新数据的时候都是线性收敛，外推之后，出现超线性收敛，加速极其明显。
关于DIIS，可以见Pulay[1980]的原始文章，起初作者是对几何优化问题考虑组合外推，DIIS是一种通用的算法思想，实际计算中并不是采用外推D_err,
而是[F(D),D]，见Pulay[1982],一种标准的收敛判断模式
（四）分子对称性的考虑
      对于原子数超过10的分子，考虑对称性是必须必的事情，起初的想法是通过将基函数映射到各自的原子，依据原子的对称排布，来确定是否属于同一
等价类，同一等价类的积分值相等对于高对称性分子，比如苯，计算效果可以加速几个档次，但程序化的困难在于如何确定这个等价类，如何对原子的拓扑连接关系进行分类？对此，不知大家有何高见？
（五）特征值求解的迭代算法
      大家已经发现，分子积分的积分速度相对提高了很多，但求解时间依然很长，每轮特征值求解所花费的时间平均是10-20ms，对于NH3的HF/STO-3G，
实际上RHF，并不需要非占据轨道的计算，如果只需迭代求解出最低的几个所需的轨道，那效果就会加速2倍，直观上如此。
davidson【1975】的文章讲述了这件事情，想要高速这是一个重点，也不是几日就能出效果，有时间再慢慢琢磨吧
以上是这段时间的一些简单回顾，希望对大家有所帮助

平辉正 · 发表于 Post on 2015-6-5 20:16:20

正要准备开始学写程序！很有参考价值

zyj19831206 · 发表于 Post on 2015-10-14 09:44:10

楼主，你是怎么学习量化编程的，能够给大家写一写你的学习过程吗？

nagami · 发表于 Post on 2015-10-15 14:56:08

zyj19831206 发表于 2015-10-14 09:44
楼主，你是怎么学习量化编程的，能够给大家写一写你的学习过程吗？

http://www.bccg.org:8888/wordpress/
这是讨论组的博客，或许有一些帮助

lsbp010 · 发表于 Post on 2015-11-12 09:41:41

好帖！学习学习！

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[算法与编程] 量化编程（二）:代码优化

评分 Rate

评分 Rate

浏览过的版块