matlab/python/swig/ libint 编写密度泛函计算代码

Shannon

分享 一个简易的Matlab编写的密度泛函代码，
原理来自 sob的博客： http://sobereva.com/69
原来的HF计算在 http://bbs.keinsci.com/forum.php?mod=viewthread&tid=63&extra= 这个帖子里。 太长了，于是另开新帖。
目前用python/swig调用 libint的双电子积分器已经完工。具体使用方法和原理在6楼。



大量参考了 multiwfn中的 积分代码，因为fortran和matlab语法不同，有很多改动，另外删掉了 不少格点调整的代码，以加大代码的可读性。
注释是用中文写的， 希望 糟糕的语文水平不影响大家阅读。因为积分格点基本没调整的关系，精度并不太好，DFT的积分精度会严重引想自洽场的收敛，因此SCF能量和高斯的结果有 6-8kJ/MoL的差距。
hatree-fock代码结果和高斯是一致的。
积分格点 如图所示（例子是水分子）：


代码可以直接读取高斯的输入文件，具体可参看 代码中提供的例子。

代码可以选择两种 基组 STO2G和3-21G。 3-21G只包含碳氢氧原子的基组信息，STO2G另外包含有氟原子的基组信息。各位可以往里面添加，基组位置参看文件名。

目前DFT 只支持Xalpha 这种 极其原始的 LDA方法，其他的 方面诸如BLYP，B3LYP之类的在后续版本逐渐添加。 这个难度很小，诸位自己写估计也就几十分钟了事。

superrice

THX!

北纬18°

还是提高精度吧

helpme

佩服！

夏梓青华

非常感谢分享   

Shannon

一直想把代码的双电子积分用libint计算，速度估计快好多。 但苦于C++和matlab在windows平台上几乎完全没法兼容。 cygwin上编译了无数次，每次matlab运行 mex文件 都会直接挂掉。
个人是坚定不移的拿来煮义者，如果有代码能利用是绝对不想多写一行的，  用matlab一个好处就是matlab在超高开发效率的同时 ， 能使用C/C++写的mex从而达到极快的计算速率， 同时包含JVM虚拟机，直接调用各种java代码混编 毫无障碍。 可惜， libint 上matlab实在是兼容不了， 线程不安全这个问题 非常 残念，动辄退出令编程极其不爽。
不得已，转战python平台，  现已路转粉。
python那简洁的编程风格 能治愈 业余程序猿 们 习惯乱缩进的顽疾，同时有海量的免费 代码库 能直接用。scipy/numpy 的函数 几乎能逐句翻译matlab代码。

目前代码版本的双电子积分是 用 python 和swig 包装 c++写的 libint库实现的。理论上价值不大，但可能是一个还不错的python参考代码，故分享共赏。
以后有时间再在这个python 代码基础上写出几何导数计算，libint已经自带有几何导数功能，只需要整理一下就能用。
（代码license为GPL，欢迎使用）
7z压缩包内，libint和c++代码已经在windows7 64bit/core i7 3770/cygwin 64bit/ swig 3.05/ python 2.7 32bit/mingw32 下编译好了，直接测试即可。

代码介绍如下：
1. 编译环境
我的环境是 windows7 64bit/core i7 3770/cygwin 64bit/ swig 3.05/ python 2.7 32bit/mingw32
需要用到swig 和 python 2.7 win32 版。 代码用到了numpy 和scipy，这两个没有官方的64bit安装包（非官方安装包有）。
编译需要用mingw的gcc 编译. cygwin的gcc 编译会报错。 libint同样需要mingw编译，但是需要在cygwin环境中。
cygwin中若有gcc，和mingw会冲突，清将cygwin内的gcc和g++改名掉（或者删掉），如果有同学在cygwin内装了python的话，也得改名或者删掉，以免出现冲突。
mingw的gcc和cygwin的gcc生成的库文件有很大差别， 没法兼容。

2. 编译libint
libint的编译是老样子：
./configure
make
make intall
非常简单，就是需要点耐心。
生产的库文件找到，用来第三部工作。
3. 编译python的c模块
3.1，准备接口文件 maincode.i，如下，直接调用头文件就可以了，以后也不需要改。另外，为了能让python使用指针和c++数组，需要cpointer.i和carray.i。这两个swig库文件用法可以参考 http://www.swig.org/Doc1.3/Library.html

1. /* File : interface.i */
   
