From 1d58caea77f2630b002efa751d7ac08dcb8c6a18 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 30 Nov 2015 14:54:37 +0100 Subject: [PATCH 1/6] Corrected bug in CISD selected --- plugins/CISD_selected/cisd_selection.irp.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/plugins/CISD_selected/cisd_selection.irp.f b/plugins/CISD_selected/cisd_selection.irp.f index b2178860..ad31269c 100644 --- a/plugins/CISD_selected/cisd_selection.irp.f +++ b/plugins/CISD_selected/cisd_selection.irp.f @@ -41,8 +41,8 @@ program cisd N_det = min(N_det,N_det_max) touch N_det psi_det psi_coef call diagonalize_CI - deallocate(pt2,norm_pert,H_pert_diag) - call save_wavefunction + call save_wavefunction call ezfio_set_cisd_selected_energy(CI_energy) call ezfio_set_cisd_selected_energy_pt2(CI_energy+pt2) + deallocate(pt2,norm_pert,H_pert_diag) end From 6d05c72143181692b1c55823dcdab21ee9c17ed3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 30 Nov 2015 19:36:18 +0100 Subject: [PATCH 2/6] Added CAS-SCF module --- plugins/CASSCF/EZFIO.cfg | 10 + plugins/CASSCF/H_apply.irp.f | 39 ++++ plugins/CASSCF/NEEDED_CHILDREN_MODULES | 1 + plugins/CASSCF/README.rst | 20 ++ plugins/CASSCF/casscf.irp.f | 220 +++++++++++++++++++++ src/Bitmask/bitmasks.irp.f | 8 +- src/Integrals_Bielec/map_integrals.irp.f | 9 - src/Integrals_Bielec/mo_bi_integrals.irp.f | 11 ++ src/MO_Basis/utils.irp.f | 1 - 9 files changed, 302 insertions(+), 17 deletions(-) create mode 100644 plugins/CASSCF/EZFIO.cfg create mode 100644 plugins/CASSCF/H_apply.irp.f create mode 100644 plugins/CASSCF/NEEDED_CHILDREN_MODULES create mode 100644 plugins/CASSCF/README.rst create mode 100644 plugins/CASSCF/casscf.irp.f diff --git a/plugins/CASSCF/EZFIO.cfg b/plugins/CASSCF/EZFIO.cfg new file mode 100644 index 00000000..e9e6e92e --- /dev/null +++ b/plugins/CASSCF/EZFIO.cfg @@ -0,0 +1,10 @@ +[energy] +type: double precision +doc: "Calculated CAS-SCF energy" +interface: ezfio + +[energy_pt2] +type: double precision +doc: "Calculated selected CAS-SCF energy with PT2 correction" +interface: ezfio + diff --git a/plugins/CASSCF/H_apply.irp.f b/plugins/CASSCF/H_apply.irp.f new file mode 100644 index 00000000..35c45fb6 --- /dev/null +++ b/plugins/CASSCF/H_apply.irp.f @@ -0,0 +1,39 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("CAS_SD") +print s + +s = H_apply("CAS_SD_selected_no_skip") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +print s + +s = H_apply("CAS_SD_selected") +s.set_selection_pt2("epstein_nesbet_2x2") +print s + +s = H_apply("CAS_SD_PT2") +s.set_perturbation("epstein_nesbet_2x2") +print s + + +s = H_apply("CAS_S",do_double_exc=False) +print s + +s = H_apply("CAS_S_selected_no_skip",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +print s + +s = H_apply("CAS_S_selected",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +print s + +s = H_apply("CAS_S_PT2",do_double_exc=False) +s.set_perturbation("epstein_nesbet_2x2") +print s + +END_SHELL + diff --git a/plugins/CASSCF/NEEDED_CHILDREN_MODULES b/plugins/CASSCF/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..29e39f2f --- /dev/null +++ b/plugins/CASSCF/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Generators_CAS Perturbation Selectors_full diff --git a/plugins/CASSCF/README.rst b/plugins/CASSCF/README.rst new file mode 100644 index 00000000..ceeb7477 --- /dev/null +++ b/plugins/CASSCF/README.rst @@ -0,0 +1,20 @@ +====== +CASSCF +====== + +This module is not a "real" CAS-SCF. It is an orbital optimization step done by : + +1) Doing the CAS+SD +2) Taking one-electron density matrix +3) Cancelling all active-active rotations +4) Finding the order which matches with the input MOs + + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/CASSCF/casscf.irp.f b/plugins/CASSCF/casscf.irp.f new file mode 100644 index 00000000..4e7450dc --- /dev/null +++ b/plugins/CASSCF/casscf.irp.f @@ -0,0 +1,220 @@ +program casscf + implicit none + BEGIN_DOC +! Optimize MOs and CI coefficients of the CAS + END_DOC + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer(bit_kind), allocatable :: generators_bitmask_save(:,:,:,:) + + integer :: degree, N_generators_bitmask_save, N_det_ci + double precision :: E_old, E_CI + double precision :: selection_criterion_save, selection_criterion_min_save + + integer :: N_det_old + integer :: i, j, k, l + integer :: i_bit, j_bit, i_int, j_int + integer(bit_kind), allocatable :: bit_tmp(:), cas_bm(:) + character*(64) :: label + + allocate( pt2(N_states), norm_pert(N_states),H_pert_diag(N_states) ) + allocate( generators_bitmask_save(N_int,2,6,N_generators_bitmask) ) + allocate( bit_tmp(N_int), cas_bm(N_int) ) + + PROVIDE N_det_cas + N_det_old = 0 + pt2 = 1.d0 + E_CI = 1.d0 + E_old = 0.d0 + diag_algorithm = "Lapack" + selection_criterion_save = selection_criterion + selection_criterion_min_save = selection_criterion_min + + + cas_bm = 0_bit_kind + do i=1,N_cas_bitmask + do j=1,N_int + cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,1,i)) + cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,2,i)) + enddo + enddo + + ! Save CAS-SD bitmask + generators_bitmask_save = generators_bitmask + N_generators_bitmask_save = N_generators_bitmask + + ! Set the CAS bitmask + do i=1,6 + generators_bitmask(:,:,i,:) = cas_bitmask + enddo + N_generators_bitmask = N_cas_bitmask + SOFT_TOUCH generators_bitmask N_generators_bitmask + + + ! If the number of dets already in the file is larger than the requested + ! number of determinants, truncate the wf + if (N_det > N_det_max) then + call diagonalize_CI + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = N_det_max + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + endif + + ! Start MCSCF iteration + + ! CAS-CI + ! ------ + + E_old = E_CI + + ! Reset the selection criterion + selection_criterion = selection_criterion_save + selection_criterion_min = selection_criterion_min_save + SOFT_TOUCH selection_criterion_min selection_criterion selection_criterion_factor + + ! Set the CAS bitmask + do i=1,6 + generators_bitmask(:,:,i,:) = cas_bitmask + enddo + N_generators_bitmask = N_cas_bitmask + SOFT_TOUCH generators_bitmask N_generators_bitmask + + do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_states))) > pt2_max) + N_det_old = N_det + call H_apply_CAS_SD_selected_no_skip(pt2, norm_pert, H_pert_diag, N_states) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + if (N_det > N_det_max) then + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = N_det_max + soft_touch N_det psi_det psi_coef + endif + call diagonalize_CI + call save_wavefunction + print *, '======' + print *, 'CAS-CI' + print *, '======' + print *, '' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E(CAS) = ', CI_energy + print *, 'E(CAS)+PT2 = ', CI_energy+pt2 + print *, '-----' + print *, '' + E_CI = sum(CI_energy(1:N_states)+pt2(1:N_states))/dble(N_states) + + call ezfio_set_casscf_energy(CI_energy(1)) + if (abort_all) then + exit + endif + if (N_det == N_det_old) then + exit + endif + + enddo + + ! Super-CI + ! -------- + + selection_criterion_min = 1.d-12 + selection_criterion = 1.d-12 + + ! Set the CAS bitmask + generators_bitmask = generators_bitmask_save + N_generators_bitmask = N_generators_bitmask_save + SOFT_TOUCH generators_bitmask N_generators_bitmask selection_criterion selection_criterion_min selection_criterion_factor + + N_det_ci = N_det + + call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_states) + + + do i=1,mo_tot_num + i_int = ishft(i-1,-bit_kind_shift)+1 + i_bit = j-ishft(i_int-1,bit_kind_shift)-1 + bit_tmp(:) = 0_bit_kind + bit_tmp(i_int) = ibset(0_bit_kind,i_bit) + if (iand(bit_tmp(i_int), cas_bm(i_int)) == 0_bit_kind) then + ! Not a CAS MO + cycle + endif + + do j=1,mo_tot_num + if (j == i) then + cycle + endif + j_int = ishft(j-1,-bit_kind_shift)+1 + j_bit = j-ishft(j_int-1,bit_kind_shift)-1 + bit_tmp(:) = 0_bit_kind + bit_tmp(j_int) = ibset(0_bit_kind,j_bit) + if (iand(bit_tmp(j_int), cas_bm(j_int)) == 0_bit_kind) then + ! Not a CAS MO + cycle + endif + ! Now, both i and j are MOs of the CAS. De-couple them in the DM + one_body_dm_mo(i,j) = 0.d0 + enddo + + enddo + + SOFT_TOUCH one_body_dm_mo + + double precision :: mx, ov + double precision, allocatable :: mo_coef_old(:,:) + integer, allocatable :: iorder(:) + logical, allocatable :: selected(:) + allocate( mo_coef_old(size(mo_coef,1), size(mo_coef,2)), iorder(mo_tot_num), selected(mo_tot_num) ) + mo_coef_old = mo_coef + label = "Canonical" + call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,-1) + selected = .False. + do j=1,mo_tot_num + mx = -1.d0 + iorder(j) = j + do i=1,mo_tot_num + if (selected(i)) then + cycle + endif + ov = 0.d0 + do l=1,ao_num + do k=1,ao_num + ov = ov + mo_coef_old(k,j) * ao_overlap(k,l) * mo_coef(l,i) + enddo + enddo + ov= dabs(ov) + if (ov > mx) then + mx = ov + iorder(j) = i + endif + enddo + selected( iorder(j) ) = .True. + enddo + mo_coef_old = mo_coef + do i=1,mo_tot_num + mo_coef(:,i) = mo_coef_old(:,iorder(i)) + enddo + + call save_mos + + call write_double(6,E_CI,"Energy(CAS)") + + deallocate( mo_coef_old ) + deallocate( pt2, norm_pert,H_pert_diag ) + deallocate( generators_bitmask_save ) + deallocate( bit_tmp, cas_bm, iorder ) +end diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 044fa18b..29588369 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -262,13 +262,7 @@ END_PROVIDER logical :: exists integer :: j,i integer :: i_hole,i_part,i_gen - PROVIDE ezfio_filename -!do j = 1, N_int -! inact_bitmask(j,1) = xor(generators_bitmask(j,1,1,1),cas_bitmask(j,1,1)) -! inact_bitmask(j,2) = xor(generators_bitmask(j,2,1,1),cas_bitmask(j,2,1)) -! virt_bitmask(j,1) = xor(generators_bitmask(j,1,2,1),cas_bitmask(j,1,1)) -! virt_bitmask(j,2) = xor(generators_bitmask(j,2,2,1),cas_bitmask(j,2,1)) -!enddo + n_inact_orb = 0 n_virt_orb = 0 if(N_generators_bitmask == 1)then diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 3b737f5d..43c207d5 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -424,15 +424,6 @@ integer*8 function get_mo_map_size() get_mo_map_size = mo_integrals_map % n_elements end -subroutine clear_mo_map - implicit none - BEGIN_DOC - ! Frees the memory of the MO map - END_DOC - call map_deinit(mo_integrals_map) - FREE mo_integrals_map -end - BEGIN_TEMPLATE subroutine dump_$ao_integrals(filename) diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 1fa303b5..83f0ce05 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -503,3 +503,14 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to enddo END_PROVIDER + + +subroutine clear_mo_map + implicit none + BEGIN_DOC + ! Frees the memory of the MO map + END_DOC + call map_deinit(mo_integrals_map) + FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti + FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map +end diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 3f4a260a..0d8ef5fa 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -98,7 +98,6 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) call write_time(output_mo_basis) mo_label = label - SOFT_TOUCH mo_coef mo_label end subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label) From 4f3c07f54e4837f735036112c488ed98c197ad96 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 30 Nov 2015 20:57:41 +0100 Subject: [PATCH 3/6] LinearAlgebra --- plugins/QmcChem/save_for_qmcchem.irp.f | 1 + src/Utils/LinearAlgebra.irp.f | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/plugins/QmcChem/save_for_qmcchem.irp.f b/plugins/QmcChem/save_for_qmcchem.irp.f index 4b028a7c..c8ddb4d9 100644 --- a/plugins/QmcChem/save_for_qmcchem.irp.f +++ b/plugins/QmcChem/save_for_qmcchem.irp.f @@ -1,6 +1,7 @@ program save_for_qmc read_wf = .True. TOUCH read_wf + print *, "N_det = ", N_det call write_spindeterminants if (do_pseudo) then call write_pseudopotential diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index dab9e921..e3ef0bfe 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -83,7 +83,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do i=1,n if ( D(i) < 1.d-12 ) then - stop 'Linear dependence in basis set' + D(i) = 0.d0 else D(i) = 1.d0/dsqrt(D(i)) endif From 665638ebc0f18d8aae65cb5980149e356e11bb91 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Dec 2015 12:34:24 +0100 Subject: [PATCH 4/6] configure OK with ZeroMQ --- install/scripts/install_f77zmq.sh | 1 + scripts/compilation/qp_create_ninja.py | 6 +++--- scripts/module/module_handler.py | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/install/scripts/install_f77zmq.sh b/install/scripts/install_f77zmq.sh index ff2417fb..8357857c 100755 --- a/install/scripts/install_f77zmq.sh +++ b/install/scripts/install_f77zmq.sh @@ -15,6 +15,7 @@ function _install() make -j 8 || exit 1 mv libf77zmq.a "${QP_ROOT}"/lib || exit 1 mv libf77zmq.so "${QP_ROOT}"/lib || exit 1 + cp f77_zmq.h "${QP_ROOT}"/src/ZMQ/ cd - return 0 } diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 7171c2df..a57c7cbf 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -36,8 +36,8 @@ except ImportError: from qp_path import QP_ROOT, QP_SRC, QP_EZFIO LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") -EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") -ZMQ_LIB = join(QP_ROOT, "lib", "libzmq.a") + " " + join(QP_ROOT, "lib", "libf77zmq.a") +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" ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja") header = r"""# @@ -262,7 +262,7 @@ def ninja_ezfio_rule(): l_flag = ["export {0}='${0}'".format(flag) for flag in ["FC", "FCFLAGS", "IRPF90"]] - install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio.a") + install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio_irp.a") l_cmd = ["cd {0}".format(QP_EZFIO)] + l_flag l_cmd += ["rm -f make.config ; ninja && ln -sf {0} {1}".format(install_lib_ezfio, EZFIO_LIB)] diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py index ebbc367a..136dc8cf 100755 --- a/scripts/module/module_handler.py +++ b/scripts/module/module_handler.py @@ -88,7 +88,7 @@ def get_l_module_descendant(d_child, l_module): except KeyError: print >> sys.stderr, "Error: " print >> sys.stderr, "`{0}` is not a submodule".format(module) - print >> sys.stderr, "Check the typo (orthograph, case, '/', etc.) " + print >> sys.stderr, "Check the typo (spelling, case, '/', etc.) " sys.exit(1) return list(set(l)) From 94cf7c2b321e436158a003bf4dc06bf776cc44e0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Dec 2015 16:32:01 +0100 Subject: [PATCH 5/6] AO integrals with ZeroMQ --- plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES | 2 +- src/Integrals_Bielec/NEEDED_CHILDREN_MODULES | 2 +- src/Integrals_Bielec/ao_bi_integrals.irp.f | 219 +++++++++---------- src/Integrals_Bielec/map_integrals.irp.f | 3 +- src/ZMQ/NEEDED_CHILDREN_MODULES | 1 + src/ZMQ/README.rst | 15 ++ src/ZMQ/f77_zmq_module.f90 | 4 + src/ZMQ/zmq.irp.f | 105 +++++++++ 8 files changed, 235 insertions(+), 116 deletions(-) create mode 100644 src/ZMQ/NEEDED_CHILDREN_MODULES create mode 100644 src/ZMQ/README.rst create mode 100644 src/ZMQ/f77_zmq_module.f90 create mode 100644 src/ZMQ/zmq.irp.f diff --git a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES index 784cb0fb..85bdd3ad 100644 --- a/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES +++ b/plugins/Hartree_Fock/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Integrals_Bielec MOGuess +Integrals_Bielec MOGuess diff --git a/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES index 5888fc95..152711f3 100644 --- a/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES +++ b/src/Integrals_Bielec/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Pseudo Bitmask +Pseudo Bitmask ZMQ diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 34833f3b..53ce68e9 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -301,7 +301,7 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) double precision :: thresh thresh = ao_integrals_threshold - integer :: n_centers, i + integer :: i if (ao_overlap_abs(j,l) < thresh) then buffer_value = 0._integral_kind @@ -329,6 +329,7 @@ end BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] implicit none + use f77_zmq use map_module BEGIN_DOC ! Map of Atomic integrals @@ -345,9 +346,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer(key_kind),allocatable :: buffer_i(:) integer,parameter :: size_buffer = 1024*64 real(integral_kind),allocatable :: buffer_value(:) - integer(omp_lock_kind) :: lock - integer :: n_integrals, n_centers, thread_num + integer :: n_integrals, rc integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax integral = ao_bielec_integral(1,1,1,1) @@ -363,120 +363,61 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] endif endif - kk=1 - do l = 1, ao_num ! r2 - do j = 1, l ! r2 - jl_pairs(1,kk) = j - jl_pairs(2,kk) = l - kk += 1 - enddo - enddo - - PROVIDE progress_bar - call omp_init_lock(lock) - lmax = ao_num*(ao_num+1)/2 print*, 'Providing the AO integrals' call wall_time(wall_0) call wall_time(wall_1) call cpu_time(cpu_1) - call start_progress(lmax,'AO integrals (MB)',0.d0) - !$OMP PARALLEL PRIVATE(i,j,k,l,kk, & - !$OMP integral,buffer_i,buffer_value,n_integrals, & - !$OMP cpu_2,wall_2,i1,j1,thread_num) & - !$OMP DEFAULT(NONE) & - !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & - !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & - !$OMP ao_overlap_abs,ao_overlap,abort_here, & - !$OMP wall_0,progress_bar,progress_value, & - !$OMP ao_bielec_integral_schwartz) - - allocate(buffer_i(size_buffer)) - allocate(buffer_value(size_buffer)) - n_integrals = 0 -!$ thread_num = omp_get_thread_num() - - !$OMP DO SCHEDULE(dynamic) - do kk=1,lmax -IRP_IF COARRAY - if (mod(kk-this_image(),num_images()) /= 0) then - cycle - endif -IRP_ENDIF - if (abort_here) then - cycle - endif - if (thread_num == 0) then - progress_bar(1) = kk - endif - j = jl_pairs(1,kk) - l = jl_pairs(2,kk) - j1 = j+ishft(l*l-l,-1) - if (ao_overlap_abs(j,l) < thresh) then - cycle - endif - do k = 1, ao_num ! r1 - i1 = ishft(k*k-k,-1) - if (i1 > j1) then - exit - endif - do i = 1, k - i1 += 1 - if (i1 > j1) then - exit - endif - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then - cycle - endif - if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then - cycle - endif - !DIR$ FORCEINLINE - integral = ao_bielec_integral(i,k,j,l) - if (abs(integral) < thresh) then - cycle - endif - n_integrals += 1 - !DIR$ FORCEINLINE - call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) - buffer_value(n_integrals) = integral - if (n_integrals > 1024 ) then - if (omp_test_lock(lock)) then - call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) - n_integrals = 0 - call omp_unset_lock(lock) - endif - endif - if (n_integrals == size(buffer_i)) then - call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) - n_integrals = 0 - endif - enddo - enddo - call wall_time(wall_2) - - if (thread_num == 0) then - if (wall_2 - wall_0 > 1.d0) then - wall_0 = wall_2 - print*, 100.*float(kk)/float(lmax), '% in ', & - wall_2-wall_1, 's', map_mb(ao_integrals_map) ,'MB' - progress_value = dble(map_mb(ao_integrals_map)) - endif - endif - enddo - !$OMP END DO NOWAIT - call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) - deallocate(buffer_i) - deallocate(buffer_value) - !$OMP END PARALLEL - call omp_destroy_lock(lock) - call stop_progress - if (abort_here) then - stop 'Aborting in AO integrals calculation' + + integer(ZMQ_PTR) :: zmq_socket_rep_inproc, zmq_socket_push_inproc + zmq_socket_rep_inproc = f77_zmq_socket(zmq_context, ZMQ_REP) + rc = f77_zmq_bind(zmq_socket_rep_inproc, 'inproc://req_rep') + if (rc /= 0) then + stop 'Unable to connect zmq_socket_rep_inproc' endif -IRP_IF COARRAY - print*, 'Communicating the map' - call communicate_ao_integrals() -IRP_ENDIF COARRAY + + integer(ZMQ_PTR) :: thread(0:nproc) + external :: ao_bielec_integrals_in_map_slave, ao_bielec_integrals_in_map_collector + rc = pthread_create( thread(0), ao_bielec_integrals_in_map_collector ) + ! Create client threads + do i=1,nproc + rc = pthread_create( thread(i), ao_bielec_integrals_in_map_slave ) + enddo + + character*(64) :: message_string + + do l = ao_num, 1, -1 + rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0) + print *, l + ! TODO : error handling + ASSERT (rc >= 0) + ASSERT (message == 'get_ao_integrals') + rc = f77_zmq_send( zmq_socket_rep_inproc, l, 4, 0) + enddo + do i=1,nproc + rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0) + ! TODO : error handling + ASSERT (rc >= 0) + ASSERT (message == 'get_ao_integrals') + rc = f77_zmq_send( zmq_socket_rep_inproc, 0, 4, 0) + enddo + ! TODO terminer thread(0) + + rc = f77_zmq_unbind(zmq_socket_rep_inproc, 'inproc://req_rep') + do i=1,nproc + rc = pthread_join( thread(i) ) + enddo + + zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH) + rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull') + if (rc /= 0) then + stop 'Unable to connect zmq_socket_push_inproc' + endif + rc = f77_zmq_send( zmq_socket_push_inproc, -1, 4, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push_inproc, 0_key_kind, key_kind, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push_inproc, 0_integral_kind, integral_kind, 0) + + rc = pthread_join( thread(0) ) + print*, 'Sorting the map' call map_sort(ao_integrals_map) call cpu_time(cpu_2) @@ -1256,3 +1197,57 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) end + + +subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) + implicit none + use map_module + BEGIN_DOC + ! Parallel client for AO integrals + END_DOC + + integer, intent(in) :: j,l + integer,intent(out) :: n_integrals + integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) + real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num) + + integer :: i,k + double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + double precision :: thresh + integer :: kk, m, j1, i1 + + thresh = ao_integrals_threshold + + n_integrals = 0 + + j1 = j+ishft(l*l-l,-1) + do k = 1, ao_num ! r1 + i1 = ishft(k*k-k,-1) + if (i1 > j1) then + exit + endif + do i = 1, k + i1 += 1 + if (i1 > j1) then + exit + endif + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then + cycle + endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + cycle + endif + !DIR$ FORCEINLINE + integral = ao_bielec_integral(i,k,j,l) + if (abs(integral) < thresh) then + cycle + endif + n_integrals += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + buffer_value(n_integrals) = integral + enddo + enddo + +end diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 43c207d5..84b08715 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -247,8 +247,7 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ] print*, 'MO map initialized' END_PROVIDER -subroutine insert_into_ao_integrals_map(n_integrals, & - buffer_i, buffer_values) +subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) use map_module implicit none BEGIN_DOC diff --git a/src/ZMQ/NEEDED_CHILDREN_MODULES b/src/ZMQ/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/src/ZMQ/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ + diff --git a/src/ZMQ/README.rst b/src/ZMQ/README.rst new file mode 100644 index 00000000..9a12751d --- /dev/null +++ b/src/ZMQ/README.rst @@ -0,0 +1,15 @@ +=== +ZMQ +=== + +Socket address : defined as an environment variable : QP_RUN_ADDRESS + + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/src/ZMQ/f77_zmq_module.f90 b/src/ZMQ/f77_zmq_module.f90 new file mode 100644 index 00000000..d0f551fa --- /dev/null +++ b/src/ZMQ/f77_zmq_module.f90 @@ -0,0 +1,4 @@ +module f77_zmq + include 'f77_zmq.h' +end module + diff --git a/src/ZMQ/zmq.irp.f b/src/ZMQ/zmq.irp.f new file mode 100644 index 00000000..1577e12f --- /dev/null +++ b/src/ZMQ/zmq.irp.f @@ -0,0 +1,105 @@ +use f77_zmq + + +BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ] + implicit none + BEGIN_DOC + ! Context for the ZeroMQ library + END_DOC + zmq_context = f77_zmq_ctx_new () +END_PROVIDER + + + BEGIN_PROVIDER [ character*(128), qp_run_address ] +&BEGIN_PROVIDER [ integer, zmq_port_start ] + implicit none + BEGIN_DOC + ! Address of the qp_run socket + ! Example : tcp://130.120.229.139:12345 + END_DOC + character*(128) :: buffer + call getenv('QP_RUN_ADDRESS',buffer) + if (trim(buffer) == '') then + stop 'QP_RUN_ADDRESS environment variable not defined' + endif + + print *, trim(buffer) + integer :: i + do i=len(buffer),1,-1 + if ( buffer(i:i) == ':') then + qp_run_address = trim(buffer(1:i-1)) + read(buffer(i+1:), *) zmq_port_start + exit + endif + enddo +END_PROVIDER + + +function zmq_port(ishift) + implicit none + integer, intent(in) :: ishift + character*(8) :: zmq_port + write(zmq_port,'(I8)') zmq_port_start+ishift + zmq_port = adjustl(trim(zmq_port)) +end + + +BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_to_qp_run_socket ] + implicit none + BEGIN_DOC + ! Socket on which the qp_run process replies + END_DOC + integer :: rc + zmq_to_qp_run_socket = f77_zmq_socket(zmq_context, ZMQ_REQ) + rc = f77_zmq_connect(zmq_to_qp_run_socket, trim(qp_run_address)) + if (rc /= 0) then + stop 'Unable to connect zmq_to_qp_run_socket' + endif + integer :: i + i=4 + rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, i) + if (rc /= 0) then + stop 'Unable to set send timout in zmq_to_qp_run_socket' + endif + rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, i) + if (rc /= 0) then + stop 'Unable to set recv timout in zmq_to_qp_run_socket' + endif +END_PROVIDER + +BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_push ] + implicit none + BEGIN_DOC + ! Socket on which to push the results (1) + END_DOC + integer :: rc + character*(64) :: address + character*(8), external :: zmq_port + zmq_socket_push = f77_zmq_socket(zmq_context, ZMQ_PUSH) + address = trim(qp_run_address)//':'//zmq_port(1) + rc = f77_zmq_connect(zmq_socket_push, trim(address)) + if (rc /= 0) then + stop 'Unable to connect zmq_socket_push' + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_pull ] + implicit none + BEGIN_DOC + ! Socket which pulls the results (2) + END_DOC + integer :: rc + character*(64) :: address + character*(8), external :: zmq_port + zmq_socket_pull = f77_zmq_socket(zmq_context, ZMQ_PULL) + address = 'tcp://*:'//zmq_port(2) + rc = f77_zmq_bind(zmq_socket_pull, trim(address)) + if (rc /= 0) then + stop 'Unable to connect zmq_socket_pull' + endif + +END_PROVIDER + + + From cf12806fc72cc97fb952b818d913d56de2daed4e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Dec 2015 16:45:41 +0100 Subject: [PATCH 6/6] Forgot file --- .../ao_bielec_integrals_in_map_slave.irp.f | 99 +++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f 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 new file mode 100644 index 00000000..7aa59c0d --- /dev/null +++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f @@ -0,0 +1,99 @@ +subroutine ao_bielec_integrals_in_map_slave + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Computes a buffer of integrals + END_DOC + + integer :: j,l,n_integrals + integer :: rc + character*(8), external :: zmq_port + integer(ZMQ_PTR) :: zmq_socket_req_inproc, zmq_socket_push_inproc + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) ) + + ! Sockets + zmq_socket_req_inproc = f77_zmq_socket(zmq_context, ZMQ_REQ) + rc = f77_zmq_connect(zmq_socket_req_inproc, 'inproc://req_rep') + if (rc /= 0) then + stop 'Unable to connect zmq_socket_req_inproc' + endif + + zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH) + rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull') + if (rc /= 0) then + stop 'Unable to connect zmq_socket_push_inproc' + endif + + + + rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0) + rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0) + + do while (l > 0) + rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0) + + do j = 1, l + if (ao_overlap_abs(j,l) < ao_integrals_threshold) then + cycle + endif + call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) + rc = f77_zmq_send( zmq_socket_push_inproc, n_integrals, 4, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push_inproc, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE) + rc = f77_zmq_send( zmq_socket_push_inproc, buffer_value, integral_kind*n_integrals, 0) + enddo + rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0) + enddo + + deallocate( buffer_i, buffer_value ) + + rc = f77_zmq_disconnect(zmq_socket_req_inproc, 'inproc://req_rep') +end + + +subroutine ao_bielec_integrals_in_map_collector + use map_module + use f77_zmq + implicit none + BEGIN_DOC +! Collects results from the AO integral calculation + END_DOC + + integer :: j,l,n_integrals + integer :: rc + character*(8), external :: zmq_port + integer(ZMQ_PTR) :: zmq_socket_pull_inproc + real(integral_kind), allocatable :: buffer_value(:) + integer(key_kind), allocatable :: buffer_i(:) + + allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) ) + + zmq_socket_pull_inproc = f77_zmq_socket(zmq_context, ZMQ_PULL) + rc = f77_zmq_bind(zmq_socket_pull_inproc, 'inproc://push_pull') + if (rc /= 0) then + stop 'Unable to connect zmq_socket_pull_inproc' + endif + + n_integrals = 0 + do while (n_integrals >= 0) + + rc = f77_zmq_recv( zmq_socket_pull_inproc, n_integrals, 4, 0) + if (n_integrals > -1) then + rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind*n_integrals, 0) + rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind*n_integrals, 0) + call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) + else + rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind, 0) + rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind, 0) + endif + + enddo + + rc = f77_zmq_unbind(zmq_socket_pull_inproc, 'inproc://push_pull') + + deallocate( buffer_i, buffer_value ) +end +