10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 15:12:14 +02:00

Merge pull request #33 from TApplencourt/master

Libint
This commit is contained in:
Anthony Scemama 2016-03-07 20:47:25 +01:00
commit 064998dba7
10 changed files with 554 additions and 15 deletions

21
configure vendored
View File

@ -72,7 +72,9 @@ d_dependency = {
"ninja": ["g++", "python"],
"make": [],
"p_graphviz": ["python"],
"bats": []
"bats": [],
"eigen": [],
"libint": ["eigen"]
}
from collections import namedtuple
@ -148,23 +150,32 @@ f77zmq = Info(
# join(QP_ROOT, "src", "ZMQ", "f77zmq.h") )
p_graphviz = Info(
url='https://github.com/xflr6/graphviz/archive/master.tar.gz',
url='{head}/xflr6/graphviz/{tail}'.format(**path_github),
description=' Python library for graphviz',
default_path=join(QP_ROOT_INSTALL, "p_graphviz"))
bats = Info(
url='https://github.com/sstephenson/bats/archive/master.tar.gz',
url='{head}/sstephenson/bats/{tail}'.format(**path_github),
description=' Bash Automated Testing System',
default_path=join(QP_ROOT_INSTALL, "bats"))
libint = Info(
url='{head}/evaleev/libint/releases/download/v2.1.0-beta2/libint-2.1.0-beta2.tgz'.format(**path_github),
description=' Libint is a high-performance library for computing Gaussian integrals in quantum mechanics',
default_path=join(QP_ROOT_INSTALL, "libint"))
eigen = Info(
url='http://bitbucket.org/eigen/eigen/get/3.2.8.tar.gz',
description=' Eigen is a C++ template library for linear algebra.',
default_path=join(QP_ROOT_INSTALL, "eigen"))
d_info = dict()
for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt",
"resultsFile", "ninja", "emsl", "ezfio", "p_graphviz",
"zeromq", "f77zmq","bats" ]:
"zeromq", "f77zmq","bats","libint","eigen"]:
exec ("d_info['{0}']={0}".format(m))
def find_path(bin_, l_installed, var_for_qp_root=False):
"""Use the global variable
* l_installed

View File

@ -0,0 +1,10 @@
#!/bin/bash -x
TARGET=eigen
function _install()
{
cp -R ${BUILD} . || exit 1
}
source scripts/build.sh

View File

@ -0,0 +1,33 @@
#!/bin/bash -x
TARGET=libint
function _install()
{
cd ..
QP_ROOT=$PWD
cd -
cp -R ${BUILD} . || exit 1
cd ${TARGET}
export CXX="g++"
export CXXFLAGS=" -O3 -std=c++0x"
./configure --with-cxx-optflags
make -j 8
cd -
g++ -O2 -std=c++0x -DHAVE_CONFIG_H -I${PWD}/${TARGET}/include -I${QP_ROOT}/install/eigen -DPREP_LIBINT2_SKIP_BOOST -L${PWD}/${TARGET}/lib -lint2 -c ${QP_ROOT}/src/Integrals_Bielec/integral_bielec.cc
mv integral_bielec.o ${QP_ROOT}/src/Integrals_Bielec/
}
BUILD=_build/${TARGET}
rm -rf -- ${BUILD}
mkdir ${BUILD} || exit 1
tar -zxf Downloads/${TARGET}.tgz --strip-components=1 --directory=${BUILD} || exit 1
_install || exit 1
rm -rf -- ${BUILD} _build/${TARGET}.log
exit 0

View File

@ -0,0 +1,124 @@
program debug_libint
use libint_module
implicit none
double precision :: ao_bielec_integral
integer, allocatable :: s2bf(:)
double precision, allocatable :: buffer_int(:)
call init_libint(trim(ezfio_filename)//char(0))
integer :: nb_shell_f
nb_shell_f = get_nb_shell()
allocate(s2bf(2*nb_shell_f))
call map_shell_to_basis_function_interval(2*nb_shell_f,s2bf)
integer :: s1, s2,s3,s4
integer :: bf1,bf2,bf3,bf4
integer :: bf1_begin,bf2_begin,bf3_begin,bf4_begin
integer :: bf1_end,bf2_end,bf3_end,bf4_end
integer :: n1,n2,n3,n4
integer :: f1,f2,f3,f4,f1234
! =================== !
! Loop over the shell !
! =================== !
do s1 = 1, nb_shell_f
print*, s1, "/", nb_shell_f
bf1_begin = s2bf(2*s1-1)
bf1_end = s2bf(2*s1)
n1 = 1 + bf1_end - bf1_begin
do s2 = 1, nb_shell_f
bf2_begin = s2bf(2*s2-1)
bf2_end = s2bf(2*s2)
n2 = 1 + bf2_end - bf2_begin
do s3 = 1, nb_shell_f
bf3_begin = s2bf(2*s3-1)
bf3_end = s2bf(2*s3)
n3 = 1 + bf3_end - bf3_begin
do s4 = 1, nb_shell_f
bf4_begin = s2bf(2*s4-1)
bf4_end = s2bf(2*s4)
n4 = 1 + bf4_end - bf4_begin
! ========================== !
! Compute the shell integral !
! ========================== !
integer :: sze
sze = n1*n2*n3*n4
allocate(buffer_int(sze))
call compute_ao_bielec_integrals_shell(s1,s2,s3,s4,sze,buffer_int)
! ============================ !
! Loop over the basis function !
! ============================ !
do bf1 = bf1_begin, bf1_end
do bf2 = bf2_begin, bf2_end
do bf3 = bf3_begin, bf3_end
do bf4 = bf4_begin, bf4_end
f1 = bf1 - bf1_begin
f2 = bf2 - bf2_begin
f3 = bf3 - bf3_begin
f4 = bf4 - bf4_begin
!Get the integral from the buffer
f1234 = f1*n2*n3*n4+f2*n3*n4+f3*n4+f4 + 1;
!Compute the norm
double precision:: coef1, coef2, coef3, coef4, norm
coef1 = ao_coef_normalization_libint_factor(bf1)
coef2 = ao_coef_normalization_libint_factor(bf2)
coef3 = ao_coef_normalization_libint_factor(bf3)
coef4 = ao_coef_normalization_libint_factor(bf4)
norm = coef1*coef2*coef3*coef4
double precision:: libint, ref
!Value of itegral bf1,bf2,bf3, bf4
libint = buffer_int(f1234) * norm
!Verify with the manu's one
! ref = ao_bielec_integral(bf1,bf2,bf3,bf4)
!
! if ( (ABS(ABS(ref) - ABS(libint)).GE.1e-6) ) THEN
!
! print*, bf1,bf2,bf3,bf4
! print*,"r", ref
! print*,"l", libint
! print*,"r/l", ref/libint
! print*,"l/r", libint/ref
! print*,"n", norm
!
! call exit(1)
! end if
enddo
enddo
enddo
enddo
!Deallocate the buffer_intergral for the shell
deallocate(buffer_int)
enddo
enddo
enddo
enddo
call finalize_libint()
end debug_libint

View File

@ -183,6 +183,9 @@ def get_nb_permutation(str_):
def order_l_l_sym(l_l_sym):
l_order_mo = [i for i,_ in enumerate(l_l_sym)]
n = 1
for i in range(len(l_l_sym)):
if n != 1:
@ -192,11 +195,11 @@ def order_l_l_sym(l_l_sym):
l = l_l_sym[i]
n = get_nb_permutation(l[2])
l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n],
key=lambda x: x[2],
cmp=compare_gamess_style)
l_l_sym[i:i + n], l_order_mo[i:i+n] = zip(*sorted(zip(l_l_sym[i:i + n],l_order_mo[i:i+n]),
key=lambda x: x[0][2],
cmp=compare_gamess_style))
return l_l_sym
return l_l_sym, l_order_mo
#==========================
@ -205,8 +208,13 @@ def order_l_l_sym(l_l_sym):
l_sym_without_header = sym_raw.split("\n")[3:-2]
l_l_sym_raw = [i.split() for i in l_sym_without_header]
print len(l_l_sym_raw)
l_l_sym_expend_sym = expend_sym_l(l_l_sym_raw)
l_l_sym_ordered = order_l_l_sym(l_l_sym_expend_sym)
print len(l_l_sym_expend_sym)
l_l_sym_ordered, l_order_mo = order_l_l_sym(l_l_sym_expend_sym)
#========
#MO COEF
@ -348,6 +356,7 @@ d_rep={"+":"1","-":"0"}
det_without_header = det_raw[pos+2::]
for line_raw in det_without_header.split("\n"):
line = line_raw
@ -355,8 +364,14 @@ for line_raw in det_without_header.split("\n"):
try:
float(line)
except ValueError:
print line_raw.strip(), len(line_raw.strip())
print l_order_mo, len(l_order_mo)
line_order = [line_raw[i] for i in l_order_mo]
line= "".join([d_rep[x] if x in d_rep else x for x in line_raw])
print line.strip()
print "END_DET"

View File

@ -37,7 +37,8 @@ from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
LIB = "" # join(QP_ROOT, "lib", "rdtsc.o")
EZFIO_LIB = join(QP_ROOT, "lib", "libezfio_irp.a")
ZMQ_LIB = join(QP_ROOT, "lib", "libf77zmq.a") + " " + join(QP_ROOT, "lib", "libzmq.a") + " -lstdc++ -lrt"
ZMQ_LIB = join(QP_ROOT, "lib", "libf77zmq.a") + " " + join(QP_ROOT, "lib", "libzmq.a") + " -lstdc++ -lrt"
INT_LIB = join(QP_ROOT, "libint","lib", ".libs", "libint2.a")
ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja")
header = r"""#
@ -96,7 +97,8 @@ def ninja_create_env_variable(pwd_config_file):
l_string.append(str_)
lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB")
l_string.append("LIB = {0} {1} {2} {3}".format(LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB))
str_lib = " ".join([LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB, INT_LIB])
l_string.append("LIB = {0} ".format(str_lib))
l_string.append("")
@ -387,6 +389,8 @@ def get_l_file_for_module(path_module):
l_src.append(f)
obj = '{0}.o'.format(os.path.splitext(f)[0])
l_obj.append(obj)
elif f.lower().endswith(".o"):
l_obj.append(join(path_module.abs, f))
elif f == "EZFIO.cfg":
l_depend.append(join(path_module.abs, "ezfio_interface.irp.f"))

