下面是小卒从计算数学的角度处理的量子化学编程及相关算法方案选择，希望对想要了解并自己动手写量化的同学一点帮助。
所谓计算数学角度是指：建立模型->分析模型解的定性性质->离散方程->数值求解->解释结果->是否需要修正模型，etc
反映这个观点的思路可以见小卒上一个帖子：http://bbs.keinsci.com/forum.php?mod=viewthread&tid=1110&extra=
此帖和Sob老师讨论了Hatree-Fock方程中C的维数的问题，之后在编写自洽场迭代程序时才发现Sob老师的见解是没错的，
但Cances等人从数学物理角度处理的方式也没错，分歧只是C中要不要加入非占据轨道而已，这无关大雅，只要密度矩阵D的构造没错就行。
这段时间小卒正是在这方面缺少很多量化知识，调试走了一些弯路，在此要感谢@Shannon，@卡开发发，两位在这方面给予的帮助。。。
下面开始介绍下程序内容，如下图片展示了求解NH3分子的RHF模型的main程序，小卒的笔记本是：i3处理器，2.10GHz，比较out了。
这套数值方案大致如下，参考的全部文献放在这里：http://yunpan.cn/cwfERQqmupjZe 访问密码 26a8
1）初期单电子积分采用Obara-Saika递推，Boys函数则使用适当精度的Gauss-Legendre数值求积过渡下；
2）双电子积分采用Gill和Pople在90年代提出的HGP-PRISM算法，使用了CCTTT Path；
3）第一次组装双电子积分形成Fock矩阵，之后保存双电子积分需要时提取；
4 ) 广义本征值问题使用平方根分解，也可以Cholesky分解，但重叠矩阵正定性很弱，故前者可能更加；
5 ) 不动点迭代求解器：Roothaan算法，Level Shifting算法，DIIS算法尽量都要有
6）后期Boys函数采用建表插值，使用Taylor插值或者Chebyshev插值，见Gill[1991]
7 ) 用PRISM算法一样的模式取代单电子积分程序的Obara-Saika方法，提高速度的同时，也为解析梯度的计算做准备
8）在组装程序中加入分子对称性信息
下面逐个简单的对上面的内容做个介绍：
1.分子积分的Obara-Saika递推见Helgaker等人著作《Molecular Electronic-Structure Theory》，里面有详细的介绍，
因为CGTO的高收缩度，如果算法满足“early Contraction”效果是很好的。Obara-Saika递推不满足，而且双电子
积分递推是8-term，这不是我所追求的高效。一个简单的分子积分程序可以看看@Shannon兄的matlab代码，
http://bbs.keinsci.com/forum.php?mod=viewthread&tid=933，他使用的是MD算法，《Molecular Electronic-Structure Theory》
书中也有详细介绍，当然MD方法也无法提前收缩，故也不考虑。当你看到HGP-PRISM算法或者MD-PRISM的CCTTT Path时，会发现那效果那是妥妥的。
因为提前收缩，一般显著减少了K^4计算量的系数，文献里面有x代表这个系数值，不妨对比下各算法的数值。
2.学习双电子积分算法PRISM，可以见Gill小组的Paper：https://rsc.anu.edu.au/~pgill/papers.php ，读过后，还需要自动动手推导
一遍，这是必须的。这份PDF文件，是小卒为了学习PRISM算法而推导的资料，并做了部分整理，不妨做个参考，入口是
http://yunpan.cn/cwf2qrcApQKez 访问密码 1881；内容如下：
一、基本定义表示 二、径向势的双电子积分和库伦势的双电子积分 三、Hermite Gaussian函数的一系列结论 四、双电子积分的一个高效递推算法 五、HGP-PRISM算法CCTTT Path的公式推导 六、Boys函数的数值逼近 3.Fock矩阵的组装是一个头痛的事，可以参考下这份资料：http://www.cchem.berkeley.edu/chem195/h2o_8m.html
RHF组装，里面给出的参考书籍是：Szabo和Ostlund的《Modern Quantum Chemistry》，也很有帮助，略微看看关键点就行，不需要细读。
这部分需要量化的一些知识了，Sob老师说太注重数学化，可能难抓住物理本质，这是深刻的教训。
4.关于Cholesky分解，发现重叠矩阵的行列式比较小，对于NH3分子，此次计算的是0.16，条件数则不大，原因是S矩阵的特征值虽然都是正的，
但比较小，谱半径也较小，这样会不会造成Cholesky分解不稳定很难说，尤其对比较大的分子体系。故平方根分解比较妥当一些，缺点是计算量比较大。但也没理由放弃Cholesky分解。
5.最近在看Polay的关于DIIS算法的论文，研读下准备程序化，一般比Roothaan算法快很多。求解器对不同分子体系会表现出不同的收敛性质，
据说DIIS算法对高角动量轨道的HF方程效果差。多准备一些求解器是很要紧的。
6-7.这是这段时间准备弄的方案，为了快速计算能量导数，是需要的。
8.组装程序，加入分子对称性，这个比较难，目前的思路是需要计算原子连接的拓扑信息，之后通过这些信息来给出哪些Shell的双电子积分是一样的。
这个目前程序化方面暂时没啥思路。
以上都是些简单的介绍，这段时间写程序，真的是乏了。慢慢整了，身体吃不消这折腾 = =！。。。。
此程序花了一个月的时间，1500+行代码，用C++OOP的方式写的，可能有人会吐槽这个，其实用什么语言无所谓。重要的是算法和数据结构，
编程过程中，都是遇到问题解决问题，也不可能提前就建立架构的。这是小卒的一个编程记录，可以给大家看看：http://yunpan.cn/cwf8TEyszneWN 访问密码 2786
| 
发表于 Post on 2015-5-22 10:40:58
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
，因为后期要拓展到Monkhorst-Pack网格求和以及RPA，就不得不考虑一些特殊的问题了。
楼主 Author