计算化学公社

标题: 量化编程（三）：HF-DFT算法及代码（Fortran） [打印本页]

作者
Author: nagami 时间: 2016-1-7 11:59
标题: 量化编程（三）：HF-DFT算法及代码（Fortran）
本帖最后由 nagami 于 2016-1-23 15:00 编辑

量化编程(三)：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

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

作者
Author: nagami 时间: 2016-1-9 18:10

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

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

作者
Author: nagami 时间: 2016-1-9 22:26

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

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

作者
Author: brothers 时间: 2016-1-10 18:34

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

学的好快~佩服

作者
Author: caofan_success 时间: 2017-1-13 10:23

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

佩服！

作者
Author: 西乡新丰客 时间: 2017-1-13 20:01
您好，云盘链接失效，还有前面两个贴子的云盘链接也都失效了，能否再传下呢？谢谢

作者
Author: Jack 时间: 2017-4-16 20:53
360云盘以关闭，希望大神能再上传下资料，非常感谢！

作者
Author: wwjie 时间: 2018-1-6 01:51
您好，云盘链接失效，资料能再上传一遍吗

正在学习量化计算的程序实现，相信您的资料会对我有很大帮助

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

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

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