View File

@ -25,7 +25,7 @@ END_PROVIDER
BEGIN_DOC
! Coefficients including the AO normalization
END_DOC
double precision :: norm, norm2,overlap_x,overlap_y,overlap_z,C_A(3), c
double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
integer :: l, powA(3), nz
integer :: i,j,k
nz=100
@ -34,9 +34,11 @@ END_PROVIDER
C_A(3) = 0.d0
ao_coef_normalized = 0.d0
do i=1,ao_num
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
do j=1,ao_prim_num(i)
call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm)
@ -51,8 +53,42 @@ END_PROVIDER
enddo
ao_coef_normalization_factor(i) = 1.d0/sqrt(norm)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num) ]
implicit none
BEGIN_DOC
! Coefficients including the AO normalization
END_DOC
double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
integer :: l, powA(3), nz
integer :: i,j,k
nz=100
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
do i=1,ao_num
powA(1) = ao_l(i)
powA(2) = 0
powA(3) = 0
! Normalization of the contracted basis functions
norm = 0.d0
do j=1,ao_prim_num(i)
do k=1,ao_prim_num(i)
call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
enddo
enddo
ao_coef_normalization_libint_factor(i) = ao_coef_normalization_factor(i) * sqrt(norm)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_expo_ordered, (ao_num_align,ao_prim_num_max) ]
implicit none

