diff --git a/configure b/configure index 250ae9d0..a8f6cc9e 100755 --- a/configure +++ b/configure @@ -77,9 +77,7 @@ d_dependency = { "ninja": ["g++", "python"], "make": [], "p_graphviz": ["python"], - "bats": [], - "eigen": [], - "libint": ["eigen"] + "bats": [] } from collections import namedtuple @@ -155,32 +153,23 @@ f77zmq = Info( # join(QP_ROOT, "src", "ZMQ", "f77zmq.h") ) p_graphviz = Info( - url='{head}/xflr6/graphviz/{tail}'.format(**path_github), + url='https://github.com/xflr6/graphviz/archive/master.tar.gz', description=' Python library for graphviz', default_path=join(QP_ROOT_INSTALL, "p_graphviz")) bats = Info( - url='{head}/sstephenson/bats/{tail}'.format(**path_github), + url='https://github.com/sstephenson/bats/archive/master.tar.gz', 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","libint","eigen"]: + "zeromq", "f77zmq","bats"]: 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 deleted file mode 100755 index a77ac1f6..00000000 --- a/install/scripts/install_eigen.sh +++ /dev/null @@ -1,10 +0,0 @@ -#!/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 deleted file mode 100755 index 384ad8b9..00000000 --- a/install/scripts/install_libint.sh +++ /dev/null @@ -1,33 +0,0 @@ -#!/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/ocaml/Libint.ml b/ocaml/Libint.ml deleted file mode 100644 index 0b88ab41..00000000 --- a/ocaml/Libint.ml +++ /dev/null @@ -1,85 +0,0 @@ -open Core.Std -open Qptypes - -let molecule = lazy ( - let atom_list = - match Input.Nuclei.read () with - | Some data -> Input.Nuclei.to_atom_list data - | None -> failwith "No coordinate found" - and data = - match Input.Electrons.read () with - | Some data -> data - | None -> failwith "No electrons found" - in - { Molecule.nuclei = atom_list; - Molecule.elec_alpha = data.Input.Electrons.elec_alpha_num; - Molecule.elec_beta = data.Input.Electrons.elec_beta_num - } -) - -let write_xyz_file libint_dir = - assert (Sys.is_directory_exn libint_dir); - let filename = - Filename.concat libint_dir "xyz" - and molecule = - Lazy.force molecule - in - Out_channel.with_file filename ~f:(fun oc -> - Molecule.to_xyz molecule - |> Out_channel.output_string oc - ) - -let write_basis_file libint_dir = - assert (Sys.is_directory_exn libint_dir); - let filename = - Filename.concat libint_dir "basis.g94" - and molecule = - Lazy.force molecule - in - let text = - let rec substitute accu i = function - | ele :: tail -> - let pattern = - Printf.sprintf "Atom %d" i - in - let new_string = - String.substr_replace_first accu ~pattern ~with_:(Printf.sprintf "%s 0" ele) - in - substitute new_string (i+1) tail - | [] -> accu - in - let accu = - let b = - match Input.Ao_basis.read () with - | Some data -> Input.Ao_basis.to_basis data - | None -> failwith "No AO basis" - in - Basis.to_string ~fmt:Gto.Gaussian b - and atom_names = - List.map molecule.Molecule.nuclei ~f:(fun x -> Element.to_string x.Atom.element) - in - substitute accu 1 atom_names - in - Out_channel.with_file filename ~f:(fun oc -> - Out_channel.output_string oc text - ) - - -let write_files ezfio_filename = - assert (Sys.is_directory_exn ezfio_filename); - - let libint_dir = - Filename.concat ezfio_filename "libint" - in - - let () = - match Sys.is_directory libint_dir with - | `Yes -> () - | `No -> Unix.mkdir libint_dir - | `Unknown -> failwith ("Unable to tell if "^libint_dir^" exists.") - in - - write_xyz_file libint_dir; - write_basis_file libint_dir - - diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 3dc1e0f0..596c947a 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -7,9 +7,14 @@ s.set_selection_pt2("epstein_nesbet_2x2") #s.unset_openmp() print s +#s = H_apply("FCI_PT2") +#s.set_perturbation("epstein_nesbet_2x2") +#s.unset_openmp() +#print s + s = H_apply_zmq("FCI_PT2") s.set_perturbation("epstein_nesbet_2x2") -#s.unset_openmp() +s.unset_openmp() print s s = H_apply("FCI_no_skip") diff --git a/plugins/Full_CI/micro_pt2.irp.f b/plugins/Full_CI/micro_pt2.irp.f index d78a942d..1d4d7eaa 100644 --- a/plugins/Full_CI/micro_pt2.irp.f +++ b/plugins/Full_CI/micro_pt2.irp.f @@ -24,6 +24,8 @@ subroutine run_wf integer(ZMQ_PTR) :: zmq_to_qp_run_socket print *, 'Getting wave function' + zmq_context = f77_zmq_ctx_new () + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() call zmq_get_psi(zmq_to_qp_run_socket, 1) diff --git a/plugins/Hartree_Fock/debug_libint.irp.f b/plugins/Hartree_Fock/debug_libint.irp.f deleted file mode 100644 index 7984d62c..00000000 --- a/plugins/Hartree_Fock/debug_libint.irp.f +++ /dev/null @@ -1,115 +0,0 @@ - program debug_libint - - implicit none - double precision :: ao_bielec_integral - - double precision, allocatable :: buffer_int(:) - - 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 - - PROVIDE has_libint - ! =================== ! - ! Loop over the shell ! - ! =================== ! - - do s1 = 1, shell_num - - print*, s1, "/", shell_num - - bf1_begin = shell_idx(1,s1) - bf1_end = shell_idx(2,s1) - n1 = 1 + bf1_end - bf1_begin - - do s2 = 1, shell_num - - bf2_begin = shell_idx(1,s2) - bf2_end = shell_idx(2,s2) - n2 = 1 + bf2_end - bf2_begin - - do s3 = 1, shell_num - - bf3_begin = shell_idx(1,s3) - bf3_end = shell_idx(2,s3) - n3 = 1 + bf3_end - bf3_begin - - do s4 = 1, shell_num - - bf4_begin = shell_idx(1,s4) - bf4_end = shell_idx(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)) >= 1.e-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 diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f index 38cd3c55..0d49be89 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -65,8 +65,17 @@ subroutine run_pt2(N_st,energy) threshold_selectors = 1.d0 threshold_generators = 0.999d0 - N_det_generators = lambda_mrcc_pt2(0) - do i=1,N_det_generators + N_det_generators = lambda_mrcc_pt2(0) + N_det_cas + do i=1,N_det_cas + do k=1,N_int + psi_det_generators(k,1,i) = psi_ref(k,1,i) + psi_det_generators(k,2,i) = psi_ref(k,2,i) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_ref_coef(i,k) + enddo + enddo + do i=N_det_cas+1,N_det_generators j = lambda_mrcc_pt2(i) do k=1,N_int psi_det_generators(k,1,i) = psi_non_ref(k,1,j) diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index df2b67a0..449d262c 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -25,11 +25,17 @@ print s -s = H_apply_zmq("mrcc_PT2") +s = H_apply("mrcc_PT2") s.energy = "ci_electronic_energy_dressed" s.set_perturbation("epstein_nesbet_2x2") s.unset_openmp() print s +#s = H_apply_zmq("mrcc_PT2") +#s.energy = "ci_electronic_energy_dressed" +#s.set_perturbation("epstein_nesbet_2x2") +#s.unset_openmp() +#print s + END_SHELL diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index b2304818..1c2e8b74 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -274,7 +274,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge enddo do i_state=1,Nstates - ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state) + ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) enddo do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) @@ -285,16 +285,20 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge enddo enddo call omp_set_lock( psi_ref_lock(i_I) ) - do l_sd=1,idx_alpha(0) - k_sd = idx_alpha(l_sd) - do i_state=1,Nstates - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) - if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then - delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) - else - delta_ii_(i_state,i_I) = 0.d0 - endif - enddo + do i_state=1,Nstates + if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + enddo + else + delta_ii_(i_state,i_I) = 0.d0 + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + enddo + endif enddo call omp_unset_lock( psi_ref_lock(i_I) ) enddo diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index 6beca8f2..d3b6c28f 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -24,6 +24,20 @@ use bitmasks enddo enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_ref_coef_inv, (psi_det_size,n_states) ] + implicit none + BEGIN_DOC + ! 1/psi_ref_coef + END_DOC + integer :: i, i_state + do i_state=1,N_states + do i=1,N_det_ref + psi_ref_coef_inv(i,i_state) = 1.d0/psi_ref_coef(i,i_state) + enddo + enddo + END_PROVIDER diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 41709c8f..cc1c8aa8 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -38,7 +38,6 @@ 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" -INT_LIB = join(QP_ROOT, "install", "libint","lib", ".libs", "libint2.a") ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja") header = r"""# @@ -97,7 +96,7 @@ def ninja_create_env_variable(pwd_config_file): l_string.append(str_) lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB") - str_lib = " ".join([LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB, INT_LIB]) + str_lib = " ".join([LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB]) l_string.append("LIB = {0} ".format(str_lib)) l_string.append("") diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index a3f3600b..d7cd9c95 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -345,7 +345,7 @@ def save_ezfio_provider(path_head, dict_code_provider): path = "{0}/ezfio_interface.irp.f".format(path_head) l_output = ["! DO NOT MODIFY BY HAND", - "! Created by $QP_ROOT/scripts/ezfio_interface.py", + "! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py", "! from file {0}/EZFIO.cfg".format(path_head), "\n"] diff --git a/scripts/ezfio_interface/ezfio_generate_provider.py b/scripts/ezfio_interface/ezfio_generate_provider.py index 6cd919dc..89fdfa03 100755 --- a/scripts/ezfio_interface/ezfio_generate_provider.py +++ b/scripts/ezfio_interface/ezfio_generate_provider.py @@ -22,6 +22,7 @@ BEGIN_PROVIDER [ %(type)s, %(name)s %(size)s ] logical :: has PROVIDE ezfio_filename + %(test_null_size)s call ezfio_has_%(ezfio_dir)s_%(ezfio_name)s(has) if (has) then call ezfio_get_%(ezfio_dir)s_%(ezfio_name)s(%(name)s) @@ -44,6 +45,7 @@ END_PROVIDER def __repr__(self): self.set_write() + self.set_test_null_size() for v in self.values: if not v: msg = "Error : %s is not set in EZFIO.cfg" % (v) @@ -54,20 +56,31 @@ END_PROVIDER return self.data % self.__dict__ + def set_test_null_size(self): + if "size" not in self.__dict__: + self.__dict__["size"] = "" + if self.size != "": + self.test_null_size = "if (size(%s) == 0) return\n" % ( self.name ) + else: + self.test_null_size = "" + def set_write(self): self.write = "" - if self.type in self.write_correspondance: - write = self.write_correspondance[self.type] - output = self.output - name = self.name + if "size" in self.__dict__: + return + else: + if self.type in self.write_correspondance: + write = self.write_correspondance[self.type] + output = self.output + name = self.name - l_write = ["", - " call write_time(%(output)s)", - " call %(write)s(%(output)s, %(name)s, &", - " '%(name)s')", - ""] + l_write = ["", + " call write_time(%(output)s)", + " call %(write)s(%(output)s, %(name)s, &", + " '%(name)s')", + ""] - self.write = "\n".join(l_write) % locals() + self.write = "\n".join(l_write) % locals() def set_type(self, t): self.type = t.lower() diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index 212c1648..9c7a1386 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -180,7 +180,6 @@ let run check_only ezfio_filename = | None -> "vi" in - Libint.write_files (!Ezfio.ezfio_filename); match check_only with | true -> () | false -> diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 436f092d..90f997ed 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -401,7 +401,16 @@ class H_apply_zmq(H_apply): H_pert_diag(k) = 0.d0 norm_psi(k) = 0.d0 enddo - """ + """ + self.data["copy_buffer"] = """ + do i=1,N_det_generators + do k=1,N_st + pt2(k) = pt2(k) + pt2_generators(k,i) + norm_pert(k) = norm_pert(k) + norm_pert_generators(k,i) + H_pert_diag(k) = H_pert_diag(k) + H_pert_diag_generators(k,i) + enddo + enddo + """ def set_selection_pt2(self,pert): H_apply.set_selection_pt2(self,pert) @@ -416,3 +425,4 @@ class H_apply_zmq(H_apply): select_max(i_generator) = 0.d0 endif """ + diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index 28513597..abb251db 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -307,14 +307,14 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) end -subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) +subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,i_generator,N_st,task_id) use f77_zmq implicit none BEGIN_DOC ! Push PT2 calculation to the collector END_DOC integer(ZMQ_PTR), intent(in) :: zmq_socket_push - integer, intent(in) :: N_st + integer, intent(in) :: N_st, i_generator double precision, intent(in) :: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st) integer, intent(in) :: task_id integer :: rc @@ -343,6 +343,12 @@ subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) stop 'error' endif + rc = f77_zmq_send( zmq_socket_push, i_generator, 4, ZMQ_SNDMORE) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_send( zmq_socket_push, i_generator, 4, 0)' + stop 'error' + endif + rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) if (rc /= 4) then print *, irp_here, 'f77_zmq_send( zmq_socket_push, task_id, 4, 0)' @@ -358,7 +364,7 @@ subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) ! endif end -subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,N_st,n,task_id) +subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n,task_id) use f77_zmq implicit none BEGIN_DOC @@ -368,7 +374,7 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,N_st,n,task_id) integer, intent(in) :: N_st double precision, intent(out) :: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st) integer, intent(out) :: task_id - integer, intent(out) :: n + integer, intent(out) :: n, i_generator integer :: rc n=0 @@ -386,7 +392,11 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,N_st,n,task_id) rc = f77_zmq_recv( zmq_socket_pull, pt2(1), 8*N_st, 0) if (rc /= 8*N_st) then - print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, pt2(1,1) , 8*N_st, 0)' + print *, '' + print *, '' + print *, '' + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, pt2(1) , 8*N_st, 0)' + print *, rc stop 'error' endif @@ -402,6 +412,12 @@ subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,N_st,n,task_id) stop 'error' endif + rc = f77_zmq_recv( zmq_socket_pull, i_generator, 4, 0) + if (rc /= 4) then + print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, i_generator, 4, 0)' + stop 'error' + endif + rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) if (rc /= 4) then print *, irp_here, 'f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)' diff --git a/src/Determinants/H_apply_zmq.template.f b/src/Determinants/H_apply_zmq.template.f index c492a739..c1f6ceed 100644 --- a/src/Determinants/H_apply_zmq.template.f +++ b/src/Determinants/H_apply_zmq.template.f @@ -10,6 +10,7 @@ subroutine $subroutine($params_main) $decls_main + integer :: i integer :: i_generator double precision :: wall_0, wall_1 integer(omp_lock_kind) :: lck @@ -26,6 +27,9 @@ subroutine $subroutine($params_main) integer(ZMQ_PTR) :: zmq_socket_pair integer(ZMQ_PTR) :: zmq_to_qp_run_socket + double precision, allocatable :: pt2_generators(:,:), norm_pert_generators(:,:) + double precision, allocatable :: H_pert_diag_generators(:,:) + call new_parallel_job(zmq_to_qp_run_socket,'$subroutine') zmq_socket_pair = new_zmq_pair_socket(.True.) @@ -37,24 +41,26 @@ subroutine $subroutine($params_main) call add_task_to_taskserver(zmq_to_qp_run_socket,task) enddo - integer(ZMQ_PTR) :: collector_thread - external :: $subroutine_collector - rc = pthread_create(collector_thread, $subroutine_collector) + allocate ( pt2_generators(N_states,N_det_generators), & + norm_pert_generators(N_states,N_det_generators), & + H_pert_diag_generators(N_states,N_det_generators) ) - !$OMP PARALLEL DEFAULT(private) - !$OMP TASK PRIVATE(rc) - rc = omp_get_thread_num() - call $subroutine_slave_inproc(rc) - !$OMP END TASK - !$OMP TASKWAIT + PROVIDE nproc N_states + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i) & + !$OMP SHARED(zmq_socket_pair,N_states, pt2_generators, norm_pert_generators, H_pert_diag_generators, n, task_id, i_generator) & + !$OMP num_threads(nproc+1) + i = omp_get_thread_num() + if (i == 0) then + call $subroutine_collector() + integer :: n, task_id + call pull_pt2(zmq_socket_pair, pt2_generators, norm_pert_generators, H_pert_diag_generators, i_generator, size(pt2_generators), n, task_id) + else + call $subroutine_slave_inproc(i) + endif !$OMP END PARALLEL - integer :: n, task_id - call pull_pt2(zmq_socket_pair, pt2, norm_pert, H_pert_diag, N_st, n, task_id) - - rc = pthread_join(collector_thread) - call end_zmq_pair_socket(zmq_socket_pair) call end_parallel_job(zmq_to_qp_run_socket,'$subroutine') @@ -62,6 +68,7 @@ subroutine $subroutine($params_main) $copy_buffer $generate_psi_guess + deallocate ( pt2_generators, norm_pert_generators, H_pert_diag_generators) end subroutine $subroutine_slave_tcp(iproc) @@ -169,7 +176,7 @@ subroutine $subroutine_slave(thread, iproc) endif call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,1) - call push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,N_st,task_id) + call push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,i_generator,N_st,task_id) enddo @@ -186,7 +193,7 @@ subroutine $subroutine_collector use f77_zmq implicit none BEGIN_DOC -! Collects results from the selection +! Collects results from the selection in an array of generators END_DOC integer :: k, rc @@ -194,7 +201,7 @@ subroutine $subroutine_collector integer(ZMQ_PTR), external :: new_zmq_pull_socket integer(ZMQ_PTR) :: zmq_socket_pull integer*8 :: control, accu - integer :: n, more, task_id + integer :: n, more, task_id, i_generator integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket integer(ZMQ_PTR) :: zmq_to_qp_run_socket @@ -202,22 +209,25 @@ subroutine $subroutine_collector zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_socket_pull = new_zmq_pull_socket() - double precision, allocatable :: pt2(:,:), norm_pert(:,:), H_pert_diag(:,:) - allocate ( pt2(N_states,2), norm_pert(N_states,2), H_pert_diag(N_states,2)) + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + double precision, allocatable :: pt2_result(:,:), norm_pert_result(:,:), H_pert_diag_result(:,:) + allocate (pt2(N_states), norm_pert(N_states), H_pert_diag(N_states)) + allocate (pt2_result(N_states,N_det_generators), norm_pert_result(N_states,N_det_generators), & + H_pert_diag_result(N_states,N_det_generators)) - pt2 = 0.d0 - norm_pert = 0.d0 - H_pert_diag = 0.d0 + pt2_result = 0.d0 + norm_pert_result = 0.d0 + H_pert_diag_result = 0.d0 accu = 0_8 more = 1 do while (more == 1) - call pull_pt2(zmq_socket_pull, pt2, norm_pert, H_pert_diag, N_states, n, task_id) + call pull_pt2(zmq_socket_pull, pt2, norm_pert, H_pert_diag, i_generator, N_states, n, task_id) if (n > 0) then do k=1,N_states - pt2(k,2) = pt2(k,1) + pt2(k,2) - norm_pert(k,2) = norm_pert(k,1) + norm_pert(k,2) - H_pert_diag(k,2) = H_pert_diag(k,1) + H_pert_diag(k,2) + pt2_result(k,i_generator) = pt2(k) + norm_pert_result(k,i_generator) = norm_pert(k) + H_pert_diag_result(k,i_generator) = H_pert_diag(k) enddo accu = accu + 1_8 call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) @@ -234,9 +244,10 @@ subroutine $subroutine_collector socket_result = new_zmq_pair_socket(.False.) - call push_pt2(socket_result, pt2(1,2), norm_pert(1,2), H_pert_diag(1,2), N_states,0) + call push_pt2(socket_result, pt2_result, norm_pert_result, H_pert_diag_result, i_generator, & + N_states*N_det_generators,0) - deallocate ( pt2, norm_pert, H_pert_diag) + deallocate (pt2, norm_pert, H_pert_diag, pt2_result, norm_pert_result, H_pert_diag_result) call end_zmq_pair_socket(socket_result) diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index deba43c5..7c1f43c2 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -315,6 +315,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun double precision, intent(inout) :: u_in(dim_in,N_st) double precision, intent(out) :: energies(N_st) + integer :: sze_8 integer :: iter integer :: i,j,k,l,m logical :: converged @@ -334,6 +335,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun double precision :: to_print(2,N_st) double precision :: cpu, wall + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, Wt, y, h, lambda call write_time(iunit) @@ -362,12 +364,15 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun enddo write(iunit,'(A)') trim(write_buffer) + integer, external :: align_double + sze_8 = align_double(sze) + allocate( & kl_pairs(2,N_st*(N_st+1)/2), & - W(sze,N_st,davidson_sze_max), & + W(sze_8,N_st,davidson_sze_max), & Wt(sze), & - U(sze,N_st,davidson_sze_max), & - R(sze,N_st), & + U(sze_8,N_st,davidson_sze_max), & + R(sze_8,N_st), & h(N_st,davidson_sze_max,N_st,davidson_sze_max), & y(N_st,davidson_sze_max,N_st,davidson_sze_max), & lambda(N_st*davidson_sze_max)) @@ -381,39 +386,52 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! ============== - k_pairs=0 - do l=1,N_st - do k=1,l - k_pairs+=1 - kl_pairs(1,k_pairs) = k - kl_pairs(2,k_pairs) = l - enddo - enddo - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & - !$OMP Nint,dets_in,u_in) & - !$OMP PRIVATE(k,l,kl,i) - - - ! Orthonormalize initial guess - ! ============================ - - !$OMP DO - do kl=1,k_pairs - k = kl_pairs(1,kl) - l = kl_pairs(2,kl) - if (k/=l) then - overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) - overlap(l,k) = overlap(k,l) - else - overlap(k,k) = u_dot_u(U_in(1,k),sze) - endif - enddo - !$OMP END DO - !$OMP END PARALLEL + if (N_st > 1) then - call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) + k_pairs=0 + do l=1,N_st + do k=1,l + k_pairs+=1 + kl_pairs(1,k_pairs) = k + kl_pairs(2,k_pairs) = l + enddo + enddo + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & + !$OMP Nint,dets_in,u_in) & + !$OMP PRIVATE(k,l,kl) + + + ! Orthonormalize initial guess + ! ============================ + + !$OMP DO + do kl=1,k_pairs + k = kl_pairs(1,kl) + l = kl_pairs(2,kl) + if (k/=l) then + overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) + overlap(l,k) = overlap(k,l) + else + overlap(k,k) = u_dot_u(U_in(1,k),sze) + endif + enddo + !$OMP END DO + !$OMP END PARALLEL + + call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) + + else + + overlap(1,1) = u_dot_u(U_in(1,1),sze) + double precision :: f + f = 1.d0 / dsqrt(overlap(1,1)) + do i=1,sze + U_in(i,1) = U_in(i,1) * f + enddo + + endif ! Davidson iterations ! =================== @@ -473,7 +491,10 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(k,i,l,iter2) SHARED(U,W,R,y,iter,lambda,N_st,sze) do k=1,N_st + !$OMP DO do i=1,sze U(i,k,iter+1) = 0.d0 W(i,k,iter+1) = 0.d0 @@ -484,7 +505,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun enddo enddo enddo + !$OMP END DO enddo + !$OMP END PARALLEL ! Compute residual vector ! ----------------------- diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 4ab9deda..9bcc95f9 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1608,12 +1608,11 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) integer :: i,j,k,l, jj,ii integer :: i0, j0 - integer, allocatable :: shortcut(:), sort_idx(:) - integer(bit_kind), allocatable :: sorted(:,:), version(:,:) + integer, allocatable :: shortcut(:,:), sort_idx(:,:) + integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) integer(bit_kind) :: sorted_i(Nint) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi - double precision :: local_threshold ASSERT (Nint > 0) @@ -1621,104 +1620,83 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) ASSERT (n>0) PROVIDE ref_bitmask_energy davidson_criterion - allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n)) + allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) v_0 = 0.d0 - call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) + call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) + call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) + !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i)& + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version) allocate(vt(n)) Vt = 0.d0 !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0) - do sh2=1,sh + do sh=1,shortcut(0,1) + do sh2=sh,shortcut(0,1) exa = 0 do ni=1,Nint - exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2))) + exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) end do if(exa > 2) then cycle end if - do i=shortcut(sh),shortcut(sh+1)-1 - org_i = sort_idx(i) - local_threshold = threshold_davidson - dabs(u_0(org_i)) + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) if(sh==sh2) then endi = i-1 else - endi = shortcut(sh2+1)-1 + endi = shortcut(sh2+1,1)-1 end if do ni=1,Nint - sorted_i(ni) = sorted(ni,i) + sorted_i(ni) = sorted(ni,i,1) enddo - do j=shortcut(sh2),endi - org_j = sort_idx(j) - if ( dabs(u_0(org_j)) > local_threshold ) then - ext = exa - do ni=1,Nint - ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j))) - end do - if(ext <= 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - vt (org_i) = vt (org_i) + hij*u_0(org_j) - vt (org_j) = vt (org_j) + hij*u_0(org_i) - endif + do j=shortcut(sh2,1),endi + org_j = sort_idx(j,1) + ext = exa + do ni=1,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + end do + if(ext <= 4) then + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + vt (org_i) = vt (org_i) + hij*u_0(org_j) + vt (org_j) = vt (org_j) + hij*u_0(org_i) endif enddo enddo enddo enddo - !$OMP END DO - - !$OMP CRITICAL - do i=1,n - v_0(i) = v_0(i) + vt(i) - enddo - !$OMP END CRITICAL - - deallocate(vt) - !$OMP END PARALLEL - - call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) - allocate(vt(n)) - Vt = 0.d0 + !$OMP END DO NOWAIT !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0) - do i=shortcut(sh),shortcut(sh+1)-1 - org_i = sort_idx(i) - local_threshold = threshold_davidson - dabs(u_0(org_i)) - do j=shortcut(sh),i-1 - org_j = sort_idx(j) - if ( dabs(u_0(org_j)) > local_threshold ) then - ext = 0 - do ni=1,Nint - ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j))) - end do - if(ext == 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - vt (org_i) = vt (org_i) + hij*u_0(org_j) - vt (org_j) = vt (org_j) + hij*u_0(org_i) - end if + do sh=1,shortcut(0,2) + do i=shortcut(sh,2),shortcut(sh+1,2)-1 + org_i = sort_idx(i,2) + do j=shortcut(sh,2),i-1 + org_j = sort_idx(j,2) + ext = 0 + do ni=1,Nint + ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) + end do + if(ext == 4) then + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + vt (org_i) = vt (org_i) + hij*u_0(org_j) + vt (org_j) = vt (org_j) + hij*u_0(org_i) end if end do end do enddo - !$OMP END DO + !$OMP END DO NOWAIT !$OMP CRITICAL - do i=1,n + do i=n,1,-1 v_0(i) = v_0(i) + vt(i) enddo !$OMP END CRITICAL + deallocate(vt) !$OMP END PARALLEL diff --git a/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES index 6aace6f2..152711f3 100644 --- a/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES +++ b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Pseudo Bitmask ZMQ Integrals_Monoelec +Pseudo Bitmask ZMQ diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index b6446a40..b7c75fb8 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -289,18 +289,11 @@ 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) @@ -375,136 +368,25 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer(ZMQ_PTR) :: zmq_to_qp_run_socket character*(32) :: task - if (.not.has_libint) then + call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') - call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') + do l=1,ao_num + write(task,*) l + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + enddo - do l=1,ao_num - write(task,*) l - call add_task_to_taskserver(zmq_to_qp_run_socket,task) - enddo - - integer(ZMQ_PTR) :: collector_thread - external :: ao_bielec_integrals_in_map_collector - rc = pthread_create(collector_thread, ao_bielec_integrals_in_map_collector) - - !$OMP PARALLEL DEFAULT(private) - !$OMP TASK PRIVATE(i) - i = omp_get_thread_num() + PROVIDE nproc + !$OMP PARALLEL DEFAULT(private) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_bielec_integrals_in_map_collector(i) + else call ao_bielec_integrals_in_map_slave_inproc(i) - !$OMP END TASK - !$OMP TASKWAIT - !$OMP END PARALLEL + endif + !$OMP END PARALLEL - rc = pthread_join(collector_thread) + call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals') - call end_parallel_job(zmq_to_qp_run_socket, 'ao_integrals') - - else - - double precision, allocatable :: buffer_int(:) - - PROVIDE has_libint - - 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, shell_num - - print*, s1, "/", shell_num - - bf1_begin = shell_idx(1,s1) - bf1_end = shell_idx(2,s1) - n1 = 1 + bf1_end - bf1_begin - - do s2 = 1, shell_num - - bf2_begin = shell_idx(1,s2) - bf2_end = shell_idx(2,s2) - n2 = 1 + bf2_end - bf2_begin - - do s3 = 1, shell_num - - bf3_begin = shell_idx(1,s3) - bf3_end = shell_idx(2,s3) - n3 = 1 + bf3_end - bf3_begin - - do s4 = 1, shell_num - - bf4_begin = shell_idx(1,s4) - bf4_end = shell_idx(2,s4) - n4 = 1 + bf4_end - bf4_begin - - ! ========================== ! - ! Compute the shell integral ! - ! ========================== ! - integer :: sze - sze = n1*n2*n3*n4 - allocate(buffer_int(sze)) - allocate(buffer_i(sze)) - allocate(buffer_value(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 - -! if (bf1 > bf3) cycle -! if (bf2 > bf4) cycle -! if (bf1 > bf2) cycle - - !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 - - n_integrals += 1 - buffer_value(n_integrals) = buffer_int(f1234) * norm - call bielec_integrals_index(bf1,bf2,bf3,bf4,buffer_i(n_integrals)) - - enddo - enddo - enddo - enddo - - !Deallocate the buffer_intergral for the shell - deallocate(buffer_int, buffer_i, buffer_value) - if (n_integrals >= 0) then - call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) - endif - - enddo - enddo - enddo - enddo - - endif print*, 'Sorting the map' call map_sort(ao_integrals_map) diff --git a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f index 6102d119..f15376b5 100644 --- a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f +++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f @@ -67,6 +67,8 @@ end + + subroutine ao_bielec_integrals_in_map_slave(thread,iproc) use map_module use f77_zmq @@ -107,7 +109,7 @@ subroutine ao_bielec_integrals_in_map_slave(thread,iproc) call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, 0) enddo call compute_ao_integrals_jl(l,l,n_integrals,buffer_i,buffer_value) - call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_integrals) + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id) enddo @@ -127,7 +129,7 @@ subroutine pull_integrals(zmq_socket_pull, n_integrals, buffer_i, buffer_value, BEGIN_DOC ! How the collector pulls the computed integrals END_DOC - integer(ZMQ_PTR), intent(out) :: zmq_socket_pull + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull integer, intent(out) :: n_integrals integer(key_kind), intent(out) :: buffer_i(*) real(integral_kind), intent(out) :: buffer_value(*) diff --git a/src/Integrals_Bielec/integral_bielec.cc b/src/Integrals_Bielec/integral_bielec.cc deleted file mode 100644 index e63345f9..00000000 --- a/src/Integrals_Bielec/integral_bielec.cc +++ /dev/null @@ -1,247 +0,0 @@ -/* - * 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 "C++11 support is required" -#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 deleted file mode 100644 index 2cedf283..00000000 --- a/src/Integrals_Bielec/libint_module.f90 +++ /dev/null @@ -1,51 +0,0 @@ -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 diff --git a/src/Integrals_Bielec/q_package.ezfio_config b/src/Integrals_Bielec/q_package.ezfio_config index 35cba011..8c64c790 100644 --- a/src/Integrals_Bielec/q_package.ezfio_config +++ b/src/Integrals_Bielec/q_package.ezfio_config @@ -4,5 +4,3 @@ work save empty logical -libint - empty logical diff --git a/src/Integrals_Bielec/qp_ao_ints.irp.f b/src/Integrals_Bielec/qp_ao_ints.irp.f index f932df0f..c60b4e5d 100644 --- a/src/Integrals_Bielec/qp_ao_ints.irp.f +++ b/src/Integrals_Bielec/qp_ao_ints.irp.f @@ -8,7 +8,8 @@ program qp_ao_ints call switch_qp_run_to_master - PROVIDE zmq_context + zmq_context = f77_zmq_ctx_new () + ! Set the state of the ZMQ zmq_state = 'ao_integrals' diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 615ed127..789bc9ea 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -3,10 +3,14 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] BEGIN_DOC ! Pseudo-potential END_DOC + ao_pseudo_integral = 0.d0 if (do_pseudo) then - ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local - else - ao_pseudo_integral = 0.d0 + if (pseudo_klocmax > 0) then + ao_pseudo_integral += ao_pseudo_integral_local + endif + if (pseudo_kmax > 0) then + ao_pseudo_integral += ao_pseudo_integral_non_local + endif endif END_PROVIDER diff --git a/src/Integrals_Monoelec/shell.irp.f b/src/Integrals_Monoelec/shell.irp.f deleted file mode 100644 index 6972bbc3..00000000 --- a/src/Integrals_Monoelec/shell.irp.f +++ /dev/null @@ -1,45 +0,0 @@ -use libint_module -BEGIN_PROVIDER [ logical, has_libint ] - implicit none - BEGIN_DOC - ! If true, use libint - END_DOC - character :: y - call getenv('QP_LIBINT', y) - if (y=='1') then - has_libint = .True. - call init_libint(trim(ezfio_filename)//char((0))) - else - PROVIDE ezfio_filename - has_libint = .False. - endif - -END_PROVIDER - -BEGIN_PROVIDER [ integer, shell_num ] - implicit none - BEGIN_DOC - ! Number of shells - END_DOC - if (has_libint) then - shell_num = get_nb_shell() - else - stop 'shell_num not implemented without libint' - endif - -END_PROVIDER - - -BEGIN_PROVIDER [ integer, shell_idx, (2,shell_num) ] - implicit none - BEGIN_DOC - ! Contains the 1st and last AO index in each shell - END_DOC - if (has_libint) then - call map_shell_to_basis_function_interval(2*shell_num,shell_idx) - else - stop 'shell_idx not implemented without libint' - endif -END_PROVIDER - - diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 13138499..11de1270 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -73,6 +73,10 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j + if (n < 2) then + return + endif + allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) @@ -151,6 +155,10 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j, k + if (n < 2) then + return + endif + call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) !$OMP PARALLEL DEFAULT(NONE) & diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index a5904183..91a61a43 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -295,18 +295,6 @@ BEGIN_PROVIDER [ integer, nproc ] !$OMP END PARALLEL END_PROVIDER -BEGIN_PROVIDER [ integer, iproc_save, (nproc) ] - implicit none - BEGIN_DOC - ! iproc_save(i) = i-1. Used to start threads with pthreads. - END_DOC - integer :: i - do i=1,nproc - iproc_save(i) = i-1 - enddo - -END_PROVIDER - double precision function u_dot_v(u,v,sze) implicit none @@ -324,6 +312,7 @@ double precision function u_dot_v(u,v,sze) t3 = t2+t2 t4 = t3+t2 u_dot_v = 0.d0 + !DIR$ VECTOR ALWAYS do i=1,t2 u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) + & u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i) @@ -359,6 +348,7 @@ double precision function u_dot_u(u,sze) ! u_dot_u = u_dot_u+u(i)*u(i) ! enddo + !DIR$ VECTOR ALWAYS do i=1,sze u_dot_u = u_dot_u + u(i)*u(i) enddo diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f index d730f612..0bbba56e 100644 --- a/src/ZMQ/utils.irp.f +++ b/src/ZMQ/utils.irp.f @@ -181,14 +181,14 @@ function new_zmq_pair_socket(bind) endif endif - rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 0, 4) + rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 1, 4) if (rc /= 0) then - stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 0, 4)' + stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_SNDHWM, 1, 4)' endif - rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 0, 4) + rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 1, 4) if (rc /= 0) then - stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 0, 4)' + stop 'f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_RCVHWM, 1, 4)' endif rc = f77_zmq_setsockopt(new_zmq_pair_socket, ZMQ_IMMEDIATE, 1, 4) @@ -229,16 +229,11 @@ function new_zmq_pull_socket() stop 'Unable to set ZMQ_LINGER on pull socket' endif - rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_RCVHWM,100000,4) + rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_RCVHWM,1000,4) if (rc /= 0) then stop 'Unable to set ZMQ_RCVHWM on pull socket' endif - rc = f77_zmq_setsockopt(new_zmq_pull_socket,ZMQ_IMMEDIATE,1,4) - if (rc /= 0) then - stop 'Unable to set ZMQ_IMMEDIATE on pull socket' - endif - rc = f77_zmq_bind(new_zmq_pull_socket, zmq_socket_pull_tcp_address) if (rc /= 0) then print *, 'Unable to bind new_zmq_pull_socket (tcp)', zmq_socket_pull_tcp_address @@ -279,7 +274,7 @@ function new_zmq_push_socket(thread) stop 'Unable to set ZMQ_LINGER on push socket' endif - rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,100,4) + rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,1000,4) if (rc /= 0) then stop 'Unable to set ZMQ_SNDHWM on push socket' endif @@ -361,6 +356,8 @@ subroutine end_zmq_pull_socket(zmq_socket_pull) stop 'error' endif + call sleep(1) ! see https://github.com/zeromq/libzmq/issues/1922 + rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4) if (rc /= 0) then stop 'Unable to set ZMQ_LINGER on zmq_socket_pull' diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 1ced9e1d..78ed973d 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -155,7 +155,7 @@ function run_all_1h_1p() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2284994316618 1.e-4 + eq $energy -76.2288648023833 1.e-4 } @@ -166,7 +166,7 @@ function run_all_1h_1p() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483708495521 + run_HF h2o_pseudo.ezfio -16.9483703905461 } @test "FCI H2O VDZ pseudo" {