diff --git a/configure b/configure index b217cbea..013ff541 100755 --- a/configure +++ b/configure @@ -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 diff --git a/install/scripts/install_eigen.sh b/install/scripts/install_eigen.sh new file mode 100755 index 00000000..a77ac1f6 --- /dev/null +++ b/install/scripts/install_eigen.sh @@ -0,0 +1,10 @@ +#!/bin/bash -x + +TARGET=eigen + +function _install() +{ + cp -R ${BUILD} . || exit 1 +} + +source scripts/build.sh diff --git a/install/scripts/install_libint.sh b/install/scripts/install_libint.sh new file mode 100755 index 00000000..384ad8b9 --- /dev/null +++ b/install/scripts/install_libint.sh @@ -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 + + diff --git a/plugins/Hartree_Fock/debug_libinit.irp.f b/plugins/Hartree_Fock/debug_libinit.irp.f new file mode 100644 index 00000000..586f393e --- /dev/null +++ b/plugins/Hartree_Fock/debug_libinit.irp.f @@ -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 diff --git a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py index 0dc99029..a140ec1e 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py @@ -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" + diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index d089e76b..df2a1dff 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -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")) diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index aa805093..67ccf03c 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -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 diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index eb443701..f2c365ae 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -4,6 +4,7 @@ double precision function ao_bielec_integral(i,j,k,l) ! integral of the AO basis 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) diff --git a/src/Integrals_Bielec/integral_bielec.cc b/src/Integrals_Bielec/integral_bielec.cc new file mode 100644 index 00000000..1ea9f160 --- /dev/null +++ b/src/Integrals_Bielec/integral_bielec.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + +// Libint Gaussian integrals library +#include + +/*** ================ ***/ +/*** 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& shells); +std::vector map_shell_to_basis_function(const std::vector& shells); +std::vector map_basis_function_to_shell(const std::vector& 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 shells_global; +std::vector shell2bf; +std::vector bf2shell; +static libint2::TwoBodyEngine *engine_pointer; + +// ___ _ +// | ._ _|_ _ ._ ._ _. | _|_ ._ _ _|_ o _ ._ +// _|_ | | |_ (/_ | | | (_| | | |_| | | (_ |_ | (_) | | +// + +size_t nbasis(const std::vector& shells) { + size_t n = 0; + for (const auto& shell: shells) + n += shell.size(); + return n; +} + +size_t max_nprim(const std::vector& shells) { + size_t n = 0; + for (auto shell: shells) + n = std::max(shell.nprim(), n); + return n; +} + +int max_l(const std::vector& shells) { + int l = 0; + for (auto shell: shells) + for (auto c: shell.contr) + l = std::max(c.l, l); + return l; +} + + +std::vector map_shell_to_basis_function(const std::vector& shells) { + std::vector result; + result.reserve(shells.size()); + + size_t n = 0; + for (auto shell: shells) { + result.push_back(n); + n += shell.size(); + } + + return result; +} + +std::vector map_basis_function_to_shell(const std::vector& shells) { + + std::vector 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 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 (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 &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 &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]; + }; diff --git a/src/Integrals_Bielec/libint_module.f90 b/src/Integrals_Bielec/libint_module.f90 new file mode 100644 index 00000000..2cedf283 --- /dev/null +++ b/src/Integrals_Bielec/libint_module.f90 @@ -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