2. %module maincode
   
3. %{
   
4. #include "maincode.h"
   
5. %}
   
6. %include "maincode.h"
   
7. %include "cpointer.i"
   
8. %include "carrays.i"
   
9. %array_class(double, doubleArray);
   
10. %pointer_class(double, doublep);

复制代码

3.2 准备头文件编译的代码类似函数库，不需要main函数，吧main改成别的名字，然后写个头文件，例如这样的：

1. int compute(double *Result,double * contrcoef1,double * alphascale1,
   
2.         int contrdepth1input,double * A,double * contrcoef2,double * alphascale2,
   
3.         int contrdepth2input,double * B,double * contrcoef3,double * alphascale3,
   
4.         int contrdepth3input,double * C,double * contrcoef4,double * alphascale4,
   
5.         int contrdepth4input,double * D,int am0i,int am1i,int am2i,int am3i);
   
6.         

复制代码

有了swig之后，python可以调用 含有数组和指针的c++函数的。
头文件写好之后，用swig生成 wrap 的c++代码文件：
swip -c++ -python XXXXXXXX.i
3.3 写setup.py，用distutils安装
如题， 直接用gcc编译非常麻烦，推荐用python自带的distutils 安装。
setup.py

1. from distutils.core import setup, Extension
   
2. 
   
3. module1 = Extension(<font color="#ff0000">'_maincode'</font>, sources=['maincode_wrap.cxx',
   
4. 
   
5. 'maincode.cpp'], include_dirs=['libint/include/libint2'],library_dirs=['libint/lib'],libraries=['int2'])
   
6. 
   
7. setup (name = 'maincode',
   
8. 
   
9. version = '1.0',
   
10. 
    
11. description = 'Simple maincode from SWIG tutorial',
    
12. 
    
13. ext_modules = [module1])
    

复制代码

maincode前面那个_符号很重要，没写会导致很难发现的错误。 另外，使用cygwin 的gcc在此处会爆出 NoneType" object is not iterable 这样的错误。 安装：
python setup.py build_ext --compiler=mingw32 --inplace
inplace 的用处是将编译好的pyd文件放在当前文件夹下。

4. matlab调用python
matlab 也可以调用python，直接import py.package名.* 就可以用了。 和调用java类似。
但是自己编译的c模块，用matlab有时会载入失败。 这个代码中是用 .mat 文件在python部分和matlab部分之间传递变量的。python 的scipy可以直接读取.mat文件，然后locals().update(XX)将读取的变量加入python的namespace里面。

ximi1986

赞一个，学习学习，一直也想用matlab编写个~~~不过C++和matlab在windows上无法兼容？Matlab R2014a 不是有code 功能直接生成C程序把？？或者 从C中 调用Matlab引擎？另外，我记得哪本书说没有大循环的时候少用mex,因为matlab 内建函数效率是比较高~~~~还有 神马 并行、显卡计算都用上会快点~~

Shannon

ximi1986 发表于 2015-7-1 15:43
赞一个，学习学习，一直也想用matlab编写个~~~不过C++和matlab在windows上无法兼容？Matlab R2014a 不是有c ...

visual studio好像可以编译mex，这个不知道能否兼容。 用 cygwin 或者mingw的gcc的话，我试了很久，都是没错一运行mex就退出。 但是用matlab coder， 也就是那个生成C程序的东东，生成的c程序gcc编译可以通过，完全不会有任何问题。 估计问题出现在libint代码的一些原因。
量化程序内全是超级大的多重循环…for循环效率这个问题，不仅matlab很慢，python也很慢。DFT积分估计可以移到显卡上，架构上非常适合这个工作，速度大约会飞快（单精度是个问题，而且我手头也没NV的显卡可以玩）。

ximi1986

Visual Studio 2012就可以,我在用，还是比较稳定。另外要装Fortran 2013和 vs 2012配套

最爱喵星人

请问怎么调用libint2的计算函数呢？
我看手册说通过Libint_eri_t 的结构体和libint2_build_eri[n][n][n][n]函数来计算双电子积分。请问怎么把基组信息喂给libint2呢？

Zork

代码就是要够乱，看的才能提神。

granvia

赞lz  matlab最大缺点就是不免费

413

大佬，编译libint-2.0.3-stable版本成功后为啥生成的是libint2.h的头文件，而非libint.h？？？

CCCDDD

你好，是不是少了一个方程的文件