量化编程（三）：HF-DFT算法及代码（Fortran）

nagami

量化编程(三)：HF-DFT算法及代码（Fortran）

文/Nagami   2016-1-7

程序说明：用于学习求解Hartree-Fock方程及Kohn-Sham方程的数值算法
编程语言：Fortran 90
编程环境：Linux/Gfortran，Windows/IVF （皆可）
核心代码量：4000行左右
源程序下载：托管于自由软件社区Github ：https://github.com/xiangyunbin/quantum_scf

一、前言

小编的前几个帖子里，曾用C++写了一个简单的HF框架。后续发现程序存在潜在的BUG，虽然分子积分采用的是MD-PRISM-CCTTT算法，但实现效率很差（用C++的Hash索引），同时，程序也不够美观，影响阅读。之后参考了其它的一些源代码，比如ergoSCF，molgw，etc，经过几次失败的Fortran编程之后，逐渐去粗取精，基本的框架方案大致形成，小编觉得:C++的泛型编程和Fortran的向量化编程，各有优缺点，对与数值编程哪一种更佳，不可定论。


一个基本的对比数据：
笔记本：Core i3，2.10GHz
上一个帖子VS2012的C++代码计算NH3的RHF/STO-3G，DIIS-SCF收敛后的时间为320ms，现在Fortran版本SCF，IVF下跑60ms就可完成。可见效率提高很明显
计算了C6H6的RHF/STO-3G，大约花16s时间，其中15s是在计算双电子积分，Rys积分没有对多中心做过优化，只简单处理了同中心轨道正交的问题。同时注意Rys积分的精度问题

框架方案：
1）对读入的分子信息进行类型划分，按类型读入基组
2）基组通过数组形式存储，便于Fortran的向量化计算
3）对分子积分的递推数据进行拆分归类，哪些只跟原子类型有关，哪些同时与原子类型和坐标位置有关。
4）单电子积分采用Obara-Saika方法，双电子积分采用Rys Quadrature，向量化方式计算分子积分
5）双电子积分按[ij|kl]的对称性进行存储
6）DFT的格点配了Gill的G0方案
7）Becke的单位分解来处理全域积分

进一步可研究的问题：
1）并行算法，及并行化
1）等价类划分处理分子对称性，如何用于双电子积分的计算，及双电子积分的旋转变换规则
2）解析导数：现配备的分子积分方案，对解析导数的计算留了后路
3）HF方程求解的直接极小化方法，见Eric Cances的文章，极为有趣
4）几何优化：带解析导数的拟牛顿方法
从数值角度看，这里的1）和3)最有研读的价值，小编后续有时间，会关注这两点


二、程序模块说明

2-1 m_atoms  :
读入原子信息，做相关的处理，带输出

2-2 m_shells :
读入基组信息，做shell pair数据等相关的处理，带输出

2-3 m_basis :
构建基组信息，及gauss乘积数据，单电子积分程序；带输出

2-4 m_rysquad :
rys数值积分方案：带移位雅克比多项式的修正矩方法处理Rys积分、Stieltjes方法处理非经典权Gauss积分程序

2-4 m_rysroot :
Gamess的Rys积分nodes和weights的1-5阶拟合数据

2-5 m_eri:
Rys双电子积分计算及对称性存储

2-6 m_hamilton：
Fock矩阵组装及广义本征值求解器，采用tred2和tqli的经典程序

2-7 m_dft_grid:
MultiExp的径向程序，Gill的SG0网格,Becke的单位分解

2-8 m_dft_xc:
密度函数计算，在DFTxclib.F中提取对应的泛函


三、算法及文献整理
https://yunpan.cn/c3YXcgvWxbq7T  访问密码 da81




四、小编寄语

DFT部分如大家自己有兴趣，可以再扩展。小编将不着眼于量化的应用，只略涉及量化中的部分算法。

为了大家在不同平台下自由学习，全部以源代码的方式做的，没有添加任何库。

以后的重点仍就是并行算法和基本数值方法的学习，如果大家有好的量化编程的点子欢迎与本人分享。

此外，小编能力泛泛，对量化有理解不到位之处在所难免，大家对程序的正确与否，应该有自己的判断。

2016新年即将来临之际，希望这份资料对大家在量化中的学习有所帮助。


