计算化学公社

标题: matlab/python/swig/ libint 编写密度泛函计算代码 [打印本页]

作者
Author: Shannon 时间: 2015-4-11 01:29
标题: matlab/python/swig/ libint 编写密度泛函计算代码
本帖最后由 Shannon 于 2015-7-2 00:48 编辑

分享一个简易的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代码结果和高斯是一致的。
积分格点如图所示（例子是水分子）：
(, 下载次数 Times of downloads: 109)

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

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

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

作者
Author: superrice 时间: 2015-4-11 01:34
THX!

作者
Author: 北纬18° 时间: 2015-4-11 06:19
还是提高精度吧

作者
Author: helpme 时间: 2015-4-12 14:49
佩服！

作者
Author: 夏梓青华 时间: 2015-5-10 13:56
非常感谢分享

作者
Author: Shannon 时间: 2015-7-1 11:24
本帖最后由 Shannon 于 2015-7-2 00:48 编辑

一直想把代码的双电子积分用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

/* File : interface.i */
%module maincode
%{
#include "maincode.h"
%}
%include "maincode.h"
%include "cpointer.i"
%include "carrays.i"
%array_class(double, doubleArray);
%pointer_class(double, doublep);

复制代码

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

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

复制代码

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

from distutils.core import setup, Extension
module1 = Extension(<font color="#ff0000">'_maincode'</font>, sources=['maincode_wrap.cxx',
'maincode.cpp'], include_dirs=['libint/include/libint2'],library_dirs=['libint/lib'],libraries=['int2'])
setup (name = 'maincode',
version = '1.0',
description = 'Simple maincode from SWIG tutorial',
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里面。

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

作者
Author: Shannon 时间: 2015-7-1 22:21

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的显卡可以玩）。

作者
Author: ximi1986 时间: 2015-7-2 09:34
Visual Studio 2012就可以,我在用，还是比较稳定。另外要装Fortran 2013和 vs 2012配套

作者
Author: 最爱喵星人 时间: 2019-3-27 20:39
请问怎么调用libint2的计算函数呢？
我看手册说通过Libint_eri_t 的结构体和libint2_build_eri[n][n][n][n]函数来计算双电子积分。请问怎么把基组信息喂给libint2呢？

作者
Author: Zork 时间: 2019-5-20 10:17
代码就是要够乱，看的才能提神。

作者
Author: granvia 时间: 2019-5-21 13:12
赞lz matlab最大缺点就是不免费

作者
Author: 413 时间: 2019-10-12 15:12
大佬，编译libint-2.0.3-stable版本成功后为啥生成的是libint2.h的头文件，而非libint.h？？？

作者
Author: CCCDDD 时间: 2022-3-1 15:01
你好，是不是少了一个方程的文件

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