View File

@ -4,6 +4,7 @@ double precision function ao_bielec_integral(i,j,k,l)
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC
integer,intent(in) :: i,j,k,l
integer :: p,q,r,s
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
@ -288,11 +289,18 @@ end
subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
implicit none
use map_module
use libint_module
BEGIN_DOC
! Compute AO 1/r12 integrals for all i and fixed j,k,l
END_DOC
! include 'Utils/constants.include.F'
! integer, intent(in) :: j,k,l,sze
! real(integral_kind), intent(out) :: buffer_value(sze)
!
! call compute_ao_bielec_integrals_libint(j,k,l,sze,buffer_value)
include 'Utils/constants.include.F'
integer, intent(in) :: j,k,l,sze
real(integral_kind), intent(out) :: buffer_value(sze)

View File

@ -0,0 +1,247 @@
/*
* This file is a part of Libint.
* Copyright (C) 2004-2014 Edward F. Valeev
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see http://www.gnu.org/licenses/.
*
*/
#if __cplusplus <= 199711L
# error "Hartree-Fock test requires C++11 support"
#endif
// standard C++ headers
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <chrono>
#include <stdlib.h>
// Libint Gaussian integrals library
#include <libint2.hpp>
/*** ================ ***/
/*** Exposed Function ***/
/*** ================ ***/
extern "C"
{
void init_libint(char ezfio_filename[]);
void finalize_libint();
int nb_shell();
void map_shell_to_basis_function_interval(int sze, int* out_val);
double ao_bielec_integral(int bf1f, int bf2f, int bf3f, int bf4f);
void compute_ao_bielec_integrals_jkl(int bf2, int bf3, int bf4, int sze, double* values);
void compute_ao_bielec_integrals_shell(int s1, int s2, int s3, int s4, int sze, double* values);
}
using libint2::Shell;
/*** ================= ***/
/*** Internal Function ***/
/*** ================= ***/
size_t nbasis(const std::vector<Shell>& shells);
std::vector<size_t> map_shell_to_basis_function(const std::vector<Shell>& shells);
std::vector<size_t> map_basis_function_to_shell(const std::vector<Shell>& shells);
/*** ================ ***/
/*** Exposed Function ***/
/*** ================ ***/
void init_libint(char ezfio_filename[]);
void finalize_libint();
int nb_shell();
void map_shell_to_basis_function_interval(int sze, int* out_val);
double ao_bielec_integral(int bf1f, int bf2f, int bf3f, int bf4f);
void compute_ao_bielec_integrals_jkl(int i, int j, int k, int sze, double* values);
void compute_ao_bielec_integrals_shell(int s1, int s2, int s3, int s4, int sze, double* values);
/*** =============== ***/
/*** Global Variable ***/
/*** =============== ***/
std::vector<Shell> shells_global;
std::vector<size_t> shell2bf;
std::vector<size_t> bf2shell;
static libint2::TwoBodyEngine<libint2::Coulomb> *engine_pointer;
// ___ _
// | ._ _|_ _ ._ ._ _. | _|_ ._ _ _|_ o _ ._
// _|_ | | |_ (/_ | | | (_| | | |_| | | (_ |_ | (_) | |
//
size_t nbasis(const std::vector<Shell>& shells) {
size_t n = 0;
for (const auto& shell: shells)
n += shell.size();
return n;
}
size_t max_nprim(const std::vector<Shell>& shells) {
size_t n = 0;
for (auto shell: shells)
n = std::max(shell.nprim(), n);
return n;
}
int max_l(const std::vector<Shell>& shells) {
int l = 0;
for (auto shell: shells)
for (auto c: shell.contr)
l = std::max(c.l, l);
return l;
}
std::vector<size_t> map_shell_to_basis_function(const std::vector<Shell>& shells) {
std::vector<size_t> result;
result.reserve(shells.size());
size_t n = 0;
for (auto shell: shells) {
result.push_back(n);
n += shell.size();
}
return result;
}
std::vector<size_t> map_basis_function_to_shell(const std::vector<Shell>& shells) {
std::vector<size_t> result;
result.reserve(nbasis(shells));
size_t n = 0;
for (auto shell: shells) {
for (auto i=0; i!=shell.size(); ++i){
result.push_back(n);
}
n++;
}
return result;
}
// _ _
// |_ ._ _ _ _ _| _|_ ._ _ _|_ o _ ._
// |_ >< |_) (_) _> (/_ (_| | |_| | | (_ |_ | (_) | |
// |
void init_libint(char ezfio_filename[]){
/*** =========================== ***/
/*** initialize molecule ***/
/*** =========================== ***/
std::string xyz_path = ezfio_filename + std::string("/libint/xyz");
// read geometry from a filename
std::ifstream input_file(xyz_path);
std::vector<libint2::Atom> atoms = libint2::read_dotxyz(input_file);
/*** =========================== ***/
/*** create basis set ***/
/*** =========================== ***/
std::string basis_path = ezfio_filename + std::string("/libint");
setenv("LIBINT_DATA_PATH", basis_path.c_str(), 1);
libint2::BasisSet shells("basis", atoms);
shells_global = shells;
for(auto& shell: shells_global)
for(auto& contraction: shell.contr)
contraction.pure = false;
// initializes the Libint integrals library ... now ready to compute
libint2::init();
// construct the electron repulsion integrals engine
engine_pointer = new libint2::TwoBodyEngine<libint2::Coulomb> (max_nprim(shells_global), max_l(shells_global), 0);
shell2bf = map_shell_to_basis_function(shells_global);
bf2shell = map_basis_function_to_shell(shells_global);
}
void finalize_libint(){
libint2::finalize(); // done with libint2
}
int nb_shell(){
return shells_global.size();
}
void map_shell_to_basis_function_interval(int sze, int* out_val) {
size_t n = 1;
for (auto i=0; i<shells_global.size() ; i++) {
out_val[2*i] = n;
n += shells_global[i].size();
out_val[2*i+1] = n-1;
}
}
double ao_bielec_integral(int bf1f, int bf2f, int bf3f, int bf4f){
auto bf1 = bf1f-1;
auto bf2 = bf2f-1;
auto bf3 = bf3f-1;
auto bf4 = bf4f-1;
// construct the electron repulsion integrals engine
libint2::TwoBodyEngine<libint2::Coulomb> &engine = *engine_pointer;
auto s1 = bf2shell[bf1];
auto n1 = shells_global[s1].size();
auto f1 = bf1-shell2bf[s1];
auto s2 = bf2shell[bf2];
auto n2 = shells_global[s2].size();
auto f2 = bf2-shell2bf[s2];
auto s3 = bf2shell[bf3];
auto n3 = shells_global[s3].size();
auto f3 = bf3-shell2bf[s3];;
auto s4 = bf2shell[bf4];
auto n4 = shells_global[s4].size();
auto f4 = bf4- shell2bf[s4];
// std::cout << "o: compute shell set {" << s1 << "," << s2 <<"," << s3 <<"," << s4 << "} ... ";
const auto* buf_1234 = engine.compute(shells_global[s1], shells_global[s2], shells_global[s3], shells_global[s4]);
// std::cout << "done" << std::endl;
auto f1234 = f1*n2*n3*n4+f2*n3*n4+f3*n4+f4;
auto result = buf_1234[f1234];
return result;
};
void compute_ao_bielec_integrals_shell(int s1, int s2, int s3, int s4, int sze, double* values){
libint2::TwoBodyEngine<libint2::Coulomb> &engine = *engine_pointer;
const auto* buf_1234 = engine.compute(shells_global[s1-1], shells_global[s2-1], shells_global[s3-1], shells_global[s4-1]);
for(auto i=0; i!=sze; i++)
values[i] = buf_1234[i];
};