五、程序输入与输出说明
注：
1.目前测试了H2,CH4,NH3,C6H6计算的能量与Gaussian是一致的，但密度矩阵和本征值与经典RHF模型不同，这里的求解模型参考Libint2的SCF，与如下书籍模型一致
P.G. Ciarlet, J.L. Lions, C. Le Bris Handbook of Numerical Analysis. Special Volume_ Computational Chemistry Volume 10 2003
具体不一样的原因，小编也是茫然。如大家知道，可以分享下结论。
2.main_grid.f90是用于测试Becke单位分解对全局积分的精度
PS：最近本想把如上本征值不一致的问题调试下，用之前的代码跑了下结果发现是正确的，匪夷所思。之前上传的代码按理本征值和密度矩阵也是OK的，大家可以试一试，是否如此。

计算甲烷分子CH4的RHF能量，RHF/STO-3G
Main程序

!*
	!* Copyright (C) 2016 XiangYunBin <1084066694@qq.com>
	!*
	program main
	use m_definitions
	use m_atoms
	use m_shells
	use m_basis
	use m_eri
	use m_hamilton
	use m_dft_grid
	use m_dft_xc
	implicit none
	character(len=100)                     ::  basis_path,basis_name
	character(len=100)                     ::  geom_path,filename
	type(geom)                                ::  mol
	type(shell_set)                            ::  shells
	type(basis_set)                           ::  basis
	type(basis_quartet_set)               ::  basisqs
	!==
	integer,parameter                       ::  iter_max = 25
	real(dp),parameter                      ::  threshold = 1.0d-8
	integer                                        ::  nbas,ne,ndocc
	integer                                        :: iter
	real(dp)                                      ::  time0,time1,time2,time3,time4
	real(dp),allocatable                      ::  S(:,:),K(:,:),V(:,:)
	real(dp),allocatable                      ::  H(:,:),F(:,:)
	real(dp),allocatable                      ::  D(:,:),S_minus_sqrt(:,:)
	real(dp)                                      ::  ehf,enuc,ehf_old,err
	logical                                        ::  print_
	call cpu_time(time0)
	filename = "out.txt"
	open(stdout,file = filename,status = "old")
	geom_path = "ch4.txt"
	call readatoms(geom_path,mol,print_)
	basis_path = "/home/xiang/gfort/scf/basis"
	basis_name = "STO-3G"
	call read_shell_set(basis_path,basis_name,mol,shells,print_)
	call init_basis_set(mol,shells,basis,print_)
	nbas = basis%nbf
	allocate( S( nbas, nbas),&
	& K( nbas, nbas),&
	& V( nbas, nbas),&
	& H(nbas,nbas),&
	& F(nbas,nbas),&
	& D(nbas,nbas), &
	& S_minus_sqrt(nbas,nbas)&
	)
	！计算重叠矩阵和动能矩阵，并输出
	call do_overlap_kinetic_matrix(mol,shells,basis,S,K,print_) ！计算核吸引能矩阵，并输出
	call do_nucleus_matrix(mol,shells,basis,V,print_)
	call cpu_time(time1)
	！计算双电子积分，并按对称性输出，可观察0值很多，轨道正交导致
	call precalc_eri(shells,basis,basisqs,print_)
	call cpu_time(time2)
	call nucleus_nucleus_energy(mol,enuc)
	call seutp_S_minus_sqrt(S,S_minus_sqrt,print_)
	ne = sum(mol%zatom(:))
	ndocc = ne/2
	H = K + V
	call generalized_eigen_solver(H,S_minus_sqrt,D,ndocc)
	ehf = sum( D * (H + H) )
	write(stdout,'(10x,a,17x,a)'),"Energy", "Error"
	！简单SCF自洽迭代
	do iter=1,iter_max
	ehf_old = ehf
	call setup_fock_matirx(D,basis,basisqs,F)
	F = F + H
	call generalized_eigen_solver(F,S_minus_sqrt,D,ndocc)
	ehf = sum( D * (H + F) )
	err = (ehf - ehf_old)/ehf
	write(stdout,'(2f12.6)'),ehf + enuc,err
	if( abs(err) < threshold ) exit
	if( iter == iter_max-1 ) stop 'SCF : Not Convergence! '
	end do
	call cpu_time(time3)
	write(stdout,'(/)')
	write(stdout,'(a,2x,f12.6)'),"1e Hamiltonian",(time1-time0)
	write(stdout,'(a,2x,f12.6)'),"Pre ERI :",(time2-time1)
	write(stdout,'(a,2x,f12.6)'),"SCF :",(time3-time2)
	write(stdout,'(a,2x,f12.6)'),"Total Time :",(time3-time0)
	call destroy_geom(mol)
	call destroy_shell_set(shells)
	call destroy_basis_set(basis)
	call destroy_basis_quartet_set(basisqs)
	deallocate(S,K,V,H,F,D,S_minus_sqrt)
	end program

输出结果文件，out.txt：

=========== Atom List ============

	NAtoms :   5
	NAtomic Type :   2
	Element  Number
	H        4
	C        1
	Input orientation:
	Element  TypeIndex    Coordination
	-------------------------------------------
	C          2        0.000000    0.000000    0.000000
	H          1        0.000000    0.000000    2.057912
	H          1        1.940218    0.000000   -0.685971
	H          1       -0.970110   -1.680278   -0.685971
	H          1       -0.970110    1.680278   -0.685971
	=========== Shell List ============
	Number of shells :   3
	Basis set : STO-3G
	The information of shells below :
	ishell=  1  mom=  0  ishell_type=  1  bool_sp=F
	alpha=      coeff=
	0.168855    0.083474
	0.623914    0.267839
	3.425251    0.276934
	ishell=  2  mom=  0  ishell_type=  2  bool_sp=F
	alpha=      coeff=
	3.530512    0.816191
	13.045096    2.618880
	71.616837    2.707814
	ishell=  3  mom=  1  ishell_type=  2  bool_sp=T
	alpha=      coeff=
	0.222290    0.161536    0.085276
	0.683483    0.214036    0.538304
	2.941249   -0.160017    0.856044
	=========== Basis List ============
	Total number of basis functions :   9
	The information of basis functions below :
	1   iatom=  1   ishell=  2   ibas=  1     am=  0
	bool_sp=F   ibf_left=  1   ibf_right=  1
	2   iatom=  1   ishell=  3   ibas=  2     am=  1
	bool_sp=T   ibf_left=  2   ibf_right=  5
	3   iatom=  2   ishell=  1   ibas=  3     am=  0
	bool_sp=F   ibf_left=  6   ibf_right=  6
	4   iatom=  3   ishell=  1   ibas=  4     am=  0
	bool_sp=F   ibf_left=  7   ibf_right=  7
	5   iatom=  4   ishell=  1   ibas=  5     am=  0
	bool_sp=F   ibf_left=  8   ibf_right=  8
	6   iatom=  5   ishell=  1   ibas=  6     am=  0
	bool_sp=F   ibf_left=  9   ibf_right=  9
	=========== Overlap Matrix ============
	1        1.000000    0.248362    0.000000    0.000000    0.000000    0.062421    0.062421    0.062421    0.062421
	2        0.248362    1.000000    0.000000    0.000000    0.000000    0.491330    0.491330    0.491330    0.491330
	3        0.000000    0.000000    1.000000    0.000000    0.000000    0.000000    0.440874   -0.220437   -0.220437
	4        0.000000    0.000000    0.000000    1.000000    0.000000    0.000000    0.000000   -0.381808    0.381808
	5        0.000000    0.000000    0.000000    0.000000    1.000000    0.467618   -0.155873   -0.155873   -0.155873
	6        0.062421    0.491330    0.000000    0.000000    0.467618    1.000000    0.169697    0.169697    0.169697
	7        0.062421    0.491330    0.440874    0.000000   -0.155873    0.169697    1.000000    0.169697    0.169697
	8        0.062421    0.491330   -0.220437   -0.381808   -0.155873    0.169697    0.169697    1.000000    0.169697
	9        0.062421    0.491330   -0.220437    0.381808   -0.155873    0.169697    0.169697    0.169697    1.000000
	=========== Kinitic Matrix ============
	1       15.891122   -0.085890    0.000000    0.000000    0.000000   -0.010490   -0.010490   -0.010490   -0.010490
	2       -0.085890    0.472250    0.000000    0.000000    0.000000    0.107188    0.107188    0.107188    0.107188
	3        0.000000    0.000000    1.477728    0.000000    0.000000    0.000000    0.246818   -0.123409   -0.123409
	4        0.000000    0.000000    0.000000    1.477728    0.000000    0.000000    0.000000   -0.213750    0.213750
	5        0.000000    0.000000    0.000000    0.000000    1.477728    0.261790   -0.087263   -0.087263   -0.087263
	6       -0.010490    0.107188    0.000000    0.000000    0.261790    0.760032   -0.005517   -0.005517   -0.005517
	7       -0.010490    0.107188    0.246818    0.000000   -0.087263   -0.005517    0.760032   -0.005517   -0.005517
	8       -0.010490    0.107188   -0.123409   -0.213750   -0.087263   -0.005517   -0.005517    0.760032   -0.005517
	9       -0.010490    0.107188   -0.123409    0.213750   -0.087263   -0.005517   -0.005517   -0.005517    0.760032
	=========== Nucleus Matrix ============
	1      -35.596459   -4.674975   -0.000000    0.000000   -0.000000   -1.195826   -1.195826   -1.195825   -1.195825
	2       -4.674975   -7.078563   -0.000000    0.000000   -0.000000   -2.961970   -2.961970   -2.961968   -2.961968
	3       -0.000000   -0.000000   -7.036713    0.000000    0.000000    0.000000   -2.456383    1.228192    1.228192
	4        0.000000    0.000000    0.000000   -7.036712    0.000000    0.000000    0.000000    2.127289   -2.127289
	5       -0.000000   -0.000000    0.000000    0.000000   -7.036713   -2.605388    0.868463    0.868462    0.868462
	6       -1.195826   -2.961970    0.000000    0.000000   -2.605388   -4.970755   -0.922493   -0.922492   -0.922492
	7       -1.195826   -2.961970   -2.456383    0.000000    0.868463   -0.922493   -4.970755   -0.922492   -0.922492
	8       -1.195825   -2.961968    1.228192    2.127289    0.868462   -0.922492   -0.922492   -4.970754   -0.922492
	9       -1.195825   -2.961968    1.228192   -2.127289    0.868462   -0.922492   -0.922492   -0.922492   -4.970754
	=========== ERI DATA ============
	1      1      1      1        3.541948
	2      1      1      1        0.576921    0.000000    0.000000    0.000000
	2      1      2      1        0.111888    0.000000    0.000000    0.000000    0.000000    0.020296    0.000000    0.000000
	0.000000    0.000000    0.020296    0.000000    0.000000    0.000000    0.000000    0.020296
	2      2      1      1        0.854682    0.000000    0.000000    0.000000    0.000000    0.852436    0.000000    0.000000
	0.000000    0.000000    0.852436    0.000000    0.000000    0.000000    0.000000    0.852436
	... ...
	6      6      6      5        0.084274
	6      6      6      6        0.774606
	=========== S^(-1/2) ============
	1        1.029974   -0.180891   -0.000000    0.000000   -0.000000    0.026095    0.026095    0.026095    0.026095
	2       -0.180891    1.582802    0.000000    0.000000    0.000000   -0.344500   -0.344500   -0.344500   -0.344500
	3       -0.000000    0.000000    1.173497    0.000000   -0.000000    0.000000   -0.330733    0.165366    0.165366
	4        0.000000    0.000000    0.000000    1.173497    0.000000   -0.000000    0.000000    0.286423   -0.286423
	5       -0.000000    0.000000   -0.000000    0.000000    1.173497   -0.350795    0.116932    0.116932    0.116932
	6        0.026095   -0.344500    0.000000   -0.000000   -0.350795    1.273019   -0.027781   -0.027781   -0.027781
	7        0.026095   -0.344500   -0.330733    0.000000    0.116932   -0.027781    1.273019   -0.027781   -0.027781
	8        0.026095   -0.344500    0.165366    0.286423    0.116932   -0.027781   -0.027781    1.273018   -0.027781
	9        0.026095   -0.344500    0.165366   -0.286423    0.116932   -0.027781   -0.027781   -0.027781    1.273018
	Energy        Error
	-34.121400   -0.790624
	-40.791844    0.122981
	-39.528968   -0.023838
	-39.762067    0.004381
	-39.720417   -0.000783
	-39.727881    0.000140
	-39.726540   -0.000025
	-39.726780    0.000005
	-39.726737   -0.000001
	-39.726744    0.000000
	-39.726743   -0.000000
	-39.726743    0.000000
	1e Hamiltonian      0.012000
	Pre ERI :      0.048003
	SCF :      0.004000
	Total Time :      0.064003