View File

@ -0,0 +1,51 @@
module libint_module
use iso_c_binding
implicit none
interface
subroutine init_libint(str) bind(c,name='init_libint')
import :: c_char
character(len=1,kind=C_char), dimension(*), intent(in) :: str
end subroutine
integer(c_int) function get_nb_shell() bind(c,name='nb_shell')
import :: c_int
end function
subroutine finalize_libint() bind(c,name='finalize_libint')
end subroutine
subroutine map_shell_to_basis_function_interval(sze, out_val) bind(c,name='map_shell_to_basis_function_interval')
import :: c_ptr
import :: c_int
integer(c_int), INTENT(IN), value :: sze
integer(c_int), INTENT(OUT) :: out_val(sze)
end subroutine
real(c_double) function ao_bielec_integral_libint(i,j,k,l) bind(c,name='ao_bielec_integral')
import :: c_int
import :: c_double
integer(c_int), value :: i
integer(c_int), value :: j
integer(c_int), value :: k
integer(c_int), value :: l
end function
subroutine compute_ao_bielec_integrals_shell(i,j,k,l,sze,values) bind(c,name='compute_ao_bielec_integrals_shell')
import :: c_ptr
import :: c_int
import :: c_double
integer(c_int), value :: i
integer(c_int), value :: j
integer(c_int), value :: k
integer(c_int), value :: l
integer(c_int), INTENT(IN), value :: sze
real(c_double), INTENT(OUT) :: values(sze)
end subroutine
end interface
end module libint_module