苯分子C6H6out.txt文件数据：（关掉Hamilton矩阵输出，关掉ERI数据输出）

=========== Atom List ============
NAtoms :  12
NAtomic Type :   2

Element  Number
   H        6
   C        6

Input orientation:
Element  TypeIndex    Coordination
-------------------------------------------
  C          2       -0.851586    1.662638    0.000000
  C          2        1.784884    1.662638    0.000000
  C          2        3.103025    3.944954    0.000000
  C          2        1.784657    6.228707   -0.002268
  C          2       -0.851170    6.228556   -0.003175
  C          2       -2.169462    3.945370   -0.001285
  H          1       -1.890482   -0.136986    0.000850
  H          1        2.823289   -0.137345    0.002494
  H          1        5.181119    3.945106    0.001191
  H          1        2.824385    8.027991   -0.002381
  H          1       -1.890765    8.028105   -0.004970
  H          1       -4.247405    3.945729   -0.001625



=========== Shell List ============

Number of shells :   3

Basis set : STO-3G                                                                                             

The information of shells below :
ishell=  1  mom=  0  ishell_type=  1  bool_sp= F
  alpha=      coeff=
     0.168855    0.083474
     0.623914    0.267839
     3.425251    0.276934


ishell=  2  mom=  0  ishell_type=  2  bool_sp= F
  alpha=      coeff=
     3.530512    0.816191
    13.045096    2.618880
    71.616837    2.707814


ishell=  3  mom=  1  ishell_type=  2  bool_sp= T
  alpha=      coeff=
     0.222290    0.161536    0.085276
     0.683483    0.214036    0.538304
     2.941249   -0.160017    0.856044


          Energy                 Error
  1    -192.49084140           -0.87523432
  2    -239.50781075            0.10624375
  3    -226.23873500           -0.03091082
  4    -227.99979991            0.00408571
  5    -227.88083852           -0.00027607
  6    -227.88963455            0.00002041
  7    -227.88932904           -0.00000071
  8    -227.88939969            0.00000016
  9    -227.88940592            0.00000001
10    -227.88940771            0.00000000


1e Hamiltonian      0.234002
Pre ERI :     15.163297
SCF :      0.390002
Total Time :     15.787301

brothers

犀利....大神能说一下这个写了多久吗？我估算一下我们之间到底差了几个等级。。。。

nagami

brothers 发表于 2016-1-8 20:41
犀利....大神能说一下这个写了多久吗？我估算一下我们之间到底差了几个等级。。。。

过奖了，断断续续，2个月应该有的

nagami

brothers 发表于 2016-1-8 20:41
犀利....大神能说一下这个写了多久吗？我估算一下我们之间到底差了几个等级。。。。

最早接触量化编程是这个时候：http://emuch.net/bbs/viewthread.php?tid=8000420

brothers

nagami 发表于 2016-1-9 22:26
最早接触量化编程是这个时候：http://emuch.net/bbs/viewthread.php?tid=8000420

学的好快~佩服

caofan_success

nagami 发表于 2016-1-9 22:26
最早接触量化编程是这个时候：http://emuch.net/bbs/viewthread.php?tid=8000420

佩服！

西乡新丰客

您好，云盘链接失效，还有前面两个贴子的云盘链接也都失效了，能否再传下呢？谢谢

Jack

360云盘以关闭，希望大神能再上传下资料，非常感谢！

wwjie

您好，云盘链接失效，资料能再上传一遍吗
正在学习量化计算的程序实现，相信您的资料会对我有很大帮助

清微

请问一下，1.在您的代码里写的是geom_path = "ch4.txt"，再算其它分子(如c6h6)时,需要在main.f90里面修改吧？2.在windows下运行，必须先有一个out文件，否则会出错：Fortran runtime error: Cannot open file 'out.txt': No such file or directory，请问这是为什么？

清微

请问一下，1.在您的代码里写的是geom_path = "ch4.txt"，再算其它分子(如c6h6)时,需要在main.f90里面修改吧？2.在windows下运行，必须在当前目录下先有一个out文件，否则会出错：Fortran runtime error: Cannot open file 'out.txt': No such file or directory，请问这是为什么？